Spin-orbit splittings in Si/SiGe quantum wells

We present a calculation of the wavevector-dependent subband level splitting from spin-orbit coupling in Si/SiGe quantum wells. We first use the effective-mass approach, where the splittings are parameterized by separating contributions from the Rashba and Dresselhaus terms. We then determine the parameters by fitting tight-binding numerical results obtained using the quantitative nanoelectronic modeling tool, NEMO-3D. We describe the relevant parameters as a function of applied electric field and well width in our numerical simulations. For a silicon membrane, we find the bulk Rashba parameter to be linear in field, $\alpha = \alpha^1E_z$ with $\alpha^1 \simeq 2\times$ 10 $^{-5}$nm$^{-2}$. The dominant contribution to the spin-orbit splitting is from Dresselhaus-type terms, and the magnitude for a typical flat SiGe/Si/SiGe quantum well can be as high as 1$\mu$eV.

We present a calculation of the wavevector-dependent subband level splitting from spin-orbit coupling in Si/SiGe quantum wells. We first use the effective-mass approach, where the splittings are parameterized by separating contributions from the Rashba and Dresselhaus terms. We then determine the parameters by fitting tight-binding numerical results obtained using the quantitative nanoelectronic modeling tool, NEMO-3D. We describe the relevant parameters as a function of applied electric field and well width in our numerical simulations. For a silicon membrane, we find the bulk Rashba parameter to be linear in field, α = α 1 Ez with α 1 ≃ 2×10 −5 nm −2 . The dominant contribution to the spin-orbit splitting is from Dresselhaus-type terms, and the magnitude for a typical flat SiGe/Si/SiGe quantum well can be as high as 1µeV.

I. I. INTRODUCTION
Silicon is a leading candidate material for spin-based quantum information processing 1,2 . Its spin-orbit coupling (SOC) is relatively weak and the hyperfine coupling can be eliminated by isotopic purification 3 . This means that the spin lifetimes should be long.
One way to measure a spin lifetime is to use electron spin resonance (ESR) on two-dimensional electron gases (2DEGs) in silicon quantum wells (QWs) 4,5,6 , where the D'yakonov-Perel' 7 mechanism accounts for the relaxation. SOC may also be measured directly, using EST 8,9 or photocurrents 10 . In order to compare with experiment, however, we need the wavevector-dependent SOC Hamiltonian, which must be calculated atomistically. That is the aim of this paper. We shall focus on Si layers grown in the [001] direction with Si x Ge 1−x layers on either side.
Most calculations of the DP relaxation have used the Rashba Hamiltonian 11 , which is of the form (1) σ i are the Pauli matrices and k i are the two-dimensional wavevector components. N is the number of atomic Si layers in the well. v is the valley degree of freedom. We focus on the lowest electric subband, in which case the valley degree of freedom is two-valued 12 and v is a two-by-two matrix. α (E z ) depends (in lowest order) linearly on E z , the external electric field, and often only this term is kept. E z is of order 1 − 5 × 10 7 V/m in heterostructures or MODFET devices 13 . The large magnitude of the field makes it important to examine the linearity assumption, and that is one of the purposes of this paper. The Rashba is generally thought of as a bulk effect. However, de Andrade e Silva et al. pointed out that surface effects may also be important 14 , and this assumption should also be re-evaluated. In addition, when a detailed treatment of the surface effects was done by Nestoklon et al. in the absence of an applied field E z , they showed that in this case one obtains a term if N is odd 15,16 . This is a Dresselhaus-like term 17 in that it arises from inversion asymmetry. The lowestorder term in E z for β (E z , N, v) is a constant, so that if N is odd then we get a surface-induced spin splitting even in the absence of an external field. The symmetry considerations refer to an ideal freestanding Si layer, or a layer that is sandwiched between two identical Si x Ge 1−x layers that are treated in the virtual crystal approximation. Real Si x Ge 1−x layers have substitutional disorder that destroys all symmetries.
The purpose of this paper is to determine the functions α (E z , N, v) and β (E z , N, v) both for the ideal, freestanding case and for the case of a well confined in a Si/Si x Ge 1−x heterostructure. These functions determine the spin properties of electrons in Si quantum wells. We shall focus particularly on the question of which term dominates for quantum wells with realistic values of E z and N.
The numerical tight-binding calculations are performed with NEMO-3D 18 on nanoHUB.org computational resources 19 . In NEMO-3D, atoms are represented explicitly in the sp 3 d 5 s * tight-binding model, and the valence force field (VFF) method is employed to minimize strain 20 . NEMO-3D enables the calculation of localized states on a QW and their in-plane dispersion relation with a very high degree of precision, allowing to extract the splittings along the in-plane directions in k space. Note that earlier calculations of α and β have used the virtual-crystal approximation for the interfaces. Since both of these quantities are dominated by atomic-scale interface effects, this is a rather crude approach. This paper is organized as follows: Sec. II discuss the symmetry operations of a silicon membrane. In Sec. III, expressions for the SOC in a δ-functional effective mass approach are given. We present a qualitative picture in Sec. IV. Sec. V contains the numerical results for ideal Si QWs and Sec. VI, the ones for a SiGe/Si/SiGe het-erostructure. We conclude in Sec. VII with a summary of the results obtained.

