Hyperfine spectroscopy in a quantum-limited spectrometer

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×10². We measure two-pulse ESEEM on two systems: erbium ions coupled to ¹⁸³W nuclei in a natural-abundance CaWO₄ crystal and bismuth donors coupled to residual ²⁹Si nuclei in a silicon substrate that was isotopically enriched in the ²⁸Si 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.

1 Introduction

Electron paramagnetic resonance (EPR) spectroscopy provides a set of versatile tools to study the magnetic environment of unpaired electron spins (Schweiger and Jeschke, 2001). EPR spectrometers rely on the inductive detection of the spin signal by a three-dimensional microwave resonator tuned to the spin Larmor frequency. While concentration sensitivity is the main concern for dilute samples available in macroscopic volumes (Song et al., 2016), there are also cases in which the absolute spin detection sensitivity matters, motivating research towards alternative detection methods to measure smaller and smaller numbers of spins. Electrical (Elzerman et al., 2004; Veldhorst et al., 2014; Morello et al., 2010; Pla et al., 2012), optical (Wrachtrup et al., 1993; Jelezko et al., 2004), and scanning-probe-based (Rugar et al., 2004; Baumann et al., 2015) detection methods of magnetic resonance have reached sufficient sensitivity to detect individual electron spins.

In parallel, recent results have shown that the inductive detection method can also be pushed to much higher absolute sensitivity than previously achieved, using planar micro-resonators (Narkowicz et al., 2008; Artzi et al., 2015) and micro-helices (Sidabras et al., 2019). Superconducting resonators (Wallace and Silsbee, 1991; Benningshof et al., 2013; Sigillito et al., 2014) are particularly useful in that context since they combine low mode volume and narrow linewidth κ. Inductive-detection spectrometers relying on a superconducting planar micro-resonator combined with a Josephson parametric amplifier (JPA), cooled down to millikelvin temperatures (Bienfait et al., 2015; Eichler et al., 2017; Probst et al., 2017), have achieved a sensitivity of 10 spin$/\sqrt{\mathrm{Hz}}$ for detecting Hahn echoes emitted by donors in silicon (Ranjan et al., 2020b). A particular feature of these quantum-limited spectrometers is that quantum fluctuations of the microwave field play an important role. First, the system output noise is governed by these quantum fluctuations, with negligible thermal noise contribution. Second, quantum fluctuations also impact spin dynamics by triggering spontaneous emission of microwave photons at a rate ${\mathrm{\Gamma }}_{\mathrm{P}}=\mathrm{4}{g}^{\mathrm{2}}/\mathit{\kappa }$, g being the spin–photon coupling (Bienfait et al., 2016; Eichler et al., 2017; Ranjan et al., 2020a). This Purcell effect forbids T₁ from becoming prohibitively long since it is at most equal to ${\mathrm{\Gamma }}_{\mathrm{P}}^{-\mathrm{1}}$, making spin detection with a reasonable repetition rate possible even at the lowest temperatures.

Hahn echoes are the simplest pulse sequence used in EPR spectroscopy, useful to determine the electron-spin density as well as the spin Hamiltonian parameters and their distribution. The richness of EPR comes from the ability to characterize the local magnetic environment of the electron spins, often consisting of a set of nuclear spins or of other electron spins. For that, hyperfine spectroscopy is required, which uses more elaborate pulse sequences and requires larger detection bandwidth. Previous hyperfine spectroscopy measurements with superconducting micro-resonators include the electron–nuclear double resonance detection of donors in silicon (Sigillito et al., 2017) and the electron-spin-echo envelope modulation (ESEEM) of erbium ions by the nuclear spin of yttrium in a Y₂SiO₅ crystal (Probst et al., 2015).

Here, we demonstrate that hyperfine spectroscopy is compatible with quantum-limited EPR spectroscopy despite its additional requirements in terms of pulse complexity and bandwidth, by measuring ESEEM in two model electron-spin systems. We measure the ESEEM of erbium ions coupled to ¹⁸³W nuclei in a scheelite crystal (CaWO₄) with a simple two-pulse sequence and get quantitative agreement with a simple dipolar interaction model. We also measure the ESEEM of bismuth donors in silicon caused by ²⁹Si nuclei using two-, three-, and five-pulse sequences (Schweiger and Jeschke, 2001; Kasumaj and Stoll, 2008). Compared to other ESEEM measurements on donors in silicon (Witzel et al., 2007; Abe et al., 2010), ours are performed in an isotopically purified sample having a 100 times lower concentration in ²⁹Si (500 ppm) than natural abundance. As a result, the dominant hyperfine interactions in the ESEEM signal are very low (on the order of 100 Hz) and have to be detected at low magnetic fields (around 0.1 mT). These results bring quantum-limited EPR spectroscopy one step closer to real-world applications.

2 ESEEM spectroscopy: theory

2.1 Phenomenology

We start by briefly discussing the ESEEM phenomenon. Consider an ensemble of electron spins placed in a magnetic field B₀. The spin ensemble linewidth Γ is broadened by a variety of mechanisms: spatial inhomogeneity of the applied field B₀, local magnetic fields generated by magnetic impurities throughout the sample, and spatially inhomogeneous strain or electric fields. One prominent way to mitigate the effect of this inhomogeneous broadening is the spin-echo sequence (also called Hahn echo, or two-pulse echo). It consists of a π∕2 pulse at time t=0 and a π pulse after a delay τ (see Fig. 1a). This π pulse reverses the evolution of the phase of the precessing magnetic dipoles, which leads at a later time 2τ to their refocusing and the emission of a microwave pulse (the echo) of amplitude V_2p(τ).

Figure 1 Sequences used for two-pulse (a), three-pulse (b), and five-pulse (c) ESEEM measurements.

In general, V_2p(τ) decays monotonically; it can however also display oscillations. Such ESEEM was first observed by Mims et al. (1961) and Rowan et al. (1965) for Ce³⁺ ions in a CaWO₄ crystal and was interpreted as being caused by the dipolar interaction of the electronic spin of the Ce³⁺ ions with the ¹⁸³W nuclear spins of the crystal. The oscillation frequencies appearing in the ESEEM pattern are related to the nuclear-spin Larmor frequencies and to their coupling to the electron spin. As such, ESEEM measurements provide spectroscopic information on the nature of the nuclear-spin bath and its density, and ESEEM spectroscopy has become an essential tool in advanced EPR (Schweiger and Jeschke, 2001; Mims et al., 1990). ESEEM has also been observed for individual spins measured optically, in particular for individual nitrogen vacancy (NV) centers in diamond coupled to a bath of ¹³C nuclear spins (Childress et al., 2006). A more complete theory of ESEEM is presented in Mims (1972). Our goal here is to provide a simple picture of the physics involved, as well as to introduce useful formulas and notations.

2.2 Two-spin 1∕2 model

We follow the analysis in Schweiger and Jeschke (2001) of the model case depicted in Fig. 2a. An electron spin of $S=\mathrm{1}/\mathrm{2}$, with an isotropic g tensor, is coupled to a proximal nuclear spin of $I=\mathrm{1}/\mathrm{2}$. Both are subject to a magnetic field B₀ applied along z. The system Hamiltonian is

$$\begin{array}{}\text{(1)}& {\mathbf{H}}_{\mathrm{0}}={\mathbf{H}}_{\mathrm{e}}+{\mathbf{H}}_{\mathrm{n}}+{\mathbf{H}}_{\mathrm{hf}},\end{array}$$

where H_e=ω_SS_z (H_n=ω_II_z) is the Zeeman Hamiltonian of the electron (nuclear) spin with Larmor frequency ω_S (ω_I), and H_hf is the electron–nuclear hyperfine interaction, which includes their dipole–dipole coupling and may include a Fermi contact term as well. We assume that ω_S is much larger than the hyperfine interaction strength, in which case terms proportional to the S_x and S_y operators can be neglected. This secular approximation leads to a hyperfine Hamiltonian of the form ${\mathbf{H}}_{\mathrm{hf}}=A{\mathbf{S}}_z{\mathbf{I}}_z+B{\mathbf{S}}_z{\mathbf{I}}_x$, with the expressions for A and B depending on the details of the hyperfine interaction (Schweiger and Jeschke, 2001).

Overall, the system Hamiltonian is

$$\begin{array}{}\text{(2)}& {\mathbf{H}}_{\mathrm{0}}={\mathit{\omega }}_S{\mathbf{S}}_z+{\mathit{\omega }}_I{\mathbf{I}}_z+A{\mathbf{S}}_z{\mathbf{I}}_z+B{\mathbf{S}}_z{\mathbf{I}}_x.\end{array}$$

Figure 2 ESEEM model system for an electron spin of $S=\mathrm{1}/\mathrm{2}$ and nuclear spin of $I=\mathrm{1}/\mathrm{2}$ with ${\mathit{\omega }}_I,A,B>\mathrm{0}$. (a) Nuclear spin (purple) subject to external field B₀ and dipole field (blue) of a nearby electron spin (green) located at relative position r. (b) Energy diagram showing the electron transitions (green), the nuclear transitions (purple), and the (normally forbidden) electro-nuclear transitions (orange). The energy levels $|\mathrm{1}\rangle ,\mathrm{\dots },|\mathrm{4}\rangle$ are labeled according to the eigenstates of the Zeeman basis. (c) Quantization axes ω_↑ and ω_↓ due to mixing of the nuclear states, which results in inclination of the quantization axis from z by the angles η_↑ and η_↓, respectively. (d) EPR spectrum showing the electron transitions (green) and the electro-nuclear transitions (orange) as well as the relation of these electron spin resonance (ESR) transitions to the nuclear frequencies ω_↑ and ω_↓ (purple).

