Introduction

Information processing directly utilizing the quantum nature of matter is among the most pursued technologies in this era. Among the physical platforms racing for this goal, magnetic molecules have been regarded as a promising candidate because of the vastly available chemistry-based bottom-up schemes to handle them1,2,3. To be specific, molecular systems can be easily assembled to reach higher nuclearities and Hilbert space dimensions4, and their energy levels can be tailored via chemical modification. This field has been attracting increasing attention, and its development, prospects and remaining challenges have been thoroughly discussed1,5,6,7,8,9,10,11,12. While an increasing number of candidate molecules showing sufficiently long spin relaxation times have been synthesized, pushing the limit to the millisecond scale13, fullerenes are uniquely advantageous because these all-carbon molecules have little intramolecular hyperfine interaction, which helps preserving the spin coherence, and the near-spherical symmetry of the cage configurations allows small zero-field splittings (ZFS)14,15,16,17. As a result, various studies related to fullerene-based quantum information processing (QIP) have been reported, studying the ultrafast control of nuclear spin qubits assisted by electron spin18,19, the entanglement between electron and nuclear spins20,21, electron spin dipolar coupling22,23 and the atomic clock transition24.

While most researches concerning QIP based on magnetic molecules presented different ways to use the spin as a qubit, the multilevel nature of S > 1/2 systems has not been thoroughly exploited yet. It has been proposed that the quantum searching algorithm can be implemented in high-spin molecular magnets25. Since then, a variety of high-spin systems with which researchers demonstrated coherent manipulations addressing the allowed transitions have been reported14,26,27,28,29,30,31,32. Nevertheless, the full control of a high-spin system requires superposing the spin states between nonadjacent levels.

Photoexcitation has been used to generate S > 1/2 systems for a long history, with many examples showing quantum entanglement or long coherence times19,33,34,35,36. Recent researches have demonstrated that a photogenerated radical pair can be utilized to perform quantum teleportation in a donor–acceptor stable radical system37. Herein, we use the photoexcitation of C70 to generate an initialized electron spin triplet and coherently manipulate it by pulsed electron paramagnetic resonance (EPR) experiments. We demonstrate the capability to arbitrarily control the system by preparing two superposition states, one regarding the forbidden transition and the other involving the full Hilbert space basis set. With the proved accessibility of all the three sublevels, we illustrate the interference between the quantum phases in the superposition states of this electron spin qutrit. Combining the interference with the concept of time proportional phase increment (TPPI)20, we further demonstrate that a spin qutrit can undergo an incommensurate evolution.

Results

Resonance spectrum and condition selecting

We used a glassy-state solid solution of C70 in d8-toluene, prepared as described in Methods. It was cooled to 5 K and irradiated by a beam of 532 nm laser. The molecules were excited from the singlet ground state (S0) to a singlet excited state (S1), after which a triplet state (T1) is reached through intersystem crossing38. The involved energy levels are illustrated in Fig. 1. The splitting of the T1 state in a magnetic field can be modeled by the spin Hamiltonian

$$\begin{array}{*{20}{c}} {\hat H = \mu _{\mathrm{B}}{\mathbf{B}}_0^{\mathrm{T}} \cdot \overline{\overline g} \cdot {\hat{\mathbf S}} + {\hat{\mathbf S}}^{\mathrm{T}} \cdot \overline{\overline D} \cdot {\hat{\mathbf S}}} \end{array}$$
(1)

where μB is the Bohr magneton, \(\overline{\overline g}\) is the g-tensor, B0 is the external magnetic field, \({\hat{\mathbf S}}\) is the spin operator and \(\overline{\overline D}\) is the ZFS tensor. By simulating the echo detected field sweep (EDFS) spectrum in Fig. 1 (inset) with EasySpin (5.2.25)39, an isotropic \(\overline{\overline g}\) with g = 2.0037 and a \(\overline{\overline D}\) with \(|D| = 152\,{\mathrm{MHz}}\) and \(|E| = 50.4\,{\mathrm{MHz}}\) are derived, which agree with previously reported values40. The spectrum was taken at \(f_0 = 9.2505\,{\mathrm{GHz}}\) and the central field was \(B_0 = hf_0/(\mu _{\mathrm{B}}g) = 3298.5\,{\mathrm{G}}\), where h is the Planck constant. The spin density of the triplet state is calculated and shown in Fig. 1. Its slightly elliptical shape rationalizes the moderate ZFS.

Fig. 1: Energy levels and EDFS spectrum.
figure 1

