Supplementary Information for All-optical coherent population trapping with defect spin ensembles in silicon carbide

Divacancy defects in silicon carbide have long-lived electronic spin states and sharp optical transitions. Because of the various polytypes of SiC, hundreds of unique divacancies exist, many with spin properties comparable to the nitrogen-vacancy center in diamond. If ensembles of such spins can be all-optically manipulated, they make compelling candidate systems for quantum-enhanced memory, communication, and sensing applications. We report here direct all-optical addressing of basal plane-oriented divacancy spins in 4H-SiC. By means of magneto-spectroscopy, we fully identify the spin triplet structure of both the ground and the excited state, and use this for tuning of transition dipole moments between particular spin levels. We also identify a role for relaxation via intersystem crossing. Building on these results, we demonstrate coherent population trapping -a key effect for quantum state transfer between spins and photons- for divacancy sub-ensembles along particular crystal axes. These results, combined with the flexibility of SiC polytypes and device processing, put SiC at the forefront of quantum information science in the solid state.


SAMPLE DESIGN
We designed a sample structure that gave a high amount of PL and PLE emission, optimized PL(E) collection efficiency, while also geometrically separating the weak divacancy emission from the much stronger control laser beams (Fig. S1). Gold coatings of 100 nm thickness on the front and back of the sample act as mirrors that trap the applied laser beams in the sample. This enhances the path length of the laser in the sample, giving a higher number of divacancy defects that contribute to the emission signal. The laser beams enter the sample at the top-right corner (window of 0.3 mm by 0.2 mm without gold coating), while the sample has a 45 • forward tilt with respect to the laser beam. Due to refraction at the interface this gives propagation in the sample that is near-parallel to the c-axis, with a deviation that is just large enough for trapping the light on a non-overlapping zig-zag trajectory. We measured a 10% intensity loss at each reflection. In this manner laser light propagates until it reaches the left facet of the sample (no gold coating), denoted by an x, where most of the laser light still does not leave the sample due to total internal reflection. A small fraction does exit because the interfaces are not perfectly smooth. This approach give a suppression of 10 6 in laser light reaching our detectors. Conversely, divacancy emission is in random directions and for a large part on trajectories that do leave the sample when reaching the left facet. This includes trajectories that first get reflected inside the sample, and this is enhanced by total internal reflection at the trapezoidal shape on the right of the sample. This gives PL and PLE emission that mainly exits the sample on the left, in a cone of directions, facilitating efficient detection. We apply additional spectral filtering before our detectors for reaching an overall suppression of over 10 10 for light from the excitation lasers. Lasers are incident on the sample from the right, entering through a window of 0.3 mm by 0.2 mm.
Light propagates as shown, hardly exiting the sample at the marker x due to total internal reflection.
This results in a factor 10 6 suppression for laser light reaching PL(E) signal detectors. Light emitted by the divacancies does exit the sample, on the left. After additional spectral filtering it is collected by a single photon counter (or spectrometer). The trapezoidal shape of the sample, combined with gold coating on the front and back, further enhances this collection efficiency.