Because of the BS_zI_x term, the nuclear spin is subjected to an effective magnetic field whose direction (and magnitude) depend on the electron-spin state |↑_e〉 or |↓_e〉. Its eigenstates therefore depend on the electron-spin state, so that transitions become allowed between all the spin system energy levels |1〉−|4〉, leading to the ESEEM phenomenon. Relevant parameters are the electron-spin-state-dependent angles between the effective magnetic field seen by the nuclear spin and the quantization axis z:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& {\mathit{\eta }}_{\uparrow }=\arctan {\displaystyle \frac{B}{A+\mathrm{2}{\mathit{\omega }}_I}}\\ & {\mathit{\eta }}_{\downarrow }=\arctan {\displaystyle \frac{B}{A-\mathrm{2}{\mathit{\omega }}_I}},\end{array}\end{array}$$

and the electron-spin-dependent nuclear-spin frequencies:

$$\begin{array}{ll}{\displaystyle }{\mathit{\omega }}_{\uparrow }& {\displaystyle }=\left({\mathit{\omega }}_I+{\displaystyle \frac{A}{\mathrm{2}}}\right)\cos {\mathit{\eta }}_{\uparrow }-{\displaystyle \frac{B}{\mathrm{2}}}\sin {\mathit{\eta }}_{\uparrow }\\ {\displaystyle }{\mathit{\omega }}_{\downarrow }& {\displaystyle }=\left({\mathit{\omega }}_I-{\displaystyle \frac{A}{\mathrm{2}}}\right)\cos {\mathit{\eta }}_{\downarrow }-{\displaystyle \frac{B}{\mathrm{2}}}\sin {\mathit{\eta }}_{\downarrow }.\end{array}$$

When η_↑, η_↓ are close to equal, only the nuclear-spin-preserving transitions are allowed; this occurs either when B=0 (due to a specific orientation of the dipolar field or a purely isotropic hyperfine coupling), or when B≠0 but ω_I≫A (very weak-coupling limit) or ω_I≪A (very strong coupling limit). On the contrary, when the direction of the effective magnetic field seen by the nuclear spin is electron-spin dependent, all transitions become allowed. This occurs when B≠0 and ${\mathit{\omega }}_I\simeq \pm A/\mathrm{2}$.

2.3 Multi-pulse ESEEM

Because of the level structure shown in Fig. 2, and assuming for simplicity microwave pulses so short that their bandwidth is much larger than ${\mathit{\omega }}_{\uparrow ,\downarrow }$, microwave pulses at the electron-spin frequency ω_S excite the allowed transitions $|\mathrm{1}\rangle \leftrightarrow |\mathrm{3}\rangle$ and $|\mathrm{2}\rangle \leftrightarrow |\mathrm{4}\rangle$ but also the normally forbidden $|\mathrm{1}\rangle \leftrightarrow |\mathrm{4}\rangle$ and $|\mathrm{2}\rangle \leftrightarrow |\mathrm{3}\rangle$, leading to coherence transfer between the levels and to beatings. Note that for simplicity we assume that the microwave pulses are ideal and so short that their bandwidth is much larger than ω₁₂ and ω₃₄.

It is then possible to compute analytically the effect of a two-pulse echo sequence consisting of an instantaneous ideal π∕2 pulse and an instantaneous ideal π pulse (see Fig. 1), disregarding any decoherence. The resulting echo amplitude (Schweiger and Jeschke, 2001) is given by

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{V}_{\mathrm{2}\mathrm{p}}\left(\mathit{\tau }\right)& =\mathrm{1}-{\displaystyle \frac{k}{\mathrm{4}}}[\mathrm{2}-\mathrm{2}\cos \left({\mathit{\omega }}_{\uparrow }\mathit{\tau }\right)-\mathrm{2}\cos \left({\mathit{\omega }}_{\downarrow }\mathit{\tau }\right)\\ & +\cos \left(\left({\mathit{\omega }}_{\uparrow }-{\mathit{\omega }}_{\downarrow }\right)\mathit{\tau }\right)+\cos \left(\left({\mathit{\omega }}_{\uparrow }+{\mathit{\omega }}_{\downarrow }\right)\mathit{\tau }\right)],\end{array}\end{array}$$

with

$$\begin{array}{}\text{(5)}& k={\left[{\displaystyle \frac{B{\mathit{\omega }}_I}{{\mathit{\omega }}_{\uparrow }{\mathit{\omega }}_{\downarrow }}}\right]}^{\mathrm{2}}.\end{array}$$

The spin-echo amplitude is modulated by a function whose frequency spectrum and amplitude contain information about the nuclear-spin Larmor frequency ω_I as well as its hyperfine coupling (A,B) to the electron spin. The modulation contrast $\mathrm{0}\le k\le \mathrm{1}$ is maximal when transitions $|\mathrm{1}\rangle -|\mathrm{4}\rangle$ and $|\mathrm{2}\rangle -|\mathrm{3}\rangle$ are maximally allowed, corresponding to ${\mathit{\omega }}_I\simeq A/\mathrm{2}$.

The above results are exact, as long as the secular approximation is valid and the pulses are ideal. In the weak-coupling limit A, B≪ω_I, ${\mathit{\omega }}_{\uparrow }\simeq {\mathit{\omega }}_{\downarrow }\simeq {\mathit{\omega }}_I$ so that $V_{\mathrm{2}\mathrm{p}}\left(\mathit{\tau }\right)=\mathrm{1}-\frac{k}{\mathrm{4}}[\mathrm{3}-\mathrm{4}\cos ({\mathit{\omega }}_I\mathit{\tau })+\cos (\mathrm{2}{\mathit{\omega }}_I\mathit{\tau }\left)\right]$, with $k=(B/{\mathit{\omega }}_I{)}^{\mathrm{2}}\ll \mathrm{1}$. In this limit, the echo modulation spectrum directly yields the nuclear-spin Larmor frequency and also contains components at twice this frequency. Note however that, in practice, the π pulse bandwidth is always finite, because of the resonator bandwidth or limited pulse power; this sets a limit to the range of detectable modulation frequencies.

The electron spin is often coupled to N nuclear spins, with N>1. Since all nuclear-spin subspaces can be diagonalized separately, the total ESEEM modulation is simply given by the product of each nuclear-spin modulation V_2p, l(τ), l being the nuclear-spin index. Taking also into account that the electron spin is also subject to decoherence processes, modeled, for instance, by an exponential decay with time constant T₂, the echo envelope is

$$\begin{array}{}\text{(6)}& V_{\mathrm{2}\mathrm{p}}^{\prime }\left(\mathit{\tau }\right)=\exp (-\mathrm{2}\mathit{\tau }/T_{\mathrm{2}})\prod _{l=\mathrm{1}}^N{V}_{\mathrm{2}\mathrm{p},l}\left(\mathit{\tau }\right).\end{array}$$

The modulation pattern $V_{\mathrm{2}\mathrm{p}}^{\prime }\left(\mathit{\tau }\right)$ yields quantitative information about the nature and coupling of the nuclear spins surrounding the electron spin whose echo is measured and is therefore a useful tool in EPR spectroscopy. When the environmental nuclei have a certain probability p to be of a given isotope with a nuclear spin of $I=\mathrm{1}/\mathrm{2}$, and a probability 1−p to be of an isotope with I=0, the above formulas are straightforwardly modified (Rowan et al., 1965) by writing

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{V}_{\mathrm{2}\mathrm{p},l}\left(\mathit{\tau }\right)& =\mathrm{1}-{\displaystyle \frac{p{k}_l}{\mathrm{4}}}[\mathrm{2}-\mathrm{2}\cos \left({\mathit{\omega }}_{\uparrow ,l}\mathit{\tau }\right)-\mathrm{2}\cos \left({\mathit{\omega }}_{\downarrow ,l}\mathit{\tau }\right)\\ & +\cos \left(\left({\mathit{\omega }}_{\uparrow ,l}-{\mathit{\omega }}_{\downarrow ,l}\right)\mathit{\tau }\right)+\cos \left(\left({\mathit{\omega }}_{\uparrow ,l}+{\mathit{\omega }}_{\downarrow ,l}\right)\mathit{\tau }\right)].\end{array}\end{array}$$

The echo signal $V_{\mathrm{2}\mathrm{p}}^{\prime }\left(\mathit{\tau }\right)$ is the sum of terms that have the general form $p^L{\prod }_{l=\mathrm{1}}^{l=L}{k}_l\cos \left({\mathit{\omega }}_{\mathit{\mu },l}\mathit{\tau }\right)$, where l runs over a subset of L nuclei and $\mathit{\mu }=\uparrow ,\downarrow$. If p≪1, this expression is well approximated by keeping only the L=1 terms, which then yields

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{V}_{\mathrm{2}\mathrm{p}}\left(\mathit{\tau }\right)& \simeq \mathrm{1}-\sum _{l=\mathrm{1}}^{l=N}{\displaystyle \frac{p{k}_l}{\mathrm{4}}}[\mathrm{2}-\mathrm{2}\cos \left({\mathit{\omega }}_{\uparrow ,l}\mathit{\tau }\right)-\mathrm{2}\cos \left({\mathit{\omega }}_{\downarrow ,l}\mathit{\tau }\right)\\ & +\cos \left(\left({\mathit{\omega }}_{\uparrow ,l}-{\mathit{\omega }}_{\downarrow ,l}\right)\mathit{\tau }\right)+\cos \left(\left({\mathit{\omega }}_{\uparrow ,l}+{\mathit{\omega }}_{\downarrow ,l}\right)\mathit{\tau }\right)].\end{array}\end{array}$$

