# Spin-orbit enhancement in Si/SiGe heterostructures with oscillating Ge
concentration

Benjamin D. Woods  Department of Physics, University of Wisconsin-Madison,
Madison, WI 53706, USA  M. A. Eriksson  Department of Physics, University of
Wisconsin-Madison, Madison, WI 53706, USA  Robert Joynt  Department of
Physics, University of Wisconsin-Madison, Madison, WI 53706, USA  Mark Friesen
Department of Physics, University of Wisconsin-Madison, Madison, WI 53706, USA

###### Abstract

We show that Ge concentration oscillations within the quantum well region of a
Si/SiGe heterostructure can significantly enhance the spin-orbit coupling of
the low-energy conduction-band valleys. Specifically, we find that for Ge
oscillation wavelengths near $\lambda=1.57~{}\text{nm}$, a Dresselhaus spin-
orbit coupling is produced that is over an order of magnitude larger than what
is found in conventional Si/SiGe heterostructures without Ge concentration
oscillations. We also provide a detailed explanation for this resonance
phenomenon. This involves the Ge concentration oscillations producing
wavefunction satellite peaks a distance $2\pi/\lambda$ away in momentum space
from each valley, which then couple to the opposite valley through Dresselhaus
spin-orbit coupling. Our results indicate that the enhanced spin-orbit
coupling can enable fast spin manipulation within Si quantum dots using
electric dipole spin resonance in the absence of micromagnets. Indeed, our
calculations yield a Rabi frequency
$\Omega_{\text{Rabi}}/B>500~{}\text{MHz/T}$ near the optimal Ge oscillation
wavelength $\lambda=1.57~{}\text{nm}$

##  I Introduction

Following the seminal work of Loss and DiVincenzo [1], quantum dots in
semiconductors have emerged as a leading candidate platform for quantum
computation [2, 3, 4, 5]. Gate-defined quantum dots in silicon [6, 7] are
particularly attractive due to their compatibility with the microelectronics
fabrication industry. Moreover, in contrast with GaAs, for which coherence
times are limited by unavoidable hyperfine interactions with nuclear spins
[8], isotropic enrichment dramatically suppresses these interactions in Si,
enabling long coherence times [9].

While recent progress in Si quantum dots has been quite promising, many of the
leading qubit architectures rely on synthetic spin-orbit coupling arising from
micromagnets [10, 11, 12, 13], leading to challenges for scaling up to systems
with many dots. An alternative approach is to use intrinsic spin-orbit
coupling for qubit manipulation, for example, through the electric dipole spin
resonance (EDSR) mechanism [14, 15]. While this possibility has been
considered for Ge and Si hole-spin qubits, where the degeneracy of the
p-orbital dominated valence band leads to strong spin-orbit coupling [16, 17],
the weak spin-orbit coupling of the Si conduction band appears unconducive for
electron-spin qubits.

Figure 1: (a) Schematic of the Si/SiGe heterostructure considered here,
consisting of a Si-dominated quantum well sandwiched between Si0.7Ge0.3
barrier regions. (b) Ge concentration profile along the growth ($z$) direction
of a Wiggle Well. Ge concentration oscillations of wavelength $\lambda$ inside
of the quantum well region lead to spin-orbit coupling enhancement for a
proper choice of $\lambda$. (c) Ge concentration profile of a “convetional”
Si/SiGe quantum well, for comparison.

In this work, we show how the spin-orbit coupling in Si/SiGe quantum well
heterostructures can be enhanced by more than an order of magnitude by
incorporating Ge concentration oscillations inside the quantum well, leading
to the possibility of exploiting intrinsic spin-orbit coupling in Si quantum
dots for fast gate operations. Fig. 1(a) shows a schematic of the system which
consists of a Si-dominated quantum well region sandwiched between
$\text{Si}_{0.7}\text{Ge}_{0.3}$ barrier regions. In contrast to
“conventional” Si/SiGe quantum wells, the quantum well regions contains a
small amount of Ge with concentration oscillations of wavelength $\lambda$, as
shown in Fig. 1(b). For comparison, Fig. 1(c) shows the Ge concentration
profile of a conventional Si/SiGe quantum well. Previous works [18, 19] have
studied such a structure, which has been named the Wiggle Well, and found that
the periodic Ge concentration leads to an enhancement of the valley splitting.
Here, we develop a theory of spin-orbit coupling within such structures and
show that the periodic nature of the device, along with the underlying diamond
crystal structure and degeneracy of the Si $z$ valleys, also gives rise to an
enhancement in spin-orbit coupling. Importantly, we find that the wavelength
$\lambda$ must satisfy a resonance condition to give rise to this spin-orbit
coupling enhancement. As discussed in detail in Sec. IV, this involves a two
step process that can be summarized as follows. First, the periodic potential
produced by the Ge concentration oscillations produces wavefunctions
satellites a distance $2\pi/\lambda$ away in momentum space from each valley.
Then, a satellite of a given valley couples strongly to the opposite valley
through Dresselhaus spin-orbit coupling, provided that the satellite-valley
separation distance in momentum space is $4\pi/a$, corresponding to the
condition $\lambda=1.57~{}\text{nm}$.

From the outset, it is important to remark that the spin-orbit coupling
introduced by the Ge concentration oscillations is fundamentally distinct from
the spin-orbit coupling of conventional Si/SiGe quantum wells. For a given
subband of a conventional Si/SiGe quantum well immersed in a vertical electric
field, the $C_{2v}$ point group symmetry of the system allows for both Rashba
and “Dresselhaus-like” linear-$\mathbf{k}_{\parallel}$ terms of the form [20,
21],

$H_{\text{SO}}=\alpha\left(k_{y}\bar{\sigma}_{x}-k_{x}\bar{\sigma}_{y}\right)+\beta\left(k_{x}\bar{\sigma}_{x}-k_{y}\bar{\sigma}_{y}\right),$
(1)

where $\bar{\sigma}_{j}$ are the Pauli matrices acting in (pseudo)spin space
and $\alpha$ and $\beta$ are the Rashba and Dresselhaus coefficients,
respectively, of the subband. The presence of Rasbha spin-orbit coupling is
unsurprising due to the structural asymmetry provided by the electric field
[22], while the presence of the Dresselhaus-like term is initially surprising
since the diamond lattice of Si/SiGe quantum wells possesses bulk inversion
symmetry [23]. However, these systems still support a Dresselhaus-like term of
the same form $\beta(k_{x}\bar{\sigma}_{x}-k_{y}\bar{\sigma}_{y})$, due to the
broken inversion symmetry caused by the quantum well interfaces [24, 20, 21,
25]. This is in stark contrast to the true Dresselhaus spin-orbit coupling in
III-V semiconductors, where the asymmetry of the anion and cation in the unit
cell leads to bulk inversion asymmetry [23]. Importantly, we find in Sec.
III.2 that the spin-orbit coupling of the Wiggle Well does not rely upon the
presence of an interface. Rather, it is an intrinsic property of a bulk system
with Ge concentration oscillations. In this sense, the spin-orbit coupling
investigated here is more akin to the true Dresselhaus spin-orbit coupling of
III-V semiconductors than the Dresselhaus-like spin-orbit coupling of
conventional Si/SiGe quantum wells brought about by interfaces. Indeed, the
only requirement for linear-$\mathbf{k}_{\parallel}$ Dresselhaus spin-orbit
coupling in a Wiggle Well with an appropriate $\lambda$ is confinement in the
growth direction (even symmetric confinement), to allow for the formation of
subbands. For simplicity in the remainder of this work, we simply refer to
this form of spin-orbit coupling as Dresselhaus.

The rest of this article is organized as follows. In Sec. II we describe our
model used to study the quantum well heterostructure. Sec. III then presents
our numerical results for the spin-orbit coefficients. This also includes the
calculation of the EDSR Rabi frequency and studies the impact of alloy
disorder on the spin-orbit coefficients. In Sec. IV we provide an extensive
explanation of the mechanism behind the spin-orbit coupling enhancement.
Finally, we conclude in Sec. V.

##  II Model

In this section, we outline the model used to study our Si/SiGe
heterostructure along with the methods used to calculate the spin-orbit
coefficients. In Sec. II.1, we describe the tight binding model used to model
generic SiGe alloy systems. Next, in Sec. II.2 we employ a virtual crystal
approximation to impart translation invariance in the plane of the quantum
well, allowing us to reduce the problem to an effective 1D Hamiltonian
parameterized by in-plane momentum $\mathbf{k}_{\parallel}$. In Sec. II.3, we
expand the model around $\mathbf{k}_{\parallel}=0$ to separate out the
Hamiltonian components that give rise to Rashba and Dresselhaus spin-orbit
coupling, respectively, and we explain the important differences between the
two components. Finally, in Sec. II.4, we transform the Hamiltonian into the
subband basis, which allows us to obtain expressions for the Rashba and
Dresselhaus spin-orbit coefficients in each subband.

###  II.1 Model of SiGe alloys

To study the spin-orbit physics of our system we use the empirical tight
binding method [26], where the electronic wavefunction is written as a linear
combination of atomic orbitals,

$\ket{\psi}=\sum_{n,j,\nu,\sigma}\ket{nj\nu\sigma}\psi_{nj\nu\sigma}.$ (2)

Here,
$\bra{\boldsymbol{r}}\ket{nj\nu\sigma}=\phi_{\nu}\left(\boldsymbol{r}-\boldsymbol{R}_{n,j}\right)\ket{\sigma}$
is an atomic orbital centered at position $\boldsymbol{R}_{n,j}$,
corresponding to atom $j$ of atomic layer $n$ along the growth direction
$[001]$, $\nu$ is a spatial orbital index, and $\ket{\sigma}$ is a two
component spinor with $\sigma=\uparrow,\downarrow$ indicating the spin of the
orbital. We use an sp3d5s* basis set with 20 orbitals per atom, onsite spin-
orbit coupling, nearest-neighbor hopping, and strain. Note that nearest-
neighbor sp3d5s* tight binding models are well established for accurately
describing the electronic structure of semiconductor materials over a wide
energy range [27]. Explicitly, $\nu$ is a spatial orbital index from the set
including $s$, $s^{*}$, $p_{i}$ ($i=x,y,z$), and $d_{i}$
($i=xy,yz,zx,z^{2},x^{2}-y^{2}$) orbitals, which are meant to model the outer-
shell orbitals of individual Si and Ge atoms that participate in chemical
bonding. Additionally, these orbitals possesses certain spatial symmetries
that, combined with the diamond crystal structure of the SiGe alloy, dictate
the forms of the nearest neighbor couplings, as first explained in the work of
Slater and Koster [26]. The free parameters of the tight binding model
(including onsite orbital energies, nearest-neighbor hopping energies, strain
parameters, etc.) are then chosen such that the band structure of the system
agree as well as possible with experimental and/or ab-initio data. In this
work, we use the tight binding model and parameters of Ref. [28], which allows
for the modeling of strained, random SiGe alloys with any Ge concentration
profile.

The Hamiltonian of an arbitrary SiGe alloy takes the form,

$H^{mi,nj}_{\mu\sigma,\nu\sigma}=\delta_{mi}^{nj}\left[\delta_{\mu\sigma}^{\nu\sigma^{\prime}}\left(\varepsilon_{\nu}^{(nj)}+V_{n}\right)+\delta_{\sigma}^{\sigma^{\prime}}C_{\mu\nu}^{(nj)}+S_{\mu\sigma,\nu\sigma^{\prime}}^{(nj)}\right]\\\
+\delta_{\sigma}^{\sigma^{\prime}}\left(\delta_{m}^{n+1}T^{(n)}_{i\mu,j\nu}+\delta_{m}^{n-1}T^{(m)\dagger}_{i\mu,j\nu}\right),$
(3)

where
$H^{mi,nj}_{\mu\sigma,\nu\sigma^{\prime}}=\matrixelement{mi\mu\sigma}{H}{nj\nu\sigma^{\prime}}$
and $\delta$ equals $1$ if its subscripts match its superscripts and $0$
otherwise. The first line in Eq. (3) contains intra-atomic terms, where
$\varepsilon_{\nu}^{(nj)}$ is the onsite energy of orbital $\nu$, for atom $j$
in atomic layer $n$, $V_{n}$ is the potential energy due to the vertical
electric field, $C_{\mu\nu}^{(nj)}$ accounts for onsite energy shifts and
couplings caused by strain, and $S_{\mu\sigma,\nu\sigma^{\prime}}^{(nj)}$
accounts for spin-orbit coupling. The matrix $C_{\mu\nu}^{(nj)}$ is determined
by the deformation of the lattice due to strain as detailed in Ref. [28] and
arises from changes in the onsite potential of the atom due to the
displacement of its neighbors. In addition, spin-orbit coupling $S^{(nj)}$ is
an intra-atomic coupling between $p$ orbitals [29] and is the only term in Eq.
(3) that does not conserve spin $\sigma$. (See Appendix B for the explicit
form of $S^{(nj)}$.) Note that the superscripts $(nj)$, which index the atoms,
are needed here because the intra-atomic terms depend on whether an atom is Si
or Ge, as well as the local strain environment. The second line in Eq. (3)
contains inter-atomic terms describing the hopping between atoms on adjacent
atomic layers, where $T^{(n)}$ in the hopping matrix from atomic layer $n$ to
atomic layer $n+1$. Nearly all elements of $T^{(n)}$ are zero, with non-zero
hoppings occurring only between nearest neighbor atoms. A non-zero hopping
matrix element $T^{(n)}_{i\mu,j\nu}$ then depends on three things: (1) the
orbital indices $\mu$ and $\nu$, (2) the types of atoms involved, and (3) the
direction and magnitude of the vector $\mathbf{R}_{n+1,i}-\mathbf{R}_{n,j}$
connecting the atoms. We then use the Slater-Koster table in Ref. [26] along
with the parameters of Ref. [28] to calculate $T^{(n)}_{i\mu,j\nu}$. We note
that strain affects the hopping elements by altering the direction and length
of the nearest neighbor vectors (i.e., the crystalline bonds) [30, 28].

In this work, we let the Ge concentration vary between layers, as shown in
Figs. 1(b) and 1(c), but assume it to be uniform within a given layer. Note
that the large difference in Ge concentration between the barrier and well
regions results in a large conduction band offset that traps electrons inside
the quantum well. This occurs naturally in the tight binding model of Eq. (3)
because $\varepsilon^{(nj)},C^{(nj)}$, and $T^{(n)}_{ij}$ are different for Si
and Ge atoms. Finally, we point the reader to Appendix A for a description of
the lattice constant dependence on strain.

###  II.2 Virtual crystal approximation and pseudospin transformation

While the Hamiltonian in Eq. (3) provides an accurate description of SiGe
alloys, it lacks translation invariance when alloy disorder is present. This
makes the model computationally expensive to solve, and it obscures the
physics of the spin-orbit enhancement coming from the averaged effects of the
inhomogeneous Ge concentration profile. We therefore employ a virtual crystal
approximation where the Hamiltonian matrix elements are replaced by their
value averaged over all alloy realizations. Specifically, we define a virtual
crystal Hamiltonian $H_{\text{VC}}$ with elements
$(H_{\text{VC}})^{mi,nj}_{\mu\sigma,\nu\sigma^{\prime}}=\left<H^{mi,nj}_{\mu\sigma,\nu\sigma^{\prime}}\right>$,
where $\left<\dots\right>$ indicates an average over all possible alloy
realizations. The Hamiltonian is then translation invariant within the plane
of the quantum well. In addition, it is useful to move beyond the original
orbital basis, where the spin is well-defined, to a pseudospin basis defined
by

$\ket{nj\bar{\nu}\bar{\sigma}}=\sum_{\nu\sigma}\ket{nj\nu\sigma}U^{(n)}_{\nu\sigma,\bar{\nu}\bar{\sigma}}$
(4)