Energy levels involved in the experiments with the molecular structure and calculated spin density of the photoexcited triplet of C70. The inset is the measured and simulated EDFS spectrum. The vertical dashed lines mark the fields in which the measured (see next section) portion of molecules signal. Boxes mark the range in which the molecules are affected by the preparation, transformation and reversion pulses.

Full size image

With high-field approximation we denote the sublevels with MS = +1, 0 and −1 by |+〉, |0〉 and |−〉, respectively. In the EDFS spectrum, the downward and upward peaks correspond to resonant emission (|0〉 → |−〉) and absorption (|0〉 → |+〉), respectively. The centrosymmetric shape indicates equal population differences regarding transitions |0〉 → |±〉. Therefore, the photoexcitation initialized the system to a pseudopure state. With equal populations of |±〉, it can be treated equivalently as a pure |0〉 state even though the purity cannot be confirmed within EPR experiments. This is more explicitly evidenced by the density matrix tomography of the initial state described in the next section. The favouring of |0〉 is assumable since the triplet is reached via intersystem crossing from S1.

The principal axes of the \(\overline{\overline D}\) tensor in different molecules were randomly oriented with respect to B0. In the following procedure of quantum state manipulation, it is important that a well-defined portion of the molecules is addressed by the microwave pulses. Selected were those with their principal axes of \(\overline{\overline D}\) tilting about 40° from B0, corresponding to a Dzz component of 60 MHz. This portion of molecules signaled at B0±32.1 G, as marked by the vertical dashed lines in Fig. 1, and can be addressed by microwave frequencies \(f_ \pm = f_0 \pm 90\,{\mathrm{MHz}}\) (corresponding to transitions |0〉 → |∓〉, respectively). These frequencies proved suitable for our experiments because they were neither too close so that the signals would be too weak, nor too far apart so that they would fall out of the resonator bandwidth. Note that since the photoexcitation states of a fullerene can be complicated, the agreeing measured EDFS spectrum and simulation do not guarantee the uniformity of excitation. Nevertheless, some uncertainty regarding photoexcitation can be accepted because our concern is to pick out the molecules with a narrow range of Dzz and find the two frequencies that can address them.

The longitudinal (T1) and transverse (T2) relaxation times have been reported to depend on the orientation of the molecules relative to the magnetic field41. We studied this dependence in our system. To excite different portions of molecules, a series of measurements were conducted at a fixed frequency, f0, and different magnetic fields. The inversion-recovery method gave T1 = 3 – 5 ms, as shown in Supplementary Fig. 1. The Hahn-echo method gave T2 = 13 – 20 μs, as shown in Supplementary Fig. 2. The field dependence is within a factor of 1.6 and not considered as very significant compared to the 95% confidence intervals of the fittings. The fact that T1 ≫ T2 and T2 is much longer than typical pulse durations enables us to prepare, manipulate and measure superposition states coherently before the system dephases or apparently relaxes towards thermal equilibrium, regardless of molecular orientation.

According to earlier findings by our group, the Rabi frequency of a high-spin system with large ZFS is also dependent on the orientation32,42. However, the ZFS in this system is too small to have any obvious effect. The orientation dependence of the Rabi frequency was studied by measuring Rabi oscillations at f0 in different magnetic fields. As shown in Supplementary Fig. 3, no significant variation was observed.

Throughout the following experiments, the working frequency of the microwave bridge was fixed at f0, and pulses were generated at f± via an arbitrary wave generator (AWG). The signals were detected at f0 off-resonantly. Details are provided in Methods.

Superposition state preparation and tomography

To demonstrate the capability to reach an arbitrary superposition state and thus qualify the system as a qutrit, we prepared the system into \(\left| {\psi _1} \right\rangle = (| + \rangle + | - \rangle )/\sqrt 2\), with its nonadjacent components not directly addressable by the allowed transitions, and \(\left| {\psi _2} \right\rangle = (| + \rangle + |0\rangle + | - \rangle )/\sqrt 3\), with its components covering the full Hilbert space basis set, and validated the preparation by density matrix tomography, as detailed below. The superposition states were prepared as

$$\left| {\psi _1} \right\rangle = P_y^{0 - }\left( {{\uppi}} \right)P_{ - y}^{0 + }\left( \frac{{\uppi}}{2}\right) \left| 0 \right\rangle$$
(2)
$$\left| {\psi _2} \right\rangle = P_y^{0 - }\left( {\frac{{\uppi} }{2}} \right)P_{ - y}^{0 + }\left( \alpha \right)\left| 0 \right\rangle$$
(3)

