Linear hyperfine tuning of donor spins in silicon using hydrostatic strain

We experimentally study the coupling of Group V donor spins in silicon to mechanical strain, and measure strain-induced frequency shifts which are linear in strain, in contrast to the quadratic dependence predicted by the valley repopulation model (VRM), and therefore orders of magnitude greater than that predicted by the VRM for small strains $|\varepsilon|<10^{-5}$. Through both tight-binding and first principles calculations we find that these shifts arise from a linear tuning of the donor hyperfine interaction term by the hydrostatic component of strain and achieve semi-quantitative agreement with the experimental values. Our results provide a framework for making quantitative predictions of donor spins in silicon nanostructures, such as those being used to develop silicon-based quantum processors and memories. The strong spin-strain coupling we measure (up to 150~GHz per strain, for Bi-donors in Si), offers a method for donor spin tuning --- shifting Bi donor electron spins by over a linewidth with a hydrostatic strain of order $10^{-6}$ --- as well as opportunities for coupling to mechanical resonators.

The exploitation of donor spins in silicon as qubits typically requires their incorporation into nano-and micro-electronic devices.This has been used to demonstrate single-shot readout of a single 31 P donor spin using a tunnel-coupled silicon single-electron transistor (SET) [21,22], and to create hybrid devices in which donor spins are coupled to superconducting resonators [23][24][25] to develop interfaces between microwave photons and solid-state spins.In both cases, the use of metal-oxide-semiconductor (MOS) nanostructures [26], or patterned superconducting films on silicon [27] involves a combination of materials with coefficients of thermal expansion that differ by up to an order of magnitude [28][29][30][31][32].The presence of strain in the silicon environment around the donor spin is therefore pervasive when studying such nanodevices at cryogenic temperatures.Furthermore, factors such as optimising spinresonator coupling or spin-readout speed motivate the placement of donors close to features such as SETs [33] or resonator inductor wires [34] where strain is maximal.
Strain modifies the band structure of silicon [35,36], as has been shown, for example, to contribute to the confinement of single electrons in silicon under nanoscale aluminium gates [37,38].The donor electron wavefunction is also modified by strain: following the valley repopulation model (VRM) developed by Wilson & Feher [22] within the framework of effective mass theory, an applied uniaxial strain lifts the degeneracy of the six silicon valleys leading to a mixture of the donor ground state, 1s(A 1 ), with the first excited state, 1s(E).In this excited state, the hyperfine coupling between the donor electron and nuclear spin is zero, therefore the VRM predicts a quadratic reduction in A as a function of uniaxial strain, as well as a strain-induced anisotropic contribution to the electron g-factor.Strain-induced perturbations of the donor hyperfine coupling have been observed for P-donor spins in 28 Si epilayers, grown on SiGe to yield built in strains of order 10 −3 [40] -piezoelectric contacts on such material have been used to modulate this built-in strain to shift the electron spin resonance frequency by up to ∼ 400 kHz [41].
In this Letter, we report the observation of a straininduced shift in the hyperfine coupling of group V donors in silicon which is linear (rather than quadratic), and therefore orders of magnitude greater than that predicted by the valley repopulation model of Wilson and Feher [22] for small strains (|ε| < ∼ 10 −5 ).We present experimental studies showing strain-tuning of the coherent evolution of each of the group V donor spins ( 31 P, 75 As, 121 Sb, and 209 Bi), extracting the strain-induced shifts of the hyperfine coupling and electron spin g-factor for each, and corroborate the results with a combination of both tight binding and density functional theory calculations which reveal the crucial role of hydrostatic strain in this novel mechanism [42].In addition to providing essential insights for the use of donor spins in nano-and micronscale quantum devices, our results provide a method for linear tuning of the donor hyperfine interaction with coupling strengths of up to 150 GHz/strain for 209 Bi donor spins.
The spin Hamiltonian for a group V donor in the presence of an external magnetic field B = B 0 z is: where g e and g n are, respectively, the electronic and nuclear g-factors, µ B and µ N are the Bohr and nuclear magnetons, Ŝ and Î are the electronic and nuclear spin operators.The Fermi contact hyperfine interaction strength, A = 1.4754GHz in Si:Bi, can be expressed as: where ψ(0) represents the amplitude of the electronic wavefunction at the nucleus and g 0 = 2.0023 is the free electron g-factor.The eigenstates of this Hamiltonian describe a Hilbert space of dimension (2S + 1)(2I + 1), with S = 1/2 and I determined by the nuclear spin species, illustrated for the case of 209 Bi:Si in Fig. 1(A).Transitions between these eigenstates obeying the selection rule (∆m S , ∆m I ) = (±1, 0) in the high field limit can be driven and detected using pulsed electron spin resonance (ESR) [43].
We use samples of isotopically enriched 28 Si doped with Bi, Sb, As and P (see Ref [44] for more details), mounted with crystal orientation shown in Fig. 1(B).The sample is situated inside a dielectric microwave ESR resonator in an Oxford Instruments CF935 liquid helium flow cryostat, and is held between two PEEK rods whose end faces have been CNC-milled to match the profile of the sample allowing it to be rotated with respect to the magnetic field.Using calibrated masses, a uniaxial stress is applied to the sample perpendicular to the [110] face, via the upper rod which extends outside the cryostat.The resulting strain tensor can be derived from the generalised form of Hooke's law for anisotropic materials and the compliance matrix for silicon [45]: in the ([110],[1 10],[001]) coordinate system per kg of applied mass, ε 11 = −1.45× 10 −5 /kg, ε 22 = 9.02 × 10 −7 /kg, ε 33 = 5.24 × 10 −6 /kg, and ε i =j = 0.While the VRM predicts frequency shifts only from uniaxial strain, we shall see that the new mechanism presented here arises from hydrostatic strain ε hs = (ε 11 + ε 22 + ε 33 )/3.In our setup, we estimate a strain per unit mass ε hs = −2.78× 10 −6 /kg.
We use a home-built pulsed ESR spectrometer [46]  π → τ → echo [47], with τ = 15 µs and a π pulse duration of 130 ns.The time-domain Hahn echo signals (top of Fig. 2) obtained while systematically increasing the applied strain are then Fourier transformed to yield the strain-induced shifts in spin resonance frequency [44].
First, we observe in Fig. 2 that the Bi donor ESR transition can be shifted by more than a linewidth (in 28 Si) for strains of order 10 −5 (uniaxial) or 10 −6 (hydrostatic).We fit the frequency-domain echo signals to a Voigt profile, and then plot the centre-frequency shifts as a function of strain, for each of the ten allowed ESR transitions (see Fig. 1(C)).Strikingly, the ESR frequency of each transition shows a linear dependence on strain, rather than the expected quadratic dependence.In Fig. 3, we plot the experimentally determined ∂f /∂ 11 for each transition against the first-order sensitivity of each transition frequency to the isotropic hyperfine coupling ∂f /∂A.Remarkably, all 10 points fall on a single line, demonstrating that the dominant effect we observe in Si:Bi is a straininduced shift in the isotropic hyperfine coupling which is linear in strain, and equivalent to ∂A/∂ε 11 = 5.4 ± 0.3 GHz or ∂A/∂ε hs = 28.2 ± 1.6 GHz.Multivalley effective mass theory (EMT) has been successful in describing many aspects of the donor electron wavefunction [48][49][50], including close agreement between theory and experimental measurements of the Stark effect [51] and predictions of exchange coupling between neighbouring donors [52].Within this framework, the wavefunction is expanded in terms of Bloch functions concentrated around the six degenerate [100] conduction band minima (valleys) such that ψ = 6 µ=1 α µ F µ φ µ where µ indexes over the valleys in the basis [+x, −x, +y, −y, +z, −z], F µ is a hydrogenlike envelope function, and φ µ is the valley Bloch function.The donor impurity potential breaks the symmetry of the crystal and induces a coupling between the valleys, leading to a valley-orbit splitting of the 1s-like donor state into three sub-levels.The ground state is singly degenerate with A 1 symmetry and has α µ = 1/ √ 6, while one of the excited states is doubly degenerate with E symmetry and has . The valley repopulation model assumes that uniaxial strain applied along a valley axis results in the corresponding pair of valley energies being decreased or increased for compressive or tensile strain, respectively [22].This modification of the valley energies results in a redistribution of the amplitude of each valley contributing to the ground state, which can be represented under strain as an admixture of the 1s(A 1 ) and 1s(E), resulting in a quadratic reduction of A as a function of uniaxial strain.At our maximum applied (A) The gradient of the strain-induced frequency shifts (df /d ) for each of the ten Si:Bi ESR transitions are shown as a function of the first-order sensitivity of each transition to the hyperfine coupling (∂f /∂A).The linear relationship confirms the observed strain-induced shifts in Si:Bi result from tuning the hyperfine coupling, A, with a gradient ∂A/∂ 11 = 5.4 ± 0.3 GHz or, equivalently, ∂A/∂ hs = 28.2±1.6GHz.(B) Calculations showing the relative change in hyperfine coupling strength A/A0 as a function of strain, comparing tight-binding (TB) and valley repopulation (VRM) models.To mimic the experiment, we show TB calculations of A/A0 under uniaxial stress along [110] (purple curve and circles), which produces a hydrostatic component of strain (ε hs , top axis) in addition to a uniaxial component along [110] (ε11, bottom axis).This behaviour can be understood by comparing with TB calculations for pure hydrostatic stress (blue curve and hexagons), plotted on the same axis of ε hs , as well as calculations from the VRM (red dotted curve), plotted on the same axis of ε11.The arrows indicate the relevant axes for each trace.
strain of ε 11 = −1.45× 10 −5 , the VRM predicts a reduction in A of 1.9 kHz, while we measure a reduction in A of 78 kHz -this discrepancy is even more pronounced for smaller strains.Therefore, in addition to predicting a different functional form of the dependence of A against strain, the VRM predicts shifts which are approximately two orders of magnitude smaller than what we measure in this strain regime, implying that another physical mechanism must dominate the changes to the structure of the donor electron wavefunction we observe.
In order to understand these trends, we have computed the bound states of bismuth impurities in silicon using the sp 3 d 5 s * tight-binding (TB) model of Ref. 2. This model reproduces the variations of the band structure of bulk silicon under arbitrary strains in the whole first Brillouin zone.The impurity is described by a Coulomb tail and by a correction of the orbital energies of the bismuth atom (similar to a central cell correction in the effective mass approximation) [4].The TB ratio A/A 0 between the strained (A) and unstrained (A 0 ) hyperfine interaction strengths is plotted in Fig. 3(B) under uniaxial stress along [110], as a function of the resulting uniaxial [110] and hydrostatic strains.Surprisingly, and in agreement with the experiments performed here, A/A 0 behaves linearly with small strain, and this trend can be assigned to the effects of the hydrostatic stress.Although not predicted by the VRM, the existence of a linear hydrostatic term is compatible with the symmetries of the system [42].A symmetry analysis indeed suggests that, to second order in the strains ε ij in the cubic axis set: A fit to the TB data yields K = 29.3,L = −9064 and N = −225.L mostly results from the coupling of the 1s(A 1 ) with the 1s(E) state by the uniaxial strain.The TB L is close to the VRM L = −2Ξ 2 u /(9∆ 2 ) = −9720 [22], where Ξ u = 8.6 eV is the uniaxial deformation potential of the conduction band of silicon and ∆ = 41 meV is the splitting between the 1s(A 1 ) and the 1s(E) state of the Bi impurity.The quadratic shear term N is usually negligible with respect to L. K = ∂(A/A 0 )/∂ε hs results from the coupling of the 1s(A 1 ) with the 2s(A 1 ) state (and higher A 1 states, since hydrostatic strain preserves the symmetry of the system) due to the change of the shape and depth of the central cell correction under strain.A/A 0 is dominated by this hydrostatic term at small strain, as evidenced in Fig. 3(B).The TB K = 29.3 is larger than the experimental K = 19.1.At variance with L (which mostly depends on a deformation potential of the silicon matrix), K indeed depends on details of the potential near the impurity, which must be specifically accounted for in the TB model in order to reach quantitative accuracy [44].In order to better capture the central cell correction around the bismuth impurity, we also performed first principles calculations using density functional theory (DFT) to describe the atomic relaxations not accounted for by our TB calculations [44].The DFT calculations further corroborate the linear dependence of the hyperfine coupling on hydrostatic strain (for hs ≤ 10 −3 ), and predict a coefficient K = 17.5, in good agreement with our experiments.Full details concerning the models and calculations can be found in Ref [44].
To test our model further and explore the expected anisotropy of a g-factor coupling to strain, we extend our study over a range of magnetic field orientations (as defined in Fig. 1(B)) and for the other group V donors: 31 P, 75 As, and 121 Sb.In all cases we find the observed ESR transition frequency shifts f are linear as a function of hydrostatic strain hs , with the resulting coupling strengths (df /d hs ) summarised in Fig. 4 and Table SI, along with values predicted from tight-binding calculations and full data sets [44].While we find no significant anisotropy in Si:Bi, the data from Si:Sb, Si:As, and Si:P display strain effects which clearly depend on the magnetic field orientation, attributed to a strain-induced anisotropic electronic g-factor.Following Wilson & Feher [22], our model for this anisotropy includes a term accounting for the effect of valley repopulation, and another accounting for the effect of spin-orbit coupling in the sheared lattice.Fits of this model to our experimental data reproduce the predicted strength of both of these effects to within a factor of two [44].
Through experiments and calculations, we have demonstrated that hydrostatic strain in silicon leads to a strong, linear tuning of the hyperfine interaction in group V donors, through coupling between the 1s(A 1 ) and 2s(A 1 ) states.The ability to shift the ESR transition frequencies by over a linewidth with hydrostatic strain in the order of 10 −6 opens up new possibilities for conditional "A-gate" control of donors as well as coupling to mechanical resonators.In addition, these insights will be crucial in supporting the design of quantum memories and processors based on donors in silicon, enabling the Supplementary Materials: Linear hyperfine tuning of donor spins in silicon using hydrostatic strain
A cylindrical aluminium plate rests on the floor of the cryostat (see Fig. S1).A rod made from the plastic PEEK screws into this plate and extends into the centre of the sapphire resonator, providing a bottom support for the sample.A second PEEK rod, which is supported radially by the resonator structure but is free to move along its axial direction, holds the sample in place from above.
Each echo is averaged 300 times and the spins are reset after each cycle with a 5 ms flash of 50 mW above-bandgap 1047 nm laser light through an optical window in the cryostat.
The time-domain Hahn echo signals can be expressed as a periodic oscillation at the frequency of the detuning between the fixed microwave drive frequency and the strain-shifted transition frequency, multiplied by an envelope function given by the inverse Fourier transform of the frequency-domain transition spectral lineshape.Then, by the convolution theorem, the Fourier transform of such a signal results in a frequency-domain spectral peak centred at the detuning frequency.The strain-induced detuning is then extracted by fitting a Voigt profile to this spectral peak.See Ref. [S1] Chapter 5 for more details.
In the valley repopulation model (VRM) [S10], A(ε ) is expected to be quadratic with small ε (the changes in A being exclusively driven by the loss of symmetries).In the TB approximation, A(ε ) indeed describes a parabola for weak stress, but centered on some ε > 0. Therefore, A(ε ) appears to behave almost linearly with small compressive ε .
To understand this trend, it is very instructive to split the strain into a hydrostatic component, an uniaxial component, and a shear component.The hydrostatic component εhs is defined as ε xx = ε yy = ε zz = ε hs .It accounts for the changes in the total volume Ω (∆Ω/Ω = 3ε hs ).The uniaxial and shear components account for the changes in symmetries (at constant volume).The uniaxial component εuni is defined as where: For uniaxial stress along [110] [Eqs.(S6)], where: A(ε) is plotted as a function of εhs , εuni , and εshear in Fig. S4.A(ε) shows a quadratic behavior as a function of ε uni and ε shear .It does, however, behave linearly as a function of ε hs .The trends at small ε < 0 evidenced in Fig. S3 can, therefore, be ascribed to the effects of hydrostatic strains on the hyperfine coupling constant.As a matter of fact, A(ε ) = A(ε hs ) + A(ε 001 uni ) for uniaxial [001] stress, and A(ε ) = A(ε hs ) + A(ε 110 uni ) + A(ε 110 shear ) for uniaxial [110] stress, showing the relevance of this decomposition.

