Microscopic model of spin flip-flop processes in rare-earth-ion-doped crystals

Flip-flop processes due to magnetic dipole-dipole interaction between neighbouring ions in rare-earth-ion-doped crystals is one of the mechanisms of relaxation between hyperfine levels. Modeling of this mechanism has so far been macroscopic, characterized by an average rate describing the relaxation of all ions. Here however, we present a microscopic model of flip-flop interactions between individual nuclear spins of dopant ions. Every ion is situated in a unique local environment in the crystal, where each ion has different distances and a unique orientation relative to its nearest neighbors, as determined by the lattice structure. Thus, each ion has a unique flip-flop rate and the collective relaxation dynamics of all ions in a bulk crystal is a sum of many exponential decays, giving rise to a distribution of rates rather than a single average decay rate. We employ this model to calculate flip-flop rates in Pr$^{3+}$:Y$_2$SiO$_5$ and show experimental measurements of population decay of the ground state hyperfine levels at $\sim$2 K. We also present a new method to measure rates of individual transitions from hole burning spectra that requires significantly fewer fitting parameters in theoretical rate equations compared to earlier work. Furthermore, we measure the effect of external magnetic field on the flip-flop rates and observe that the rates slow down by two orders of magnitude in a field of 5 - 10 mT.

Even though hyperfine lifetimes can be as long as seconds or much more, relaxation can be a problem in many of the above applications. It results in decreased absorption depth, leading to lower efficiency of echoes in quantum memories, degrading of spectral filters and decreased gate fidelity in quantum computing. In general, hyperfine relaxation can occur either via lattice vibrations mediated by phonons (Spin-Lattice Relaxation) or via interactions with neighbouring spins (Spin-Spin Relaxation). Spin-Lattice Relaxation processes are well understood [21,22] and the mechanism relevant at cryogenic temperatures is the Direct process, whose rate increases proportional to the temperature (∝ T ) and square of magnetic field (∝ B 2 ) [23]. However, various experiments at cryogenic temperatures have demonstrated a decrease in relaxation rates with the application of a magnetic field [24][25][26][27]. The mechanisms responsible for relaxation in such cases are not phonon related but magnetic dipole interactions between dopants. They are known as flipflop interactions whereby two nearby ions exchange their spins via magnetic dipole-dipole interaction and the interaction strength wanes with distance r as r −6 . Studies in Kramers ions like Er 3+ and Nd 3+ have used a macroscopic model to explain spectral hole decay due to flipflop process by taking a single average rate to be related to the dopant concentration [27] and an average ion-ion distance for the ions in the crystal [26], resulting in a rate R ∝ n 2 r 6 . However, it has also been reported that this mechanism can lead to non-exponential decays [28] and the focus of this work is to develop a model that captures this effect.
Each dopant ion in the crystal is randomly placed in the crystalline structure such that it experiences a different magnetic environment and has different distances to and orientations of its nearest neighbouring dopants. Therefore, each ion relaxes with a unique rate and when the relaxation dynamics is studied in a bulk crystal, we see a sum of many exponential decays. In this work, we use a numerical simulation of a host crystal to create a distribution of ion-ion distances. The flip-flop rate between all the pairs of ions is then calculated using Fermi's Golden rule. The shape of this distribution of flip-flop rates mimics that of r −6 , where r is the ion-ion distance. We compare the model to experimental measurements of population decay of hyperfine levels in Pr 3+ : Y 2 SiO 5 . The experiments are done using an alternative method to measure rates of individual transitions using hole-burning spectra. An earlier work [29] used hole-burning spectra in Pr 3+ :YAlO 3 to fit 21 parameters to theoretical rate equations. We reduce the number of fitting parameters to 3 by initializing the ground state population in one of the hyperfine levels in a narrow spectral region. Three additional parameters are used to describe the effect of small magnetic fields between 5-10 mT on the flip-flop rates on each of the transitions.
The paper is structured as follows: We first introduce the relaxation pathways for flip-flop interactions considered and enumerate the steps in simulating a distribution of flip-flop rates in Section II. We then explain the new experimental method used to measure population decay in Pr 3+ : Y 2 SiO 5 in Section III. In Section IV, we compare the experiments with simulations to extract the flip-flop rates and also show that the distribution of rates arises from a distribution of ion-ion distance. Lastly, we conclude with some comments on further additions to the microscopic model.