One limitation of the previous pulse sequence is that the modulation envelope can only be measured up to a time of order T₂ due to electron-spin decoherence, which may be too short for appreciable spectral resolution. This limitation can be overcome by the three-pulse echo sequence shown in Fig. 1b. It consists of a π∕2 pulse applied at t=0 followed, after a time τ chosen such that τ<T₂, by a second π∕2 pulse. After a variable delay T, a third π∕2 pulse is applied, leading to the emission of a stimulated echo at time $t=T+\mathrm{2}\mathit{\tau }$. The interest of this sequence is that the first pair of π∕2 pulses generates nuclear-spin coherence that can survive up to the nuclear-spin coherence time T_2, n which is in general much longer than T₂ (and close to the electron energy spin relaxation time T₁). An analytical formula can be derived for the three-pulse echo amplitude in the ideal pulse approximation (Schweiger and Jeschke, 2001):

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{V}_{\mathrm{3}\mathrm{p}}\left(T\right)& =\exp (-T/T_{\mathrm{2},\mathrm{n}})\exp (-\mathrm{2}\mathit{\tau }/T_{\mathrm{2}})\\ & \mathit{\{}\mathrm{1}-{\displaystyle \frac{k}{\mathrm{4}}}[\left[\mathrm{1}-\cos {\mathit{\omega }}_{\downarrow }\mathit{\tau }\right]\left[\mathrm{1}-\cos {\mathit{\omega }}_{\uparrow }(T+\mathit{\tau })\right]\\ & +\left[\mathrm{1}-\cos {\mathit{\omega }}_{\uparrow }\mathit{\tau }\right]\left[\mathrm{1}-\cos {\mathit{\omega }}_{\downarrow }(T+\mathit{\tau })\right]\left]\mathit{\right\}}.\end{array}\end{array}$$

Contrary to two-pulse ESEEM, three-pulse echo modulation as a function of T only contains the ${\mathit{\omega }}_{\downarrow },{\mathit{\omega }}_{\uparrow }$ frequency components, and not their sum or difference; that is, in the weak-coupling limit A, B≪ω_I, only the nuclear-spin Larmor frequency ω_I appears in the spectrum. Another difference is that the modulation pattern and amplitude depend on τ; in particular, its amplitude is zero whenever ${\mathit{\omega }}_{\downarrow ,\uparrow }\mathit{\tau }=\mathrm{2}\mathit{\pi }n$ with n integer (“blind spots”).

For weakly coupled nuclei, the modulation amplitude of three-pulse ESEEM can be enhanced by up to 1 order of magnitude by using a more complex pulse sequence known as five-pulse ESEEM (Schweiger and Jeschke, 2001; Kasumaj and Stoll, 2008) and shown in Fig. 1. The analytical formula for the five-pulse echo amplitude V_5p is given in the Supplement.

Equation (6), with proper modification to take into account contributions of different pathways, can be applied to the three- and five-pulse ESEEM to treat coupling to multiple nuclear spins. The details are shown in Sect. S3.3 of the Supplement.

2.4 Fictitious spin model

The electronic spins that we consider in this work involve an unpaired electron with spin of $S_{\mathrm{0}}=\mathrm{1}/\mathrm{2}$ either located around or trapped by an ionic defect, which itself can possess a non-zero nuclear spin of I₀. These two spins of the defect are strongly coupled and therefore form a multi-level system, which can nevertheless be mapped to an effective, fictitious, spin-1∕2 model as explained below (Schweiger and Jeschke, 2001), to which the model of Sect. 2.3 can be applied.

The system spin Hamiltonian writes

$$\begin{array}{}\text{(10)}& {\mathbf{H}}_{\mathrm{ion}}={\mathit{\beta }}_{\mathrm{e}}{\mathit{B}}_{\mathbf{0}}\cdot {\mathbf{g}}_{\mathrm{e}}\cdot {\mathbf{S}}_{\mathrm{0}}+{\mathbf{S}}_{\mathrm{0}}\cdot {\mathbf{A}}_{\mathrm{0}}\cdot {\mathbf{I}}_{\mathrm{0}}.\end{array}$$

Here, β_e is the electron Bohr magneton, g_e is the (possibly anisotropic) gyromagnetic tensor, and A₀ the hyperfine tensor. The nuclear Zeeman interaction of the defect system, being small compared to the hyperfine interaction in the range of magnetic fields explored here, is neglected from the Hamiltonian.

This multi-level electron-spin system is coupled to other nuclear spins in the lattice, giving rise to ESEEM. Consider a nuclear spin at a lattice site j, defined by its location r_j with respect to the electron spin. The nuclear Zeeman Hamiltonian is ${\mathbf{H}}_j={\mathit{\omega }}_I{\mathbf{I}}_{j,z}$, with ω_I=g_nβ_nB₀, g_n being the nuclear g factor and β_n the nuclear magneton. Its hyperfine coupling to the electron-spin system is described by the Hamiltonian:

$$\begin{array}{}\text{(11)}& {\mathbf{H}}_{j,\mathrm{hf}}={\mathbf{S}}_{\mathrm{0}}\cdot {\mathbf{A}}_j\cdot {\mathbf{I}}_j,\end{array}$$

with

$$\begin{array}{}\text{(12)}& {\mathbf{A}}_j={\mathbf{A}}_{j,\mathrm{cf}}+{\mathbf{A}}_{j,\mathrm{dd}}.\end{array}$$

This hyperfine tensor consists of a Fermi contact term, ${\mathbf{A}}_{j,\mathrm{cf}}=\frac{\mathrm{2}}{\mathrm{3}}{\mathit{\mu }}_{\mathrm{0}}{\mathit{\beta }}_{\mathrm{e}}{g}_{\mathrm{n}}{\mathit{\beta }}_{\mathrm{n}}{\mathbf{g}}_{\mathrm{e}}{\left|\mathit{\psi }\left({\mathit{r}}_j\right)\right|}^{\mathrm{2}}$, and a dipole–dipole term, ${\mathbf{A}}_{j,\mathrm{dd}}=\frac{\mathrm{3}{\mathit{\mu }}_{\mathrm{0}}}{\mathrm{4}\mathit{\pi }{\left|{\mathit{r}}_j\right|}^{\mathrm{5}}}{\mathit{\beta }}_{\mathrm{e}}{\mathit{\beta }}_{\mathrm{n}}{g}_{\mathrm{n}}\left[{\mathit{r}}_j^{\mathrm{2}}{\mathbf{g}}_{\mathrm{e}}-\mathrm{3}\left({\mathbf{g}}_{\mathrm{e}}\cdot {\mathit{r}}_j\right){\mathit{r}}_j\right]$, ψ(r_j), being the electron wave function at the nuclear-spin location.

The Hamiltonian H_ion (Eq. 10) can be diagonalized, yielding 4I₀+2 energy levels. It is in general possible to isolate two levels |α〉 and |β〉 that are coupled by an electron spin resonance (ESR)-allowed transition and are resonant or quasi-resonant with the microwave cavity, with a transition frequency ω_S. If these two levels are sufficiently separated in energy from other levels of H_ion, they define a fictitious $S=\mathrm{1}/\mathrm{2}$ system. Writing the total Hamiltonian ${\mathbf{H}}_{\mathrm{ion}}+{\mathbf{H}}_j+{\mathbf{H}}_{j,\mathrm{hf}}$ restricted to this two-dimensional subspace yields

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\mathbf{H}}_{\mathrm{0}}& ={\mathit{\omega }}_S{\mathbf{S}}_z+\left({\mathit{\omega }}_I+{\displaystyle \frac{m_S^{\mathit{\alpha }}+m_S^{\mathit{\beta }}}{\mathrm{2}}}{\mathbf{A}}_{j,zz}\right){\mathbf{I}}_{j,z}\\ & +{\displaystyle \frac{m_S^{\mathit{\alpha }}+m_S^{\mathit{\beta }}}{\mathrm{2}}}{\mathbf{A}}_{j,zx}{\mathbf{I}}_{j,x}\\ & +\left(m_S^{\mathit{\alpha }}-m_S{}^{\mathit{\beta }}\right)\left({\mathbf{A}}_{j,zz}{\mathbf{S}}_z{\mathbf{I}}_{j,z}+{\mathbf{A}}_{j,zx}{\mathbf{S}}_z{\mathbf{I}}_{j,x}\right),\end{array}\end{array}$$

where $m_S^{\mathit{\alpha },\mathit{\beta }}=\langle \mathit{\alpha },\mathit{\beta }\left|{\mathbf{S}}_{\mathrm{0},z}\right|\mathit{\alpha },\mathit{\beta }\rangle$.

