# Decoherence induced by dipole-dipole couplings between atomic species in rare-earth ion-doped Y$_{2}$SiO$_{5}$

Charlotte Pignol, Antonio Ortu, Louis Nicolas, Virginia D'Auria, Sébastien

Tanzilli, Thierry Chanelière, Mikael Afzelius, Jean Etesse

► To cite this version:

Charlotte Pignol, Antonio Ortu, Louis Nicolas, Virginia D'Auria, Sébastien Tanzilli, et al.. Decoherence induced by dipole-dipole couplings between atomic species in rare-earth ion-doped Y$_{2}$SiO$_{5}$. Physical Review B, 2024, 110 (21), pp.214208. 10.1103/PhysRevB.110.214208. hal-04668718

HAL Id: hal-04668718

https://hal.science/hal-04668718v1

Submitted on 7 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Decoherence induced by dipole-dipole couplings between atomic species in rare-earth ion-doped Y₂SiO₅.

C. Pignol,¹ A. Ortu,²,³ L. Nicolas,⁴ V. D'Auria,¹ S. Tanzilli,¹ T. Chanelière,⁵ M. Afzelius,⁴ and J. Etesse¹

¹Université Côte d'Azur, CNRS, Institut de Physique de Nice, 06200 Nice, France

²Institut für Quantenoptik und Quanteninformation,
Österreichische Akademie der Wissenschaften, Technikerstrasse 21a, 6020 Innsbruck, Austria

³Universität Innsbruck, Institut für Experimentalphysik, Technikerstrasse 25, 6020 Innsbruck, Austria

⁴Department of Applied Physics, University of Geneva, CH-1211 Geneva, Switzerland

⁵Université Grenoble Alpes, CNRS, Grenoble INP, Institut Néel, 38000, Grenoble, France

(Dated: August 6, 2024)

Rare-earth ion doped crystals are state-of-the-art platforms for processing quantum information, particularly thanks to their excellent optical and spin coherence properties at cryogenic temperatures. Experimental observations have shown that the application of a static magnetic bias field significantly improves the coherence times in the rare-earth ions ensemble, but only a few studies have focused on its the dependency as a function of both magnetic field direction and amplitude. This is especially true for magnetic field amplitudes under the mT, and for low magnetic dipole moment ions. In this paper, we investigate the relationship between the magnetic field parameters and the decoherence caused by magnetic dipole-dipole coupling with the nearest neighbors nuclear spins in the crystal. The primary non-Kramers rare-earth ions investigated here are europium and praseodymium, but we also extend our study to the ytterbium Kramers ion due to its low magnetic dipole in the mT range. We perform theoretical investigations and simulations of the energy structure and coherence time evolution and identify good correspondences between experimental and simulated spin echo data. This work allows us to pinpoint the most relevant decoherence mechanisms in the considered magnetic field regime, and to predict favorable magnetic configurations.

## I. INTRODUCTION

Quantum networks offer exciting prospects such as unconditionally secure communication through quantum key distribution [1, 2], enhanced computational power via distributed quantum computing [3, 4], and ultra-sensitive measurements with quantum sensors [5, 6]. Furthermore, they hold promise for fundamental scientific investigations, including testing the foundations of quantum theory [7, 8]. The foreseen challenges for creating operational quantum networks include (i) maintaining quantum states over distances above hundreds of kilometers despite propagation losses in fibers [9], (ii) managing decoherence and noise [10], and (iii) developing reliable protocols for quantum operations and information transfer [11].

This is where quantum memories set themselves as crucial elements. They can act as quantum repeaters to extend the distance over which quantum information can be propagated [12–14], to delay the information transmission for protocols that require synchronization [15], as well as to enable single photon generation [16] and perform operations on qubits [17]. Many platforms have been envisioned to implement quantum memories, such as cold atoms [18], hot atomic vapors [19, 20], solid-state artificial atoms [21], and all-optical storage loops [22]. Amongst all these possibilities, rare-earth (RE) ion doped crystals are recognized as state-of-the-art candidates. Moreover, such material also allow, among others, for photonic information processing [23, 24], single-spin manipulation for microwave processing [25] as well as opto-microwave transduction [26].

Their principal interest lies in their intrinsically long optical and spin coherence times [27], allowing to perform storage of up to an hour [28]. Furthermore, their large optical inhomogeneous broadening offers great multiplexing capabilities which allow handling thousands of modes in the spectro-temporal domain [29]. One means to achieve long spin coherence time is to operate in a very specific magnetic regime called Zero First-Order Zeeman (ZEFOZ) point (also referred to as clock transitions), where the magnetic sensitivity to the environmental fluctuations is minimal [30]. The downsides of working in this regime are that it usually requires to operate in the 100 mT-10 T amplitude range, and to adjust the field direction with very high precision [31]. In contrast, recent experiments have shown that the application of a weak magnetic field (< 10 mT) already suffices to strongly affect the spin coherence times of the ensembles [32, 33]. Depending on the RE species, such effect can lead either to an increase [32–34] or to a decrease [35] of the spin coherence time.

Even if each RE ion ensemble behaves differently in different magnetic field regimes, we however postulate that the underlying mechanisms that lead to spin decoherence are the same for all species and stem from magnetic dipole-dipole coupling [36]. Several approaches have been used to explore this question, using either averaging techniques, perturbative treatments, and stochastic models, that enabled the models to be effective within a specific frame of study [37–40]. More recently, a microscopic model based on dipole-dipole interaction in erbium allowed to replicate complex decay patterns [41] and selective spin bath ion addressing [42]. This ap-
proach, directly linked to cluster correlation expansion (CCE) methods [39, 43, 44], can however not be used with all ions in all magnetic field regimes. Inspired by this approach, we present here a general and predictive model that can be applied to any ion species in any host matrix, allowing one to estimate the effect of magnetic dipole-dipole interaction on coherence time.

The main novelty of our approach is that it does not rely on particular coupling or dynamic hypothesis but on an exact simulation of spin dynamics in interaction with a few (~5) nearest neighbors nuclear spins. To this end, we adopt a microscopic approach and theoretically investigate the impact of magnetic dipole-dipole couplings between the host matrix ions (that possess a nuclear spin) and the RE ion on the coherence time. In our study, we focus on three RE ions: europium $^{151}$Eu$^{3+}$ (Eu), ytterbium $^{171}$Yb$^{3+}$ (Yb) and praseodymium $^{141}$Pr$^{3+}$ (Pr). Here, we consider the yttrium orthosilicate matrix ($Y_2SiO_5$, hereafter noted YSO), in which the RE ions substitute Y$^{3+}$ ions. Such a crystal is commonly used due to its low nuclear spin density [31, 34, 45]. We adopt three different approaches for the three RE at stake for confrontation with our numerical model. For europium, we perform an experiment based on the atomic frequency comb (AFC) protocol [46] to estimate the spin coherence time. For ytterbium, we use already published experimental data. Finally, for praseodymium, we only perform numerical simulations to predict a favorable magnetic configuration. As stated above, the considered magnetic field amplitude ranges from the µT to the mT.

This article is organized as follows. First, we present the theoretical model of the problem in Sec. II, by introducing the model itself in Sec. II A, and its possible alternatives and approximations in Sec. II B. We then present our simulations and findings in Sec. III. We compare our simulations with experimental results for Eu in Sec. III A, then for Yb in Sec. III B. We continue by performing coherence time predictions for the Pr ion, as a function of magnetic field amplitude and orientation in Sec. III C. In Sec. IV A, we summarize and explicit the relationship between the RE-Y dipole-dipole interaction, the magnetic field and the coherence time. Finally, we present in Sec. IV B the directions for future investigations.

## II. MODEL

### A. Simulated system and Hamiltonian

We show in this article that the decoherence mechanisms in the RE ion-doped crystal under study can be investigated by focusing on a Hilbert space of a rather small dimension restricted to the central RE ion plus a few nearest neighbors nuclear spins in the bath. For that, the concentration in doping RE ions has to be low enough to be able to neglect the magnetic dipole-dipole coupling between themselves. Because of their strong electronic

moment, Kramers ions need a much lower concentration to satisfy this condition than non-Kramers ions. We can cite as an example erbium-doped crystals at 10 ppm [47], meanwhile a 1000 ppm concentration is viable for an europium-doped crystal [48].

Without interactions amongst the RE ions, the crystal is then considered as constituted of randomly distributed independent and identical cells containing one RE ion surrounded by yttrium (Y) ions. We simulate one of these cells, but because of computation power limitations we restrict the number of Y ions inside to a maximum of seven Y ions. The Y ions are added to the system from the closest to the RE ion to the furthest, and their positions are taken from published computations made from the crystal lattice parameters [41].

YSO has a monoclinic structure and belongs to the $C_{2h}^6$ space group. Crystals are usually cut along the ($D_1, D_2, b$) axes that correspond to the polarization extinction axes [49]. The RE ions can substitute the Y ions in two different optical sites, I and II [50]. Each optical site consists of two magnetic sub-sites that are related to each other via the $C_2$ symmetry axis (which coincides with the b axis). For a magnetic field direction in the ($D_1, D_2$) plane or along b, these magnetic sub-sites become equivalent. Most of the simulations presented in the following involve fields in this plane.

#### 1. The rare-earth ion

The RE ion Hamiltonians of non-Kramers (such as Eu and Pr) and Kramers ions (such as Yb) are different. They lead to the energetic structures displayed in Fig. 1. Eu- and Pr- doped YSO have an Hamiltonian of the form [45, 51]:

$$ \hat{H}_{\text{nK}} = \hat{\mathbf{I}} \cdot \bar{Q} \cdot \hat{\mathbf{I}} + \mathbf{B}_{\text{DC}} \cdot \bar{M} \cdot \hat{\mathbf{I}} + \mathbf{B}_{\text{DC}} \cdot \bar{Z} \cdot \mathbf{B}_{\text{DC}}, \quad (1) $$

with the static magnetic field vector $\mathbf{B}_{\text{DC}}$ expressed in the ($D_1, D_2, b$) basis as:

$$ \mathbf{B}_{\text{DC}} = B_{\text{DC}} \begin{pmatrix} \cos(\varphi) \sin(\theta) \\ \sin(\varphi) \sin(\theta) \\ \cos(\theta) \end{pmatrix}. \quad (2) $$

In Eq. (1), $\hat{\mathbf{I}}$ is the nuclear spin vector, $\bar{Q}$ the quadrupolar tensor, $\bar{M}$ the linear Zeeman tensor and $\bar{Z}$ the quadratic Zeeman tensor. The quadrupolar interaction is present for spins $I \ge 1$, and both Pr and Eu ions have an $I = 5/2$ spin. The Zeeman interaction lifts the degeneracy of the quadrupolar levels by an amount dependent on the magnetic field orientation and amplitude. The order of magnitude of this effect is ~1 to 10 kHz/mT for Eu and ~10 to 100 kHz/mT for Pr, as illustrated in Fig. 1. Then, the quadratic Zeeman tensor has (for instance with Eu) a magnitude on the order of 100 Hz/mT$^2$ [31]. This means that the induced dipole is below 100 Hz/mT in the magnetic field amplitude range considered here, far
<figure><img src="image_2.png" /><figcaption>FIG. 1: Schematic energy levels of Eu:YSO, Pr:YSO and Yb:YSO with the energetic structure of the spin nuclear bath $\hat{H}_Y$ for two Y ions. The green arrows highlight the spin transitions under study. Blue, red, and purple arrows in the top inset indicate spin dynamics mechanisms discussed in the body of the text.</figcaption></figure>