where $\bar{\nu}$ are the hybridized orbital states,
$\bar{\sigma}=\Uparrow,\Downarrow$ are the pseudospins, and $U^{(n)}$ is the
transformation matrix for layer $n$. Full details of this basis transformation
can be found in Appendix B. Here, we mention three important features of Eq.
(4). First, the basis transformation diagonalizes the onsite spin-orbit
coupling, making the spin-orbit physics more transparent. Second, the
pseudospin states represent linear combinations of orbitals including both
spin $\uparrow$ and spin $\downarrow$. Third, the transformation matrix
$U^{(n)}$ can be shown to satisfy

$U^{(n)}=\begin{cases}U^{(0)},&n\in\mathbb{Z}_{\text{even}},\\\
U^{(1)},&n\in\mathbb{Z}_{\text{odd}}.\end{cases}$ (5)

This alternating structure for the transformation matrix is crucial to the
results that follow, and results from the presence of two sublattices in the
diamond crystal structure of Si, as shown in Fig. 2(a).

Making use of the virtual crystal approximation, we now convert our 3D
Hamiltonian into an effective 1D Hamiltonian. To begin, we note that the
virtual crystal Hamiltonian takes the simplified form,

$\left(H_{\text{VC}}\right)^{mi,nj}_{\bar{\mu}\bar{\sigma},\bar{\nu}\bar{\sigma}^{\prime}}=\delta_{mi\bar{\sigma}}^{nj\bar{\sigma}^{\prime}}\left[\delta_{\bar{\mu}}^{\bar{\nu}}\left(\bar{\varepsilon}_{\bar{\nu}}^{(n)}+V_{n}\right)+\bar{C}_{\bar{\mu}\bar{\nu}}^{(n)}\right]\\\
+\delta_{m}^{n+1}\bar{T}^{(n)}_{i\bar{\mu}\bar{\sigma},j\bar{\nu}\bar{\sigma}^{\prime}}+\delta_{m}^{n-1}\bar{T}^{(m)\dagger}_{i\bar{\mu}\bar{\sigma},j\bar{\nu}\bar{\sigma}^{\prime}},$
(6)

where
$\left(H_{\text{VC}}\right)^{mi,nj}_{\bar{\mu}\bar{\sigma},\bar{\nu}\bar{\sigma}^{\prime}}=\matrixelement{mi\bar{\mu}\bar{\sigma}}{H_{\text{VC}}}{nj\bar{\nu}\bar{\sigma}^{\prime}}$,
$\bar{\varepsilon}_{\bar{\nu}}^{(n)}$ is the onsite energy of pseudospin
orbital $\bar{\nu}$, and $\bar{C}^{(n)}$ and $\bar{T}^{(n)}$ are the onsite
strain and hopping matrices, respectively, transformed into the pseudospin
basis and averaged over alloy realizations. Note that
$\bar{\varepsilon}^{(n)}$ includes contributions from diagonalizing the onsite
spin-orbit coupling. Importantly, $\bar{\varepsilon}^{(n)}$ and
$\bar{C}^{(n)}$ maintain a dependence on the layer index $n$ due to the non-
uniform Ge concentration profile. In contrast, the dependence of intra-atomic
terms on the intra-layer atom index $j$ has vanished due to the virtual
crystal approximation. Moreover, the translation invariance of the virtual
crystal approximation implies that hopping matrix elements between any two
layers only depend upon the relative position of atoms, i.e.,
$\bar{T}^{(n)}_{i\bar{\mu}\bar{\sigma},j\bar{\nu}\bar{\sigma}^{\prime}}=\bar{T}^{(n)}_{\bar{\mu}\bar{\sigma},\bar{\nu}\bar{\sigma}^{\prime}}(\boldsymbol{R}_{n+1,i}-\boldsymbol{R}_{n,j})$.
We therefore introduce in-plane momentum $\boldsymbol{k}_{\parallel}$ as a
good quantum number and Fourier transform our Hamiltonian. To do so, we define
the basis state

$\ket{\boldsymbol{k}_{\parallel}n\bar{\nu}\bar{\sigma}}=\frac{1}{\sqrt{N_{\parallel}}}\sum_{j}e^{i\boldsymbol{k_{\parallel}}\cdot\boldsymbol{R}_{nj}}\ket{nj\bar{\nu}\bar{\sigma}},$
(7)

where $\boldsymbol{k}_{\parallel}=(k_{x},k_{y})$ and $N_{\parallel}$ is the
number of atoms within each layer. The Hamiltonian has matrix elements

$\widetilde{H}^{mn}_{\bar{\mu}\bar{\sigma},\bar{\nu}\bar{\sigma}^{\prime}}(\boldsymbol{k}_{\parallel})=\delta_{m\bar{\sigma}}^{n\bar{\sigma}^{\prime}}\left[\delta_{\bar{\mu}}^{\bar{\nu}}\left(\bar{\varepsilon}_{\bar{\nu}}^{(n)}+V_{n}\right)+\bar{C}_{\bar{\mu}\bar{\nu}}^{(n)}\right]\\\
+\delta_{m}^{n+1}\widetilde{T}_{\bar{\mu}\bar{\sigma},\bar{\nu}\bar{\sigma}^{\prime}}^{(n)}\left(\boldsymbol{k}_{\parallel}\right)+\delta_{m}^{n-1}\widetilde{T}_{\bar{\mu}\bar{\sigma},\bar{\nu}\bar{\sigma}^{\prime}}^{(m)\dagger}\left(\boldsymbol{k}_{\parallel}\right),$
(8)

where
$\widetilde{H}^{mn}_{\bar{\mu}\bar{\sigma},\bar{\nu}\bar{\sigma}^{\prime}}(\boldsymbol{k}_{\parallel})=\matrixelement{\boldsymbol{k}_{\parallel}m\bar{\mu}\bar{\sigma}}{H_{\text{VC}}}{\boldsymbol{k}_{\parallel}n\bar{\nu}\bar{\sigma}^{\prime}}$,
and
$\widetilde{T}_{\bar{\mu}\bar{\sigma},\bar{\nu}\bar{\sigma}^{\prime}}^{(n)}\left(\boldsymbol{k}_{\parallel}\right)$
is the Fourier transformed hopping matrix given by

$\widetilde{T}_{\bar{\mu}\bar{\sigma},\bar{\nu}\bar{\sigma}^{\prime}}^{(n)}\left(\boldsymbol{k}_{\parallel}\right)=\sum_{l=1}^{2}e^{-i\boldsymbol{k_{\parallel}}\cdot\boldsymbol{r}_{l}^{(n)}}\bar{T}_{\bar{\mu}\bar{\sigma},\bar{\nu}\bar{\sigma}^{\prime}}^{(n)}(\boldsymbol{r}_{l}^{(n)}),$
(9)

where $\boldsymbol{r}_{l}^{(n)}$ is a nearest neighbor vector from a reference
atom in layer $n$ to one of its nearest neighbors in layer $n+1$. For a
diamond lattice, each atom has only has two such bonds, as indicated in Fig.
2(a). Note that the Hamiltonian matrix elements vanish between states with
different momenta due to translational invariance. Hence we obtain an
effective 1D Hamiltonian as a function of $\boldsymbol{k}_{\parallel}$.

Figure 2: (a) Diamond crystal structure of silicon. The dashed lines outline
the conventional unit cell of the face centered cubic lattice. Both red and
blue atoms are silicon but belong to different sublattices. Notice that the
vectors connecting an atom to its four nearest neighbors are fundamentally
different for the red and blue atoms, giving rise to the alternating hopping
structure shown in (b). (b) Effective 1D tight-binding chain, with hopping
matrix terms alternating between
$\widetilde{T}_{+}^{(n)}(\boldsymbol{k}_{\parallel})$ and
$\widetilde{T}_{-}^{(n)}(\boldsymbol{k}_{\parallel})$. Note that each site has
20 orbitals, and only the forward hopping terms are shown. Onsite and backward
hopping terms are not shown. For a SiGe alloy in the virtual crystal
approximation, the two-sublattice structure is retained, but the atoms are
replaced by virtual atoms, with averaged properties consistent with the Ge
concentration of a given layer.

An important feature of the Fourier transformed hopping matrix
$\widetilde{T}^{(n)}(\boldsymbol{k}_{\parallel})$ is that it depends on the
layer index $n$ for two reasons. First, the inhomogeneous Ge concentration
along the growth axis causes the hopping parameters to change slightly from
layer to layer. Second, and more importantly, the diamond crystal structure is
composed of two interleaving face-centered cubic sublattices which each
contribute an inequivalent atom to the primitive unit cell. This is
illustrated in Fig. 2(a) where the atoms belonging to the two sublattices are
colored red and blue, respectively. Indeed, the atoms for
$n\in\mathbb{Z}_{\text{even}}$ and $n\in\mathbb{Z}_{\text{odd}}$ belong to
sublattice 1 and 2, respectively, and have different nearest neighbor vectors.
It is therefore useful to define

$\widetilde{T}^{(n)}\left(\boldsymbol{k}_{\parallel}\right)=\begin{cases}\widetilde{T}_{+}^{(n)}\left(\boldsymbol{k}_{\parallel}\right)&n\in\mathbb{Z}_{\text{even}}\\\
\widetilde{T}_{-}^{(n)}\left(\boldsymbol{k}_{\parallel}\right)&n\in\mathbb{Z}_{\text{odd}}\end{cases}$
(10)

as the hopping matrices for the two sublattices. We stress that the dependence
of $\widetilde{T}_{+}^{(n)}(\boldsymbol{k}_{\parallel})$ and
$\widetilde{T}_{-}^{(n)}(\boldsymbol{k}_{\parallel})$ on the layer index $n$
is due to the inhomogeneous Ge concentration profile, and that
$\widetilde{T}^{(n)}_{+}(\boldsymbol{k}_{\parallel})$ and
$\widetilde{T}^{(n)}_{-}(\boldsymbol{k}_{\parallel})$ differ due to the
diamond crystal structure having two inequivalent atoms in its primitive unit
cell. We can therefore visualize the system, for any given
$\boldsymbol{k}_{\parallel}$, as a 1D, multi-orbital tight binding chain, as
shown in Fig. 2(b), where the hopping terms alternate in successive layers.

###  II.3 Expansion around $\mathbf{k}_{\parallel}=0$

Our goal is to understand the spin-orbit physics of low-energy conduction band
states near the Fermi level. In strained Si/SiGe quantum wells, these derive
from the two degenerate valleys near the Z-point of the strained Brillouin
zone [31]. Therefore, the low-energy states have small
$|\mathbf{k}_{\parallel}|$, and we can understand the spin-orbit physics by
expanding the Fourier transformed hopping matrices
$\widetilde{T}_{\pm}^{(n)}\left(\boldsymbol{k}_{\parallel}\right)$ to linear
order. We find that

$\widetilde{T}_{\pm}^{(n)}(\mathbf{k}_{\parallel})=\widetilde{T}_{0}^{(n)}+\widetilde{T}_{R}^{(n)}(\mathbf{k}_{\parallel})\pm\widetilde{T}_{D}^{(n)}(\mathbf{k}_{\parallel})+\mathcal{O}(\mathbf{k}_{\parallel}^{2}),$
(11)

where $\widetilde{T}_{0}^{(n)}$ is the hopping matrix for
$\mathbf{k}_{\parallel}=0$, and $\widetilde{T}_{R}^{(n)}$ and
$\widetilde{T}_{D}^{(n)}$ contain the linear $\mathbf{k}_{\parallel}$
corrections. These hopping matrix components are found to be

$\displaystyle\widetilde{T}_{0}^{(n)}$
$\displaystyle=\Omega^{(n)}\bar{\sigma}_{0},$ (12)
$\displaystyle\widetilde{T}_{R}^{(n)}(\mathbf{k}_{\parallel})$
$\displaystyle=\Phi^{(n)}\left(k_{y}\bar{\sigma}_{x}-k_{x}\bar{\sigma}_{y}\right),$
(13) $\displaystyle\widetilde{T}_{D}^{(n)}(\mathbf{k}_{\parallel})$
$\displaystyle=\Phi^{(n)}\left(k_{x}\bar{\sigma}_{x}-k_{y}\bar{\sigma}_{y}\right),$
(14)

where $\Omega^{(n)}$ and $\Phi^{(n)}$ are real-valued $10\times 10$ matrices,
and $\bar{\sigma}_{j}$ are the Pauli matrices acting on pseudospin space, with
$j=0,x,y,z$.

There are several features to remark on in Eqs. (12)-(14). First, the
momentum-spin structure of $\widetilde{T}_{R}^{(n)}$ and
$\widetilde{T}_{D}^{(n)}$ have the familiar forms of Rasbha and Dresselhaus
spin-orbit coupling; hence, we apply the subscript labels $R$ and $D$. Second,
the hopping matrices for the two sublattices, $\widetilde{T}_{+}^{(n)}$ and
$\widetilde{T}_{-}^{(n)}$, differ by the sign in front of
$\widetilde{T}_{D}^{(n)}$, while the sign of $\widetilde{T}_{R}^{(n)}$ is
sublattice independent. As noted above, this alternating sign is a consequence
of the two sublattices of the diamond crystal being inequivalent. As we shall
show in Sec. IV.3, this alternating hopping structure of
$\widetilde{T}_{D}^{(n)}$ is key to explaining the mechanism behind the
enhanced spin-orbit coupling. In Appendix C, we also provide a symmetry
argument for why the system has this particular alternating hopping structure
for the Rashba and Dresselhaus terms. Third, for the special case of
$\mathbf{k}_{\parallel}=0$, the two sublattice become equivalent (i.e., there
is no even/odd structure due to the vanishing of $\widetilde{T}_{R}^{(n)}$ and
$\widetilde{T}_{D}^{(n)}$ at $\mathbf{k}_{\parallel}=0$). Furthermore, the
pseudospin sectors are uncoupled and equivalent at $\mathbf{k}_{\parallel}=0$,
which implies that all eigenstates of our Hamiltonian are doubly degenerate at
$\mathbf{k}_{\parallel}=0$ as is expected since the system has time-reversal
symmetry. Fourth, the dependence of $\Omega^{(n)}$ and $\Phi^{(n)}$ on the
layer index $n$ arises only from the inhomogeneous Ge concentration profile.
However, note that neither matrix vanishes if the Ge concentration is uniform.
For interested readers, we provide expressions for $\Omega^{(n)}$ and
$\Phi^{(n)}$ in Appendix D, for the case of a pure Si structure. Finally, we
note that the diagonal elements of $\Phi^{(n)}$ all vanish.

To linear order in $\mathbf{k}_{\parallel}$, the Hamiltonian then takes the
compact form

$H=H_{0}^{(z)}\bar{\sigma}_{0}+H_{R}^{(z)}\left(k_{y}\bar{\sigma}_{x}-k_{x}\bar{\sigma}_{y}\right)\\\
+H_{D}^{(z)}\left(k_{x}\bar{\sigma}_{x}-k_{y}\bar{\sigma}_{y}\right)+\mathcal{O}(\mathbf{k}_{\parallel}^{2}),$
(15)

where $H_{0}^{(z)}$ (which we call the subband Hamiltonian) describes the
physics at $\mathbf{k}_{\parallel}=0$, and $H_{R}^{(z)}$ and $H_{D}^{(z)}$
describe the linear $\mathbf{k}_{\parallel}$ perturbations arising from
$\widetilde{T}_{R}^{(n)}$ and $\widetilde{T}_{D}^{(n)}$, respectively. These
take the form