Discussion
Although not predicted by the valley repopulation model, the existence of a ∝ ε hs term in A is allowed by symmetries [S11].Indeed, a symmetry analysis suggests that, to second order in ε ij : where K = ∂A/∂ε hs , L, M and N are constants.A fit to the tight-binding data yields: Eq. ( S12) then simplifies into: We find L −M because there is no sizable non-linearity in the dependence of A on hydrostatic strain in the investigated range |ε hs | < 10 −4 (no ∝ ε 2 hs term above when ε xx = ε yy = ε zz ).Also note that the effects of the quadratic shear terms are usually negligible with respect to the effects of the quadratic uniaxial terms (N L).The quadratic L term is mostly due to the coupling of the 1s(A 1 ) with the 1s(E) states of the impurity under uniaxial strain.The VRM of Ref. S10 actually suggests L = −2Ξ 2 u /(9∆ 2 ), where Ξ u = 8.6 eV is the uniaxial deformation potential of the conduction band of silicon and ∆ is the splitting between the 1s(A 1 ) and the 1s(E) state.For Bi (∆ = 41 meV), the VRM predicts L = −9720, in close agreement with the TB data.The linear hydrostatic term is due, on the other hand, to the coupling of the 1s(A 1 ) with 2s(A 1 ) state (and possibly higher s states with the same symmetry).This coupling results from the variations of the on-site correction on the bismuth impurity, and from the variations of the bismuth-silicon interactions under hydrostatic strain -in other words, from the variations of the depth and shape of the "central cell correction" [S12] not accounted for by the VRM.
Indeed, the total potential on the bismuth atom catches contributions from the tails of the atomic potentials of the neighboring silicon atoms, and these contributions depend on the silicon-bismuth bond lengths.As a matter of fact, the present TB model includes a strain-dependent correction for the energy E iµ of orbital µ ≡ s, p, d, s * of atom i [S2]: where the sum runs over the nearest neighbors j of atom i, d ij is the distance between atoms i and j, and d 0 ij is the relaxed bond length.E 0 iµ is the energy of the orbital in the reference, unstrained system and α iµ characterizes the deepening of the potential under strain.The interactions between bismuth and the nearest neighbor silicon atoms scale, on the other hand, as (d ij /d 0 ij ) n , with n close to 2. In the present TB model, the parameters α iµ and the exponents n of bismuth are the same as for silicon.Although this choice is a safe first guess, it can only provide a semi-quantitative description of the dependence of A on hydrostatic strain.This is why the TB K = 29.3 is significantly larger than the experimental K = 19.1.We may, in the spirit of Eqs.(S2), lump all corrections to the TB model into the α iµ .Therefore, we tentatively set: and adjust ∆α on the experimental K.This yields ∆α = 3.48 eV for bismuth.
We have repeated the same procedure for P, As and Sb.We give in Table SI the on-site parameters U s = U s * , U p , U d and ∆α of each impurity, as well as the values of K, L and N .The model for P and As [S13] is based on a simple Coulomb tail V (R i ) = −e 2 /(κR i ) instead of Eq. (S1).
The dependence of the binding energy E b of As impurities on the hydrostatic pressure P has been measured by Holland and Paul (dE b /dP −0.05 meV/kbar) [S14] and by Samara and Barnes (dE b /dP −0.1 meV/kbar) [S15].The electron hence gets more loosely bound to the impurity under pressure (or equivalently under compressive hydrostatic strain).This is consistent with the decrease of the hyperfine coupling constant A reported here (K < 0).Samara and Barnes explain the decrease of E b under pressure by the variations of the effective masses and dielectric constant κ.We point out, though, that there is also a significant contribution from the variations of the central cell correction.The variations of effective masses and dielectric constant actually make little contribution to K. In the simplest effective mass approximation, the wave function of the electron bound to the donor is indeed Ψ( r) = e −r/a B /( √ πa , where a B = h2 κ/(m * e 2 ) is the Bohr radius.Hence, |Ψ( 0)| 2 = 1/(πa 3 B ), so that: Ab-initio calculations within density functional theory (see next section) give (1/κ)(∂κ/∂ε hs ) = 0.78, (1/m * )(∂m * /∂ε hs ) = −0.19 for the longitudinal mass, and (1/m * )(∂m * /∂ε hs ) = 1.54 for the transverse mass.Therefore, the variations of the masses and dielectric constant are expected to have little net effect on the hyperfine coupling constant.

DENSITY FUNCTIONAL THEORY CALCULATIONS
In order to strengthen the above interpretation, we have also performed first principles calculations using density functional theory (DFT) with the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional [S16] and the projector-augmented wave method [S17] in the Vienna Ab-initio Simulation Package (VASP) [S18], following the methodology described in Ref. [S19].The calculations were carried out on one Bi impurity in a 1728-atom supercell, with a 250 eV plane-wave cutoff energy.DFT describes the central cell correction around the bismuth impurity from first principles and captures the atomic relaxations not accounted for by TB calculations -due to its accuracy in the immediate vicinity of the donor it could be expected to provide a good description of the variation of the hyperfine coupling with strains.Nevertheless, due to the finite size of the supercell, DFT misses the long range Coulomb tail of the potential [S20], which (along with over-delocalization arising from the self-interaction error in PBE) contributes to an significant underestimation in the absolute value the hyperfine interaction (1102 MHz).
The ab-initio hyperfine coupling shows the expected linear dependence on hydrostatic strain over the entire range explored here (up to ε xx = ε yy = ε zz = 2 × 10 −3 ).We extract a coefficient K = 17.5, in good agreement with the experimental data.The ab-initio quadratic term is L = −1.17× 10 4 , as obtained from a fit to the calculated data with ε zz ≤ 10 −3 and ε xx = ε yy = 0. Deviations are observed for higher strains, consistent with the higher-order terms present in the VRM.The equilibrium bismuth-silicon bond length (2.651 Å) is significantly larger than the Si-Si bond length (2.367 Å in the bulk).All Bi-Si and Si-Si bonds are simply scaled by the hydrostatic strain, to within better than 0.001 Å.
These data provide further support for the linear dependence of the hyperfine parameter on the hydrostatic component of strain, and illustrate the complementary strengths of these two computational approaches (DFT and tight-binding) in modelling the behaviour of donors in silicon.

MODELING G-FACTOR ANISOTROPY
The ellipsoidal shape of the Si conduction band minima in k-space results in differing effective masses for Bloch waves parallel and perpendicular to the valley axis at these points [S21].This leads to a single-valley g-factor which is anisotropic, with g 2 µ = g 2 cos 2 φ + g 2 ⊥ sin 2 φ, where g and g ⊥ are the g-factors parallel and perpendicular to the valley axis, and φ is the angle between the magnetic field and the valley axis.In the case of a donor, the g-factor is found by summing over the relative contribution from each valley such that g donor = 6 µ=1 α 2 µ g µ .For the unstrained 1s(A 1 ) donor ground state, which is an equal superposition of all six equivalent valleys, this summation leads to a cancellation of the anisotropy, leaving g 0 = 1 3 g + 2 3 g ⊥ .Under strain, the valleys repopulate, breaking this cancellation symmetry.In our system with θ defined as in figure S1, the resulting g-factor anisotropy can be modeled using the VRM [S22] by: where: g and g ⊥ are the parallel and perpendicular g-factors, Ξ u is the uniaxial deformation potential, ∆ is the donordependent 1s(A 1 )-1s(E) splitting, and k 1 and k 2 are defined in Eq.S11.
It is also known that the g-factor of a single valley is changed in the presence of shear strain by the coupling with the opposite valley.Following Wilson & Feher [S22], the effective Hamiltonian for this mechanism can be written [S23]: where C is a coefficient involving spin-orbit coupling matrix elements [S22].Rewriting the corresponding g-factor contribution in terms of ε 11 results in a second anisotropic term: where: Then, the derivative of the ESR transition frequency with respect to the uniaxial strain reads: ) ∂f /∂g can be calculated for each transition in the same manner as ∂f /∂A by solving the spin Hamiltonian while varying the value of g.We use a linear least squares regression to fit this model to the experimental data.We introduce three fitting parameters characterising the strength of the different effects for each donor: K = ∂(A/A 0 )/∂ hs , β VRM and β shear .The results of these fits are reported in table SII and figure 4 in the main text for Si:Sb, Si:As, and Si:P.The absence of a clear anisotropy for Si:Bi could be explained by the relatively small magnitude of the predicted g-factor effects in comparison with the absolute shifts due to the modified hyperfine interaction.
β VRM × 10 −3 (exp.) β shear × 10 −3 (exp.)We compare the extracted parameters with theoretical predictions for β VRM and β shear [Eqs.(S19) and (S22)] calculated using C = 0.44 [S22], ∆ values from Ref. S8, and g − g ⊥ values from Ref. [S22].The measured strengths of the g-factor effects agree with theory within approximately a factor of 2. Tight binding simulations with the present model (see also Ref. S24) predict that the very small g − g ⊥ is approximately an order of magnitude larger than given by Ref. [S22].It is interesting to note that the experimental data sit between Ref. [S22] and the TB predictions.

