Hyperfine spectroscopy in a quantum-limited spectrometer

Abstract We report measurements of electron-spin-echo envelope modulation (ESEEM) performed at millikelvin temperatures in a custom-built high-sensitivity spectrometer based on superconducting micro-resonators. The high quality factor and small mode volume (down to 0.2 pL) of the resonator allow us to probe a small number of spins, down to 5×102 . We measure two-pulse ESEEM on two systems: erbium ions coupled to 183W nuclei in a natural-abundance CaWO4 crystal and bismuth donors coupled to residual 29Si nuclei in a silicon substrate that was isotopically enriched in the 28Si isotope. We also measure three- and five-pulse ESEEM for the bismuth donors in silicon. Quantitative agreement is obtained for both the hyperfine coupling strength of proximal nuclei and the nuclear-spin concentration.

Section S1 -Details of the Model for ESEEM of Bi:Si In this Section, we give the details of the Si:Bi Hamiltonian, as a supplementary note to Section III B of the main text.
The Hamiltonian of the bismuth donor spins and the 29 Si nuclear spins is 5 where H Bi is the Hamiltonian for the Bismuth electron and nuclear spin, H Si is the Zeeman energy and the dipole-dipole interaction of the 29 Si nuclear spins, and H hf is the hyperfine interaction between the Bismuth electron spin and the 29 Si nuclear spins. The Hamiltonian of the Bi center spins is (Ma et al., 2015) H Bi = g e β e S 0,z + A Bi S 0 · I 0 , (S2) 10 where the notations are the same as in main text. The electron gyromagnetic ratio g e β e = 1.76 × 10 11 S −1 T −1 and the Fermi contact hyperfine coupling A Bi /2π = 1.4754 GHz. Note that the Zeeman energy of the nuclear spin is neglected (as explained in the main text). The internal Hamiltonian of the 29 Si nuclear spin bath is where ω I = g n β n B 0 denotes the Larmor frequency of the 29 Si nuclear spin with the gyromagnetic ratio g n β n = −5.319 × 10 7 s −1 T −1 , the spin operator I j denotes the nuclear spin at position r j , and r ij = r i − r j .

15
The hyperfine interaction between the electron spin and the 29 Si nuclear spins is where the hyperfine coupling tensor includes both the Fermi contact term A cf,j (a scalar) and the dipolar interactionĀ dd,j , i.e.,Ā j = A cf,j 1 +Ā dd,j . The dipolar hyperfine coupling tensor is A dd,j = µ 0 g e β e g n β n 4πr 3 j 1 − 3r j r j r 2 j ≡ A dd,j 1 − 3r j r j r 2 j , decaying with a cubic power of the distance. The Fermi contact interaction is proportional to the electron density |ψ(r j )| 2 of the Bismuth donor at the nuclear spin position r j (Feher, 1959;Kohn, 1957;de Sousa and Sarma, 2003), A cf,j = 2 3 µ 0 g e β e g n β n |ψ(r j )| 2 .
The electron density |ψ(r j )| 2 = 2 3 η|f (r j )| 2 [cos(k 0 x j ) + cos(k 0 y j ) + cos(k 0 z j )] 2 , where k 0 = 0.85 2π aSi (with a Si = 0.543 nm) is the wavenumber of the conduction band minimum, η ≈ 180 is the charge density 20 on each site, and the envelope function is taken as the Kohn-Luttinger wave function form (de Sousa and Sarma, 2003) f j (r) = 1 π(sa) 2 (sb) exp − z 2 (sb) 2 + x 2 + y 2 (sa) 2 , with a = 2.51 nm and b = 1.44 nm being the characteristic lengths for hydrogenic impurities in Si and the scaling factor s = 0.64 for bismuth (Hale and Mieher, 1969). The hyperfine coupling strength |Ā j | is mostly << MHz. : Section S2 -Applying the fictitious spin-1/2 model to Bi:Si coupled to 29 Si spins In this section, we provide detailed justifications for applying the fictitious model in Section II D of main text to Bi:Si coupled to 29 Si spins.
The eigenstates |±, m of the Bismuth Hamiltonian H Bi in Eq. (S2) are 5 |±, m = cos with tan θ m = √ 25−m 2 m+(1+δ)geβeB0/ABi , and the corresponding eigenenergies are For the experimental condition |g e β e B 0 | |A Bi |, the eigenenergies in Eq. (S5) can be approximated as and the level splitting For a field about 1 Gauss (a typical value in our experiments), (2π) −1 g e β e B 0 /10 ∼ 300 kHz. The hyperfine coupling to the 29 Si nuclear spins can be written in the basis of the eigenstates of H Bi . The non-vanishing 10 matrix elements of the electron spin operators are The diagonal matrix elements cause the frequency shift due to the coupling to the nuclear spins,which is ∼ |Ā j |/2, and the 15 off-diagonal elements cause the coupling between neighboring energy levels ∼ |Ā j |/4. For a nuclear spin with relatively weak hyperfine coupling, namely, the hyperfine-induced mixing between different energy levels of the Bismuth center can be neglected. Then the interaction Hamiltonian becomes where as defined in Eq. (S6a). This Hamiltonian is diagonal in the Bismuth center eigenstates, corresponding to the pure dephasing model (or secular approximation).