$\displaystyle\begin{split}\matrixelement{m\bar{\mu}}{H_{0}^{(z)}}{n\bar{\nu}}=&~{}\delta_{m}^{n}\left[\delta_{\bar{\mu}}^{\bar{\nu}}\left(\bar{\varepsilon}_{\bar{\nu}}^{(n)}+V_{n}\right)+\bar{C}_{\bar{\mu}\bar{\nu}}^{(n)}\right]\\\
&+\delta_{m}^{n+1}\Omega_{\bar{\mu}\bar{\nu}}^{(n)}+\delta_{m}^{n-1}\Omega^{(m)T}_{\bar{\mu}\bar{\nu}}\end{split},$
(16) $\displaystyle\matrixelement{m\bar{\mu}}{H_{R}^{(z)}}{n\bar{\nu}}=$
$\displaystyle\delta_{m}^{n+1}\Phi^{(n)}_{\bar{\mu}\bar{\nu}}+\delta_{m}^{n-1}\Phi^{(m)T}_{\bar{\mu}\bar{\nu}},$
(17) $\displaystyle\matrixelement{m\bar{\mu}}{H_{D}^{(z)}}{n\bar{\nu}}=$
$\displaystyle(-1)^{n}\left(\delta_{m}^{n+1}\Phi^{(n)}_{\bar{\mu}\bar{\nu}}-\delta_{m}^{n-1}\Phi^{(m)T}_{\bar{\mu}\bar{\nu}}\right),$
(18)

where the superscript $z$ indicates that only the orbital degrees of freedom
in the $z$ direction (i.e., the growth direction) are acted upon. Hence the
momentum $\mathbf{k}_{\parallel}$ and pseudospin $\bar{\sigma}$ indices are
both dropped in Eqs. (16)-(18). Also note that the alternating $\pm$ factor in
front of $\widetilde{T}_{D}^{(n)}$ in Eq. (11) is reflected in the $(-1)^{n}$
factor in Eq. (18).

###  II.4 Transformation to the subband basis

The largest term in Hamiltonian (15) (by far) is the subband Hamiltonian
$H_{0}^{(z)}$. The eigenstates of $H_{0}^{(z)}$ are referred to as the orbital
subbands of the quantum well, including two distinct valley states per
subband. The subband and valley states, in turn, serve as a natural basis for
representing the Hamiltonian, since the lateral confinement associated with a
quantum dot barely perturbs this subband designation, although disorder may
cause hybridization of the valley states. It is therefore the properties of
the individual subbands that largely determine the properties of quantum dot
states, including their spin-orbit behavior.

To perform a subband basis transformation, we define $\ket{\varphi_{\ell}}$ as
the $\ell^{\text{th}}$ eigenstate of $H_{0}^{({z})}$ with energy $E_{\ell}$.
(Here, for convenience, we include both subband and valley states in the set
$\\{\ell\\}$.) Generically we can write

$\ket{\varphi_{\ell}}=\sum_{n\bar{\nu}}\ket{n\bar{\nu}}Q_{n\bar{\nu},\ell},$
(19)

where $Q$ is an orthogonal matrix, defined such that

$H_{0}^{({z})}\ket{\varphi_{\ell}}=E_{\ell}\ket{\varphi_{\ell}}$ (20)

for each $\ell$. Using these eigenstates as a basis, the Hamiltonian can then
be expressed as

$H=\bar{\Lambda}\sigma_{0}+\bar{\alpha}\left(k_{y}\bar{\sigma}_{x}-k_{x}\bar{\sigma}_{y}\right)\\\
+\bar{\beta}\left(k_{x}\bar{\sigma}_{x}-k_{y}\bar{\sigma}_{y}\right)+\mathcal{O}(\mathbf{k}_{\parallel}^{2}),$
(21)

where $\bar{\Lambda}$, $\bar{\alpha}$ and $\bar{\beta}$ are real-symmetric
matrices acting in subband-space with elements

$\displaystyle\bar{\Lambda}_{\ell\ell^{\prime}}=$
$\displaystyle\delta_{\ell\ell^{\prime}}E_{\ell},$ (22)
$\displaystyle\bar{\alpha}_{\ell\ell^{\prime}}=$
$\displaystyle\matrixelement{\varphi_{\ell}}{H_{R}^{(z)}}{\varphi_{\ell^{\prime}}},$
(23) $\displaystyle\bar{\beta}_{\ell\ell^{\prime}}=$
$\displaystyle\matrixelement{\varphi_{\ell}}{H_{D}^{(z)}}{\varphi_{\ell^{\prime}}}.$
(24)

The matrix elements $\bar{\alpha}_{\ell\ell^{\prime}}$ and
$\bar{\beta}_{\ell\ell^{\prime}}$ are referred to as the Rashba and
Dresselhaus spin-orbit coupling coefficients, respectively [16]. The diagonal
elements are of particular importance since they determine the linear
dispersion of a given subband near $\mathbf{k}_{\parallel}=0$. Indeed, the
diagonal elements $\bar{\alpha}_{\ell\ell}$ and $\bar{\beta}_{\ell\ell}$
themselves are often referred to in the literature as the Rasbha and
Dresselhaus spin-orbit coupling coefficients, respectively, and are typically
denoted simply as $\alpha$ and $\beta$. Furthermore, we focus on the diagonal
elements of the ground ($\ell=0$) and excited ($\ell=1$) valley states
corresponding to the lowest orbital subband, which we henceforth refer to as
simply the ground and excited valley states. These represent the lowest-energy
conduction subbands, which are nearly degenerate due to the wide separation of
the two degenerate $z$ valleys within the Brillouin zone of Si [32], thus
playing a dominating role in the physics of Si spin qubits. In some cases, we
may also be interested in the Dresselhaus coupling between the ground and
excited valleys, since the valley states are much closer in energy than the
orbitally excited subbands. Finally, we mention that confinement in the growth
direction is a crucial ingredient for obtaining nonzero values of
$\bar{\alpha}_{\ell\ell}$ and $\bar{\beta}_{\ell\ell}$. While this latter fact
is not obvious from the structure of the Hamiltonian, it can be shown to be
true, using the fact that $\Phi^{(n)}$ has vanishing diagonal elements and the
structure of the $\Omega^{(n)}$ and $\Phi^{(n)}$ matrices described in
Appendix D.

###  II.5 Summary of calculation procedure

To conclude this section, we present a brief summary of the procedure used in
a typical calculation, like those reported in Sec. III. First, we specify a Ge
concentration profile as a function of layer index $n$. Second, we construct
the subband Hamiltonian $H_{0}^{(z)}$ using Eq. (16). Third, we diagonalize
the subband Hamiltonian to obtain a set of eigenstates
$\left\\{\ket{\varphi_{\ell}}\right\\}$. Finally, we construct the
$H_{R}^{(z)}$ and $H_{D}^{(z)}$ matrices in Eqs. (17) and (18) and calculate
the matrix elements in Eqs. (23) and (24), which yields the spin-orbit
coupling coefficients.

##  III Numerical results

In this section, we present our numerical results. We first present in Sec.
III.1 spin-orbit coupling results for a “conventional” Si/SiGe quantum well
system without Ge oscillations included in the quantum well region. These
results serve as a baseline for comparison. Next, we provide spin-orbit
coupling results in Sec. III.2 for the Si/SiGe quantum well system with Ge
oscillations included in the quantum well region, namely a Wiggle Well. In
this case, we observe significant enhancement of the Dresselhaus spin-orbit
coupling for appropriate Ge concentration oscillation wavelengths $\lambda$.
In Sec. III.3, we show that the enhanced spin-orbit coupling allows for fast
Rabi oscillations using EDSR. Finally, we study in Sec. III.4 the impact of
alloy disorder on the spin-orbit coupling.

###  III.1 Spin-orbit coupling in Si/SiGe quantum wells without Ge
concentration oscillations

We first calculate spin-orbit coupling for the conventional Si/SiGe quantum
well shown in Fig. 1(c). We assume barrier regions with a uniform Ge
concentration of $n_{\text{Ge,bar}}=30\%$ and a quantum well width of
$L_{z}\approx 20~{}\text{nm}$, consisting of an even number of atomic layers;
the latter value was chosen to ensure that the wavefunctions have negligible
weight at the bottom barrier except in the limit of very weak electric fields.
In addition, the interfaces are given a nonzero width of
$L_{\text{int}}\approx 0.95~{}\text{nm}$ ($7$ atomic layers), in which the Ge
concentration linearly interpolates between values appropriate for the barrier
and well regions. Such finite-width interfaces occur in realistic devices, and
are known to significantly impact important properties of the quantum well
such as the valley splitting [33].

Figure 3: Diagonal Dresselhaus $\bar{\beta}_{\ell\ell}$ (a) and Rashba
$\bar{\alpha}_{\ell\ell}$ (b) spin-orbit coupling coefficients as a function
of vertical electric field $F_{z}$, for the ground (solid blue,
$\bar{\beta}_{00}$ and $\bar{\alpha}_{00}$) and excited (dashed red,
$\bar{\beta}_{11}$ and $\bar{\alpha}_{11}$) valleys of the “conventional”
Si/SiGe quantum well shown in Fig. 1(c). Notice that the Dresselhaus
coefficients $\bar{\beta}_{\ell\ell}$ are $\sim 7$ times larger in magnitude
than the Rashba coefficients $\bar{\alpha}_{\ell\ell}$.

The diagonal Dresselhaus and Rashba spin-orbit coupling coefficients are
calculated for the ground (solid blue, $\bar{\beta}_{00}$ and
$\bar{\alpha}_{00}$) and excited (dashed red, $\bar{\beta}_{11}$ and
$\bar{\alpha}_{11}$) valley states, and are plotted in Figs. 3(a) and 3(b) as
a function of vertical electric field $F_{z}$. At zero field, $F_{z}=0$, the
diagonal spin-orbit coupling coefficients in Fig. 3 all vanish, as consistent
with the system being inversion symmetric [20, 21]. By turning on the electric
field we break the structural inversion symmetry, and the resulting spin-orbit
coefficients vary linearly over the entire field range considered here. Note
again that Dresselhaus spin-orbit coupling in SiGe requires the presence of a
broken structural inversion symmetry [24, 20, 21, 25], in contrast with GaAs,
which requires a broken bulk inversion symmetry [23]. As consistent with
previous studies [34, 20, 35, 36], the Dresselhaus coefficients of the ground
and excited valleys are found to be approximately opposite in sign, with the
Dresselhaus coefficient of the excited valley being slightly smaller in
magnitude. Additionally, the diagonal Rashba matrix elements are seen to be
much smaller in magnitude than the Dresselhaus elements. Here, we find that
$|\bar{\beta}_{\ell\ell}|\approx 7|\bar{\alpha}_{\ell\ell}|$ for both low-
energy valleys ($\ell=0,1$). In addition, the magnitudes of the spin-orbit
elements are found to quantitatively agree with Ref. [20], in the large
electric field regime. (The system studied in Ref. [20] assumed a narrower
quantum well, resulting in different behavior at low electric fields.) We
note, however, that the Rashba coefficients of the ground and excited valley
states were found to have opposite signs in Ref. [20], while we find them to
have the same sign here. We have numerically confirmed our results for the
diagonal spin-orbit matrix elements using the computational scheme of Ref.
[20], which is unrelated to our scheme, summarized in Eqs. (23) and (24).

###  III.2 Spin-orbit coupling in Si/SiGe quantum wells with Ge concentration
oscillations

We now calculate spin-orbit coefficients for the Wiggle Well geometry shown in
Fig. 1(b), with a sinusoidally varying Ge concentration ranging from
$n_{\text{Ge}}=0\%$ to a maximum amplitude of $n_{\text{Ge}}=10\%$, and an
oscillation wavelength of $\lambda$. These parameters were chosen to match
those of an experimental device reported in Ref. [18]. Here, the system
parameters, $n_{\text{Ge,bar}}$, $L_{z}$, and $L_{\text{int}}$ are the same as
before, and we apply an electric field of $F_{z}=10~{}\text{mV/nm}$.

Figure 4: Diagonal Dresselhaus $\bar{\beta}_{\ell\ell}$ (a) and Rashba
$\bar{\alpha}_{\ell\ell}$ (b) spin-orbit matrix elements as a function of Ge
oscillation wavelength $\lambda$ for an electric field of strength
$F_{z}=10~{}\text{mV/nm}$. Coefficients for the ground and excited valleys are
shown in blue (solid, $\bar{\beta}_{00}$ and $\bar{\alpha}_{00}$) and red
(dashed, $\bar{\beta}_{11}$ and $\bar{\alpha}_{11}$), respectively. A wide
resonance centered at $\lambda\approx 1.57~{}\text{nm}$ occurs in (a),
corresponding to a dramatically enhanced Dresselhaus spin-orbit coupling. At
the center of the resonance, $|\bar{\beta}_{\ell\ell}|$ is $\sim 15$ times
larger than the “conventional” Si/SiGe system at the same electric field. (See
Fig. 3(a).) Narrow resonances for the Dresselhaus and Rashba coefficients at
small $\lambda$ values are shown in (c) and (d). The corresponding Ge
concentration profile is shown in Fig. 1(b).

The resulting diagonal Dresselhaus $\bar{\beta}_{\ell\ell}$ and Rashba
$\bar{\alpha}_{\ell\ell}$ spin-orbit coefficients for the ground (solid blue,
$\bar{\beta}_{00}$ and $\bar{\alpha}_{00}$) and excited (dashed red,
$\bar{\beta}_{11}$ and $\bar{\alpha}_{11}$) valley states are shown in Figs.
4(a) and 4(b) as a function of the Ge oscillation wavelength $\lambda$. We
find that the spin-orbit coefficients vary nontrivially with the choice of
oscillation wavelength $\lambda$, exhibiting resonances where the coefficients
are significantly enhanced. In particular, we observe a broad resonance in
Fig. 4(a), centered at $\lambda\approx 1.57~{}\text{nm}$, in which the
Dresselhaus spin-orbit coefficient reaches the values
$|\bar{\beta}_{00}|\approx 365~{}\mu\text{eV}\cdot\text{nm}$ for the ground
valley and $|\bar{\beta}_{11}|\approx 357~{}\mu\text{eV}\cdot\text{nm}$ for
the excited valley. Note that this resonance corresponds to the same
wavelength as the long-period Wiggle Well predicted in Refs. [18, 19] to
enhance the valley splitting. (Note that a slightly different $\lambda$ value
was predicted in Refs. [18, 19], due to the valley minima residing at slightly
different momenta $\pm k_{0}$ compared to our model.) We comment on the
relation between the spin-orbit enhancement and predicted valley-splitting
enhancement in Sec. IV.3. Similar to the conventional Si/SiGe quantum well
shown in Fig. 4, the ground and excited valleys have approximately opposite
Dresselhaus coefficients for the range of $\lambda$ values considered here,
except at very small $\lambda$. At the center of the resonance, however, the
diagonal Dresselhaus coefficients $\bar{\beta}_{\ell\ell}$ for the Wiggle Well
are $\sim 15$ times larger than those of the conventional system at the same
electric field strength. (See Fig. 3(a) at $F_{z}=10~{}\text{mV/nm}$.)
Moreover, the resonance appears very broad, with a full-width-at-half-maximum
of $\Delta\lambda=0.55~{}\text{nm}$. In contrast to the diagonal Dresselhaus
coefficients $\bar{\beta}_{\ell\ell}$, the diagonal Rashba coefficients
$\bar{\alpha}_{\ell\ell}$ show no resonance features except at small
wavelengths. Indeed, for $\lambda\gtrsim 0.5~{}\text{nm}$, the Rashba
coefficients are essentially independent of $\lambda$ with
$\bar{\alpha}_{00}\approx-3.7~{}\mu\text{eV}\cdot\text{nm}$ and
$\bar{\alpha}_{11}\approx-3.3~{}\mu\text{eV}\cdot\text{nm}$ for the ground and
excited valleys, respectively. Comparing this to the results of the
conventional system in Fig. 3(b) for the same electric field, we see that the
Rashba elements are nearly identical in the two cases. Evidently, the Rashba
coefficients are unaffected by the Ge concentrations oscillations for all but
the shortest $\lambda$ values.