at 9 . 7 6 FIG. 1 .
FIG. 1. (A)Energy levels of Si:Bi spin eigenstates as a function of magnetic field, with the ten allowed ESR transitions at 9.7 GHz highlighted and labelled according to their high-field nuclear spin projection mI .(B) Schematic of experimental setup.The silicon single crystal sample is mounted between two engineered plastic rods with the ability to apply compressive stress.θ, the angle of the applied magnetic field to the [001] direction, can be varied by rotating the sample.(C) Observed linear frequency shifts for each of the ten allowed ESR transitions shown in panel (A), as a function of strain ε11, with θ = 30 • .

FIG. 2 .
FIG. 2. Electron spin echo signals in the frequency domain measured in 28 Si:Bi as a function of compressive strain (shown in terms of the uniaxial strain ε11 and hydrostatic strain ε hs ) arising from the applied stress in our experiment.Timedomain echoes for zero strain and ε11 = 1.4 × 10 −5 are shown as insets.Data shown is from the mI = −1/2 transition with θ = 45 • , taken at T = 8 K.
FIG. 3.(A) The gradient of the strain-induced frequency shifts (df /d ) for each of the ten Si:Bi ESR transitions are shown as a function of the first-order sensitivity of each transition to the hyperfine coupling (∂f /∂A).The linear relationship confirms the observed strain-induced shifts in Si:Bi result from tuning the hyperfine coupling, A, with a gradient ∂A/∂ 11 = 5.4 ± 0.3 GHz or, equivalently, ∂A/∂ hs = 28.2±1.6GHz.(B) Calculations showing the relative change in hyperfine coupling strength A/A0 as a function of strain, comparing tight-binding (TB) and valley repopulation (VRM) models.To mimic the experiment, we show TB calculations of A/A0 under uniaxial stress along [110] (purple curve and circles), which produces a hydrostatic component of strain (ε hs , top axis) in addition to a uniaxial component along [110] (ε11, bottom axis).This behaviour can be understood by comparing with TB calculations for pure hydrostatic stress (blue curve and hexagons), plotted on the same axis of ε hs , as well as calculations from the VRM (red dotted curve), plotted on the same axis of ε11.The arrows indicate the relevant axes for each trace.