smaller than the ~1 to 10 kHz/mT induced by $\bar{M}$. Consequently, in the investigated magnetic field range (μT-mT), we can reasonably neglect the quadratic Zeeman interaction term. The energy structure created by the $\hat{H}_{nK}$ Hamiltonian is schematically displayed in Fig. 1 left. The optical electronic transition addressed is a $4f \leftrightarrow 4f$ transition ($^3H_4 \leftrightarrow ^1D_2$ for Pr at 605.98 nm and $^5D_0 \leftrightarrow ^7F_0$ for Eu at 580.04 nm), for the crystallographic site I of both systems. The spin transitions under study are the radiofrequency (RF) $|\pm 3/2\rangle_g \leftrightarrow |\pm 1/2\rangle_g$ 10.2 MHz transition for Pr and the $|\pm 5/2\rangle_g \leftrightarrow |\pm 3/2\rangle_g$ 46.2 MHz transition for Eu (green arrows in Fig. 1 left). The Yb Kramers ion Hamiltonian is [52]:

$$ \hat{H}_K = \hat{\mathbf{I}} \cdot \bar{\mathcal{A}} \cdot \hat{\mathbf{S}} + \mu_B \mathbf{B}_{DC} \cdot \bar{\mathcal{g}} \cdot \hat{\mathbf{S}} + \mu_n \cdot g_n \cdot \mathbf{B}_{DC} \cdot \hat{\mathbf{I}}. \quad (3) $$

The Yb electronic 1/2 spin $\hat{\mathbf{S}}$ is coupled to its nuclear 1/2 spin $\hat{\mathbf{I}}$ through the hyperfine interaction tensor $\bar{\mathcal{A}}$. It creates a set of four non-degenerate levels when doping YSO, that are separated by a frequency splitting wide enough to selectively address one transition optically. $\mu_B$ and $\mu_n$ are the Bohr magneton and the nuclear magneton, respectively, while $\bar{\mathcal{g}}$ and $g_n$ are the electronic Zeeman interaction tensor and nuclear Zeeman interaction coefficient, respectively.

The energetic structure of Yb:YSO is schematically displayed on the right of Fig. 1. The electronic optical transition for Yb is the $^2F_{7/2}(0) \leftrightarrow ^2F_{5/2}(0)$ at 978.85 nm for the crystallographic site II. The spin transition that we consider here is the $|2_g\rangle \leftrightarrow |4_g\rangle$ microwave (MW) transition at 2.497 GHz (green arrow on the right

of Fig. 1), whose spin decoherence was investigated recently [35].

## 2. The bath

The RE ion is considered to be surrounded by N yttrium ions, sole ions of the YSO matrix that possess a nuclear spin with non-negligible abundance (we neglect the contribution of $^{29}$Si, only present at 4.7% [53, 54]). The nuclear spin bath has an Hamiltonian of the form:

$$ \hat{H}_Y = \sum_{k=1}^{N} \hat{H}_{Y_k}^Z + \sum_{i \neq j=1}^{N} \hat{H}_{Y_i - Y_j}^{dd}. \quad (4) $$

The first term in this expression is the sum of the Zeeman Hamiltonians of each Y ion labeled k:

$$ \hat{H}_{Y_k}^Z = \gamma_Y \mathbf{B}_{DC} \cdot \hat{\mathbf{I}}_{1/2}^k, \quad (5) $$

where $\gamma_Y = 2$ kH/mT is the magnetic dipole moment of yttrium and $\hat{\mathbf{I}}_{1/2}^k$ is the 1/2 spin vector operator. The second term $\hat{H}_{Y_i-Y_j}^{dd}$ is the dipole-dipole interaction of all the Y ions with each other. The general expression of such an Hamiltonian for two ions 1 and 2 whose positions are linked with the vector $\mathbf{r}_{12}$ is [55]:

$$ \hat{H}_{1-2}^{dd} = -\frac{\mu_0 h}{4\pi r_{12}^3} (3(\hat{\mathbf{m}}_1 \cdot \mathbf{r}_{12}) \otimes (\hat{\mathbf{m}}_2 \cdot \mathbf{r}_{12}) - \hat{\mathbf{m}}_1 \otimes \hat{\mathbf{m}}_2), \quad (6) $$
where $\hat{\mathbf{m}}_1$ and $\hat{\mathbf{m}}_2$ are the magnetic moments of the ions, $h$ the Planck constant, $\mu_0$ the vacuum magnetic permeability and $r_{12} = \|\mathbf{r}_{12}\|$. The Y-Y dipole-dipole interaction lifts the degeneracy of each Y Zeeman level. As they depend on the relative distance and orientation of each pair of ions in the crystalline matrix, they create a manifold of energy levels, whose energy separation is on the order of a few Hz that is different for each level pair. The energetic structure associated with this Hamiltonian is displayed in the center of Fig. 1, for two Y ions. Total spin projection configurations have different energies due to the dipole-dipole coupling, and this ultimately leads to an apparent energetic broadening of the central spin transition, as will be illustrated later. Then, we define

$$ \hat{H}_0 = \hat{H}_{\text{RE}} + \hat{H}_{\text{Y}} \quad (7) $$

the diagonalized Hamiltonian defining the energetic structure on which the dynamical evolution will next be considered. We chose to include the Y-Y dipole-dipole interaction in this diagonalized Hamiltonian, leading to a lift of degeneracy in the bath manifolds as illustrated in Fig. 1. Such a choice is motivated by the fact that we want to focus on the effect of the bath as a whole on the RE ion through the RE-Y dipole-dipole couplings as we will introduce in the next paragraph. Notice that choosing to include Y-Y dipole-dipole interaction in the diagonalized part does not mean that we only consider diagonal interaction terms, as sometimes done in the literature [56].

### 3. Dynamics

Now that we have defined the energetic structure of the system, we focus on its dynamics. To this extent, we add a non-diagonal interaction Hamiltonian $\hat{H}_{\text{int}}$ written in the eigenbasis of $\hat{H}_0$ to define the system total Hamiltonian:

$$ \hat{H} = \hat{H}_0 + \hat{H}_{\text{int}}. \quad (8) $$

This interaction Hamiltonian consists of two contributions:

$$ \hat{H}_{\text{int}} = \hat{H}_{\text{RE}}^{\text{RF}} + \sum_{k=1}^{N} \hat{H}_{\text{RE}-\text{Y}_k}^{\text{dd}}. \quad (9) $$

The first term $\hat{H}^{\text{RF}}$ represents the resonant radiofrequency or microwave coupling via Zeeman interaction, which reads

$$ \hat{H}_{\text{nK}}^{\text{RF}} = \mathbf{B}_{\text{AC}} \cdot \bar{M} \cdot \hat{\mathbf{I}}_{\text{RE}} \quad (10) $$

for Eu and Pr, and

$$ \hat{H}_{\text{K}}^{\text{RF}} = \mu_B \mathbf{B}_{\text{AC}} \cdot \bar{g} \cdot \hat{\mathbf{S}} \quad (11) $$

for Yb. These Hamiltonians drive the spin transitions (green arrows in Fig. 1) during the echo sequences and the spin rephasing protocols, via the application of an oscillating magnetic field $\mathbf{B}_{\text{AC}}$ resonant with the spin transition.

The second interaction term $\sum \hat{H}_{\text{RE}-\text{Y}_k}^{\text{dd}}$ describes the interaction of the central RE ion with the Y ions labeled by $k$. This Hamiltonian is also of the form (6), but now with different species in interaction. This contribution is the one that leads to decoherence in the sequences considered here.

## B. Framework for the simulations

In this section, we precisely define how decoherence is modeled, and discuss approximations that are to be taken in the previously defined Hamiltonians. We also compare them with different approximations usually considered in the literature for similar treatment, and explicit why most of them cannot be taken here.

### 1. Approximations in our model

In order to understand which approximations can be made on our system, it is instructive to first rewrite the general dipole-dipole interaction Hamiltonian (6) as [55]:

$$ \hat{H}_{1-2}^{\text{dd}} = -\frac{\mu_0 h}{4\pi r_{12}^3} (\hat{A} + \hat{B} + \hat{C} + \hat{C}^\dagger + \hat{D} + \hat{D}^\dagger), \quad (12) $$

with:

$$ \begin{aligned} \hat{A} &= (1 - 3 \cos^2(\theta_{12})) \hat{m}_1^z \hat{m}_2^z, \\ \hat{B} &= -\frac{1}{4} (1 - 3 \cos^2(\theta_{12})) (\hat{m}_1^+ \hat{m}_2^- + \hat{m}_1^- \hat{m}_2^+), \\ \hat{C} &= -\frac{3}{2} \sin(\theta_{12}) \cos(\theta_{12}) e^{-i\varphi_{12}} (\hat{m}_1^+ \hat{m}_2^z + \hat{m}_1^- \hat{m}_2^+), \\ \hat{D} &= -\frac{3}{4} \sin^2(\theta_{12}) e^{-2i\varphi_{12}} \hat{m}_1^+ \hat{m}_2^+. \end{aligned} $$

In this expression, $\theta_{12}$ and $\varphi_{12}$ are the polar angles of the $\mathbf{r}_{12}$ vector that connects ions 1 and 2 with respect to the z-axis and $\hat{m}_k^\pm = \hat{m}_k^x \pm i\hat{m}_k^y$ the ladder operator of spin $k$.

The first term $\hat{A}$ adds to each energy level of the system a different frequency offset. This term is sometimes stamped as ‘z-coupling’ and can induce single ion inhomogeneous broadening of transitions [38]. The $\hat{B}$ term describes the mechanism through which the two interacting ions exchange their spin states, often referred to as flip-flop mechanism. For two identical ions, this mechanism is energy conserving, however it is not conserving for the RE-Y interaction. The $\hat{C}$ part allows one ion to change its spin state while the other remains unchanged. Finally, the $\hat{D}$ term is the double spin flip mechanism. Note that the conjugate interaction terms are included in $\hat{C}^\dagger$ and $\hat{D}^\dagger$. In our case, some of the mechanisms
grasped with this picture, sketched in the inset of Fig. 1, and highlighted with blue, red, and purple arrows, consist in:

- Blue arrows: the central ion can spin-flip, due to either flip-flop with one or several ions in the bath ($\hat{B}$ term in $\hat{H}_{\text{RE}-Y_k}^{\text{dd}}$, dashed blue arrow) or to an off-resonant single spin-flip ($\hat{m}_{Y_k}^z \hat{m}_{\text{RE}}^{\pm}$ terms in $\hat{H}_{\text{RE}-Y_k}^{\text{dd}}$, solid blue arrow);

- Red arrows: the bath ions can spin-flip, due to either off-resonant central ion coupling ($\hat{m}_{\text{RE}}^z \hat{m}_{Y_k}^{\pm}$ terms in $\hat{H}_{\text{RE}-Y_k}^{\text{dd}}$) or to off-resonant bath coupling ($\hat{C}$ and $\hat{C}^\dagger$ terms in $\hat{H}_{Y_i-Y_j}^{\text{dd}}$);