II. MICROSCOPIC MODEL FOR FLIP-FLOP INTERACTION
In this section, we first explain the relaxation pathways and different strengths of each pathway considered in the model, with the example of Pr 3+ :Y 2 SiO 5 . We set up the magnetic dipole-dipole interaction Hamiltonian for a pair of ions and explain the use of Fermi's rule to calculate the flip-flop rate between them. Then, we enumerate the steps in simulating flip-flop interactions and specify the parameters used for Pr 3+ :Y 2 SiO 5 .
Figure 1(i) shows the three hyperfine levels in the electronic ground state 3 H 4 Pr 3+ :Y 2 SiO 5 and the twelve relaxation pathways considered in the simulations. Each hyperfine level is doubly degenerate but the degeneracy is lost in the presence of an external magnetic field, giving rise to six levels in total. It should be noted that conventional labels for the hyperfine levels are ± 1 2 g , ± 3 2 g , ± 5 2 g but each level is in reality an admixture of all six hyperfine wave functions. So we instead use the labels ± |a , ± |b , ± |c . For six levels, we could expect fifteen unique flip-flop transitions. But we do not consider the transitions where only the parity changes, for example transitions of the type + |a ↔ − |a since we do not measure these individually in our experiments. Hence, we have twelve pathways in total. In the experiments described in Section III, it is seen that the strongest interaction is between ions occupying ± |a and ± |b (indicated by solid double-sided arrows), while ± |c couples weakly to the other two levels (indicated by dashed and dotted double-sided arrows). Figure 1(ii) shows the dependence of strength of flip-flop interaction on the hyperfine level occupied by ions. For example, ions initialized in ± |a and ± |b flip-flop strongly to give fairly mixed populations (blue and red circles) while ions in ± |c (purple circles) flip-flop with either of the other two levels with less likelihood. Figure 1(iii) visualizes how the interaction strength scales with distance as r −6 , thus closely-lying neighbours in a crystal interact strongly (shown as solid ovals) while ions far away from each other show weaker interaction (shown as dashed and dotted ovals).
The rate for ion 'i' to flip from |x to |y due to interactions with its 'j' neighbours initially in the state |y is calculated using Fermi's Golden Rule : It is worth noting that calculation of the matrix elements | y i ⊗ x j H ij dd x i ⊗ y j | for all pair of ions makes our model 'microscopic', setting it apart from previous similar works, where this was taken as an average value and related to the concentration of dopants in Ref. [26,27].
Wavefunctions of hyperfine levels x i and y j are the eigenstates of the spin Hamiltonian and they depend on the external magnetic field B. They are calculated using the following equation for the spin Hamiltonian, as used in [30] : The crystallographic axes of the crystal D 1 D 2 b form the common frame of reference for the above calculations.Ĩ is the vector of nuclear spin operatorsĨ x ,Ĩ y ,Ĩ z and B is the magnetic field vector. M is the effective Zeeman tensor and Q is the effective quadrupole tensor, defined as follows : Each of the above matrices is transformed into the frame D 1 D 2 b using rotation matrices with appropriate Euler angles : R k = R(α, β, γ). The two terms on the right-hand side of Equation (2) are evaluated according to [22]: B.M.Ĩ = g pq B pĨq andĨ.Q.Ĩ = Q pqĨpĨq where p, q = x, y, z and the usual summation rules are to be observed whenever a suffix occurs twice.
H ij dd is the Hamiltonian for magnetic dipole -dipole interaction between an ion 'i' and a neighbouring ion 'j' [22].
where each of the suffixes p, q, s, t take the values x, y, z. r ij is the vector connecting the two ions.
The last factor in Fermi's golden rule in Equation (1) is f (E), the density of initial and final states for transitions between two levels in the continuum of initial and final states x i ⊗ y j and y i ⊗ x j respectively. The form of density of states we use is f (E) = 1 πh where Γ hom and Γ xy are the homogeneous and inhomogeneous linewidths of the transition |x ↔ |y . Γ hom is a function of external magnetic field B and κ xy (B) is a phenomenological addition to describe the increase in inhomogeneous linewidths in the presence of a magnetic field. The details of derivation of density of states is given in Appendix VI A. In principle, all pairs of ions are spectrally separated by a different value in the distribution of spin inhomogeneous broadening Γ xy . However, here we focus on the microscopic effect of distances between the ions being different and take an average value for Γ xy .
In brief, the simulation steps required for calculating flip-flop rates are : 1. A small sphere of a host crystal is simulated, where ions are placed according to the crystal lattice structure [31]. It is doped with a rare-earth ion with the specified concentration. Alternatively, one could also assume a continuous random distribution function of ions to determine the position of nearest N th neighbour, as done in [32]. More details about modelling the host crystal can be found in [33]. An ion 'i' is picked in the sphere and nearest neighbours 'j' are found. Nuclear wave functions + a i .. − c i and + a j .. − c j are calculated to be eigenstates of the spin Hamiltonian in Equation (2) and depend on the orientation of the ion in the crystal and the magnetic field.
2. The dipole-dipole interaction Hamiltonian for ion 'i' due to interaction with neighbours 'j' is calculated according to Equation (5).
3. Flip-flop rates for the transitions between all hyperfine levels are calculated using Fermi's rule in Equation (1).
In Pr 3+ : Y 2 SiO 5 , only the ions in site 1 corresponding to the 3 H 4 → 1 D 2 transition at 606 nm were used. The radius of sphere used was 100 nm and the flip-flop rate of ion 'i' was calculated due to the interaction with its twenty nearest neighbors. Pr 3+ has a nuclear spin 5 2 , thusĨ is a (3 × 1) vector where each element is a (6 × 6) matrix. The eight basic molecules in a unit cell of Y 2 SiO 5 have four different directions so for any ion 'i', the tensors M and Q in Equation [2] have one of the four orientations. Values for all the parameters in Equation (3) and (4) were taken from Raman Heterodyne Spectroscopy measurements done in Ref. [30]. The magnetic field was directed along the crystal axis b. Homogeneous linewidths Γ hom = 1 πT2 were taken Ref. [34,35] where the spin coherence time T 2 was measured to be 0.5 ms with zero magnetic field and 6 ms in the presence of magnetic field of 2mT. It does not change appreciably even up to 100 mT, so 6 ms was used for the data with a field between 5-10 mT. The values for all the individual transitions have not been measured, so the same was used for all. Furthermore, our experiments do not distinguish between the rates of the form of + |a → + |b from + |a → − |b , − |a → + |b or − |a → − |b , so we in the following sections sum and average the rates such that only three effective rates R ab , R bc and R ac were obtained for each ion 'i'. The details of reducing twelve rates down to three are described in Appendix VI B. After this reduction, the model contains six unknowns: the three inhomogeneous spin linewidths Γ ab , Γ bc , Γ ac and the factors describing their magnetic field dependence κ ab , κ bc , κ ac used in the density of states f(E).