25
As will be demonstrated in Sec. , the contribution of a relatively strongly coupled 29 Si nuclear spin (|Ā j |/2 |g e β e B 0 |/10) to the ESEEM signal in our experiments is negligible. Thus the pure dephasing model in Eq. (S8) is justified.
The microwave pulses in general can induce many transitions in the Bi center as along as the transitions are within the bandwidth of the microwave cavity (164 kHz for the high-Q resonator). Since the interference between the transitions whose frequencies are not near-degenerate would cause oscillation much faster than the time resolution of the experiments and the transitions with near-degenerate frequencies do not share an eigenstate of the Bi center (and therefore has no interference), we can reduce the Bi center to independent two-level systems ("fictitious spin-1/2") coupled to the Si spin bath.
The "fictitious spin-1/2" for the transition |+, m ↔ |−, m − 1 has the Hamiltonian where the Overhauser field along the z direction is The transition frequency or effective Larmor frequency of the fictitious spin-1/2 Using a fictitious spin-1/2 to represent the transition |−, m − 1 ↔ |+, m and define δ m ≡ (α m − α m−1 )/2 andᾱ m ≡ (α m + α m−1 )/2, we obtain the Hamiltonian For a single nuclear spin (denoted as spin-j), the conditional Hamiltonian is Note that we have dropped the ω S -term (by working in the rotating reference frame) and for convenience chosen the coordinate 15 system such that A j,zy = 0. In the basis of the eigenstates ofH m , the effective Hamiltonian can be written as wherẽ with θ = arcsin (δ m A j,zx /2ω I ).
Section S3 -Analytical solutions of ESEEM envelopes for the fictitious spin-1/2 model In this Section, we provide details on the analytical solutions for the 2p-, 3p-and 5p-ESEEM (Kasumaj and Stoll, 2008), as a supplementary note on how the curves in Figs. 9-11 of Sec. VB in main text are calculated.

25
: Section S3.1 -Formula In the n-pulse ESEEM experiment, a sequence of n pulses R j ∈ {R x π/2 , R y π/2 , R x π , R y π } are applied at t j (j = 1, 2, . . . n), and then the signal V np is measured at the echo time t. In this Section, we take the pulses as ideal, i.e., R x/y θ = exp −iθS x/y . With the pulse sequences described in Sec. VB of main text, the evolution for the 2p-, 3p-, and 5p-ESEEM is in turn 5 We assume the Bi center is initially in the state ρ S = |−, m − 1 −, m − 1| and the nuclear spin bath is in the maximally mixed state ρ B . The spin coherence at the echo time is Section S3.2 -Exact formula for a single nuclear spin For 2p-ESEEM, the modulation amplitude due to the j-th 29 Si spin is where , and (the same as Eq. 6 of main text)