- Purple arrows: the bath ions can flip-flop, due to resonant bath coupling ($\hat{B}$ term in $\hat{H}_{Y_i-Y_j}^{\text{dd}}$) or to central ion-mediated interaction ($\hat{m}_{\text{RE}}^z \hat{m}_{Y_i}^+$ with $\hat{m}_{\text{RE}}^z \hat{m}_{-Y_j}^-$ involved respectively in $\hat{H}_{\text{RE}-Y_i}^{\text{dd}}$ and $\hat{H}_{\text{RE}-Y_j}^{\text{dd}}$).

<figure><img src="image_3.png" /><figcaption>FIG. 2: Eigenenergies of the bath $\hat{H}_Y$ (given by Eq.(4)) in the case $N=2$, as a function of the magnetic field magnitude.</figcaption></figure>

In order to assess the strength of the different contributions in the Y bath interaction Hamiltonian $\hat{H}_{Y_i-Y_j}^{\text{dd}}$, we plot in Fig. 2 the eigenvalues of $\hat{H}_Y$ for the first two Y neighbor ions in blue, and by removing the non-resonant contributions $\hat{C}$ and $\hat{D}$ in dashed red. We see that as soon as the ambient magnetic field $B_{\text{DC}}$ is larger than a few micro-Tesla, the Zeeman splitting is already sufficient to prevent non-resonant processes to occur via direct dipole-dipole coupling. Therefore, in the following, we neglect these contributions and restrict the study only to z-coupling ($\hat{A}$ term) and flip-flop ($\hat{B}$ term) processes within the bath alone: this is the only hypothesis that we make in the following.

In the next section, we relate our work to considerations found in the literature, where hypothesis were taken to conduct estimations of coherence times, and show why we cannot map them to our study.

## 2. *Usual approximations not made here*

If we focus on the central ion, an usual and efficient way of estimating decay rates in the case of a system coupled to a quasi-continuum bath is to rely on the Fermi golden rule [42, 57]. This has been for instance used to estimate flip-flop rates of Y spins in the bath [58]. This approximation is frequent in models that follow a stochastic approach. The Y ions influence is often modeled as a mean fluctuating magnetic field emerging from the Y-Y dynamics [38, 58], as a sudden jump model of Y spin changes [37], or as random noise [59]. All three models act as a perturbation that drives the decoherence on the central spin. Here we aim at estimating something slightly different: the central ion is coupled in a discrete manner to each state of the bath, with highly non-constant couplings between all levels. Moreover, no quasi-continuum hypothesis can be done in our case due to the purely deterministic landscape of the surrounding Y spins for each RE ion. Discrepancies between the Fermi golden rule and our model are addressed in App. B1.

Another very common hypothesis done on the systems is the frozen central spin, in which spin flips of the central ion are neglected (mechanisms depicted with a blue arrow in Fig. 1). In this case, the decoherence is estimated through the sole bath dynamics during an echo sequence. Interestingly, in this regime, explicit dynamics formula can be given by relying on methods such as the CCE [39, 43, 44]. This technique relies on the decomposition of the bath into clusters of increasing size, allowing to infer the explicit dynamics within each cluster and to estimate with more and more precision the bath effect on the central spin when increasing the cluster size (and therefore the expansion order). For the simplest *Loschmidt* echo decay, the decoherence over time is given by the overlap of the bath wavefunction with itself after free evolution whether the central spin is in its ground or excited state [60]. In the present study, this approach cannot grasp the full dynamics of the system due to the low nuclear dipole moments of the central ions considered here. Indeed, the scaling difference between the Y nuclear moment ($\sim 2$ kHz/mT) and Eu and Pr ($\sim 10-100$ kHz/mT) is not large enough to prevent multi-spin contributions for flipping the central ion. In the low field regime, we will show that this is indeed the case. For Yb, the zero-field splitting however prevents these events and central ion flips can be neglected at low fields.

Recently, a mean-dipole model has also been developed to describe the echo decay dynamics in erbium-doped YSO [41, 47]. Here, the frozen central spin hypothesis is also made, and no bath interaction is considered. We show in App. B2a that this approach is exactly a first-order CCE.
<figure><img src="image_4.png" /><figcaption>FIG. 3: Sequence under study for simulating single-ion coherence decays: after creation of an initial coherence with a radiofrequency pulse, the coherence (16) is plotted as a function of time (coherence decay).</figcaption></figure>

### 3. The sequence under consideration

Finally, we define the sequence that we will consider to assess a coherence time. Our approach aims at modeling decoherence via the interaction of the central ion with a bath. As it is a single-spin decoherence process, no inhomogeneity between RE ions of the ensembles is considered. Therefore, a simple free-induction decay (FID)-type sequence is used to simulate the effective decoherence of the central RE spin, in the same philosophy as when considering the previously mentioned Loschmidt echo. We however keep in mind that single-ion static broadening can be induced by sole z-coupling due to the $\hat{A}$ term in Eq. (12). This term effectively leads to frequency shifts of spin-preserving transitions that beat due to balanced statistical mixture of the initial bath population. To ensure that such effect is not limiting the coherence time, we can plot Hahn echo sequence and witness that similar coherence time is obtained (see App. E). However, computation resources required for simulating such a sequence is very demanding, therefore we focus on FID sequences whenever possible.

The considered sequence is depicted in Fig. 3. Initially, only the ground hyperfine level is populated and Zeeman (for Eu and Pr) and superhyperfine energy levels are in a balanced mixed state. A coherence is then created in the system with an RF or MW pulse on the relevant spin transition of the central RE ion (green arrows in Fig. 1). This is illustrated in Fig. 3 with the representation of the absolute value of the density matrix before and after the RF pulse for one Pr ion and two Y bath ions. Note that in the subsequent plots we discard this time interval, such that we only plot the following decay, that we will from now on call the coherence decay (time after $t = 0$ in Fig. 3). Next, the system is left to evolve under the

sole drive of the dipole-dipole interaction term of $\hat{H}_{\text{int}}$. The evolution of the state vector can then be explicitly written at each instant:

$$|\psi(t)\rangle = \sum_n \sum_{k,l} \alpha_k(0) p_{n,l} p_{kl}^* e^{i\lambda_l t} |n\rangle \quad (13)$$

$$= \sum_n \alpha_n(t) |n\rangle, \quad (14)$$

where the $\alpha_n(t)$ are the time-dependent components of the state vector of the system along the eigenvector $|n\rangle$ of $\hat{H}_0$. The elements $p_{i,j}$ and $\lambda_l$ are respectively the coefficients of the base change matrix $P$ and of the diagonal matrix $D$, such that, in its matrix form, the total Hamiltonian reads:

$$H = PDP^{-1}, \quad (15)$$

where the interaction Hamiltonian $\hat{H}_{\text{int}}$ only consists in the dipole-dipole interaction (no RF drive during the FID).

The coherence decay intensity is then read optically at any time $t$ and computed as:

$$I(t) = |\text{Tr}[\hat{\rho}(t)\hat{\mu}^{\text{opt}}]|^2, \quad (16)$$

where $\hat{\mu}^{\text{opt}}$ is the optical moment of the RE ion. This formula emulates an optical readout of the coherence, as it is usually performed experimentally via the Raman heterodyne scattering (RHS) protocol [61]. Note that here we assume that the hyperfine population is maintained by the transfer pulse for optical readout (the density matrix in the spin transition is assumed to be the same as in the optical one). We however estimate that such consideration is valid for reproducing decay curves in qualitative agreement with experimental data.

Injecting expression (13) in (16) gives the coherence decay in the pure state case:

$$\begin{aligned} I(t) &= |\langle\psi(t)|\hat{\mu}^{\text{opt}}|\psi(t)\rangle|^2 \\ &= \left| \sum_{\substack{n,m \\ k,l,k',l'}} \mu_{n,m}^{\text{opt}} \alpha_k^*(0) \alpha_{k'}(0) p_{n,l}^* p_{k,l} p_{m,l'} p_{k',l'}^* e^{i(\lambda_{l'} - \lambda_l)t} \right|^2 \\ &= \left| \sum_{l,l'} g_{l,l'} e^{i(\lambda_{l'} - \lambda_l)t} \right|^2. \end{aligned} \quad (17)$$

This expression shows that the frequencies appearing in the dynamics of the coherence decay are the interval between the eigenvalues of $\hat{H}$ noted $\Delta\lambda_{l,l'} = |\lambda_l - \lambda_{l'}|$ weighted by a function related to the coupling between the eigenvectors (coefficients $p_{i,j}$) and to the initial coherences $\alpha_k^*(0)\alpha_{k'}(0)$. We will use this formula to visualize relevant weighted frequencies appearing in the dynamics of the echo. From now on, frequencies $\Delta\lambda_{l,l'}$ are called the dynamical frequencies.

Finally, when a coherence time shall be inferred and
displayed for a given coherence decay trace, we follow the procedure detailed in App. A. Such a procedure is required due to the very changing nature of the decay temporal profile with respect to magnetic field properties.

For the following sections, we perform our simulations on the three RE ions presented above. By identifying the dynamical frequencies at stake, we will for each case determine the predominant decoherence mechanism among the ones presented in Fig. 1.

### III. RESULTS

#### A. Europium

The first system under study is europium. We confront in the following our numerical simulations with experimentally acquired data points with the spin-wave AFC protocol, with a setup showed in Fig. 4.

<figure><img src="image_5.png" /><figcaption>FIG. 4: Experimental setup for the AFC spin-wave storage in Eu:YSO. See body of the text for details.</figcaption></figure>

A 1000 ppm Eu:YSO crystal is cooled to 4 K in a closed-cycle cryostat. Population manipulations are performed to prepare a 3 MHz wide atomic frequency comb on the $|\pm 5/2\rangle_g \leftrightarrow |\pm 1/2\rangle_e$ transition with a preparation beam (gray Prep beam in Fig. 4), allowing to reach a double-pass optical depth of 6. Then the spin-wave storage of classical pulses of light is performed in the ensemble, by sending an optical input pulse on the $|\pm 5/2\rangle_g \leftrightarrow |\pm 1/2\rangle_e$ transition (orange In beam in Fig. 4), followed by a 15 µs-long optical control pulse on the $|\pm 3/2\rangle_g \leftrightarrow |\pm 1/2\rangle_e$ transition (blue Ctrl beam in Fig. 4). Note that all beams propagate along a direction close to the *b* axis (offset angles serve for spatial filtering) and are all polarized along $D_1$. Spin rephasing is then performed by applying an RF XY-4 sequence on the 46.2 MHz $|\pm 5/2\rangle_g \leftrightarrow |\pm 3/2\rangle_g$ transition. The RF field is applied thanks to a coil wrapped around the crystal. Details about pulse shape, preparation procedure and pulse sequences can be found in [32] and its Sup. Mat., for which the same experimental setup and sequence was used. An additional static magnetic field of varying amplitude is applied along $D_1$. It is important to mention that in the experimental sequence, the RF pulses compensate for the effect of the inhomogeneous spin decoherence, such that the plotted decay curves are limited

by single-ion decoherence, which can be compared to the model discussed in this paper.

The decays of the spin-wave storage of classical pulses of light are presented with orange dots in Fig. 5a for bias magnetic field amplitudes of 1, 10, 40 and 60 µT. These curves display a clear dependency on the applied external field, as well as modulation patterns.