III. EXPERIMENTS
Relaxation between spin levels has been studied in many different ways, for example using methods that combine optical spectral hole burning and RF fields resonant with a hyperfine transition [36][37][38]. A method to extract rate constants for individual transitions using only hole burning spectra has been used in [29] but it requires many fitting parameters for each rate equation to be able to keep track of the initial population of any ion that was excited during the hole burning. For example, Pr 3+ :YAlO 3 has three hyperfine levels in the ground and excited states. Thus, a laser at a single frequency on the 3 H 4 → 1 D 2 transition can excite nine different transitions or classes of ions. So the method in Ref. [29] required 21 independent fitting parameters (18 initial spin populations and 3 rates). Here, we present an alternative method to measure individual transition rates by initializing population in one hyperfine level (or, equivalently in a single class) within a narrow spectral region and tracking the decay of this state-specific hyperfine population versus its neighbouring spectral background. This method can be advantageous for measurements in rareearth ions with more than one ground hyperfine level, where there are multiple classes of ions since the number of parameters for initial spin population are reduced due to initialization.
We now describe the steps in experiments. We first create a transmission window using spectral tailoring techniques as described in Ref. [39] and initialize the population in one of the ground state hyperfine levels within a spectral region of 1 MHz inside the window. This enables coupling of the laser to a single class of ions and appropriate selection of spectral background range enables us to monitor only this class of ions rather than all the nine classes. Initial population conditions and evolution for all classes of ions are explained in Appendix [VI C]. By probing the ions at different intervals of time, we recorded decay curves for each of the levels, up to 2700s. The absorption structure was erased and the population was reset using a strong frequency scanning pulse after the last readout. The transmission window was then recreated. Experiments were also carried out in the presence of an external magnetic field in the range 5-10 mT, along the crystal axis 'b'. For each experiment, the field was turned on after the step of population initialization. For a given hyperfine level, the population decays at the same rate (within ± 5-10 %) in the range 5-10 mT. Thus, we take the average of the decay for each hyperfine level for this range of magnetic fields.
Transitions used for evaluating the population and an example of the absorption structure after the initialization process within 1 MHz region are shown in Figure 2. The optical transitions labelled as '1', '2', '3' in (i) have higher oscillator strength than other transitions, thus the corresponding absorption peaks in (ii) show high absorption and are used for data analysis. Three spectra in blue, red and purple (inset) show the absorption spectrum after initializing ions in ± |a ,± |b and ± |c respectively. Evaluation of population is done in two steps. First, a slope is subtracted across the width of each peak since the background on either side might be different on the low and high frequency sides of the peak. This can be seen, for example in peak '1' in Figure 2 under the peaks labelled '1' is summed up to obtain the population ± |a and the same is done for peaks '2' and '3' to obtain populations in ± |b and ± |c respectively. Background absorption level is indicated with black arrows at 2 MHz, 12.2 MHz and 38.9 MHz in Figure 2. More details can be found in Appendix VI E.
All experiments were done in a Pr 3+ :Y 2 SiO 5 crystal with 0.05% concentration and dimensions 10mm x 10mm x 0.8mm along D 1 ,D 2 ,b axes respectively. The crystal was placed inside a liquid helium bath cryostat and cooled down to ∼2 K. The light source was a dye laser tuned to the 3 H 4 → 1 D 2 transition in Pr 3+ :Y 2 SiO 5 at 606 nm and was locked to an ultra low expansion glass cavity using the Pound-Drever-Hall locking technique, reducing the linewidth to sub-kHz. A schematic of the experimental setup is shown in Figure 3. All the pulses were shaped using an arbitrary waveform generator and two AOMs. A half-wave plate in combination with a polar- izer aligns the polarization of the light to the D 2 axis of the crystal with an absorption coefficient measured to be 40 cm −1 . The optical power of light for burning pulses was about 20 mW. The readout probe had sufficiently low power such that the same absorption structure could be read up to 100 times without disturbing the population. This was checked by reducing the power until the change in absorption after 100 readouts was within shotto-shot fluctuations. A collimated 1 mm diameter beam, propagating along the b axis (0.8 mm) of the crystal was used. More details about experiments are described in Appendix VI D.