15
For the 3p-ESEEM, the modulation amplitude is given by where and V ↓ 3p,j (τ, T ) is obtained by exchanging ↑ and ↓.

Section S3.3 -Contributions of multiple nuclear spins
For multiple nuclear spins, if they are taken as independent (with interactions between the nuclear spins neglected), the ESEEM signal is obtained by applying the product rule for each pathway followed by average over different pathways ((Kasumaj and Stoll, 2008)). For the 2p/3p/5p-ESEEM signals are in turn Above we have considered the ESEEM signal due to the transition |−, m−1 ↔ |+, m . The transitions |−, m ↔ |+, m−1 can be considered similarly. We use V

Section S3.4 -Weighting factors of different Bi spin transitions
To take into account the contributions of all possible transitions, we assume the Bi center spin is initially randomly populated in the lower manifold of the eigenstates, i.e., with the probabilities m P m = 1. The coherence of different transitions is summed up as with the weighting factor W m,± accounting for the initial probabilities P m and the different amplitudes for different transitions 20 given the microwave pulse spectra and the distribution of the hyperfine coupling A Bi due to the strain in the Si layer.
To determine the relative transition weights, we resort to a complete model of the experiment, based on the physical parameters, namely the control pulse amplitude and temporal profile, the Rabi frequencies distribution, and the repetition time Γ −1 rep . We use a simulation code that was purposely written. For each class of fictitious spin-1/2's with given Larmor frequency ω S and Rabi frequency Ω R , the program integrates the Bloch equations to compute the time dependence of the spin density 25 matrix ρ(t). The initial conditions take into account the Rabi-frequency-dependent spin relaxation T 1 because of the Purcell effect, by taking ρ(0) = |x x| (the spin polarized along the x-axis). The simulation results are furthermore averaged over the strain-induced distribution of the 209 Bi hyperfine coupling σ A (A Bi ), which results in the distribution of the transition frequency σ S (ω S ), and over the distribution of the Rabi frequency σ R (Ω R ). The Bi hyperfine coupling distribution σ A (A Bi ) is taken to be flat since the inhomogeneous broadening (∼50 MHz) is two orders of magnitude larger than the cavity bandwidth (∼160 kHz). 30 -1 1 1 10 10 -1 10 -4 100 1000   The Rabi frequency distribution σ R (Ω R ) is computed using finite-element modelling of the AC field spatial profile generated by running a constant current through the resonator inductance wire. The AC field map and the resulting distribution σ R (Ω R ) are shown in Fig. S1. To calibrate the AC pulse amplitude, we rely on the fact that due to the AC field inhomogeneity and to the Purcell spin relaxation, the measured T 1 in an inversion recovery sequence is in fact amplitude-dependent. We then adjust the pulse amplitude in the simulation so that the simulated inversion recovery sequence reproduces the same relaxation curve 5 as measured. The weighting factors thus derived from simulating the experimental data are given in Table. S1.