To simulate this behavior, we use the previously introduced model and sequence, with an hypothesis on the initial coherence creation pulse: in the experimental case, the initial coherence is indeed created through optical excitation followed by an optical transfer pulse, each of the two pulses having very different parameters. This is in contrast with our simulations, where all calculations and simulations are done on the spin transitions. Therefore, we chose to mimic this by using a 200 µs FWHM adiabatic pulse on the spin transition, chirped by 50 kHz. Such a choice essentially conditions the depth of the modulations at the Zeeman frequency but not the overall shape of the coherence decay [62, 63]. The simulated curves are shown in blue, by considering the first five Y neighbors. Even if the simulated sequence is much simpler than the actual full spin-wave protocol, we can observe two important correspondences.

The first agreement is the global decay time, which is well reproduced: an increase in the amplitude of the magnetic field leads to a longer coherence time, from less than a ms at low field to 5 ms at 60 µT. Here we remind that no empirical decay is introduced manually in the shape of the coherence decay (17) (no $e^{-(t/\tau)^\alpha}$ term introduced): it is solely due to path interference in the $\hat{H}_0$ manifold driven by the $\hat{H}_{int}$ Hamiltonian. Therefore, our simulations always present long-duration coherence revivals, that are pushed further and further by considering more and more bath ions.

The second observation is the faster modulation present in the dynamics, at a frequency of the Eu Zeeman splitting (≈ 14 kHz/mT), that we also reproduce numerically. The contrast of this modulation also diminishes when increasing the magnetic field amplitude. This is mainly due to the decrease of the probability of having a central spin flip during the coherence creation with an increasing field, as we showed in a previous publication [48] (both Zeeman splittings of the spin transitions are actually involved in the evolution). Moreover, the experimental decay pattern strongly depends on the rephasing pulse position, that we do not take into account due to the simplified sequence considered here [48]. It can also be observed that the mismatch between the curves is higher at low magnetic fields, which we attribute to a difficulty to precisely estimate experimentally the applied magnetic field for such low values. For quick comparison with the first-order CCE approximation, we also show the curves obtained with the first 30 Y neighbors with black dashed lines in Fig. 5a. The simulated CCE curves clearly display longer decays, as they do not grasp the main decoherence mechanisms for this RE, as we will see below.
<figure><img src="image_6.png" /><figcaption>FIG. 5: a) Experimental (orange dots) and simulated (blue lines) echo decay curves for Eu:YSO for the first $N = 5$ neighboring Y ions, under a magnetic field oriented along $D_1$ set to an amplitude of $B_{DC} = 1, 10, 40$ and $60$ µT. For comparison, decays obtained with first order CCE are shown in dash black (Hahn echo) and solid black (Loschmidt echo). b) Simulation of the dynamical frequencies involved in the echo decays as a function of the magnetic field amplitude for Eu surrounded by two Y ions. The color scale represents the normalized *g* weights (see Eq. (17)) associated with each frequency. Blue arrows indicate Eu spin flip processes, red arrows indicate Y spin flips, and purple arrows indicate Y flip-flop processes. Inset: plotted frequencies without color scale.</figcaption></figure>

It should also be pointed that in the experimental echo decays, one revival stands precisely at the yttrium Zeeman frequency (3.6 ms at 40 µT and 6 ms at 60 µT), which is not reproduced by our simulations. This discrepancy leads us to investigate more in details the different frequencies involved in the dynamics, as well as their relative weights. Figure 5b shows the frequencies $\Delta\lambda_{l,l'}$ and their respective frequency weights $g_{l,l'}$ of Eq. (17) (color scale), as a function of the magnetic field amplitude between 0 and $60$ µT. In the inset, all the dynamical frequencies $\Delta\lambda_{l,l'}$ are plotted without the color scale. The frequency weights show the different processes at play depending on the magnetic field amplitude. Up to $10$ µT, all Eu spin flip (blue arrows), Y spin flip (red arrows), and Y flip-flop processes (purple arrows) are non-negligible. Beyond $10$ µT and up to $40$ µT, as the frequency gap between the Eu Zeeman level increases, the Eu-Y dipole-dipole interactions are not sufficient to induce Eu spin-flips anymore. Above $40$ µT, Y spin-flips are forbidden for the same reason, and only Y flip-flops contribute. Among these mechanisms, only the Y spin flips are taken into account in the first-order CCE, explaining the previously mentioned mismatch.

Above $40$ µT, predominant dynamical frequencies stay mostly constant with the field, leading to a constant characteristic time in the echo decays (see Fig. 5a). It should also be pointed that Fig. 5b is plotted with two Y ions, whereas the echo decays of Fig. 5a are simulated with five Y ions. Therefore the actual amount of dynamical frequencies giving rise to the decay is higher than the one shown here. However, this graph indicates which process

is relevant in each field regime. To further strengthen this interpretation, we propose in App. C to show population dynamics that reveal the aforementioned spin flips. The last necessary validation concerns the pertinence of the approach consisting in simulating a simple FID process with respect to a full echo sequence. As mentioned in Sec. IIB 1, the $\hat{A}$ term in the dipole-dipole coupling creates a level-dependent energy shift: if the system is in an initial mixed state, it will dephase quickly as it is the case for an inhomogeneously broadened medium. This occurs even if we consider a single ion interacting with the bath. Therefore, we should ensure that the limiting decay time (essentially independent of the magnetic field) is long enough as compared to the simulated one. To this end, we plot in Fig. 5a the decays associated with Loschmidt echoes [60], in which static decoherence is not compensated by a rephasing pulse as in either Hahn echo or spin-wave AFC sequences. The curves are represented in solid black. They reveal that the simulated coherence decay begins to be similar to Loschmidt echoes above $40$ µT, showing that simulation at higher fields would require to introduce rephasing pulses that also counteract single-ion inhomogeneous dephasing. We confirm the validity of our approach at fields up to $60$ µT in App. E by comparing a Hahn echo sequence with our simulated sequence, and witnessing a good correspondence.
<figure><img src="image_7.png" /><figcaption>FIG. 6: a) Yb echo decays at various magnetic field amplitudes. Simulations for one Yb ion and 7 Y ions (in blue) and experimental data (in orange) [35]. b) Dynamical frequencies involved in the echo dynamics as a function of the magnetic field amplitude. The color scale represents the *g* weight (see Eq. (17)) associated with each frequency (normalized). Inset: frequencies involved in the echo dynamics, without color scaling. The red arrows point out components associated with Y spin flips processes, and purple arrows points components associated with Y flip-flops processes. In both figures, the blue dashed line represents the 2.1 kHz/mT Y Zeeman frequency found in the experimental echo decay revivals of Fig. 6a.</figcaption></figure>

## B. Ytterbium

The second RE ion investigated is ytterbium (Yb:YSO). Recent studies have revealed long coherence times in particular magnetic configurations [35, 64], and we propose here to confront our simulations with previously published experimental echo decays. Experimental data points and details about the used sequence can all be found in [35]. The experimental measurement of the coherence is performed using a Hahn echo sequence and detected through coherent RHS techniques. Both the experiments and our simulations are performed on the 2.497 GHz transition as stated before, for the optical site II (see Fig. 1).

In our first simulations, we aim at replicating the experimental echo decays recorded at varying magnetic field amplitude along the *D*₁ direction shown in Fig. 6a in orange. We perform our simulations with seven Y ions and the decays are represented in blue in the figure. For magnetic fields above ~ 60 µT, we obtain a good correspondence in terms of initial decay time between the curves. For magnetic field amplitudes below 60 µT, the simulation predicts a longer decay time than experimentally measured. This difference can be explained by the difficulty to precisely estimate the magnetic field experimentally (in the same way as for europium in the previous section), and by the fact that in this regime the superhyperfine interaction might not be the limiting process anymore. In addition to the overall good match in the decay characteristic time, we also notice that the initial horizontal slope is well reproduced by our simulations. This is explained by the fact that in our approach, the decay is composed of multiple sinusoidal contributions,

each with a zero slope at zero delay. This results in Gaussian-like decays, that are in good agreement with experimental data, as illustrated in Fig. 6a.

<figure><img src="image_8.png" /><figcaption>FIG. 7: Dependence of the simulated decay time in Yb:YSO as a function of the magnetic field amplitude. *T*₂ red dots values fitted from the simulated coherence decays of Fig. 6a. Black line is a fit of red dots, with a law $T_2 = A_1/B_{D_1}$ and $A_1 = 1.36 \cdot 10^{-7}$ s.T.</figcaption></figure>

The second point highlighted by these curves is the absence of revivals in the simulated decays as compared to experimental ones, in which a clear revival close to the Y Zeeman splitting uncoupled frequency (2.1 kHz/mT) is visible. Two possible explanations are proposed. On one hand, the mismatch could come from the initial state, that may not be in reality a perfectly balanced mixed state of the levels inside the Y manifold (see the $\alpha_k(0)\alpha'_k(0)$ term in Eq. (17)). The consequence for this is that different frequency components involved in the dynamics can interfere in a different way given their different weights. On the other hand, we do not simulate an
<figure><img src="image_9.png" /><figcaption>FIG. 8: a) Simulated coherence time (color map, in ms) as a function of magnetic field angles $\varphi$ and $\theta$ (see Eq. (2)) for one Pr ion and three Y ions and for magnetic field amplitudes of $B_{DC} = 0.07, 0.35$ and $1$ mT. Simulations are performed for one magnetic subsite. Echo decays associated with colored symbols in Fig. (a) are represented on: b) for the dot (orange), c) for the cross (purple), and d) for the triangle (green). For comparison, all echo decay simulations were also performed with five Y ions (black curves).</figcaption></figure>

optical pulse to transfer the coherence from one ground hyperfine Yb level to an excited hyperfine level. Instead we assume that we have a perfect population transfer from the $|2\rangle_g$ state to the $|2\rangle_e$ state, similarly as for europium (see Sec. III A). This idealized population transfer in our simulations may also affect the dynamics.

Then, similarly as for Eu, we associate the frequencies observed in the dynamics with spin processes. To this extent, we plot in Fig. 6b the evolution of the dynamical frequencies $\Delta\lambda_{l,\nu}$ as a function of the magnetic field amplitude for one Yb ion and two Y ions, still for the same field orientation (aligned along $D_1$). Their respective weights $g$ are normalized and represented with a color scale. The Y Zeeman frequency 2.1 kHz/mT is plotted as a blue dashed line, and it clearly appears that such a frequency is actually not present in the spin dynamics, corroborating the observations made with the simulated decays in Fig. 6a. Frequencies with the strongest weights stand at lower frequencies, at 0.67 kHz/mT and 0.3 kHz/mT. To associate each curve with a physical process, we plot in the inset of Fig. 6b all dynamical frequencies (unweighted, as in Fig. 5b). In this inset, three groups of curves can be identified. The highest energy branches correspond to double Y spin flips (2 branches

that are indistinguishable on the figure), then a group of 8 branches (3 visible on the figure) correspond to single Y spin flips. Eventually, at lower energies, 4 branches (3 visible on the figure) are associated with Y flip-flops. Notice here that these branches are actually field-dependent due to the linear dependence of the magnetic dipole of Yb with the field, on the contrary to the previous section with Eu. This feature is addressed in detail in the Supp. Mat. of [35]. In our simulated decays, we see that the only components with non-negligible weights are the Y flip-flop components, driven at a frequency essentially governed by the Yb-Y dipole-dipole coupling (one order of magnitude larger than simple Y-Y dipole-dipole interaction).