II. SYMMETRY
We give here for clarity the symmetry operations of this system in the ideal case, since they are far from obvious. We stress that this is only an expanded discussion of the analysis already given by Nestoklon etal. 15 . The lattice considered as a bulk sample has the diamond structure with a tetragonal distortion due to the Si x Ge 1−x layers: the [001] axis along the growth direction is compressed relative to the in-plane [100] and [010] axes. In a symmorphic lattice this would simply reduce the point group symmetry from the cubic group O h , with 48 operations, to the tetragonal group D 4h , with 16. (Recall that a symmorphic space group is one that is generated by translations and by rotations and reflections about a point. A nonsymmorphic lattice requires combined operations such as screw axis and glide plane operations in its generating set.) For the diamond lattice, a nonsymmorphic lattice, there is no point group, but the factor group still has 16 operations, realized as follows. There are 4 proper rotations through 180 • using the 2-fold [100], [010], and [001] axes passing through the origin, which may be taken at the position of any atom. In addition there are two ±90 • improper rotations about the [001] axis: the rotation is followed by a reflection in the (001) plane. There are also reflections with respect to the (110) and (110) planes. Any of these 8 operations may be combined with an inversion through a point midway between any pair of nearest neighbors. When we consider a layer, all 16 operations might appear to preserve the positions of the interfaces, whose presence would therefore not reduce the symmetry. However, this turns out to be not quite the case.
For a layer with an odd number N of atomic layers, take the origin (0, 0, 0) at an atom in the central plane. Then the other atoms in this plane lie at the points (a x /2) [n 1 (1, 1, 0) + n 2 (1, −1, 0)], where a x is the in-plane lattice constant and n 1 and n 2 denote any integer. Each atom in the central plane is accompanied by another one in the plane z = a z /4 shifted from it by (a x /4, a x /4, a z /4) , where a z is the lattice constant in the z-direction. The atoms in the plane at z = −a z /4 are shifted with respect to the atoms in the central plane by (a x /4, −a x /4, −a z /4) . The atoms in the planes at z = ±a z /2 are shifted with respect to the atoms in the central plane by (a x /2, 0, ±a z /2) . The positions of other atoms can be found by translating this layer of thickness a z by integer multiples of (0, 0, a z ) , so to understand the symmetry properties of the full layer it is sufficient to consider a layer with N = 5. Simple reflection in the z = 0 plane is not a symmetry, since it interchanges the z = a z /4 and z = −a z /4 layers, whose in-plane shift is (a x /2) (0, 1, 0) . The z = 0 plane is also not a glide plane, since following the reflection by the translation of (a x /2) (0, 1, 0) to restore the z = a z /4 and z = −a z /4 layers would change the z = ±a z /2 layers. The 6 rotations using x-, y-, and z-axes passing through the origin are readily seen to be symmetries of the layer, as are the reflections through the (110) and 110 planes. They take the point (x, y, z) into the points {(x, y, z) , (x, −y, −z) , (−x, y, −z) , (−x, −y, z) , (−y, x, −z) , (y, −x, −z) , (y, x, z) , (−y, −x, z)} which is the group D 2d . This is a true point group, and the space group is therefore symmorphic. Spin-orbit effects come from terms linearly proportional to σ i , the spin operators that transform as pseudovectors, while the electric field E = (0, 0, E z ) and in-plane momentum k = (k x , k y , 0) transform as vectors, the same as the coordinates. Our interest here is in combinations of these three quantities. Using the above list of operations, we find in zeroth order in E z that there is one invariant term in the Hamiltonian of the form k x σ x −k y σ y . In first order in E z , there is only the Rashba term E z (k x σ y − k y σ x ) (which is of course invariant under all isometries). Either of these terms can be multiplied by any even function of E z , which is invariant under all the operations. Thus we find that when N is an odd number, β (N, E z , v) is an even function of E z and α (E z , N, v) is an odd function of E z .
For even N, take the origin (0, 0, 0) at the center of a bond between atoms at positions (a x /8, a x /8, a z /8) and (−a x /8, −a x /8, −a z /8) . The origin is then a center of inversion. The rotation through 180 • about the 110 axis is a symmetry operation.
The (110) axis is a screw axis since the 180 • rotation about this axis must be accompanied by a translation through (a x /4, a x /4, 0) to be a symmetry operation. The same is true for the 180 • rotation about the (001) axis. The 8 operations of the factor group obtained by combining these operations take the point (x, y, z) into the points {(x, y, z) , (−x, −y, −z) , (−y, −x, −z) , (y, x, z) , Modulo translations, this is isomorphic to the group D 2h . Because of the appearance of the translations, this is not a true point group and the space group is not symmorphic.
Its action in the Hilbert space reduces in many cases to projective rather than faithful representations of D 2h . However, this does not affect the symmetry analysis of the Hamiltonian. In zeroth order in E z , the group does not allow any combination of terms of the form k i σ j , since all such terms change sign under inversion. In first order, we again have the Rashba term E z (k x σ y − k y σ x ) . Again, multiplication of this expression by any even function of E z is permissible. Thus we find that when N is an even number, α (N, E z , v) is an odd function of E z and β (E z , N, v) is an odd function of E z .

III. VALLEY AND SPIN ORBIT COUPLING
We take into account the above symmetry arguments and employ a δ-functional approach for the interfaceinduced valley mixing 21 .
In the envelope-function picture 22 , the Bloch periodic functions (valleys) at the interfaces |± =u ±k0 e ±ik0L mix, resulting in valley split- , with ∆ and φ v being phenomenological parameters, and k 0 the conduction band minima. |Φ(z i )| 2 is the zero-field value of the envelope function at the interfaces. Nestoklon et al. 15 extend this work by introducing valley-orbit and spin-orbit mixing, as a spin-dependent reflection of the wavefunction at the interfaces, situated at z=z u,d (see Fig. 1(a)).
Following their approach, we introduce next H D and H R as δ-functional perturbations in the lowest spinor-valleys functions at the interfaces. We consider the hight-symmetry directions x ′ [110] and y ′ [110], along which the spin eigenstates are parallel (see Fig. 1 (b)). In this rotated basis, we have: As mentioned in the text, a QW with N even is isomorphic with the D 2h , which contains full inversion symmetry. Under the D 2h operations, both H R and H D change sign: Taking into account these symmetry arguments, we introduce a δ-functional SOC, where the effective potential is at both interfaces situated in z u and z d , as depicted in Fig. 1(a). The advantage of the δ-functional is that it allows to consider the symmetry arguments in a simple way, but we shall recover the common notation of (1) and (2). The Hamiltonian in the δ-functional approach can be expressed in terms of the Pauli matrices in the valley-space,ŝ i : where we have introduced a i and b i as parameters that determine the strength of Rashba and Dresselhaus SOC. We shall express those in terms of the commonly used α and β. Eq. (3) can be divided in two contributions: the intra-valley SOC, which is Eq. (3) with i = 0, and the inter-valley SOC, for i = x, y: We first consider the intra-valley contribution, H SO1 (i = 0). s 0 is then the identity matrix in the valley space, and thus, the spin-mixing parameters, a 0 and b 0 are real. Integrating (3) we recover the usual notation: where we have identified α and β in terms of the value of the envelope functions at the interfaces, Note that the parity is contained in the Dresselhaus term, which would be non-zero for N odd and E z = 0 (or |Φ(z u )| 2 =|Φ(z d )| 2 )|). The Hamiltonian of Eq. (4) mixes spins of the same valley, which we denote as {|+ } and {|− }. Remember that in this basis, the two lowest eigenstates have even or odd parity 21 , ϕ e/o = e −iφv /2 ; ±ηe iφv/2 , with η = sgn{cos (k 0 N a z /4)} determining the parity of the ground state, and φ v defined in 21 . We include the spin degree of freedom to diagonalize Eq. (4), and find the eigen-vectors in the basis , so the eigen-vectors are eigenstates ofσ x ′ or toσ y ′ , depending on the phase φ 1 given by the direction of the in-plane wave vector, k. Along thex ′ orŷ ′ directions, the spin 'up' and 'down' states split in both valleys by the same We consider next the valley-mixing SOC, H SO2 , which is Eq. (3) with i = x, y.
The parameters, b x =b y =b=|b|e iφ β and a x =a y =a=|a|e iφα (a z =b z =0) are now complex, as they mix different valleys, and deter-mine the valley-spin coupling strength. As previously noted by Nestoklon et al., a Si QW possesses mirror rotation operation S 4 , resulting in a relative phase change of φ α (φ β ) for the Rashba (Dresselhaus) interaction at either interface. Combined with the absence (existence) of inversion center for even (odd) N , a change of sign is also observed in the Rashba (both Rashba and Dresselhaus) terms. Taking into account these symmetry arguments, the valley-mixing Hamiltonian has the shape: The valley-spin orbit Hamiltonian in the bi-spinor basis for the lowest two valleys will thus consist of a 4×4 matrix, which can be expressed in terms of the Pauli matrices for the valleys s i and for the spin, σ i as: with: where φ 0 = k 0 L and φ 0i = φ 0 − φ i . We note that translation of the vector r by a three-dimensional Bravaislattice vector a results in multiplication of the Bloch functions |± by the factors exp (±ik 0 L) 23 , and thus a phase φ 0 =k 0 L appears in the valley-mixing terms. with where we have defined: Hence, we can write the SOC Hamiltonian in a compact way, merging Eq. (4) and Eq. (8): Along k x ′ (k y ′ ), the eigenvectors are eigenstates ofσ y ′ (σ x ′ ). The inter-valley term i = z, has a relative change of sign for the splittings of the spin 'up' and 'down' states in either valley, ∆ 2 ∝ ±2k(β z ± α z ), so there is a valleydependent spin splitting, as depicted in the inset of Fig.  2: |ε ↑ − ε ↓ |=|∆ 2 ± ∆ 1 |. From our numerical results, we observe that in general, |∆ 2 | > |∆ 1 |, causing a reversed symmetry in the spin structure in the lowest two valleys. We also observe that higher order terms contribute to the valley-mixing SOC, as well as intra-valley SOC in the heterostructure case. Hence, we have generalized Eq. (9) by labelling the order n of the interaction, so far considered to zero order: α (n) ∝ (kα (0) ) 2 /|ε n − ε 0 |. Note that the numerical results presented in this work correspond to α = n α (n) and β = n β (n) .