Equation (13) maps the more complex system to the simple model of Sect. 2.2. Compared to Eq. (2), two differences appear. First, the hyperfine interaction parameters A, B are rescaled by the effective longitudinal magnetization difference $(m_S^{\mathit{\alpha }}-m_S^{\mathit{\beta }})$ which depends on the two levels considered. Second, when the average longitudinal magnetization of the two levels $(m_S^{\mathit{\alpha }}+m_S^{\mathit{\beta }})$ is non-zero, the nuclear spin sees an extra Zeeman contribution which may be tilted with respect to the z axis. Once taken into account, these corrections, the analysis, and formulas of Sect. 2.3 remain valid.

3 Spin systems

3.1 Erbium-doped CaWO₄

The first system investigated consists of erbium Er³⁺ ions doped into a CaWO₄ matrix, substituting Ca²⁺. The crystal has a tetragonal body-centered structure (see Fig. 3) with lattice constants $a=b=\mathrm{0.524}$ nm and c=1.137 nm. Rare-earth ions with an odd number of electrons such as Er³⁺ have a ground state consisting of two levels that are degenerate in zero magnetic field and separated from other levels by an energy scale equivalent to several tens of Kelvin due to the crystalline electric field and the spin–orbit interaction. This pair of electronic levels is known as a Kramers doublet and forms an effective $S_{\mathrm{0}}=\mathrm{1}/\mathrm{2}$ electron-spin system, with a spin Hamiltonian H_Er (Abragam and Bleaney, 2012) whose form is given by Eq. (10).

Due to the S4 site symmetry in which rare-earth ions are found in CaWO₄, the g tensor is diagonal in the crystallographic frame with $g_{xx}=g_{yy}=\mathrm{8.38}$ and g_zz=1.247 (Antipin et al., 1968) (x, y, z corresponding to a, b, c). Of all erbium atoms, 77 % are from an isotope that has nuclear spin of I₀=0 and therefore no contribution from the hyperfine term in Eq. (10). Their energy levels are shown in Fig. 3 for B₀ applied in the (a,b) plane.

The remaining 23 % are from the ¹⁶⁷Er isotope with $I_{\mathrm{0}}=\mathrm{7}/\mathrm{2}$. Its hyperfine coupling tensor to the Er³⁺ electron spin is diagonal, with coefficients $A_{xx}=A_{yy}=\mathrm{873}$ MHz and A_zz=130 MHz. The 16 eigenfrequencies of the ¹⁶⁷Er spin Hamiltonian are also shown in Fig. 3, again for B₀ applied in the (a,b) plane. In the high-magnetic-field limit of $B_{\mathrm{0}}\gg A_{\mathrm{Er}}/\left(g_{\mathrm{Er}}{\mathit{\beta }}_{\mathrm{e}}\right)$, the eigenstates are simply described by $|\pm ,m_I\rangle$, ± describing the electron-spin quantum number $m_S=\pm \mathrm{1}/\mathrm{2}$ and m_I the nuclear-spin quantum number. For B₀<100 mT, as is the case in the measurements described here, this limit is only approximate, but we will nevertheless use the high-field state vectors as labels for the lower-field eigenstates. The strongest EPR-allowed transitions are the m_I-preserving transitions. In the following, we will apply the fictitious spin model with $|\mathit{\alpha },\mathit{\beta }\rangle$ = $|\pm ,m_I\rangle$.

The CaWO₄ matrix also contains nuclear spins. Indeed, the ¹⁸³W isotope has a spin of $I=\mathrm{1}/\mathrm{2}$ with nuclear g factor g_n=0.235 (corresponding to a gyromagnetic ratio of 1.8 MHz T⁻¹) and is present in a p=0.13 abundance, whereas the other tungsten isotopes are nuclear-spin free. The interaction of the ¹⁸³W atoms with the erbium ions gives rise to the ESEEM studied below. Because the 4f electron wave function is mainly located on the Er³⁺ ion, the contact hyperfine with the nuclear spins of the lattice is expected to be negligibly small. We therefore model the hyperfine interaction with ¹⁸³W by the dipole–dipole term in Eq. (12).

Figure 3 Structure and energy diagram of erbium ions in CaWO₄. (a) Crystal structure with oxygen atoms hidden for clarity. Erbium atoms are in substitution of the calcium. The crystal has a rotational symmetry around the c axis. A fraction p=0.13 of the W atoms are of the ¹⁸³W isotope, with a nuclear spin of 1∕2. (b) Energy level diagram of the I=0 erbium isotopes (dashed black line) and of the ¹⁶⁷Er isotope (solid black lines) with $I=\mathrm{7}/\mathrm{2}$ for B₀ applied perpendicular to the c axis. Vertical red lines indicate the value of B₀ for which an allowed EPR transition becomes resonant with the 4.372 GHz frequency of our detection resonator (see main text, Sect. 4). Four red arrows indicate the values of B₀ at which ESEEM data were measured.

3.2 Bismuth donors in silicon

The other system considered is the bismuth donor in silicon. Bismuth, as an element of the fifth column, substitutes in the silicon lattice by making four covalent bonds with neighboring atoms, leaving one unpaired electron that can be weakly trapped by the hydrogenic potential generated by the Bi⁺ ion, whose spin gives rise to the resonance signal (see Fig. 4a). The donor wave function ψ(r) has a complex structure that extends over ≈1.5 nm in the silicon lattice (Kohn and Luttinger, 1955; Feher, 1959) (see the Supplement). As for Er:CaWO₄, the donor spin Hamiltonian H_Bi is given by Eq. (10). However, in this case, the g tensor g_e1 is isotropic with g_e=2, and the hyperfine tensor A_Bi1 with the nuclear spin of $I_{\mathrm{0}}=\mathrm{9}/\mathrm{2}$ of the bismuth atom is also isotropic, with $A_{\mathrm{Bi}}/\mathrm{2}\mathit{\pi }=\mathrm{1.4754}$ GHz.

The eigenstates of H_Bi have simple properties because of its isotropic character. Denoting m_S (m_I) the eigenvalue of S_z,0 (I_z,0), we note that $m=m_I+m_S$ is a good quantum number since H_Bi commutes with ${\mathbf{S}}_{z,\mathrm{0}}+{\mathbf{I}}_{z,\mathrm{0}}$ (Mohammady et al., 2010), z being the direction of B₀. States with equal m are hybridized by H_Bi. States $|m=\mathrm{5}\rangle$ and $|m=-\mathrm{5}\rangle$, corresponding to $|m_S=+\mathrm{1}/\mathrm{2},m_I=\mathrm{9}/\mathrm{2}\rangle$ and $|m_S=-\mathrm{1}/\mathrm{2},m_I=-\mathrm{9}/\mathrm{2}\rangle$, are non-degenerate and are thus also eigenstates of H_Bi. States with $\left|m\right|\le \mathrm{4}$ belong to nine two-dimensional subspaces spanned by $|m_S=+\mathrm{1}/\mathrm{2},m_I=m-\mathrm{1}/\mathrm{2}\rangle ,|m_S=-\mathrm{1}/\mathrm{2},m_I=m+\mathrm{1}/\mathrm{2}\rangle$ within which the two eigenstates of H_Bi are given by $|\pm ,m\rangle =a_m^{\pm }|\pm \frac{\mathrm{1}}{\mathrm{2}},m\mp \frac{\mathrm{1}}{\mathrm{2}}\rangle +b_m^{\pm }|\mp \frac{\mathrm{1}}{\mathrm{2}},m\pm \frac{\mathrm{1}}{\mathrm{2}}\rangle$, with values of $a_m^{\pm },b_m^{\pm }$ that can be determined analytically (Mohammady et al., 2010).

Contrary to the erbium case, the measurements of bismuth donor spins are performed in the low-field limit |g_eβ_eB₀|≪|A_Bi|, in which the eigenstates are fully hybridized. In this limit, a useful approximate expression for the eigenenergy of level $|\pm ,m\rangle$ is

$$\begin{array}{}\text{(14)}& E_m^{\pm }\approx -{\displaystyle \frac{A_{\mathrm{Bi}}}{\mathrm{2}}}\pm {\displaystyle \frac{\mathrm{5}{A}_{\mathrm{Bi}}}{\mathrm{2}}}\pm {\displaystyle \frac{m{g}_{\mathrm{e}}{\mathit{\beta }}_{\mathrm{e}}\left|{\mathit{B}}_{\mathrm{0}}\right|}{\mathrm{10}}}.\end{array}$$

The magnetic-field dependence of the $|\pm ,m\rangle$ energy levels is shown in Fig. 4b for B₀<1 mT. Note in particular that the separation between neighboring hyperfine levels is given by $E_m^{\pm }-E_{m-\mathrm{1}}^{\pm }\approx \pm \frac{g_{\mathrm{e}}{\mathit{\beta }}_{\mathrm{e}}{B}_{\mathrm{0}}}{\mathrm{10}}=\pm \mathrm{2}\mathit{\pi }\times \mathrm{2.8}{B}_{\mathrm{0}}$ GHz.