Another important difference of this system with Eu is the constant weights of the involved dynamical frequencies, because of the linearly increasing magnetic dipole of Yb with the field. This implies that augmenting the field amplitude increases the dynamical frequencies at stake while maintaining their weight constant. The consequence for this is that the frequency map remains identical whatever the field up to a global scaling factor. In other words, changing the field by a factor $K$ amounts to changing time by a factor $1/K$. To quantify this ef-
fect more explicitly, we extract from our four simulated coherence decay curves (Fig. 6a) a coherence time $T_2$ by fitting them with a stretched exponential $\exp[-(t/T_2)^\beta]$. The resulting $T_2$ are plotted in Fig. 7 with red dots and are well fitted with a $1/B$ scaling law, as illustrated by the black solid line.

### C. Praseodymium

In this section, we will show numerical predictions of spin coherence times for Pr:YSO (see Fig. 1 for the energetic structure). As stated in Sec. II A, we will focus on the spin transition $|\pm 1/2\rangle_g \leftrightarrow |\pm 3/2\rangle_g$ at 10.2 MHz.

We plot in Fig. 8a the coherence time as a function of the magnetic field angles $\varphi \in [0, 2\pi]$ and $\theta \in [0, \pi]$ (see Eq. (2) for angles definition), for magnetic field amplitudes of $B_{DC} = 0.07, 0.35$ and 1 mT. In this graph, due to simulation duration constraints, we limit the number of Y ions in the system to 3, and the simulation time range to 6.5 ms. We recall that the coherence time is inferred with the method presented in App. A. The first information revealed by these plots is the high anisotropy of the coherence time with the field orientation: at magnetic fields above 0.35 mT, there is a factor of ten between low and high coherence time regions ($\approx 0.5$ to 5 ms). This feature is to be linked with the high anisotropy of the magnetic dipole of Pr:YSO, with a factor of more than 4 difference along the three directions of space [45]. The second information is that the coherence time increases when the magnetic field amplitude (up to 0.35 mT here), as witnessed experimentally [36, 65]. However, the plots show that above this value, increasing the magnetic field does not lead to a coherence time extension anymore. This observation is in agreement with simulations performed with Eu:YSO in Sec. III A, in the regime where Y spin flip-flops are the only remaining decoherence mechanism: both REs have a field amplitude independent magnetic dipole (as long as the quadratic Zeeman contribution remains negligible).

In order to further investigate the improvement in the coherence time, we plot the decays in three configurations indicated by colored symbols in Fig. 8a. These are represented in Fig. 8b for the orange dot, Fig. 8c for the purple cross, and Fig. 8d for the green triangle. Notice that for these particular curves, we push the simulation duration up to 8 ms to visualize the time dependence in more details. The curves reveal the non-monotonous nature of the decay and display, similarly as for Eu:YSO (see Sec. III A), a strong component at the Pr Zeeman frequencies. In the same manner, the depth of the modulation is directly conditioned by the initial state, that we however don't aim at adjusting here.

Then, in order to assess if limiting the simulation to the first three neighboring Y ions is sufficient, we plot the echoes obtained with five ions in black on top of the three Y ions curves. Adding more Y ions does not change the hierarchy of decay times (highest coherence regions identified

in Fig. 8a still lead to the longest coherence times with 5 Y ions) but modifies the decay shape. In particular, strong revivals appearing after a few milliseconds are blurred by the introduction of additional accessible energies in the system. Notice that this is also the case for Yb:YSO, in which strong revival can appear in decay simulations if the number of surrounding Y taken into account is insufficient.

In Fig. 8c, we also clearly identify beatings at the yttrium Zeeman frequency (particularly visible at 0.35 mT). The frequency components of the corresponding decays are shown in Fig 9a, and correspond to spectral structures that have also been measured experimentally in Pr:YSO [36], and associated with single or multiple Y spin flips by the authors. We confirm here this interpretation numerically.

Finally, in order to gain insight into the shape of the region of long coherence times (sinusoidal like shape in Fig. 8a), we plot the average of the norm of the four gradients $||S_1||$ of the transitions involved in the four distinct spin transitions $|\pm 1/2\rangle_g \leftrightarrow |\pm 3/2\rangle_g$ in Fig. 9b. In this plot, we consider both magnetic subsites when calculating the average $||S_1||$ because when $\theta \notin \{0, \pi/2\}$, both subsites are nonequivalent. The concordance of the shapes for subsite 1 shows that the sinusoidal pattern indeed mostly originates from a weaker sensitivity of the transitions to external fluctuations. However, the longest coherence times simulated in Fig. 8a are not always superimposed with the minimum of $||S_1||$, indicating that a better optimization of the coherence time could be achieved with our model. Experimentally, we expect the decays to consist of contributions of both subsites. Therefore, a simple working point to be identified in these plots is at the crossing of the two low $||S_1||$ zones in Fig. 9b. This point corresponds to $\theta = \pi/2$ (field in the $(D_1, D_2)$ plane), where the two magnetic subsites are equivalent by symmetry, with a corresponding decay plotted in Fig. 8b. We focus on this region by also plotting the coherence time as a function of $\varphi$ and the field amplitude in Fig. 9c, showing that it is for $\varphi$ close to $\pi$ (field along $D_1$) that the coherence time is the highest. It also gives the value of the field ($\sim 200$ µT) for which Y spin flip-flops remain the only decoherence mechanism. This is five times higher than for Eu, in agreement with the strength of the RE-Y dipole-dipole coupling as we will see in the next section.

## IV. DISCUSSION

### A. Impact of the Y ions on the coherence times

Both the simulations and the experimental data presented previously showed that the coherence times in RE-Y systems are highly dependent on the magnetic field amplitude and orientation, as witnessed for Pr in Fig. 8a and for Yb in App. D, Fig. 14a. Due to the anisotropy of the magnetic moments of the RE ions, as the mag-
<figure><img src="image_10.png" /><figcaption>FIG. 9: a) Fourier transform of the echo decays presented in Fig. 8c. The red arrows indicate the frequencies corresponding to the 2.1 kHz/mT yttrium nuclear magnetic moment. b) Mean $||S_1||$ (MHz/T) of the four transitions $|\pm 1/2\rangle_g \leftrightarrow |\pm 3/2\rangle_g$ for each magnetic subsite as a function of the magnetic field direction, for an amplitude of 1 mT. c) Coherence time (color map, in ms) as a function of magnetic field angle $\phi$ (in rad) and amplitude for one Pr ion and three Y ions in the ($D_1, D_2$) plane.</figcaption></figure>

netic field orientation changes, so does the strength of the RE ion effective magnetic dipole and consequently, the dipole-dipole RE-Y interaction strength.

In the case of Yb, as the magnetic dipole linearly increases with the applied magnetic field, the coupling with the surrounding Y ions increases accordingly. This leads to the contribution of dynamical frequencies of increasing frequencies with the field, with weights that do not vanish when increasing the amplitude of the field [35]. On the contrary, for a given magnetic field orientation, the magnetic dipoles of the Eu and Pr ions are constant regardless of the magnetic field amplitude. The constant strength of the Eu-Y and Pr-Y dipole-dipole interactions result in a different impact on the whole RE-Y dynamics, depending on the field amplitude. This leads us to identify three magnetic field regimes for the Eu and Pr ions.

a. Very low field regime: When both the RE Zeeman splitting $\delta_{RE}$ and the Y Zeeman splitting $\delta_Y$ are smaller than the RE-Y magnetic dipole-dipole coupling, we are in the very low field regime. In this regime, RE-Y interactions can cause RE spin flips (blue arrows in Fig. 1), Y spin flips (red arrows in Fig. 1) and Y flip-flops (purple arrows in Fig. 1). This causes the dynamics of the system to present many frequencies with non-negligible weights $g$ (see Eq. (17)), which may also span a larger frequency range than for higher fields. Both effects lead to fast decays, and are observed in Fig. 5 for Eu, for magnetic field amplitudes below 10 $\mu$T.

b. Low field regime: When the Y Zeeman splittings $\delta_Y$ are of the same order of magnitude than the dipole-dipole RE-Y coupling, we enter the low magnetic field regime where the RE-Y couplings are not able to induce RE spin flips anymore. We saw that this regime is situated between 10 $\mu$T and 40 $\mu$T for Eu (field along $D_1$).

c. Intermediate field regime: Finally, when both $\delta_{RE}$ and $\delta_Y$ are larger than the RE-Y coupling, we enter the intermediate magnetic field regime. Dipole-dipole couplings are not able to cause any spin flips but can only drive Y flip-flops. If the magnetic dipole of the RE ion does not scale with the magnetic field amplitude, such as for Eu and Pr, the coherence time will become mostly independent on the magnetic field amplitude. Indeed, the frequency splittings created by the dipole-dipole interactions are independent of the field amplitude, and consequently the spectral range explored in the system dynamics will remain unchanged. This spectral range being narrower than the spectral range explored in low and very low magnetic field regimes, the intermediate magnetic field regime is the one that provides the longest coherence times among the three regimes investigated here.

## B. Perspectives

Our work has allowed us to identify decoherence mechanisms in different magnetic field regimes in RE-doped YSO. Simulations have helped to identify favorable regimes, by mainly focusing on fields applied in the ($D_1, D_2$) plane for symmetry reasons. A deeper investigation of the interplay between the two magnetic subsites could potentially help pushing the coherence times even further (see for instance high coherence time zones of subsite 1 in Fig. 8). The sequence that we have used to simulate the coherence times shall also be discussed: in the very low field regime, only a few hypotheses can be done on the total Hamiltonian, therefore full dynamics simulation must be
done for estimating the coherence time. Limited computation power ultimately limits the complexity of the sequences that are simulated, as well as the size of the Hilbert space (and therefore the number of Y ions). In the very low field regime, we have seen that the coherence decays could reproduce well the experimental echo decays, indicating that taking into account only a few Y ions is sufficient. However, we predict that if one pushes the field further, the simplistic sequence used here will not reflect anymore actual coherence times. This is mainly due to the single-ion z-coupling contribution discussed at the end of Sec. III A. However, entering the intermediate field regimes also means that spin flips are forbidden, and therefore that frozen central spin approximation can be made. Given that first-order CCE does not allow to simulate ion bath flip-flop, we believe that higher-order CCE would be a viable and efficient way of simulating the system in such regime. It would also allow to focus more attention on the spin manipulation during spin rephasing sequences: imperfections in spin inversions can cause the rephasing sequence to be less efficient than expected [48], therefore a precise determination of the magnetic field regime and inversion pulse parameter is still to be precisely investigated.

Such investigation would also allow us to simulate more complex sequences such as dynamical decoupling sequences, routinely applied to counteract spectral diffusion, that is usually attributed to Y flip-flops in the intermediate field regime [32, 33, 66]. In this regime, it is therefore instructive to illustrate the range of frequencies involved in such flip-flop dynamics. In our picture, we recall that the Y flip-flop mechanism gives rise to a lift of degeneracy of equal total Y spin projection states (see Fig. 1). If we plot the spread of frequencies spanned by the distribution of the first ten Y nearest neighbors, we obtain the histogram of Fig. 10 for the configuration 5 spins "up" and 5 spin "down", leading to a total of $C_{10}^5 = 252$ configurations. This distribution fits well with a Gaussian distribution with FWHM of 30 Hz. Links between this energy spread and spectral diffusion witnessed in non-Kramers ions deserve further investigation, but such a value is already in good agreement with typical linewidths considered in spectral diffusion models (35 Hz found for instance in [33]).