Here, we denote the unitary transformation implemented by a pulse addressing sublevels i and j (\(i,j = + ,0, -\)) at phase k (\(k = x,y, - x, - y\), or ϕ for an arbitrary phase angle) with tipping angle θ as \(P_k^{ij}\left( \theta \right)\), and \(\alpha = {\mathrm{arccos}}(1/3)\). Figure 2(a) illustrates the preparation and tomography procedure, where the latter contains a transformation part followed by the measurement. After laser initialization, the system was manipulated with the preparation pulses. After that, a sequence of operations was applied to determine the real and imaginary parts of each density matrix term. In general, for a 3-by-3 density matrix,

$$\begin{array}{*{20}{c}} {{\boldsymbol{\rho }}= \left( {\begin{array}{*{20}{c}} {\rho ^{ + + }} & {\rho ^{0 + }} & {\rho ^{ + - }} \\ {\overline {\rho ^{0 + }} } & {\rho ^{00}} & {\rho ^{0 - }} \\ {\overline {\rho ^{ + - }} } & {\overline {\rho ^{0 - }} } & {\rho ^{ - - }} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {\rho ^{ + + }} & {S_x^{0 + } - {\mathrm{i}}S_y^{0 + }} & {S_x^{ + - } - {\mathrm{i}}S_y^{ + - }} \\ {S_x^{0 + } + {\mathrm{i}}S_y^{0 + }} & {\rho ^{00}} & {S_x^{0 - } - {\mathrm{i}}S_y^{0 - }} \\ {S_x^{ + - } + {\mathrm{i}}S_y^{ + - }} & {S_x^{0 - } + {\mathrm{i}}S_y^{0 - }} & {\rho ^{ - - }} \end{array}} \right)} \end{array}$$
(4)

where the bar on the top means complex conjugate, and all the variables on the right side of the latter equal sign are real numbers.

Fig. 2: Tomography of the density matrix.
figure 2

a Pulse sequence for the preparation and density matrix tomography of \(\left| {\psi _1} \right\rangle\) and \(\left| {\psi _2} \right\rangle\). b–e Real and imaginary parts of the density matrices of (b, c) \(\left| {\psi _1} \right\rangle\) and (d, e) \(\left| {\psi _2} \right\rangle\), with deviations from the ideal values represented by the dashed yellow bars.

Full size image

The prepared state’s density matrix ρ was first transformed to ρt, in which the to-be-measured part appears in the diagonal. Then ρt is measured in one of its 2-level subspace magnetizations, \(\hat M^{0 + }\) at f−, or \(\hat M^{0 - }\) at f+, defined as in Supplementary Note 1.

The result would be

$$\begin{array}{*{20}{c}} {{\mathrm{Tr}}\left( {\hat M^{0 + }{\boldsymbol{\rho }}_{\mathrm{t}}} \right) = \left( {{\boldsymbol{\rho }}_t^{ + + } - {\boldsymbol{\rho }}_t^{00}} \right)/2} \end{array}$$
(5)

or

$$\begin{array}{*{20}{c}} {{\mathrm{Tr}}\left( {\hat M^{0 - }{\boldsymbol{\rho }}_{\mathrm{t}}} \right) = \left( {{\boldsymbol{\rho }}_t^{00} - {\boldsymbol{\rho }}_t^{ - - }} \right)/2} \end{array}$$
(6)

When determining the diagonal terms, the transformation part was null, leaving ρt = ρ. ρ++, ρ00, and \(\rho ^{ - - }\)can be solved with two measurements and the condition Tr(ρ) = 1. To determine \(S_x^{0 + }\), \(S_y^{0 + }\), \(S_x^{0 - }\) and \(S_y^{0 - }\), a π/2 pulse at f− or f+ was needed in the transformation part to transfer each of them to the diagonal terms of ρt. \(S_x^{ + - }\) and \(S_y^{ + - }\) cannot be transferred to the diagonal terms by a single pulse. They were instead transferred by a sequence of two pulses. The off-diagonal terms were measured as described in Supplementary Note 2. To determine the density matrix of a prepared state, ten experiments in total were run.

The measurement was done in three steps: (i) a delay for about 4T2 to let the system decohere, (ii) a selective π/2 pulse at either f+ or f−, and (iii) recording the free induction decay (FID) thereafter. The signed intensity of the component in the FID signal with the same frequency as the selective pulse was taken as the result.

A combination of short preparation and transformation pulses and long measuring pulses was used to optimize the results. The preparation and transformation pulses were 25.5 – 51.0 ns square pulses, which means they had bandwidths up to 78 MHz. This converted to the width of the resonance magnetic field is 28 G. Molecules that contributed to the EDFS spectrum in this range, as boxed in Fig. 1, were efficiently affected by the preparation and transformation pulses. Note that the two intervals are symmetric about B0, therefore correspond to the two transitions in the same portion of molecules. Although not very selective, they do not overlap to complicate the situation. The measuring pulses were 2500 ns square pulses. With a bandwidth of 0.8 MHz, they ensured that the recorded FID come from a very selective portion of molecules, corresponding to intervals narrower than the lines in the figure. To minimize the loss of fidelity, the delay between all the preparation and transformation pulses were kept within 2 ns, so that the decoherence affecting the off-diagonal terms happened mainly in the duration of the pulses.