FIG. 4 .
FIG.4.Extracted linear fit gradients df /d 11 for each transition for all four donors under consideration as a function of the angle of B0 w.r.t. the crystal θ.For Si:Sb, Si:As, and Si:P, these fits are overlaid with a model taking into account the linear shift of hyperfine interaction strength A as well as an anisotropic g-factor as a function of 11.
FIG. S1. a) Schematic of experimental setup showing Bruker ESR resonator mounted inside liquid Helium flow cryostat.The sample is held in place in the centre of the sapphire resonator by two engineered plastic rods.b) Magnified view of sample inside strain mount.

FULL 5 FIG. S5 .
FIG. S5.Full dataset for Si:Bi showing frequency shifts for all ten ESR transitions as a function of strain ε11 and θ.Linear fits are shown for each transition.

FIG. S6 .
FIG. S6.Full dataset for Si:Sb showing frequency shifts for all six ESR transitions as a function of strain ε11 and θ.Linear fits are shown for each transition.
FIG. S7.Full dataset for Si:As showing frequency shifts for all four ESR transitions as a function of strain ε11 and θ.Linear fits are shown for each transition.

FrequencyFIG. S8 .
FIG. S8.Full dataset for Si:P showing frequency shifts for both ESR transitions as a function of strain ε11 and θ.Linear fits are shown for each transition.

TABLE SI .
Coulomb tail Us = Us * (eV) Up (eV) U d (eV) ∆α (eV) Nature of the Coulomb tail, on-site corrections U [Eq. (S2)] and ∆α [Eq.(S16)], and value of K, L and N for P, As, Sb and Bi donors in silicon.