By pushing the field amplitude further, the quadratic Zeeman contribution (third term in Eq. (1)) becomes non negligible anymore, and we believe that our approach can provide new paths for the comprehension of mechanisms such as the 'frozen core' effect [27, 45].

V. CONCLUSION

By performing a small-scale exact treatment of the spin interactions in RE ion-doped crystals in the µT to mT magnetic field amplitude range, we have shown that we could understand and simulate the processes that induce decoherence and experimentally limit the coherence

<figure><img src="image_11.png" /><figcaption>FIG. 10: Histogram of the energy distribution of the $I_z = 0$ manifold of the first ten Y ions (blue histogram, occurrence in the vertical axis). The energy spread is due to Y-Y magnetic dipole coupling (see second term in Eq. (4)), and fits with a Gaussian of FWHM 30 Hz (solid red).</figcaption></figure>

times. We have validated our approach with a side-to-side comparison between simulations and experimental results for both Eu and Yb, and although slight differences in the dynamics have been witnessed, similar coherence decay times have been found in the magnetic field amplitude range of 0-100 µT, with no adjustable parameter. Remarkably, in Yb:YSO the $1/B$ dependency of the experimental echo decay time could be well reproduced numerically. Following this, we have used our simulator to make coherence time predictions as a function of magnetic field amplitude and orientation in a Pr:YSO system, and have identified a favorable set of magnetic field parameters that minimize the decoherence induced by the surrounding Y ions.

Our results have led to the identification of different magnetic field regimes, in which different spin flip processes induced by the RE-Y dipole-dipole interactions are predominant.

ACKNOWLEDGMENTS

The authors thank A. Holzäpfel, H. de Riedmatten, F. Appas, P. Goldner and A. Tiranov for fruitful discussions and important feedbacks on the simulations. JE acknowledges the National Research Agency (ANR) through the project WAQUAM (ANR-21-CE47-0001-01), France 2030 National Research Program on Quantum Technologies project Qmemo (ANR-22-PETQ-0010), Univ. Côte d'Azur CSI through the project FITRALP, Doeblin Federation through the project STALAR and DEP from Académie 2 Univ. Côte d'Azur. VDA acknowledges the
Institut Universitaire de France. MA acknowledges fund-
ing from the Swiss FNS NCCR program Quantum Sci-
ence Technology (QSIT), European Unio Horizon 2020
research and innovation program within the Flagship on

Quantum Technologies through GA 820445 (QIA) and
under the Marie Skłodowska-Curie program through GA
675662 (QCALL).

[1] N. Gisin, G. Ribordy, W. Tittel, and H. Zbinden, Quan-
tum cryptography, *Rev. Mod. Phys.* **74**, 145 (2002).

[2] A. Ekert and R. Renner, The ultimate physical limits of privacy, *Nature* **507**, 443 (2014).

[3] J. I. Cirac, A. K. Ekert, S. F. Huelga, and C. Macchiavello, Distributed quantum computation over noisy channels, *Phys. Rev. A* **59**, 4249 (1999).

[4] A. S. Cacciapuoti, M. Caleffi, F. Tafuri, F. S. Cataliotti, S. Gherardini, and G. Bianchi, Quantum internet: Networking challenges in distributed quantum computing, *IEEE Network* **34**, 137 (2020).

[5] P. Kómár, E. M. Kessler, M. Bishof, L. Jiang, A. S. Sørensen, J. Ye, and M. D. Lukin, A quantum network of clocks, *Nature Physics* **10**, 582 (2014).

[6] D. Gottesman, T. Jennewein, and S. Croke, Longer-
baseline telescopes using quantum repeaters, *Phys. Rev.
Lett.* **109**, 070503 (2012).

[7] B. Hensen, H. Bernien, A. E. Dréau, A. Reiserer, N. Kalb, M. S. Blok, J. Ruitenberga, R. F. L. Vermeulen, R. N. Schouten, C. Abellán, W. Amaya, V. Pruneri, M. W. Mitchell, M. Markham, D. J. Twitchen, D. Elkouss, S. Wehner, T. H. Taminiau, and R. Hanson, Loophole-free Bell inequality violation using electron spins separated by 1.3 kilometres, *Nature* **526**, 682 (2015).

[8] W. Rosenfeld, D. Burchardt, R. Garthoff, K. Redeker, N. Ortegel, M. Rau, and H. Weinfurter, Event-ready bell test using entangled atoms simultaneously closing detection and locality loopholes, *Phys. Rev. Lett.* **119**, 010402 (2017).

[9] T. Inagaki, N. Matsuda, O. Tadanaga, M. Asobe, and H. Takesue, Entanglement distribution over 300 km of fiber, *Opt. Express* **21**, 23241 (2013).

[10] A. Reiserer, N. Kalb, M. S. Blok, K. J. M. van Bem-
melen, T. H. Taminiau, R. Hanson, D. J. Twitchen, and
M. Markham, Robust quantum-network memory using
decoherence-protected subspaces of nuclear spins, *Phys.
Rev. X* **6**, 021040 (2016).

[11] W. Kozlowski, A. Dahlberg, and S. Wehner, Designing a quantum network protocol, in *Proceedings of the 16th International Conference on Emerging Networking Experiments and Technologies*, CoNEXT '20 (Association for Computing Machinery, New York, NY, USA, 2020) p. 1–16.

[12] H.-J. Briegel, W. Dür, J. I. Cirac, and P. Zoller, Quan-
tum repeaters: The role of imperfect local operations
in quantum communication, *Phys. Rev. Lett.* **81**, 5932
(1998).

[13] S. L. N. Hermans, M. Pompili, H. K. C. Beukers, S. Baier, J. Borregaard, and R. Hanson, Qubit teleportation between non-neighbouring nodes in a quantum network, *Nature* **605**, 663 (2022).

[14] Z.-Q. Zhou, D.-L. Chen, M. Jin, L. Zheng, Y.-Z. Ma, T. Tu, A. Ferrier, P. Goldner, C.-F. Li, and G.-C. Guo, A transportable long-lived coherent memory for light pulses, *Science Bulletin* **67**, 2042 (2022).

[15] W. Zhang, D.-S. Ding, Y.-B. Sheng, L. Zhou, B.-S. Shi, and G.-C. Guo, Quantum secure direct communication with quantum memory, *Phys. Rev. Lett.* **118**, 220501 (2017).

[16] N. H. Lindner and T. Rudolph, Proposal for pulsed on-demand sources of photonic cluster state strings, *Phys. Rev. Lett.* **103**, 113602 (2009).

[17] I. Schwartz, D. Cogan, E. R. Schmidgall, Y. Don, L. Gantz, O. Kenneth, N. H. Lindner, and D. Ger-shoni, Deterministic generation of a cluster state of entangled photons, *Science* **354**, 434 (2016), https://www.science.org/doi/pdf/10.1126/science.aah4758.

[18] T. Chanelière, D. N. Matsukevich, S. D. Jenkins, S.-Y. Lan, T. A. B. Kennedy, and A. Kuzmich, Storage and retrieval of single photons transmitted between remote quantum memories, *Nature* **438**, 833 (2005).

[19] G. Buser, R. Mottola, B. Cotting, J. Wolters, and P. Treutlein, Single-photon storage in a ground-state vapor cell quantum memory, *PRX Quantum* **3**, 020349 (2022).

[20] K. T. Kaczmarek, P. M. Ledingham, B. Brecht, S. E. Thomas, G. S. Thekkadath, O. Lazo-Arjona, J. H. D. Munns, E. Poem, A. Feizpour, D. J. Saunders, J. Nunn, and I. A. Walmsley, High-speed noise-free optical quantum memory, *Phys. Rev. A* **97**, 042316 (2018).

[21] R. Stassi and F. Nori, Long-lasting quantum memories: Extending the coherence time of superconducting artificial atoms in the ultrastrong-coupling regime, *Phys. Rev. A* **97**, 033823 (2018).

[22] M. Bouillard, G. Boucher, J. F. Ortas, B. Pointard, and R. Tualle-Brouri, Quantum storage of single-photon and two-photon fock states with an all-optical quantum memory., *Physical review letters* **122 21**, 210501 (2018).

[23] I. Craiciu, M. Lei, J. Rochman, J. G. Bartholomew, and A. Faraon, Multifunctional on-chip storage at telecommunication wavelength for quantum networks, *Optica* **8**, 114 (2021).

[24] A. Kinos, L. Rippe, D. Serrano, A. Walther, and S. Kröll, High-connectivity quantum processor nodes using single-ion qubits in rare-earth-ion-doped crystals, *Phys. Rev. A* **105**, 032603 (2022).

[25] Z. Wang, L. Balembois, M. Rančić, E. Billaud, M. Le Dantec, A. Ferrier, P. Goldner, S. Bertaina, T. Chanelière, D. Esteve, D. Vion, P. Bertet, and E. Flurin, Single-electron spin resonance detection by microwave photon counting, *Nature* **619**, 276 (2023).

[26] J. R. Everts, M. C. Berrington, R. L. Ahlefeldt, and J. J. Longdell, Microwave to optical photon conversion via fully concentrated rare-earth-ion crystals, *Phys. Rev. A* **99**, 063830 (2019).

[27] M. Zhong, M. P. Hedges, R. L. Ahlefeldt, J. G. Bartholomew, S. E. Beavan, S. M. Wittig, J. J. Longdell, and M. J. Sellars, Optically addressable nuclear spins in a solid with a six-hour coherence time, *Nature* **517**, 177 (2015).
[28] Y. Ma, Y.-Z. Ma, Z.-Q. Zhou, C.-F. Li, and G.-C. Guo, One-hour coherent optical storage in an atomic frequency comb memory, *Nature Communications* **12**, 2381 (2021).

[29] M. Businger, L. Nicolas, T. S. Mejia, A. Ferrier, P. Goldner, and M. Afzelius, Non-classical correlations over 1250 modes between telecom photons and 979-nm photons stored in 171yb3+:y2sio5, *Nature Communications* **13**, 6438 (2022).

[30] G. Heinze, C. Hubrich, and T. Halfmann, Stopped light and image storage by electromagnetically induced transparency up to the regime of one minute, *Phys. Rev. Lett.* **71**, 033601 (2013).

[31] M. Zhong, *Development of Persistent Quantum Memories*, PhD Thesis, Laser Physics Centre, Research School of Physics and Engineering, The Australian National University (2017).

[32] A. Ortu, A. Holzäpfel, J. Etesse, and M. Afzelius, Storage of photonic time-bin qubits for up to 20 ms in a rare-earth doped crystal, *npj Quantum Information* **8**, 29 (2022).

[33] A. Holzäpfel, J. Etesse, K. T. Kaczmarek, A. Tiranov, N. Gisin, and M. Afzelius, Optical storage for 0.53 s in a solid-state atomic frequency comb memory using dynamical decoupling, *New Journal of Physics* **22**, 063009 (2020).

[34] R. W. Equall, Y. Sun, R. L. Cone, and R. M. Macfarlane, Ultraslow optical dephasing in Eu³⁺:Y₂SiO₅, *Phys. Rev. Lett.* **72**, 2179 (1994).