Because of the hybridization, all transitions that satisfy $\left|\mathrm{\Delta }m\right|=\mathrm{1}$ are to some extent EPR allowed at low field; i.e., they have a non-zero matrix element of operator S_0,x. In this work, we particularly focus on the 18 $\left|\mathrm{\Delta }m\right|=\mathrm{1}$ transitions that are in the ≃7 GHz frequency range at low magnetic fields $|+,m\rangle \leftrightarrow |-,m-\mathrm{1}\rangle$ and $|-,m\rangle \leftrightarrow |+,m-\mathrm{1}\rangle$, as shown in Fig. 4c. The $|-,m\rangle \leftrightarrow |+,m+\mathrm{1}\rangle$ and $|-,m+\mathrm{1}\rangle \leftrightarrow |+,m-\mathrm{1}\rangle$ transitions are degenerate in frequency for $-\mathrm{4}\le m<\mathrm{4}$ as seen from Eq. (14), which results in only 10 different transition frequencies (see Figs. 4b–c and 8a).

The most abundant isotope of silicon is ²⁸Si, which is nuclear-spin free. The lattice also contains a small percentage p of ²⁹Si atoms that have a nuclear spin of $I=\mathrm{1}/\mathrm{2}$ and give rise to the ESEEM. The g factor of ²⁹Si is $g_{\mathrm{n}}=-\mathrm{1.11}$, yielding a gyromagnetic ratio of 8.46 MHz T⁻¹.

The donor–²⁹Si hyperfine interaction is given by Eq. (12). Due to the spatial extent of the electron wave function, the Fermi contact term is not negligible and needs to be taken into account together with the dipole–dipole coupling (Hale and Mieher, 1969); more details can be found in the Supplement.

The restriction of the total system Hamiltonian to each of the 18 ESR-allowed transitions of the bismuth donor manifold can be mapped onto the fictitious spin-1∕2 model of Sect. 2.4. Note however that the hyperfine term $\left|A_j\right|$ can take values up to ∼1 MHz for proximal nuclear spins, which is comparable to or larger than the frequency difference between hyperfine states of the bismuth donor manifold at low field as explained above. The validity of the fictitious spin-1∕2 model in this context will be discussed in Sect. 5.

Figure 4 Structure and energy diagram of bismuth donors in silicon. (a) Silicon crystal structure, showing a substitutional bismuth atom coupled to nearby ²⁹Si nuclear spins. The donor electron is trapped around the Bi⁺ ion and its wave function covers many lattice sites. (b) Energy levels of the bismuth donor for B₀<1 mT. (c) Schematic representation of the allowed transitions (black and grey arrows) between the bismuth donor energy levels in the low field limit.

4 Experimental setup and samples

The EPR spectrometer has been described in detail in Bienfait et al. (2015) and Probst et al. (2017) and is shown schematically in Fig. 5a. It is built around a superconducting micro-resonator of frequency ω_r consisting of a planar interdigitated capacitor shunted by an inductor, directly patterned on the crystal. We detect the spins that are located in the immediate vicinity of the resonator inductance. Note that the microwave B₁ field generated by the inductance is spatially inhomogeneous. If the spin location is broadly distributed, this can make the application of control pulses with a well-defined Rabi angle problematic (Ranjan et al., 2020a). As explained below, the resonator is more strongly coupled to the measurement line than in Bienfait et al. (2015) to increase the measurement bandwidth as requested for ESEEM spectroscopy.

Figure 5 Experimental setup and samples. (a) Schematic of the low-temperature EPR spectrometer. The LC resonator is inductively coupled to electron spins, which are coupled to a nuclear-spin bath that causes the ESEEM. The spins are probed by sequences of microwave pulses at the resonator frequency ${\mathit{\omega }}_{\mathrm{r}}=\mathrm{1}/\sqrt{\text{LC}}$. Reflected pulses as well as the echo signal are routed to a parametric amplifier, further amplified at 4 K, and finally demodulated and digitized at room temperature. (b, c) Design of the LC resonator used for the detection of erbium ion spins, with a 725 µm long, 5 µm wide inductor. It is patterned out of a 100 nm thick niobium film deposited on top of a CaWO₄ substrate bulk doped with Er³⁺ ions. (d, e) Design of the LC resonator used for the detection of bismuth donor spins, with a 100 µm long, 0.5 µm wide inductor. It is patterned out of a 100 nm thick aluminum film deposited on top of a silicon substrate isotopically enriched in ²⁸Si, in which bismuth ions were implanted at a 50–100 nm depth.

The sample is mounted in a copper sample holder thermally anchored at the mixing chamber of a dilution refrigerator. A DC magnetic field B₀ is applied parallel to the sample surface and along the resonator inductance. The resonator is coupled capacitively to an antenna, which is itself connected to a microwave measurement setup in reflection. To minimize heat load, the coaxial cables between 4 K and 10 mK are in superconducting niobium–titanium (NbTi). To suppress thermal noise, the input line is heavily attenuated at low temperatures. Microwave pulses for driving the spins are sent to the resonator input, and their reflection or transmission, together with the echo signal emitted by the spins, is fed into a superconducting JPA, either of the flux-pumped type (Zhou et al., 2014) or of the Josephson traveling-wave parametric amplifier (JTWPA) type (Macklin et al., 2015). Further microwave amplification takes place at 4 K with a high-electron-mobility transistor (HEMT) from Low Noise Factory, and then at room temperature, before homodyne demodulation which yields the two signal quadratures [I(t),Q(t)]. The echo-containing quadrature signal is integrated to yield the echo amplitude A_e. Such a setup was shown to reach sensitivities of the order of 10²–10³ spin/$\sqrt{\mathrm{Hz}}$ (Bienfait et al., 2015; Eichler et al., 2017; Probst et al., 2017).

Because of the small resonator mode volume and high quality factor, little microwave power is needed to drive the spins. The exact amount depends on the resonator geometry, as conveniently expressed by the power-to-field conversion factor $\mathit{\alpha }=B_{\mathrm{1}}/\sqrt{P_{\mathrm{in}}}$. In the experiments reported here, the maximum microwave power used to drive the spins is on the order of 10 nW. At this power, the superconducting pre-amplifiers saturate; however, they recover rapidly enough (within a few microseconds) to amplify the much weaker subsequent spin echoes. Flux-pumped JPAs are moreover switched off during the control pulses by pulsing the pump tone, whereas the JTWPA was kept on all the time. All microwave powers reaching the 4 K HEMT are low enough that neither saturation nor damage are to be expected at this stage.

The erbium-doped sample (from Scientific Materials) was prepared by mixing erbium oxide with calcium and tungsten oxides before crystal growth, yielding a uniform Er concentration of 6×10¹⁷ cm⁻³ (50 ppm) throughout the sample. For resonator fabrication, the bulk crystal was cut and polished to a thin rectangular sample with dimensions $\mathrm{0.4},\mathrm{mm}\times \mathrm{3}\mathrm{mm}\times \mathrm{6},\mathrm{mm}$ parallel to $a\times b\times c$ axes. The resonator was patterned out of a 100 nm thick (sputtered) niobium (Nb) layer, using a design similar to that shown in Bienfait et al. (2015). More specifically, 15 interdigitated fingers on either side of a 720 µm×5  µm inductive wire form an LC resonator, corresponding to a detection volume of V_Er∼20 pL. In the absence of magnetic field, the resonance frequency is ${\mathit{\omega }}_{\mathrm{r}}/\mathrm{2}\mathit{\pi }=\mathrm{4.323}$ GHz. Its total quality factor of 8×10³ is set both by the internal losses, characterized by the energy loss rate ${\mathit{\kappa }}_{\mathrm{i}}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{5}}$ s⁻¹, and by its coupling to the measurement line ${\mathit{\kappa }}_{\mathrm{C}}=\mathrm{3}\times {\mathrm{10}}^{\mathrm{6}}$ s⁻¹. For this geometry, the power-to-field factor is α=1.7 T W${}^{-\mathrm{1}/\mathrm{2}}$.

The bismuth donors have been implanted at ≈100 nm depth with a peak concentration of 8×10¹⁶ cm⁻³ in a silicon sample. They lie in a 700 nm thick silicon epilayer enriched in the nuclear-spin-free ²⁸Si isotope (nominal concentration of 99.95 %), grown on top of a natural-abundance silicon sample. The resonator is patterned out of a 50 nm thick aluminum film. It has the same geometry as reported in Probst et al. (2017), with a 100 µm long, 500 nm wide inductor, and a detection volume of 0.2 pL. Its frequency ${\mathit{\omega }}_{\mathrm{r}}/\mathrm{2}\mathit{\pi }=\mathrm{7.370}$ GHz is only slightly below the zero-field splitting of unperturbed Bi:Si donors at $\mathrm{5}{A}_{\mathrm{Bi}}/\left(\mathrm{2}\mathit{\pi }\right)=\mathrm{7.37585}$ GHz (Wolfowicz et al., 2013). The resonator internal loss is given by ${\mathit{\kappa }}_{\mathrm{i}}=\mathrm{3}\times {\mathrm{10}}^{\mathrm{5}}$ s⁻¹. The coupling to the measurement line can be tuned at will by modifying the length of a microwave antenna that capacitively couples the measurement waveguide to the on-chip resonator via the copper sample holder (Bienfait et al., 2015; Probst et al., 2017). For the experiments reported below, we used two settings: one for which the resonator was overcoupled (κ_C1=10⁷ s⁻¹), corresponding to a loaded quality factor of $Q_{\mathrm{1}}=\mathrm{4}\times {\mathrm{10}}^{\mathrm{3}}$, and one for which the coupling was closer to critical (κ_C2=10⁶ s⁻¹), corresponding to a loaded quality factor of $Q_{\mathrm{2}}=\mathrm{3.4}\times {\mathrm{10}}^{\mathrm{4}}$. In the low-Q case, square microwave pulses were used, of duration ≃100 ns similar to the cavity field damping time. In the high-Q case, shaped pulses were used (Probst et al., 2019) so that the intra-cavity field was a square pulse of 1 µs without any ringing. In some experiments, we additionally used a train of π pulses (Carr–Purcell–Meiboom–Gill sequence; CPMG), which generated extra echoes for significant gain in signal-to-noise ratio. More details on the pulse sequences used, the phase cycling scheme, and the repetition time, will be given in the following sections, together with experimental results. For this geometry, the power-to-field factor is α=9 T W${}^{-\mathrm{1}/\mathrm{2}}$ for the low-Q case, and α=21 T W${}^{-\mathrm{1}/\mathrm{2}}$ for the high-Q case.