Figure 5: Same as Fig. 4 except with an electric field strength of
$F_{z}=2~{}\text{mV/nm}$. Notice that the magnitude of
$\bar{\beta}_{\ell\ell}$ at the center of the resonances are nearly the same
as in Fig. 4, but that the resonance features are narrower than in Fig. 4.

To illustrate the resonance features at small wavelengths, we plot the
Dresselhaus and Rashba coefficients over selected, narrow ranges of $\lambda$
in Figs. 4(c) and 4(d). Here we observe a resonance in the Dresselhaus
coefficient $\bar{\beta}_{\ell\ell}$ centered at $\lambda=a/2\approx
0.27~{}\text{nm}$, where $a$ is the size of the cubic unit cell in Fig. 2(a)
with an amplitude about twice that of the $\lambda=1.57~{}\text{nm}$
resonance. This $\lambda$ value corresponds to having a Ge concentration
profile that alternates on every other atomic layer, essentially transforming
the diamond lattice of Si into a zinc blende lattice. In this case, the system
can be thought of as a III-V semiconductor with the cation and anion
corresponding to different concentrations of Ge. Interestingly, we find that
the ground and excited valleys have the same sign of $\bar{\beta}_{\ell\ell}$
here, in contrast to the resonance centered at $\lambda\approx
1.57~{}\text{nm}$. We also observe a resonance in the Rashba coefficient
$\bar{\alpha}_{\ell\ell}$ near $\lambda=0.46~{}\text{nm}$ as shown in Fig.
4(d). However, note that the magnitude of this resonance peak is still more
than an order of magnitude smaller than the $\bar{\beta}_{\ell\ell}$
resonances. Indeed, the magnitude of this resonance peak is even smaller than
the Dresselhaus coefficients $\bar{\beta}_{\ell\ell}$ shown in Fig. 3(a) for
the conventional Si/SiGe system at the same electric field strength. Finally,
we note that these small-$\lambda$ resonances have a much narrower width than
the $\bar{\beta}_{\ell\ell}$ resonance at $\lambda=1.57~{}\text{nm}$, which
has important consequences for practical applications.

To study the tunability of the spin-orbit physics, we also calculate the spin-
orbit coefficients of a Wiggle Well for a weaker vertical electric field.
These results are shown in Fig. 5, which is the same system as Fig. 4, except
with $F_{z}=2~{}\text{mV/nm}$. The results are qualitatively similar to the
stronger electric field case. However, there exists an important quantitative
similarity and difference between the two cases, which we address in the
following two paragraphs.

The remarkable similarity is that the Dresselhaus coefficients
$\bar{\beta}_{\ell\ell}$ at the center of the resonances are nearly identical
in the two cases. Note that this is true for both the
$\lambda=1.57~{}\text{nm}$ and $\lambda=0.27~{}\text{nm}$ resonances.
Evidently, at the center of the resonances, the electric field plays a minor
role. This is in stark contrast with the Dresselhaus coefficients for the
conventional system, where the diagonal Dresselhaus coefficients are
proportional to $F_{z}$ as shown in Fig. 3(a). This highlights the important
difference first discussed in the introduction between spin-orbit coupling in
conventional Si/SiGe quantum wells and Wiggle Wells, namely that the
conventional system requires the presence of an interface and structural
asymmetry, while the Wiggle Well fundamentally does not. In the latter case,
the spin-orbit coupling is an intrinsic property of the bulk system. Indeed,
we have checked that the Dresselhaus spin-orbit coupling persists in a Wiggle
Well of wavelength $\lambda=1.57~{}\text{nm}$, in the absence of an interface,
by calculating the spin-orbit coefficients for a system without barriers, but
instead immersed in a harmonic potential, $V_{n}=V_{0}z_{n}^{2}$, obtaining
similar results.

The main difference between the two cases is that the resonance features are
narrower in $\lambda$ space for the weak electric field of Fig. 5, compared to
the strong electric field of Fig. 4. This represents an important advantage of
strong electric fields for the Wiggle Well system, since it can be challenging
to grow heterostructures with perfect oscillation periods. We conclude that
stronger electric fields provide more reliable access to the enhanced spin-
orbit coupling provided by Ge concentration oscillations, since the control of
$\lambda$ does not need to be as precise during the growth process. The reason
for narrower resonance features in weaker electric field will be explained in
Sec. IV.3. Finally, we note that stronger electric fields also provide larger
valley splittings [37].

###  III.3 Electric dipole spin resonance

Figure 6: Electric dipole spin resonance (EDSR) Rabi frequency
$\Omega_{\text{Rabi}}$ as a function of magnetic field $B$ and Ge
concentration oscillation wavelength $\lambda$ for a quantum dot with
confinement energy $\hbar\omega=1~{}\text{meV}$ and an in-plane AC electric
field amplitude $F_{x}=10^{-2}~{}\text{mV/nm}$. The vertical electric field is
$F_{z}=10~{}\text{mV/nm}$.

Besides being of interest from a purely scientific standpoint, the presence of
spin-orbit coupling can be exploited to perform gate operations within a
quantum computation context. In particular, electric dipole spin resonance
(EDSR) is a powerful technique to manipulate individual spins through all-
electrical means [14, 15]. Here, we calculate the EDSR Rabi frequency for a
single electron in a quantum dot embedded in a Wiggle Well. As shown in Ref.
[15], applying an AC, in-plane electric field of amplitude $F_{x}$ with
frequency $\omega_{d}$ across a quantum dot with spin-orbit coupling leads to
an effective AC in-plane magnetic field. For valley $\ell$, the magnitude of
this effective AC magnetic field is given by

$B_{\text{eff}}(t)=\frac{2eF_{x}}{\hbar\omega}\left(\frac{\omega_{d}}{\omega}\right)\frac{\sqrt{\bar{\beta}^{2}_{\ell\ell}+\bar{\alpha}^{2}_{\ell\ell}}}{g\mu_{B}}\sin\left(\omega_{d}t\right),$
(25)

where $\hbar\omega$ is the level spacing characteristic of the dot’s harmonic
confinement potential, $g$ is the g-factor, and $\mu_{B}$ is the Bohr
magneton. Note that for
$|\bar{\beta}_{\ell\ell}|\gg|\bar{\alpha}_{\ell\ell}|$, as is the case for
Si/SiGe systems, the effective AC magnetic field is parallel to the in-plane
electric field $F_{x}$. Applying a static, out-of-plane magnetic field $B$,
the effective Hamiltonian for the quantum dot restricted to the orbital ground
state is then
$H_{\text{eff}}(t)=\frac{1}{2}g\mu_{B}\left[B\sigma_{z}+B_{\text{eff}}(t)\sigma_{x}\right]$.
For a system initialized in the spin-$\uparrow$ ground valley state and driven
at the resonance frequency $\omega_{d}=g\mu_{B}B/\hbar$ set by the external
magnetic field, the probability of finding the electron in the
spin-$\downarrow$ state is given by
$P_{\downarrow}(t)=\sin^{2}(\Omega_{\text{Rabi}}t/2)$ [38], where the Rabi
frequency is found to be

$\Omega_{\text{Rabi}}=\frac{eF_{x}g\mu_{B}\sqrt{\bar{\beta}^{2}_{\ell\ell}+\bar{\alpha}^{2}_{\ell\ell}}}{\hbar(\hbar\omega)^{2}}B.$
(26)

The EDSR Rabi frequency $\Omega_{\text{Rabi}}$ of a quantum dot in a Wiggle
Well is plotted in Fig. 6 as a function of magnetic field and Ge oscillation
wavelength $\lambda$ for the realistic parameters of
$\hbar\omega=1~{}\text{meV}$ and $F_{x}=10^{-2}~{}\text{mV/nm}$. We see that
for oscillation wavelengths near the middle of the the spin-orbit resonance at
$\lambda\approx 1.57~{}\text{nm}$, we can obtain Rabi frequencies of $100$’s
of MHz for moderate magnetic field strengths. Indeed, a Rabi frequency of
$\Omega_{\text{Rabi}}=1~{}\text{GHz}$ is achieved at the center of the
resonance for a magnetic field $B\approx 1.5~{}\text{T}$. We therefore
conclude that including Ge concentration oscillations in Si/SiGe quantum wells
enables a dramatic speedup of single qubit gates using EDSR.

###  III.4 Impact of alloy disorder

SiGe is a random alloy, and the resulting alloy disorder is known to have an
effect on quantum well properties such as valley splitting, both for
conventional quantum wells [33] and Wiggle Wells [18, 19]. It is therefore
important to explore the effects of alloy disorder on spin-orbit coupling.
Recall from Sec. II.2 that we have employed a virtual crystal approximation in
our model, which averages over all possible alloy realizations. While
providing tractability to our calculations, this approximation ignores the
fluctuations arising from the random nature of the Ge atom arrangements in the
SiGe alloy. In this section, we explain how Ge concentration fluctuations can
be reintroduced into our model, to explore the effects of alloy disorder.

A full 3D calculation including alloy disorder is computationally expensive
due to the loss of translation invariance, and is beyond the scope of the this
work. We can still, however, include the effects of alloy disorder
approximately within our 1D effective model by allowing for fluctuations in
the Ge concentration in each atomic layer. This is accomplished using the
procedure described in Ref. [33], which can be summarized as follows. First,
we assume a dot of radius
$a_{\text{dot}}=\sqrt{\hbar/m_{\parallel}\omega}=20~{}\text{nm}$ in the plane
of the quantum well, where $m_{\parallel}=0.19m_{e}$ is the in-plane effective
mass and $\hbar\omega=1~{}\text{meV}$ is the orbital excitation energy
characterizing the parabolic confinement of the dot. We then calculate the
effective Ge concentration $n_{\text{Ge},n}^{\text{eff}}$ in layer $n$ of our
disordered system by counting the number of Ge atoms within our dot. Here, the
probability of any given atom in layer $n$ being a Ge atom is
$n_{\text{Ge},n}$, where $n_{\text{Ge},n}$ is the average germanium
concentration throughout the entire layer $n$, and the number of atoms in our
dot is $N_{\text{eff}}=4\pi a_{\text{dot}}^{2}/a^{2}\approx 17100$, where
$a=0.543~{}\text{nm}$ is the cubic lattice constant of Si. The effective Ge
concentration $n_{\text{Ge},n}^{\text{eff}}$ in layer $n$ can then be drawn
from the distribution
$n_{\text{Ge},n}^{\text{eff}}=N_{\text{eff}}^{-1}\text{Binom}(N_{\text{eff}},n_{\text{Ge},n})$,
where $\text{Binom}(n,p)$ is the binomial distribution with $n$ trials and
probability of success $p$. In the limit of $N_{\text{eff}}\rightarrow\infty$,
the resulting, randomized effective Ge concentration
$n_{\text{Ge},n}^{\text{eff}}$ approaches the ideal Ge concentration
$n_{\text{Ge},n}$, but for smaller dots, fluctuations from this ideal limit
become more pronounced. We then calculate the spin-orbit coefficients for our
effective 1D model as in previous sections but with the Ge concentration
profile given by $n_{\text{Ge},n}^{\text{eff}}$ instead of $n_{\text{Ge},n}$.
Note that this method of including alloy disorder in the 1D effective model
was shown in Ref. [33] to yield valley splitting distributions in good
agreement with 3D calculations.

Figure 7: (a) Distribution of the diagonal Dresselhaus coefficients
$\bar{\beta}_{\ell\ell}$ of the Wiggle Well when including alloy disorder as
described in the main text. System parameters are $\lambda=1.57~{}\text{nm}$
and $F_{z}=10~{}\text{mV/nm}$. The data includes $1000$ alloy realizations. We
see that the spin-orbit coupling enhancement is robust against alloy
fluctuations with $86\%$ of alloy realizations having
$|\bar{\beta}_{00}|>200~{}\mu\text{eV}\cdot\text{nm}$. (b) Scatter plot of the
diagonal [$\bar{\beta}_{00}$ (blue) and $\bar{\beta}_{11}$ (red)] and off-
diagonal ($\bar{\beta}_{01}$) spin-orbit coefficients of all alloy
realizations. This result indicates that the main effect of the alloy disorder
is to mix the ground and excited valley states.

We now calculate the diagonal Dresselhaus coefficients
$\bar{\beta}_{\ell\ell}$ for the same Si/SiGe system as Fig. 4, with Ge
concentration oscillations of wavelength $\lambda=1.57~{}\text{nm}$, but now
with Ge concentration fluctuations included. Note that this $\lambda$ value
corresponds to the center of the main resonance peak in Fig. 4(a). The
distribution of the spin-orbit coefficients is plotted in Fig. 7(a) for $1000$
random-alloy realizations, for the ground (blue, $\bar{\beta}_{00}$) and
excited (red, $\bar{\beta}_{11}$) valleys. Unsurprisingly, we see that the
alloy fluctuations affect the spin-orbit coefficients, with the ground valley
coefficient spanning the range
$-377<\bar{\beta}_{00}<355~{}\mu\text{eV}\cdot\text{nm}$. The distributions
are highly peaked, however, near the bare values
$\bar{\beta}_{00}=-365~{}\mu\text{eV}\cdot\text{nm}$ and
$\bar{\beta}_{11}=357~{}\mu\text{eV}\cdot\text{nm}$ of the disorder-free
system (see Fig. 4). Indeed, we find that $86\%$ of the alloy realizations
have $|\bar{\beta}_{00}|>200~{}\mu\text{eV}\cdot\text{nm}$. We therefore
conclude that the spin-orbit enhancement arising from the Ge concentration
oscillations within Wiggle Wells is robust against alloy fluctuations.

We gain further insight into the effects of the alloy fluctuations by studying
the inter-valley spin-orbit coefficient $\bar{\beta}_{01}$. Figure 7(b)
presents a scatter plot showing both the diagonal ($\bar{\beta}_{00}$ and
$\bar{\beta}_{11}$) and off-diagonal ($\bar{\beta}_{01}$) coefficients of all
$1000$ alloy realizations. Interestingly, all realizations yield coefficients
that land in a narrow semi-circular region of parameter space. This is an
indication that the alloy disorder is essentially mixing the original ground
and excited valley states of the system without disorder. Note that
perturbations arising from higher orbital subbands are insignificant, except
for inducing a small width to the distribution. This is not unexpected since
the energy difference between the ground and excited valley states is
$<1~{}\text{meV}$, while the lowest orbital excitation energy is $\gtrsim
20~{}\text{meV}$.

Figure 8: (a) Distribution of the diagonal Dresselhaus coefficients
$\bar{\beta}_{\ell\ell}$ of a conventional Si/SiGe system with a pure Si
quantum well when including alloy disorder as described in the main text. The
electric field is $F_{z}=10~{}\text{mV/nm}$. (b) Same as (a) except that
$n_{\text{Ge},n}=5\%$ throughout the well region. Insets show the
$n_{\text{Ge},n}$ profiles. We see that inclusion of Ge inside of the well
region widens the distribution of $\bar{\beta}_{\ell\ell}$ but does not on
average increase its magnitude.

