Coherent electrical control of a single high-spin nucleus in silicon

TL;DR

This paper demonstrates the first experimental coherent control of a single high-spin nucleus in silicon using localized electric fields, enabling scalable quantum information processing without oscillating magnetic fields.

Contribution

It introduces a novel all-electrical control method for a single high-spin nucleus in silicon, supported by a microscopic theoretical model.

Findings

Achieved nuclear spin dephasing time of 0.1 seconds.

Controlled a single antimony nucleus using electric fields.

Theoretical model explains quadrupole interaction modulation.

Abstract

Nuclear spins are highly coherent quantum objects. In large ensembles, their control and detection via magnetic resonance is widely exploited, e.g. in chemistry, medicine, materials science and mining. Nuclear spins also featured in early ideas and demonstrations of quantum information processing. Scaling up these ideas requires controlling individual nuclei, which can be detected when coupled to an electron. However, the need to address the nuclei via oscillating magnetic fields complicates their integration in multi-spin nanoscale devices, because the field cannot be localized or screened. Control via electric fields would resolve this problem, but previous methods relied upon transducing electric signals into magnetic fields via the electron-nuclear hyperfine interaction, which severely affects the nuclear coherence. Here we demonstrate the coherent quantum control of a single antimony (spin-7/2) nucleus, using localized electric fields produced within a silicon nanoelectronic device. The method exploits an idea first proposed in 1961 but never realized experimentally with a single nucleus. Our results are quantitatively supported by a microscopic theoretical model that reveals how the purely electrical modulation of the nuclear electric quadrupole interaction, in the presence of lattice strain, results in coherent nuclear spin transitions. The spin dephasing time, 0.1 seconds, surpasses by orders of magnitude those obtained via methods that require a coupled electron spin for electrical drive. These results show that high-spin quadrupolar nuclei could be deployed as chaotic models, strain sensors and hybrid spin-mechanical quantum systems using all-electrical controls. Integrating electrically controllable nuclei with quantum dots could pave the way to scalable nuclear- and electron-spin-based quantum computers in silicon that operate without the need for oscillating magnetic fields.

Figures