5 Results

5.1 Erbium-doped CaWO₄

5.1.1 Spectroscopy

Figure 6 shows a spectrum comprising a series of microwave transmission measurements recorded on a vector network analyzer, measured at 100 mK, as a function of the magnetic field B₀ applied along the b crystal axis (Probst et al., 2020). Note that compared to Fig. 5a, the resonator is coupled to the measurement line in a hanger geometry (Day et al., 2003), so that its resonance appears as a dip in the amplitude transmission coefficient $\left|S_{\mathrm{21}}\right|$ (see Fig. 6). The nine red lines indicate the values of B₀ at which the calculated Er³⁺ ion transitions are equal to ω_r (see Fig. 3b). Avoided level crossings are observed, which indicate a strong coupling of the resonator to the erbium transitions. Several additional anti-crossings and discontinuities are visible above 40 mT. These are attributed to ytterbium impurities ¹⁷¹Yb and ¹⁷³Yb and magnetic flux vortices penetrating the resonator.

Noticeable in the spectrum at 37 mT is the large anti-crossing attributed to the highly concentrated I=0 erbium isotopes. Here, the high-cooperativity regime (C>30) is reached between the electronic spins and the resonator (Kubo et al., 2010; Probst et al., 2013). Typical linewidth of $\mathrm{\Gamma }/\mathrm{2}\mathit{\pi }\sim \mathrm{20}$ MHz is observed. The coupling strength is also observed to be different for the eight ¹⁶⁷Er transitions, which are labeled according to their corresponding nuclear-spin projections m_I. This is explained by the partial polarization of the ground-state hyperfine levels of ¹⁶⁷Er³⁺ at millikelvin temperatures (see Fig. 3b).

Figure 6 Spectroscopy of ${\mathrm{Er}}^{\mathrm{3}+}:{\mathrm{CaWO}}_{\mathrm{4}}$. Transmission coefficient $\left|S_{\mathrm{21}}\right|\left(\mathit{\omega }\right)$ at 100 mK as a function of the magnetic field B₀ applied along to the a crystalline axis, around 4.323 GHz. Vertical red lines indicate the expected erbium transitions either for the I=0 isotopes (dashed) or the $I=\mathrm{7}/\mathrm{2}$ isotope (solid). Red arrows indicate the field at which the ESEEM data are measured.

5.1.2 Two-pulse ESEEM

Four values of B₀ were selected for investigating ESEEM, indicated by the arrows in Fig. 6; the first, second, and fourth corresponding to electronic-spin transitions of ¹⁶⁷Er, and the third one to the I=0 isotopes. The two-pulse echo sequence of Fig. 1a was implemented with square pulses of 1 µs duration applied at the resonator input, with double amplitude for the second pulse. Note that due to the B₁ spatial inhomogeneity combined with the homogeneous spin distribution throughout the crystal, the spread of Rabi frequency is too large to observe a well-defined nutation signal. The Rabi angle is therefore not well defined, and the echo is the average of different rotation angles.

The control pulses driving the spins are filtered by the resonator bandwidth $\mathit{\kappa }/\mathrm{2}\mathit{\pi }\simeq \mathrm{600}$ kHz, corresponding to a field decay time $\mathrm{2}{\mathit{\kappa }}^{-\mathrm{1}}=\mathrm{3.3}$ µs. The repetition time between echo sequences was 1 s, close to the spin relaxation time T₁∼1–2 s measured by saturation recovery on the transitions studied. The echo signal was averaged 10 times with phase cycling of the π pulse to improve signal-to-noise ratio and to remove signal offsets.

Figure 7 shows the two-pulse echo integrated amplitude A_e as a function of τ for each of the four Er transitions investigated (Probst et al., 2020). A clear envelope modulation signal is observed, together with an overall damping. Here, we are interested only in the modulation pattern; a detailed study of the coherence time T₂ will be provided elsewhere. Qualitatively, we observe that the modulation frequency increases with B₀ and the modulation amplitude overall decreases with B₀, as expected from the discussion in Sect. 2. A Fourier transform of the I=0 data (see Fig. 7b) shows the ESEEM spectrum. Well-resolved peaks are observed in the 5–100 kHz range, distributed around the ¹⁸³W bare Larmor frequency ω_W.

A very rough estimate of the number of erbium ions contributing to the signal is [Er]V_Erκ∕Γ, which is 2.5×10⁸ for the I=0 data, and 10⁷ for each ¹⁶⁷Er transition.

Figure 7 Two-pulse ESEEM on Er:CaWO₄. (a) Integrated echo area A_e as a function of the inter-pulse delay τ for four values of B₀ corresponding to different transitions. Open circles are measurements, and solid lines are the results of the ESEEM calculations as explained in Sect. S5.1. (b) Measured (open red circles) and computed (solid line) fast Fourier transform of the I=0 data. The dashed blue line shows the Larmor frequency of ¹⁸³W nuclei in free space.

5.1.3 Comparison with the model

We compute the echo envelope $V_{\mathrm{2}\mathrm{p}}^{\prime }\left(\mathit{\tau }\right)$ described in Sect. 2.3, with the nearest 1000 coupled tungsten nuclei (N=1000) and a natural ¹⁸³W abundance of 14.4 % (p=0.144). The hyperfine interaction is taken to be purely dipolar, as already explained (Guillot-Noël et al., 2007; Car et al., 2018). The fitting proceeds by assigning an initial “guess” to six free parameters, then minimizing using the limited-memory Broyden–Fletcher–Goldfarb–Shanno algorithm for bound-constrained optimization (L-BFGS-B) (Byrd et al., 1995). Three of these parameters $\left(\left|B_{\mathrm{0}}\right|,\mathit{\varphi },\mathit{\theta }\right)$ describe the applied magnetic field:

$${\displaystyle }{\mathit{B}}_{\mathrm{0}}=\left|{\mathit{B}}_{\mathrm{0}}\right|\left[\sin \mathit{\theta }\cos \mathit{\varphi }\mathit{x}+\sin \mathit{\theta }\sin \mathit{\varphi }\mathit{y}+\cos \mathit{\theta }\mathit{z}\right].$$

Here, θ is the angle of the field relative to the crystal c axis (z) and ϕ is the angle relative to the a axis (x) in the a–b plane (x–y plane). The other three parameters (C, T₂, n) account for the echo envelope decay:

$${\displaystyle }\mathrm{1}{A}_{\mathrm{e}}\left(\mathit{\tau }\right)=V_{\mathrm{2}\mathrm{p}}\left(\mathit{\tau }\right)\cdot C\exp {\left(-{\displaystyle \frac{\mathrm{2}\mathit{\tau }}{T_{\mathrm{2}}}}\right)}^n,$$

where C represents the signal magnitude, T₂ represents the coherence time, and $n\in [\mathrm{1},\mathrm{2}]$ accounts for non-exponential decay. To determine the global minimum of the fit, the minimization is repeated 200 times with randomly seeded initial values for the six parameters, bounded within the known uncertainty of the applied magnetic field B₀, signal strength C, and coherence time T₂. This approach reveals single local minima for each fitted parameter within the bounded range, with the variance of the 200 outcomes determining the uncertainty for each parameter. In particular, it yields precise values for the angles $\mathit{\theta }=\mathrm{91.47}\pm \mathrm{0.01}$^∘ and $\mathit{\varphi }=\mathrm{90.50}\pm \mathrm{0.01}$^∘. The result of this fitting is presented in Fig. 7a, overlaid on the data for the I=0 transition at 37 mT. Only the decay parameters (C, T₂, n) and magnetic-field magnitude $\left|B_{\mathrm{0}}\right|$ are left free when fitting the other three transitions in Fig. 7a. This was done for consistency between datasets, and because the I=0 data yield the most accurate values for ϕ and θ due to the low decoherence rate. The fits yield coherence times T₂ varying between 40 and 400 µs, depending on the transition considered. Good agreement was also reached between the fitted and expected (pre-calibrated) field magnitudes.

Note that good fits to the data are also achieved by including only the nearest 100 tungsten nuclei, although noticeable deviations between the data and fit are observed with any amount less. The dimensionless “anisotropic hyperfine interaction parameter” ρ described in the seminal publication on ESEEM (Rowan et al., 1965) is not required here. This parameter was introduced with the earliest attempts of ESEEM fitting, likely to compensate for the low number of simulated nuclear spins (typically 10 nearest nuclei or less) and was interpreted as an account for a potential distortion of the local environment caused by dopant insertion. Finally, a consideration of the spectral components presented in Fig. 7b helps to more clearly identify the difference between the fit and the data. In particular, the high-frequency components of the fitted model are not present experimentally due to the filtering effect of the superconducting resonance (260 kHz half width at half maximum; HWHM). This high-Q resonator greatly reduces the bandwidth of the radio frequency (RF) field absorbed by the coupled Er−¹⁸³W system and further limits the bandwidth of the detected echo signal.