[35] L. Nicolas, M. Businger, T. Sanchez Mejia, A. Tiranov, T. Chanelière, E. Lafitte-Houssat, A. Ferrier, P. Goldner, and M. Afzelius, Coherent optical-microwave interface for manipulation of low-field electronic clock transitions in 171Yb3+:Y2SiO5, *npj Quantum Information* **9**, 21 (2023).

[36] E. Fraval, M. Sellars, A. Morrison, and A. Ferris, Pr-Y interaction in Pr³⁺:Y₂SiO₅, *J. Lumin.* **107**, 347 (2004).

[37] J. R. Klauder and P. W. Anderson, Spectral diffusion decay in spin resonance experiments, *Phys. Rev.* **125**, 912 (1962).

[38] W. B. Mims, Phase Memory in Electron Spin Echoes, Lattice Relaxation Effects in CaWO₄: Er, Ce, Mn, *Phys. Rev.* **168**, 370 (1968).

[39] W. Yang and R.-B. Liu, Quantum many-body theory of qubit decoherence in a finite-size spin bath. ii. ensemble dynamics, *Phys. Rev. B* **79**, 115320 (2009).

[40] Y.-Z. Ma, Y.-C. Lv, T.-S. Yang, Y. Ma, Z.-Q. Zhou, C.-F. Li, and G.-C. Guo, Monte carlo simulation of the nuclear spin decoherence process in eu³⁺:y2sio5 crystals, *Phys. Rev. B* **107**, 014310 (2023).

[41] B. Car, J.-L. Le Gouët, and T. Chanelière, Superhyperfine induced photon-echo collapse of erbium in y₂SiO₅, *Phys. Rev. B* **102**, 115119 (2020).

[42] B. Car, L. Veissier, A. Louchet-Chauvet, J.-L. Le Gouët, and T. Chanelière, Optical study of the anisotropic erbium spin flip-flop dynamics, *Phys. Rev. B* **100**, 165107 (2019).

[43] W. M. Witzel and S. Das Sarma, Quantum theory for electron spin decoherence induced by nuclear spin dynamics in semiconductor quantum computer architectures: Spectral diffusion of localized electron spins in the nuclear solid-state environment, *Phys. Rev. B* **74**, 035322 (2006).

[44] W. Yang and R.-B. Liu, Quantum many-body theory of qubit decoherence in a finite-size spin bath, *Phys. Rev. B* **78**, 085315 (2008).

[45] E. Fraval, M. J. Sellars, and J. J. Longdell, Method of extending hyperfine coherence times in pr³⁺:y2sio5, *Phys. Rev. Lett.* **92**, 077601 (2004).

[46] M. Afzelius, C. Simon, H. de Riedmatten, and N. Gisin, Multimode quantum memory based on atomic frequency combs, *Phys. Rev. A* **79**, 052329 (2009).

[47] B. Car, Selective optical addressing of nuclear spins through superhyperfine interaction in rare-earth doped solids, *Physical Review Letters* **120** (2018).

[48] J. Etesse, A. Holzäpfel, A. Ortu, and M. Afzelius, Optical and spin manipulation of non-kramers rare-earth ions in a weak magnetic field for quantum memory applications, *Phys. Rev. A* **103**, 022618 (2021).

[49] C. Li, C. Wyon, and R. Moncorge, Spectroscopic properties and fluorescence dynamics of er/sup 3+/ and yb/sup 3+/ in y/sub 2/sio/sub 5/, *IEEE Journal of Quantum Electronics* **28**, 1209 (1992).

[50] I. Kurkin and K. Chernov, Epr and spin-lattice relaxation of rare-earth activated centres in y2sio5 single crystals, *Physica B+C* **101**, 233 (1980).

[51] E. Z. Cruzeiro, J. Etesse, A. Tiranov, P.-A. Bourdel, F. Fröwis, P. Goldner, N. Gisin, and M. Afzelius, Characterization of the hyperfine interaction of the excited ⁵D₀ state of eu³⁺: y₂sio₅, *Phys. Rev. B* **97**, 094416 (2018).

[52] A. Tiranov, A. Ortu, S. Welinski, A. Ferrier, P. Goldner, N. Gisin, and M. Afzelius, Spectroscopic study of hyperfine properties in ¹⁷¹Yb³⁺: y₂sio₅, *Phys. Rev. B* **98**, 195110 (2018).

[53] R. W. Equall, R. L. Cone, and R. M. Macfarlane, Homogeneous broadening and hyperfine structure of optical transitions in Pr³⁺:Y₂SiO₅, *Phys. Rev. B* **52**, 3963 (1995).

[54] F. Könz, Y. Sun, C. W. Thiel, R. L. Cone, R. W. Equall, R. L. Hutcheson, and R. M. Macfarlane, Temperature and concentration dependence of optical dephasing, spectral-hole lifetime, and anisotropic absorption in Eu³⁺:Y₂SiO₅, *Phys. Rev. B* **68**, 085109 (2003).

[55] A. Abragam, *Principles of Nuclear Magnetism* (Oxford University press, 1961).

[56] B. Herzog, Transient nuclear induction and double nuclear resonance in solids, *Physical Review* **103**, 148 (1956).

[57] H. Syed, A. Kinos, C. Shi, L. Rippe, and S. Kröll, Microscopic model of spin flip-flop processes in crystals doped by rare-earth ions, *Phys. Rev. B* **106**, 184109 (2022).

[58] B. Gong, T. Tu, Z.-Q. Zhou, X.-Y. Zhu, C.-F. Li, and G.-C. Guo, Environment spectrum and coherence behaviours in a rare-earth doped crystal for quantum memory, *Scientific Reports* **7**, 18030 (2017).

[59] W. Yang, W.-L. Ma, and R.-B. Liu, Quantum many-body theory for electron spin decoherence in nanoscale nuclear spin baths, *Reports on Progress in Physics* **80**, 016001 (2016).

[60] D. Rossini, P. Facchi, R. Fazio, G. Florio, D. A. Lidar, S. Pascazio, F. Plastina, and P. Zanardi, Bang-bang control of a qubit coupled to a quantum critical spin bath, *Phys. Rev. A* **77**, 052112 (2008).

[61] J. Mlynek, N. C. Wong, R. G. DeVoe, E. S. Kintzer, and R. G. Brewer, Raman heterodyne detection of nuclear magnetic resonance, *Phys. Rev. Lett.* **50**, 993 (1983).

[62] J. c. v. Minář, N. Sangouard, M. Afzelius, H. de Riedmatten, and N. Gisin, Spin-wave storage using chirped control fields in atomic frequency comb-based quantum memory, *Phys. Rev. A* **82**, 042309 (2010).
[63] I. Roos and K. Mølmer, Quantum computing with an inhomogeneously broadened ensemble of ions: Suppression of errors from detuning variations by specially adapted pulses and coherent population trapping, Phys. Rev. A 69, 022321 (2004).

[64] A. Ortu, A. Tiranov, S. Welinski, F. Fröwis, N. Gisin, A. Ferrier, P. Goldner, and M. Afzelius, Simultaneous coherence enhancement of optical and microwave transitions in solid-state electronic spins, Nature Materials 17, 671 (2018).

[65] E. Fraval, M. J. Sellars, and J. J. Longdell, Method of extending hyperfine coherence times in $^{3+}$:y$_2$si$_5$, Phys. Rev. Lett. 92, 077601 (2004).

[66] M. F. Pascual-Winter, R.-C. Tongning, T. Chanelière, and J.-L. Le Gouët, Spin coherence lifetime extension in $^{3+}$:yag through dynamical decoupling, Phys. Rev. B 86, 184301 (2012).

# Appendix A: Coherence time automated evaluation

To infer a coherence time from our simulations, we automatically detect the end of the initial decay by choosing an arbitrary threshold condition for which we consider that the coherence is lost. As exposed in Sec. III A and Sec. III C, the non-Kramers ions coherence decay curves display highly contrasted modulations at the RE ion Zeeman frequency that cause the coherence to drop to zero periodically. However, such strong modulation cannot be seen as a loss of coherence due to the very strong revival of the coherence. With that consideration, using only a threshold condition on these curves is not sufficient to obtain a relevant characteristic decay time. We instead are interested on detecting the time at which the echo envelope goes below our threshold. Consequently, for non-Kramers ions we first apply a low-pass filter on the echo decay to get rid of the RE Zeeman frequencies if they are present. This is illustrated here in the case of Pr:YSO in Fig. 11.

Once the signal envelope is obtained, we discard the first 100 µs and last 15 % of the signal, as the filtered signal in these time ranges frequently displays significant deviations from the unfiltered signal and causes errors on the detection (the grey areas on the figure 11 are the discarded parts of the echo). This is a direct consequence of the filtering function and unrelated to our simulations. To illustrate this time cutoff, for a 8 ms simulation the minimum decay time is 100 µs, and the maximum detected decay time is 6.5 ms.

Unfortunately, despite these precautions the simulation results may not resemble to an echo decay from which a coherence time could be deduced. For these curves, a coherence time of 0 ms is attributed by default (see e.g. zones with 0 coherence time in Fig. 8a).

However, for all valid decays, we choose a threshold condition for which we consider our coherence lost. The threshold value is an arbitrary choice based on visual evidence to have the most accurate decay time for the largest variety of cases. The chosen value for our threshold condition is an echo intensity going below 9 %.

# Appendix B: Fermi Golden rule, mean dipole and CCE computations details

## 1. Fermi Golden rule

Instead of computing exactly the coefficients in Eq. (17), the Fermi golden rule provides a faster and simpler method to derive them and express a decay rate. Several assumptions are needed in order to unfold the calculations. First, we assume that only one level is populated at the initial time. Secondly, the perturbation induced by the Hamiltonian $\hat{H}_{\text{RE-Y}}^{\text{dd}}$ must be weak compared to $\hat{H}_0$. The probability for $\hat{H}_{\text{RE-Y}}^{\text{dd}}$ to change the
<figure><img src="image_12.png" /><figcaption>FIG. 11: Comparison between the unfiltered coherence decay (in blue) in Pr:YSO for 0.7 mT, $\varphi = 2.95$, $\theta = 2.8$ and 3Y with the same echo after filtering out the Pr Zeeman frequency (in red). The gray areas around 0 and 7 ms are excluded from the decay time detection procedure. The black cross points out the end of the echo decay detected with the chosen threshold condition ($y = 0.09$).</figcaption></figure>

state of the system from the eigenstate $|\varphi_i\rangle_g$ to $|\varphi_f\rangle_e$ is:

$$ \mathcal{P}_{if} = |\langle\varphi_f|\hat{H}_{\text{RE-Y}}^{\text{dd}}|\varphi_i\rangle|^2, \quad (B1) $$

which develops to

$$ \mathcal{P}_{if} = 4(\hat{H}_{\text{int}})_{if}^2 \left| \frac{\sin\left[\frac{(\omega_f - \omega_i)t}{2\hbar}\right]}{(\omega_f - \omega_i)} \right|^2, \quad (B2) $$

with $\omega_f$ and $\omega_i$ the eigenfrequencies of the final and the initial state.

If the state $|\varphi_i\rangle$ is connected to $n$ levels, the total probability for the system to experience a state change is then:

$$ \mathcal{P} = \sum_{k=1}^{n} \mathcal{P}_{ik}. \quad (B3) $$