For completeness, we also calculate the distribution of Dresselhaus spin-orbit
coefficients in conventional Si/SiGe quantum wells that include alloy
fluctuations. These results are shown in Fig. 8(a) for a system with no Ge in
the well region and electric field $F_{z}=10~{}\text{mV/nm}$. Here, the Ge
concentration profile is shown in the inset. Again, we obtain
$\bar{\beta}_{\ell\ell}$ distributions centered at the same values as the
disorder-free system. (See Fig. 3, with $F_{z}=10$ mV/nm.) Here, however, the
spread is significantly narrower than the Wiggle Well results shown in Fig.
7(a), with the full-width-at-half-maximum being $\sim
10~{}\mu\text{eV}\cdot\text{nm}$. This is because no alloy fluctuations occur
in the well region where the majority of the wavefunction resides. We also
consider a system with $n_{\text{Ge}}=5\%$ distributed uniformly throughout
the well region, taking into account the effects of random alloy fluctuations.
[Here, we do not include intentional Ge concentration oscillations; the
resulting Ge concentration profile is shown in the inset of Fig. 8(b).] Note
that this system contains the same amount of Ge in the well region as the
Wiggle Well system shown in Fig. 7. Unsurprisingly, we find that the
distribution of the $\bar{\beta}_{\ell\ell}$ coefficients spreads
considerably, compared to the results in Fig. 8(a), due to the presence of
alloy fluctuations inside of the well. Importantly, however, a uniform Ge
concentration in the well region does not on average increase the spin-orbit
coupling, in contrast to the effect on valley splitting. Indeed, we find that
$\left<|\bar{\beta}_{00}|\right>=26~{}\mu\text{eV}\cdot\text{nm}$ for the data
in Fig. 8(a), while
$\left<|\bar{\beta}_{00}|\right>=22~{}\mu\text{eV}\cdot\text{nm}$ for the data
in Fig. 8(b). This highlights that fact that in order to enhance spin-orbit
coupling, it is not enough to simply include Ge in the well region, but rather
the Ge concentration must oscillate with the appropriate wavelength $\lambda$.
We note that this result is consistent with the experimental observation that
the $g$-factor measured in uniform Si1-xGex alloys is only slightly altered by
changing the Ge concentration [39, 40].

##  IV Mechanism behind the spin-orbit enhancement

Having shown through simulations that the inclusion of Ge concentration
oscillations of an appropriate wavelength can significantly enhance the spin-
orbit coupling, the natural question is what mechanism leads to this
enhancement? In this section, we provide an explanation of this mechanism. To
begin, we describe in Sec. IV.1 a simplified version of our model that allows
for easier understanding of the spin-orbit enhancement mechanism. Next, we
study in Sec. IV.2 the real-space representation of the ground valley
wavefunction in the absence and presence of Ge concentration oscillations. We
find that a real-space picture is inadequate in explaining the spin-orbit
enhancement. We therefore study in Sec. IV.3 the structure of the wavefunction
in momentum space. We show how the combination of the oscillating potential
produced by the Ge concentration oscillations and the selection rules of
Dresselhaus spin-orbit coupling leads to the spin-orbit coupling enhancement.

###  IV.1 Simplified Model

To focus on the essential physics for the spin-orbit coupling enhancement, we
use in this section a model for SiGe alloys that is slightly simplified
compared to the model presented in Sec. II and used in our numerical
calculations in Sec. III. In this model, Ge atoms are assumed to be identical
to Si atoms, except for their orbitals being shifted up in energy by a
constant, i.e.
$\bar{\varepsilon}_{\bar{\nu}}^{(\text{Ge})}=\bar{\varepsilon}_{\bar{\nu}}^{(\text{Si})}+E_{\text{Ge}}$
where $E_{\text{Ge}}=0.8~{}\text{eV}$ is the extra energy of every Ge orbital.
This is meant to capture at the simplest level that inclusion of Ge increases
the energy of the conduction band minima. In particular, the chosen value
produces a band offset of $0.24~{}\text{eV}$ between a pure Si region and a
barrier region with a nominal 30% Ge concentration. Note, however, that the
precise value is not important since we are only using this simplified model
to understand the spin-orbit enhancement mechanism, leaving quantitative
questions to the more accurate model of Sec. II. In addition, we also neglect
the effects of strain, such that $\bar{C}^{(n)}\rightarrow 0$, since they are
not crucial in understanding the spin-orbit enhancement mechanism. With these
simplifications, the addition of Ge is equivalent to adding a term to the
potential energy $V$ of a pure Si system. For simplicity, we define a new
potential $\mathcal{V}_{n}=eF_{z}z_{n}+E_{\text{Ge}}n_{\text{Ge},n}$, that
includes both the electric field $F_{z}$ and the energy shift from the Ge
concentration $n_{\text{Ge},n}$ of the layer, and we let all orbitals energies
take values appropriate for Si:
$\bar{\varepsilon}_{\bar{\nu}}^{(n)}\rightarrow\bar{\varepsilon}_{\bar{\nu}}^{(\text{Si})}$.
Here, $z_{n}=na/4$ is the $z$ coordinate of atomic layer $n$, and
$a=0.543~{}\text{nm}$ is the cubic lattice constant of Si. Importantly, in
this simplified model, the onsite orbital energies and hopping matrices all
lose their dependence on the layer index $n$,
$\left\\{\bar{\varepsilon}_{\bar{\nu}}^{(n)},\widetilde{T}^{(n)}_{0},\widetilde{T}^{(n)}_{R}(\mathbf{k}_{\parallel}),\widetilde{T}^{(n)}_{D}(\mathbf{k}_{\parallel}),\Omega^{(n)},\Phi^{(n)}\right\\}\rightarrow\left\\{\bar{\varepsilon}_{\bar{\nu}}^{(\text{Si})},\widetilde{T}_{0},\widetilde{T}_{R}(\mathbf{k}_{\parallel}),\widetilde{T}_{D}(\mathbf{k}_{\parallel}),\Omega,\Phi\right\\}$.
The Hamiltonian components then take the simplified forms

$\displaystyle\begin{split}\matrixelement{m\bar{\mu}}{H_{0}^{(z)}}{n\bar{\nu}}=&~{}\delta_{m\bar{\mu}}^{n\bar{\nu}}\left(\bar{\varepsilon}_{\bar{\nu}}^{(\text{Si})}+\mathcal{V}_{n}\right)\\\
&+\delta_{m}^{n+1}\Omega_{\bar{\mu}\bar{\nu}}+\delta_{m}^{n-1}\Omega^{T}_{\bar{\mu}\bar{\nu}}\end{split},$
(27) $\displaystyle\matrixelement{m\bar{\mu}}{H_{R}^{(z)}}{n\bar{\nu}}=$
$\displaystyle\delta_{m}^{n+1}\Phi_{\bar{\mu}\bar{\nu}}+\delta_{m}^{n-1}\Phi^{T}_{\bar{\mu}\bar{\nu}},$
(28) $\displaystyle\matrixelement{m\bar{\mu}}{H_{D}^{(z)}}{n\bar{\nu}}=$
$\displaystyle(-1)^{n}\left(\delta_{m}^{n+1}\Phi_{\bar{\mu}\bar{\nu}}-\delta_{m}^{n-1}\Phi^{T}_{\bar{\mu}\bar{\nu}}\right).$
(29)

Note that the numerical values of the $\Omega$ and $\Phi$ matrices used in the
above equations are specified in Appendix D.

###  IV.2 Real space wavefunctions

Let us now study the effects of the Ge concentration oscillations on the
wavefunctions of the subband Hamiltonian $H_{0}^{(z)}$. To do so, we first
calculate the ground valley wavefunction of a conventional Si/SiGe quantum
well that has the Ge concentration profile shown as the blue line in Fig.
9(a). The ground valley wavefunction $|\psi|^{2}$ of the subband Hamiltonian
$H_{0}^{(z)}$ is shown in red in Fig. 9(a), where the state is pushed up
against the barrier-well interface by an electric field of strength
$F_{z}=5~{}\text{mV/nm}$. The wavefunction also exhibits fast oscillations
characteristic of the superposition of the two valley minima [32]. We next
consider the same system, except with Ge concentration oscillations of
wavelength $\lambda=1.62~{}\text{nm}$ included in the well region, as shown by
the blue line in Fig. 9(b). Comparing the two cases, we see that the
wavefunction is suppressed in regions of high Ge concentration, as consistent
with the fact that the conduction band minima (and thus the local potential
energies) are higher in energy in those regions. However, these qualitative
observations do not directly explain the enhancement of Dresselhaus spin-orbit
coupling observed in Figs. 4 and 5.

Figure 9: Comparison of the ground valley wavefunction for a “conventional”
Si/SiGe system with a pure Si well region (a) and a Wiggle Well (b) with
oscillation wavelength $\lambda=1.62~{}\text{nm}$. Blue lines show the Ge
concentration profiles while red filled in curves are the wavefunctions
$|\psi|^{2}$. Note that we sum over orbital indices. The electric field is
$F_{z}=5~{}\text{mV/nm}$ in both (a) and (b). We see in (b) that the
wavefunction is suppressed in regions of high Ge concentration, consistent
with the fact that the conduction band minima are higher in energy in those
regions.

###  IV.3 Momentum space wavefunctions

We gain a clearer understanding of the effects of the Ge concentration
oscillations by studying the momentum space representation of the ground
valley. To begin, let us first investigate the band structure of Si at
$\mathbf{k}_{\parallel}=0$, which is shown as the orange dashed lines in Fig.
10. A central feature of the band structure is the presence of two degenerate
valleys near zero energy that give rise to the ground and excited valley
states in a quantum well. In addition, the valley minima are located only a
short distance of $0.17(2\pi/a)$ away from the Brillouin zone edge at
$k_{z}=\pm 2\pi/a$. Interestingly, we notice that at $k_{z}=\pm 2\pi/a$, there
are no band anti-crossings (only crossings). This is an indication that the
Brillouin zone can be enlarged for $\mathbf{k}_{\parallel}=0$. Indeed, for the
special case of $\mathbf{k}_{\parallel}=0$, the two sublattice hopping
matrices become equal,
$\widetilde{T}_{+}(\mathbf{0})=\widetilde{T}_{-}(\mathbf{0})$, and the
primitive unit cell of the 1D chain shown in Fig. 2(b) reduces from two sites
to one. The Brillouin zone should therefore extend to $k_{z}=\pm 4\pi/a$
instead of $k_{z}=\pm 2\pi/a$. Here, we will refer to the Brillouin zone that
extends to $k_{z}=\pm 4\pi/a$ as the extended zone, while the zone extending
only to $k_{z}=\pm 2\pi/a$ as the conventional zone.

Figure 10:  The band structure of Si for $\boldsymbol{k}_{\parallel}=0$. The
band structure for the conventional Brillouin zone extending to $k_{z}=\pm
2\pi/a$ is shown by dashed orange lines. The presence of band crossings at the
zone edge $k_{z}=\pm 2\pi/a$ indicates that the Brillouin zone can be
enlarged. Indeed, the equivalence of the two sublattices of Si for
$\mathbf{k}_{\parallel}=0$ allows us to extend the Brillouin zone to
$k_{z}=\pm 4\pi/a$. The band structure for the extended zone is shown as solid
black lines. See the main text for details on the green, blue, and red points.

To calculate the band structure in the extended zone, we define a plane wave
basis state with momentum $k_{z}$ as

$\ket{k_{z}\bar{\nu}}=\frac{1}{\sqrt{N_{z}}}\sum_{n}e^{ik_{z}z_{n}}\ket{n\bar{\nu}},$
(30)

where $N_{z}$ is the number of sites in the 1D chain and $-4\pi/a<k_{z}\leq
4\pi/a$, where $k_{z}=8\pi n/(N_{z}a)$. The subband Hamiltonian $H_{0}^{(z)}$
then has matrix elements,

$\begin{split}\matrixelement{k_{z}\bar{\mu}}{H_{0}^{(z)}}{k_{z}^{\prime}\bar{\nu}}=&\delta_{k_{z}}^{k_{z}^{\prime}}\Big{(}\delta_{\bar{\mu}}^{\bar{\nu}}\bar{\varepsilon}_{\bar{\nu}}^{(\text{Si})}+e^{-ik_{z}\frac{a}{4}}\Omega_{\bar{\mu}\bar{\nu}}+e^{ik_{z}\frac{a}{4}}\Omega^{T}_{\bar{\mu}\bar{\nu}}\Big{)}\\\
&+\delta_{\bar{\mu}}^{\bar{\nu}}~{}\widetilde{\mathcal{V}}(k_{z}-k_{z}^{\prime})\end{split},$
(31)

where $\widetilde{\mathcal{V}}(k_{z}-k_{z}^{\prime})$ is the Fourier transform
of the potential $\mathcal{V}$ and is given by
$\widetilde{\mathcal{V}}(q_{z})=N_{z}^{-1}\sum_{n}\exp(-iq_{z}z_{n})\mathcal{V}_{n}$.
Notice that in the absence of the potential
[$\widetilde{\mathcal{V}}(q_{z})\rightarrow 0$], $k_{z}$ is a good quantum
number for $H_{0}^{(z)}$, and Eq. (31) represents a Bloch Hamiltonian. The
spectrum of $H_{0}^{(z)}$ for $\widetilde{\mathcal{V}}(q_{z})=0$ is shown as
the solid black lines in Fig. 10, which we refer to as the extended band
structure. Within the conventional zone, we see that the extended band
structure aligns perfectly with half of the conventional (orange) bands, while
the other half of the bands have been removed from the conventional zone and
instead reside in the regions between $k_{z}=|2\pi/a|$ and $k_{z}=|4\pi/a|$.
Furthermore, we observe that the conventional bands that do not match with the
extended band structure can be made to match by shifting them by a reciprocal
lattice vector $G_{z}=\pm 4\pi/a$ of the lattice containing two sites. For
example, the green point at $k_{z}=-\pi/a$ in Fig. 10 gets shifted to the blue
point at $k_{z}=3\pi/a$, which coincides with a band of the extended band
structure. This is expected since the points $k_{z}$ and $k_{z}+G_{z}$ are
equivalent from the point of view of the conventional zone [41].

We stress that the extended band structure contains more information than the
conventional band structure. Indeed, the extended scheme clarifies which
states can be coupled by a given Fourier component $q_{z}$ of the potential
$\widetilde{\mathcal{V}}(q_{z})$. For example, let us consider what Fourier
component of the potential could couple the states of the conventional band
structure marked by the green and red points at $k_{z}=-\pi/a$ and
$k_{z}=-1.3\pi/a$, respectively, in Fig. 10. Looking at the conventional band
structure, one may initially believe the $\widetilde{\mathcal{V}}(-0.3\pi/a)$
Fourier component could couple the states. Looking at these states in the
extended band structure (red and blue points in Fig. 10), however, we
immediately see that these states are instead coupled by the
$\widetilde{\mathcal{V}}(3.7\pi/a)$ Fourier component. This additional
information will be crucial in understanding the enhanced spin-orbit coupling
mechanism below.

Figure 11: Wavefunction profiles of the ground valley in the plane wave
representation. Panels (a) and (b) correspond to Ge profiles (a) and (b) in
Fig. 9, respectively. Both states show main peaks centered at the valley
minima, $k_{z}\approx\pm 0.83(2\pi/a)$, as shown in Fig. 10. As shown in (b),
however, the Ge concentration oscillations produce wavefunction satellites a
distance $2\pi/\lambda$ away from the main peaks. The wave vector from a main
peak to the opposite outer satellite is $4\pi/a$ for this choice of $\lambda$.
Since Dresselhaus spin-orbit coupling connects locations in reciprocal space
separated by $4\pi/a$, spin-orbit enhancement occurs when a satellite peak is
separated from a main peak by $4\pi/a$.

Returning to our comparison of the systems with and without Ge concentration
oscillations, we now plot the ground valley wavefunctions in the plane wave
representation. These are shown in Figs. 11(a) and 11(b) and correspond to the
real-space wavefunctions shown in Figs. 9(a) and 9(b), respectively. That is,
we plot
$(\sum_{\bar{\nu}}|\bra{k_{z}\bar{\nu}}\ket{\varphi_{\ell}}|^{2})^{-1/2}$ as a
function of $k_{z}$ where $\ket{\varphi_{\ell}}$ is the ground valley
wavefunction. As shown in Fig. 11(a), the ground valley wavefuction of the
conventional Si/SiGe heterostructure consists of two peaks centered at
$k_{z}\approx\pm 0.83(2\pi/a)$. These coincide with the conduction band minima
in Fig. 10(a) as expected. Note that the positions of these minima in the band
structure of Si depends on the exact tight binding parameters used, and
differs slightly from other models in the literature. As shown in Fig. 11(b),
the ground valley wavefuction in the presence of the Ge oscillations still has
its main peaks but also includes surrounding satellite features. The location
of these satellites with respect to the the central peaks is determined by the
wavelength of the Ge oscillations, with a peak-satellite separation of
$2\pi/\lambda$, as shown in Fig. 11(b).