IV. QUALITATIVE PICTURE
Beyond the symmetry arguments, we can also analyze the functions α (N, E z , v) and β (E z , N, v) in qualitative terms. Let us first note that there are two distinct regimes for these functions considered in the E z − N plane. In the weak-field (WF), thick well regime, the wavefunction for the electrons in the lowest electric subband is spread throughout the well. The strong-field (SF), thin-well regime is reached when the wavefunction is confined near one interface and does not feel the other. In both the ideal and sandwich cases, the potential is rather flat in the interior of the well over a region whose extent ∼ N and the confinement comes from relatively sharp interfaces. In this case, placing the classical turning point in the middle of the well shows that the dividing line between the two regimes is described by m l is the longitudinal mass. We need to consider the two sides of this line separately. 1. N dependence of α (E z , N, v) . For α the parity of N is not important. Let us define the lowest order term in E z for α (E z , N, v) as α 1 (N, v) E z (k x σ y − k y σ x ) . At first sight, the Rashba effect appears to be a bulk effect and therefore we expect α 1 (N, v) to be independent of N. However, the Ehrenfest theorem implies that the expectation value of E z , which is proportional to the mean force, must vanish for any wavefunction bound in the z-direction. Thus in a continuum effective mass approximation the lowest-order term must vanish even though it is allowed by symmetry. Only when we put in interface effects and other atomic-scale effects will this term emerge. We shall assume that the extent of the interface in the z-direction is independent of N. If this is the case, then the probability to find the electron at the interface in the WF regime is ∼ 1/N, and we may expect α (E z , N, v) to be a decreasing function of N in the WF regime. In the SF regime (large N ) α becomes independent of N for fixed E z since we only add layers that are unoccupied.
2. E z dependence of α (E z , N, v) . We have seen that at small E z (WF) the dependence on E z is linear. At large E z for fixed N (SF) the wavefunction is increasingly squeezed onto the interface and we may expect some continued increase in α. Hence, we can set: (a) N odd. This is the only case for which β (E z = 0, N, v) = 0. In the WF regime this fieldindependent term may be considered as a perturbation in 1/N, since it is zero for even N and the adding of an additional layer to make N odd is the same as adding a term to the Hamiltonian whose matrix elements vanish as 1/N. So we expect an initial decrease in the term as a function of N. Again, α should approach a constant at large N and fixed E z for the same reasons as in 1.
(b) N even. β (E z = 0, N, v) = 0. In the WF regime the field-independent term should converge to the result for even N as N increases, since they differ by terms of order 1/N. The same holds for the SF regime.
(a) N odd. There is a constant term but no strong dependence on E z in the WF regime. In the SF regime the wavefunction is strongly confined to the interface. If we consider just a two-layer interface, there is a very strong orthorhombic anisotropy: the [110] and [110] directions are different, since the nearest-neighbor bond is in one of the two directions. Hence β (E z , N, v) can be expected to be large, so β (E z , N, v) should increase strongly at large E z with fixed N.
(b) N even. In the WF regime the symmetry is very important and β (E z , N, v) is linear in E z . Again, in the SF regime we expect to converge to the odd N result.
Finally we consider the v dependence. The valley splitting ∆ v vanishes in the effective-mass continuum approximation -it is due to interface effects. For E z = 0, the eigenfunctions must be even in z: ψ + (r) = F (z) φ ( r) cos k 0 z or odd in z: ψ − ( r) = F (z) φ ( r) sin k 0 z, where F (z) is an even, slowly-varying envelope function and φ ( r) is an even (in z) Bloch function. k 0 is the wavector of the conduction-band minima. For a well with smooth surfaces such as we consider here, E v has the order of magnitude ∼ 1 − 10 meV and oscillates with thickness on the scale ∆N = π/2k 0 a z and is proportional to 1/N, as expected for an interface effect in the WF regime. In the SF regime, E z = 0 and the eigenstates are no longer even or odd. ∆ v saturates for large N at fixed E z , and its overall magnitude increases with E z , also as expected as the wavefunction is squeezed onto an interface.
The oscillations with N arise in the following way. Let V (z) be the confining potential and V (k) its Fourier transform: V (k) = exp (ikz) V (z) dz. If we apply lowest-order degenerate perturbation theory for states in the two valleys, we find ∆ v = d 3 r F 2 (z) φ 2 ( r) e 2ik0z V (z) ∼ |V (2k 0 )| . As we change N, V (z) has variation on the scale z = 4N a, the separation between the two interfaces. |V (2k 0 )| then has constructive interference when 2k 0 × 4 (∆N ) a = 2π or ∆N = π/4k 0 a. This ignores Umklapp, which will be present in the actual system, but is absent in the tightbinding approximation. It turns out that the dependence on the valley index can be quite dramatic. The valley states differ substantially right at the interface, where much of the spin-orbit effect arises. By the same token, we can expect the same oscillations with N that are seen in ∆ v to be present in ∆ SO , the spin-orbit energy splitting.
In actual heterostructures, the interfaces are not sharp. Ge is substituted for Si on randomly chosen sites, which will generally mean that the penetration length of the wavefunctions into the barriers varies randomly in x and y. All symmetries are violated by the disorder and it no longer makes clear sense to speak of even and odd N. Furthermore, the free electrons come from dopants that create an electric field in the structure, so that the E z = 0 limit is not accessible. It will probably be very difficult to observe the parity dependences that are predicted from symmetry arguments and also the oscillations. However, it may be possible to observe such effects in free-standing membranes.