If we now make the assumption that all the $f$ levels are spaced evenly, such that $\omega_f - \omega_i = k\delta$, and if we consider that the level $i$ is connected with a similar strength ($\hat{H}_{\text{RE-Y}}^{\text{dd}}$)$_{if} = \nu$ to all $f$ levels, (B2) and (B3) give:

$$ \mathcal{P} = 4\nu^2 \sum_k \frac{\sin^2\left[\frac{k\delta t}{2\hbar}\right]^2}{(k\delta)}. \quad (B4) $$

Given that:

$$ \lim_{a \to 0} \sum_{k=-\infty}^{k=+\infty} \frac{\sin^2(ak)^2}{k} = \pi a, \quad (B5) $$

we obtain the following change of state probability per unit of time:

$$ \mathcal{P} = 2\pi \frac{\nu^2}{\delta}. \quad (B6) $$

This result is known as Fermi's golden rule.

In our case, several of the assumptions needed to derive Eq. (B6) do not hold. First, the coefficients of $\hat{H}_{\text{RE-Y}}^{\text{dd}}$ cannot be considered as constants for all pairs of levels. The RE-Y dipole-dipole interactions that generates them is strongly dependent on the Y ions distance and orientation, causing levels associated with a spin change of the closest Y ion to be strongly connected compared to farther ions. On a second hand, the energy levels cannot be considered as evenly spaced. This is directly illustrated with Fig. 10, that shows the Gaussian-like distribution of the Y manifold.

## 2. Frozen central ion models

### a. Model

In this paragraph, we detail the link existing between two models that can be found in the literature: the mean dipole moment model [41, 42] and the CCE model [39, 43, 44]. We also illustrate why they cannot be applied to our simulation cases.

Among others, the CCE method helps to estimate the time dependence of arbitrary echo sequences. In the case of a Hahn echo, it can be written of the form [43]:

$$ \mathcal{E}(t) = \text{Tr}(\hat{\rho}_0 e^{-i\hat{H}^{-}t} e^{-i\hat{H}^{+}t} e^{i\hat{H}^{-}t} e^{i\hat{H}^{+}t}), \quad (B7) $$

where $\hat{H}^{-}$ (resp. $\hat{H}^{+}$) is the Hamiltonian of the bath (Y spins in our case) when the central spin is in its ground (resp. excited) state, and $\hat{\rho}_0$ is the initial spin bath state. The dimension of the bath Hamiltonian makes the calculations quickly intractable if the number of ions is higher than a few. There comes the CCE into play, which consists in expressing the total echo as a product of contribution of individual clusters. In particular, at order 1, each cluster contains a single Y bath ion, and the total echo reads:

$$ \mathcal{E}(t) = \tilde{\mathcal{E}}_{\{1\}}(t)\tilde{\mathcal{E}}_{\{2\}}(t)\cdots\tilde{\mathcal{E}}_{\{N\}}(t), \quad (B8) $$

where

$$ \tilde{\mathcal{E}}_{\{k\}}(t) = \text{Tr} \left( (\hat{\rho}_k)_0 e^{-i\hat{H}_k^{-}t} e^{-i\hat{H}_k^{+}t} e^{i\hat{H}_k^{-}t} e^{i\hat{H}_k^{+}t} \right) \quad (B9) $$

is the echo computed with the sole yttrium ion labeled $k$ in the bath and $(\hat{\rho}_k)_0$ its initial state. The Hamiltonians that shall be considered in expression (B9) are then

$$ \hat{H}_{k}^{\pm} = \hat{H}_{Y_{k}}^{Z} + \hat{H}_{\text{RE}^{\pm}-Y_{k}}^{\text{dd}} \quad (B10) $$

Here, the bath Hamiltonian is solely given by $\hat{H}_{Y_k}^Z$, as at order 1 only single bath ions must be considered. The interaction Hamiltonian (9) now reduces to the dipole-dipole coupling with yttrium ion labeled $k$, for the central ion in its ground ($\hat{H}_{\text{RE}-Y_k}^{\text{dd}}$) or excited ($\hat{H}_{\text{RE}+Y_k}^{\text{dd}}$)
<figure><img src="image_13.png" /><figcaption>FIG. 12: Eigenvalues $\lambda_l$ of the total Hamiltonian (8) as a function of the magnetic field amplitude, for one europium ion and the first two yttrium neighboring ions, for a field oriented along $D_1$, for a) $|\pm 1/2\rangle_g$, b) $|\pm 3/2\rangle_g$ and c) $|\pm 5/2\rangle_g$. The two quickly separating branches stem from the Eu Zeeman contribution (separating $|+k/2\rangle$ from $|-k/2\rangle$), while the four-level manifold contain the two Y contributions. The mean values are subtracted for each hyperfine level for clarity.</figcaption></figure>

state. These Hamiltonians are simply the restriction to the subspaces of central spin down or up of the magnetic dipole-dipole coupling, and are given by

$$ \hat{H}_{\text{RE}\pm-Y_k}^{\text{dd}} = \langle \pm | \hat{H}_{\text{RE}-Y_k}^{\text{dd}} | \pm \rangle = -\gamma_Y \cdot \mathbf{I} \cdot \mathbf{B}_{\text{RE}\pm}^{\{k\}} \quad (\text{B11}) $$

This turns to selecting only the diagonal blocks of the dipole-dipole interaction, in agreement with the hypothesis of a frozen central spin. Notice that this hypothesis is strictly equivalent to taking the mean value of the central spin dipole moment in the interaction Hamiltonian, allowing to define an equivalent magnetic field $\mathbf{B}_{\text{RE}\pm}^{\{k\}}$ radiated by the central ion at the position of the yttrium ion labeled $k$, as done in [41, 47]. With this in mind, we can re-write (B10) as:

$$ \hat{H}_{k}^{\pm} = -\gamma_{Y} \cdot \mathbf{I} \cdot \mathbf{B}_{\text{DC}} - \gamma_{Y} \cdot \mathbf{I} \cdot \mathbf{B}_{\text{RE}\pm}^{\{k\}} := -\gamma_{Y} \cdot \mathbf{I} \cdot \mathbf{B}_{\text{eff}\pm}^{\{k\}} \quad (\text{B12}) $$

By injecting this expression in (B9) and taking initially perfectly mixed Y states $(\hat{\rho}_k)_0 = 1/2|-\rangle\langle-|+1/2|+\rangle\langle+|$, we find the expression

$$ \tilde{\mathcal{E}}_{\{k\}}(t) = 1 - \frac{\sin^2(\theta_k)}{2} [1 - \cos(\Delta_k^- t)] [1 - \cos(\Delta_k^+ t)], \quad (\text{B13}) $$

with $\theta_k$ the angle between the two effective magnetic fields

$$ \sin^2(\theta_k) = 1 - \mathbf{B}_{\text{eff}+}^{\{k\}} \cdot \mathbf{B}_{\text{eff}-}^{\{k\}} \quad (\text{B14}) $$

and frequencies defined as

$$ \Delta_k^{\pm} = \gamma_Y ||\mathbf{B}_{\text{eff}\pm}^{\{k\}}||. \quad (\text{B15}) $$

Expression (B13) coincides with the one that was used in [41, 42] in the mean field model. Finally, the first order CCE of the echo is computed with (B8).

## b. Invalidity of the approximation in the weak field regime

In order to assess the validity of the frozen ion approximation, we plot the eigenvalues $\lambda_l$ of $\hat{H}$ (8) in the case of Eu:YSO for the three hyperfine ground states, as a function of the magnetic field amplitude, for a field oriented along the $D_1$ axis. Such eigenvalues are the ones dictating the dynamics of the evolution on the energetic ladder defined by $\hat{H}_0$. We remind here that the dynamical frequencies are given by the intervals $\Delta\lambda_{l,l'} = \lambda_l - \lambda_{l'}$. We also plot the same eigenvalues calculated with the mean dipole hypothesis in dashed black lines for assessing its validity. What can be directly witnessed is that in the magnetic field amplitude regime considered here, the frozen spin approximation is not always relevant, as strong energy mismatch (up to ~100 Hz for transitions of equivalent energy gap) are witnessed. Another important observation is the very different behavior of the hyperfine levels: while $|\pm 1/2\rangle_g$ actively spin-flips with the neighboring Y ions, $|\pm 5/2\rangle_g$ remains quite protected from this phenomenon, even at fields below 10 µT.

# Appendix C: Population dynamics in Eu

We show in this appendix the population dynamics of level $|\pm 3/2\rangle_g$ in Eu:YSO as a function of the magnetic field amplitude, to highlight spin flip processes at stake at different field amplitudes. To this extent we consider a Eu ion in interaction with the first two neighboring Y ions, and place the system initially in an eigenstate of $\hat{H}_0$. Therefore, on the contrary to the sequence used to simulate coherence evolution over time in the body of the paper, no RF pulse is applied to create superpositions between the hyperfine states. We then let the system to evolve under dipole-dipole interaction only and plot the population of each level versus time. This is shown in Fig. 13 for a field applied along $D_1$, with a magnitude ranging from 1 to 60 µT. The population is shown with
<figure><img src="image_14.png" /><figcaption>FIG. 13: Population dynamics of the $|\pm 3/2\rangle_g$ level of Eu:YSO as a function of the magnetic field amplitude (field along $D_1$).</figcaption></figure>

a colorscale, and reveals that at fields lower than 10 µT essentially all states end up being populated over time. Above 10 µT, we see that the Eu spin flip begins to be forbidden, and mostly Y spin flips and Y flip-flops remain. Then, above 40 µT, only two levels are populated, corresponding to Y-Y flip flops. This observation corroborates the observations done in coherence with simple population simulations.

## Appendix D: Coherence map of Yb

In order to get a general picture of the coherence dependence of Yb:YSO, we plot in Fig. 14a the echo decay times (see App. A for its estimation), as a function of the magnetic field direction and amplitude for one Yb ion and the first three neighboring Y ions.

The magnetic dipole of Yb scales linearly with the magnetic field amplitude, therefore Yb has a zero-field ZEFOZ point, where the coupling with the surrounding Y ions is minimal [64]. From that point, the coherence time

gets longer than the simulated evolution time (10 ms), causing the automated decay time detection to fail (black region on the plot). Past that zone, the simulated decay time decreases at higher magnetic field amplitudes because of the increased Yb-Y couplings. We also notice that the coherence time is significantly improved along a specific direction. It corresponds to the direction for which the first order transition gradient $||S_1||$ is minimal in the $(D_1, D_2)$ plane, at $\phi = 55.9^\circ$. This has also been witnessed experimentally [35].

## Appendix E: Relevance of the sequence considered

Here, we compare the simulated decay for europium at 60 µT (case of Fig. 5a) with a decay generated for a full simulation of a Hahn echo sequence. Both simulations are shown in Fig. 14b, and even if slight deviations from one curve to the other can be seen, the overall shape, modulation pattern and characteristic decay time are similar. This curve comforts us in the relevance of our approach, even at such field, for which single ion inhomogeneous decoherence is close to be limiting our decay time.
<figure><img src="image_15.png" /><figcaption>FIG. 14: a) Coherence time in millisecond as a function of magnetic field angle, for one Yb ion and 3 Y ions. b) Decays simulated with the simple FID sequence (used in the paper, solid blue), and with a Hahn sequence (orange dots), for 3 Y ions.</figcaption></figure>

<figure><img src="image_16.png" /></figure>