IV. RESULTS AND DISCUSSION
Results of population decay are shown in Figure 4. We will first describe the relaxation dynamics in Figure  4(i-ii), in the absence of an applied field. The figures show the experimental data with error bars indicating the weighted standard deviation of three data sets taken for decay after initializing the populations in each hyperfine level ± |a (blue), ± |b (red) and ± |c (purple). There was no external magnetic field applied but previous measurements indicate that there is a residual field <0.2mT in our cryostat. Population decay for the first 5 ms is negligible so a moving average is performed up to this point and subsequently, the population is normalized with respect to this point. Decay from ± |c is slower than ± |a or ± |b so one can expect R ac and R bc to be lower than R ab . In other words, ions in ± |c flip-   flop with those in either of the other levels at a much slower rate. As described earlier in Section I, a single average value for each of R ab , R bc , R ac is typically used to describe the relaxation of all ions in the crystal. Thus, all ions relax bi-exponentially with N = 1 in Equations (12), (13), (14) in Appendix VI C. As an example of this macroscopic' model, we attempted a biexponential fit to our data, plotted using dashed-dotted lines in Figure 4(i). The best fit obtained for ± |a , ± |b , ± |c respectively was 0.52e −t/5.52 + 0.48e −t/2193 (blue), 0.48e −t/5.52 +0.52e −t/2193 (red) and 0.03e −t/5.52 + 0.97e −t/2193 (purple), where time t is in seconds. While these curves fit well to many data points, several data points do not follow the fits especially ± |a and ± |b .
Each decay curve obtained in Figure 4 is in fact, an average of many exponential decays of different ions within the 1 MHz peak shown in Figure 2(ii). Each ion may have a different flip-flop rate for a given transition, depending on its position and orientation in the crystal. In the microscopic model, the effective decay is instead an average of the bi-exponential decay of many ions in the crystal, shown as the solid colored (blue, red, purple) lines in Figure 4(ii). These are the simulations which evaluate population according to steps detailed in Section II and they match the experimental data quite well. The solid colored lines in Figure 4(ii) and (iii) show the fits from simulation of our microscopic model in the absence and presence of magnetic field respectively. In addition to the list of steps in simulations described in Section II, a few more steps were followed in order to be able to compare the simulations with the experiments: 1. Using the rates R ab , R bc , R ac , population decay in levels ± |a ,± |b ,± |c is calculated using respectively, Equations (12), (13), (14) derived in Ap-pendix VI C.