V. RESULTS FOR IDEAL CASE
In this section we get the tight-binding results for the free-standing layer. Our approach to determining α and β will be to compute ∆ SO along [110] and [110] using NEMO-3D for free-standing layers with varying thickness N and applied electric field E z . To discriminate the Rashba and Dresselhaus contributions, we also compute the expectation value of σ x ′ and σ y ′ , which determines φ 1 of Eq. (6) and thus, the sign of (β − α).  In Fig. 2 we show the four energy eigenvalues of the lowest electric subband as a function of wavevector in the [110] direction for a 8 nm QW, which corresponds to N = 60, and E z = 0. The valley splitting is about 9 meV, which is much larger than the spin-orbit splitting for all k. It is seen that ∆ SO is linear in k at small k. Extracting ∆ 2 and ∆ 1 (see diagram in the inset of Fig. 2) we can then calculate (β + α) and (β − α) from the slopes of these lines, which allows us to determine α i (E z , N ) and β i (E z , N ) (see Eq. (9)).
We consider first an ideal Si QW, were the symmetry arguments can be directly applied. We impose hard wall boundary conditions at both interfaces and perform tight-binding calculations to obtain the eigenvalues and eigenvectors. As mentioned above, we expect the Rashba contribution to be linear in E z , α 0 (E z , N ) ≈ α 1 0 E z , and insensitive to N in the SF regime. The inset of Fig. 3 shows that this is the case: The intra-valley contribution of the Rashba SOC, α 1 0 , is plotted as a function of N . For large N , α 0 converges to a value and is insensitive to any changes of N , as corresponds to the SF regime. From this data, we find α 1 0 ≃ 2 × 10 −5 nm 2 in the SF limit. The inset shows α 0 as a function of E z : for the regime of fields considered here, it appears to be linear. For the Dresselhaus contribution, we consider first the N even case. We do not observe any SOC related terms at E z = 0, as expected: We recall that any term of the form k j σ i is violated under D 2h symmetry operations. We expect the terms to be linear to lowest order of E z , which is indeed the case: β 0(z) ≃ β 1 0(z) E z . The results for the intra-valley are shown in Fig. 4. The SF value for the intra-valley is β 1 0 ≃8×10 −5 nm −2 . We note an abrupt change in β 0 accompanied by a parity flip (depicted in the inset of Fig. 4), an event that has already been noted in literature 24,25,26,27 . This reveals that β 0 is more sensitive to higher energy contributions than α 0 .
The inter-valley parameters are plotted next. We observe again a linear behavior with electric field for both terms, as expected from our qualitative arguments.
Recall that under the D 2h operations, both the Rashba (for any N ) and the Dresselhaus SOC (for even N ) transform in a similar manner, hence a similar behavior is expected with E z . We observe that both contributions exhibit oscillations as a function of N . These oscillations are related to the valley-splitting oscillations, already observed in literature 15,21,23,27,28 . For large N , the interaction converges to the SF limit value: These results lead to two significant conclusions: (i) the Rashba inter-and intra-valley contributions are of the same order of magnitude, which would indicate that the bulk part of α predominates over the interface part in the SF regime. On the contrary, β 1 z is more than one order of magnitude larger than β 0 z , indicating that it is a pure interface effect, also sensitive to higher order energy levels. (ii) We also find from Figs. 3-5 that the dividing line between the SF and WF regime is at E z N 3 ≃ 10 11 .
Next, we consider the Dresselhaus contribution for N odd. Fig. 6 shows β 0 as a function of E z for different N . We observe that for large N , the parity effect is less apparent, and β 0 becomes independent of N, as predicted from our qualitative arguments. In the WF regime, β z is not as sensitive as β 0 to electric fields. We note that β 0 (E z = 0) is non-zero, as shown in the inset of Fig. 6. It presents strong oscillations with N in the WF regime, due to mixing with higher energy states, and vanishes in the SF regime, as expected in the bulk limit for silicon. The overall β 0 (E z = 0) ∼ 1/N dependence is also evident in this curve. The inter-valley Dresselhaus coupling constant for N odd is shown in Fig. 7. For large N , β z becomes linear with E z , as it corresponds to the bulk limit. The behavior is quite similar to β 0 .

