Coherent manipulation and quantum phase interference in a fullerene-based electron triplet molecular qutrit

High-spin magnetic molecules are promising candidates for quantum information processing because their intrinsic multiplicity facilitates information storage and computational operations. However, due to the absence of suitable sublevel splittings, their susceptibility to environmental disturbances and limitation from the selection rule, the arbitrary control of the quantum state of a molecular electron multiplet has not been realized. Here, we exploit the photoexcited triplet of C70 as a molecular electron spin qutrit with pulsed electron paramagnetic resonance. We prepared the system into 3-level superposition states characteristic of a qutrit and validated them by the tomography of their density matrices. To further elucidate the coherence of the operation and the nature of the system as a qutrit, we demonstrated the quantum phase interference in the superposition. The interference pattern is further interpreted as a map of possible evolution paths in the space of phase factors, representing the quantum nature of the 3-level system.


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 them [1][2][3] . To be specific, molecular systems can be easily assembled to reach higher nuclearities and Hilbert space dimensions 4 , 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 discussed 1,[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 scale 13 , fullerenes are uniquely advantageous because these allcarbon 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 spin 18,19 , the entanglement between electron and nuclear spins 20,21 , electron spin dipolar coupling 22,23 and the atomic clock transition 24 .
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 magnets 25 . Since then, a variety of high-spin systems with which researchers demonstrated coherent manipulations addressing the allowed transitions have been reported 14,[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 times 19,[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 system 37 . Herein, we use the photoexcitation of C 70 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.

Resonance spectrum and condition selecting
We used a glassy-state solid solution of C 70 in d 8 -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 (S 0 ) to a singlet excited state (S 1 ), after which a triplet state (T 1 ) is reached through intersystem crossing 38 . The involved energy levels are illustrated in Fig. 1. The splitting of the T 1 state in a magnetic field can be modeled by the spin Hamiltonian where μ B is the Bohr magneton, g is the g-tensor, B 0 is the external magnetic field,Ŝ is the spin operator and D is the ZFS tensor. By simulating the echo detected field sweep (EDFS) spectrum in Fig. 1  With high-field approximation we denote the sublevels with M S = +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 S 1 .
The principal axes of the D tensor in different molecules were randomly oriented with respect to B 0 . 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 D tilting about 40°from B 0 , corresponding to a D zz component of 60 MHz. This portion of molecules signaled at B 0 ±32.1 G, as marked by the vertical dashed lines in Fig. 1, and can be addressed by microwave frequencies f ± ¼ f 0 ± 90 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 D zz and find the two frequencies that can address them.
The longitudinal (T 1 ) and transverse (T 2 ) relaxation times have been reported to depend on the orientation of the molecules relative to the magnetic field 41 . We studied this dependence in our system. To excite different portions of molecules, a series of measurements were conducted at a fixed frequency, f 0 , and different magnetic fields. The inversion-recovery method gave T 1 = 3 -5 ms, as shown in Supplementary Fig. 1. The Hahn-echo method gave T 2 = 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 T 1 ≫ T 2 and T 2 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 orientation 32,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 f 0 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 f 0 , and pulses were generated at f ± via an arbitrary wave generator (AWG). The signals were detected at f 0 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 ψ 1 j i ¼ ðjþi þ jÀiÞ= ffiffi ffi 2 p , with its nonadjacent components not directly addressable by the allowed transitions, and ψ 2 j i ¼ ðjþi þ j0i þ jÀiÞ= ffiffi ffi 3 p , 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 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 ij k θ ð Þ, and α ¼ 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, 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. 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, M 0þ at f − , orM 0À at f + , defined as in Supplementary Note 1. The result would be When determining the diagonal terms, the transformation part was null, leaving ρ t = ρ. ρ ++ , ρ 00 , and ρ ÀÀ can be solved with two measurements and the condition Tr(ρ) = 1. To determine S 0þ x , S 0þ y , S 0À x and S 0À y , 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 4T 2 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 B 0 , 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 j0i, and the 2-level superposition states ψ þ ¼ ðjþi þ j0iÞ= ffiffi ffi Þj0i. The results were Note that ρ 0 as measured is not an appropriate density matrix, this is because the results TrM 0 ± ρ 0 À Á have to be defined ∓ 1=2 to determine the instrument responsivities at f − and f + , while the offdiagonal 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 j0i state with high quality.
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 derivative 19 . The most significant deviation from ideal values occurred with ρ þÀ 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 ψ 1 j i was twice as long as any other. Supplementary Fig. 4 shows the decay of jρ þÀ 1 j and jρ 0þ þ j, measured by the tomography method, with an increased delay after the preparation of the state. The coherence times are 2.4(5) μs for jρ 0þ þ j and 0.3(1) μs for jρ þÀ 1 j. 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 T 2 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 T 2 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 ψ 1 j i and ψ 2 j i, which is inspired by TPPI aiming at illustrating the phase evolution of an entangled system 20 . The pulse sequence was as illustrated in Fig. 3(a). The preparation pulses at f ± were shifted in phase by ϕ 0 ∓ 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 M 0+ , was measured. It was observed that ϕ 0+ and ϕ 0− interfere to make M 0+ , 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 ψ 1 j i and ψ 2 j i. In the ideal case, M 0þ 2 ϕ 0À ; ϕ 0þ À Á ¼ À 4 9 cos ϕ 0þ À 1 9 cos ϕ 0À À 4 9 cos ϕ 0þ þ ϕ 0À À Á (14) the derivation of which is given in Supplementary Note 3.
Note that the reversion and detection parts altogether can be perceived as the measuring of an observable, which is M 0+ 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 k 0± they represent, and are symmetric about (0,0) by principle. Components with k 0À ; k 0þ ð Þ¼± ð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 k 0À ; k 0þ ð Þ¼ð 0; ± 1Þ and (±1,0) (B and C) marks the superposition of |0〉 with |−〉, and |0〉 with |+〉. The perfect ψ 1 j i ¼ ðjþi þ jÀiÞ= ffiffi ffi 2 p will have A1 only, and ψ 2 j i ¼ ðjþi þ j0i þ jÀiÞ= ffiffi ffi 3 p 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 f 0 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 ψð0Þ j i¼ c þ jþi þ c 0 j0i þ c À jÀi of this qutrit is prepared by P 0À 0 θ 2 ð ÞP 0þ 0 θ 1 ð Þ, then its phase evolution Y.-X. Wang et al.
observed in a rotating frame with frequency f 0 is The derivation is given by Supplementary Note 4. The latter is a state that can also be prepared directly by À2π f0ÀfÀ ð Þ t θ 1 ð Þ. Therefore, the phase evolution can be viewed as going along a linear path in the ϕ 0À ; ϕ 0þ À Á plane, namely Along this line, M 0þ ϕ 0À ðtÞ; ϕ 0þ ðtÞ À Á 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 t 0À ; t 0þ ð Þwith Here, t 0± 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 M 0+ would fluctuate in this process. Nevertheless, the area aside from this line are not meaningless. With a proper ratio r ¼ f 0 À f À ð Þ = f þ À f 0 ð Þ >0 any point is reachable. The setting of r can be done by simply shifting B 0 (and f 0 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 ψ j i ¼ c þ jþi þ c 0 j0i þ c À jÀi can be made into any ψ 0 j i ¼ c þ expðiϕ 1 Þjþi þ c 0 j0i þ c À expðiϕ 2 ÞjÀi (ϕ 1 ,ϕ 2 > 0), by setting and B 0 ¼ hf 0 =ðμ B gÞ, and waiting for 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 d Summarized evolution paths with different r, each represented by a horizontal section of this graph. The vertical axis is spanned evenly by θ ¼ 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.
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 C 70 to demonstrate the coherent manipulation of a 3-level quantum system. Superposition states ψ 1 i j ¼ ðjþi þ jÀiÞ= ffiffi ffi 2 p and ψ 2 j i ¼ ðjþi þ j0i þ jÀiÞ= ffiffi ffi 3 p , 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 algorithm 25 . 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 ΔM S 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 D zz . 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 Mn 2+ and Gd 3+ 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 thoroughly 32,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.

Sample preparation
The sample was prepared by dissolving high purity C 70 (from 1-Material) in d 8 -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 highpower 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 T 2 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 T 1 measurements were done by the "inversion-recovery" method, and the Rabi oscillations were measured by nutation experiments. The intensities in T 1 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 B 0 = 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 GHz. The pulses addressing the selected portion of molecules are of frequencies f 0 ± 90 MHz, making the corresponding signals also ± 90 MHz off-resonant with respect to f 0 . A combination of wide-band preparation, transformation and reversion pulses and narrowband 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.