2.
Steps [1][2][3] from the list in Section II are repeated for i = 2, 3...N , where N is the number of ions in the sphere. An average decay of 'N ' ions gives a single decay curve describing the decay of all ions in the crystal. These are the solid colored lines in Figure 4(iii).
3. All of the above steps are repeated for data with an external magnetic field.
The experiments show little difference between the decay from ± |a and ± |b , indicating that ions occupying these states have the strongest magnetic dipole-dipole interaction. This is shown by the blue and red solid lines almost overlapping with each other in 4(i). The optimized values of spin inhomogeneous linewidths, Γ ab , Γ bc and Γ ac were found to be 0.618, 3.309 and 2.664 kHz respectively. We choose ions in a sphere of radius 100 nm for the simulations, thus the fitted linewidths represent the local spin inhomogeneity and can be less than the measured values Γ ab = 50.5 kHz, Γ bc =75.4 kHz [30] in a bulk crystal. The optimization is fairly insensitive to Γ ac and the relaxation is predominantly governed by the rates R ab and R bc . We now try to understand why the fitted values of Γ ab and Γ bc differ by a factor of ∼4.4. A possible contribution to spin inhomogeneity is inhomogeneity in the g-tensor which stem from strains or defects [40]. Local inhomogeneity in spin could also be due to magnetic dipole-dipole interactions between a Pr ion with its neighbouring Pr ions of the type given by Equation 5. If the hyperfine wavefunctions ± |a , ± |b and ± |c were composed of pure ± 1 2 g , ± 3 2 g and ± 5 2 g states, then the shift in hyperfine frequencies due to interaction between a pair of Pr ions scales linearly with the quantum number m I = 1 2 , 3 2 , 5 2 . Thus, the frequency shift of one ion in a pair occupying ± 5 2 g and ± 3 2 g is five times larger than that of another ion in a pair occupying ± 1 2 g and ± 3 2 g . Furthermore, the effect of external magnetic field is largest on ± |c since it undergoes a larger Zeeman shift compared to ± |a , as seen in Figure 8 (ii) and (vi) in Appendix VI E. This could explain why Γ bc is ∼4.4 x Γ ab even though ± |a , ± |b and ± |c are actually an admixture of the pure hyperfine states.
For the experiments with magnetic field in Figure  4(iii), the field is put on after the population initialization step and it takes a few seconds for the field to ramp up to the set value. Thus, the simulation evolves the population until the dotted line at 4.6 seconds assuming there is no external field and normalizes the data so that the population at the time corresponding to the dotted line in (iii) equals the population in (ii). Some of the data points in ± |c in (iii) before the dotted line show population greater than 1. This is an experimental artefact and occurs because the peak corresponding to these points, '3' in Figure [2] split in the presence of field due to nuclear Zeeman effect and thus a different spectral region is chosen for evaluating population before and after the peak has split. This is shown in detail in Appendix [VI E].
After the dotted line, it is assumed that the magnetic field has reached the set value and the simulation evolves the population by including the phenomenological terms, κ ab , κ bc and κ ac introduced earlier in Fermi's rule [1]. The optimized values were found to be respectively, 2.6, 3.6 and 1.5 with a field between 5-10 mT. Figure 5 shows the effect of magnetic field on the calculated rates, where (i),(ii) and (iii) show the histogram of R ab , R bc and R ac respectively with and without a magnetic field. R ab (with no applied field) is spread over a distribution ranging from 10 −4 − 10 2 Hz and peaks at 10 −2 Hz. R bc is slower, ranging from 10 −6 − 1 Hz and peaks at 10 −4 Hz, with no field while R ac is slowest, ranging from 10 −7 − 1 Hz and peaks at 2 x 10 −6 Hz. All three rates slow down by two orders of magnitude with a field of 5 -10 mT. The distribution of rates shown in (i) -(iii) follow from the distribution of r −6 N N shown in (iv) and the inset shows the distribution of r N N , where r is the distance to any of the twenty closest neighbours of any ion considered in the simulations.
To understand why the rates slow down in a magnetic field, one can infer from Equation (1) that the cause could either be evolution of matrix elements in the dipoledipole interaction term or a change in density of states f (E). While the matrix elements do not change appreciably with a small field of 5 -10 mT, the density of states changes drastically due to the decrease in homogeneous linewidth by more than a factor of ten, as measured in Ref. [35] and this is attributed to minimizing spin flips of the neighbouring Y ions. In the absence of an external field, the magnetic field experienced by the core Y ions is due to the local Pr ion, which is of the order of ∼ 0.1 mT and, a change in the spin state of Pr flips the spin state of Y ions. Thus, dephasing of Pr ions is dominated by neighbouring Y flips in the core. When the external field significantly exceeds the field due to the local Pr ion, such flips are minimized. Another factor contributing to the change in density of states is the increase in the spin inhomogeneous linewidth, characterized by the fitting parameters κ ab , κ bc and κ ac . A linear increase in spin inhomogeneous linewidths has also been reported in Nd 3+ :Y 2 SiO 5 [27] and in erbium doped glass fibers [41]. Measurements of spin linewidths as a function of magnetic field has partly been done in some Kramers ions [42] and similar measurements in Pr 3+ :Y 2 SiO 5 may shed more light on this explanation but such data is unavailable at this point.
We conclude this section by noting that there are two conditions that need to be satisfied for two Pr ions to flip-flop : they need to be close to each other in the crystal and they also need to be spectrally close in the spin inhomogeneous profile. In our model, we take an average value for the spin inhomogeneity and model the ion-ion distance as a distribution. One could also model the spin inhomogeneity as a distribution, for example by including the effect of the local magnetic field around each Pr ion. The term B in Equation (2) could be replaced by B total = B ext + B local so that each ion has a unique Spin Hamiltonian, resulting in a distribution of Zeeman frequencies of Pr ion.