To validate the pulses as quantitative operations and estimate the level of inaccuracy, we conducted tomography on ρ0, ρ+ and ρ−, namely the density matrices of the initial state, which is supposed to be \(|0\rangle\), and the 2-level superposition states \(\left| {\psi _ + } \right\rangle = (| + \rangle + |0\rangle )/\sqrt 2 = P_{ - y}^{0 + }\left( {{\uppi} /2} \right)|0\rangle\) and \(\left| {\psi _ - } \right\rangle = (|0\rangle + | - \rangle )/\sqrt 2 = P_y^{0 - }\left( {{\uppi} /2} \right)|0\rangle\). The results were

$$\begin{array}{*{20}{c}} {{\boldsymbol{\rho }}_0 = \left( {\begin{array}{*{20}{c}} 0 & {0.02 + 0.01{\mathrm{i}}} & {0.01 + 0.03{\mathrm{i}}} \\ {0.02 - 0.01{\mathrm{i}}} & 1 & { - 0.008 + 0.007{\mathrm{i}}} \\ {0.01 - 0.03{\mathrm{i}}} & { - 0.008 - 0.007{\mathrm{i}}} & 0 \end{array}} \right)} \end{array}$$
(7)
$$\begin{array}{*{20}{c}} {{\boldsymbol{\rho }}_ + = \left( {\begin{array}{*{20}{c}} {0.48} & {0.43 + 0.06{\mathrm{i}}} & { - 0.07 + 0.03{\mathrm{i}}} \\ {0.43 - 0.06{\mathrm{i}}} & {0.53} & { - 0.01 - 0.02{\mathrm{i}}} \\ { - 0.07 - 0.03{\mathrm{i}}} & { - 0.01 + 0.02{\mathrm{i}}} & { - 0.01} \end{array}} \right)} \end{array}$$
(8)
$$\begin{array}{*{20}{c}} {{\boldsymbol{\rho }}_ - = \left( {\begin{array}{*{20}{c}} {0.04} & { - 0.003 + 0.02{\mathrm{i}}} & { - 0.02 - 0.002{\mathrm{i}}} \\ { - 0.003 - 0.02{\mathrm{i}}} & {0.46} & {0.39 + 0.07{\mathrm{i}}} \\ { - 0.02 + 0.002{\mathrm{i}}} & {0.39 - 0.07{\mathrm{i}}} & {0.50} \end{array}} \right)} \end{array}$$
(9)

Note that ρ0 as measured is not an appropriate density matrix, this is because the results \({\mathrm{Tr}}\left( {\hat M^{0 \pm }{\boldsymbol{\rho }}_0} \right)\) have to be defined \(\mp 1/2\) to determine the instrument responsivities at f− and f+, while the off-diagonal terms cannot be strictly zero. Nevertheless, the norms of off-diagonal terms not exceeding 0.035 indicate from another point of view that the initial state is a pseudopure \(|0\rangle\) state with high quality.

ρ+ and ρ− agree quite well with the ideal density matrices with levels of inaccuracy similar with ρ0, validating the control over the two transitions separately. The fidelities, defined as

$$\begin{array}{*{20}{c}} {F\left( {{\boldsymbol{\rho }}_i,{\boldsymbol{\sigma }}_i} \right) = \left[ {{\mathrm{Tr}}\left( {\sqrt {\sqrt {{\boldsymbol{\sigma }}_i} {\boldsymbol{\rho }}_i\sqrt {{\boldsymbol{\sigma }}_i} } } \right)} \right]^2} \end{array}$$
(10)

where σi are the ideal density matrices, are \(F\left( {{\boldsymbol{\rho }}_ + ,{\boldsymbol{\sigma }}_ + } \right) = 0.93\) and \(F\left( {{\boldsymbol{\rho }}_ - ,{\boldsymbol{\sigma }}_ - } \right) = 0.87\).

The density matrices of \(\left| {\psi _1} \right\rangle\) and \(\left| {\psi _2} \right\rangle\) are determined to be