The second step in explaining the enhanced spin-orbit coupling requires a
mechanism for coupling different regions of the Brillouin zone. We first
express the matrix elements of the Rashba $H_{R}^{(z)}$ and Dresselhaus
$H_{D}^{(z)}$ Hamiltonian components in the plane wave basis, giving

$\displaystyle\matrixelement{k_{z}\bar{\mu}}{H_{R}^{(z)}}{k_{z}^{\prime}\bar{\nu}}$
$\displaystyle=\delta_{k_{z}}^{k_{z}^{\prime}}\left(e^{-ik_{z}\frac{a}{4}}\Phi_{\bar{\mu}\bar{\nu}}+e^{ik_{z}\frac{a}{4}}\Phi^{T}_{\bar{\mu}\bar{\nu}}\right),$
(32)
$\displaystyle\matrixelement{k_{z}\bar{\mu}}{H_{D}^{(z)}}{k_{z}^{\prime}\bar{\nu}}$
$\displaystyle=\delta_{|k_{z}-k_{z}^{\prime}|}^{4\pi/a}\Big{(}e^{-ik_{z}\frac{a}{4}}\Phi_{\bar{\mu}\bar{\nu}}-e^{ik_{z}\frac{a}{4}}\Phi^{T}_{\bar{\mu}\bar{\nu}}\Big{)}.$
(33)

The key features to notice here are the selection rules between $k_{z}$ and
$k_{z}^{\prime}$: while $H_{R}^{(z)}$ conserves $k_{z}$, $H_{D}^{(z)}$ couples
states with momenta differing by $4\pi/a$. The latter result is obtained by
Fourier transforming the $(-1)^{n}$ factor in Eq. (29), which itself is a
manifestation of the alternating sign in front of $\widetilde{T}_{D}^{(n)}$ in
Eq. (11). To our knowledge, these selection rules and their relation to Rashba
and Dresselhaus spin-orbit couplings have not been noticed previously,
although we speculate that they could be deduced from group theory methods,
such as the method of invariants [16, 42]. We emphasize, however, that the
extended zone scheme is key to obtaining such results, since in the
conventional zone scheme, $k_{z}$ values differing by $4\pi/a$ are equivalent,
and cannot yield a selection rule like Eq. (33). We also note that Eqs. (32)
and (33) apply to any system with a diamond crystal structure.

The Dresselhaus momentum selection rule in Eq. (33) indicates that the spin-
orbit coupling will be enhanced when wavefunction peaks associated with two
different valleys are separated by $4\pi/a$. In Fig. 11(b), we see that this
can only occur for coupling between a central valley peak and the outer
satellite associated with the opposite valley. Since the central valley peaks
are located at $k_{z}\approx\pm 0.83(2\pi/a)$, we see that the resonance
condition for the oscillation wavelength to enhance the Dresselhaus spin-orbit
is given by $\lambda_{\text{res}}=2.94a=1.59~{}\text{nm}$. Uncoincidentally,
this wavelength value is very close to the middle of the resonance of our
numerical result in Fig. 4(a) where the Dresselhaus spin-orbit coupling is
significantly enhanced. Note that similar considerations apply to the excited
valley.

As noted in Sec. III.2, the resonance wavelength $\lambda_{\text{res}}$ is the
same as the long-period Wiggle Well, which was predicted to enhance the valley
splitting in Refs. [18] and [19]. The understanding developed in the present
work allows us to elucidate the relation between these oscillation-induced
spin-orbit and valley-splitting enhancements. First, we note that the valley
splitting cannot be enhanced in a long wavelength Wiggle Well unless the
sublattice symmetry at $\mathbf{k}_{\parallel}=0$ is broken by alloy disorder
[19]. The valley splitting enhancement mechanism therefore also involves a
two-step process: first, the Ge concentration oscillations cause wavefunction
satellites to form in momentum space, similar to Fig. 11(b); then the
satellite peaks couple opposite valleys through alloy disorder, which breaks
the sublattice symmetry. This is precisely analogous to the spin-orbit
enhancement mechanism described above, with the Dresselhaus coupling mechanism
replaced by alloy disorder.

Our understanding of the spin-orbit enhancement mechanism also explains the
narrower resonance features observed in Fig. 5 for a weak vertical electric
field $F_{z}=2~{}\text{mV/nm}$, as compared to the wider resonance observed in
Fig. 4 for a strong electric field $F_{z}=10~{}\text{mV/nm}$. Essentially, a
weaker electric field produces wavefunctions with narrower features in the
plane-wave representation than the those from a strong electric field. These
thinner features make it more difficult to satisfy the Dresselhaus resonance
condition since the narrower peaks need to be situated more precisely in
$k_{z}$ space.

Finally, we comment that the spin-orbit enhancement mechanism relies
fundamentally on the degeneracy of the two $z$ valleys in the band structure
of strained Si. Such a situation could not occur if, for example, Si was a
direct bandgap semiconductor with a single non-degenerate valley at the
$\Gamma$ point, since the key coupling in Fig. 11(b) occurs between the
central peak of one valley and the outer satellite of the opposite valley.
Interestingly, in this case the valley degeneracy of Si can be considered as
beneficial, while in other scenarios it is often considered to be problematic.
In addition, the spin-orbit enhancement mechanism has similarities to holes in
semiconductors, for which the enhancement arises due to degeneracy at the
valence band edge. In contrast to the valence band case, however, which has
degenerate bands at the same momentum, the degenerate $z$ valleys in Si
require a periodic potential in the form of Ge concentration oscillations to
enhance the coupling strength.

##  V Conclusions

We have shown that the inclusion of periodic Ge concentrations oscillations
within the quantum well region of a Si/SiGe heterostructure leads to enhanced
spin-orbit coupling when the oscillation wavelength $\lambda$ is properly
chosen. Specifically, we find that the Dresselhaus spin-orbit coupling
coefficient is enhanced by over an order of magnitude when $\lambda\approx
1.57~{}\text{nm}$, as shown in Figs. 4 and 5. We have provided a detailed
explanation for this behavior: the Ge concentration oscillations produce
wavefuction satellites in momentum space which can couple strongly to the
valley minima through Dresselhaus spin-orbit coupling provided that the
satellite-valley separation is approximately $4\pi/a$ in the extended
Brillouin zone as shown in Fig. 11. Importantly, the resonance in Fig. 4 is
quite wide in $\lambda$-space, which has the important implication that the
Wiggle Well structure should allow for rather large growth errors in the Ge
concentration profile while maintaining the enhanced spin-orbit effect.
Additionally, our results indicate that the spin-orbit enhancement is robust
against alloy disorder, as shown in Fig. 7.

Enhancement of both the Dresselhaus and Rashba coefficients at smaller
$\lambda$ values have also been found in Figs. 4(c) and 4(d), although these
resonances are much narrower in width than the $\lambda\approx
1.57~{}\text{nm}$ resonance, making such structures more challenging to
fabricate. Assuming that the Wiggle Well with $\lambda\approx
0.27~{}\text{nm}$ can be practically realized, however, this period is quite
attractive as it could provide the enhanced Dresselhaus spin-orbit coupling
studied here along with a huge deterministic valley splitting [18, 19].

With regards to possible applications, the enhanced spin-orbit coupling of the
Wiggle Well indicates that EDSR can be used for fast, electrically-driven
manipulations of single-spin, Loss-DiVincenzo qubits without the use of
micromagnets. Indeed, a fast, spin-orbit driven EDSR capability is one of the
main attractive features of hole-spin qubits [43, 44, 17, 45, 46], and has
recently also attracted interest in Si electron-spin qubits [47]. This
possibility is supported by our calculations in Sec. III, where it was shown
that an EDSR Rabi frequency of $\Omega_{\text{Rabi}}/B>500~{}\text{MHz/T}$ can
be obtained near the optimal Ge oscillation wavelength
$\lambda=1.57~{}\text{nm}$. Finally, we mention that the enhanced spin-orbit
coupling between the valleys may also be used to drive fast singlet–triplet
rotations near the valley-Zeeman hotspot [48].

## Acknowledgements

Research was sponsored in part by the Army Research Office (ARO) under Award
Numbers W911NF-17-1-0274 and W911NF-22-1-0090. The views, conclusions, and
recommendations contained in this document are those of the authors and are
not necessarily endorsed nor should they be interpreted as representing the
official policies, either expressed or implied, of the Army Research Office
(ARO) or the U.S. Government. The U.S. Government is authorized to reproduce
and distribute reprints for Government purposes notwithstanding any copyright
notation herein.

##  Appendix A Lattice constants in strained Si/SiGe quantum wells

To determine the lattice constants of the strained Si/SiGe heterostructure, we
use pseudomorphic boundary conditions [49, 31], where the in-plane lattice
constant $a_{\parallel}$ throughout the system is given by the relaxed lattice
constant of the Si0.7Ge0.3 barrier regions. This then also sets the lattice
spacing along the growth direction as described below.

We take the relaxed lattice constant of a Si1-$x$Ge$x$ alloy as

$a_{o}(x)=(1-x)a_{\text{Si}}+xa_{\text{Ge}},$ (34)

where $a_{\text{Si}}=0.5431~{}\text{nm}$ and
$a_{\text{Ge}}=0.5657~{}\text{nm}$ are the relaxed lattice constants of Si and
Ge, respectively. Our structure therefore has an in-plane lattice constant
$a_{\parallel}=a_{o}(0.3)=0.5499~{}\text{nm}$ throughout the entire system.
Atomic layers in our system with Ge concentration $n_{Ge,n}\neq 0.3$ are
therefore strained. Explicitly, the in-plane strain $\epsilon_{\parallel,n}$
of layer $n$ is

$\epsilon_{\parallel,n}=\frac{a_{\parallel}-a_{o}(n_{Ge,n})}{a_{o}(n_{Ge,n})}.$
(35)

For a bulk Si1-$x^{\prime}$Ge$x^{\prime}$ alloy under biaxial stress
perpendicular to the $[001]$, we have strains
$\epsilon_{xx}=\epsilon_{yy}=\epsilon_{\parallel}$ and
$\epsilon_{zz}=\epsilon_{\perp}$ that are related by [49]

$\epsilon_{\perp}(x^{\prime})=-2\frac{C_{12}(x^{\prime})}{C_{11}(x^{\prime})}\epsilon_{\parallel},$
(36)

where $C_{12}$ and $C_{11}$ are elastic constants that depend on the Ge
concentration $x^{\prime}$. For pure Si (x = 0), we have
$C_{11}(0)=165.8~{}\text{GPa}$ and $C_{12}(0)=63.9~{}\text{GPa}$, while for
pure Ge (x = 1), we have $C_{11}(1)=131.8~{}\text{GPa}$ and
$C_{12}(1)=48.3~{}\text{GPa}$ [28]. For simplicity, we assume that the elastic
constants vary linearly with the Ge concentration For a well region composed
of Si1-$x^{\prime}$Ge$x^{\prime}$, we would then have the out-of-plane lattice
constant,

$a_{\perp}(x^{\prime})=\Big{(}1+\epsilon_{\perp}(x^{\prime})\Big{)}a_{o}(x^{\prime}).$
(37)

In our system with its inhomogeneous Ge concentration profile, we take the
lattice spacing between layers $n$ and $n+1$ with Ge concentrations $n_{Ge,n}$
and $n_{Ge,n+1}$, respectively, as $a_{\perp}^{(n+1,n)}/4$, where

$a_{\perp}^{(n+1,n)}=\frac{1}{2}\Big{(}a_{\perp}(n_{Ge,n+1})+a_{\perp}(n_{Ge,n})\Big{)},$
(38)

is the average of the out-of-plane lattice constants expected for strained
regions with Ge concentration $n_{Ge,n}$ and $n_{Ge,n+1}$, respectively.

##  Appendix B Pseudospin basis transformation details

In this appendix, we provide details of the pseudospin basis introduced in
Sec. II.2 of the main text. The pseudospin basis is defined as

$\ket{\boldsymbol{k}_{\parallel}n\bar{\nu}\bar{\sigma}}=\sum_{\nu\sigma}\ket{\boldsymbol{k}_{\parallel}n\nu\sigma}U^{(n)}_{\nu\sigma,\bar{\nu}\bar{\sigma}},$
(39)

where $\bar{\nu}$ and $\bar{\sigma}$ label the new orbitals with
$\bar{\sigma}=\Uparrow,\Downarrow$ being a pseudospin label, $\nu$ and
$\sigma$ are indices of the original basis with $\sigma=\uparrow,\downarrow$
simply denoting spin, and $U^{(n)}$ is the transformation matrix of layer $n$.

Our first requirement of the new basis is that it diagonalizes the onsite
spin-orbit coupling. Following Chadi [29], spin-orbit coupling is taken to be
an intra-atomic (onsite) coupling between p orbitals and enters into the
Hamiltonian as the matrix $S$. The explicit form of the $S$ matrix of atom $j$
in atomic layer $n$ in the p-orbital subspace
$\left\\{\ket{p_{z}\uparrow},\ket{p_{x}\downarrow},\ket{p_{y}\downarrow},\ket{p_{z}\downarrow},\ket{p_{x}\uparrow},\ket{p_{y}\uparrow}\right\\}$
is

$S^{(nj)}=\frac{\Delta^{(nj)}_{\text{SO}}}{3}\begin{pmatrix}0&-1&i&0&0&0\\\
-1&0&i&0&0&0\\\ -i&-i&0&0&0&0\\\ 0&0&0&0&1&i\\\ 0&0&0&1&0&-i\\\
0&0&0&-i&i&0\\\ \end{pmatrix},$ (40)

where $\Delta_{\text{SO}}^{(nj)}$ is the spin-orbit energy. All matrix
elements involving s, s*, and d orbitals are set to zero in $S$ and are not
shown. Note that spin-orbit coupling, in principle, does exist between
d-orbitals, but is much smaller than the p-orbital couplings and is typically
neglected [27]. Also note that the spin-orbit energy depends on if the atom is
Si or Ge, but the form of the spin-orbit coupling matrix is independent of
atom type. It turns out that the $S$ matrix is diagonalized by the eigenstates
of total angular momentum [16]. Within the p-orbital subspace, these states
are given by

$\displaystyle\left|p_{1}\Uparrow\right>$
$\displaystyle=\frac{i}{\sqrt{2}}\left(\left|p_{x}\downarrow\right>-i\left|p_{y}\downarrow\right>\right),$
(41) $\displaystyle\left|p_{2}\Uparrow\right>$
$\displaystyle=\frac{1}{\sqrt{6}}\left(2\left|p_{z}\uparrow\right>-\left|p_{x}\downarrow\right>-i\left|p_{y}\downarrow\right>\right),$
(42) $\displaystyle\left|p_{3}\Uparrow\right>$
$\displaystyle=\frac{1}{\sqrt{3}}\left(\left|p_{z}\uparrow\right>+\left|p_{x}\downarrow\right>+i\left|p_{y}\downarrow\right>\right),$
(43) $\displaystyle\left|p_{1}\Downarrow\right>$
$\displaystyle=\frac{i}{\sqrt{2}}\left(\left|p_{x}\uparrow\right>+i\left|p_{y}\uparrow\right>\right),$
(44) $\displaystyle\left|p_{2}\Downarrow\right>$
$\displaystyle=\frac{1}{\sqrt{6}}\left(2\left|p_{z}\downarrow\right>+\left|p_{x}\uparrow\right>-i\left|p_{y}\uparrow\right>\right),$
(45) $\displaystyle\left|p_{3}\Downarrow\right>$
$\displaystyle=\frac{1}{\sqrt{3}}\left(\left|p_{z}\downarrow\right>-\left|p_{x}\uparrow\right>+i\left|p_{y}\uparrow\right>\right),$
(46)