V. CONCLUSION
We have presented a method to model microscopic effects of flip-flop interactions between individual ions in a rare-earth-ion doped crystal. We have simulated a random doping based on the crystal structure of the host, where the position and orientation of all ions is known. Every dopant ion is situated in a unique position and orientation with respect to its neighbours so the ion-ion distance is a distribution and the flip-flop rate of any ion with its neighbours is different owing to this distribution. We apply this model to experiments of population decay of ground state hyperfine levels in Pr 3+ :Y 2 SiO 5 . The experimental method used is an alternative to methods used in earlier works. The collective relaxation dynamics of all ions probed in the crystal is an average sum of many exponential decays of different ions. Thus, the flip-flop rate between two hyperfine levels is a distribution of rates rather than one average rate describing the dynamics of all ions.
The fastest rate is R ab between the levels ± |a and ± |b , whose distribution has a peak at 10 −2 Hz while R bc and R ac have a peak at 10 −4 and 2 x 10 −6 Hz respectively, in the presence of a residual field <0.2 mT. All the rates decrease by 2 orders of magnitude upon applying an external field of 5 -10 mT and the reason could be a combination of an order of magnitude decrease in the spin homogeneous linewidths and an increase in spin inhomogeneous linewidths [26,27]. An improvement to the model could be to include the effect of differences in the local magnetic field around each dopant ion. Nonetheless, our model serves as a general tool to calculate other kinds of interactions at the microscopic level. It could be used to study the dynamics of other rare-earth ions in different materials as well.