5.2 Bismuth donor sample

5.2.1 Spectroscopy

Given the resonator frequency ω_r, four bismuth donor resonances should be observed when varying B₀ between 0 and 1 mT, as seen in Fig. 8a. Figure 8b shows an echo-detected field sweep, measured at 12 mK: the integrated area A_e of echoes obtained with a sequence shown in Fig. 1a with τ=50 µs pulse separation is plotted as a function of B₀  (Probst et al., 2020). Instead of showing well-separated peaks as in the erbium case, echoes are observed for all fields below 1 mT, with a maximum close to 0.1 mT, and extends in particular down to B₀=0 mT . This is the sign that each of the expected peaks is broadened and overlaps with neighboring transitions. Close to zero field, the echo amplitude goes down by a factor of 2 on a scale of ∼0.1 mT, before showing a sharp increase at exactly zero field. These zero-field features are not currently understood, but they are reproducible as confirmed by the measurements at B₀<0, which are approximately symmetric to the B₀>0 data as they should be.

Line broadening was reported previously for bismuth donors in silicon in related experiments (Bienfait et al., 2015; Probst et al., 2017) and was attributed to the mechanical strain exerted by the aluminum resonator onto the silicon substrate due to differential thermal contractions between the metal and the substrate. At low strain, A_Bi depends linearly on the hydrostatic component of the strain tensor ${\mathit{ϵ}}_{\mathrm{hs}}=({\mathit{ϵ}}_{xx}+{\mathit{ϵ}}_{yy}+{\mathit{ϵ}}_{zz})/\mathrm{3}$ with a coefficient $\mathrm{d}{A}_{\mathrm{Bi}}/\mathrm{d}{\mathit{ϵ}}_{\mathrm{hs}}/\left(\mathrm{2}\mathit{\pi }\right)=\mathrm{28}$ GHz (Mansir et al., 2018). Quantitative understanding of the line shape was achieved in a given sample geometry based on this mechanism (Pla et al., 2018), using a finite-element modeling to estimate the strain profile induced upon sample cooldown. A similar modeling was performed for the Bi sample reported here (see Fig. 8d). Based on the typical strain distribution $\left|{\mathit{ϵ}}_{\mathrm{hyd}}\right|\sim \mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$ and on the hyperfine to strain coefficient $\mathrm{d}{A}_{\mathrm{Bi}}/\mathrm{d}{\mathit{ϵ}}_{\mathrm{hs}}/\left(\mathrm{2}\mathit{\pi }\right)=\mathrm{28}$ GHz, we expect the zero-field splitting 5A_Bi∕(2π) to have a spread of ∼50 MHz, which would indeed result in complete peak overlap in the B₀<1 mT region, as observed in Fig. 8b.

This broadening has two consequences worth highlighting. First, the bismuth donor echo signals can be measured down to B₀=0 mT, which otherwise is generally impossible in X-band spectroscopy. Here, this is enabled by the large hyperfine coupling of the Bi:Si donor, combined with strain-induced broadening. This makes it possible to detect ESEEM caused by very weakly coupled nuclear spins, which requires low magnetic fields as explained in Sect. 2. Second, at a given magnetic field, the spin-echo signal contains contributions from several overlapping EPR transitions. This last point is best understood from Fig. 8c, which shows how several classes of bismuth donors, each with different hyperfine coupling A_Bi, may have transitions resonant with ω_r. We will assume in the following that the inhomogeneous distribution of A_Bi is so broad that each of the 10 A_Bi values for which one bismuth donor transition is resonant with ω_r at fixed B₀ is equally probable, which is likely to be valid for B₀<1 mT.

Figure 8 (a) EPR-allowed transitions of a bismuth donor in silicon for $\mathrm{0}<B_{\mathrm{0}}<\mathrm{1}$ mT. The dashed red line denotes the resonator frequency ω_r. The spectrum is for an unstrained donor for which the frequency at B₀=0 is 5A_Bi∕(2π). (b) Echo-detected field sweep. The echo integral A_e is plotted versus B₀. (c) Frequency of all 18 bismuth donor transitions that may contribute to the echo signal at a given field (here, B₀=0.23 mT). This is made possible by the strain-induced spread in A_Bi between different donors. (d) Hydrostatic component of strain in silicon simulated using COMSOL.

5.2.2 Two-pulse ESEEM

Two-pulse echoes are measured with the pulse sequence shown in Fig. 1, which consists of a square π∕2_X pulse of duration 50 ns followed by a square π_Y pulse of duration 100 ns after a delay τ. Note that due to the donor spatial location in a shallow layer below the surface and to the strain shifting of their Larmor frequency (Pla et al., 2018), the Rabi frequency is more homogeneous than in the erbium-doped sample, and Rabi rotations with a well-defined angle can be applied (Pla et al., 2018; Probst et al., 2017). To increase the signal-to-noise ratio, a CPMG sequence of 198 π pulses separated by 10 µs are used following the echo sequence (Probst et al., 2017). The curves are repeated 20 times, with a delay of 2 s in between to enable spin relaxation of the donors. All the resulting echoes are then averaged. Phase cycling is performed by alternating sequences with opposite phases for the π∕2 pulses and subtracting the resulting echoes. The data are obtained in the low-Q configuration (see Sect. 4).

Figure 9 Two-pulse ESEEM of Bi:Si donors. (a) Echo integral A_e versus inter-pulse delay τ for a two-pulse echo sequence for varying magnetic field B₀. Dots indicate experimental data; lines indicate results of the model (see text), assuming a concentration in ²⁹Si of $p=\mathrm{4.4}\times {\mathrm{10}}^{-\mathrm{4}}$. The curves are vertically shifted for clarity. (b) Fit residue χ² for different ²⁹Si relative abundance p. The best fit is obtained for $p=\mathrm{4.4}\pm \mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$, in agreement with the specified value.

Figure 9 shows the integral of the averaged echoes A_e(τ) as a function of τ for various values of B₀ (Probst et al., 2020). At a non-zero field, A_e(τ) shows B₀-dependent oscillations on top of an exponential decay with time constant T₂=2.6 ms. Similar decay times were measured on the same chip with another resonator (Probst et al., 2017), and are attributed to a combination of donor–donor dipolar interactions and magnetic noise from defects at the sample surface.

Figure 10 Amplitude of the Fourier transform of the experimental (a) and theoretical (b) two-pulse Bi:Si donor ESEEM data.

In the subsequent discussion, we concentrate on the ESEEM pattern. To analyze the data, each curve was divided by a constant exponential decay with 2.6 ms time constant, mirrored at t=0, and Fourier transformed (see Fig. 10). Only two peaks are observed. Their frequencies vary linearly with B₀, and are found to be approximately 8 and 16 kHz mT⁻¹. This is in good agreement with the gyromagnetic ratio of ²⁹Si (8.46 kHz mT⁻¹); the presence of the second peak at twice this value is expected as explained in Sect. 2 for the two-pulse ESEEM in the weak-coupling limit. The oscillation amplitude goes down with B₀, again as expected from the model put forward in Sect. 2.

A rough estimate of the number of donors contributing to the measurements shown in Fig. 9 can be obtained by comparison with (Probst et al., 2017). Given the nearly identical resonator geometry, and assuming identical strain broadening in both samples, the ratio of the number of donors involved in both measurements is simply given by the ratio of resonator bandwidths. For the low-Q configuration, such as the two-pulse echo of Fig. 9, this corresponds to $\simeq \mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$ dopants; in the high-Q configuration (see the three- and five-pulse data in the next paragraph), this number is reduced to $\simeq \mathrm{5}\times {\mathrm{10}}^{\mathrm{2}}$ dopants.

5.2.3 Three- and five-pulse ESEEM

The spectral resolution provided by the measurement protocol is limited because of the finite electron coherence time T₂. As discussed in Sect. 2.3, this can be overcome by three- or five-pulse ESEEM.

Figure 11 (a) Three-pulse (blue circles) and five-pulse (green circles) ESEEM signals of Bi:Si donors at B₀=0.1 mT. Black lines are simulations assuming a ²⁹Si concentration of $p=\mathrm{4.4}\times {\mathrm{10}}^{-\mathrm{4}}$. (b) Imaginary and (c) real part of the Fourier transform of the five-pulse ESEEM data. The spectrum only contains a peak at 850 Hz, which is the ²⁹Si nuclei Larmor frequency at this field.