where $\Uparrow,\Downarrow$ are pseudospin labels. We find
$\ket{p_{1}\bar{\sigma}}$ and $\ket{p_{2}\bar{\sigma}}$ both have an S matrix
eigenvalue of $\Delta^{(nj)}_{\text{SO}}/3$, while $\ket{p_{3}\bar{\sigma}}$
has an eigenvalue of $-2\Delta^{(nj)}_{\text{SO}}/3$. We also define
pseudospin s and s* orbitals as

$\displaystyle\left|s\Uparrow\right>=\ket{s\uparrow},$ (47)
$\displaystyle\left|s^{*}\Uparrow\right>=\ket{s^{*}\uparrow},$ (48)
$\displaystyle\left|s\Downarrow\right>=\ket{s\downarrow},$ (49)
$\displaystyle\left|s^{*}\Downarrow\right>=\ket{s^{*}\downarrow},$ (50)

and pseudospin d-orbitals as

$\displaystyle\left|d_{1}\Uparrow\right>=\frac{i}{\sqrt{2}}\left(\left|d_{zx}\downarrow\right>-i\left|d_{yz}\downarrow\right>\right),$
(51)
$\displaystyle\left|d_{2}\Uparrow\right>=\frac{1}{\sqrt{6}}\left(2\left|d_{z^{2}}\uparrow\right>-\left|d_{zx}\downarrow\right>-i\left|d_{yz}\downarrow\right>\right),$
(52)
$\displaystyle\left|d_{3}\Uparrow\right>=\frac{1}{\sqrt{3}}\left(\left|d_{z^{2}}\uparrow\right>+\left|d_{zx}\downarrow\right>+i\left|d_{yz}\downarrow\right>\right),$
(53) $\displaystyle\left|d_{4}\Uparrow\right>=\ket{d_{xy}\uparrow},$ (54)
$\displaystyle\left|d_{5}\Uparrow\right>=\ket{d_{x^{2}-y^{2}}\uparrow},$ (55)
$\displaystyle\left|d_{1}\Downarrow\right>=\frac{i}{\sqrt{2}}\left(\left|d_{zx}\uparrow\right>+i\left|d_{yz}\uparrow\right>\right),$
(56)
$\displaystyle\left|d_{2}\Downarrow\right>=\frac{1}{\sqrt{6}}\left(2\left|d_{z^{2}}\downarrow\right>+\left|d_{zx}\uparrow\right>-i\left|d_{yz}\uparrow\right>\right),$
(57)
$\displaystyle\left|d_{3}\Downarrow\right>=\frac{1}{\sqrt{3}}\left(\left|d_{z^{2}}\downarrow\right>-\left|d_{zx}\uparrow\right>+i\left|d_{yz}\uparrow\right>\right)$
(58) $\displaystyle\left|d_{4}\Downarrow\right>=\ket{d_{xy}\downarrow},$ (59)
$\displaystyle\left|d_{5}\Downarrow\right>=\ket{d_{x^{2}-y^{2}}\downarrow}.$
(60)

Note that all pseudospin s, s*, and d-orbitals are trivially eigenstates of
the spin-orbit matrix $S$ with eigenvalue $0$. In addition, note that for
$j=1,2,3$, $\ket{d_{j}\Uparrow}$ and $\ket{d_{j}\Downarrow}$ are found from
$\ket{p_{j}\Uparrow}$ and $\ket{p_{j}\Downarrow}$, respectively, by letting
$p_{x},p_{y},p_{z}\rightarrow d_{zx},d_{yz},d_{z^{2}}$. The other pseudospin d
orbitals and s orbitals are trivially related to the original basis. Note that
we adopt these altered d-orbital pseudospin states even though they posses no
spin-orbit coupling such that we obtain the pseudospin structure of the
hopping matrices in Eqs. (12 \- 14) of the main text. Failure to adopt these
d-orbital pseudospin states would result in coupling between the pseudospin
sectors for $\mathbf{k}_{\parallel}=0$.

Secondly, we require the pseudospin basis to transform the Hamiltonian in such
a way that the minimum unit cell (in the absence of an inhomogenous potential
$V_{n}$) decreases from two sites to one site for $\mathbf{k}_{\parallel}=0$.
In other words, the Fourier transformed hopping matrix, which is introduced in
Eq. (9) of the main text, must become site independent for
$\mathbf{k}_{\parallel}=0$. Naively adopting the orbitals defined in Eqs. (41
\- 60) for every site does not fulfill this requirement. However, we find this
requirement to be fulfilled if we adopt the orbitals defined in Eqs. (41 \-
60) if the $\ket{p_{1}\Uparrow}$, $\ket{p_{1}\Downarrow}$,
$\ket{d_{1}\Uparrow}$, $\ket{d_{1}\Downarrow}$, $\ket{d_{4}\Uparrow}$, and
$\ket{d_{4}\Downarrow}$, orbitals are multiplied by $(-1)^{n}$, where $n$ is
the site index in the $1D$ chain. In other words, we flip the sign of these
select orbitals on every other site. As stated in the main text, this
alternating structure for the transformation matrix occurs due to the presence
of two sublattices in the diamond crystal structure of Si as shown in Fig. 2
(a).

For clarity, we now provide the explicit form of $U^{(n)}$. We write $U^{(n)}$
as

$U^{(n)}=\begin{pmatrix}U_{s}^{(n)}&0&0\\\ 0&U_{p}^{(n)}&0\\\
0&0&U_{d}^{(n)},\end{pmatrix}$ (61)

where $U_{s}^{(n)}$, $U_{p}^{(n)}$, and $U_{d}^{(n)}$ are matrix blocks which
describe how the $s$, $p$, and $d$ orbitals transform, respectively. Here,
$U_{s}^{(n)}=I_{4\times 4}$ is just identity. The p block is given by

$U_{p}^{(n)}=\begin{pmatrix}0&\frac{2}{\sqrt{6}}&\frac{1}{\sqrt{3}}&0&0&0\\\
0&0&0&\frac{i(-1)^{n}}{\sqrt{2}}&\frac{1}{\sqrt{6}}&\frac{-1}{\sqrt{3}}\\\
0&0&0&\frac{(-1)^{n+1}}{\sqrt{2}}&\frac{-i}{\sqrt{6}}&\frac{i}{\sqrt{3}}\\\
0&0&0&0&\frac{2}{\sqrt{6}}&\frac{1}{\sqrt{3}}\\\
\frac{i(-1)^{n}}{\sqrt{2}}&\frac{-1}{\sqrt{6}}&\frac{1}{\sqrt{3}}&0&0&0\\\
\frac{(-1)^{n}}{\sqrt{2}}&\frac{-i}{\sqrt{6}}&\frac{i}{\sqrt{3}}&0&0&0\\\
\end{pmatrix},$ (62)

where the columns correspond to the pseudospin orbitals in the order
$\\{p_{1}\Uparrow,p_{2}\Uparrow,p_{3}\Uparrow,p_{1}\Downarrow,p_{2}\Downarrow,p_{3}\Downarrow\\}$,
and the rows correspond to the “standard” orbitals in the order
$\\{p_{z}\uparrow,p_{x}\uparrow,p_{y}\uparrow,p_{z}\downarrow,p_{x}\downarrow,p_{y}\downarrow\\}$.
Finally, the d block is given by

$U_{d}^{(n)}=\begin{pmatrix}0&\frac{2}{\sqrt{6}}&\frac{1}{\sqrt{3}}&0&0&0&0&0&0&0\\\
0&0&0&0&0&\frac{i(-1)^{n}}{\sqrt{2}}&\frac{1}{\sqrt{6}}&\frac{-1}{\sqrt{3}}&0&0\\\
0&0&0&0&0&\frac{(-1)^{n+1}}{\sqrt{2}}&\frac{-i}{\sqrt{6}}&\frac{i}{\sqrt{3}}&0&0\\\
0&0&0&(-1)^{n}&0&0&0&0&0&0\\\ 0&0&0&0&1&0&0&0&0&0\\\
0&0&0&0&0&0&\frac{2}{\sqrt{6}}&\frac{1}{\sqrt{3}}&0&0\\\
\frac{i(-1)^{n}}{\sqrt{2}}&\frac{-1}{\sqrt{6}}&\frac{1}{\sqrt{3}}&0&0&0&0&0&0&0\\\
\frac{(-1)^{n}}{\sqrt{2}}&\frac{-i}{\sqrt{6}}&\frac{i}{\sqrt{3}}&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&(-1)^{n}&0\\\ 0&0&0&0&0&0&0&0&0&1\\\ \end{pmatrix},$ (63)

where the columns correspond to the pseudospin orbitals in the order
$\\{d_{1}\Uparrow,d_{2}\Uparrow,d_{3}\Uparrow,d_{4}\Uparrow,d_{5}\Uparrow,d_{1}\Downarrow,d_{2}\Downarrow,d_{3}\Downarrow,d_{4}\Downarrow,d_{5}\Downarrow\\}$,
and the rows correspond to the “standard” orbitals in the order
$\\{d_{z^{2}}\uparrow,d_{zx}\uparrow,d_{yz}\uparrow,d_{xy}\uparrow,d_{x^{2}-y^{2}}\uparrow,d_{z^{2}}\downarrow,d_{zx}\downarrow,d_{yz}\downarrow,d_{xy}\downarrow,d_{x^{2}-y^{2}}\downarrow\\}$.
Notice the $(-1)^{n}$ factors that flips the signs of the
$\ket{p_{1}\Uparrow}$, $\ket{p_{1}\Downarrow}$, $\ket{d_{1}\Uparrow}$,
$\ket{d_{1}\Downarrow}$, $\ket{d_{4}\Uparrow}$, and $\ket{d_{4}\Downarrow}$
orbitals on every odd site. This allows for a unit cell containing only one
site for $\mathbf{k}_{\parallel}=0$ as described above.

##  Appendix C Symmetry argument for the hopping structure of the Rashba and
Dresselhaus Hamiltonian components

In Sec. II.3 of the main text, we found the Fourier transformed hopping matrix
$\widetilde{T}^{(n)}_{\pm}(\mathbf{k}_{\parallel})$ to have the form

$\widetilde{T}^{(n)}_{\pm}(\mathbf{k}_{\parallel})=\widetilde{T}^{(n)}_{o}+\widetilde{T}^{(n)}_{R}(\mathbf{k}_{\parallel})\pm\widetilde{T}^{(n)}_{D}(\mathbf{k}_{\parallel})+\mathcal{O}(\mathbf{k}_{\parallel}^{2}),$
(64)

where $\widetilde{T}^{(n)}_{R}(\mathbf{k}_{\parallel})$ and
$\widetilde{T}^{(n)}_{D}(\mathbf{k}_{\parallel})$ are the Rashba and
Dresselhaus hopping matrices, respectively, and are given by

$\displaystyle\widetilde{T}^{(n)}_{R}(\mathbf{k}_{\parallel})$
$\displaystyle=\Phi^{(n)}\left(k_{y}\sigma_{x}-k_{x}\sigma_{y}\right),$ (65)
$\displaystyle\widetilde{T}^{(n)}_{D}(\mathbf{k}_{\parallel})$
$\displaystyle=\Phi^{(n)}\left(k_{x}\sigma_{x}-k_{y}\sigma_{y}\right),$ (66)

with $\Phi^{(n)}$ being a real-valued $10\times 10$ matrix with vanishing
diagonal elemental. Importantly, the sign in front of the Rashba hopping
matrix $\widetilde{T}^{(n)}_{R}$ is site independent, while the sign of the
Dresselhaus hopping matrix $\widetilde{T}^{(n)}_{D}$ changes sign between
every site, as indicated by the $\pm$ in Eq. (64). This can be understood as
originating from the diamond crystal structure of the Si by the following
symmetry argument; Let us consider the case of pure Si such that
$\left\\{\Phi^{(n)},\widetilde{T}^{(n)}_{R},\widetilde{T}^{(n)}_{D}\right\\}\rightarrow\left\\{\Phi,\widetilde{T}_{R},\widetilde{T}_{D}\right\\}$
all lose their dependence on the layer index. Next, note that under a $C_{4}$
rotation about the $z$-axis (growth axis), we have
$\\{k_{x},k_{y},\sigma_{x},\sigma_{y}\\}\rightarrow\\{k_{y},-k_{x},\sigma_{y},-\sigma_{x}\\}$.
This leaves invariant the Rashba term in Eq. (65) and flips the sign of the
Dresselhaus term in Eq. (66). Finally, performing the same $C_{4}$ rotation on
our diamond crystal structure in Fig. 2 (a) transforms red atoms into blue
atoms and vice versa in the sense that the nearest neighbor vectors of even
and odd atomic layers swap. In other words, the two sublattices of the Si
lattice swap. This in turn, swaps $\widetilde{T}_{+}$ and $\widetilde{T}_{-}$
within the tight binding chain. Clearly then, the Rashba hopping term
$\widetilde{T}_{R}$, being invariant under a $C_{4}$ rotation, should contain
the part common to $\widetilde{T}_{+}$ and $\widetilde{T}_{-}$. In contrast,
the Dresselhaus hopping term $\widetilde{T}_{D}$ should contain the part which
is different between $\widetilde{T}_{+}$ and $\widetilde{T}_{-}$, since it
flips sign under a $C_{4}$ rotation. This then explains the $\pm$ in front of
$\widetilde{T}^{(n)}_{D}$ in Eq. (64).

##  Appendix D $\Omega$ and $\Phi$ matrices

In Sec. II.3 of the main text, we introduced the $10\times 10$ matrices
$\Omega^{(n)}$ and $\Phi^{(n)}$ as components of the hopping matrices in Eq.
(12 \- 14). These matrices can further decomposed into the block forms,

$\displaystyle\Omega^{(n)}$
$\displaystyle=\begin{bmatrix}\Omega^{(n)}_{00}&\Omega^{(n)}_{01}&0\\\
-\Omega_{01}^{(n)T}&\Omega^{(n)}_{11}&0\\\
0&0&\Omega^{(n)}_{22}\end{bmatrix},$ (67) $\displaystyle\Phi^{(n)}$
$\displaystyle=\begin{bmatrix}\Phi_{00}^{(n)}&\Phi_{01}^{(n)}&\Phi_{02}^{(n)}\\\
\Phi_{01}^{(n)T}&\Phi_{11}^{(n)}&\Phi_{12}^{(n)}\\\
-\Phi_{02}^{(n)T}&\Phi_{12}^{(n)T}&0\end{bmatrix},$ (68)

where the shape of the diagonal blocks are $5\times 5$, $4\times 4$, and
$1\times 1$, respectively, and the diagonal block matrices satisfy
$\Omega_{ii}^{(n)T}=\Omega_{ii}^{(n)}$ and
$\Phi_{ii}^{(n)T}=-\Phi_{ii}^{(n)}$. Note that this implies that all diagonal
elements of $\Phi^{(n)}$ are zero. Here, the orbital ordering used is
$\left\\{s,s^{*},p_{1},d_{2},d_{3},p_{2},p_{3},d_{1},d_{4},d_{5}\right\\}$.
Generically, these matrices depend on the layer index $n$ due to Ge
concentration changing from layer to layer. In the case of a uniform Ge
concentration, however, this layer dependence goes away,
$\left\\{\Omega^{(n)},\Phi^{(n)}\right\\}\rightarrow\left\\{\Omega,\Phi\right\\}$.
In the particular case of an unstrained Si system, the $\Omega$ matrix blocks
(in eV) are given by