A. Density of states for transitions between levels in two different continuum of states
In this section, we give details of how Fermi's rule is applied to the case of a flip-flop transition between two homogeneously broadened levels centered around different energies. We first start with Fermi's rule for a transition between two discrete levels and then extend it to the case of transition between a level in a continuum of states with finite width to another continuum of states. Finally, we apply this to flip-flop transitions, where each state is a two-level system.
Let us start with two levels |1 and |2 with energies E 1 and E 2 respectively. The transition rate from a discrete state with energy E 1 in |1 to E 2 in |2 under a perturbation H is given by Fermi's rule : Assuming the matrix element | 1|H |2 | is independent of the energies E 1 and E 2 , the total transition rate from the continuum of states in |1 → |2 is given by integrating R E1→E2 over the density of states for both the initial and final energies : Further, assume that both ρ 1 and ρ 2 have a normalized Lorentzian lineshape centered around 1 and 2 with homogeneous HWHM (half width at half maxima) ∆ 1 and ∆ 2 respectively such that, ρ l (E) = 1 π ∆ l ∆ 2 l +(E− l ) 2 , where l = 1, 2. The rate is then : For the case of flip-flop transitions, perturbation H is the magnetic dipole-dipole Hamiltonian between two ions 'i' and 'j' : H ij dd . |1 is the two-level system x i ⊗ y j , which flips to |2 i.e., y i ⊗ x j due to H ij dd , where |x , |y are the wavefunctions of the hyperfine levels + |a , − |a ... − |c (shown in Figure 2). In the equation 7 above, 1 and 2 are the energies of x i ⊗ y j and y i ⊗ x i respectively and ∆ 1 and ∆ 2 are their respective homogeneous HWHM of Lorentzian lineshapes. The lineshape of such a transition is given by the convolution of the two Lorentzian lineshapes of the individual levels and is another Lorentzian with HWHM ∆ 1 +∆ 2 , centered at 1 + 2 . From Ref. [34,35], we take 2∆ 1 = 1 πhT2 , where T 2 is the spin coherence time at zero magnetic field and it is equal to 0.5 ms (6 ms with a field of 5-10 mT). Considering ∆ 1 = ∆ 2 , HWHM is ∆ = 1 πhT2 , which is just the homogeneous linewidth of the transition |x ↔ |y . Also, ( 1 − 2 ) = hΓ xy , where Γ xy is the spin inhomogeneity .
Writing Equation 7 in terms of frequencies and replacing ∆ 1 + ∆ 2 by hΓ hom and ( 1 − 2 ) by hΓ xy , we get the final expression : Thus we can use the last term on the right of the above equation as the form of density of states f (E). Similar expressions have been used to describe transitions between broadened states in different quantum wells, as described in Section 3.3 of Ref. [43]. In Equation 1 in the main text, the homogeneous linewidths are functions of magnetic field and we also have a phenomenological factor κ xy to describe the increase in inhomogeneous linewidths in the presence of magnetic field.