Section S3.5 -Comparison with experiments
For comparison with experimental data, the curves in Figs. 9-11 of main text are obtained by first calculating V 2p/3p/5p using Eq. (S22) for each nuclear spin configuration (with random positions of the 29 Si nuclei) and then averaging over different nuclear spin spatial configurations.
10 Section S4 -Justification of independent bath spin approximation In the calculations above, the nuclear spins in the bath are taken as independent, i.e., the interactions between the 29 Si spins are neglected. In this Section, we will provide a justification of the approximation by calculating the Bi spin coherence with the cluster correlation expansion (CCE) (Yang andLiu, 2008, 2009;Zhao et al., 2012). The leading order of expansion (CCE-1) corresponds to the independent bath spin approximation used in the previous Section. The numerical calculation shows that 15 the CCE-1 is a good approximation for the timescale relevant to the experimental studies of the ESEEM.
We check the contributions of various orders of the nuclear spin correlations to the ESEEM signals using the fictitious spin-1/2 model in Eq. (S10) and the CCE method (Yang andLiu, 2008, 2009;Zhao et al., 2012). The interactions between the 29 Si nuclear spins are included in H Si . The CCE method allows us to consider in a recursive way the dynamics of the nuclear spin bath due to correlations of different sizes (CCE-n for n-spin irreducible correlation).
. CCE calculation of the Hahn echo signal of a Bi donor spin coupled to a 29 Si nuclear spin bath. The CCE-2 and CCE-3 results have excluded the CCE-1 effects. (a) Decoherence for two different random bath configurations (denoted as "A" and "B"), calculated up to CCE-2 or CCE-3 (both with CCE-1 contributions excluded). B0 = 1 G. The number of bath spins N = 2000. The negligible difference between CCE-2 and CCE-3 indicates that the CCE-2 has already converged. (b) Comparison of the CCE results with the experimental ESEEM signal for B = 1 G. For the "CCE-1 w/ decay", an overall exponential decay is included as explained in the main text.
Considering the 2p-ESEEM for example, the spin coherence is expanded as withṼ C 2p defined as the irreduible correlation of cluster C excluding the irreducible correlations of all sub-clusters, i.e., where V C 2p is the center spin coherence under coupling to the cluster C in the spin bath (the bath spins outside the cluster dropped).
As shown in Fig. S2, in the timescale relevant to the ESEEM signals in the experiments, the echo is affected mainly by the single-spin dynamics (CCE-1) and the contributions from pair dynamics (CCE-2) and higher order correlations in the bath are 10 negligible. The decoherence due to CCE-2 and CCE-3 occurs at timescales of ∼ 100 ms [ Fig. S2(a)], which are much longer than the experimental time regime (∼ 1 ms).
Thus it is well justified to take the nuclear spins in the bath as independent of each other.
Section S5 -Effect of a strongly coupled Si-29 nuclear spin In this Section, we show that the strongly coupled 29 Si nuclear spins have negligible contributions to the ESEEM signals. 15 Therefore the calculated ESEEM signals presented in Figs. 9-11 in the main text are those from 29 Si with hyperfine couplings < 20 kHz.
In the calculations, we have assumed the pure dephasing model (the secular approximation) in which the transitions between different eigenstates of the Bi center spin due to the hyperfine coupling to the 29 nuclear spins are taken as negligible. This approximation is well justified if the hyperfine coupling is much less than the energy splitting between different Bi center spin 20 states [Eq. (S7)]. When the coupling is strong the secular and the CCE-1 approximation become invalid. Using simulations that take into account exactly the effect of the strongly coupled nuclear spin, we show that a strongly coupled spin has negligible effects on the ESEEM signal in the timescale considered in the experiments. A strongly coupled spin contributes only fast oscillations in the signal, which would vanish if we take into account the inhomogeneous broadening effects. In addition, the influence of the strongly coupled spins on the distant weakly coupled spins are also negligible.
For a 29 Si with hyperfine coupling 200 kHz, the state mixing due to the hyperfine coupling enables many transitions and the interference between these transitions will cause rather complicated and fast oscillations in the spin echo signal. This is 5 seen in Figs. S7. However, in such a strong coupling case, the ESEEM frequencies depend sensitively on the local Overhauser field on the Bi electron spin. Since the Overhauser field has a large inhomogeneous broadening (∼ 0.5 MHz), the ensemble average over the nuclear spin state thermal distribution leads to a rapid decay of the signal (decay in < 1 µs). We estimate that about 10% Bi centers have one or more 29 Si with coupling > 200 kHz in the proximity, which would contribute to a fast initial decay of the total echo signal in < 1 µs by about 10%. In the experiments, the echo signal is measured at times much greater 10 than µs (the first time point is ∼ 1 ms). As shown in Figs. S5 and S7, the ESEEM amplitude of a nuclear spin with coupling strength between 20 kHz and 200 kHz is much less than 1%. And for nuclear spins with the relatively weak hyperfine coupling (2π) −1 |Ā j | < 100 kHz, the pure dephasing model produces results with negligible errors of the modulation frequencies from the exact solution (Figs. S5-S6). Furthermore, the systematic numerical studies (Fig. S9) show that a nearby Si nuclear spin with coupling < 200 kHz has negligible effects on the ESEEM due to other distant nuclear spins.