$$\begin{array}{*{20}{c}} {{\boldsymbol{\rho }}_1 = \left( {\begin{array}{*{20}{c}} {0.50} & { - 0.05 + 0.02{\mathrm{i}}} & {0.29} \\ { - 0.05 - 0.02{\mathrm{i}}} & {0.04} & { - 0.02 + 0.04{\mathrm{i}}} \\ {0.29} & { - 0.02 - 0.04{\mathrm{i}}} & {0.45} \end{array}} \right)} \end{array}$$
(11)
$$\begin{array}{*{20}{c}} {{\boldsymbol{\rho }}_2 = \left( {\begin{array}{*{20}{c}} {0.37} & {0.27 + 0.03{\mathrm{i}}} & {0.28 + 0.04{\mathrm{i}}} \\ {0.27 - 0.03{\mathrm{i}}} & {0.31} & {0.29 + 0.01{\mathrm{i}}} \\ {0.28 - 0.04{\mathrm{i}}} & {0.29 - 0.01{\mathrm{i}}} & {0.33} \end{array}} \right)} \end{array}$$
(12)

as shown in Fig. 2(b–e). The fidelities are \(F\left( {{\boldsymbol{\rho }}_1,{\boldsymbol{\sigma }}_1} \right) = 0.77\) and\(F\left( {{\boldsymbol{\rho }}_2,{\boldsymbol{\sigma }}_2} \right) = 0.89\). Since \({\mathrm{Tr}}\left( {{\boldsymbol{\rho }}_i^2} \right) < 1\) (\(i = + , - ,1,2\)), all the four prepared states were apparently mixed.

While previous reports of manipulating the electron triplet state in a molecular system are rare, this fidelity is considerably better than that in a study with similar material, where a pair of nuclear spins are made into a Bell state using the photoexcited triplet of a fullerene derivative19. The most significant deviation from ideal values occurred with \(\rho _1^{ + - }\). This might have resulted mainly from that the nonadjacent superposition decayed much faster than the adjacent superpositions and that the π pulse used to prepare \(\left| {\psi _1} \right\rangle\) was twice as long as any other. Supplementary Fig. 4 shows the decay of \(|\rho _1^{ + - }|\) and \(|\rho _ + ^{0 + }|\), measured by the tomography method, with an increased delay after the preparation of the state. The coherence times are 2.4(5) μs for \(|\rho _ + ^{0 + }|\) and 0.3(1) μs for \(|\rho _1^{ + - }|\). The imperfection of the pulses might have been another source of error. In this experiment, preparation and transformation pulses were chosen short to minimalize decoherence, but short pulses are disadvantageous in their broad bandwidth. To further improve the fidelity, a finer balance between the two factors might be needed.

With the validated ability to prepare an arbitrary superposition state, we studied the T2 relaxation of superposition states with different proportions of |+〉, |0〉 and |-〉 with the pulse sequence “preparation-τ-π-τ-echo”. The decreasing echo intensities are plotted versus 2τ and fitted with exponential decay, as shown in Supplementary Fig. 5. The T2 of different superposition states showed no significant difference with respect to the 95% confidence interval of the fitting. The π pulses and echoes herein were at f−, and the same conclusion is expected for f+.

Quantum phase interference

To further validate the coherent control over this 3-level system and elucidate the system’s characteristic as a qutrit, we conducted an experiment to unveil the quantum phase interference of the phase factors in \(\left| {\psi _1} \right\rangle\) and \(\left| {\psi _2} \right\rangle\), which is inspired by TPPI aiming at illustrating the phase evolution of an entangled system20. The pulse sequence was as illustrated in Fig. 3(a). The preparation pulses at f± were shifted in phase by \(\phi ^{0 \mp }\) respectively, generating superposition states with the same amplitudes and different phase factors. Then the reversion pulses, which are reversed in order, have opposite tipping angles with respect to the preparation pulses and fixed in phases were applied. Afterwards, the 2-level subspace magnetization M0+, was measured. It was observed that ϕ0+ and ϕ0− interfere to make M0+, as a function of the two phases, exhibit a pattern characteristic of the superposition proportion. Figure 3(b, c) show these patterns of the prepared states (left) and simulated results (right) for \(\left| {\psi _1} \right\rangle\) and \(\left| {\psi _2} \right\rangle\). In the ideal case,

$$\begin{array}{*{20}{c}} {M_1^{0 + }\left( {\phi ^{0 - },\phi ^{0 + }} \right) = - \frac{1}{2}\cos \left( {\phi ^{0 + } + \phi ^{0 - }} \right)} \end{array}$$
(13)
$$\begin{array}{*{20}{c}} {M_2^{0 + }\left( {\phi ^{0 - },\phi ^{0 + }} \right) = - \frac{4}{9}\cos \phi ^{0 + } - \frac{1}{9}\cos \phi ^{0 - } - \frac{4}{9}\cos \left( {\phi ^{0 + } + \phi ^{0 - }} \right)} \end{array}$$
(14)

