# Valley population of donor states in highly strained silicon

B Voisin, K.S H Ng, J Salfi, M Usman, J.C Wong, A Tankasala, B.C Johnson, J.C Mccallum, L Hutin, B Bertrand, et al.

► To cite this version:

B Voisin, K.S H Ng, J Salfi, M Usman, J.C Wong, et al. Valley population of donor states in highly strained silicon. Mat.Quant.Tech., 2022, 2 (2), pp.025002. 10.1088/2633-4356/ac5d1d. hal-04204454

HAL Id: hal-04204454

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

Submitted on 26 Sep 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Materials for Quantum Technology

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

CrossMark

PAPER

## Valley population of donor states in highly strained silicon

OPEN ACCESS

RECEIVED

6 October 2021

REVISED

20 February 2022

ACCEPTED FOR PUBLICATION

11 March 2022

PUBLISHED

6 April 2022

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.

Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

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

B Voisin¹,²,* , K S H Ng³,⁴, J Salfi³,⁵ , M Usman⁶,⁷ , J C Wong⁸, A Tankasala⁹,

B C Johnson⁶,¹⁰ , J C McCallum⁶,¹⁰ , L Hutin¹¹ , B Bertrand¹¹, M Vinet¹¹,

N Valanoor⁸ , M Y Simmons³,⁹ , R Rahman³,⁹ , L C L Hollenberg⁶ and S Rogge³,*,¹⁰

1 Silicon Quantum Computing, Sydney, NSW 2052, Australia

2 School of Physics, Sydney, The University of New South Wales, Sydney, NSW 2052, Australia

3 Centre for Quantum Computation and Communication Technology, School of Physics, The University of New South Wales, Sydney, NSW 2052, Australia

4 5. Physikalisches Institut and Center for Integrated Quantum Science and Technology, Universität Stuttgart, 70569 Stuttgart, Germany

5 Department of Electrical and Computer Engineering, University of British Columbia, Vancouver, BC V6T 1Z4, Canada

6 Centre for Quantum Computation and Communication Technology, School of Physics, The University of Melbourne, Parkville, VIC 3010, Australia

7 School of Computing and Information Systems, Faculty of Engineering and Information Technology, The University of Melbourne, Parkville, 3010, Victoria

8 School of Materials Science and Engineering, The University of New South Wales, Sydney, NSW 2052, Australia

9 Electrical and Computer Engineering Department, Purdue University, West Lafayette, Indiana, United States of America

10 Centre for Quantum Computation and Communication Technology, School of Engineering, RMIT University, Melbourne, VIC 3001, Australia

11 Université Grenoble Alpes, CEA, LETI, 38000 Grenoble, France

* Authors to whom any correspondence should be addressed.

E-mail: benoit.voisin@unsw.edu.au and s.rogge@unsw.edu.au

Keywords: donors, strain silicon, valley population, STM

Supplementary material for this article is available online

## Abstract

Strain is extensively used to controllably tailor the electronic properties of materials. In the context of indirect band-gap semiconductors such as silicon, strain lifts the valley degeneracy of the six conduction band minima, and by extension the valley states of electrons bound to phosphorus donors. Here, single phosphorus atoms are embedded in an engineered thin layer of silicon strained to 0.8% and their wave function imaged using spatially resolved spectroscopy. A prevalence of the out-of-plane valleys is confirmed from the real-space images, and a combination of theoretical modelling tools is used to assess how this valley repopulation effect can yield isotropic exchange and tunnel interactions in the *xy*-plane relevant for atomically precise donor qubit devices. Finally, the residual presence of in-plane valleys is evidenced by a Fourier analysis of both experimental and theoretical images, and atomistic calculations highlight the importance of higher orbital excited states to obtain a precise relationship between valley population and strain. Controlling the valley degree of freedom in engineered strained epilayers provides a new competitive asset for the development of donor-based quantum technologies in silicon.

## 1. Introduction

Strain is a modification of a crystal's lattice constant, which alters the electronic band structure and its degeneracies by breaking crystal symmetries. This effect, which can be captured by first principles methods [1, 2], can modify the band energies, lift degeneracies, and reduce inter-band scattering [3]. These properties have been used over the decades in microelectronics to enhance the mobility of nanoscaled MOSFETs [3, 4] as well as the luminescence of solar cells or 2D material systems [5, 6]. By extension, strain also alters the quantum states of electrons bound to single dopants, defects or quantum dots [7] which can be isolated and controlled in nanoscale electronic devices, since these states are built upon band states of the host material. Strain is invariably present to some degree in nanoelectronic devices for classical and quantum technologies, notably because of the material stack, of the presence of interfaces and nearby doping. Therefore, it is essential to develop
<figure><img src="image_4.png" /></figure>