15
Considering these different contributions of Si nuclear spins of different hyperfine couplings, it is justified to assume the pure-dephasing model and consider only the contributions of those Si nuclear spins that have couplings weaker than a certain cut-off (chosen as |Ā j | ≤ 20 kHz).

Section S5.1 -CCE-1 for the multi-level central system
Strong hyperfine interaction would cause the mixing between the Bismuth center eigenstates. To consider the effect of a 20 strongly coupled 29 Si nuclear spin, we take the electron spin S 0 , the 209 Bi nuclear spin I 0 , and the j 0 "strongly coupled" 29 Si nuclear spins (denoted as I j for 1 ≤ j ≤ j 0 ) as a hybrid center spin system. The hybrid center system can be diagonalized, with 2×10×2 j0 eigenstates, separated into two manifold {|±, m }. The corresponding eigenenergies of the hybrid center spin system is denoted as E ± m . See Fig. S3 for an example of the eigenenergies of a hybrid spin system with one strongly coupled 29 Si nuclear spin.  Figure S3. Eigenenergies of a hybrid spin system composed of the Bi donor electron spin S0, the 209 Bi nuclear spin I0, and a "strongly coupled" 29 Si nuclear spin at position r1 = (1.8, 0.4, −0.9) nm. The Fermi-contact hyperfine coupling to the Si spin is A fc,1 /(2π) ≈ 250 kHz. The magnetic field along the z-axis B0 = 1 G. The x axis denotes the state label for the upper and lower eigenstates (the ± manifolds).
The Hamiltonian of the hybrid center system plus the rest N − j 0 weakly coupled 29 Si nuclear spins can be written as where the dipole-dipole interaction between the nuclear spins in the bath is neglected because it is weak and would have effects only on higher order CCE. By assumption, the hyperfine couplings |Ā j | (for j > j 0 ) are much smaller than the energy differences among the eigenstates of the central system (otherwise they would have been absorbed into the hybrid central system). Therefore, we have the pure dephasing model with the bath Hamiltonians conditioned on the states of the central system (for the conciseness, we label (η, m) as k) where the Si nuclear spin bath Hamiltonian H Si excludes the strongly coupled 29 Si nuclear spins I j (for j ≤ j 0 ) and neglects the dipolar interactions between the Si nuclear spins. 5 In the echo experiment, a sequence of n pulses R j ∈ {R x π/2 , R x π , R y π/2 , R y π } are applied at t j (j = 1, 2, . . . n), and then the signal V 3p is measured at time t. The π/2 and π control pulses are taken as R x/y π/2 ≡ e −i(ΩRS 0,x/y +HBi)τp/2 and R x/y π = e −i(2ΩRS 0,x/y +HBi)τp/2 , where τ p is the duration of a π pulse and Ω R the Rabi frequency. Note that the pulse durations are assumed much shorter than the timescales of dynamics of the 29 Si bath spins and thus the rotation transform does not include H Si . If the system starts with a certain state |k 0 ⊗ |J at t = 0 (in which |k 0 denotes a certain eigenstate |−, m and |J a certain bath state), the state at time t is k1,k2,...,kn C k0,k1,...kn e −iφ k 1 ,k 2 ,...kn |k n |J k1,k2,...kn , where |k j denotes an eigenstate |±, m , the phase the coefficients of the eigenstates C k0,k ≡ k n |R n |k n−1 · · · k 2 |R 2 |k 1 k 1 |R 1 |k 0 , and the bath state with the shorthand notation k = (k 1 , k 2 , . . . , k n ). The echo signal is In general, the nuclear spin state overlap J k |J k is computed using the CCE (Yang andLiu, 2008, 2009;Zhao et al., 2012), and in our current case for the Hamiltonian in Eq. (S25), CCE-1 gives the exact solution. The result of V 2p/3p/5p is further averaged over different initial states |k 0 and |J in the thermal distribution.