the derivation of which is given in Supplementary Note 3.

Fig. 3: Quantum phase interference.
figure 3

a The pulse sequence. b, c Measured (left) and simulated (right) patterns of M0+ as a function of ϕ0+ and ϕ0− for \(\left| {\psi _1} \right\rangle\) and \(\left| {\psi _2} \right\rangle\). d, e The 2D Fourier transformations of interference patterns (b) and (c), with solid bars showing the experimental data and dashed ones the simulation results. For clarity, only low-frequency data are shown and the heights are normalized.

Full size image

Note that the reversion and detection parts altogether can be perceived as the measuring of an observable, which is M0+ transformed by a unitary operator representing the reversion pulses, as shown in Supplementary Note 3. Therefore, although phase factors in a superposition state are not observables by themselves, they do exhibit this internal interference regarding a well-defined observable.

The 2D Fourier transforms (2D-FFT) in Fig. 3(d, e) show the components in the patterns’ periodicity. The solid-line bars represent data from experiments and dashed-line bars give calculated results, with the magnitudes normalized. They are indexed with the dimensionless frequencies k0± they represent, and are symmetric about (0,0) by principle. Components with \(\left( {k^{0 - },k^{0 + }} \right) = \pm (1,1)\) (A1 and A2) are an indication of the superposition of |−〉 and |+〉 states, which corresponds to a spin-forbidden transition and does not exist in 2-level systems. Components with \(\left( {k^{0 - },k^{0 + }} \right) = (0, \pm 1)\) and (±1,0) (B and C) marks the superposition of |0〉 with |−〉, and |0〉 with |+〉. The perfect \(\left| {\psi _1} \right\rangle = (| + \rangle + | - \rangle )/\sqrt 2\) will have A1 only, and \(\left| {\psi _2} \right\rangle = (| + \rangle + |0\rangle + | - \rangle )/\sqrt 3\) will have A2, B, and C with relative heights 4:4:1, as indicated by Eqs. (13, 14). These were reasonably reproduced by the experiment. The deviation can be attributed to the same reasons as those discussed for tomography.

In the common practice of pulsed magnetic resonance, the physical process is studied in a rotating frame with the working frequency f0 to simplify the description. However, speaking of a qudit in which a number of transitions are simultaneously addressable with different frequencies, the phases of the superposition of different components evolve at different rates. Suppose a superposition state \(\left| {\psi (0)} \right\rangle = c_ + | + \rangle + c_0|0\rangle + c_ - | - \rangle\) of this qutrit is prepared by \(P_0^{0 - }\left( {\theta _2} \right)P_0^{0 + }\left( {\theta _1} \right)\), then its phase evolution observed in a rotating frame with frequency f0 is

$$\begin{array}{lll}\left|{\psi\left(t\right)}\right\rangle&=&\hat R\left(t\right)\hat E\left(t\right)\left|\psi \left( 0\right)\right\rangle\\ &&=c_+\exp\left( {{\mathrm{i}}\cdot2{{\uppi}}\left({f_0-f_-}\right)t}\right)|+\rangle\\ && +c_0\left|0 \right\rangle+c_-\exp\left( {{\mathrm{i}}\cdot 2{{\uppi}}\left({f_+-f_0}\right)t}\right)|-\rangle\end{array}$$
(15)

The derivation is given by Supplementary Note 4. The latter is a state that can also be prepared directly by \(P_{2{{\uppi}} \left( {f_ + - f_0} \right)t}^{0 - }\left( {\theta _2} \right)P_{ - 2{{\uppi}} \left( {f_0 - f_ - } \right)t}^{0 + }\left( {\theta _1} \right)\). Therefore, the phase evolution can be viewed as going along a linear path in the \(\left( {\phi ^{0 - },\phi ^{0 + }} \right)\) plane, namely