tools that are able to probe the impact of strain at the nanoscale in order to understand and mitigate, but also possibly control and leverage strain in new architecture designs [8–10]. In particular, strain lifts the valley and heavy/light hole degeneracies of silicon [11, 12] (see figure 1(a)), with direct implications for exchange and hyperfine interactions [11, 13–15] as well as dipole and possible quadrupolar moments [10, 16] of single dopants, all core properties that enable quantum information processing.

The sensitivity of electron spin resonance spectroscopy has been used to detect shifts in the hyperfine frequencies induced by strain on dopants in silicon [17–20]. However, this technique averages the impact of strain over microscopic areas and over a large number of dopants, while the local impact of strain in a nanostructure and on a single quantum state is of interest for quantum applications based on individual dopants. Instead, microscopy techniques are well suited to study strain effects at the nanoscale. Whilst scanning transmission electron microscopy [21] and atom probe tomography are able to detect single atoms in nanoscale devices [22, 23], scanning tunneling microscopy (STM) presents the advantage of giving direct access to wave function information [24]. This technique was recently used to probe the valley population of single donors in natural silicon [25, 26]. Here, we extend this work to spatially resolve the wave function of an electron bound to a single phosphorus atom embedded in a thin layer of highly strained silicon compatible with atomic precision lithography [27].

In this work, we report that the complex oscillating pattern which was observed in these images in the unstrained case due to the presence of all six valleys [25, 26] vanishes because of the valley repopulation effect. Using a combination of theoretical models for the exchange interaction, we show how this dominant *z*-valley donor wave function can be used in the context of atomically precise qubit devices to obtain robust exchange and tunnel couplings for any in-plane orientation and down to 5 nm separation. A Fourier transform analysis highlights the remaining presence of an *xy*-valley population, which is unexpected according to effective mass theories for such large strain [11]. The presence of this remaining in-plane valley population observed experimentally is in agreement with atomistic tight-binding (TB) simulations [28], which notably take into account high orbital excited states [29, 30] and have been benchmarked against *ab initio* density functional theory [14]. Our work provides a quantitative understanding of the effect of strain on valley population, relevant for the development of robust two-donor coupling in silicon.
## 2. Sample fabrication and structural analysis

The sample fabrication procedure achieves three main objectives. The native oxide needs to be removed without dewetting the epilayer nor relaxing the strain, a flat surface must be obtained to enable STM spectroscopy, and dopants must be embedded in the originally undoped strained epilayer to induce a single electron tunneling framework on isolated ones. The fabrication of our sample starts from a commercial 300 mm wafer from Soitec with a 10 nm-thick layer of highly strained silicon [31, 32]. This silicon layer was grown on a Si0.8Ge0.2 substrate to induce the desired amount of strain, and subsequently transferred on a silicon oxide substrate according to the Smart Cut process™. Importantly, the strain should be preserved throughout the fabrication procedure which is achieved by implementing a low-temperature fabrication budget.