$\displaystyle\Omega_{00}$ $\displaystyle=\begin{pmatrix}-3.73&-2.78&0&0&0\\\
-2.78&-9.03&0&0&0\\\ 0&0&-0.73&-1.71&2.42\\\ 0&0&-1.71&1.63&2.03\\\
0&0&2.42&2.03&0.19\end{pmatrix},$ (69) $\displaystyle\Omega_{11}$
$\displaystyle=\begin{pmatrix}0.73&0&1.71&-2.42\\\ 0&0.73&-2.42&-1.71\\\
1.71&-2.42&1.24&0\\\ -2.42&-1.71&0&1.24\\\ \end{pmatrix},$ (70)
$\displaystyle\Omega_{22}$ $\displaystyle=3.06,$ (71)
$\displaystyle\Omega_{01}$ $\displaystyle=\begin{pmatrix}2.74&1.94&0&-2.59\\\
2.89&2.05&0&-0.90\\\ 2.15&-3.04&0.19&0\\\ -2.07&-1.60&0.70&-2.92\\\
-1.60&0.94&-0.99&-2.07\end{pmatrix},$ (72)

and the $\Phi$ matrix blocks (in eV $\cdot$ Å) are given by

$\displaystyle\Phi_{00}$ $\displaystyle=\begin{pmatrix}0&0&3.23&1.43&-2.03\\\
0&0&3.40&0.50&-0.70\\\ -3.23&-3.40&0&-1.25&-0.89\\\ -1.43&-0.50&1.25&0&1.72\\\
2.03&0.70&0.89&-1.72&0\end{pmatrix},$ (73) $\displaystyle\Phi_{11}$
$\displaystyle=\begin{pmatrix}0&3.57&-0.15&-0.10\\\ -3.57&0&-0.1&0.15\\\
-0.15&-0.1&0&-1.16\\\ -0.10&-0.15&1.16&0\\\ \end{pmatrix},$ (74)
$\displaystyle\Phi_{02}$
$\displaystyle=\begin{pmatrix}0&0&2.66&1.72&2.43\end{pmatrix},$ (75)
$\displaystyle\Phi_{12}$
$\displaystyle=\begin{pmatrix}-1.54&2.17&2.98&0\end{pmatrix}.$ (76)

Note that these $\Omega$ and $\Phi$ are precisely the matrices used in the
simplified model of Sec. IV, where Ge atoms are treated as Si atoms with
orbitals shifted up by a constant energy $E_{\text{Ge}}$.

## References

  * Loss and DiVincenzo [1998] D. Loss and D. P. DiVincenzo, Quantum computation with quantum dots, Phys. Rev. A 57, 120 (1998).
  * Hanson  _et al._ [2007] R. Hanson, L. P. Kouwenhoven, J. R. Petta, S. Tarucha, and L. M. K. Vandersypen, Spins in few-electron quantum dots, Rev. Mod. Phys. 79, 1217 (2007).
  * Burkard  _et al._ [2021] G. Burkard, T. D. Ladd, J. M. Nichol, A. Pan, and J. R. Petta, Semiconductor spin qubits, arXiv:2112.08863 (2021).
  * Chatterjee  _et al._ [2021] A. Chatterjee, P. Stevenson, S. D. Franceschi, A. Morello, N. P. de Leon, and F. Kuemmeth, Semiconductor qubits in practice, Nature Reviews Physics 3, 157 (2021).
  * Zhang  _et al._ [2018] X. Zhang, H.-O. Li, G. Cao, M. Xiao, G.-C. Guo, and G.-P. Guo, Semiconductor quantum computation, National Science Review 6, 32 (2018).
  * Zwanenburg  _et al._ [2013] F. A. Zwanenburg, A. S. Dzurak, A. Morello, M. Y. Simmons, L. C. L. Hollenberg, G. Klimeck, S. Rogge, S. N. Coppersmith, and M. A. Eriksson, Silicon quantum electronics, Rev. Mod. Phys. 85, 961 (2013).
  * Gyure  _et al._ [2021] M. F. Gyure, A. A. Kiselev, R. S. Ross, R. Rahman, and C. G. V. de Walle, Materials and device simulations for silicon qubit design and optimization, MRS Bulletin 46, 634 (2021).
  * Petta  _et al._ [2005] J. R. Petta, A. C. Johnson, J. M. Taylor, E. A. Laird, A. Yacoby, M. D. Lukin, C. M. Marcus, M. P. Hanson, and A. C. Gossard, Coherent manipulation of coupled electron spins in semiconductor quantum dots, Science 309, 2180 (2005).
  * Assali  _et al._ [2011] L. V. C. Assali, H. M. Petrilli, R. B. Capaz, B. Koiller, X. Hu, and S. Das Sarma, Hyperfine interactions in silicon quantum dots, Phys. Rev. B 83, 165301 (2011).
  * Pioro-Ladrière  _et al._ [2008] M. Pioro-Ladrière, T. Obata, Y. Tokura, Y.-S. Shin, T. Kubo, K. Yoshida, T. Taniyama, and S. Tarucha, Electrically driven single-electron spin resonance in a slanting Zeeman field, Nat. Phys. 4, 776 (2008).
  * Xue  _et al._ [2022] X. Xue, M. Russ, N. Samkharadze, B. Undseth, A. Sammak, G. Scappucci, and L. M. K. Vandersypen, Quantum logic with spin qubits crossing the surface code threshold, Nature 601, 343 (2022).
  * Noiri  _et al._ [2022] A. Noiri, K. Takeda, T. Nakajima, T. Kobayashi, A. Sammak, G. Scappucci, and S. Tarucha, Fast universal quantum gate above the fault-tolerance threshold in silicon, Nature 601, 338 (2022).
  * Mills  _et al._ [2021] A. R. Mills, C. R. Guinn, M. J. Gullans, A. J. Sigillito, M. M. Feldman, E. Nielsen, and J. R. Petta, Two-qubit silicon quantum processor with operation fidelity exceeding 99, arXiv:2111.11937 (2021).
  * Golovach  _et al._ [2006] V. N. Golovach, M. Borhani, and D. Loss, Electric-dipole-induced spin resonance in quantum dots, Phys. Rev. B 74, 165319 (2006).
  * Rashba [2008] E. I. Rashba, Theory of electric dipole spin resonance in quantum dots: Mean field theory with gaussian fluctuations and beyond, Phys. Rev. B 78, 195302 (2008).
  * Winkler [2003] R. Winkler,  _Spin-Orbit Coupling Effects in Two-Dimensional Electron and Hole Systems_ (Springer, 2003).
  * Terrazos  _et al._ [2021] L. A. Terrazos, E. Marcellina, Z. Wang, S. N. Coppersmith, M. Friesen, A. R. Hamilton, X. Hu, B. Koiller, A. L. Saraiva, D. Culcer, and R. B. Capaz, Theory of hole-spin qubits in strained germanium quantum dots, Physical Review B 103, 125201 (2021).
  * McJunkin  _et al._ [2021] T. McJunkin, B. Harpt, Y. Feng, M. Losert, R. Rahman, J. P. Dodson, M. A. Wolfe, D. E. Savage, M. G. Lagally, S. N. Coppersmith, M. Friesen, R. Joynt, and M. A. Eriksson, SiGe quantum wells with oscillating Ge concentrations for quantum dot qubits, arXiv:2112.09765 (2021).
  * Feng and Joynt [2022] Y. Feng and R. Joynt, Enhanced valley splitting in Si layers with oscillatory Ge concentration, Phys. Rev. B 106, 085304 (2022).
  * Nestoklon  _et al._ [2008] M. O. Nestoklon, E. L. Ivchenko, J.-M. Jancu, and P. Voisin, Electric field effect on electron spin splitting in $\mathrm{Si}\mathrm{Ge}/\mathrm{Si}$ quantum wells, Phys. Rev. B 77, 155328 (2008).
  * Prada  _et al._ [2011] M. Prada, G. Klimeck, and R. Joynt, Spin–orbit splittings in Si/SiGe quantum wells: from ideal Si membranes to realistic heterostructures, New Journal of Physics 13, 013009 (2011).
  * Bychkov and Rashba [1984] Y. A. Bychkov and E. I. Rashba, Oscillatory effects and the magnetic susceptibility of carriers in inversion layers, Journal of Physics C: Solid State Physics 17, 6039 (1984).
  * Dresselhaus [1955] G. Dresselhaus, Spin-orbit coupling effects in zinc blende structures, Phys. Rev. 100, 580 (1955).
  * Ivchenko  _et al._ [1996] E. L. Ivchenko, A. Y. Kaminski, and U. Rössler, Heavy-light hole mixing at zinc-blende (001) interfaces under normal incidence, Phys. Rev. B 54, 5852 (1996).
  * Vervoort  _et al._ [1997] L. Vervoort, R. Ferreira, and P. Voisin, Effects of interface asymmetry on hole subband degeneracies and spin-relaxation rates in quantum wells, Phys. Rev. B 56, R12744 (1997).
  * Slater and Koster [1954] J. C. Slater and G. F. Koster, Simplified LCAO method for the periodic potential problem, Phys. Rev. 94, 1498 (1954).
  * Jancu  _et al._ [1998] J.-M. Jancu, R. Scholz, F. Beltram, and F. Bassani, Empirical ${\mathrm{spds}}^{*}$ tight-binding calculation for cubic semiconductors: General method and material parameters, Phys. Rev. B 57, 6493 (1998).
  * Niquet  _et al._ [2009] Y. M. Niquet, D. Rideau, C. Tavernier, H. Jaouen, and X. Blase, Onsite matrix elements of the tight-binding hamiltonian of a strained crystal: Application to silicon, germanium, and their alloys, Phys. Rev. B 79, 245201 (2009).
  * Chadi [1977] D. J. Chadi, Spin-orbit splitting in crystalline and compositionally disordered semiconductors, Phys. Rev. B 16, 790 (1977).
  * Ren  _et al._ [1982] S. Y. Ren, J. D. Dow, and D. J. Wolford, Pressure dependence of deep levels in GaAs, Phys. Rev. B 25, 7661 (1982).
  * Schäffler [1997] F. Schäffler, High-mobility Si and Ge structures, Semiconductor Science and Technology 12, 1515 (1997).
  * Friesen  _et al._ [2007] M. Friesen, S. Chutia, C. Tahan, and S. N. Coppersmith, Valley splitting theory of $\mathrm{Si}\mathrm{Ge}/\mathrm{Si}/\mathrm{Si}\mathrm{Ge}$ quantum wells, Phys. Rev. B 75, 115318 (2007).
  * Wuetz  _et al._ [2021] B. P. Wuetz, M. P. Losert, S. Koelling, L. E. A. Stehouwer, A.-M. J. Zwerver, S. G. J. Philips, M. T. Mądzik, X. Xue, G. Zheng, M. Lodari, S. V. Amitonov, N. Samkharadze, A. Sammak, L. M. K. Vandersypen, R. Rahman, S. N. Coppersmith, O. Moutanabbir, M. Friesen, and G. Scappucci, Atomic fluctuations lifting the energy degeneracy in Si/SiGe quantum dots, arXiv:2112.09606 (2021).
  * Nestoklon  _et al._ [2006] M. O. Nestoklon, L. E. Golub, and E. L. Ivchenko, Spin and valley-orbit splittings in $\mathrm{Si}\mathrm{Ge}/\mathrm{Si}$ heterostructures, Phys. Rev. B 73, 235334 (2006).
  * Veldhorst  _et al._ [2015] M. Veldhorst, R. Ruskov, C. H. Yang, J. C. C. Hwang, F. E. Hudson, M. E. Flatté, C. Tahan, K. M. Itoh, A. Morello, and A. S. Dzurak, Spin-orbit coupling and operation of multivalley spin qubits, Phys. Rev. B 92, 201401 (2015).
  * Ruskov  _et al._ [2018] R. Ruskov, M. Veldhorst, A. S. Dzurak, and C. Tahan, Electron $g$-factor of valley states in realistic silicon quantum dots, Phys. Rev. B 98, 245424 (2018).
  * Boykin  _et al._ [2004] T. B. Boykin, G. Klimeck, M. A. Eriksson, M. Friesen, S. N. Coppersmith, P. von Allmen, F. Oyafuso, and S. Lee, Valley splitting in strained silicon quantum wells, Applied Physics Letters 84, 115 (2004).
  * Townsend [2012] J. S. Townsend,  _A Modern Approach to Quantum Mechanics_ (University Science Books, 2012).
  * Jantsch  _et al._ [2002] W. Jantsch, Z. Wilamowski, N. Sandersfeld, M. Mühlberger, and F. Schäffler, Spin lifetimes and g-factor tuning in si/sige quantum wells, Physica E: Low-dimensional Systems and Nanostructures 13, 504 (2002).
  * Malissa  _et al._ [2004] H. Malissa, W. Jantsch, M. Mühlberger, F. Schäffler, Z. Wilamowski, M. Draxler, and P. Bauer, Anisotropy of g-factor and electron spin resonance linewidth in modulation doped sige quantum wells, Applied Physics Letters 85, 1739 (2004).
  * Kittel and McEuen [2019] C. Kittel and P. McEuen,  _Introduction to solid state physics_ (John Wiley & Sons, 2019).
  * Li and Dery [2011] P. Li and H. Dery, Spin-orbit symmetries of conduction electrons in silicon, Phys. Rev. Lett. 107, 107203 (2011).
  * Maurand  _et al._ [2016] R. Maurand, X. Jehl, D. Kotekar-Patil, A. Corna, H. Bohuslavskyi, R. Laviéville, L. Hutin, S. Barraud, M. Vinet, M. Sanquer, and S. D. Franceschi, A cmos silicon spin qubit, Nature Communications 7, 13575 (2016).
  * Watzinger  _et al._ [2018] H. Watzinger, J. Kukučka, L. Vukušić, F. Gao, T. Wang, F. Schäffler, J.-J. Zhang, and G. Katsaros, A germanium hole spin qubit, Nature Communications 9, 3902 (2018).
  * Wang  _et al._ [2022] K. Wang, G. Xu, F. Gao, H. Liu, R.-L. Ma, X. Zhang, Z. Wang, G. Cao, T. Wang, J.-J. Zhang, D. Culcer, X. Hu, H.-W. Jiang, H.-O. Li, G.-C. Guo, and G.-P. Guo, Ultrafast coherent control of a hole spin qubit in a germanium quantum dot, Nature Communications 13, 206 (2022).
  * Kloeffel and Loss [2013] C. Kloeffel and D. Loss, Prospects for spin-based quantum computing in quantum dots, Annual Review of Condensed Matter Physics 4, 51 (2013).
  * Gilbert  _et al._ [2022] W. Gilbert, T. Tanttu, W. H. Lim, M. Feng, J. Y. Huang, J. D. Cifuentes, S. Serrano, P. Y. Mai, R. C. C. Leon, C. C. Escott, K. M. Itoh, N. V. Abrosimov, H.-J. Pohl, M. L. W. Thewalt, F. E. Hudson, A. Morello, A. Laucht, C. H. Yang, A. Saraiva, and A. S. Dzurak, On-demand electrical control of spin qubits, arXiv:2201.06679 (2022).
  * Jock  _et al._ [2022] R. M. Jock, N. T. Jacobson, M. Rudolph, D. R. Ward, M. S. Carroll, and D. R. Luhman, A silicon singlet–triplet qubit driven by spin-valley coupling, Nature Communications 13, 641 (2022).
  * Van de Walle and Martin [1986] C. G. Van de Walle and R. M. Martin, Theoretical calculations of heterojunction discontinuities in the Si/Ge system, Phys. Rev. B 34, 5621 (1986).