B). We use Eq. ([36](#E36)) to relate the computed strains to EFG and numerically compute the eigenenergies of Eq. [2](#E2) and show the estimated quadrupole splitting $f_{Q}$ in the main text (Fig. 4C). The variation in quadrupole splitting shows a good qualitative correspondence with the variation of the shear strain, reflecting the larger contribution of $S_{44}$ to the quadrupole splitting. Values found for the strain-induced quadrupole splitting range between 90 and $160\text{\,}\mathrm{kHz}$ throughout the area where the donor is estimated to be. This result represents an outstanding qua

G) on the $m_{I}=+5/2\leftrightarrow+7/2$ ($\Delta m_{I}=\pm 1$) transition to extract the pure dephasing time $T_{2n+}^{*}(+5/2\leftrightarrow+7/2)=92(8)$ ms, which corresponds to an NER broadening (full width at half maximum) $\Gamma_{n}=\ln{2}/(\pi T_{2n+}^{*})=2.4(2)$ Hz. The $m_{I}=+3/2\leftrightarrow+7/2$ ($\Delta m_{I}=\pm 2$) transition has shorter $T_{2n+}^{*}(+3/2\leftrightarrow+7/2)=28(1)$ ms (Fig. 2H). Both values, while extremely long in absolute terms, are noticeably shorter than the $T_{2n+}^{*}\approx 250-600$ ms measured on the *31*P nucleus in two other similar devices [[16](

C,D) from the nuclear quadrupole interaction requires breaking the $\mathrm{T}_{d}$ symmetry of the silicon crystal. This can be caused by uniaxial strain ($\epsilon_{zz}$), which lowers the symmetry to $\mathrm{D}_{2d}$, or shear strain ($\epsilon_{xx}-\epsilon_{zz}$), which lowers it to $\mathrm{C}_{2v}$. Even in the absence of strain, an electric field can break the $\mathrm{T}_{d}$ symmetry by polarizing the atomic bonds.

B). We use Eq. ([36](#E36)) to relate the computed strains to EFG and numerically compute the eigenenergies of Eq. [2](#E2) and show the estimated quadrupole splitting $f_{Q}$ in the main text (Fig. 4C). The variation in quadrupole splitting shows a good qualitative correspondence with the variation of the shear strain, reflecting the larger contribution of $S_{44}$ to the quadrupole splitting. Values found for the strain-induced quadrupole splitting range between 90 and $160\text{\,}\mathrm{kHz}$ throughout the area where the donor is estimated to be. This result represents an outstanding qua

G) on the $m_{I}=+5/2\leftrightarrow+7/2$ ($\Delta m_{I}=\pm 1$) transition to extract the pure dephasing time $T_{2n+}^{*}(+5/2\leftrightarrow+7/2)=92(8)$ ms, which corresponds to an NER broadening (full width at half maximum) $\Gamma_{n}=\ln{2}/(\pi T_{2n+}^{*})=2.4(2)$ Hz. The $m_{I}=+3/2\leftrightarrow+7/2$ ($\Delta m_{I}=\pm 2$) transition has shorter $T_{2n+}^{*}(+3/2\leftrightarrow+7/2)=28(1)$ ms (Fig. 2H). Both values, while extremely long in absolute terms, are noticeably shorter than the $T_{2n+}^{*}\approx 250-600$ ms measured on the *31*P nucleus in two other similar devices [[16](

C,D) from the nuclear quadrupole interaction requires breaking the $\mathrm{T}_{d}$ symmetry of the silicon crystal. This can be caused by uniaxial strain ($\epsilon_{zz}$), which lowers the symmetry to $\mathrm{D}_{2d}$, or shear strain ($\epsilon_{xx}-\epsilon_{zz}$), which lowers it to $\mathrm{C}_{2v}$. Even in the absence of strain, an electric field can break the $\mathrm{T}_{d}$ symmetry by polarizing the atomic bonds.

b,d). This, however, does not yet break the $\mathrm{T}_{d}$ symmetry. An EFG is obtained by further introducing strain. The $S$-tensor that links EFG to strain has two unique components, $S_{11}$ (uniaxial) and $S_{44}$ (shear). We conducted a first-principles, density functional theory (DFT) calculation to extract $S_{11}=2.4\times 10^{22}$ V/m2 and $S_{44}=6.1\times 10^{22}$ V/m2 [[14](#bib.bib14)]. With a finite-element numerical model we computed the strain profile in our device, as caused by the different thermal expansions of Si and Al upon cooling to cryogenic temperatures [[19](#bib.b

Tables

Table 1. TABLE S1.: Rabi frequency coefficients for NER and NMR transitions. Transitions are labeled by their secondary spin quantum number m I subscript 𝑚 𝐼 m_{I} . The coefficient α m I − 1 ↔ m I subscript 𝛼 ↔ subscript 𝑚 𝐼 1 subscript 𝑚 𝐼 \alpha_{m_{I}-1\leftrightarrow m_{I}} governs the NER transition between neighboring levels with Δ ​ m I = ± 1 Δ subscript 𝑚 𝐼 plus-or-minus 1 \Delta m_{I}=\pm 1 , β m I − 2 ↔ m I subscript 𝛽 ↔ subscript 𝑚 𝐼 2 subscript 𝑚 𝐼 \beta_{m_{I}-2\leftrightarrow m_{I}} governs the NER transition between next nearest neighboring levels with Δ ​ m I = ± 2 Δ subscript 𝑚 𝐼 plus-or-minus 2 \Delta m_{I}=\pm 2 (hence the last transition in the table not existing), and ζ m I − 1 ↔ m I subscript 𝜁 ↔ subscript 𝑚 𝐼 1 subscript 𝑚 𝐼 \zeta_{m_{I}-1\leftrightarrow m_{I}} governs the NMR transition between neighboring levels with Δ ​ m I = ± 1 Δ subscript 𝑚 𝐼 plus-or-minus 1 \Delta m_{I}=\pm 1 .

$m_I$	$7/2$	$5/2$	$3/2$	$1/2$	$-1/2$	$-3/2$	$-5/2$
${\alpha }_{m_I-1\leftrightarrow m_I}$	$\sqrt{63}$	$\sqrt{48}$	$\sqrt{15}$	0	$\sqrt{15}$	$\sqrt{48}$	$\sqrt{63}$
${\beta }_{m_I-2\leftrightarrow m_I}$	$\sqrt{21}/2$	$\sqrt{45}/2$	$\sqrt{60}/2$	$\sqrt{60}/2$	$\sqrt{45}/2$	$\sqrt{21}/2$	-
${\zeta }_{m_I-1\leftrightarrow m_I}$	$\sqrt{7}/2$	$\sqrt{12}/2$	$\sqrt{15}/2$	$\sqrt{16}/2$	$\sqrt{15}/2$	$\sqrt{12}/2$	$\sqrt{7}/2$

Table 2. TABLE S2.: Simulation parameters used in the COMSOL device model. Parameters taken directly from the COMSOL material library are not listed here. The top section of the table lists the different thermal expansion coefficients. Thermal deformation in the two cooldown steps is captured by the average coefficient of thermal expansion in this temperature range for all materials present. The bottom section lists all applied gate voltages, see Fig.  11 for labeling of the gates.

Parameter[Source]	Symbol	Value
Thermal expansion $850\mathrm{°C}$ to $400\mathrm{°C}$
Silicon [?]	$\overline{\alpha }$	$4.198\text{×}{10}^{-6}{K}^{-1}$
Silicon oxide [?]	$\overline{\alpha }$	$0.4\text{×}{10}^{-6}{K}^{-1}$
$\alpha$-Quartz [?]	${\overline{\alpha }}_{\perp }$	$23.62\text{×}{10}^{-6}{K}^{-1}$
	${\overline{\alpha }}_{\parallel }$	$13.33\text{×}{10}^{-6}{K}^{-1}$
Thermal expansion $400\mathrm{°C}$ to $200\mathrm{mK}$
Silicon [?]	$\overline{\alpha }$	$2.28\text{×}{10}^{-6}{K}^{-1}$
Silicon oxide [?]	$\overline{\alpha }$	$0.446\text{×}{10}^{-6}{K}^{-1}$
$\alpha$-Quartz [?,?]	${\overline{\alpha }}_{\perp }$	$14.186\text{×}{10}^{-6}{K}^{-1}$
	${\overline{\alpha }}_{\parallel }$	$7.824\text{×}{10}^{-6}{K}^{-1}$
Aluminum [?,?]	$\overline{\alpha }$	$21.43\text{×}{10}^{-6}{K}^{-1}$
Aluminum oxide [?]	$\overline{\alpha }$	$4.98\text{×}{10}^{-6}{K}^{-1}$
Applied gate voltages
Top gate	$V_{\mathrm{TG}}$	$1.881V$
Right barrier	$V_{\mathrm{RB}}$	$0.3421V$
Left barrier	$V_{\mathrm{LB}}$	$0.4293V$
Donor gate front left	$V_{\mathrm{DFL}}$	$0.4823V$
Donor gate front right	$V_{\mathrm{DFR}}$	$0.3915V$
Donor gate back left	$V_{\mathrm{DBL}}$	$0.4177V$
Donor gate back right	$V_{\mathrm{DBR}}$	$0.4059V$
Plunger gate	$V_{\mathrm{PL}}$	$-0.0032V$
Source	$V_{\mathrm{SRC}}$	$0.0002V$

Equations

$$\hat{\mathcal{H}}/h={\gamma }_n{B}_0{\hat{I}}_z+\sum _{\alpha ,\beta \in \{x,y,z\}}{Q}_{\alpha \beta }{\hat{I}}_{\alpha }{\hat{I}}_{\beta },$$

$${\hat{\mathcal{H}}}_{m_I-1\leftrightarrow m_I}^{\mathrm{NER}}\left(t\right)/h=$$

$$+\delta Q_{yz}\left({\hat{I}}_y{\hat{I}}_z+{\hat{I}}_z{\hat{I}}_y\right)],$$

$${\hat{\mathcal{H}}}_{m_I-2\leftrightarrow m_I}^{\mathrm{NER}}\left(t\right)/h=$$

$$+\delta Q_{xy}\left(I_x{I}_y+I_y{I}_x\right)].$$

$$\hat{\mathcal{H}}={\gamma }_e{B}_0{\hat{S}}_z+{\gamma }_n{B}_0{\hat{I}}_z+A\vec{S}\cdot \vec{I}+{\hat{\mathcal{H}}}_Q,$$

$$\hat{\mathcal{H}}={\gamma }_n{B}_0{\hat{I}}_z+{\hat{\mathcal{H}}}_Q.$$

$${\hat{\mathcal{H}}}_Q=Q_{xx}{\hat{I}}_x^2+Q_{yy}{\hat{I}}_y^2+Q_{zz}{\hat{I}}_z^2+Q_{yz}({\hat{I}}_y{\hat{I}}_z+{\hat{I}}_z{\hat{I}}_y)+Q_{xz}({\hat{I}}_x{\hat{I}}_z+{\hat{I}}_z{\hat{I}}_x)+Q_{xy}({\hat{I}}_x{\hat{I}}_y+{\hat{I}}_y{\hat{I}}_x),$$

$$Q_{\alpha \beta }=\frac{e{q}_n{\mathcal{V}}_{\alpha \beta }}{2I(2I-1)h},$$

$$E_{m_I}={\gamma }_n{B}_0{m}_I+⟨m_I⟩{\hat{\mathcal{H}}}_Q{m}_I,$$

$$f_{m_I-1\leftrightarrow m_I}=-{\gamma }_n{B}_0+\left(m_I-\frac{1}{2}\right)\left(Q_{xx}+Q_{yy}-2Q_{zz}\right),$$

$$f_Q=f_{m_I\leftrightarrow m_I+1}-f_{m_I-1\leftrightarrow m_I}=Q_{xx}+Q_{yy}-2Q_{zz}.$$

$${\hat{\mathcal{H}}}_Q^{\prime }=\frac{e{q}_n{\mathcal{V}}_{z^{\prime }{z}^{\prime }}}{4I(2I-1)h}\left[3{\hat{I}}_{z^{\prime }}^2-I^2+\frac{{\mathcal{V}}_{x^{\prime }{x}^{\prime }}-{\mathcal{V}}_{y^{\prime }{y}^{\prime }}}{{\mathcal{V}}_{z^{\prime }{z}^{\prime }}}\left({\hat{I}}_{x^{\prime }}^2-{\hat{I}}_{y^{\prime }}^2\right)\right],$$

$$\hat{\mathcal{H}}(t)=\cos (2\pi ft)\delta {\hat{\mathcal{H}}}_Q=\cos (2\pi ft){\sum }_{\alpha ,\beta \in \{x,y,z\}}\delta Q_{\alpha \beta }{\hat{I}}_{\alpha }{\hat{I}}_{\beta }.$$

$${\hat{I}}_x$$

$${\hat{I}}_y$$

$${\hat{I}}_z$$

$$f_{m_I-1\leftrightarrow m_I}^{\text{Rabi, NER}}=|⟨m_I-1⟩\delta {\hat{\mathcal{H}}}_Q{m}_I|=|\delta Q_{xz}⟨m_I-1⟩{\hat{I}}_x{\hat{I}}_z+{\hat{I}}_z{\hat{I}}_x{m}_I+\delta Q_{yz}⟨m_I-1⟩{\hat{I}}_y{\hat{I}}_z+{\hat{I}}_z{\hat{I}}_y{m}_I|.$$

$${\hat{I}}_x{\hat{I}}_z+{\hat{I}}_z{\hat{I}}_x=\left[\begin{array}{cccccccc}0& \sqrt{63}& 0& 0& 0& 0& 0& 0\\ \sqrt{63}& 0& \sqrt{48}& 0& 0& 0& 0& 0\\ 0& \sqrt{48}& 0& \sqrt{15}& 0& 0& 0& 0\\ 0& 0& \sqrt{15}& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& -\sqrt{15}& 0& 0\\ 0& 0& 0& 0& -\sqrt{15}& 0& -\sqrt{48}& 0\\ 0& 0& 0& 0& 0& -\sqrt{48}& 0& -\sqrt{63}\\ 0& 0& 0& 0& 0& 0& -\sqrt{63}& 0\end{array}\right],$$

$$f_{m_I-1\leftrightarrow m_I}^{\text{Rabi, NER}}={\alpha }_{m_I-1\leftrightarrow m_I}\sqrt{\delta Q_{xz}^2+\delta Q_{yz}^2}$$

$${\alpha }_{m_I-1\leftrightarrow m_I}=\frac{1}{2}\left|2m_I-1\right|\sqrt{I(I+1)-m_I(m_I-1)}.$$

$$f_{m_I-2\leftrightarrow m_I}^{\text{Rabi, NER}}=|\delta Q_{xx}⟨m_I-2⟩{\hat{I}}_x^2{m}_I+\delta Q_{yy}⟨m_I-2⟩{\hat{I}}_y^2{m}_I+\delta Q_{xy}⟨m_I-2⟩{\hat{I}}_x{\hat{I}}_y+{\hat{I}}_y{\hat{I}}_x{m}_I|.$$

$${\hat{I}}_x^2=\frac{1}{4}\left[\begin{array}{cccccccc}7& 0& \sqrt{84}& 0& 0& 0& 0& 0\\ 0& 19& 0& \sqrt{180}& 0& 0& 0& 0\\ \sqrt{84}& 0& 27& 0& \sqrt{240}& 0& 0& 0\\ 0& \sqrt{180}& 0& 31& 0& \sqrt{240}& 0& 0\\ 0& 0& \sqrt{240}& 0& 31& 0& \sqrt{180}& 0\\ 0& 0& 0& \sqrt{240}& 0& 27& 0& \sqrt{84}\\ 0& 0& 0& 0& \sqrt{180}& 0& 19& 0\\ 0& 0& 0& 0& 0& \sqrt{84}& 0& 7\end{array}\right],$$

$$f_{m_I-2\leftrightarrow m_I}^{\text{Rabi, NER}}={\beta }_{m_I-2\leftrightarrow m_I}\sqrt{(\delta Q_{xx}-\delta Q_{yy}{)}^2+4\delta Q_{xy}^2}$$

$${\beta }_{m_I-2\leftrightarrow m_I}=\frac{1}{4}\sqrt{(I-m_I-7)(I-m_I-6)(I-m_I+1)(I-m_I+2)}$$

$$f_{m_I-1\leftrightarrow m_I}^{\text{Rabi, NMR}}=|⟨m_I-1⟩{\gamma }_n{B}_1{\hat{I}}_x{m}_I|,$$

$$f_{m_I-1\leftrightarrow m_I}^{\text{Rabi, NMR}}=|{\gamma }_n|B_1{\zeta }_{m_I-1\leftrightarrow m_I}$$

$${\zeta }_{m_I-1\leftrightarrow m_I}=\frac{1}{2}\sqrt{I(I+1)+m_I(m_I-1)}.$$

$$P_{\mathrm{flip}}=\frac{{\gamma }^2}{{\Omega }^2}{\sin }^2\Omega t,$$

$$\frac{dV({\vec{r}}_0;V_{\mathrm{PL}},V_{\mathrm{DFL}},...)}{d{V}_{\mathrm{PL}}}\cdot \Delta V_{\mathrm{PL}}+\frac{dV({\vec{r}}_0;V_{\mathrm{PL}},V_{\mathrm{DFL}},...)}{d{V}_{\mathrm{DFL}}}\cdot \Delta V_{\mathrm{DFL}}=0$$

$$\displaystyle-\frac{\Delta V_{\rm PL}}{\Delta V_{\rm DFL}}=\frac{dV(\vec{r}_{0};V_{\rm PL},V_{\rm DFL},...)}{dV_{\rm DFL}}\bigg{/}\frac{dV(\vec{r}_{0};V_{\rm PL},V_{\rm DFL},...)}{dV_{\rm PL}}.$$

$$s_{\mathrm{DBL}}=-\Delta V_{\mathrm{PL}}/\Delta V_{\mathrm{DBL}}=\tan (\theta ),$$

$$\displaystyle\sigma_{\rm DFL}=(s_{\rm DFL}^{2}+1)\cdot 2\pi\cdot$1\text{\,}\mathrm{\SIUnitSymbolDegree}$/$360\text{\,}\mathrm{\SIUnitSymbolDegree}$.$$

Full text

Nuclear spins are highly coherent quantum objects. In large ensembles, their control and detection via magnetic resonance is widely exploited, e.g. in chemistry, medicine, materials science and mining. Nuclear spins also featured in early ideas [1] and demonstrations [2] of quantum information processing. Scaling up these ideas requires controlling individual nuclei, which can be detected when coupled to an electron [3, 4, 5]. However, the need to address the nuclei via oscillating magnetic fields complicates their integration in multi-spin nanoscale devices, because the field cannot be localized or screened. Control via electric fields would resolve this problem, but previous methods [6, 7, 8] relied upon transducing electric signals into magnetic fields via the electron-nuclear hyperfine interaction, which severely affects the nuclear coherence. Here we demonstrate the coherent quantum control of a single antimony (spin-7/2) nucleus, using localized electric fields produced within a silicon nanoelectronic device. The method exploits an idea first proposed in 1961 [9] but never realized experimentally with a single nucleus. Our results are quantitatively supported by a microscopic theoretical model that reveals how the purely electrical modulation of the nuclear electric quadrupole interaction, in the presence of lattice strain, results in coherent nuclear spin transitions. The spin dephasing time, 0.1 seconds, surpasses by orders of magnitude those obtained via methods that require a coupled electron spin for electrical drive. These results show that high-spin quadrupolar nuclei could be deployed as chaotic models, strain sensors and hybrid spin-mechanical quantum systems using all-electrical controls. Integrating electrically controllable nuclei with quantum dots [10, 11] could pave the way to scalable, nuclear- and electron-spin-based quantum computers in silicon that operate without the need for oscillating magnetic fields.

Nuclear Magnetic Resonance (NMR) relies on the presence of a static magnetic field, $B_{0}$, that separates the energy levels of the nuclear spins, and a radio-frequency (RF) oscillating magnetic field, $B_{1}$, that induces transitions between such levels. Magnetic fields cannot easily be confined or screened at the nanoscale. Therefore, identical nuclear spins within large regions would all respond to the same signal, preventing the spins from being individually addressed. Electric fields, instead, can be efficiently routed and confined within highly complex nanoscale devices, a prime example being the sophisticated interconnects found in modern silicon computer chips. These observations suggest that an ideal route to scale up nuclear spin based quantum devices would involve the use of RF electric fields for spin control.

A theoretical idea in this direction was proposed by Bloembergen as early as 1961 [9]: for nuclei with spin $I>1/2$ and nonzero electric quadrupole moment $q_{n}$, a resonant electric field induces nuclear spin transitions by modulating the nuclear quadrupole interaction, if the nuclei are placed in solids that lack point inversion symmetry at the lattice site. In bulk ensembles, the static shift of the NMR frequency by a DC electric field, named Linear Quadrupole Stark Effect (LQSE), was observed in the 1960s [12]. The resonant version of LQSE, called Nuclear Electric Resonance (NER) was only demonstrated more recently [13] in a bulk GaAs crystal.

We report here the demonstration of NER and coherent electrical control of a single 123Sb nucleus in silicon. The discovery that this nucleus could be electrically controlled was, in fact, fortuitous. The 123Sb atom possesses a nuclear spin $I=7/2$ with electric quadrupole moment $q_{n}=-0.69$ barn. Depending on its electrochemical potential relative to a nearby electron reservoir, an electron (with spin $S=1/2$) can optionally be bound to the nucleus. The atom was implanted in a metal-oxide-semiconductor nanostructure [14] (Fig. 1A) fabricated on isotopically-enriched 28Si, similar to those developed for 31P spin qubits [4, 15, 16]. The structure comprises a single-electron transistor for electron spin readout [17], electrostatic gates to control the donor electrochemical potential, and a broadband on-chip microwave antenna [18]. The antenna is nominally terminated by a short circuit, in order to obtain maximum current at its tip and produce strong oscillating magnetic fields to control both the electron (at $\approx 40$ GHz) and the nuclear (at $\approx 10$ MHz) spins of the donor. In this device, however, an electrostatic discharge damaged the short circuit termination (Fig. 1A). The small gap in the termination has a low enough impedance at 40 GHz to allow current flow for electron spin resonance, but at $\approx 10$ MHz it produces solely an RF electric field. Once we realized that NER was possible, we changed to using the electrostatic gates fabricated right above the donor, which have an even stronger effect.

We focus here on the 123Sb donor in its ionized state; the removal of the donor-bound electron precludes any interpretation of the data involving modulation of hyperfine fields [6, 8].

In nanoscale Si devices, the aluminum gates can cause significant lattice strain at low temperatures, due to the different thermal contraction of Al and Si [19]. Lattice strain creates an electric field gradient (EFG) $\mathcal{V}_{\alpha\beta}=\partial^{2}V/\partial\alpha\partial\beta$ ($V$ is the electric potential and $\alpha,\beta=x,y,z$) at the nuclear site [20, 21] (Fig. 1B), which produces a static nuclear quadrupole interaction $Q_{\alpha\beta}=eq_{n}\mathcal{V}_{\alpha\beta}/[2I(2I-1)]$, resulting in a quadrupole splitting $f_{Q}$ of the nuclear resonance frequencies (Fig. 1D), making all transitions individually addressable.

The application of an RF electric field of amplitude $E_{1}$ modulates the nuclear quadrupole energies by $\delta Q_{xz}$ and $\delta Q_{yz}$, and induces transitions between nuclear states at a rate $f^{\textrm{Rabi,NER}}_{m_{I}-1\leftrightarrow m_{I}}\propto|\delta Q_{xz}\langle m_{I}-1|\hat{I}_{x}\hat{I}_{z}+\hat{I}_{z}\hat{I}_{x}|m_{I}\rangle|$. Notably, the transition rate is predicted to be zero for the $m_{I}=-1/2\leftrightarrow+1/2$ transition [see Eq. 15 in [14]]. Figure 2A shows the experimental NER spectrum for $\Delta m=\pm 1$ transitions, containing six sharp resonances separated by $f_{Q}=66$ kHz, with the $m_{I}=-1/2\leftrightarrow+1/2$ absent.

The quadrupole interaction is quadratic in the spin operators. Therefore, transitions between spin states that differ by $\Delta m_{I}=\pm 2$ are allowed, to first-order, and occur at a rate $f^{\textrm{Rabi, NER}}_{m_{I}-2\leftrightarrow m_{I}}\propto|\delta Q_{xx}\langle m_{I}-2|\hat{I}_{x}^{2}|m_{I}\rangle|$ [see Eq. 19 in [14]]. Importantly, all $\Delta m_{I}=\pm 2$ transitions have nonzero rate (Fig. 2B). This allowed us to “jump over” the NER-forbidden $m_{I}=-1/2\leftrightarrow+1/2$ transition, reaching the $m_{I}=-1/2$ state via the $m_{I}=-1/2\leftrightarrow+3/2$ transition. From there, the 3 remaining $\Delta m_{I}=\pm 1$ transitions between states with negative $m_{I}$ could be accessed.

Figure 2C,D presents the observed transition rates between each pair of states, in excellent agreement with the predicted trends from NER theory. For the $\Delta m_{I}=\pm 1$ transitions, the NMR Rabi frequencies would follow $f^{\textrm{Rabi,NMR}}_{m_{I}-1\leftrightarrow m_{I}}=|\gamma_{n}B_{1}\langle m_{I}-1|\hat{I}_{x}|m_{I}\rangle|$, which is notably maximal for $m_{I}=-1/2\leftrightarrow+1/2$ transition. The $\Delta m_{I}=\pm 2$ NMR transitions are forbidden, to first-order. These results prove decisively that our experiments do not constitute a form of magnetic resonance.

As observed in earlier experiments on 31P [22, 16], the nuclear spins of ionized donors in 28Si have exceptional quantum coherence properties. We performed a Ramsey experiment (Fig. 2G) on the $m_{I}=+5/2\leftrightarrow+7/2$ ($\Delta m_{I}=\pm 1$) transition to extract the pure dephasing time $T_{2n+}^{*}(+5/2\leftrightarrow+7/2)=92(8)$ ms, which corresponds to an NER broadening (full width at half maximum) $\Gamma_{n}=\ln{2}/(\pi T_{2n+}^{*})=2.4(2)$ Hz. The $m_{I}=+3/2\leftrightarrow+7/2$ ($\Delta m_{I}=\pm 2$) transition has shorter $T_{2n+}^{*}(+3/2\leftrightarrow+7/2)=28(1)$ ms (Fig. 2H). Both values, while extremely long in absolute terms, are noticeably shorter than the $T_{2n+}^{*}\approx 250-600$ ms measured on the 31P nucleus in two other similar devices [16], fabricated on the same 28Si wafer. Since the 31P nucleus has zero quadrupole moment, this suggests that the 123Sb coherence may be affected by electrical noise [23], in a way that the 31P is not. Therefore, the 123Sb nucleus could become a useful tool for spectroscopy of very slow electrical noise. Nonetheless, our dephasing time remains 2 orders of magnitude longer than that observed in 31P when adding a hyperfine-coupled electron, $T_{2n0}^{*}\approx 430-570$ $\mu$s [16], and 3 orders of magnitude longer than the $T_{2}^{*}=64$ $\mu$s observed in a Tb nucleus in a single-atom magnet [6]. This observation highlights the benefit of a pure electrical control mechanism that does not rely upon hyperfine interactions.

We measured the Rabi frequencies of the $\Delta m_{I}=\pm 1$ and $\Delta m_{I}=\pm 2$ NER transitions as a function of the amplitude of the RF voltage applied to the gate, finding $g_{\rm E,1}=34.21(3)$ Hz/mV (Fig. 3A) and $g_{\rm E,2}=1.995(4)$ Hz/mV (Fig. 3B). These transition rates show that NER is a weak effect, but due to the long nuclear spin coherence in 28Si we were able to perform high-fidelity Rabi flops persisting for tens of milliseconds (Fig. 2E,F).

In addition to driving nuclear spin transitions with an RF voltage, we were able to Stark shift the resonance frequencies using an additional DC voltage $\Delta V$ on the gates (Fig. 3C,D). All NER frequencies shifted according to $\delta f=\Delta f_{Q}\Delta Vf(m_{I},\Delta m_{I})$, with $\Delta f_{Q}=9.9(3)$ Hz/mV and $f(m_{I},\Delta m_{I})$ a factor of order unity that represents the matrix element of the electric quadrupole interaction between the initial and final state of each transition [14].

The results reported here constitute the first observation of coherent, purely electrical control of a single nuclear spin. Achieving this in silicon is, at first sight, remarkable: no effect of electric fields on nuclear spins has ever been observed in a non-polar, non-piezoelectric material in the absence of a hyperfine-coupled electron. To gain a microscopic understanding of this phenomenon, we conjectured that our results are a form of LQSE [12]. Resonant transitions between nuclear levels induced by electric fields (NER) require that the crystal does not possess point inversion symmetry at the atomic site [9], as is indeed the case for silicon. However, the static shift of the resonance lines (Fig. 3C,D) from the nuclear quadrupole interaction requires breaking the $\mathrm{T}_{d}$ symmetry of the silicon crystal. This can be caused by uniaxial strain ($\epsilon_{zz}$), which lowers the symmetry to $\mathrm{D}_{2d}$, or shear strain ($\epsilon_{xx}-\epsilon_{zz}$), which lowers it to $\mathrm{C}_{2v}$. Even in the absence of strain, an electric field can break the $\mathrm{T}_{d}$ symmetry by polarizing the atomic bonds.

The larger and charged donor atom introduces a local lattice distortion, displacing its four coordinating Si atoms by 0.2 Å, and polarizes the charge density along the bonds (Fig. 4b,d). This, however, does not yet break the $\mathrm{T}_{d}$ symmetry. An EFG is obtained by further introducing strain. The $S$-tensor that links EFG to strain has two unique components, $S_{11}$ (uniaxial) and $S_{44}$ (shear). We conducted a first-principles, density functional theory (DFT) calculation to extract $S_{11}=2.4\times 10^{22}$ V/m2 and $S_{44}=6.1\times 10^{22}$ V/m2 [14]. With a finite-element numerical model we computed the strain profile in our device, as caused by the different thermal expansions of Si and Al upon cooling to cryogenic temperatures [19, 21] (Fig. 1B). Finally, we triangulated the most likely location of the 123Sb+ donor by combining the implantation depth profile with a modeling of the relative capacitive coupling between the donor and different pairs of control gates, extracted from the experimental charge stability diagrams [14]. Combining these three pieces of information, we arrived at a spatial map of quadrupole splittings $f_{Q}$ (Fig. 4C), which shows good agreement between the models and the experiment around the predicted location of the donor under study.

The effect of electric fields, both static (LQSE) and dynamic (NER), can be understood as arising from the single unique component of the $R$-tensor, $R_{14}$ [14]. By combining a finite-element model of the electric field in the device, the estimated 123Sb+ donor position, and the experimental values of LQSE and NER Rabi frequencies, we extracted $R_{14}=1.7\times 10^{12}$ m*-1* [14]. The strength of this coupling is comparable to prior bulk measurements of LQSE on 75As in GaAs [24]. This can be understood by observing that, while the Sb+-Si bond has a weaker ionic character compared to the Ga-As bond, $R_{14}$ scales with atomic number, leading to a similar overall value. Since our model agrees with the experiment within a factor of order unity, and no alternative explanation comes within orders of magnitude of the results [14], we conclude that we have observed the manifestation of LQSE and NER in a single nuclear spin in silicon.

Our results have significant consequences for the development of nuclear spin based quantum computers, and the design of nanoscale quantum devices. The Hilbert space of the $I=7/2$ 123Sb nucleus has 8 dimensions. It can encode the equivalent of 3 quantum bits of information, allowing simple quantum algorithms [25] or quantum error correction codes [26], all solely using electric fields. The donor electron and nuclear spins combined form a “flip-flop” qubit [10], controllable by electric dipole spin resonance. This scheme normally requires a magnetic antenna to reset the nuclear state in the appropriate qubit subspace. This need for could be removed completely by using an electrically-drivable high-spin nucleus like 123Sb. A recent result showed that lithographic quantum dots in silicon can be entangled with nuclear spins, and that the nuclear coherence can be preserved while shuttling the electron between different dots [11]. Electron spin qubits in silicon can be coherently controlled by electric fields with high speed and high fidelity [27]. Adding the ability to electrically control quadrupolar nuclei paves the way to quantum computer architectures that integrate fast electron spin qubits with long-lived nuclear quantum memories, while fully exploiting the controllability and scalability of silicon metal-oxide-semiconductor devices, without the complication of routing RF magnetic fields within the device.

The experimental validation of a microscopic model of the relation between strain and quadrupole splitting, obtained in a functional silicon electronic device, suggests the use of quadrupolar nuclei as single-atom probes of local strain, which has a key role in enhancing the performance of ultra-scaled transistors [28]. The observation of a large quadrupole splitting $f_{Q}=66$ kHz in a high-spin nucleus creates a platform in which to study quantum chaotic dynamics in a single particle [29]. This has further applications in quantum information science, e.g. because of the remarkable analogies between chaotic spin models and digital quantum simulations [30]. Although the strain in the present device is static, our work allows us to predict the nuclear Rabi frequencies that would arise from time-dependent strain [14]. A dynamical strain $\approx 5\times 10^{-8}$ would cause a Rabi frequency of 10 Hz, comparable to both the inhomogeneous nuclear linewidth $\Gamma_{n}\approx 2.4$ Hz, and to the linewidth $\Gamma_{m}$ of high-quality silicon mechanical resonators in the MHz range [31]. Therefore, it is conceivable to achieve the strong-coupling limit of cavity-quantum electrodynamics between a single nuclear spin and a macroscopic mechanical oscillator, adding a novel spin-mechanical coupling pathway to the toolbox of hybrid quantum systems for quantum information processing and precision sensing [32].

Methods

Fabrication

The device was fabricated on a <100> p-type silicon wafer, with a 900 nm thick epitaxial layer of isotopically purified 28Si on top (concentration of residual 29Si is 730 ppm). Metallic leads for the SET are formed using optical lithography and phosphorus diffusion. The substrate is subsequently covered with a 250 nm thick field oxide, with a small central window (10 x 20 $\mu$m) containing a high quality, thermally grown layer of SiO2 with a thickness of 8 nm. Using a combination of standard optical and electron beam lithography techniques, the device is fabricated on this thin oxide window. First, a small (90 x 100 nm) window is defined, through which ${}^{123}\mathrm{Sb}$ ions are implanted at an energy of 8 keV and fluence of 2 x 1011/cm2, corresponding to an average of 14 donors in the implantation window. Donors are activated using a rapid thermal anneal at 1000 C. Next, in two electron beam lithography steps, the gates forming the SET, the donor gates, and the microwave antenna are created, using thermally evaporated aluminum and lift off, with aluminum oxide as gate dielectric. Ohmics to the n-doped SET leads are formed using optical lithography, thermally evaporated aluminium and lift off, followed by a forming gas anneal. A detailed step-by-step process flow is given in [14].

Experimental setup

The sample is cooled to a temperature of  20 mK in a dilution refrigerator (Bluefors BF-LD400) with superconducting magnet. During measurements, arbitrary waveform generators (Signadyne M3201A and M3300A) were used to tune the donor electrochemical potential, generate NER pulses, and IQ-modulate the microwave signals generated by a vector microwave source (Keysight E8267D). The SET current was amplified with a transimpedance amplifier (FEMTO DLPCA-200 in combination with Stanford Instruments SIM911), and subsequently measured with a digitizer (Signadyne M3300A). Full details of the experimental setup including wiring schematic can be found in [14].

Nuclear spin readout

The nuclear spin state is measured through the electron spin state, by using ESR pulses conditional on the nuclear spin eigenstates. Electron spin readout is achieved by spin-to-charge conversion through spin-dependent tunneling onto an SET and subsequent detection of the change in charge occupation of the donor (see [14] for details).

An NER pulse has a probability $P_{\mathrm{flip}}$ of flipping the nuclear spin between two spin states. Both nuclear spin states are measured after each NER pulse, and are compared to the previous measurement. This reveals if the NER pulse has flipped the nuclear spin. This process is repeated $N_{\mathrm{repetitions}}$ times, and the resulting nuclear flips $N_{\mathrm{flips}}$ gives an estimate of the flip probability $P_{\mathrm{flip}}\approx N_{\mathrm{flips}}/(N_{\mathrm{repeititions}}-1)$.

Theoretical modelling

The spin Hamiltonian of the 123Sb nucleus takes the form:

[TABLE]

where $h=6.626\times 10^{-34}$ J/Hz is the Planck constant, $\gamma_{n}=-5.553$ MHz/T is the nuclear gyromagnetic ratio and $B_{0}=1.496$ T. $\hat{I}_{x}$, $\hat{I}_{y}$ and $\hat{I}_{z}$ are the 8-dimensional operators describing the $x,y,z$ projections of the $I=7/2$ spin. In the presence of an RF electric field of amplitude $E_{1}$, the $\Delta m_{I}=\pm 1$ transitions are driven by an additional Hamiltonian term of the form:

[TABLE]

The $\Delta m_{I}=\pm 2$ transitions are driven by a term of the form:

[TABLE]

A detailed derivation of the matrix elements responsible for driving the $\Delta m_{I}=\pm 1$ and $\Delta m=_{I}\pm 2$ NER transitions is given in [14].

A finite element model is used to compute the strain and electric fields in the silicon layer near the donor position, using the COMSOL multiphysics software. The donor position has been triangulated by comparing simulated gate-to-donor coupling strengths with the experimentally observed strength, combined with the donor implantation profile (see Extended Data Fig.4 and [14]).

Kohn-Sham density functional theory (DFT) has been employed to calculate the S-tensor components that describe the impact of strain on the electric field gradient (EFG). To this end, 64- and 512-atom supercells were strained using the Projector Augmented-Wave (PAW) formalism [33] with a plane wave basis, as implemented in the Vienna Ab-Initio Simulation Package (VASP) [34, 35, 36]. The electric field response tensor is estimated by comparing the data points from the DC LQSE (Fig. 3C,D) and the $\Delta m_{I}=\pm 1$ (Fig. 2C) $\Delta m_{I}=\pm 2$ (Fig. 2D) Rabi frequencies with the simulated electric fields at the triangulated donor position. The final $R_{14}$ was found by minimizing the normalized residuals of the three separate $R_{14}$ estimates. Full theoretical modelling details can be found in [14].

Data availability: All data necessary to evaluate the claims of this paper is provided in the main manuscript or the supplementary information. Raw data files can be obtained from the corresponding author upon request.

Extended data

Supplementary Materials: Coherent electrical control of a single high spin nucleus in silicon

Contents

@starttoc

toc

S1 Antimony in silicon

Antimony (Sb) is a group-V atom which behaves as a donor when placed as a substitutional dopant in a silicon crystal. The neutral D0 state of the Sb donor has five valence electrons, four of which form covalent bonds with neighbouring silicon atoms, while the fifth remains loosely-bounded to the donor with a binding energy of 42.7 meV  [?]. This outer electron can be removed by tuning the electrochemical potential above the Fermi level of a nearby charge reservoir, resulting in an ionized D+ state of the donor.

The Sb isotope used in this experiment, ${}^{123}\mathrm{Sb}$, has a nuclear spin quantum number $I=7/2$, and therefore has eight ($2I+1$) nuclear eigenstates. All core electrons of the ${}^{123}\mathrm{Sb}$ donor, as well as the four electrons participating in the covalent bonds with the silicon lattice are paired in singlet states. Therefore, the spin Hamiltonian of the D0 state can only include the spin of the loosely-bound outer electron and the nuclear spin. In frequency units, the D0 spin Hamiltonian is

[TABLE]

where $\gamma_{e}=$27.97\text{,}\mathrm{GHz}\text{,}{\mathrm{T}}^{-1} is the electron gyromagnetic ratio, $\gamma_{n}=$-5.55\text{\,}\mathrm{MHz}\text{\,}{\mathrm{T}}^{-1} is the nuclear gyromagnetic ratio, $B_{0}$ is the static magnetic field, $A$ is the hyperfine interaction, equal to $101.52\text{\,}\mathrm{MHz}$ in bulk  [?] and assumed to be isotropic, $\vec{I}=\{\hat{I}_{x},\hat{I}_{y},\hat{I}_{z}\}$ ($\vec{S}=\{\hat{S}_{x},\hat{S}_{y},\hat{S}_{z}\}$) are the spin operators describing the $x,y,z$ projections of the nuclear (electron) spin, and $\hat{\mathcal{H}}_{Q}$ is the nuclear quadrupole Hamiltonian (Sec. S2).

In the limit where $\gamma_{e}B_{0}\gg A$, as is the case in our experiment ($B_{0}=$1.496\text{,}\mathrm{T}$$), the hyperfine interaction $A\vec{S}\cdot\vec{I}$ may be approximated as $A\hat{S}_{z}\hat{I}_{z}$. In this case, the nuclear and electron spin eigenstates are separable, and given by the tensor product of the eigenstates of the nuclear $\hat{I}_{z}$ and electron $\hat{S}_{z}$ spin operators. The hyperfine interaction between electron and nucleus introduces a dependence of the electron spin resonance (ESR) frequency on the nuclear spin state. This dependence is exploited to measure the nuclear spin state via the electron spin (Sec. S5).

If the D0 donor’s outer electron is removed, the resulting D+ donor has a considerably simpler Hamiltonian

[TABLE]

The nuclear Zeeman interaction separates the eight energy levels, resulting in seven identical nuclear transition frequencies (Fig. 1C). A nonzero quadrupole interaction separates these nuclear transition frequencies, allowing addressing of individual nuclear transitions (Fig. 1D). Throughout the remainder of the text, we will refer to the D+ donor Hamiltonian unless explicitly stated otherwise.

S2 Nuclear electric resonance and nuclear quadrupole interaction

A Nuclear quadrupole interaction

Assuming perfect spherical symmetry of the nuclear charge distribution, its electrical response to its environment is described by the Coulomb potential of a point charge. However, deviations from spherical symmetry require correction terms to this simple model. The lowest-order, non-zero correction term to the point charge model is the nuclear quadrupole interaction, which captures oblate or prolate deviations from spherical symmetry. This interaction is quadratic in the nuclear spin operators, and scales linearly with the EFG of the local electric field [?].

A general Hamiltonian $\hat{\mathcal{H}}_{Q}$ of the nuclear quadrupole interaction may be written as

[TABLE]

with $\hat{I}_{\alpha}$ as defined in S1, $Q_{\alpha\beta}$ with $\alpha,\beta\in\{x,y,z\}$, determines the strength of the different terms and is given by

[TABLE]

where $e$ is the elementary charge, $h$ is Planck’s constant, $I$ is the nuclear spin quantum number, and $q_{n}$ is the quadrupole moment, a constant specific for each nucleus, which parametrizes the degree of spherical asymmetry of the nuclear wavefunction. For 123Sb, $q_{n}=-0.69$ barn  [?,?], i.e. $-0.69\text{\times}{10}^{-28}\text{\,}\mathrm{m}\text{\,}{\mathrm{}}^{2}$. The EFG tensor components $\mathcal{V}_{\alpha\beta}=\partial^{2}V(x,y,z)/\partial\alpha\partial\beta$ are defined as the second partial derivatives of the electric potential $V(x,y,z)$ experienced by the nucleus. Note that energy is given in units of frequency throughout this work, hence the division in Eq. (4) by $h$. Furthermore, all Sternheimer anti-shielding factors  [?] are incorporated into $\mathcal{V}_{\alpha\beta}$. Finally, the EFG tensor is traceless and symmetric, as $V$ experienced at the nucleus obeys the Laplace equation $\nabla^{2}V=0$ and its partial derivatives commute. Therefore, $Q_{\alpha\beta}=Q_{\beta\alpha}$.

Two important aspects of $\hat{\mathcal{H}}_{Q}$ are worth emphasizing here. First, this interaction can only be non-zero for nuclei having $I\geq 1$. The quadrupole interaction was therefore irrelevant to the extensive body of earlier work on identical device structures containing a ${}^{31}\mathrm{P}$ donor with $I=1/2$. Second, as $q_{n}$ has a very small value, the EFG necessary to generate a significant quadrupole interaction needs to be large. Typically, only a microscopic mechanism in a crystal lattice or molecule, such as the distortion of covalent bonds in the vicinity of the nucleus, creates a significant EFG. An in-depth discussion of the microscopic origin of the EFG at the position of the nucleus in our experiment, and its dependence on electric field, may be found in Sec. S7. For the following, it suffices to know that such an EFG is present, and that it can be changed by varying an external electric field.

B Estimating nuclear quadrupole interaction parameters

As all experiments are performed at $B_{0}=$1.496\text{,}\mathrm{T}$$, the level splitting due to the Zeeman interaction is of order $8\text{\,}\mathrm{MHz}$. This is more than a factor 100 larger than the experimentally observed quadrupole shift, which is visible as a splitting of $66\text{\,}\mathrm{kHz}$ in the spectrum. This justifies treating the quadrupole interaction as a perturbation to the Zeeman interaction.

Eigenstates $\ket{m_{I}}$ of the unperturbed system are labelled with the secondary spin quantum number $m_{I}$, ranging from $-I$ to $I$ in steps of one. Applying standard first-order perturbation theory to find the quadrupole interaction correction to the energy, i.e.

[TABLE]

results in a transition frequency between neighboring levels given by

[TABLE]

where the transition frequency $f_{m_{I}-1\leftrightarrow m_{I}}$ corresponds to a transition between states $\ket{m_{I}-1}$ and $\ket{m_{I}}$.

The quadrupole splitting $f_{Q}$ introduced in the main text is now defined as the first-order splitting between two subsequent spectral lines

[TABLE]

Thus far, the quadrupole interaction Hamiltonian has been considered in the laboratory frame ($z$-axis along $B_{0}$, $x$-axis in plane with device substrate). However, since the EFG is a real, traceless, and symmetric tensor, a set of principal axes $x^{\prime}$, $y^{\prime}$, $z^{\prime}$ exists for which the tensor is diagonal. The corresponding diagonalized Hamiltonian $\hat{\mathcal{H}}_{Q}^{\prime}$ is given by

[TABLE]

with the pricipal axes chosen such that $|\mathcal{V}_{x^{\prime}x^{\prime}}|\leq|\mathcal{V}_{y^{\prime}y^{\prime}}|\leq|\mathcal{V}_{z^{\prime}z^{\prime}}|$. It is standard practice to define an asymmetry parameter $\eta=\frac{\mathcal{V}_{x^{\prime}x^{\prime}}-\mathcal{V}_{y^{\prime}y^{\prime}}}{\mathcal{V}_{z^{\prime}z^{\prime}}}$ to characterize the deviation from cylindrical symmetry ($0\leq\eta\leq 1$, $\eta=0$ corresponding to perfect cylindrical symmetry). In summary, the five independent degrees of freedom of the quadrupole Hamiltonian may now be understood as an overall interaction strength $Q$ (the prefactor to $I_{z^{\prime}}^{2}$ in Eq. (8)), an asymmetry $\eta$, and three angles defining the orientation of the principal axes.

We stress that the measured quadrupole splitting $f_{Q}$, as introduced in the main text and defined by Eq. (7), has no one-to-one correspondence to the ‘intrinsic’ quadrupole interaction strength $Q$. This is due to a strong dependence of the spectral line splitting on the value of the asymmetry parameter as well as the angle between external magnetic field and the quadrupole’s principal $z^{\prime}$-axis. Both of these quantities are unknown in the absence of any a-priori knowledge of the orientation and strength of the EFG. Given this, $f_{Q}$ may take on any value in the range $-2Q\leq f_{Q}\leq 2Q$. Extracting the five independent quadrupole degrees of freedom requires measuring the nuclear spectrum for different orientations of the external magnetic field, which is beyond the scope of the current work. Consequently, based on the available experimental data alone we cannot comment on the strength, asymmetry or orientation of the quadrupole interaction in the system. However, the observed splitting $f_{Q}$ in the spectrum gives a conservative estimate of a quadrupole interaction strength of at least $\frac{3eq_{n}\mathcal{V}_{z^{\prime}z^{\prime}}}{4I(2I-1)h}=$33\text{,}\mathrm{kHz}$$.

C Nuclear electric resonance

Resonant electric driving of transitions between nuclear spin states can be achieved by modulation of the quadrupole interaction. We postulate (see Sec. S7 for a thorough introduction and discussion) the existence of a modulation of the EFG components at the nucleus. This modulation has an amplitude $\delta\mathcal{V}_{\alpha\beta}$ (independent of frequency) and is of the form $\cos{\left(2\pi ft\right)}\delta\mathcal{V}_{\alpha\beta}$, where $f$ is the drive frequency. This EFG modulation causes a modulation of the quadrupole interaction (Eq. (4)), which we shall refer to as $\delta Q_{\alpha\beta}$. The time-dependent addition $\hat{\mathcal{H}}(t)$ to the static Hamiltonian $\hat{\mathcal{H}}$ (given by Eq. (2)) is

[TABLE]

Without prior knowledge of the exact EFG modulation at the nucleus, the absolute value of the Rabi frequencies of the different transitions cannot be predicted. Crucially, however, a clear and distinctive prediction about the relationship between the Rabi frequencies for the different nuclear electric resonance (NER) transitions is possible. For ${}^{123}\mathrm{Sb}$ with $I=7/2$, the spin matrices are given by

[TABLE]

Transitions between neighboring states $\ket{m_{I}-1}$ and $\ket{m_{I}}$ require the first off-diagonal matrix elements to be non-zero. Since $\hat{\mathcal{H}}(t)$ contains only products of the spin operators, the NER transitions between neighboring states $\ket{m_{I}-1}$ and $\ket{m_{I}}$ are only possible for terms involving the products $\hat{I}_{x}\hat{I}_{z}$, $\hat{I}_{z}\hat{I}_{x}$, $\hat{I}_{y}\hat{I}_{z}$, and $\hat{I}_{z}\hat{I}_{y}$. More formally, when $f$ equals the transition frequency for the transition between $\ket{m_{I}-1}$ and $\ket{m_{I}}$ ($\Delta m_{I}=\pm 1$), the resonant NER Rabi frequency $f^{\textrm{Rabi, NER}}_{m_{I}-1\leftrightarrow m_{I}}$ is given by

[TABLE]

The operator $\hat{I}_{x}\hat{I}_{z}+\hat{I}_{z}\hat{I}_{x}$ is given by

[TABLE]

and the purely imaginary operator $\hat{I}_{y}\hat{I}_{z}+\hat{I}_{z}\hat{I}_{y}$ has matrix elements with equal modulus. The relationship between the $\Delta m_{I}=\pm 1$ Rabi frequencies is determined by the relative magnitudes of the first off-diagonal matrix elements. In particular, the matrix element corresponding to the $\ket{-1/2}\leftrightarrow\ket{1/2}$ transition is zero and therefore this transition cannot be driven via NER. Taking the modulus of the matrix elements of Eq. (13) yields the Rabi frequencies $f^{\textrm{Rabi,NER}}_{m_{I}-1\leftrightarrow m_{I}}$ of the $\Delta m_{I}=1$ transitions, which are given by

[TABLE]

with $\alpha_{m_{I}-1\leftrightarrow m_{I}}$ given by

[TABLE]

Values of $\alpha_{m_{I}-1\leftrightarrow m_{I}}$ for the different transitions are given in Table S1.

Transitions between next-nearest-neighboring levels $\ket{m_{I}-2}\leftrightarrow\ket{m_{I}}$ ($\Delta m_{I}=\pm 2$) are allowed in first order for quadratic interactions such as the quadrupole interaction. Again, inspecting the spin matrices reveals that matrix elements coupling such states must stem from operator products $\hat{I}^{2}_{x}$, $\hat{I}^{2}_{y}$, $\hat{I}_{x}\hat{I}_{y}$, and $\hat{I}_{y}\hat{I}_{x}$, as only such products result in non-zero matrix elements on the second off-diagonal. Explicitly, the Rabi frequency for such transitions, $f^{\textrm{Rabi, NER}}_{m_{I}-2\leftrightarrow m_{I}}$, is given by

[TABLE]

The operator $\hat{I}^{2}_{x}$ is given by

[TABLE]

showing the coupling matrix elements on the second off-diagonal. The operator $\hat{I}^{2}_{y}$ has elements on the second off-diagonal with the same magnitude but opposite sign, whereas the operator $\hat{I}_{x}\hat{I}_{y}+\hat{I}_{y}\hat{I}_{x}$ is purely imaginary and its second off-diagonal elements have twice the strength of those of the other operators. This leads to a Rabi frequency given by

[TABLE]

with $\beta_{m_{I}-2\leftrightarrow m_{I}}$ given by

[TABLE]

Values of $\beta_{m_{I}-2\leftrightarrow m_{I}}$ for the different transitions are given in Table S1.

The NER Rabi frequencies for the $\Delta m_{I}=\pm 1$ and $\Delta m_{I}=\pm 2$ transitions, as given by Eq. (15) and Eq. (19) respectively, are only rescaled by the particulars of the EFG modulation. Importantly, the *relationship *between them presents a trend unique to electric drive via the nuclear quadrupole interaction.

These predictions for NER should be compared against a similar prediction for nuclear magnetic resonance (NMR). In the case of NMR, the driving of transitions is caused by an interaction with a time-varying magnetic field of the form $\cos{(2\pi ft)\gamma_{n}B_{1}\hat{I}_{x}}$. For transitions between neighboring levels, the NMR Rabi frequency is now determined by

[TABLE]

which leads to

[TABLE]

with $\zeta_{m_{I}-1\leftrightarrow m_{I}}$ given by

[TABLE]

In case an additional driving component along $\hat{I}_{y}$ is present, $B_{1}$ is replaced by the modulus of the total field, but $\zeta_{m_{I}-1\leftrightarrow m_{I}}$ remains unchanged. Values of $\zeta_{m_{I}-1\leftrightarrow m_{I}}$ for the different transitions are given in Table S1. Importantly, transitions between next nearest neighboring states ($\Delta m_{I}=\pm 2$) are to first order forbidden under NMR.

This experiment resulted in two key observations, as shown in the main text. First, the Rabi frequencies of the transitions between neighboring states accurately follow the trend predicted for NER in Eq. (15) (including the absence of the $\ket{-1/2}\leftrightarrow\ket{1/2}$ transition, which is predicted to have zero matrix element), and does not follow the prediction of Eq. (23) for NMR. Second, transitions between next-nearest-neighboring states, which are first-order forbidden for NMR, accurately follow the prediction of Eq. (19) for NER. These combined observations enable the main conclusion of this work: the demonstration of coherent electrical control of the ${}^{123}\mathrm{Sb}$ nuclear spin via NER, without transducing the RF electric field into a magnetic field.

S3 Device fabrication

A Fabrication protocol

Devices are fabricated using standard microelectronic fabrication techniques. The substrate is a p-type <100> silicon (Si) wafer (resistivity 10-$20\text{\,}\mathrm{\SIUnitSymbolOhm}\text{\,}\mathrm{cm}$), with a low pressure chemical vapor deposition (LPCVD) grown epitaxial layer of isotopically enriched ${}^{28}\mathrm{Si}$ (concentration of residual ${}^{29}\mathrm{Si}$ = 730 ppm) of $900\text{\,}\mathrm{nm}$ thickness on top.

The main fabrication steps are:

Definition of negative optical alignment markers via tetramethylammonium hydroxide (TMAH) etch of the silicon, using optical lithography and a wet thermally grown masking oxide. The masking oxide is subsequently removed using buffered hydrogen fluoride (HF) etching. 2. 2.

Creation of p-doped regions in between the n-doped metallic leads of step 4, defined using optical lithography and a wet thermally grown masking oxide. This step is designed to suppress spurious leakage currents flowing at the silicon / silicon dioxide (Si/SiO2) interface between the n-doped leads. Doping is achieved via thermal diffusion of boron (B). The masking oxide is subsequently removed using buffered HF etching. 3. 3.

Two-step thermal oxidation process to repair defects and drive B into the silicon, resulting in a final wet thermal oxide of $200\text{\,}\mathrm{nm}$ thickness. 4. 4.

Creation of n-doped metallic leads, defined by optical lithography using the thermally grown masking oxide of the previous step. Doping is achieved via thermal diffusion of phosphorus (P). Masking oxide is subsequently removed via buffered HF etching. 5. 5.

Single-step thermal oxidation to drive P into the silicon, resulting in a wet thermal oxide of $200\text{\,}\mathrm{nm}$ thickness which serves as thick field oxide. 6. 6.

Etching of a central window of 20\text{\,}\mathrm{\SIUnitSymbolMicro m}$\times$40\text{\,}\mathrm{\SIUnitSymbolMicro m} in the field oxide using buffered HF etching, and subsequent growth of an $8\text{\,}\mathrm{nm}$ thick, high quality dry thermal oxide. The p-doped regions and n-doped leads of steps 2 and 4 extend for $\sim$2\text{,}\mathrm{\SIUnitSymbolMicro m}$$ underneath the high quality gate oxide. 7. 7.

Definition of positive alignment markers via optical lithography. Markers are created by electron beam evaporation of $15\text{/}75\text{\,}\mathrm{nm}$ titanium/platinum (Ti/Pt) and lift-off using warm N-Methyl-2-pyrrolidone (NMP). 8. 8.

Definition of positive electron beam alignment markers by patterning a poly(methyl methacrylate) (PMMA) mask aligned to the optical alignment markers using electron beam lithography (EBL). Markers are created by electron beam evaporation of $15\text{/}75\text{\,}\mathrm{nm}$ titanium/platinum (Ti/Pt) and lift-off using warm NMP. 9. 9.

Definition of a 90\text{\,}\mathrm{nm}$$\times$$90\text{\,}\mathrm{nm} implantation window, patterned in a PMMA mask using EBL. 10. 10.

Low-dose implantation of ${}^{123}\mathrm{Sb}$ (see Sec. S3 B for details), followed by a rapid thermal anneal for $5\text{\,}\mathrm{s}$ at $1000\text{\,}\mathrm{\SIUnitSymbolCelsius}$ to achieve donor activation and repair implantation damage to the silicon lattice. 11. 11.

Definition of ohmic contacts to the n-doped metallic leads via optical lithography; ohmic contact is achieved by etching through the field oxide using buffered HF, evaporation of $200\text{\,}\mathrm{nm}$ aluminum (Al), and a $15\text{\,}\mathrm{min}$ forming gas (5% hydrogen in nitrogen) anneal at $400\text{\,}\mathrm{\SIUnitSymbolCelsius}$. 12. 12.

Definition of all gates necessary to control and read out the ${}^{123}\mathrm{Sb}$ donor, including gates to tune the electrochemical potential of the donor, a single electron transistor (SET) for read-out, and a microwave antenna. This is achieved via two steps of standard EBL using PMMA masks, thermal evaporation of Al and lift-off using warm NMP. The first layer consists of $20\text{\,}\mathrm{nm}$ of Al, and contains the SET barrier gates and two of the donor tuning gates. The second layer consists of $40\text{\,}\mathrm{nm}$ of Al, and contains the top-gate of the SET, the two remaining donor tuning gates, and the microwave antenna. The two layers are electrically isolated by the native oxide formed on the first Al layer upon exposure to air. 13. 13.

Forming gas anneal for $15\text{\,}\mathrm{min}$ at $400\text{\,}\mathrm{\SIUnitSymbolCelsius}$, to passivate traps in the gate oxide. 14. 14.

Dicing and packaging, including bonding (see also Sec. S4).

B ${}^{123}\mathrm{Sb}$ implantation parameters

Standard modeling using SRIM/TRIM software was used to predict the ${}^{123}\mathrm{Sb}$ dopant profile. The implantation energy was chosen such that the peak of the ${}^{123}\mathrm{Sb}$ dopant profile is located $2\text{\,}\mathrm{nm}$ below the SiO2/Si interface. Such a shallow implantation was chosen to maximize the donor’s exposure to static strain upon cool-down of the sample (see Sec. S7 A 1 for details on strain calculations). Given a gate oxide thickness of $8\text{\,}\mathrm{nm}$, the optimal implantation energy was found to be $8\text{\,}\mathrm{keV}$; the simulated vertical doping profile is shown in Fig. 1. A low dose of $2\text{\times}{10}^{11}\text{\,}{\mathrm{cm}}^{-2}$ was implanted. This corresponds to an average of 14 donors per implantation window, i.e. each device is expected to have a few donors tunnel-coupled to the single-electron transistor.

C Microwave antenna

The antenna that was originally designed to perform NMR and ESR is an on-chip coplanar waveguide coplanar waveguide (CPW) terminated by a short circuit. Due to the short circuit termination, the electric field at the tip of the antenna will be minimized, whilst a time-varying magnetic field is generated by the current flowing through the short.

In an effort to maximize the magnetic field $B_{1}$ at the donor site, multiple steps were taken to ensure the donor is located as close as possible to the antenna. This includes placing the implantation window no further than $150\text{\,}\mathrm{nm}$ away from the antenna and making the shorting wire as narrow as possible. The microwave antenna is made of aluminum (see step 12 of the fabrication procedure outlined in Sec. S3) and was designed to have a length of $500\text{\,}\mathrm{nm}$ on each side of the CPW, a thickness of $40\text{\,}\mathrm{nm}$ and a width of $50\text{\,}\mathrm{nm}$.

However, we have noticed across multiple devices that a short of these dimensions is not robust enough to withstand our standard experimental procedures such as connecting and disconnecting the microwave line, or applying high-power RF pulses (of the order of $-10$ to [math] dBm at the antenna). Upon inspection with a scanning electron microscope (SEM), the short typically appears ‘molten’ (Fig. 1A), which suggests that the cross-sectional area of the short is too small to sustain high-power microwaves or current spikes.

One can approximate the behaviour of a broken antenna, to first order, by modelling it as a series RC circuit. The electric field produced is proportional to the voltage at the capacitor, and the magnetic field is proportional to the current flowing through the circuit. Thus, the electric field will be at a relatively constant value for all frequencies until it starts rolling off at higher frequencies. The capacitor strongly attenuates any low-frequency current and thus causes the magnetic field to increase with frequency. At high frequencies, the impedance of the broken part of the antenna decreases to the point where a reasonable amount of current can flow. This all points to the antenna acting as an electric antenna at low frequencies ($\sim$8\text{,}\mathrm{MHz} needed for NER) and a (albeit poor) magnetic antenna at higher frequencies ($\sim$40\text{\,}\mathrm{GHz} needed for ESR).

To study the electromagnetic response in some depth, the width and location of the antenna breaks were extracted from the SEM (Fig. 1A) image and added to the original antenna design. We simulated the broken antenna using CST Studio Suite, and extracted the electric and magnetic fields as a function of frequency at the expected donor site (Sec. S7 A 2). Furthermore, the original, intact antenna was simulated in CST Studio Suite. This allowed the definition of a relative antenna performance, in which the electric and magnetic fields generated by the broken antenna at the donor site are divided by those of the working antenna.

The simulation results show a general trend that the magnetic field is reduced for a broken antenna and the electric field is enhanced, which match the experimental observations. When analyzing the simulation results, we find that the simulations for a broken antenna suggest that the magnetic field is on the order of $60\text{\,}\mathrm{dB}$ lower than that of a working antenna at ESR frequencies. A significantly reduced magnetic field at ESR frequencies has indeed been observed, but comparison to previous functional devices in our group suggests that the magnetic field is $\sim 30\text{\,}\mathrm{dB}$ weaker. The difference between experimentally observed values and simulation is likely due to the simplifications (exact shape of the discontinuity, grain structure, etc.) that were made in the CST model for ease of simulation.

The antenna simulations in the range of the nuclear resonance frequencies ($\sim 8\text{\,}\mathrm{MHz}$) show a highly-attenuated magnetic field ($\sim 100\text{\,}\mathrm{dB}$ attenuation), which would inhibit any coherent driving of the NMR transitions. Instead in this frequency range, the structure acts as a good electric antenna with $\sim$80\text{,}\mathrm{dB}$$ higher electric fields, enabling driving via NER. This explains how we were able to drive transitions via NER instead of NMR with a structure that was designed to be a broad-band magnetic antenna (Fig. 9). As outlined in the main text, the validity of this interpretation is confirmed by our observation of NER using a nearby electrostatic gate as an electrical antenna.

S4 Experimental setup

Device packaging consists of a custom-made printed circuit board (PCB) with microwave and low-frequency lines, positioned around a central cutout accommodating the fabricated chip. The PCB is mounted to a copper enclosure with SK, i.e. $2.92\text{\,}\mathrm{mm}$, (rated to $40\text{\,}\mathrm{GHz}$) and MMCX (rated to $6\text{\,}\mathrm{GHz}$) connectors. The device is glued to the enclosure and Al wire bonded to the PCB. The copper enclosure is mounted on a gold-plated copper cold-finger bolted to the mixing chamber of a Bluefors BF-LD400 cryogen-free dilution refrigerator at a base temperature of $12\text{\,}\mathrm{mK}$. The static magnetic field in the experiment is produced by a superconducting magnet kept in persistent mode. The magnet is equipped with a low-drift persistent mode switch that results in a typical magnetic field drift of less than 50 ppb/hour.

Three different types of control lines are present in the experimental apparatus. A single high-frequency microwave line is used for ESR in the $\sim$40\text{,}\mathrm{GHz}$$ regime, and has an inner/outer DC block at room temperature and a $10\text{\,}\mathrm{dB}$ attenuator at the $4\text{\,}\mathrm{K}$ stage of the dilution refrigerator. The coaxial microwave line has a silver-plated copper-nickel inner conductor, a copper-nickel outer conductor and a PTFE dielectric. Six radio-frequency coaxial lines are used for static and dynamic tuning of the donor electrochemical potential, driving via NER, and for the readout signal from the SET. These lines have a graphite coating on the dielectric to reduce triboelectric noise effects [?] and are low-pass filtered, with a $145\text{\,}\mathrm{MHz}$ cut-off. Three more lines of a Constantan loom, low-pass filtered to a $20\text{\,}\mathrm{Hz}$ cut-off, are used for the static electrical tuning of the SET.

The electronic setup, as depicted in Fig. 3, includes two instruments as DC sources, a Stanford Research Systems (SRS) SIM928 and National Instruments (NI) PXIe-4322. We use a total of eight arbitrary waveform generator (AWG) channels between the Keysight M3201A and M3300A modules, which are bandwidth limited to $200\text{\,}\mathrm{MHz}$, and a 100 MSPS digitizer channel of the M3300A is used to record SET current traces. These modules have on-board field-programmable gate arrays (FPGAs) where a variety of custom code is implemented, including an in-house direct digital synthesis (DDS) system allowing generation and sequencing of sinusoidal and linear chirp pulses. These are used as the IQ modulation inputs for the Keysight E8267D PSG vector microwave source, which is used to perform ESR.

To measure the nuclear flipping probability from NER, two microwave signals with ESR frequencies corresponding to the two nuclear states of interest need to be measured (Sec. S5). Single-sideband modulation, either upper or lower, is used to address the relevant ESR transitions. To avoid simultaneously driving both transitions with leakage of the alternate sideband, the microwave carrier frequency is offset from the center of the two transition frequencies. The SET current passes through the following chain of amplifiers and filters: FEMTO DLPCA-200 transimpedance amplifier, SRS SIM911 BJT amplifier and a SRS SIM965 low-pass filter. The Python-based QCoDeS data acquisition framework  [?] was used to control the instruments and perform measurements.

S5 Nuclear resonance measurement scheme and analysis

A Electron initialization, control, and readout

In the experiments outlined in this paper, the D0 donor’s outer electron is used as a tool to measure the nuclear spin state  [?]. In the presence of a static magnetic field $B_{0}$ along the $z$-axis, the Zeeman interaction causes an energy splitting $\gamma_{e}B_{0}$ between a spin-up $\ket{\uparrow}$ and spin-down $\ket{\downarrow}$ electron projected along the $z$-axis. If the electrochemical potential of a D0 donor is tuned such that a $\ket{\uparrow}$ electron is above the Fermi level, and a $\ket{\downarrow}$ electron below the Fermi level, the characteristic electron tunnel time from donor to SET is much shorter for a $\ket{\uparrow}$ electron ($\tau_{\ket{\uparrow}\rightarrow\mathrm{SET}}$) than a $\ket{\downarrow}$ electron ($\tau_{\ket{\downarrow}\rightarrow\mathrm{SET}}$). Note that $\tau_{\ket{\downarrow}\rightarrow\mathrm{SET}}$ is still finite, primarily due to thermal broadening of the states that are occupied in the SET. Therefore, if we keep these electrochemical potentials fixed for a duration $t$, where $\tau_{\ket{\uparrow}\rightarrow\mathrm{SET}}<t<\tau_{\ket{\downarrow}\rightarrow\mathrm{SET}}$, only a $\ket{\uparrow}$ electron will likely have tunneled onto the SET, an effect known as spin-dependent tunneling  [?]. Capacitive coupling between a donor electron and the SET enables the SET to be in Coulomb blockade exclusively when the donor is in the D0 charge state. Consequently, the tunneling of an electron from the donor onto the SET can be measured as an onset of SET current $I_{\mathrm{SET}}$ which, combined with spin-dependent tunneling, enables readout of the donor electron spin state.

Initialization of a $\ket{\downarrow}$ electron is performed using the electrochemical configuration described above, this time using the property that the electron tunnel time from SET to D+ donor is much shorter for a $\ket{\downarrow}$ electron. Readout and initialization of the donor’s outer electron can therefore be combined by remaining in the spin-dependent tunneling configuration for a fixed duration while measuring the SET current. Within this time, a $\ket{\uparrow}$ electron will likely tunnel onto the SET, followed by a $\ket{\downarrow}$ electron tunneling back onto the donor, which is measured as a brief increase in SET current. A $\ket{\downarrow}$ electron, on the other hand, will likely not tunnel out, and no SET current measured. In both cases, the final electron is in the $\ket{\downarrow}$ state. A more detailed description of readout and initialization of the donor electron is described in Ref.  [?,?].

Microwave-frequency linear chirp pulses are applied to adiabatically flip the electron spin via ESR. During the ESR pulse, the electrochemical potential of the donor is lowered well below the Fermi level, to ensure that the electron remains bound to the donor. The chirp pulse frequency range is chosen such that the ESR frequency is close to the center frequency, and the frequency ramp rate (chirpyness) is chosen to be below the adiabatic limit set by the Landau-Zener equation  [?], and above the limit set by the homogeneous linewidth. The main advantage of adiabatic pulses is that they give a high inversion fidelity without the need for a well-calibrated pulse amplitude, and are insensitive to slow resonance frequency drifts.

B Manipulation and readout of the nuclear state

The nuclear resonance measurements always focus on two chosen states out of the eight nuclear spin states of ${}^{123}\mathrm{Sb}$. We assume the nucleus to be initialized in one of these two relevant spin states. The general pulse sequence for measurements of the nucleus consists of two stages (Fig. 5). The first stage is nuclear spin manipulation, which we performed on the ionized donor (D+). A nearby AC gate applies oscillating electric fields that can drive Rabi oscillations between two nuclear spin states via NER.

After completing the nuclear spin manipulation, the final state is measured via the electron. To this end, a $\ket{\downarrow}$ electron is initialized, followed by an ESR chirp pulse that flips the electron only if the nucleus is in the first of the two states. The electron spin state is then read out by spin-dependent tunneling, where an electron tunneling event indicates that the electron was successfully flipped, and hence that the nucleus is in the nuclear state corresponding to the applied ESR frequency. This process of electron initialization, nuclear-state dependent flipping, and electron readout is then also performed for the ESR frequency corresponding to the other nuclear state. Since the electron readout method does not have a 100% readout fidelity, the two electron reads are repeated $N_{\mathrm{shots}}\approx 10$ times and the proportion of electron reads with tunneling event is counted for each nuclear state.

The nuclear state is successfully measured if exactly one of the two nuclear reads has a fraction of tunneling events above a preset threshold (Fig.  6A). A subsequent repetition of the pulse sequence will have a probability of flipping the nuclear spin to the other measured state, in which case the nuclear state with a high fraction of tunneling events will be switched. By repeating the pulse sequence $N_{\mathrm{repetitions}}$ times and counting the number of nuclear flips $N_{\mathrm{flips}}$, the resulting flip probability is $P_{\mathrm{flip}}=\dfrac{N_{\mathrm{flips}}}{N_{\mathrm{repetitions}}-1}$ (Fig.  6B). Occasionally, both electron spin-up fractions are either above or below the threshold. This could be caused by the nuclear state randomly flipping to one of the other six nuclear states, or simply due to the electron readout method not having a 100% fidelity. As a result, the nuclear state, and consequently $P_{\mathrm{flip}}$ could not be determined, and thus the data point is discarded. To account for effects such as slow drifts, multiple repetitions of the varied parameters were measured, and the final $P_{\mathrm{flip}}$ is taken as the average over the repeated measurements.

C Nuclear spin initialization

Throughout a measurement, the nuclear spin is operated in a two-dimensional subspace spanned by the two relevant spin states. However, two causes were found to result in undesired flipping of the nuclear spin out of this subspace. The first cause is due to the $\mathrm{D}^{0}$ and $\mathrm{D}^{+}$ donor having different spin eigenbases, and therefore a non-zero overlap between distinct nuclear spin eigenstates. The repeated ionization/neutralization that occurs during readout therefore yields a nonzero probability of the nuclear spin to transition to a different spin state. This effect has previously been observed in 31P (Ref.  [?]), but is enhanced for ${}^{123}\mathrm{Sb}$ due to the quadrupole interaction. The second cause is a flip-flop decay process when the donor is in the neutral state, where a $\ket{m_{I},\uparrow}$ state can decay into the $\ket{m_{I}+1,\downarrow}$ state. While this decay mechanism conserves the total angular momentum, the nuclear spin $z$-projection has increased. Consequently, the nucleus is likelier to settle at high $m_{I}$ values, which is in agreement with our experimental observations.

Two methods were employed to initialize the nuclear spin into a target state $\ket{m_{I}^{\mathrm{target}}}$. The first method is based on using calibrated NER $\pi$-pulses according to the following scheme

Determine the nuclear state by measuring the electron flipping probability for each of the eight ESR frequencies. Stop if the nucleus is in the state $\ket{m_{I}^{\mathrm{target}}}$. 2. 2.

Apply a sequence of NER $\pi$-pulses on an ionized donor to flip the nucleus to $\ket{m_{I}^{\mathrm{target}}}$. The pulse sequence consists of a minimal number of pulses, which can be a combination of $\Delta m_{I}=\pm 1$ and $\Delta m_{I}=\pm 2$ pulses. 3. 3.

Repeat from step 1.

While this method can be highly successful, it suffers from two significant drawbacks. The first drawback is that it is only possible to reach $\ket{m_{I}^{\mathrm{target}}}$ if all the $\pi$-pulses connecting the current and target state have been calibrated. The second drawback is that the nuclear resonance frequencies can slowly drift over time, and as a result the $\pi$-pulses need occasional recalibration.

Over the course of an experiment, a second nuclear spin initialization scheme was found and utilized that did not suffer from either of these drawbacks. The scheme relied on spin pumping using the flip-flop transitions, in which both the nuclear and electron spin states flipped (Fig. 7). Using high-power microwave chirp pulses (henceforth referred to as a flip-flop pulse), a transition between two flip-flop states could be addressed with a fidelity of $\sim 0.5$. For an initial $\ket{\downarrow}$ electron, application of a flip-flop pulse followed by a donor ionization results in a 50% chance of decreasing the nuclear spin ($\Delta m_{I}=-1$), being limited by decoherence. By repeatedly cycling through all flip-flop transition frequencies $f^{FF}_{\uparrow m_{I}-1\leftrightarrow\downarrow m_{I}}$ where the initial nuclear spin state $\ket{m_{I}}$ is above the target state ($m_{I}>m_{I,\mathrm{target}}$), the higher nuclear states are effectively pumped down to $\ket{m_{I}^{\mathrm{target}}}$ (Fig. 7C). A similar sequence is used for all nuclear spin states $\ket{m_{I}}$ below the target state $m_{I}<m_{I}^{\mathrm{target}}$, with the modification of starting with an initial $\ket{\uparrow}$ electron by applying two ESR pulses with frequencies $f^{\mathrm{ESR}}_{m_{I}}$ and $f^{\mathrm{ESR}}_{m_{I}+1}$. In this case, the flip-flop transition will increase the nuclear spin ($\Delta m_{I}=+1$), and cycling through the flip-flop transition frequencies effectively pumps all lower nuclear states up to $\ket{m_{I}^{\mathrm{target}}}$ (Fig. 7B).

One of the main advantages of this initialization scheme is that it does not require knowledge of the initial nuclear state, and instead cycles through the different flip-flop frequencies to pump the nuclear spin into a target state. This is an entirely open-loop control sequence to initialize the nuclear spin. The use of chirp pulses removes the requirement of an accurately-tuned flip-flop frequency, thereby negating the need to periodically recalibrate the flip-flop pulses. An additional benefit is that the flip-flop transition frequencies can be accurately estimated from the ESR frequencies alone, and thus do not require knowledge of nuclear transition frequencies.

S6 Ancillary measurements

A Slope in $\Delta m_{I}=\pm 2$ Rabi frequencies

Both the $\Delta m_{I}=\pm 1$ and $\Delta m_{I}=\pm 2$ NER Rabi frequencies should be symmetric about the center transition (Table. S1). While this symmetry has been largely observed (Fig. 2, C and D), the $\Delta m_{I}=\pm 2$ Rabi frequencies exhibit a small negative slope with increasing frequency. This slight asymmetry either has a physical cause not captured by the model Hamiltonian of the ionized nucleus, or is caused by an experimental artifact external to the donor, such as a frequency-dependent attenuation of the gate line in this frequency range.

To discriminate between these two possible causes, we reduced the static magnetic field $B_{0}$ from $1.496\text{\,}\mathrm{T}$ to $1.376\text{\,}\mathrm{T}$ (Decreasing transition frequencies by $20f_{Q}$), and remeasured the Rabi frequencies of the outer transitions $\ket{3/2}\leftrightarrow\ket{7/2}$ and $\ket{-7/2}\leftrightarrow\ket{-3/2}$ (Fig. 8A) keeping the driving amplitude constant. The two measured Rabi frequencies are both higher than their counterparts at higher magnetic field, thereby excluding that the effect is associated with specific nuclear transitions. Furthermore, by scaling the Rabi frequencies by their coefficients $\beta_{m_{I}-2\leftrightarrow m_{I}}$ (Table. S1), all Rabi frequencies follow a straight line with a gradient of $14.1(3)\text{\,}\mathrm{Hz}\text{\,}{\mathrm{MHz}}^{-1}$ (Fig. S6 AB). This is strongly suggestive of a dependence of the observed Rabi frequencies on the electric drive frequency, but not due to an underlying microscopic cause. The likely cause is a frequency dependence of the gate line attenuation, as room-temperature measurements of identical transmission lines have shown variations up to $1\text{\,}\mathrm{dB}\text{\,}{\mathrm{MHz}}^{-1}$ in the few-MHz regime. Although the electric drive amplitude is kept constant at the gate line input, this frequency-dependent attenuation causes the drive amplitude to vary at the donor site. This interpretation is further confirmed by the lack of a slope in the Rabi frequencies when driving with the antenna (Sec. S6 B), as the antenna transmission line (designed to operate up to $40\text{\,}\mathrm{GHz}$) does not have any significant frequency dependence in the few-MHz regime.

B Rabi oscillations driven via the microwave antenna

The NER spectrum has also been measured while driving with the damaged microwave antenna, from which the Rabi frequencies have been extracted for all $\Delta m_{I}=\pm 1$ (Fig. 9A) and $\Delta m_{I}=\pm 2$ (Fig. 9B) transitions. In both cases, the relationship between the Rabi frequencies match well with those expected for NER (Table S1). Additionally, the transition $\ket{-1/2}\leftrightarrow\ket{1/2}$ could not be driven, further confirming that the driving mechanism from the antenna in the $\sim$MHz frequency range is electric instead of magnetic.

The Rabi frequencies of the $\Delta m_{I}=\pm 2$ transitions are significantly higher than those of the $\Delta m_{I}=\pm 1$ transitions. This is opposite to the results obtained while driving with an electric gate, where the $\Delta m_{I}=\pm 1$ Rabi frequencies are much higher. This could be caused by the differing AC electric field orientation at the donor site when driving with an antenna versus with a donor gate, which will affect the ratio of $\Delta m_{I}=\pm 1$ to $\Delta m_{I}=\pm 2$ Rabi frequencies. Additionally, the $\Delta m_{I}=\pm 2$ antenna-driven Rabi frequencies do not exhibit a linear dependence on drive frequency, which is in contrast with the gate-driven Rabi frequencies (Sec. S6 A). This agrees with the estimate of a uniform transmission of electric signals (see Sec. S3 C in this frequency range, and further strengthens the above interpretation of non-uniform gate line transmission.

C Close-up of nuclear spectrum

The nuclear spectrum in Figs. 2A and 2B is presented over a wide frequency range, and does not highlight the shape of the spectral lines near the resonance frequencies. A close-up of each of the measured NER spectral lines is therefore shown in Fig. 10. Each resonance is power-broadened, and its shape is fitted to Rabi’s formula, given by

[TABLE]

where $\gamma=\pi f_{\mathrm{Rabi}}$ is the drive strength, $\Omega=\sqrt{\gamma^{2}+2\pi(f-f_{m_{I}-\Delta m_{I}\leftrightarrow m_{I}})/4}$, and $t$ is the NER pulse duration.

S7 Microscopic origin of the quadrupole interaction and nuclear electric resonance

The quadrupole interaction and NER in this device are related to the presence of a strain- and electric field-dependent EFG at the site of the ${}^{123}\mathrm{Sb}$ donor in silicon. In this section, we use finite-element models to estimate the strain and electric fields that are present in our device and to triangulate the position of the donor under study (Sec. S7 A 1). We then develop a microscopic model that relates the strain and electric fields to the EFG at the site of the donor, which directly relates to the quadrupole interaction strength. (Sec. S7 B).

By combining the strain profile from our finite element model with Density functional theory (DFT) calculations for the strain-EFG coupling, we estimate the portion of the quadrupole interaction that is due to strains arising from the thermal contraction of metallic gates in the device (Sec. S7 C). Then, combining the electric field profile from our finite element model with the measured NER shifts and Rabi frequencies, we estimate the portion of the quadrupole interaction that is due to external applied electric fields. We compare the inferred electric field EFG coupling to historical measurements of linear quadrupole stark effect (LQSE) in bulk materials (Sec. S7 D). Concluding this section, we use our finite element model to quantitatively rule out other conceivable physical mechanisms that could lead to NER (Sec. S7 E).

A Finite elements simulations of the device

1 The model set up

The finite element solver COMSOL was used for simulating electrostatics and strain due to the different thermal contraction of Si and Al upon cooling down the device. To achieve this, the electrostatics package (part of the AC/DC package), and the solid mechanics package (part of the structural mechanics package) were used, respectively. The model consists of a 2\text{\,}\mathrm{\SIUnitSymbolMicro m}$\times$2\text{\,}\mathrm{\SIUnitSymbolMicro m}$\times$2\text{\,}\mathrm{\SIUnitSymbolMicro m} silicon substrate that is mechanically fixed and electrically grounded at the bottom. An $8\text{\,}\mathrm{nm}$ thick silicon dioxide layer is defined on top of the silicon, and the aluminum gates are defined on top of this layer. Following the fabrication procedure described in Section S3, the aluminum is covered by $2\text{\,}\mathrm{nm}$ of aluminum oxide, and consists of 2 separate layers with a thickness of $20\text{\,}\mathrm{nm}$ for the bottom layer and $40\text{\,}\mathrm{nm}$ for the top layer. The gate layout is matched to that of the actual device upon its imaging and is shown in Fig. 11. In the electrostatic simulations, the gate-induced two-dimensional electron gas forming the SET leads and island is modelled by a 1 nm thick metallic layer at the SiO2/Si interface, with the lateral dimensions of the top gate TG, and gaps of the width of the barrier gates LB and RB. The two SET leads and island are kept at a potential of 0 V. A tetrahedral mesh was used, adjusting the properties to ensure elements do not have dimensions larger than $2\text{\,}\mathrm{nm}$ in the implantation region of the donor.

The structural deformation during the cooldown is modeled in two simulation steps. It is assumed that the device is strain free at $850\text{\,}\mathrm{\SIUnitSymbolCelsius}$, i.e. the growth temperature of the silicon oxide. At this stage, no aluminum gates are present. In the first step, the Si/SiO2 chip is cooled down from 850 to $400\text{\,}\mathrm{\SIUnitSymbolCelsius}$. Then, the Al gates are added to the structure and, at the temperature of $400\text{\,}\mathrm{\SIUnitSymbolCelsius}$, the Al/AlxOy gates are annealed in forming gas. We assume that the $400\text{\,}\mathrm{\SIUnitSymbolCelsius}$ forming gas anneal strain-relaxes the aluminum structures at this temperature. In the second step, including the initial strain of the silicon/silicon oxide, the deformation during the cooldown from $400\text{\,}\mathrm{\SIUnitSymbolCelsius}$ to $200\text{\,}\mathrm{mK}$ is modelled. At $200\text{\,}\mathrm{mK}$ the thermal deformation has saturated and further cooldown has a negliable effect.

Silicon is an anisotropic material. When setting up the model, special care needs to be taken in correctly aligning the orientation of the aluminum gate layout with the silicon crystal axes. Our sample is fabricated on top of a (100) wafer, where [001] is the crystal direction normal to the wafer surface. The wafer was diced along the in-plane crystal directions [110] and [$\bar{1}$10] (the natural cleaving directions) and the device subsequently aligned to these directions. Consequently, the [110] and [$\bar{1}$10] directions correspond to the $z$ and $x$ axis in the COMSOL model, respectively, and this is the coordinate convention used throughout the manuscript. The stiffness matrix for silicon is commonly defined with respect to the [100], [010] and [001] crystal axes. We define a new coordinate system with respect to the crystal axis $x^{\prime\prime}$ along [010], $y^{\prime\prime}$ along [001], and $z^{\prime\prime}$ along [100]. To match our chosen lab frame ($x$, $y$, $z$) with the standard crystal orientation ($x^{\prime\prime}$, $y^{\prime\prime}$, $z^{\prime\prime}$), either the model has to be rotated by $\pi$/4 about $y$, or the stiffness tensor has to be rotated. We chose the first option, as it is more convenient to rotate electric field components and strain matrices than the stiffness tensor. In what follows, all spatial maps of quantities will be shown in the lab frame.

The aluminum oxide is assumed to be isotropic. All thermal expansion coefficients used in the strain simulations can be found in Table S2. All other material properties are taken from the COMSOL material library.

2 Donor triangulation

The donor’s electrochemical potential, which sets the ionized-to-neutral charge transition, is influenced differently by different gates, according to their capacitive coupling to the donor charge. The allows us to triangulate the location of the donor by comparing the experimental ratios of different gate capacitances with electrostatic COMSOL simulations.

Fig. 12 shows the SET current as a function of the voltages on PL and one of the donor gates, with a different donor gate in each panel. A discrete shift in the SET Coulomb peaks indicates a charge transfer, e.g. the ionization of the measured donor marked by the white dashed lines. For each gate voltage configuration along these lines, the electrostatic potential at the donor position $V(\vec{r}_{0};V_{\rm PL},V_{\rm DFL},...)$ is constant. Changing the gate voltages in Fig. 12A by $\Delta V_{\rm PL}$ and $\Delta V_{\rm DFL}$ along the dashed line, we may to first order write

[TABLE]

We identify the left hand side as the experimentally measured slope

[TABLE]

which we can equivalently define for the three other panels. To estimate the precision of extracting these slopes, we have introduced the angle $\theta$. With an uncertainty $\sigma_{\theta}=$1\text{,}\mathrm{\SIUnitSymbolDegree}$$, Gaussian error propagation leads to an estimated standard deviation $\sigma_{\rm DFL}$ given by

[TABLE]

We now triangulate the donor by finding the position $\vec{r}$ at which an electrostatic simulation of the right hand side in Eq. (27) best matches the experimentally measured slopes. We use COMSOL to compute the electrostatic potential landscape across the model for the set of gate voltages given in Table S2. Subsequently, the relevant gate voltages are varied by $10\text{\,}\mathrm{mV}$, and the electrostatic potential landscape across the model is computed again for variation of each gate voltage. This results in a spatially-varying simulated slope, which, at a given position $\vec{r}$, is defined as

[TABLE]

Here again the DFL gate is used as an example to illustrate the procedure, and the plunger gate is used as the common gate.

To compare the simulated slopes to the measured ones at each point within the model space, we use a least-squares estimate [?] that is defined by

[TABLE]

with $P_{\mathrm{triangulation}}(\vec{r})$ the probability density, which we normalize over the given volume via $\mathcal{N}$, and the summation runs over the DFL, DFR, DBL and DBR gates.

We have verified that the thickness and exact geometry of the two-dimensional electron gas has negligble influence on the triangulation probability density function. Furthermore, we have implemented a Thomas-Fermi approximation to the density of the two-dimensional electron gas in one specific geometry and found that this only negligibly alters the triangulation probability density function, and we therefore choose the less computationally expensive, metallic approximation.

The resulting triangulation probability density function is shown in Fig. 13(A, C). The planar layout of the gate structures leads to a low sensitivity of the capacitance triangulation method in the out-of-plane (vertical) $y$ direction, as reflected by a large uncertainty along this axis, visible in Fig. 13(C). However, a significant section of this region constitutes highly unlikely positions of the donor, as revealed by the donor implantation depth profile in Fig. 1. A more accurate probability distribution of the donor’s position is therefore given by

[TABLE]

with $P_{\mathrm{implantation}}(\vec{r})$ the dopant probability distribution and $\mathcal{N}$ a normalization factor. $P_{\mathrm{implantation}}(\vec{r})$ is the distribution shown in Fig. 1, a modified Gaussian distribution homogenous in the lateral dimensions $x$ and $z$ and with mean at -2.5 nm, standard deviation of 2.6 nm and slight skewness (0.54) and kurtosis (3.59). This distribution is cut off at the Si/SiO2 interface, and we have verified that donor diffusion due to the rapid thermal anneal of $5\text{\,}\mathrm{s}$ at $1000\text{\,}\mathrm{\SIUnitSymbolCelsius}$ has a negligible effect.

Further methods to narrow $P(\vec{r})$, e.g. by modeling the donor-SET tunnel coupling [?], were not adopted here.

Fig. 13(B, D) shows color maps of this final probability density function for the donor’s location. The best estimate (maximum of $P(\vec{r})$) is found at a lateral position of $(x,z)=($13\text{,}\mathrm{nm}$,$8\text{,}\mathrm{nm}$)$ and a depth of $-5\text{\,}\mathrm{nm}$ below the Si/SiO2 interface.

B Microscopic origin of the electric field gradient

As discussed in Section S2, the quadrupole interaction is generated by an EFG at the site of the ${}^{123}\mathrm{Sb}$ nucleus. This EFG is due to inhomogeneities in both the external applied electric field and the internal electric field arising from the charged electrons and nuclei comprising the host crystal. Due to the relatively large physical length scales of the gates that supply the external field, the former contribution can be ignored, as quantitatively supported by a finite element model in Sec. S7 E 2. Here, we focus entirely on a microscopic model for the distortion of the electronic charge distribution around the ${}^{123}\mathrm{Sb}$ nucleus. We propose that both strain and externally applied electric fields distort the chemical environment coordinating the ${}^{123}\mathrm{Sb}$ nucleus, and thus give rise to a non-zero EFG, providing the physical mechanism that allows us to drive NER with an AC electric field. The coupling of the externally applied electric field to the EFG also leads to a manifestation of the LQSE, in which the quadrupole splitting can be tuned with a DC electric field (see main text for details).

We first consider a generic expression that relates the total charge density of the electrons and nuclei, $\rho(\vec{r})$, to the EFG tensor, $\mathcal{V}_{\alpha\beta}$. Taking the ${}^{123}\mathrm{Sb}$ nucleus as the origin of coordinates,

[TABLE]

where $\vec{r}=(r_{x},r_{y},r_{z})$ is the position vector, $\Delta_{\alpha\beta}$ is the Kronecker delta, and $\alpha,\beta\in\{x,y,z\}$. The integration kernel coefficient to $\rho(\vec{r})$ will remain invariant under the application of strain or an electric field, the effects of which will enter through their coupling to the charge density. Expanding the charge density as a first-order functional Taylor series about its form at zero strain and zero electric field we arrive at

[TABLE]

Here $\epsilon_{\gamma\delta}$ is the strain tensor and $E_{\gamma}$ is the external applied electric field, where $\gamma,\delta\in\{x,y,z\}$ index their respective components. This suggests the simplified functional form,

[TABLE]

$\mathcal{V}^{0}_{\alpha\beta}$ is the EFG arising due to the unperturbed charge density. $S_{\alpha\beta\gamma\delta}$ is the fourth-rank gradient-elastic tensor that describes the effect of strain on the EFG at the ${}^{123}\mathrm{Sb}$ nucleus. $R_{\alpha\beta\gamma}$ is the third-rank electric field response tensor, capturing the relationship between the electric field and the EFG at the ${}^{123}\mathrm{Sb}$ nucleus. The symmetry properties of the tensors that appear in Eq. (35) suggest a more compact representation using Voigt’s notation. $\mathcal{V}_{\alpha\beta}$ is a symmetric tensor with 6 independent components that we represent with integer indices, i.e., $\mathcal{V}_{1}=\mathcal{V}_{xx},\mathcal{V}_{2}=\mathcal{V}_{yy},\mathcal{V}_{3}=\mathcal{V}_{zz},\mathcal{V}_{4}=\mathcal{V}_{yz}=\mathcal{V}_{zy},\mathcal{V}_{5}=\mathcal{V}_{xz}=\mathcal{V}_{zx},$ and $\mathcal{V}_{6}=\mathcal{V}_{xy}=\mathcal{V}_{yx}$. Similarly, because $S$ and $R$ are symmetric they can be represented by 36 and 18 independent components, respectively, rather than 81 and 27. The tetrahedral Td symmetry of the ${}^{123}\mathrm{Sb}$ donor site permits an even more dramatic reduction in the number of independent components of $S$ and $R$ that will be discussed in Sec. S7 C 1 and Sec. S7 D 1.

We only consider the three leading order contributions to the EFG in detail. It is worth noting that the higher-order mixed terms may be of interest in future experiments. Specifically, shear strains may break symmetries that lead to a permanent dipole moment that can enhance the electric field’s coupling to the EFG, relative to uniaxial strains. However, this would require measurements of samples under different or variable strains to be substantiated.

We first consider whether the strain and electric field independent contribution, $\mathcal{V}^{0}_{\alpha\beta}$, vanishes by symmetry. This can be determined by identifying whether the representation of the integrand, a product of $\rho(\vec{r})$ and the EFG kernel, includes the totally symmetric irreducible representation of the point group associated with the symmetry of the substitutional site. In the absence of strain or an external applied electric field, this point group is Td and it is evident that the totally symmetric irreducible representation does not appear, leading this contribution to vanish by symmetry. Our measurements indicate a non-zero EFG, so the symmetry of the ${}^{123}\mathrm{Sb}$ site must be effectively lowered by strain and/or an external applied electric field.

Both uniaxial and shear strains will change the positions of the silicon atoms that coordinate the ${}^{123}\mathrm{Sb}$ donor such that the Td symmetry of the substitutional site is broken. Uniaxial strain will reduce the site’s symmetry to D2d, giving rise to a linear coupling between the strain and the diagonal components of the EFG tensor. Shear strain will reduce the site’s symmetry to C2v, giving rise to a linear coupling between the strain and the off-diagonal components of the EFG tensor. These couplings correspond to the two non-trivial components of the $S$ tensor for Td symmetry, $S_{11}$ and $S_{44}$, respectively. DFT calculations of these couplings are given in Sec. S7 C 2. These calculations include the direct evaluation of the EFG as a function of strain and corroborate a linear relationship that extrapolates to zero in the absence of strain.

Even in the absence of strain, an external applied electric field can create a non-vanishing EFG. Because the silicon crystal is not piezoelectric, the effect of this electric field will not be to change the positions of the coordinating silicon atoms. We substantiate in Sec. S7 E 1 that an indirect modulation of the strain through a piezoelectric coupling in the oxide is not a significant effect. Rather, the primary effect of an electric field will be to mix the electronic ground state of the ionized ${}^{123}\mathrm{Sb}$ with excited states of opposite parity, proportional to the transition dipole moment between them. This polarization leads to a symmetry breaking in $\rho(\vec{r})$ that gives rise to a linear coupling between the electric field and the EFG through the $R$ tensor. From symmetry considerations, the only non-trivial component will be $R_{14}$, which gives rise to off-diagonal contributions to the EFG tensor. This direct coupling between the electric field and the EFG is the physical mechanism that allows us to drive NER with an external applied AC electric. That this coupling is present and of the correct order of magnitude to rationalize our observations is further substantiated by the observed shift of the nuclear spectral lines in response to a DC electric field, which results from the same mechanism. In the NER literature, the resulting shift in transition frequencies is known as the LQSE. The presence of the LQSE relies on a product of matrix elements with opposite parities under inversion, and thus it is forbidden by symmetry for nuclei that inhabit sites with inversion symmetry. Because the ${}^{123}\mathrm{Sb}$ donor does not occupy such a site, it is plausible in principle that such an electric field induced shift can explain our observations. We provide an estimate of this coupling from our measurements of the DC electric field induced shifts in Sec. S7 D 2 and a comparison to measurements of the LQSE in bulk crystals in Sec. S7 D 3.

C Strain induced quadrupole splitting

1 The gradient-elastic tensor

For Td tetrahedral symmetry the gradient-elastic tensor $S$ has two different nonzero components $S_{11}$ and $S_{44}$. Here, Voigt’s notation in the conventional crystallographic, non-rotated basis, and engineering strains, are used. In this frame (given by $x^{\prime\prime}$, $y^{\prime\prime}$, $z^{\prime\prime}$, see Sec. S7 A 1), the EFG tensor is given by

[TABLE]

Note that the relationship between shear components of the engineering and infinitesimal strain is $\gamma_{y^{\prime\prime}z^{\prime\prime}}=2\epsilon_{y^{\prime\prime}z^{\prime\prime}}$. DFT calculations have been performed (see Sec. S7 C 2) to evaluate $S_{11}$ and $S_{44}$. In the experiment, the magnetic field is applied along the [110] crystal direction, e.g. in the frame of the gradient-elastic tensor above this corresponds to $\hat{\mathcal{H}}=\gamma_{n}B_{0}(\hat{I}_{x^{\prime\prime}}+\hat{I}_{z^{\prime\prime}})/\sqrt{2}$. Computing the first order perturbative correction, as in Sec. S2 B, yields the quadrupole splitting

[TABLE]

Rotating the strain tensor by $\pi/4$ about $y$ to align the magnetic field with the $z$ axis results in the laboratory frame expression for the spectral quadrupole splitting given by

[TABLE]

2 Gradient-elastic tensor calculations

We use Kohn-Sham density functional theory (DFT) to develop an atomistic understanding of the impact that strain has on the EFG at the ionized ${}^{123}\mathrm{Sb}$ donor. We perform a series of supercell calculations to determine the manner in which the silicon atoms coordinating the ${}^{123}\mathrm{Sb}$ site relax as a function of strain. These calculations also provide us with first principles values for the $S$ tensor that avoid the need for empirically derived Sternheimer factors.

For all supercell calculations we use the Projector Augmented-Wave (PAW) formalism [?] with a plane wave basis, as implemented in the Vienna Ab-Initio Simulation Package (VASP)  [?,?,?]. Within the PAW formalism the all-electron Kohn-Sham orbitals, and their associated density, can be accessed via an explicit linear transformation on smoother pseudo orbitals that can be efficiently represented in a plane wave basis. As such, we can accurately compute properties that depend sensitively on the charge density and local potential near nuclei, while also making use of the frozen core approximation to reduce the number of orbitals that must be explicitly included in a given calculation. One such property is the EFG, which has been shown to be accurately represented within the PAW formalism when compared to both experiments and calculations using the linear augmented plane wave (LAPW) formalism  [?].

Because we are studying an ionized Sb donor, we do not need to represent its hydrogenic bound state in our calculations and we can thus work with smaller supercells than are required for neutral shallow defects  [?]. That we do not need to represent states near the conduction band edge also mitigates concerns about the impact of the band gap problem on our results. We thus find the semilocal Perdew-Burke-Ernzerhof (PBE) construction of the generalized gradient approximation to the exchange-correlation functional  [?] to be adequate. Our calculations are carried out using a plane wave cut-off of $500\text{\,}\mathrm{eV}$ for the orbitals and $1000\text{\,}\mathrm{eV}$ for the augmentation charges, with a 3$\times$3$\times$3 Gamma-centered sampling of the first Brillouin zone for 64 atom supercells, and a Gamma-only sampling for 512 atom supercells. Even though we do not need to represent the neutral donor wave function, finite size effects still impact our results because our supercell calculations are non-cubic and have a net positive charge. It is well-known that total energy and forces in such calculations are slow to converge with system size. Dipole corrections that would otherwise account for finite size errors are only implemented for cubic supercells in the software package that we are using. The values of the EFG computed for the 64 and 512 atom supercells at a given strain are within 10% of one another and that this level of accuracy suffices to demonstrate a linear trend across different strains. In the worst case, the difference between the values of the slope obtained regressing on a linear model with a zero $y$-intercept and a non-zero $y$-intercept is 0.1%.

We self-consistently relax the Si crystal geometry around the ionized ${}^{123}\mathrm{Sb}$ donor for both uniaxial strain on [100] and shear strain at 5 different values (-10*-2*, -3$\times 10^{-3}$, -10*-3*, -3$\times 10^{-4}$, and -10*-4*). Our starting point for each calculation is a pristine Si crystal with the strain applied relative to the experimental lattice constant ($5.431\text{\,}\mathrm{\SIUnitSymbolAngstrom}$) and a single substitutional ${}^{123}\mathrm{Sb}$ donor. To account for the +1 ionization state of the Sb donor, we include 4 electrons per atom in our supercells, and a homogeneous background charge is assumed to maintain the overall neutrality necessary to obtain a finite total energy. We relax the atomic positions at a fixed supercell volume until the interatomic forces are all less than $1\text{\,}\mathrm{meV}\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{-1}$, and find that in all cases the primary impact of the ionized ${}^{123}\mathrm{Sb}$ is to displace the 4 silicon atoms that coordinate it by $\sim$ $0.2\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. The next nearest neighbors are displaced by $\sim$ $0.05\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, and beyond that all displacements are below $\sim$ $0.025\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. This is consistent with a simple picture in which the ionized Sb occupies a larger volume than neutral Si by merit of its filled M and N shells. For each strain type and amplitude we compute the EFG and perform a linear regression on the computed values to extract estimates for $S_{11}$ and $S_{44}$.

For the case of uniaxial strain, we evaluate $S_{11}$ for both a model in which the Poisson ratio is 0 (i.e., unphysical “cork-like” Si) and another in which it takes on its experimental bulk value of 0.28. The resulting values of $S_{11}$ are $1.9\text{\times}{10}^{22}\text{\,}\mathrm{V}\text{\,}{\mathrm{m}}^{-2}$ and $2.4\text{\times}{10}^{22}\text{\,}\mathrm{V}\text{\,}{\mathrm{m}}^{-2}$, respectively, noting that this distinction preserves the order of magnitude of the coupling. For shear strain, we compute $S_{44}$ to be $6.1\text{\times}{10}^{22}\text{\,}\mathrm{V}\text{\,}{\mathrm{m}}^{-2}$. These values are comparable to those reported for As in Si by Franke et al., which are respectively $1.5\text{\times}{10}^{22}\text{\,}\mathrm{V}\text{\,}{\mathrm{m}}^{-2}$ and $6.0\text{\times}{10}^{22}\text{\,}\mathrm{V}\text{\,}{\mathrm{m}}^{-2}$  [?]. That the ratio of the shear component to the uniaxial component is larger for As than Sb is to be expected because Sb occupies a larger volume than As and will thus have an EFG that couples more strongly to strains that change volume.

Our calculations were carried out using 4 electron PAW potentials for Si and 15 electron PAW potentials for Sb. To assess the primary chemical contribution to the EFG, we repeated the calculations with a 5 electron PAW potential for Sb and found that the EFG decreased by 5$\%$ for both models of uniaxial strain and 2.5$\%$ for shear strain. This small change confirms that the EFG is generated primarily by the contribution to the local potential due to the charge density in the sp3-like Sb-Si bonding orbitals, and not a strain-induced distortion of the on-site contribution to the local potential due to the Sb donor’s d electrons.

The relaxed supercell geometries expose microscopic details of the symmetry breaking at the ${}^{123}\mathrm{Sb}$ site. In the absence of strain this site has the symmetry of the point group Td, which is apparent from the tetrahedral coordination of the four nearest silicon atoms. Uniaxial strain reduces the symmetry of this site to D2d, whereas shear strain reduces it to C2v. We note that this detail may become important when considering contributions to the EFG beyond first order to which both the strain and electric field are coupled, and will be discussed briefly in Sec. S7 D 3.

3 Calculation of quadrupole splitting due to strain

As $S_{44}/S_{11}\approx 2.5$, Eq (38) shows that the leading contribution to the quadrupole splitting will be the shear strain along the magnetic field aligned with the [110] axis. We use COMSOL to compute the thermal deformation and show the shear strain $\epsilon_{xx}-\epsilon_{zz}$ in the main text (Fig. 1B). We use Eq. (36) to relate the computed strains to EFG and numerically compute the eigenenergies of Eq. 2 and show the estimated quadrupole splitting $f_{Q}$ in the main text (Fig. 4C). The variation in quadrupole splitting shows a good qualitative correspondence with the variation of the shear strain, reflecting the larger contribution of $S_{44}$ to the quadrupole splitting. Values found for the strain-induced quadrupole splitting range between 90 and $160\text{\,}\mathrm{kHz}$ throughout the area where the donor is estimated to be. This result represents an outstanding quantitative agreement with the experimentally observed value of $f_{Q}=$66\text{,}\mathrm{kHz}$$, in light of the complex combination of strain modeling, DFT calculations and donor triangulation required in our modeling effort.

D Electric field induced quadrupole splitting and NER

1 The electric field response tensor

The effect of an applied electric field on the EFG can be described by a third rank tensor. In Td symmetry we only have one unique non-vanishing element $R_{14}=R_{25}=R_{36}$. Here, the first index refers to the electric field component at the donor site in the crystallographic frame (i.e. $x^{\prime\prime}=1,y^{\prime\prime}=2,z^{\prime\prime}=3$) and the second index refers to the resulting EFG component in Voigt’s notation. In the frame given by the crystal axes ($x^{\prime\prime}$, $y^{\prime\prime}$ and $z^{\prime\prime}$) the tensor acts as follows:

[TABLE]

We note that the above assumption of Td symmetry is not entirely obvious, as strain breaks this symmetry. However, to first order, the effects of electric field and strain are independent; the development of a more elaborate theory will be the topic of future research.

The coupling between an external applied electric field and the EFG at the ${}^{123}\mathrm{Sb}$ nucleus is the physical origin of both NER and the LQSE. Both AC and DC fields will act to modulate the off-diagonal components of the EFG tensor relative to the crystallographic axes of the silicon host. It is important to note that this is not the same coordinate system as the laboratory frame, which is defined relative to the $B_{0}$ field. Thus, off-diagonal modulation in the crystal frame can lead to modulation of both the off-diagonal (NER) and diagonal components (LQSE) in the laboratory frame. Future experiments could employ a vector magnet to explore these details. For example, one could imagine aligning the external applied magnetic field such that the laboratory and crystal coordinate frames coincide to suppress LQSE.

2 Electric field response tensor estimate

As for the strain in Sec. S7 C 1, we can derive an analytic expression for the shift of the quadrupole splitting with electric field, assuming a strong magnetic field applied along the [110] crystal direction. We rotate the lab frame expression for the quadrupole splitting (Eq. (7)) into the crystal frame

[TABLE]

As discussed and shown above (Eq. (39)), an electric field only gives rise to off-diagonal EFG tensor components in the crystal frame. Thus, only $\Delta\mathcal{V}_{x^{\prime\prime}z^{\prime\prime}}=R_{14}\,\Delta E_{y^{\prime\prime}}$ contributes to the spectral line shift as a function of electric field,

[TABLE]

and since the transformation between crystal and lab frames is a rotation about $y$, the same analytic form is found in the lab frame:

[TABLE]

Similarly, we can derive analytic expressions for the expected Rabi frequencies of the $\Delta m_{I}=\pm 1$ and $\Delta m_{I}=\pm 2$ transitions under NER. Applying a sinusoidal pulse to an electric gate with amplitude $\delta V$ will cause an electric field modulation at the donor position, $\delta\vec{E}$. This in turn gives rise to the modulation of the EFG, $\delta\mathcal{V}_{\alpha\beta}$, postulated in Sec. S2 C. There, we derived the lab frame representation of the NER $\Delta m_{I}=\pm 1$ transition Rabi frequencies in Eq. (15), which, upon rotation into the crystal frame, becomes

[TABLE]

Again, only off-diagonal EFG tensor components contribute in the crystal frame, and we find

[TABLE]

and for the lab frame

[TABLE]

Note that only the electric field component along the applied magnetic field contributes. Repeating the procedure for the $\Delta m_{I}=\pm 2$ transition Rabi frequencies using Eq. (19), we find

[TABLE]

for the crystal frame, and for the lab frame

[TABLE]

We can exploit these linear relationships to estimate the value of $R_{14}$. We use COMSOL to model the out-of-plane DC electric field shift $\Delta E_{y}$ for increasing the gate voltage $V_{\rm DFL}$ by $20\text{\,}\mathrm{mV}$, as shown across the device in Fig. 14A. This is the relevant electric field component for the LQSE as a function of gate voltage, as shown in Eq. (42). Furthermore, the electric drive is modelled by applying the peak voltage of a sinusoidal pulse to the DFR gate, $\delta V_{\rm DFR}=$20\text{,}\mathrm{mV}$$, as used in the experiment. The relevant electric field components determining the Rabi frequency as given by Eq. (46) and Eq. (50) are shown across the device in Fig. 14C for the $\Delta m_{I}=\pm 1$ transitions and in Fig. 14E for the $\Delta m_{I}=\pm 2$ transitions.

For the DC spectral line shift, the experiment found that $\Delta f_{Q}=$9.9\text{,}\mathrm{kHz} per applied volt on gate DFL (main text). For the electric drive, the experiment found that $f^{\mathrm{Rabi,NER}}_{5/2\leftrightarrow 7/2}=$684\text{\,}\mathrm{Hz} and $f^{\mathrm{Rabi,NER}}_{3/2\leftrightarrow 7/2}=$38.5\text{,}\mathrm{Hz} for $20\text{\,}\mathrm{mV}$ drive amplitude. Using these three experimental values and the simulated electric fields throughout the device, $R_{14}$ is the only unknown constant. Fig. [14](#F14)(B, D, F) shows the $R_{14}$ values needed to explain the measured DC spectral line shift (B), the Rabi frequencies of the $\Delta m_{I}=\pm 1$ NER transitions (D) and the Rabi frequencies of the $\Delta m_{I}=\pm 2$ NER transitions (F). At the most likely donor position (black marker), as determined by combining triangulation with the ion implantation profile, we find values of $R_{14}=$10\text{\times}{10}^{12}\text{\,}{\mathrm{m}}^{-1} (B), $R_{14}=$1.7\text{\times}{10}^{12}\text{,}{\mathrm{m}}^{-1} (D) and $R_{14}=$0.7\text{\times}{10}^{12}\text{\,}{\mathrm{m}}^{-1} (F). The order of magnitude spread of these values should be viewed as a remarkable achievement, given the complexity of the models and the absence of any free parameters in them.

We combine these results and determine a unique value of $R_{14}$ at the triangulated donor position by minimizing the sum of the square of the normalized residuals

[TABLE]

Here, $M_{i}$ is the observed measurement value, i.e. the left hand side in Eq. (42), Eq. (46) or Eq. (50). $S_{i}(\vec{r})$ is the product of all known constants on the right hand side of Eq. (42), Eq. (46) or Eq. (50) and the electric field component. $R_{14}(\vec{r})$ is the targeted minimization value. $\sigma_{M_{i}}$ is the uncertainty of the measurement value $M_{i}$. Assigning an equal measurement error to each observation, we find at the most likely donor position $R_{14}=$1.7\text{\times}{10}^{12}\text{,}{\mathrm{m}}^{-1}$$, and this is the value used in Fig. 4E.

In Sec. S7 C we explained the observed quadrupole splitting $f_{Q}$ in the spectrum as caused by static strain. However, a small part of this should be due to the LQSE, as static electric fields are present, caused by the gate voltages required for device operation. Using the static electric field from our COMSOL model and the value of $R_{14}$ we found, we estimate an electric contribution to $f_{Q}$ of $1.7\text{\,}\mathrm{kHz}$, i.e. $\sim 2.5\%$ of the experimentally observed $f_{Q}$. This suggest that indeed the largest fraction of $f_{Q}$ is caused by strain.

In what follows we aim to justify the value of $R_{14}$ on the basis of a microscopic theory.

3 Comparison to LQSE measurements and empirical microscopic theory for bulk crystals

We have conjectured that the DC electric field induced shift of the NER spectral lines is a manifestation of the LQSE, predicted in 1961 by Bloembergen: “If a nucleus or paramagnetic ion with spin greater than one-half occupies a site in a crystal lattice without inversion symmetry, an external applied electric field will cause, in general, a change in the quadrupole coupling constant or crystal-line field splitting proportional to E”  [?]. The effect was first invoked to explain contemporary observations of electric field induced shifts in the zero-field quadrupolar spectra of ionic crystals containing halogen nuclei (81Br  [?] and 35Cl  [?,?]). Subsequently, it was measured for 27Al in Al2O3  [?] and 69Ga, 71Ga, and 75As in GaAs  [?,?]. The physical mechanism by which the external applied electric field couples to the nuclear quadrupolar interaction is through a distortion of the chemical environment of a given nucleus, which gives rise to a change in the EFG at the nucleus. Typical measured proportionalities between the external applied electric field and the resulting EFG are on the order of $1\text{\times}{10}^{12}\text{\,}{\mathrm{m}}^{-1}$.

The manifestation of LQSE in our experiment is distinguished from previous measurements not only by being observed in a single nucleus, but one in which the nuclear quadrupole moment is between 2 (75As) and 8 (35Cl) times larger than any previous measurements and in a host material that is non-polar in bulk. We stress that the necessary condition for LQSE of broken point inversion symmetry is naturally fulfilled by silicon; the $T_{d}$ symmetry of Si lacks point inversion symmetry at the lattice sites. Comparisons to the GaAs measurements are the most relevant because the individual nuclei experience the same tetrahedral coordination as our ${}^{123}\mathrm{Sb}$ donor. In the preceding section (Sec. S7 D 2) we estimated a value of $R_{14}=$1.7\text{\times}{10}^{12}\text{,}{\mathrm{m}}^{-1}$$ for the ${}^{123}\mathrm{Sb}$ substitutional donor in silicon under study here. This is comparable to the values reported for 75As in bulk GaAs  [?,?], respectively $1.55\text{\times}{10}^{12}\text{\,}{\mathrm{m}}^{-1}$ and $2.0\text{\times}{10}^{12}\text{\,}{\mathrm{m}}^{-1}$.

The microscopic theory used to rationalize the bulk GaAs measurements of Gill and Bloembergen [?] decomposes the electric field response tensor into ionic and covalent contributions, $R_{14}=R_{14}^{\mathrm{ion}}+R_{14}^{\mathrm{cov}}$. The ionic contribution accounts for the electric field gradient due to an electric field induced distortion of the crystal lattice (i.e. piezoelectricity), whereas the covalent contribution accounts for the electric field gradient due to an electric field induced distortion of the valence orbitals. In the case of 75As, their theory predicts that roughly 2/3 of the measured value of $R_{14}$ is covalent, while the remaining 1/3 is ionic, i.e., $R_{14}^{\mathrm{cov}}/R_{14}\approx 2/3$ and $R_{14}^{\mathrm{ion}}/R_{14}\approx 1/3$. Because the silicon host crystal is non-polar, we anticipate that our inferred value of $R_{14}$ is almost entirely covalent, $R_{14}\approx R_{14}^{\mathrm{cov}}$. Thus, in comparing the strength of the LQSE in our experiment to bulk GaAs, we compare to Gill and Bloembergen’s covalent value of $1.1\text{\times}{10}^{12}\text{\,}{\mathrm{m}}^{-1}$. The nucleus-specific constant of proportionality between the LQSE splitting and the external applied electric field is

[TABLE]

Accounting for the difference in nuclear quadrupole moment, nuclear spin, and $R_{14}$, we find $K^{\mathrm{cov}}$ values of $5.8\text{\times}{10}^{-18}\text{\,}\mathrm{m}$ for the covalent contribution to 75As in GaAs and $2.8\text{\times}{10}^{-18}\text{\,}\mathrm{m}$ for ${}^{123}\mathrm{Sb}$ in silicon. That these effects have the same order of magnitude suggests that it is reasonable to rationalize the DC bias induced shift in the NER spectral lines with the LQSE. In Sec. S7 D 2 it is also demonstrated that this same picture can rationalize the observed $\Delta m_{I}=\pm 1$ and $\Delta m_{I}=\pm 2$ Rabi frequencies.

It is also informative to directly compare the values of $R^{\mathrm{cov}}_{14}$ in the context of the atomistic model Gill and Bloembergen provide for $R^{\mathrm{cov}}_{14}$ in GaAs. Applying second-order perturbation theory to derive a linear coupling between the external applied electric field and the LQSE shift, they arrive at the empirical expression

[TABLE]

where $\lambda\in\left[0,1\right]$ is related to the ionicity of the Ga-As bond, $eq_{at}/(4\pi\epsilon_{0})$ is the EFG at the 75As nucleus due to an electron in a 4p orbital in an atomic model (i.e., not a solid), $|\psi_{c}\rangle$ and $W_{c}$ are the electronic orbital and orbital energy associated with a bonding orbital near the conduction band edge, and $|\psi_{v}\rangle$ and $W_{v}$ are the corresponding quantities for an antibonding orbital near the valence band edge. For 75As in GaAs, $\lambda=0.4$ corresponding to $\pm 2$ effective charges on the nuclei, $q_{at}=$40.8\text{\times}{10}^{30}\text{,}{\mathrm{m}}^{-3}$$, a transition dipole moment of $1.3\text{\,}\mathrm{eV}\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, and an orbital gap of $2\text{\,}\mathrm{eV}$ are found to be in good agreement with the measured value of $R_{14}$ when combined with a model for $R_{14}^{\mathrm{ion}}$. Because the symmetry of the ${}^{123}\mathrm{Sb}$ site is the same, an empirical model for the $R_{14}\approx R_{14}^{\mathrm{cov}}$ measured in our experiment would take a form similar to Eq. (53). Our inferred value of $R^{\mathrm{cov}}_{14}$ is only $\sim 1.5$ times larger than the one that was measured for 75As. That these values are so close in value can then be rationalized in terms of these empirical parameters, namely the ionicity of the Si-Sb+ bond, the EFG at the ${}^{123}\mathrm{Sb}$ nucleus due to electrons in p orbitals, and the electronic structure of the bonding and antibonding orbitals at the defect site.

While we expect that the Si-Sb+ bonding and antibonding orbitals will have some ionic character, we do not expect that they will have a more ionic character than the relevant Ga-As orbitals. This expectation is verified by preliminary quantum chemistry calculations on cationic hydrogen-terminated Sb-doped crystalline Si clusters using the NWChem software package  [?]. Studying finite clusters allows us to circumvent the technical difficulties associated with applying the modern theory of polarization in charged supercells. These calculations indicate that the relevant transition dipole moment is on the order of $0.1\text{\,}\mathrm{eV}\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, noting that the comparable transition dipole matrix element used by Gill and Bloembergen is an order of magnitude larger. Because our inferred value of $R^{cov}_{14}$ is the same order of magnitude as for GaAs, this suggests that $q_{at}$ needs to be commensurately larger to compensate. We expect that $q_{at}$ will, in fact, be larger because it will be the EFG at the ${}^{123}\mathrm{Sb}$ nucleus due to an electron in a 5p orbital, instead of a 4p orbital. The Sternheimer effect can also be invoked for even further increases, particularly due to the positive ionization state of the ${}^{123}\mathrm{Sb}$ and the observation that this effect grows with atomic number.

It is worth noting again that the sensitivity of our inference of $R_{14}$ to details in the COMSOL model suggests that we should interpret it as an order of magnitude estimate. Given the dependencies on empirical quantities described above, the microscopic physics seems to corroborate that this estimate is quite reasonable. Further experiments and the development of a first principles theory that does not rely on empirical expressions of the form in Eq. 53 will provide even more clarity. Such a theory is a topic of ongoing work. While there are reasonably mature first principles methods for computing the electric field gradient with and without strain (see Sec. S7 C 2), the direct evaluation of the LQSE parameters has not yet been demonstrated. Such a theory would circumvent the need for a description in terms of empirical parameters like bond ionicities, atomic EFG s, and Sternheimer factors. Another item of interest for a more detailed theory is the possibility of a second-order coupling of the EFG to both the strain and electric field, due to the development of a permanent electric dipole moment under shear strains as alluded to in Sec. S7 B.

E Alternative (unlikely) sources of NER

To further corroborate the claim that our data is most convincingly explained by an electric-field induced modulation of the nuclear quadrupole tensor caused by a local distortion of the atomic bonds, we quantitatively analyze two possible alternative explanations for the origin of a dynamical EFG. Both alternatives explanations are found to be in discrepancy with the observed NER rate by many orders of magnitude.

1 Mechanical driving through SiO2 piezoelectricity

Ref.  [?] suggests that SiO2 thin films thermally grown on Si can exhibit piezoelectric properties in the first few monolayers of the SiO2 film. In the presence of a piezoelectric material, the time-dependent part of the applied electric potential will result in a periodic deformation of the sample. Consequently, such a time-dependent strain component, which relates to a time-dependent EFG (Eq. (36)), results in a mechanical component to $\delta Q_{\alpha\beta}$ in Eq. (9). This would enable electrical driving of the nuclear spin.

Ref.  [?] models the mechanical displacement by assuming that the first $3.5\text{\,}\mathrm{nm}$ of oxide consist of quartz. We therefore incorporate this quartz layer at the Si/SiO2 interface in our COMSOL model. As for the silicon, we need to adjust the orientation of the piezoelectric and stiffness matrices. The crystal structure of $\alpha$-quartz is trigonal trapezohedral. Here, the [100]Quartz, [010]Quartz and [001]Quartz quartz crystal directions do not correspond to the Cartesian axes $x^{\prime\prime}$, $y^{\prime\prime}$ and $z^{\prime\prime}$ as for the silicon ([010], [001] and [100] respectively, see Sec. S7 A 1). In the $x^{\prime\prime}-y^{\prime\prime}$ plane, the quartz unit cell has the shape of an equilateral parallelogram ($a=b=$4.91\text{,}\mathrm{\SIUnitSymbolAngstrom}), where the angle between the crystal directions [100]Quartz and [010]Quartz is $\gamma=$120\text{\,}\mathrm{\SIUnitSymbolDegree}. In this convention, the standard silicon Cartesian $x^{\prime\prime}$ and $y^{\prime\prime}$ axes correspond to the [010]Quartz and [210]Quartz crystal directions, respectively. Perpendicular to these ($\alpha=\beta=$90\text{,}\mathrm{\SIUnitSymbolDegree}) extends the [001]Quartz direction of the unit cell with a length of $c=$5.40\text{\,}\mathrm{\SIUnitSymbolAngstrom}. This is the exceptional three-fold rotation axis of $\alpha$-quartz. Ref.  [?] found that the growth direction of $\alpha$-quartz on a (100) silicon substrate, using our crystal axis, is $[210]_{\rm Quartz}\parallel[001]_{\rm Silicon}$. Hence in our model the exceptional [001]Quartz axis will be in plane and aligned with the [100]Silicon or [010]Silicon directions. We show the first orientation as an example. We tested different quartz orientations, which all gave the same order of magnitude result.

A time-varying potential $\delta V_{\rm DFR}=$20\text{,}\mathrm{mV} applied to the DFR gate results in the driving strength shown in Fig. [15](#F15). At the expected donor position, the estimated Rabi frequency on resonance is 5 to $10\text{\,}\mathrm{Hz}$, two orders of magnitude smaller than the experimentally observed value $f^{\rm Rabi,NER}_{5/2\leftrightarrow 7/2}=$684\text{\,}\mathrm{Hz}. We therefore conclude that piezoelectricity is an unlikely cause of NER.

2 Direct effect of electric gate potentials

The electric gate potentials will directly lead to an EFG at the donor site, since the produced field is spatially inhomogeneous. To estimate the scale of this effect, the resulting quadrupole splitting is calculated in the absence of any other EFG generating effects such as strain or the LQSE. Fig. 16 shows the resulting quadrupole splitting $f_{Q}$ (a), the spectral line shift per applied volt on gate DFL (b) and NER driving strength $f^{\rm Rabi,NER}_{5/2\leftrightarrow 7/2}$ for $20\text{\,}\mathrm{mV}$ driving strength on the gate DFR (c), taking the EFG directly generated by the gates at the donor site. Their strengths are in the $100\text{\,}\mathrm{\SIUnitSymbolMicro Hz}$, $1\text{\,}\mathrm{mHz}\text{\,}{\mathrm{V}}^{-1}$ and $1\text{\,}\mathrm{\SIUnitSymbolMicro Hz}$ range, respectively, as opposed to the experimentally observed values of the order of $100\text{\,}\mathrm{kHz}$, $10\text{\,}\mathrm{kHz}\text{\,}{\mathrm{V}}^{-1}$ and $1\text{\,}\mathrm{kHz}$. Although some enhancement of this EFG is expected at the ${}^{123}\mathrm{Sb}$ nucleus due to the Sternheimer anti-shielding effect, this is not expected to surpass two orders of magnitude. As the direct effect of the electric gates results in an underestimation of the experimental observations by at least 6 orders of magnitude, we conclude this effect is insignificant. This illustrates that in our experiment a significant quadrupole interaction and NER strength requires a microscopic mechanism to be at play in direct vicinity of the donor site.

S8 Spin-mechanical coupling

A fascinating prospect raised by our results is the possibility to coherently couple a single nuclear spin to a mechanical resonator. This would require that the dynamical strain caused by the zero-point fluctuations of the mechanical beam results in a nuclear Rabi frequency exceeding the inhomogenous linewidths of both the electron and the mechanical resonator. Fortunately, both systems can have exceptionally narrow linewidths, in the range of a few Hz. Below we provide an estimate of the zero-point strain in a silicon doubly-clamped mechanical oscillator, following the calculation presented in Ref.  [?] (see also Ref.  [?] for a similar calculation in the case of a singly-clamped beam).

Consider a mechanical beam of length $L$, width $w$ and thickness $t$, clamped at both its extremities. We assume that the long axis of the beam coincides with the [110] crystallographic direction. For Si, the Young’s modulus is $E=188$ GPa for stress along [110], and the mass density is $\rho=2300$ kg/m3. The moment of inertia of the beam is $I=wt^{3}/12$ and the wave number of the fundamental resonance mode is $k_{0}=4.73/L$. From this, one derives the (angular) frequency of the fundamental mode as:

[TABLE]

In practice, one will need to choose $f_{0}$ to match the nuclear Larmor frequency $\gamma_{n}B_{0}$ (we neglect the small quadrupole shift $f_{Q}$ in this context). Therefore, $f_{0}$ and the resulting choice of beam geometry are not all free parameters. Choosing for example to fix the thickness $t$ of the beam (typically to the lowest value allowed by fabrication), we can derive the length $L$ necessary to yield the mechanical resonance that matches the nuclear precession frequency:

[TABLE]

Next, we calculate the zero-point mechanical strain $\epsilon_{\rm zpf}$ in such a beam. We assume that the nucleus is placed in the middle of the beam ($z=L/2$), as close as possible to the surface (where the strain is maximum), and find:

[TABLE]

This equation shows that, at a constant frequency, $\epsilon_{\rm zpf}$ is maximized by making the beam as thin and narrow as possible, which will result in a minimal length too [Eq. (55)]. It also shows that $\epsilon_{\rm zpf}$ depends on frequency (and thus on field) to the power of $3/4$.

Now we can use our models, benchmarked against the LQSE experimental data on the 123Sb nucleus, to estimate the spin-elastic coupling. In other words, we analyze a hypothetical Nuclear Acoustic Resonance (NAR) experiment  [?] where a dynamical strain $\delta\epsilon$ at the nuclear site is caused by the zero-point strain $\epsilon_{\rm zpf}$ of the resonator. We assume that the length of the beam, the z-axis, coincides with the [110] crystallographic axis of Si. We apply a magnetic field along the [100] crystal axis. Considering the $m_{I}=5/2\leftrightarrow 7/2$ transition, the Rabi frequency is given by:

[TABLE]

where $\alpha_{5/2\leftrightarrow 7/2}=\sqrt{63}$. Taking the value $S_{44}=6.1\times 10^{22}$ V/m2 extracted from the analysis in Sec. S7 C 2, we obtain:

[TABLE]

Identifying the dynamical strain component along the beam, $\delta\epsilon_{zz}$, with the zero-point strain of the mechanical beam, we find a simple linear relation:

[TABLE]

In the context of cavity-QED, $f_{5/2\leftrightarrow 7/2}^{\rm Rabi,NAR}$ can be identified as the spin-mechanical vacuum Rabi splitting, $2g$. The interesting strong-coupling regime is achieved when $g$ exceeds both the phonon loss $\kappa$ and the qubit dephasing $\gamma$. In the notation and language of relevance to our system, this means $f_{5/2\leftrightarrow 7/2}^{\rm Rabi,NAR}\gg\Gamma_{n},\Gamma_{m}$, where $\Gamma_{n}=3.3$ Hz is the experimentally observed nuclear spin inhomogeneous linewidth, and $\Gamma_{m}=f_{0}/Q_{m}$ is the resonator linewidth, determined by its mechanical quality factor $Q_{m}$.

For the purpose of achieving strong coupling, the resonator design must be optimized considering both $\epsilon_{\rm zpf}$ and $\Gamma_{m}$. Typical Si resonators have a constant product $Q_{m}\cdot f_{0}\sim 10^{13}$ Hz determined by material properties  [?]. Therefore, $\Gamma_{m}\propto f_{0}^{2}$, whereas $g\propto f_{0}^{3/4}$ (56), indicating that the best chance to achieve strong coupling is by reducing the operating frequency until the lower bound $\Gamma_{m}\sim\Gamma_{n}$ is reached ($\Gamma_{n}$ is usually independent of frequency).

Taking for example $w=40$ nm, $t=20$ nm, and placing the 123Sb nuclear spin in a static field $B_{0}=1.5$ T, the resulting beam would need to have a resonance frequency $f_{0}=8.33$ MHz, a length $L=4.71$ $\mu$m, and would produce a zero-point strain $\epsilon_{\rm zpf}\approx 0.6\times 10^{-8}$, resulting in $f_{5/2\leftrightarrow 7/2}^{\rm Rabi,NAR}\approx 1.2$ Hz.

Assuming that the resonator has a quality factor $Q_{m}\approx 10^{6}$ at this frequency (consistent with $Q_{m}\cdot f_{0}\sim 10^{13}$ Hz), i.e. $\Gamma_{m}=8.3$ Hz, we thus see that the spin-cavity coupling $g$ is within less than a factor 10 of the dephasing rates.

Importantly, our experiment shows that the challenge in reaching strong coupling would not be set by the nuclear spin coherence, but rather by the zero-point strain and quality factor of the mechanical resonator, both of which can be significantly improved from the values used in the simple estimate above.

Combining novel high-$Q_{m}$ resonators with designs that explicitly maximize the zero-point strain, could bring the goal of coherent coupling between a single nuclear spin and a mechanical oscillator tantalizingly within reach. Achieving the single-phonon limit will be very challenging ($f_{0}\sim 10$ MHz is equivalent to $\sim 0.5$ mK phonon energy), but interesting physics can already be explored in hot cavities  [?].