10
Without loss of generality we consider the 3p-ESEEM experiment with the magnetic field B 0 = 1 G. For the 3p-ESEEM, the CCE-1 result is where U (k1,k2,k3) j is the evolution operator of the spin I j for the center spin pathway k 0 → k 1 → k 2 → k 3 .

15
For the hybrid center that contains only the Bi electron and nuclear spins (no 29 Si spins), the nonzero elements of the S 0,x operator are : Table S2. FFT I 1 I 2 Figure S4. Exact solution of the 3p-ESEEM of a Bi donor spin coupled to only one 29 Si nuclear spin. The Fermi contact and dipolar hyperfine couplings are A cf,1 = 2.4 kHz and A dd,1 = 0.6 kHz in the first case (labeled as I1), and A cf,2 = 0.8 kHz and A dd,2 = 0.5 kHz in the second case (labeled as I2). The external field B0 = 1 G. The left panel shows the ESEEM, and the right panel shows the Fourier transform of the ESEEM. The cutoff thresholds |C * k C k | < 10 −5 and 10 −4 for a pathway to be dropped are indicated in the left panel. The two choices of cutoff thresholds produce nearly identical results. Table. S2.

See numerical values in
Not all of the pathways have contributions to the ESEEM signal because of the inhomogeneous broadening and selection rules. For instance, considering a hybrid center spin system that contains only the Bi electron and nuclear spins, the phase difference accumulated for k j = (+, m) and k j = (+, m ) during time from t j−1 to t j is 5 where h z is the Overhauser field from the bath spins. For the relevant timescales (∼ 1 ms) and inhomogeneous broadening of h z (which is ∼ 0.5 MHz), the ensemble averaged phase factor would vanish unless k j = k j . In the numerical simulation, a pathway (k , k) is dropped if (1) the echo condition is not satisfied, or (2) the amplitude is too small, e.g., |C * k C k | < 10 −4 (see Fig. S4).

Section S5.2 -Exact simulation for one 29 Si spin
Taking the special case of Eq. (S28) for j 0 = 1 and N = 1, we can obtain the exact solution of the ESEEM due to a single 29 Si spin. The result for a certain initial state |k 0 is where φ (k1,k2,k3) ≡ (E k1 + E k3 )τ + E k2 T is the phase accumulated for the pathway k 0 → k 1 → k 2 → k 3 .

5
In numerical simulations we neglect the pathways that have negligible probabilities C * k 0 ,k 1 ,k 2 ,k3 C k0,k1,k2,k3 |. For instance, Fig. S4 compares the results with pathways neglected when |C * k C k | < 10 −5 or 10 −4 . Actually, for these two cases shown in the figure, there are only 4 main contributing pathways, since most pathways have no contributions in the ESEEM signal after taking into account of the inhomogeneous broadening effects (see Eq. S29). We expect that when the hyperfine coupling is comparable to the level splitting of the Bismuth donor spin within each manifold, the secular and the CCE-1 approximation become invalid (Eq. S7). Figure S5 compares the exact solution and the 15 CCE-1 approximation for various hyperfine coupling strengths. As expected, the CCE-1 agrees well with the exact solution for relatively weak coupling. The deviation of the ESEEM frequency calculated by the CCE-1 from the exact solution is shown in Fig. S6 as a function of the hyperfine coupling strength. For coupling weaker than 20 kHz, the error is less than 5%. Furthermore, the ESEEM depth (shown in Fig. S6) becomes negligible for coupling greater than 20 kHz for the field. (It should be noted that if the echo condition in Eq. (S30) is applied, strong hyperfine coupling would lead to mixing between Bi spin 20 states, causing fast oscillations with amplitude increasing with the coupling strength. Such fast oscillations, however, decay rapidly to zero due to inhomogeneous broadening. See Sec. below for more discussions.) Considering both the ESEEM frequency calculation precision and the modulation depth, a nuclear spin with hyperfine coupling weaker than 20 kHz can be well approximated by the CCE-1.
Section S5.4 -ESEEM due to a strongly coupled 29 Si spin 25 Figure. S7 show the exact simulation of the 3p-ESEEM due to a strongly coupled 29 Si spin. Two coupling strengths (∼ 100 kHz and ∼ 200 kHz) are considered. The strong hyperfine coupling induces mixing between the Bi spin states (violating the pure dephasing approximation). To show the effect of state mixing, we do not impose the echo condition Eq. (S30), and therefore need to include about 100 pathways to produce converged results. The interference between different pathways (due to the mixing among Bi spin states) induce fast and complicated modulations. The Fourier transform shows that the ESEEM 30 frequencies are spread around 280 kHz, which is approximately the splitting between Bi levels in absence of Si spins. The fact that modulation frequency is near 280 kHz confirms that the ESEEM is mostly due to the mixing of the Bi spin states. On the contrary, in the case of weakly coupled 29 Si spins where the mixing between different Bi levels is negligible, the ESEEM is due to the mixing between different 29 Si nuclear spin states, so there the ESEEM has frequencies essentially determined by the 29 Si nuclear spin Larmor frequency (see Fig. S5).