We measure three- and five-pulse ESEEM with the pulse sequence shown in Fig. 11. The high-Q configuration is chosen for which T₁=120 ms is measured (see the Supplement); shaped pulses generate an intra-cavity field in the form of a rectangular pulse of 1 µs duration with sharp rise and fall (Probst et al., 2019) despite the high resonator quality factor. The data are acquired at B₀=0.1 mT, so that ${\mathit{\omega }}_I/\mathrm{2}\mathit{\pi }\simeq \mathrm{850}$ Hz. The first blind spot for three-pulse ESEEM is thus at $\mathrm{2}\mathit{\pi }/{\mathit{\omega }}_I=\mathrm{1.2}$ ms; we chose τ=290 µs for the three-pulse echo, and ${\mathit{\tau }}_{\mathrm{1}}={\mathit{\tau }}_{\mathrm{2}}=\mathrm{290}$ µs for the five-pulse sequence. A sequence of 19 CPMG π pulses, separated by 50 µs, was used to enhance the signal-to-noise ratio. The sequences were repeated after a fixed waiting time of 100 ms between the last π pulse of one sequence and the first π∕2 pulse of the following, to enable spin relaxation. Phase cycling is used to suppress unwanted echoes (see the Supplement for the schemes; Schweiger and Jeschke, 2001; Kasumaj and Stoll, 2008). Each point is averaged over 2.5×10⁴ sequences, with a total acquisition time of 2 weeks for each curve (Probst et al., 2020).

The results are shown in Fig. 11, together with their fast Fourier transform (Probst et al., 2020). Both the three-pulse ESEEM (3PE) and five-pulse ESEEM (5PE) curves show oscillations that last 1 order of magnitude longer than the electron spin T₂ (up to 20 ms), enabling higher spectral resolution of the ESEEM signal. The 5PE curve has a higher oscillation amplitude than the 3PE by a factor of 2–3, as expected. The decay of the oscillations occurs in ∼10 ms, 1 order of magnitude faster than the stimulated echo amplitude (see the 3PE curve), suggesting that it is an intrinsic feature of the ESEEM signal, as discussed below.

The spectrum shows only one peak at the ²⁹Si frequency. This is consistent with the expression provided in Sect. 2 and the Supplement for the three- and five-pulse ESEEM, in which the terms oscillating at the sum and difference frequency are absent in contrast to the two-pulse ESEEM. The peak width is ≃100 Hz, which indicates that the nuclei contributing to the ESEEM signal have hyperfine coupling strengths A, B of at most 100 Hz. Neglecting the contact interaction term, this corresponds to ²⁹Si nuclei that are located at least ∼5 nm away from the donor spin.

The measured ESEEM spectrum of the bismuth donor sample qualitatively differs from the erbium sample, since it only contains a peak at the unperturbed silicon nuclei Larmor frequency (and at twice this frequency for the two-pulse ESEEM), instead of the many peaks observed in Fig. 7 indicating nuclear-spin contribution with vastly different hyperfine strengths. This can be qualitatively understood by examining Eq. (8). Defining N_l as the number of lattice sites with approximately the same hyperfine parameters A_l, B_l and modulation frequency ${\mathit{\omega }}_{\downarrow /\uparrow ,l}$, the component at ${\mathit{\omega }}_{\downarrow /\uparrow ,l}$ is visible in the spectrum if N_lk_lp∼1, which can only be achieved if N_lp∼1. In the case of erbium, p=0.144 so that even the sites closest to the ion (for which N_l is of order unity) may satisfy this condition for well-chosen B₀. In the bismuth donor sample where $p=\mathrm{4.4}\times {\mathrm{10}}^{-\mathrm{4}}$, this condition can only be met for N_l∼10³, and therefore for crystal sites l that are far from the donor, for which the hyperfine coupling is small, so that ${\mathit{\omega }}_{\downarrow /\uparrow ,l}\simeq {\mathit{\omega }}_I$. This is confirmed by the more quantitative modeling below.

5.2.4 Comparison with the model

As explained above, the measured echo signal results from the contribution of all 18 Bi:Si transitions because of strain broadening. To model the data, we therefore apply the fictitious spin-1∕2 model to each transition, and sum the resulting echo amplitudes weighted by their relative contribution, which we determine using numerical simulations described in the Supplement.

Moreover, as discussed in Sect. 3.2, and in contrast to the erbium case, the fictitious spin model for a given transition needs to be validated in the low-B₀ regime because the energy difference between neighboring hyperfine levels of the bismuth donor manifold ($E_m^{\pm }-E_{m-\mathrm{1}}^{\pm })/h\simeq \mathrm{0.3}$ MHz for B₀=0.1 mT is comparable to or even lower than the hyperfine coupling to some ²⁹Si nuclei. In that case, the hyperfine interaction induces significant mixing between the bismuth donor and the ²⁹Si eigenstates, and we should describe the coupled “electron spin S₀+²⁰⁹Bi nuclear spin I₀+²⁹Si nuclear spin I” system as a single 40-level quantum system.

This study is described in Sect. S5 for a ²⁹Si with strong hyperfine coupling (≥200 kHz). The state mixing makes many transitions EPR allowed, and the interference between these transitions causes fast oscillations in the spin-echo signal, as seen in Fig. S7 in the Supplement. The frequencies of these oscillations depend greatly on the local Overhauser field on the donor electron spin. Since the latter has a large inhomogeneous broadening (∼0.5 MHz), the ensemble average leads to a rapid decay of the signal (<1 µs). Given the ²⁹Si concentration, about 10 % of the donors have one or more ²⁹Si with coupling >300 kHz in the proximity, which therefore leads to a rapid decay of the total echo signal within ∼1 µs by about 10 %. In the experimental data, this fast decay is not visible because the echo signal is measured at longer times, and therefore the ESEEM signals presented in Fig. S5 are those from ²⁹Si with couplings <200 kHz.

As for spins with a coupling strength between 20 and 200 kHz, they lead to ESEEM amplitude much less than 1 %, as shown in Figs. S7–S9. For nuclear spins with a hyperfine coupling <100 kHz, the fictitious spin model produces results with negligible errors of the modulation frequencies from the exact solution (Figs. S5 and S6). Furthermore, the systematic numerical studies (Figs. S9) show that a nearby Si nuclear spin with coupling <100 kHz has little effect on the ESEEM due to other distant nuclear spins.

Considering these different contributions of Si nuclear spins of different hyperfine couplings, as discussed in the paragraph above and in more details in the Supplement, we apply the fictitious spin-1∕2 model to each EPR-allowed transition of the bismuth donor manifold, considering only Si nuclear spins that have a hyperfine coupling weaker than a certain cut-off which we choose as 20 kHz, and discarding all the others.

For each transition, we compute the hyperfine parameters that enter the fictitious spin-1∕2 model for all sites of the silicon lattice. We then generate a large number of random configurations of nuclear spins. We compute the corresponding two-, three-, or five-pulse ESEEM signal using the analytical formulas of Sect. 2.4 after discarding all nuclei whose hyperfine coupling is larger than 20 kHz. We average the signal for one configuration over all bismuth donor transitions using the weights determined by simulation and then average the results over all the configurations computed. In this way, we obtain the curves shown in Fig. 9.

We use the two-pulse echo dataset to determine the most likely sample concentration in ²⁹Si, using p as a fitting parameter. As seen in Fig. 9b, the best fit is obtained for $p=\mathrm{4.4}\pm \mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$, which is compatible with the specified $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$. The agreement is satisfactory but not perfect, as seen, for instance, in the amplitude of the short-time ESEEM oscillations which are lower in the measurements than in the simulations, particularly at larger field. Also, the peak at 2ω_I is notably broader and has a lower amplitude than in the experiment.

For the fitted value of p, the three- and five-pulse theoretical signals are also computed and found to be in overall agreement with the data, even though the decay of the ESEEM signal predicted by the model is faster than in the experiment and correspondingly the predicted ESEEM spectrum broader than the data.

6 Discussion and conclusion

We have reported two-, three-, and five-pulse ESEEM measurements using a quantum-limited EPR spectrometer on two model systems: erbium ions in a CaWO₄ matrix and bismuth donors in silicon. Whereas the erbium measurements are done in a commonly used regime of high field, the bismuth donor measurements are performed in an unusual regime of low nuclear-spin density, low hyperfine coupling, and almost zero magnetic field. Good agreement is found with the simplest analytical ESEEM models.

Having demonstrated that ESEEM is feasible in a millikelvin quantum-limited EPR spectrometer setup on two model spin systems, it is worth speculating in broader terms about its potential for real-world hyperfine spectroscopy. First, high magnetic fields are desirable for a better spectral resolution. Superconducting resonators in Nb, niobium nitride (NbN), or niobium titanium nitride (NbTiN) can retain a high quality factor up to ∼1 T (Graaf et al., 2012; Samkharadze et al., 2016; Mahashabde et al., 2020), so that quantum-limited EPR spectroscopy in the Q band can in principle be envisioned. Resonator bandwidths larger than demonstrated here are also desirable. Given, increasing κ in the Purcell regime leads to longer relaxation times T₁, this should be done with care. One option is to increase also the coupling constant g by further reduction of the resonator mode volume (Ranjan et al., 2020b). Interestingly, this provides another motivation to apply higher magnetic fields, since g is proportional to ω_r. Overall, a resonator at ${\mathit{\omega }}_{\mathrm{r}}/\mathrm{2}\mathit{\pi }\simeq \mathrm{30}$ GHz, in a magnetic field B₀≃1 T, and with a $\mathit{\kappa }/\mathrm{2}\mathit{\pi }\sim \mathrm{10}$ MHz bandwidth seems within reach, while keeping the Purcell T₁ well below 1 s. Such a high-bandwidth, high-sensitivity EPR spectrometer would be ideally suited for studying surface defects. One potential concern, however, is the power-handling capability of the resonator, as the kinetic inductance causes a non-linear response at high power.