RATE MODEL OF TWO-LASER SPECTROSCOPY
We model the spectral positions and amplitudes of PLE signals for data as in Fig. 2a in an approach that combines rate equations for the transitions between the levels |g i and |e j with solving the spin Hamiltonian (equation (1)) for the ground and excited state. After fitting this yields results as in Fig. 2e. The rate equations are the set of coupled equations for the processes and parameters that are illustrated in Fig. S2, and this system's steady-state solution determines the PLE signals from our experiments with continuous excitation. We use an approach that only describes populations on the levels (and no coherences between the levels) in order to keep the total number of fitting parameters low. We repeated our experiments with variation of laser intensities over two orders of magnitude, and obtain nominally the same values for the fitting parameters (besides the trivial increase in PLE signal). This justifies that we can neglect coherences in our modeling.
The parameters for the rate equations, and the parameters D e and E e of the excitedstate spin Hamiltonian, are obtained from fitting the model to data sets as in Fig. 2a. We find parameters D g and E g for the ground-state spin Hamiltonian that are consistent with literature, and we therefore use these as fixed values. Similarly, the g-factor is set to 2 for both ground and excited state (deriving the g-factors from our experimental data gives In the latter case, we treat all six ISC rates X g i , X e j as fitting parameters which may slowly vary with magnetic field (further discussed below). During the fitting, we apply (in a self-consistent manner) that the homogeneous linewidth for the optical transitions is governed by the total decay rate Γ j out of a level |e j , which we calculate as Γ j = G 0 + X e j . The optical excitation rates E ij are modeled in a manner similar to G z ij , Here we only need to consider ZPL transitions since the phonon-side-band excitation is negligible for laser excitation near resonance with the ZPL at the temperatures of our experiments. For excitation from a level |g i to a level |e j where C 2 is a constant, f the laser frequency and f ij the resonance frequency of the transition. The last factor in this equation describes a Lorentzian lineshape of halfwidth Γ j /2π for the resonance (we checked that all our experiments were in the regime where the laser intensities give negligible power broadening of the PLE lines). We implement this with a fit parameter E 0 , The relaxation rates between ground states levels |g i are assumed to be negligible compared to the radiative decay and excitation rates (we checked that including these as parameters gave very low rates that do not improve the fits). For our fitting we subtract from the data a constant PLE background that is present from single laser driving of the inhomogeneous ensemble.

Fitting results and discussion
The fitting yields G 0 = XXX ± XXX MHz, D e = sequent on the s field 17 . field one for diva not exp ber of fi the resu data set that are We th six glob while m for the we effec ranges o states d show a at magn by notin and |e j states b in stabl The r for the (the app ues bet ror from these IS rates in |e X XX This ide side ban pulses ( levels, a results i ij FIG. S2: Schematic of energy levels and processes of the rate-equation model. E ij is the rate of laser excitation from |g i to |e j , G ij is the radiative decay rate from |e j to |g i , X e j is the intersystem-crossing (ISC) decay rate from |e j to |s , and X g i is the ISC decay rate from |s to |g i . agreement with published values to within the experimental error 1 ). For a given magnetic field, the energy splittings between the three levels |g i and |e j , and the associated spin states, are obtained by calculating eigenvalues and eigenvectors of equation (1) of the main text. The singlet state |s is fixed. We calculate in parallel PLE signals from the P and R sub-ensembles, but for displaying the result in Fig. 2e we use orange and green coloring for the respective signal contributions. For the small offset angles θ and ϕ we had in our experiments ( Fig. 1a) the differences in eigenvalues for different divacancy orientations within the sub-ensembles P and R are negligible for the present analysis.
We start with describing the influence of the spin states on the radiative decay from population in a level |e j to a level |g i . We first analyze this for the case of zero-phonon-line (ZPL) transitions (below we will further discuss the role of radiative decay via phononsideband emission), which occurs at a rate G z ij and is given by where C 1 is a constant and d ij is the transition dipole matrix element. We can write this out as whered is the electric dipole operator, and the states |ψ e j and |ψ g i are the initial and final state for this transition, respectively. We will use the approximation that |ψ i can be treated as a product state of the spin state |g i and the orbital state |χ g i of the ground state: Similarly, we will use |ψ e j = |χ e j |e j where |χ e j is the orbital state for the excited state. This approximation is valid for the case of very low spin-orbit interaction, and applies for our system given the low atomic masses of SiC. We thus also assume that the orbital state |χ g i is identical for the three levels |g i , and that |χ e j is identical for the three levels |e j . Further, we assume that at the magnetic fields of our experiments the highly confined orbital states |χ g i and |χ e j have negligible magnetic field dependence. Since the dipole operator only works on the orbital part of a state, the expression for G z ij can now be written as and with the above assumptions | χ g i |d |χ e j | is equal for all nine |e j -|g i transitions. We can thus express G z ij as Here G z 0 is the total decay rate via ZPL emission out of the excited state, and | g i |e j | 2 is the Franck-Condon factor with respect to spin 2 .
We need to account for the fact that a significant part of the radiative decay occurs via phonon-sideband emission. The ODMR work on the divacancies in SiC 3,4 shows that after phonon-sideband emission the spin state is preserved during subsequent phonon relaxation (which occurs at a fast time scale). Also, with the above assumptions, the rate of phononsideband emission is simply proportional to the ZPL emission for all |e j -|g i transitions.
We can thus simply express the total decay rate from all the possible radiative transitions as rates G ij (we drop the superscript z) that obey and the PLE signal for a particular |e j -|g i transition is proportional with G ij and the steady-state population in |e j . We treat G 0 as a fit parameter.
In addition to this radiative decay, we include parallel, non-radiative decay paths via the singlet state |s for modeling intersystem-crossing (ISC) events. Each |e j has its own transition rate X e j into |s , and from there the spin decays with rates X g i to |g i . We cannot exclude that the ISC rates depend on the spin state, and hence on magnetic field. We investigated fitting both with constant (field-independent) and field-dependent ISC rates.
In the latter case, we treat all six ISC rates X g i , X e j as fitting parameters which may slowly vary with magnetic field (further discussed below). During the fitting, we apply (in a selfconsistent manner) that the homogeneous linewidth for the optical transitions is governed by the total decay rate Γ j out of a level |e j , which we calculate as Γ j = G 0 + X e j . The optical excitation rates E ij are modeled in a manner similar to G z ij . Here we only need to consider ZPL transitions since the phonon-sideband excitation is negligible for laser excitation near resonance with the ZPL at the temperatures of our experiments. For excitation from a level |g i to a level |e j where C 2 is a constant, f the laser frequency and f ij the resonance frequency of the transition.
The last factor in this equation describes a Lorentzian lineshape of half-width Γ j /2π for the resonance (we checked that all our experiments were in the regime where the laser intensities give negligible power broadening of the PLE lines). We implement this with a fit parameter The relaxation rates between ground state levels |g i are assumed to be negligible compared to the radiative decay and excitation rates (we checked that including these as parameters gave very low rates that do not improve the fits). For our fitting we subtract from the data a constant PLE background that is present from single laser driving of the inhomogeneous ensemble.

Fitting results and discussion
The fitting yields G 0 = 20 ± 5 MHz, D e = 0.95 ± 0.02 GHz, and E e = 0.48 ± 0.01 GHz (and we used literature values 4 D g = 1.334 GHz and E g = 18.7 MHz). The value for E 0 is always such that it gives rates E ij well below G 0 and shows no interdependence with the other parameters. For these parameters only G 0 shows a weak dependence on the various approaches to fitting the ISC rates. Attempts to fit without accounting for ISC yields PLE lines at the correct frequencies (mainly governed by the values of the D e,g and E e,g parameters) and with the correct width (mainly governed by G 0 ) in Fig. 2a. However, for that case the relative amplitudes of calculated PLE lines (at a certain magnetic field value) show poor agreement with the experimental results. On this aspect we get much better agreement in an approach with magnetic-field-independent values for the six ISC rates (Fig. S3). However, we cannot exclude that the ISC rates depend on magnetic field.
While the ISC processes are far from completely understood, recent theory work for diamond NV − centers proposes a role for spin-orbit coupling, where the orbital quantization axis is linked to the defect orientation. Consequently, this work proposes that the ISC rates depend on the spin states of the levels, and thereby on magnetic field 5 . This also indicates that for non-zero magnetic field one should in fact fit with different ISC parameters for divacancies along different crystal directions. We did not explore this latter aspect since this makes the number of fitting parameters very high and we analyzed that the results would not have statistical significance for our data set. For non-zero field we thus present ISC rates that are an average over the different defect orientations.
We thus carried out a fitting approach where we used six global ISC rates that can vary as a function of field, while maintaining single parameters D e , E e , G 0 and E 0 for the full data set. To suppress statistical fluctuations we effectively fit average ISC rates over magnetic field ranges of 10 mT width in Fig. 2a (over which the spins states do not strongly vary). We then see ISC rates that show a weak variation with field, and become constant at magnetic fields over 40 mT. This can be understood by noting that this is beyond any anticrossings in the |g i and |e j levels where the spin states mix. Here, the spin states become field independent and this results indeed in stable ISC rates from the fitting.
The results for the ISC rates are presented in Table I for the applied field near B = 0 mT and B > 40 mT (the approach with field-independent ISC rates gave values between these two cases, with a least-squares error from fitting that was 20% higher). The result with these ISC rates are also presented in Fig 2e. At low fields, these ISC rates indicate relatively strong values for ISC out of levels |e m and |e u , and preferential ISC decay into the level |g l . This identifies how single-laser excitation into the phonon sideband or driving the ZPL with spectrally broad laser pulses (which both give parallel excitation from all spin levels, and is for example applied in ODMR studies 3,4 ) results in spin polarization in the ground state.

MODELING OF COHERENT POPULATION TRAPPING
We investigated the agreement between our CPT observations (Fig. 3) and the standard theoretical description for this phenomenon 6 (unlike the PLE modeling of the previous section, this modeling does include coherences between the levels). This also allows for extracting the ground-state spin dephasing time T * 2 for our ensemble. For making the fits for panels I-V in Fig. 3, we approximate the six-level system as an effective three-level system, where we only account for the levels that participate in a laser-driven Λ configuration (Fig. 2d). For this system we determine the steady-state solution of the density-matrix master equation, with the approach described in Ref. 6. In order to relate it to the measured data from our two-laser spectroscopy technique, we fit using a sum of master equation solutions, where the summation is over a range of inhomogeneous broadening values for the optical transition. (giving that the CPT dip is always in the middle of a PLE spectral line).
In addition, we take into account the decaying power along the beam path in our multipass sample (using a decay constant that was measured), by fitting with a sum of master equation solutions with exponentially decreasing Rabi frequencies. For the measurements that show two (four) CPT dips due to small misalignment angles for the applied magnetic field, we sum the PLE from two (four) contributions. Here each contribution is calculated for one set of misalignment angles that occur within the ensemble that is addressed (representing one of the two or four orientations within ensemble P or R, respectively). In order to obtain good fits we need to account for a background signal in the PLE data, in particular for PLE lines with CPT that are on the wing of a neighboring PLE line. In the calculations we account for this with a linear-in-frequency background signal. We can thus calculate the CPT line shapes that appear in the PLE data, and the results are presented as the solid lines in Fig. 3. These results give a T * 2 value of 42±8 ns.

ANGULAR DEPENDENCE OF CPT
To find evidence for our interpretation that the splitting of the CPT dips in panels III-V of Fig. 3 is a result of minute magnetic field misalignment, we rotated the sample holder inside the cryostat in steps of 0.5 • . This is a rotation around the vertical axis of Fig. S1, and yields a combined variation of ϕ and θ offsets due to the 45 • forward tilt in our sample mounting, see Fig. 1a and Fig. S1. Results from these measurement are presented in Fig. S4a

BLEACHING AND CHARGE STATE DYNAMICS
When addressing the divacancies with lasers resonant on the zero-phonon line (ZPL), the amount of PLE signal that is detected gradually decreases over time, on a timescale of minutes to hours for the laser powers we use. When keeping a single laser fixed central on the ZPL for several hours, and subsequently scanning a single laser over the entire ZPL (in a few seconds), we see a clear dip in the PLE signal (Fig. S5). This bleaching phenomenon for SiC divacancies is an effect also seen for NV − centers in diamond, where it has been identified as a two-photon process. In case of an NV − center, it is initially negatively charged, absorbs two photons, and subsequently an electron is lost to the conduction band 7 . For the SiC divacancy the initial state is charge neutral, and two-photon absorption also causes the defect to make a transition to a different meta-stable charge state. In this charged state the wavelength for the ZPL optical transitions is shifted by over 50 nm. Applying an additional laser that excites these charged divacancies brings them back to their original charge state (again in analogy with observations on NV − centers 7 ). We found that applying a 300 µW repumping laser at 685 nm completely removes the bleaching dip after a few seconds, and we had such a repump laser permanently on during all the two-laser spectroscopy experiments.
It is important that this laser is far-removed from the ZPL and the phonon sideband for excitation of our neutral divacancy, such that it only excites divacancies in the altered charge state. This avoids that the repump laser interferes with driving optical transitions for the neutral V SiC and our two-laser spectroscopy. Notably, with selectively bleaching most of the left and right wing of the inhomogeneous linewidth of the ZPL of our neutral divacancy (the inverse of the result in Fig. S5) one can create a long-lived sub-ensemble with a strongly reduced inhomogeneity for the optical transition.