35
The rapid oscillations due to interference between different Bi spin states, however, have their frequencies sensitively depending on the external field and the Overhauser field due to the random configurations of the nuclear spin baths (the inhomogeneous broadening). Therefore, the ensemble average over the distribution of the Overhauser field h z (which has a broadening of about 0.5 MHz) would lead to rapid decay of the modulation (in µs timescales). This is indeed shown in Fig. S8. Thus, the contributions of the strongly coupled nuclear spins are not observable in the timescales relevant in the experiments (∼ms).

40
Section S5.5 -Influence of a strongly coupled 29 Si spin on weakly coupled 29 Si spins In general, a strongly coupled nuclear spin may affect the ESEEM due to a weakly coupled nuclear spin. This is because the Bi spin states can be mixed by the strong coupling. The state mixing renormalizes the effective couplings between the weakly coupled nuclear spin and the "ficticious" spin-1/2 of the Bi donor [in particular, α m and α m−1 in Eq. (S13)] and hence affects the ESEEM. To exam such an effect, we calculate the ESEEM due to a weakly coupled nuclear spin I 2 in the presence or in the 45 absence of a strongly coupled nuclear spin I 1 . The weakly coupled nuclear spin is considered using the CCE-1 and the strongly coupled one, if taken into account, is exactly considered by absorbing it into the hybrid center spin system [corresponding to Eq. (S28) with j 0 = 1 and N = 2]. Figure S9 shows that while the strongly coupled spin causes some fast, small-amplitude modulations, it has negligible effect on the ESEEM due to the weakly coupled nuclear spin. The weak effect is understandable considering that the matrix elements α m and α m−1 in Eq. (S13) are only slightly (< 1/10) modified by the hyperfine coupling to the spin I 1 even when A fc,1 is as 5 large as 200 kHz.

Section S6 -Phase cycling scheme
In the tables below is provided the phase cycling scheme for the 3-and 5-pulse ESEEM measurements.  Figure S8. Ensemble average of (a) 3p-ESEEM over 100 samples of nuclear spin states in a thermal distribution and (b) the spectrum. The locations and the hyperfine couplings are r1 = (8, 5, −9) × aSi/4, A cf,1 ≈ 207.1 kHz, and A dd,1 ≈ 2.8 kHz. The external field B0 = 1 G.
All pathways with C * k 0 ,k C k 0 ,k > 10 −4 are taken into account.  Figure S9. The ESEEM of a weakly coupled 29 Si nuclear spin with or without another nuclear spin strongly coupled to the Bi donor spin. The Fermi contact hyperfine coupling A fc,2 ≈ 2.4 kHz for the weakly coupled nuclear spin, and A fc,1 ≈ 45.5 or 207.1 kHz for the strongly coupled nuclear spin. The external field B0 = 1 G. 3 pulse ESEEM CPMG π/2 π/2 π/2 Det. π Det. +x +x +x +y +y +y -x +x +x -y +y -y -x -x +x +y +y +y +x -x +x -y +y -y Table S4. Phase cycling table for 5 pulse ESEEM with CPMG and offset removal.