VI. RESULTS FOR HETEROSTRUCTURE CASE
Having verified that the results are reasonable overall, we redo the calculations for a more realistic model of an actual heterostructure.
First, we present results on a Si x Ge 1−x /Si/Si x Ge 1−x membrane with x = 0.5. The Si layer is surrounded by 28 layers of Si x Ge 1−x on both sides. This is sufficient to avoid surface effects and to confine the wavefunction in the Si QW for the electric fields presented here. Biaxial (shear) strain is minimized using VFF, which includes precise values of the elastic constants c i , as defined elsewhere 29 . The unit cell is chosen to have 24 atomic layers along the [010] and [100] directions, for which realizations of the substitutional disorder are averaged over, as noted in our numerical data. These calculations are much more time-consuming than those for the ideal layer, so fewer results are presented. We recall that calculations of this kind have not been done previously. Earlier work used the virtual-crystal approximation for the outer layers, which artificially preserves the symmetry. The reduction in symmetry can only increase the number of possible terms in the Hamiltonian, so not only the Rashba and Dresselhaus terms exist, but in principle all terms k i σ j can exist -strictly speaking, even k x and k y are no longer good quantum numbers. However, we shall take advantage of the approximate symmetry to present the results in the same way.  Fig. 8 shows the parameters α 0 (a) and β 0 (b) as a function of electric field, E z . We note that the dramatic dependence on the parity of N is no longer present, as the disorder destroys their distinct symmetry properties. We also note that α is always non-zero, even for E z = 0. The intra-valley α 0 is non-linear in the WF regime, reaching a linear-in-E z value comparable to the one of Si QW in the SF regime (open triangles), α 1 0 ∼ 1.9×10 −5 nm 2 . This is consistent with the value obtained above for the Si QW. Fig. 8 (b) shows an overall 1/N dependency of β i for E z = 0 , consistent with the previous section results. For large N we observe that β 0 is linear with electric fields. We fit QW in the 10-30nm range and find β 1 0 ≃ 38×10 −5 nm 2 . We also find, for the SF regime, β 1 z ≃ 58×10 −5 nm 2 , with β i = β 1 i E z . We have also studied the eigenvectors for each case and find frequent parity-flips (see inset of Fig. 4) by varying N or E z . As a consequence, both β and α depend on N even in the SF regime.
Finally, we include the most experimentally relevant case: a Si x Ge 1−x /Si/Si x Ge 1−x grown on a Si 1−x Ge x substrate, with x = 0.3. The in-plane lattice constant a is now relaxed to the SiGe. We apply fixed boundary conditions for the lattice constant to the SiGe value at the bottom of the QW, and allow NEMO-3D to minimize the strain energy by varying the lattice constant along the z-axis.    9 shows the linear-in-k SOC intra-valley contributions. The splittings are linear in E z in the SF limit, although parity flips and higher energy levels results in small non-linearities. We note, however, that the zero field shows a visible splitting, even for wide QWs (N ≃140). We observe that α 0 is much smaller than in the previous samples, where silicon was the dominant component. For a typical Si/SiGe heterostructure with n s = 4 × 10 11 cm −230 we have k F ≃ 0.16nm −1 , at which the Dresselhaus-induced SOC splitting is ∆ β ≃ 1.25 µeV, whereas the Rashba splitting is only 0.02 µeV.
We show in table I the linear SOC coefficients for a 20nm Si QW (N ≃140) in the three different structures considered in this work, for a typical electric field E z = 10 7 V/m. α 1 i and β 1 i are in units of 10 −5 nm 2 . Note that the Dresselhaus is dominant in all cases.