B. Reduction of rates from 15 to 3 and optimization of parameters in simulations
In the absence of an external field, there are three ground state hyperfine levels in Pr 3+ :Y 2 SiO 5 and each splits into two in the presence of a field as shown in Figure 1. There can be fifteen unique rates due to magnetic dipole transitions in this case. Flip-flop interactions where initial and final state are the same or only change parity, for example R +|b ↔−|b are ignored since we do not measure these individually in our experiments. Thus, the simulations calculate twelve unique rates. However, our experiments are designed to measure only three rates R ±|a ↔±|b , R ±|b ↔±|c , R ±|a ↔±|c , or referred to in the main article as R ab , R bc , R ac . Each of the rates is divided by 6 since the neighbouring ion can only be in one of the six hyperfine levels. So we are left with the task of reducing the twelve rates from simulations down to three.
We divide the rates into three categories, each involving the pair of levels ± |a and ± |b , ± |b and ± |c , ± |a and ± |c . Histogram of rates involving the transitions in each pair are shown in Figure 6 (i),(iii),(v). We first consider (i).
• Transitions originating from the same level are added together. For example, The rates + |a → + |b and + |a → − |b (labelled as (I) and (II) in the Figure 6 (i) ). Similarly (III) and (IV), − |a → + |b and − |a → − |b are added together.
• Since we cannot distinguish + |a → + |b from − |a → + |b , we take the average of (I+II) and (III+IV), which is reasonable since they are very similar and in experiments, the split peaks appear to decay at the same rate (see Figure 8(iv) and (vi) for the case of ± |b and ± |c . This total rate is called as R ab . The relative difference between the two quantities, calculated as [(I+II)-(III+IV)]/Mean[R ab ] is shown in Figure  6(ii).
A similar argument is applied to (iii-iv) and (v-vi) in Figure 6, thus the rates are reduced from twelve to three.
Optimization of the parameters Γ ab , Γ bc , Γ ac , κ ab , κ bc and κ ac was done using Global Search and fminsearch functions in Matlab. The cost function, which is a measure of deviation in population between simulated model (p m ) and experiments (p e ) is minimized to have the lowest 'score' simultaneously for decay in all three hyperfine levels as well as for all values of magnetic field : For data without external magnetic field, w is the standard deviation of three data sets taken at each time point and for cases with magnetic field, w is the standard de-viation of three data sets, each taken at 5, 7 and 10 mT. In total, there are six parameters describing the decay : Γ ab , Γ bc , Γ ac , κ ab , κ bc , κ ac . The score is most sensitive to value of Γ ab . A change of ±5% in Γ ab changes the score by ∼ 10% while the same change in Γ bc or κ ab changes the score by ∼ 3%. Changing κ bc by ±5% results in a change in score by ∼ 1.5%. Thus, the relaxation is mostly governed by the fastest rate, Γ ab . It was also checked that increasing the size of the sphere or number of neighbours did not further improve the score appreciably, thus it is sufficient to calculate the effect of twenty nearest neighbours on each other and in a YSO crystal, this distance varies between 1-20 nm.

C. Rate Equations and Initial Conditions of population
In this section, we derive the rate equations for a three level system ± |a , ± |b , ± |c shown in Figure 2 and elucidate the initial conditions of population for nine classes of ions. If N a , N b , N c are populations normalized to the total population N and R ab , R bc , R ac are the flip-flop rates between the three hyperfine levels, the rate equations can be written as: f peak = 0 MHz f background = 2MHz Class Transition probed at 0 MHz n a n b n c n a n b n c I ± |a → ± |a e 1 0 The solutions for arbitrary initial conditions N a = n a , N b = n b , N c = n c are : N c (t) = n a + n b + n c 3 + 1 6σ e −t(K+σ) n a A 2 where Due to the optical inhomogenous broadening, the laser can couple to nine different transitions from the ground to excited state (shown in Figure 2) at a given frequency.
To extract the individual rates R ab , R bc , R ac from a simple hole burning spectra, one would need to follow the evolution of all classes [29] by summing over contributions from nine transitions for each of the three levels with three unknown initial conditions, thereby giving 30 unknowns. We simplify this by creating a transmission window using optical pumping and isolating peaks of ions in one hyperfine level within this window. A simulation in a six level system (with three ground state and three excited states) was performed to predict an absorption spectra after the initialization step.
An example of initializing in ± |a , with a peak at 0 MHz corresponding to ± |a → ± |a e is shown in Figure  7. It can be seen that the experimental data in purple matches quite well with the black line showing the simulation, indicating successful isolation of one class of ions from the other eight classes. We call this isolated group of ions as Class I, for which the starting conditions after the initialization are n a = 1, n b = 0, n c = 0 and the corresponding spectral background region at 2MHz has all Class I ions shelved in n c . The initial population conditions of the other eight classes for both the peak and background are charted out in Table [I]. The conditions for all classes but one, 'I' are the same for both peak and background. Thus, we can subtract background from the peak for this class only using Equation (9) to calculate population decay in ± |a . This describes the decay of an ion 'i'. To account for all ions i = 1, 2..., N in the sphere considered in the simulations, we take the average as follows: Similar considerations for peaks corresponding to ± |b → ± |b e at 14.8 MHz and ± |c → ± |c e at 36.9 MHz (peaks '2' and '3' respectively in Figure 2 are given in Tables II and III. Coupled with simulations similar to Figure 7, isolation of one class was ensured. The evolution of population N b (t) and N c (t) can also be obtained:

E. Absorption spectra with and without magnetic field
Absorption spectra from the first and last readout are shown in Figure 8. (i),(iii) and (v) show the absorption after initializing the populations in ± |a (blue), ± |b (orange) and ± |c (purple) respectively. All black traces show the absorption after the last readout at t f = 2700s. Although no external magnetic field was applied in these cases, we expect there to be a stray field less than 0.2mT. Results of similar experiments with an external field of 10mT are shown in (ii),(iv),(vi) . The first readout t 0 is at 3.8s for (ii) and (iv), 5.5s for (vi), which is also indicated by the vertical dashed line in Figure 4. Peaks '2' and '3' in (iv) and (vi) respectively split due to nuclear Zeeman effect. Although the absorption level of split peaks is roughly halved, the absorption level of their backgrounds are not the same as their counterparts in (iii) and (v). Due to this, some of the data points before the dotted line in Figure 4(iii) are greater than 1. A different choice of background could perhaps have been better but the important information for simulations is what happens to the population after the magnetic field has reached the set value at the dotted line.