The substrate was first implanted at 3.5 µm-depth with phosphorus atoms at an energy of 6 MeV and a dose of 10¹⁵ cm⁻² in order to facilitate the annealing steps required during the in situ fabrication procedure. Normally for such devices, this in situ preparation starts with a 10 h-anneal at 600 °C followed by a high-temperature flash annealing (typically above 1050 °C) to remove the native oxide [25]. However, this high temperature anneal cannot be performed here with a thin strained layer because it results in cross-hatching due to the strain variation across the wafer, which is notably unsuitable for STM lithography [33]. This high-temperature flash also incurs a high risk of dewetting the thin strained layer. Instead, we adopted a low-temperature recipe [8] where the sample is kept at 700 °C facing a silicon deposition source calibrated to a rate of 0.5 ML min⁻¹ for 1 h. This was followed with a 2 min-anneal at 800 °C to ensure a high-quality 2 × 1 surface reconstruction [see supporting information section S1 (https://stacks.iop.org/MQT/2/025002/mmedia)].

A schematic of the final device is shown in figure 1(b), with the corresponding electronic transport diagram shown in figure 1(c). Similar to previous work [25, 26, 34, 35], a tunnel current occurs from the highly doped reservoir labelled $P_1$ to an isolated donor found in the lightly doped layer (labelled $P_2$) and then to the STM tip. The temperature of the fabrication process, starting from the initial step where the native oxide is removed, is described in figure 1(d). A saturation dose of phosphorus dopants $P_1$, which will be conductive and act as a reservoir at low temperature, is deposited and incorporated at the surface. This reservoir is electrically contacted and voltage biased ($V_b$) through the silicided part of the sample under the clamp of the STM sample holder. A silicon growth of 7.5 nm follows, divided in three equivalent sub-growths of 2.5 nm. For each sub-growth, the first nanometer is grown at room temperature and the remaining 1.5 nm grown between 250 °C and 400 °C, which is designed to minimise dopant segregation whilst preserving a flat surface for STM studies [36, 37]. A second phosphorus deposition $P_2$ follows, with a lower dose of 10¹² cm⁻², designed to be able to randomly find isolated donors. A second silicon growth of 2.5 nm buries the isolated donors, followed by a 10 s-anneal at 600 °C to further flatten the surface, which is finally hydrogen passivated. A large scale topography image is shown in figure 1(e), which shows an excellent surface quality, with low defects and large terraces (larger than 20 nm) suitable for STM spectroscopy and amenable to STM lithography to fabricate atomically precise qubit devices [27, 38].

To confirm that we successfully inhibited strain relaxation by limiting the thermal and growth budget, we performed x-ray diffraction measurements through the overgrown layer. A 2θ – ω x-ray 2D map was taken along the [224] crystallographic axis for the bare sSOI wafer first, which is reconstructed in real space in figure 2(a). The lattice coordinates were corrected in order for the map maximum (yellow cross) to match the lattice constant of the sSOI strain-relaxed handling wafer in all three dimensions. The amount of germanium in the virtual substrate Si1-xGex determines the resulting amount of strain in the silicon epilayer, which can be predicted from linear elasticity theory using the lattice constant mismatch between the two materials ($a_{Ge} = 0.5658$ nm and $a_{Si} = 0.5431$ nm). The resulting in-plane $a_‖$ and out-of-plane $a_⊥$ silicon lattice constants with respect to the unstrained case are:

$$ 
\begin{gathered}
a_‖ = (1 - x)a_{Si} + x a_{Ge} \\
a_⊥ - a_{Si} = 2(c_{12}/c_{11})(a_{Si} - a_‖) \sim 0.75(a_{Si} - a_‖)
\end{gathered}
\qquad (1) $$

where $c_{11}$ and $c_{12}$ are the elastic constants of relaxed silicon [39]. Since $a_{Ge} > a_{Si}$, growing a silicon epilayer on Si1-xGex results in a tensile in-plane strain, i.e. $a_‖ > a_{Si}$ and $a_⊥ < a_{Si}$. We define $\epsilon_‖$ the in-plane strain (i.e., the relative change in the in-plane lattice constant), with $\epsilon_‖ = 0.8\%$ for $x = 0.2$. The black dashed line represents the linear relationship between the in-plane $a_‖$ and out-of-plane $a_⊥$ lattice constants following equation (1). The thin strained silicon epilayer is evidenced as an ellipsoidal signal, which is centred very close to the expected location for a silicon layer grown on a 20%-rich germanium substrate with $a_‖ = 0.5476$ nm and $a_⊥ = 0.5397$ nm (red cross). A 2D Gaussian fit of this signal (black ellipse, see supporting information section S2) yields a strain value $\epsilon_‖^{bare} = 0.76 \pm 0.18\%$ (black cross). A similar measurement and analysis was performed on the overgrown and dosed sample, as shown in figure 2(b). A similar signal can be observed for this sample close to the expected strain value, and the same fitting analysis yields $\epsilon_‖^{exp} = 0.71 \pm 0.15\%$. The reservoir doping
<figure><img src="image_5.png" /><figcaption>**Figure 2.** X-ray diffraction of thin strained silicon epilayers. (a) 2D $2\theta - \omega$ map taken along [224], reconstructed in real space, of the bare wafer. The black dashed line corresponds to the expected relationship between $a_\parallel$ and $a_\perp$ from elasticity theory, see equation (1). The yellow cross indicates the position of the maximum, which corresponds to the unstrained substrate and is used as a reference of the relaxed lattice constants. The strained layer appears as a 2D ellipse centred (black cross) very close to the expectation value for a layer grown on a $Si_{0.8}Ge_{0.2}$ wafer (red cross). The black ellipse denotes the variance of a Gaussian fit of the strain ellipse. The other diagonal features are spurious interference from the measurement setup. (b) Similar x-ray diffraction measurements and analysis were performed on the sample that was spectroscopically measured at low temperature after growth and doping. The strain has remained after overgrowth and phosphorus incorporation with very limited relaxation.</figcaption></figure>

concentration is on the order of $10^{20}-10^{21}$ cm$^{-3}$, which induces a strain of the order of 0.05% [40]. The strain obtained here is therefore much larger and can be confidently attributed to the strain originally present in the epilayer. The errors bars between the two samples are also quantitatively similar, which demonstrates that the strain is maintained during our anneal, growth and dosing procedure. We note that the SiGe growth method used to fabricate the original strained epilayer is known to result in a strain field which can vary by around 10% with a characteristic length of 1 to 2 µm due to the cross-hatch effect reflecting the misfit dislocations present in the original SiGe buffer [41]. The implications for single donors will be addressed later in the manuscript.

### 3. Spectroscopy and wave function imaging of a strained donor

The objectives of this section are to spectroscopically evidence the presence of single donors addressable in the single electron tunneling regime, before spatially resolving their charge distribution and analysing these images in the context of the strain applied to the crystal. The sample is voltage biased ($V_b$) and the STM tip is connected to the input of a trans-impedance amplifier placed at room temperature. In this configuration, donor states can be evidenced spectroscopically at 4.2 K as they yield resonances below the direct transport onset from silicon valence band states to the STM tip [25, 26, 42]. This onset can be measured away from a single donor as shown in figure 3(a). Here, it occurs at $V_b \sim 0.2$ V, which is offset compared to previous results [25, 26, 42] in unstrained silicon where it is found around $V_b \sim -0.8$ V. This offset is attributed to
<figure><img src="image_6.png" /></figure>

an electrical short between the silicon epilayer and the handling wafer on the sample holder. Two extra reso-
nances, at $V_b \sim 0.6$ V and $V_b \sim 0.2$ V respectively, can be observed when tunneling over a single dopant occurs
(red line), which correspond to adding respectively the first (called $D^0$ state) and the second electron
(called $D^-$ state) on the donor, the latter being close to the onset of direct transport from the valence band
states [25, 34]. This double resonance can be fitted to a thermal broadening model from which lever arms of
$0.065 \pm 0.009$ and $0.072 \pm 0.003$ are obtained for the $D^0$ and $D^-$ resonances, respectively, as well as a charging
energy of $14 \pm 2$ meV (see supporting information section S3) between the two charge states. This value is
low compared to the charging energy of a single phosphorus donor in bulk silicon (about 47 meV [27]). This
difference can be explained by the presence of a metallic reservoir 10 nm away [35] and the strain environment
which leads to a slight increase in the wave function envelope extent. We have performed full configuration
interaction (FCI) calculations of the $D^-$ state following an atomistic tight-binding simulation of $D^0$. These
calculations take into account both the presence of the strain in the lattice and of a screening metallic reservoir
[30], and find that the charging energy of a $4.5a_0$-deep donor is decreased to 23 meV for 20%-Ge content
strain at zero externally applied electric field, in reasonable agreement with the experimental value mentioned
above.

A double-pass technique [25] is used to image the charge distribution of the single donor evidenced in spec-
troscopy. The tip is first stabilised in topography, before setting the bias in the gap in a second pass where only
the donor density of states contribute to the tunnel current. Following a drift correction protocol performed
on the lattice topography image of the first pass, the charge distribution images are shown in figures 3(b) and
(c), respectively for the $D^0$ state taken at $V_b = 0.4$ V and for the $D^-$ state taken at $V_b = 0.2$ V. In the unstrained
case, the donor image reflects interference between waves at the valley and lattice frequencies, and the pattern
strongly depends on the lattice plane at which the donor sits [25, 26, 43, 44]. In contrast, the strained case
studied here shows that the maxima of the charge distribution are aligned along the dimer rows. Following
our previous work [26], we have computed theoretical STM images based on a multi-million atom TB simu-
lation of the donor $D^0$ state at a depth of $4.5a_0$ for several amounts of strain. At 1% (in Ge content), as shown
in figure 3(d), the image reveals a complex symmetry with a set of oscillating patterns whose maxima are distributed both along and across the dimer rows. This is characteristic of the presence of all six valleys and of lattice frequencies. Going to 5% and 20%, as shown respectively in figures 3(e) and (f), this oscillating pattern progressively vanishes leaving a signal aligned with the surface dimer rows in agreement with the experimental images. This behaviour is consistent for any dopant depth (see supporting information section S4). Two other donors measured in different sample locations show very similar features with charge distribution maxima aligned along the dimer rows (see supporting information section S5). Together with the support from theoretical modelling, this analysis suggests that the strain remains large enough anywhere on the sample to induce a significant change in the valley population of the donors.

## 4. Valley population impact on tunnel and exchange coupling

The change in the pattern observed in the STM images links to a change in the valley population. This effect has important consequences for dopant electronic properties, which we investigate theoretically in this section. We focus here on interactions based on the direct overlap between donor wave functions, such as tunnel and exchange couplings. These quantities can be very sensitive to the exact position of the dopants in the lattice because of the presence of interference at the valley frequency, an effect modulated by the anisotropy of the donor envelope that originates from the mass anisotropy in the conduction band of silicon [42, 45]. Below, we investigate the effect of the valley population on these interactions in the context of atomically precise donor qubit devices, where donors are placed within the same crystallographic plane [13, 27]. We start from the Heitler–London expression of the exchange coupling using an effective mass description of the single electron wave functions, and of the relationship between tunnel and exchange couplings [42, 45]:

$$ J(\vec{R}) = \sum_{\substack{\mu,\nu= \\ \pm\{x,y,z\}}} \alpha_{\mu}^{2} \alpha_{\nu}^{2} j_{\mu\nu}(\vec{R}) \cos(\Delta\phi_{\mu\nu}(\vec{R})) \\ t = \sqrt{JU}/2 \qquad (2) $$

where $\alpha_\mu$ is the amplitude of the valley population, $j_{\mu\nu}$ the envelope weight or overlap between the orbital part of valleys $\mu$ and $\nu$, $\Delta\phi_{\mu\nu}$ a geometric phase difference term between the two valleys which depends on the distance between the two atoms and the valley momentum, and $U$ the on-site Coulomb repulsion energy of the 2e-state which we consider to be constant. This exchange expression is in good agreement with full configuration exchange calculations down to inter-donor distances of 5 nm [30, 42] (see supporting information section S6). In equation (2), the valley population gives an extra weight to each interfering envelope term $j_{\mu\nu}$. Their impact is shown in figure 4(a), where all the tunnel coupling values corresponding to in-plane donor positions between 6 and 8 nm are plotted for two different valley populations. In the unstrained case, i.e. a z-valley population of 33% ($\alpha_z = \alpha_x = \alpha_y = 1/\sqrt{6}$), clear variations of over an order of magnitude close to the [100] axis can be observed, which are due to the sensitivity to y-valley interference, while the [110] axis is more resilient [42]. Going to a z-valley population of 80%, these large variations vanish to make the tunnel coupling isotropic in the plane, and we also note a relative increase of the maximum values. Both effects arise since the wave function is mainly concentrated in the z-valley and the $j_{zz}$ term is fully constructive, where less amplitude of the wave function is lost through destructively interfering terms in building the overlap. Making the direct interactions isotropic in the plane opens a way for quantum simulation architectures based on donor arrays with non-square geometries, where the donors are not only positioned along the [110] or [110] axis [10, 46, 47].

We further investigated the impact of strain on the exchange variations caused by small changes in atom position. For two target orientations ([100] and [110]) and two target distances (5 and 12 nm), we calculated the tunnel couplings values corresponding to all possible inter-donor positions within a 1 nm in-plane neighbourhood. The ratios of the minimum to the maximum value are plotted in figure 4(b) as a function of the total z-valley population. Along [110] the variations remain below a factor of 10 for small target distance and change in valley population. However, large variations are predicted along the [100] axis due to the inherent presence of destructive positions close by, for any target distance. Strain provides a way to overcome these variations. Overall, the impact of valley interference on tunnel coupling, and thus exchange coupling as well, is quenched for a z-valley population beyond 70% leaving variations only due to changes in the donor envelope overlap. A remaining question to address is to determine the amount of strain to achieve such population, which is the focus of the following section.
a Target distance 6-8nm

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

b Variations around a target distance

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

**Figure 4.** Valley population and tunnel coupling robustness to dopant position. (a) Normalised tunnel coupling values calculated for all the in-plane inter-donor positions between 6 and 8 nm, for a z-valley population of 33% (black points) and 80% (red points). The angle $\theta_{in}$ is defined from the [100] crystallographic direction. The tunnel coupling becomes isotropic upon strain with a reduction of the large variations close to the [100] crystallographic axis. (b) Ratios of the minimum to maximum tunnel coupling values found within a ±1 nm in-plane neighbourhood as a function of the z-valley population, for two target distances (5 and 12 nm) and two target directions ([100] and [110]). Valley interference-induced variations are quenched irrespective of the distance and orientation for a z-valley population beyond 70%.

**5. Valley population of strained donors - Fourier analysis**

The valley and lattice frequencies that are present in the donor's ground state image can be distinguished in Fourier space around the first Brillouin zone [44]. The 2D Fourier transform of the experimental images of the $D^0$ and $D^-$ states are shown in figures 5(a) and (b) respectively. A strong signal can be observed around the lattice frequencies located at $(k_x, k_y) = \pm0.5k_0$ and $(k_x, k_y) = \pm1k_0$ with $k_0 = 2\pi/a_0$ (grey dashed circles), which reflects the dominance of the z-valley population in the donor ground state and the alignment of the real-space image signal along the surface dimer rows. The same strong lattice components can also be observed in the Fourier transform of the theoretical $D^0$ image at 20%-Ge strain (figure 5(c)) to establish a clear agreement between experiment and theory. Regarding the *xy*-valleys, important components to focus on are located around $k_x, k_y \sim 0.81k_0$ (green ellipses in figures 5(a)–(c)), since they relate to *xy*-valley interference only and are independent of the lattice [44]. The $D^-$ image also provides information on the $D^0$ bound states including valleys [35]. While the *xy*-valley population is expected to be very low in such large strain environment, a clear signal can be observed around the *y*-valleys of the $D^0$ Fourier transform, as well as for both *x* and *y*-valleys of the $D^-$ experimental and the $D^0$ theoretical Fourier transforms, which establishes a residual presence of *xy*-valley population. The asymmetry between the *x* and the *y*-valley signal observed in the experimental $D^0$ Fourier transform, with the *x*-valley signal absent, is most likely due to an extrinsic cause to the dopant (lattice defect or tip effect for instance), but further investigation is needed. The presence of a residual *xy*-valley population observed experimentally and theoretically based on a TB approach is not expected from effective mass
<figure><img src="image_9.png" /></figure>

(EM) theory for such large strain [11, 17] as explained below. Following this EM framework, an energy shift of the valley states is associated with a change of the lattice constants upon strain, which can be characterised by a valley strain parameter $\chi$ linked to the deformation potential $\Xi_u \sim 8.6$ eV and the valley-orbit parameter $\Delta_c \sim 2$ meV for silicon [11, 17]:

$$ \chi = \frac{\Xi_u}{3\Delta_c} \frac{a_{Si} - a_{Ge}}{a_{Si}} \left( \frac{2c_{12}}{c_{11}} + 1 \right) x. \quad (3) $$

This energy shift of the valley states causes the donor ground state to evolve from the so-called $A_1$ state, where all the valleys are equally distributed, to the $A_{12}^\infty$ state made out z-valleys only (see supporting information section S7). The evolution of the z-valley population according to the EM model is plotted in figure 5(d) as a function of the Ge concentration (see supplementary information), and it can be compared to the valley population obtained from our TB atomistic simulations. The z-valleys quickly dominate (above 90%) for strain larger than 3.2% in the EM case, to exceed 99.7% at 20%-Ge, while the valley repopulation is much slower in the case with a remaining 4.6% of xy-population at 20%-Ge, which qualitatively match our experimental observation.

As depicted in figure 1(a), the evolution of the valley population is directly related to the energy difference between the valley states, which is a good starting point to gain a better understanding of why the change in valley population as calculated by TB is slower than by EM. Therefore, we show in figure 5(e) the evolution of the state energies as a function of strain for both EM and TB calculations. Both models show a qualitatively similar behaviour for the 1s-manifold (black lines for EM, blue dots for TB), with a large first excited state energy (about 12 meV) at zero strain which decreases to a few meV at finite strain. Both models also show an increasing splitting between the z-valley and xy-valley states (black and blue arrows for the EM and TB model, respectively) as the Ge content (hence the strain) is increased. However, the energy splitting between
these valley states is smaller for TB than it is for EM calculations for the same amount of strain (see arrows), which drives a slower valley repopulation of the ground state. In order to understand this quantitative difference between the two models, it is important to note that the EM model only considers the first 1s-like manifold, while TB takes into account a larger donor spectrum basis including higher orbital states [29, 30] (e.g., 2s, 2p and so on). In figure 5(e), some of these higher orbital states (red dots) can be seen coming down in energy with strain and crossing the xy-valley states, which can reduce the energy gap between the z and xy-valley states in turn driving the valley population. As future work, it would be relevant to extend the EM model to include higher orbital states in the calculations to see if it converges towards the TB ones, as well as to further investigate in the TB framework the interaction and possible hybridization between the 1s-like and higher orbital states when they cross around 5% Ge concentration.

## 6. Conclusion

In conclusion, we have engineered a thin epilayer of highly strained silicon to embed and image the wave functions of single phosphorus donors. A low-temperature fabrication procedure was used to maintain a high level of strain, as confirmed by x-ray diffraction, and to obtain a high quality surface for STM purposes. We evidence a clear dominance of the z-valley population in the STM image in contrast to the relaxed case, which can be used to make the exchange and tunnel coupling isotropic in the xy-plane. A Fourier analysis highlights the residual presence of xy-valley states for large strain that goes beyond predictions based on EM theory. This effect can be explained using atomistic TB simulations showing the influence of higher orbital states, not previously considered in EM calculations. These results provide the necessary understanding of strain mechanisms and impact on donor quantum states in silicon to harness their potential in view of quantum nanoelectronics. Further work will focus on optimising the comparison between experimental and theoretical images to provide a better accuracy on the strain each dopant is embedded in and its local variations across the sample. Electrically separating the top and back contacts of the sample would also provide a way to modify the vertical electric field and therefore the valley population of single donors [48], which could be directly tracked using quantum state imaging.

## Acknowledgments

BV, KN, JS, MU, AT, BCJ, JCMC, MYS, RR, LCLH and SR acknowledge support from the ARC Centre of Excellence for Quantum Computation and Communication Technology (CE170100012). JS acknowledges support from an ARC DECRA fellowship (DE160101490). This work was in part funded by the US Army Research Office (W911NF-17-1-0202). BCJ and JCMC acknowledge the AFAiiR node of the NCRIS Heavy Ion Capability for access to ion-implantation facilities. CW and NV acknowledge support from ARC Centre of Excellence Future Low Energy Electronics Technologies Grant CE170100039. Computational resources were provided by the National Computing Infrastructure (NCI) and Pawsey Supercomputing Center through National Computational Merit Allocation Scheme (NCMAS) and from the LIEF HPC-GPGPU facility hosted at the University of Melbourne with the assistance of LIEF Grant LE170100200. MYS is a director of the company Silicon Quantum Computing Pty Ltd.

## Data availability statement

The data that support the findings of this study are available upon reasonable request from the authors.

## ORCID iDs

B Voisin https://orcid.org/0000-0002-5231-1965

J Salfi https://orcid.org/0000-0001-9240-4245

M Usman https://orcid.org/0000-0003-3476-2348

B C Johnson https://orcid.org/0000-0002-2174-4178

J C McCallum https://orcid.org/0000-0002-6692-7728

L Hutin https://orcid.org/0000-0001-6429-3867

N Valanoor https://orcid.org/0000-0003-2534-5868

M Y Simmons https://orcid.org/0000-0002-6422-5888

R Rahman https://orcid.org/0000-0003-1649-823X
L C L Hollenberg https://orcid.org/0000-0001-7672-6965

S Rogge https://orcid.org/0000-0003-1078-9482

## References

[1] Richard S, Aniel F, Fishman G and Cavassilas N 2003 *J. Appl. Phys.* **94** 1795

[2] Munguia J, Bremond G, Bluet J M, Hartmann J M and Mermoux M 2008 *Appl. Phys. Lett.* **93** 102101

[3] Niquet Y-M, Delerue C and Krzeminski C 2012 *Nano Lett.* **12** 3545

[4] Yu D, Zhang Y and Liu F 2008 *Phys. Rev.* B **78** 245204

[5] Wu W, Wang L, Yu R, Liu Y, Wei S-H, Hone J and Wang Z L 2016 *Adv. Mater.* **28** 8463

[6] Deng S, Sumant A V and Berry V 2018 *Nano Today* **22** 14

[7] Hanson R, Kouwenhoven L P, Petta J R, Tarucha S and Vandersypen L M K 2007 *Rev. Mod. Phys.* **79** 1217

[8] Zhang P, Tevaarwerk E, Park B-N, Savage D E, Celler G K, Knezevic I, Evans P G, Eriksson M A and Lagally M G 2006 *Nature* **439** 703

[9] Thorbeck T and Zimmerman N M 2015 *AIP Adv.* **5** 087107

[10] Salfi J, Mol J A, Culcer D and Rogge S 2016 *Phys. Rev. Lett.* **116** 246801

[11] Koiller B, Hu X and Das Sarma S 2002 *Phys. Rev.* B **66** 115201

[12] Zwanenburg F A, Dzurak A S, Morello A, Simmons M Y, Hollenberg L C L, Klimeck G, Rogge S, Coppersmith S N and Eriksson M A 2013 *Rev. Mod. Phys.* **85** 961

[13] He Y, Gorman S K, Keith D, Kranz L, Keizer J G and Simmons M Y 2019 *Nature* **571** 371

[14] Usman M, Hill C D, Rahman R, Klimeck G, Simmons M Y, Rogge S and Hollenberg L C L 2015 *Phys. Rev.* B **91** 245209

[15] Franke D P, Hrubesch F M, Künzl M, Becker H-W, Itoh K M, Stutzmann M, Hoehne F, Dreher L and Brandt M S 2015 *Phys. Rev. Lett.* **115** 057601

[16] Asaad S *et al* 2020 *Nature* **579** 205

[17] Wilson D K and Feher G 1961 *Phys. Rev.* **124** 1068

[18] Huebl H, Stegner A R, Stutzmann M, Brandt M S, Vogg G, Bensch F, Rauls E and Gerstmann U 2006 *Phys. Rev. Lett.* **97** 166402

[19] Lo C C, Urdampilleta M, Ross P, Gonzalez-Zalba M F, Mansir J, Lyon S A, Thewalt M L W and Morton J J L 2015 *Nature Mater* **14** 490

[20] Mansir J *et al* 2018 *Phys. Rev. Lett.* **120** 167701

[21] Cooper D, Le Royer C, Béché A and Rouvière J-L 2012 *Appl. Phys. Lett.* **100** 233121

[22] Koelling S, Li A, Cavalli A, Assali S, Car D, Gazibegovic S, Bakkers E P A M and Koenraad P M 2017 *Nano Lett.* **17** 599

[23] Ishikawa R, Shibata N, Taniguchi T and Ikuhara Y 2020 *Phys. Rev. Applied* **13** 034064

[24] Yakunin A M *et al* 2007 *Nature Mater* **6** 512

[25] Salfi J, Mol J A, Rahman R, Klimeck G, Simmons M Y, Hollenberg L C L and Rogge S 2014 *Nature Mater* **13** 605

[26] Usman M, Bocquel J, Salfi J, Voisin B, Tankasala A, Rahman R, Simmons M Y, Rogge S and Hollenberg L C L 2016 *Nat. Nanotechnol.* **11** 1

[27] Fuechsle M, Miwa J A, Mahapatra S, Ryu H, Lee S, Warschkow O, Hollenberg L C L, Klimeck G and Simmons M Y 2012 *Nature Nanotech* **7** 242

[28] Usman M, Rahman R, Salfi J, Bocquel J, Voisin B, Rogge S, Klimeck G and Hollenberg L C L 2015 *J. Phys.: Condens. Matter.* **27** 154207

[29] Lansbergen G P, Rahman R, Verduijn J, Tettamanzi G C, Collaert N, Biesemans S, Klimeck G, Hollenberg L C L and Rogge S 2011 *Phys. Rev. Lett.* **107** 136602

[30] Tankasala A *et al* 2018 *Phys. Rev.* B **97** 195301

[31] Hartmann J M *et al* 2003 *Semicond. Sci. Technol.* **19** 311

[32] Hartmann J M, Abbadie A, Guinche Y, Holliger P, Rolland G, Buisson M, Defranoux C, Pierrel F and Billon T 2007 *Semicond. Sci. Technol.* **22** 354

[33] Lee W C T, McKibbin S R, Thompson D L, Xue K, Scappucci G, Bishop N, Celler G K, Carroll M S and Simmons M Y 2014 *Nanotechnology* **25** 145302

[34] Voisin B, Salfi J, Bocquel J, Rahman R and Rogge S 2015 *J. Phys.: Condens. Matter.* **27** 154203

[35] Salfi J, Voisin B, Tankasala A, Bocquel J, Usman M, Simmons M Y, Hollenberg L C L, Rahman R and Rogge S 2018 *Phys. Rev.* X **8** 031049

[36] Keizer J G, Koelling S, Koenraad P M and Simmons M Y 2015 *ACS Nano* **9** 12537

[37] Koch M, Keizer J G, Pakkiam P, Keith D, House M G, Peretz E and Simmons M Y 2019 *Nature Nanotech* **14** 137

[38] Kranz L, Gorman S K, Thorgrimsson B, He Y, Keith D, Keizer J G and Simmons M Y 2020 *Adv. Mater.* **32** 2003361

[39] Windl W, Sankey O F and Menéndez J 1998 *Phys. Rev.* B **57** 2431

[40] Lee M, Ko E and Ko D-H 2017 *J. Mater. Chem.* C **5** 9744

[41] Rouchon D, Hartmann J-M, Crisci A and Mermoux M 2008 *ECS Trans.* **16** 181

[42] Voisin B, Bocquel J, Tankasala A, Usman M, Salfi J, Rahman R, Simmons M Y, Hollenberg L C L and Rogge S 2020 *Nat. Commun.* **11** 6124

[43] Kohn W and Luttinger J M 1955 *Phys. Rev.* **98** 915

[44] Saraiva A L, Salfi J, Bocquel J, Voisin B, Rogge S, Capaz R B, Calderón M J and Koiller B 2016 *Phys. Rev.* B **93** 045303

[45] Koiller B, Hu X and Das Sarma S 2001 *Phys. Rev. Lett.* **88** 027903

[46] Georgescu I M, Ashhab S and Nori F 2014 *Rev. Mod. Phys.* **86** 153

[47] Le N H, Fisher A J, Curson N J and Ginossar E 2020 *npj Quantum Inf* **6** 24

[48] Verduijn J, Tettamanzi G C and Rogge S 2013 *Nano Lett.* **13** 1476