VII. CONCLUSIONS
We have been able to extract the anisotropic-in-k splittings due to SOC using a sp 3 d 5 s * tight binding model capable to take into account the interface effects to atomic scale. For a Si QW, we distinguish the N even and odd cases, since symmetry operations over the Dresselhaustype terms are fundamentally different: while for N odd it appears at zero-order in E z , for N even is linear in E z , to lowest order (as Rashba-type terms are). We have extracted the linear-in-E z and linear-in-k parameters for the N even case.
We also distinguish two regimes of operation for typical wells: in the weak field regime, the splittings vary strongly as a function of N . The intra-valley mixing components show roughly a 1/N behavior, whereas the valley-mixing ones present also oscillations. On the contrary, the splittings do not change with N in the strong field limit. In the case of Dresselhaus, the dependency on N is more prominent. Together with some oscillations observed for the intra-valley mixing, this would reveal that higher-energy states are also very important. We also observe a reverse spin structure in the spin-split valleys, a direct consequence of the inter-valley splitting being larger than the intra-valley.
We find consistent results for Si QWs formed in SiGe heterostructures. We have also studied the lowest subband eigenstates, and found frequent parity flips by varying N or E z , suggesting a sample-dependent SOC. In accordance with ref. 15 , we find that the energy splittings due to Dresselhaus are in general larger than Rashba type ones: the Dresselhaus parameter β is almost one order of magnitude larger than the Rashba α for the ideal (pure Si QW) case, and roughly two orders of magnitude larger for the heterostructure (Si on SiGe substrate) case. New experimental data fits show indeed that the Dresselhaus term is larger than the Rashba one 31 . We recall, however, that throughout this paper the numerical data have been obtained for QWs with flat interfaces along a main crystallographic axis. More realistic samples include a small tilted angle with respect to a high symmetry axis. The numerical data are beyond the scope of this publication, however, we expect the interface induced Dresselhaus parameter β to become smaller in tilted samples in a similar manner as valley splitting does 32 . We expect, nevertheless, that α would remain of the same order of magnitude even in realistic tilted samples. Simulations carried out on a thin membrane confirm this hypothesis, although a more extensive study would quantitatively determine the SOC parameters.