$$\begin{array}{*{20}{c}} {\left\{ {\begin{array}{*{20}{c}} {\phi ^{0 - }(t) = 2{{\uppi}} \left( {f_ + - f_0} \right)t} \\ {\phi ^{0 + }(t) = - 2{{\uppi}} \left( {f_0 - f_ - } \right)t} \end{array}} \right.} \end{array}$$
(16)

Along this line, \(M^{0 + }\left( {\phi ^{0 - }(t),\phi ^{0 + }(t)} \right)\) gives a function of time characterizing the evolution of the superposition state.

Therefore, the 2D interference pattern above can be viewed as a map of phase evolution by transforming the phase coordinates into time coordinates \(\left( {t^{0 - },t^{0 + }} \right)\) with

$$\begin{array}{*{20}{c}} {t^{0 \pm } = \frac{{\phi ^{0 \pm }}}{{2{{\uppi}} \left( {f_ \mp - f_0} \right)}}} \end{array}$$
(17)

Here, t0± can be viewed as fictitious evolution times, as if the phase evolution happened only in either of the 2-level subspaces and the phase of the third state is frozen. Apparently, the line \(t^{0 - } = t^{0 + }\) is the path of the actual evolution.

Taking Fig. 3(c) as an example, having its coordinates converted by Eq. (17) we get the upper panel of Fig. 4(a). Its diagonal marks the actual path of evolution through time \(t = t^{0 - } = t^{0 + }\). The lower panel shows how M0+ would fluctuate in this process. Nevertheless, the area aside from this line are not meaningless. With a proper ratio \(r = \left( {f_0 - f_ - } \right)/\left( {f_ + - f_0} \right) > 0\) any point is reachable. The setting of r can be done by simply shifting B0 (and f0 accordingly in the meantime). As long as f± are fixed, the portion of molecules that are addressed keeps the same. The interference pattern, rescaled in different ways, could cover cases with any positive r, as exemplified in Fig. 4(b) with r = 0.5 and 4(c) with r = 2. Generally, the evolution path of the system with any positive r can be represented by a horizontal line in Fig. 4(d). More generally, a state \(\left| \psi \right\rangle = c_ + | + \rangle + c_0|0\rangle + c_ - | - \rangle\) can be made into any \(\left| {\psi ^\prime } \right\rangle = c_ + {\mathrm{exp}}({\mathrm{i}}\phi _1)| + \rangle + c_0|0\rangle + c_ - {\mathrm{exp}}({\mathrm{i}}\phi _2)| - \rangle\) (ϕ1,ϕ2 > 0), by setting

$$\begin{array}{*{20}{c}} {f_0 = \frac{{f_ + \phi _1\, + \,f_ - \phi _2}}{{\phi _1\, + \,\phi _2}}} \end{array}$$
(18)

and \(B_0 = hf_0/(\mu _{\mathrm{B}}g)\), and waiting for

$$\begin{array}{*{20}{c}} {t = \frac{1}{{2{{\uppi}} \left( {f_ + - f_ - } \right)}}} \end{array}$$
(19)
Fig. 4: Interpretation of the quantum phase interference pattern as a map of evolution.
figure 4

a–c The interference pattern of \(\left| {\psi _2} \right\rangle\) converted to time coordinates \(\left( {t^{0 - },t^{0 + }} \right)\) with a r = 1, b r=0.5 and c r = 2. The red lines mark the path of evolution through time in each condition, along which the fluctuation of M0+ is shown in the lower panels, with solid and dashed lines giving experimental and simulated results. d Summarized evolution paths with different r, each represented by a horizontal section of this graph. The vertical axis is spanned evenly by \(\theta = {\mathrm{arctan}}(r)\). e Paths of evolution represented by curves on a torus, with blue, red and green curves for r = 0.5, 1, and 2 in the upper part, and the incommensurate case for an irrational r in the lower part. The tori are colored according to simulated data for clarity. The same color scale is shared by a–e.

Full size image

The collection of all these states, parameterized by the phases ϕ1 and ϕ2, can be represented by a torus, as shown in Fig. 4(e), where ϕ1 spans along the horizontal circles going around the torus and ϕ2 spans along the ones perpendicular to them. The path of evolution can then be viewed as a curve twining on the torus. When r is a rational number, the paths are closed, as shown in the upper part (blue, red and green curves for r = 0.5, 1 and 2, respectively). However, in the most cases where r is irrational, the paths would never return to where it started, as depicted in the lower part. This means that qutrits, unlike qubits, tend to undergo non-periodic evolution. An observable might naturally fluctuate along an incommensurate curve, although approximate periodicity in a finite period of time can be achieved if the frequencies are set precisely enough.

Discussion

Utilizing the advantages of the fullerene, we used the photoexcited triplet of C70 to demonstrate the coherent manipulation of a 3-level quantum system. Superposition states \(\left| {\psi _1} \rangle\right. = (| + \rangle + | - \rangle )/\sqrt 2\) and \(\left| {\psi _2} \right\rangle = (| + \rangle + |0\rangle + | - \rangle )/\sqrt 3\), characteristic of a qutrit, were prepared, characterized and observed by phase interference. The interference pattern is further explained as a map of evolution. This shows the capability to fully control this molecule-based S = 1 electron spin state, in both magnitudes and phases of superposition. This is a prerequisite for further applications such as Grover’s algorithm25. The loss of fidelity is mainly due to the faster decay of the nonadjacent superposition. This might be accounted for by that the dephasing of both of the two adjacent superpositions may contribute to the dephasing of this nonadjacent superposition, and that the twice-as-large ΔMS and level splitting makes it more sensitive to environmental inhomogeneity and fluctuations. Further eliminating the factors leading to decoherence is still important. In the current system and experimental setup, the main cause of decoherence is the diversity of the ZFS parameter Dzz. Other factors include field inhomogeneity and intermolecular dipole-dipole interactions. The former calls for instrumentational improvements, while the latter is only eliminated by dilution, at the price of lowering the signal-to-noise ratio. Therefore, crystalizing the sample to attain a more concentrated distribution of orientations is of interest, and chemical approaches such as embedding the fullerene molecules in rigid clathrates also seem promising.

The multilevel nature means not only the possibility of encoding more information within a molecular unit, but also the emergence of phenomena such as incommensurate evolution not present in 2-level systems. Arbitrary coherent manipulation of this S > 1/2 system is realized on a molecular and purely electron basis, which paves the way toward addressable, tailorable and scalable assemblies of multilevel quantum systems. This reminds us that high-spin magnetic molecules such as Mn2+ and Gd3+ complexes might deserve more attention. In these systems, although the initialization procedure might not be so straightforward as photoexcitation and suitable ZFS might be more difficult to find, the stability of the multiplets allows the level structure to be studied and manipulated by chemical approaches much more thoroughly32,43. It is expectable that molecular systems with multilevel nature can be exploited to a greater extent by coherent manipulation for QIP in the future.

Methods

Sample preparation

The sample was prepared by dissolving high purity C70 (from 1-Material) in d8-toluene (from J&K, 99.5% D) to form a 0.1 mg ml−1 solution. The solution was then degassed, sealed in argon atmosphere and cooled rapidly to 5 K by the liquid helium cryostat to form a glassy-state solid solution.

Equipment

The pulsed EPR experiments are done on a Bruker Elexsys E580 spectrometer with a MS-5W resonator with an optical window. The Q factor is 36, and the resonator broadening is greater than f+−f−. The pulses were applied through a 300 W high-power amplifier with 0 dB high-power attenuation. The fine control of the pulse amplitudes, frequency shifts and phases were realized with a Bruker SpinJet AWG. The photoexcitation was conducted with an MQV-350 laser from Beijing ZK Laser Co., Ltd.

Pulsed EPR experiments

All pulsed EPR experiments started with photoexcitation by a nanosecond laser pulse with 532 nm wavelength. The full width at half maximum was 7 ns and the energy is 200 mJ per pulse. The T2 measurements, as shown in Supplementary Figs. 2 and 5, were done by measuring the decay of the Hahn-echo intensity (integrated) while increasing τ in the “π/2-τ-π-τ-echo” or “preparation-τ-π-τ-echo” sequence and plotting the echo intensity against 2τ. The T1 measurements were done by the “inversion-recovery” method, and the Rabi oscillations were measured by nutation experiments. The intensities in T1 and Rabi oscillation measurements were both integrated Hahn-echo intensities.

The preparation and tomography of superposition states and the quantum phase interference experiments were done with a fixed static magnetic field B0 = 3298.5 G. The working frequency of the microwave bridge, being also the frequency about which the signals were recorded, is fixed at \(f_0 = 9.2505\,{\mathrm{GHz}}\). The pulses addressing the selected portion of molecules are of frequencies \(f_0 \pm 90\,{\mathrm{MHz}}\), making the corresponding signals also \(\pm 90\,{\mathrm{MHz}}\) off-resonant with respect to f0. A combination of wide-band preparation, transformation and reversion pulses and narrow-band measuring pulses were used to promote the precision of these experiments. The wide-band pulses were of large amplitudes and thus short durations, not exceeding 51 ns. Tens of MHz in bandwidth, they had high efficiency but poor selectivity. The narrow-band measuring pulses were 2500 ns and of low amplitudes to attain high selectivity. The designs and procedures of these experiments are elaborated in the main text.

Calculation

The spin density was calculated with B3LYP/6-31g* by Gaussian 09w. The simulation of the EDFS spectrum was done with EasySpin (5.2.25) MATLAB toolbox.

Data processing

In the tomography and the quantum phase interference experiments, all the raw data were FID curves, 2048 ns in length and 1 ns in time resolution. The fast Fourier transform was applied to extract from each FID curve the phase and intensity at the frequency of the preceding selective measuring pulse. This processing is done with MATLAB by scripts available from the authors upon reasonable request.