Electromagnetic control of valley splitting in ideal and disordered Si quantum dots

In silicon spin qubits, the valley splitting must be tuned far away from the qubit Zeeman splitting to prevent fast qubit relaxation. In this work, we study in detail how the valley splitting depends on the electric and magnetic fields as well as the quantum dot geometry for both ideal and disordered Si/SiGe interfaces. We theoretically model a realistic electrostatically defined quantum dot and find the exact ground and excited states for the out-of-plane electron motion. This enables us to find the electron envelope function and its dependence on the electric and magnetic fields. For a quantum dot with an ideal interface, the slight cyclotron motion of electrons driven by an in-plane magnetic field slightly increases the valley splitting. Importantly, our modeling makes it possible to analyze the effect of arbitrary configurations of interface disorders. In agreement with previous studies, we show that interface steps can significantly reduce the valley splitting. Interestingly, depending on where the interface steps are located, the magnetic field can increase or further suppress the valley splitting. Moreover, the valley splitting can scale linearly or, in the presence of interface steps, non-linearly with the electric field.


I. INTRODUCTION
The spin of isolated electrons trapped in silicon-based heterostructures is very promising for building high performance and scalable qubits [1]. The long relaxation time [2][3][4] and dephasing time [5][6][7] that are achieved in these qubits are due to the weak spin-orbit interaction and nuclear zero-spin isotopes. Strong coherent coupling between Si spin qubits and photons using superconducting resonators has been realized [8,9] while the fidelities demonstrated for single and two-qubit gates are steadily improving [10][11][12][13][14]. Having mentioned all these advantages, we note that the nature of the degenerate conduction band minima, known as valleys, in bulk silicon poses a significant challenge for the operation and scalability of silicon spin qubits. It can be shown that a combination of biaxial strain as well as the sharp interface potential lifts the valley degeneracy in Si heterostructures, and gives rise to two low-lying states [1]. In general, a qubit performs well only when the qubit energy splitting is well-separated from any other energy scale in the environment. In Si spin qubits, the spin couples to the valley degree of freedom due to interface-induced spinorbit interaction [15][16][17]. If the valley splitting becomes equal to the qubit Zeeman splitting, a condition known as spin-valley hotspot, the valley-spin mixing for the qubit excited state reaches its maximum and gives rise to a very fast qubit relaxation via electron-phonon interaction [3,18]. It is, therefore, of crucial importance to understand how the valley splitting behaves as a function of parameters of the system; namely, electric and magnetic fields, the quantum dot geometry and the roughness at the Si/barrier interface.
In studying the valley splitting, a suitable starting point is the effective mass theory that can be used to ob-tain the electron envelope function. This envelope function, in turn, depends on the above-mentioned system parameters, and the valley splitting can be deduced from it. As we review in Section II A, in the absence of a magnetic field, the Hamiltonian describing the full envelope function is separable. While the in-plane envelope function is trivially given by the harmonic-oscillator wave function (due to the in-plane parabolic confinement), to our knowledge, the (ground state of the) out-of-plane envelope function has only been studied and approximated via variational methods [19][20][21][22] or by setting the barrier potential to infinity [30]. However, the assumptions involved in these methods render them less accurate for higher electric fields. In this paper, we model a realistic potential profile for a SiGe/Si/SiGe quantum dot by taking into account both Si/SiGe interfaces as well as an interface between SiGe and the insulating layer hosting the gate electrodes. Within this model, we then find the exact solution for the ground state as well as excited state envelope functions for the out-of-plane electron motion.
In the presence of an in-plane magnetic field, a cyclotron motion of electrons takes place which tends to increase the electron probability amplitude at the Si/SiGe interface [23]. This effect can, in turn, modify and increase the valley splitting. The magnetic field couples in-plane to out-of-plane degrees of freedom and thus prevents us from finding the exact solution for the electron envelope function. Using the exact excited states for the out-of-plane envelope function, we find the full envelope function in the presence of a magnetic field by applying perturbation theory. We show that an in-plane magnetic field indeed slightly increases the valley splitting; up to a few Tesla, the valley splitting increases quadratically with the magnetic field. Besides this, for a quantum dot with an ideal interface, i.e, no miscuts and steps at the interface, we find that the dominant contribution to the valley splitting scales linearly with the electric field, F z .
During the experimental process of fabricating silicon heterostructures, the formation of steps and miscuts at the Si/SiGe seems to be inevitable [24,25]. It has been shown that the presence of interface steps can severely suppress the valley splitting [22,26,27]. Here we again use the exact excited states for the out-of-plane envelope function in order to perturbatively treat the effects of interface steps to the envelope function. We argue that our modeling is applicable to any configuration for the interface disorder. We first study how the interface steps suppress the valley splitting in the absence of a magnetic field. We show that the valley splitting of a disordered quantum dot can scale either sub-linearly, linearly or super-linearly with the electric field, depending on the steps' configuration. We then consider the effects of an in-plane magnetic field. While it has been speculated that the magnetic field can increase the valley splitting in the presence of interface steps [26,27], interestingly, we find that the magnetic field can both increase or further suppress the valley splitting depending on the locations of the steps. This paper is structured as follows: In Section II A we present our model and find the exact solution for the outof-plane electron motion. In Section II B we obtain the envelope function in the presence of an in-plane magnetic field for a quantum dot with an ideal Si/SiGe interface.
In Section II C we extend our model to include the interface disorders and derive the envelope function for a certain configuration of steps. In Section III we build on our findings for the envelope function to obtain and discuss the valley splitting; in Secs. III A and III B, we study how the valley splitting of an ideal quantum dot depends on the electric and magnetic field field. In III C, we consider interface disorder and calculate the valley splitting and its phase depending on the location of the step. We then investigate the role of the electric and magnetic fields in modifying the valley splitting. In Section IV we summarize and conclude the paper. The Appendices contain further details of our analysis.

II. MODEL
We consider a SiGe/Si/SiGe heterostructure grown along theẑ direction ([001]) while the silicon layer is between −d t < z < 0. Panel (a) of Figure 1 shows a schematic cross-section of the layer structure of the system. An electric field is applied alongẑ via the gates, and we consider both interfaces between Si and SiGe at z = 0 as well as z = −d t . The energy offset between the minima of the conduction band in Si and SiGe is given by U 0 = 150 meV. Moreover, we also consider the interface between the top SiGe barrier and the insulating layer that hosts the electric gates. Inside the insulating layer, we take U ∞ = ∞ which indicates the envelope function does not leak into that region. Panel (b) of Figure 1 shows the full potential alongẑ. Moreover, here we assume a generally elliptical quantum dot with harmonic in-plane confinement. We denote the radius of the quantum dot alongx by x 0 and the radius alongŷ by y 0 . (a) Schematic layered structure of a single SiGe/Si/SiGe quantum dot. The dark gray area is an insulating layer hosting the electric gates. The ±V gates are used to trap and confine a single electron in the silicon layer. (b) The electrostatic potential profile along the growth directionẑ. Fz is the out-of-plane electric field generated by the gates, U0 is the potential barrier of the SiGe layers, U∞ is the infinite potential barrier due to the insulating layer. The potential energy U (z) is defined in Eq. (2).

A. Exact envelope function in absence of a magnetic field
Within the effective mass theory and in the absence of a magnetic field, the Hamiltonian describing the envelope function reads, and its associated energy scale, in order to piecewise expressing Eq. (3) as, Herez = z/z 0 ,˜ z,n = E z,n / 0 andŨ 0 = U 0 / 0 are the normalized length, eigenenergy, and potential. The above equation, at each interval, has generally two linearly-independent solutions known as Airy functions of the first and second kind, Ai and Bi [28]. We thus find the exact solution for ψ z,n : Note that the Bi function is omitted from the solution for z < −d t . This is based on the physical ground that Bi does not give rise to a decaying behaviour inside the extended barrier layer.
In order to find the eigenenergies and determine the coefficients involved in the envelope function Eq. (7), we note that ψ z,n and its first derivative must be continuous at the boundaries between Si and SiGe, i.e. at z = −d t and z = 0. Moreover, since there is no leakage to the insulating layer, the envelope function must vanish at z = d b . By imposing these boundary conditions, we obtain the equation below from which we can numerically find all possible eigenenergies, (9) with the definitions, and, where Once the (normalized) eigenenergy˜ z,n is found, we use it to calculate the coefficients c 1 to c 5 . The coefficient N 0 is found by using the normalization of the envelope function. We note that by solving Eq. (9), we also find a set of states where the envelope function is not localized in the Si quantum well but rather in the upper SiGe barrier underneath the insulating layer. As we discuss it in Section III, the valley splitting is basically determined by the ground state localized in the Si quantum well. In the presence of a magnetic field or interface steps, we also need to take into account the excited states which have sizable overlap with the localized ground state in the Si quantum well; see Eqs. (26), (27) and (33). As such, the states that are localized underneath the insulating layer do not contribute to the behavior of the valley splitting, and we neglect them in this paper.
For the ground state of the electron motion alongẑ, we can simplify the analysis presented above and find analytic relations. As we show in Appendix A, the (normalized) ground state energy in the regime of a deep quantum well,Ũ 0 1, can be expressed up to the leading order as˜ where −r 0 −2.338 is the smallest root (in absolute value) of the Ai function. The normalized envelope function in this case is approximated by, Figure 2 we show the obtained energies for the ground state as well as first few excited states as a function of the applied electric field. In Figure 3 we also show the probability density |ψ z,n | 2 for the ground state and first two excited states for two different electric fields. In both figures, a comparison between the numerics and the analytical relations Eqs. (14) and (15) for the ground state shows a very good agreement.
We note here that the interface-induced spin-orbit interaction is neglected in our model. Consideration of this effect has been shown to be essential for explaining the valley-dependent g-factor in silicon quantum dots [17,30]. However, as noted in Ref. [30], the matrix elements involved in the spin-orbit interaction are much smaller than the valley splitting matrix element. This justifies our omission of the interface-induced spin-orbit interaction. As we will show in the next sections, the information stored in the excited states ψ z,n≥1 enables us to obtain the full envelope function ψ xyz,0 in a finite magnetic field, and also makes it possible to study realistic cases where there are steps and miscuts at the Si/SiGe interface.
B. Envelope function in the presence of an in-plane magnetic field with ideal Si/SiGe interface Let us now consider a quantum dot with an ideally flat Si/SiGe interface in the presence of an in-plane magnetic field B = (B x , B y , 0). We use a gauge for which the vector potential becomes A = (0, 0, yB x − xB y ). By substituting p z → p z − eA z (B) in Eq. (1), we arrive at the following form for the Hamiltonian describing the envelope function, where we start from the separable Hamiltonian, and treat the field-induced couplings as a perturbation, We note that the confinement frequencies and lengths alongx andŷ are modified by the magnetic field, In order to obtain the envelope function from Eq. (16), we treat H 0 (B) exactly and apply perturbation theory in n αn (10 −4 T −1 ) βn (10 −4 T −1 ) H pert (B). The ground state up to the first order perturbation then reads, where, and, where the number of relevant bound excited states in the vertical direction for F z = 15 MV/m is found to be n max = 3, see Fig. 2. Here we defined the coefficients, We numerically calculate α n and β n using the excited states ψ z,n obtained in Section II A. For a circular dot we obtain α n = β n . For an elliptical dot with realistic parameters, these coefficients remain close to each other since the confinement alongẑ in quantum dots is always stronger than the in-plane confinements. Table I shows an example for the values of α n and β n . With the set of parameters used in Table I, we find (25) are one order of magnitude larger than x 0 y 0 η. In Figure 4 we show the probability density in the x − z plane in leading order, |ψ xyz,0 (B)| 2 , as well as the first order correction, |ψ (1) xyz,0 (B)| 2 .

C. Envelope function with disordered Si/SiGe interface
So far, we have studied structures where the interface between Si and SiGe is perfectly flat and is located at xyz,0 (B)(x, y = 0, z)| 2 (1/nm 3 ) due to an in-plane magnetic field; see Eq. (25). The parameters used are the same as given in the caption of Table I. The dashed line in both panels marks the ideally flat Si/SiGe interface. z = −d t and z = 0. However, during the experimental fabrication of Si qubit nanostructures, the formation of miscuts and steps at the interfaces is highly probable. Such uncontrolled disorder can modify the valley splitting and its phase, and is considered to be the main reason that makes the valley splitting a device-dependent quantity. In Ref. [22], several configurations for the steps at the interface are considered, and in each case, the envelope function is formed from a variational ansatz that uses a smooth interpolation between the envelope functions far from the step position (i.e. envelope functions for perfect interface).
In this section, we extend our model to include stairlike interface steps parallel to theŷ axis at the upper Si/SiGe interface as depicted in Figure 6. Our main objective here is to study how these miscuts influence the quantum dot envelope wavefunction. Note that disorder could also be present at the lower SiGe/Si interface. However, since the amplitude of the envelope function is small at the lower interface, the effects of possible disorder is negligible. We also point out that other disorder configurations at the upper Si/SiGe interface can be analyzed using the same approach.
In silicon, the thickness of each atomic layer is a 0 /4 where a 0 = 0.543 nm denotes the lattice constant. This indicates that the change in the interface position due to a few miscuts is much smaller than the total thickness of the envelope function alongẑ and enables us to use perturbation theory in order to obtain the electron envelope function. We take the z-position of the interface layer that contains the quantum dot center as the reference for the position of the barrier interface (e.g., the layer within [x 2 , x 3 ] in Fig. 6), and take any change to the interface position due to the miscuts as a perturbation. We describe the disordered SiGe/Si/SiGe interface with the step potential (29) where, The Hamiltonian describing the envelope function with disordered interface at finite in-plane magnetic field can again be written in the form of Eq. (16) where, in this case, Eq. (30) is added to the perturbative part of the Hamiltonian, Eq. (18). The ground-state envelope function then reads up to the second -order perturbation with respect to the interface disorders, xyz,0 (B)  Table I. We also find γ 1,0 = 0.3784 and γ 2,0 = 0.0663; see Eq. (38).
where N 0 is a normalization constant and D (1) xyz,0 is the first-order correction due to the interface disorder that amounts to for which the coefficients shall be calculated numerically. Table II shows examples for the values of γ m,n . Since the out-of-plane confinement is much stronger than the in-plane confinement, the largest contribution comes from m = 1 and n = 0. Moreover, we observe that by taking up to 4 excited states ψ x,m , the values of γ m,n substantially decay. As such, we can set m max = 4 as a cutoff in the summation in Eq. (32). For the second-order correction due to the interface disorder, we only keep the leading-order terms to arrive at (see Apendix C for more detail), where the perturbative coefficients c 1 and c 2 are given by, We now take, FIG. 6. The electron probability density in the x − z plane, that has the same functional form as D (1) xyz,0 given by Eq. (32) in which the perturbative coefficients become, otherwise . (38) In Figure 6, we show the electron probability density in the x − z plane in the presence of interface steps. The asymmetry around x = 0 in this case is due to the change of the quantum-dot thickness due to the interface disorder. Since in our model the Si quantum well is thicker for x > 0, the peak of the probability density is also shifted towards x > 0.
In the next section, we use the envelope functions we found in this section to study and discuss how the valley splitting of a quantum dot depends on the electric and magnetic fields for an ideally flat as well as disordered Si/SiGe interfaces.

III. DISCUSSION
Within the effective mass theory, the two low-lying valley components of the quantum dot can be written as Here k 0 0.85(2π/a 0 ) describes the Bloch wave vector of the conduction band minima and u ±z (r) are the periodic parts of the Bloch functions for the ±z valleys in silicon.
We express these functions by a plane wave expansion, for which G = (G x , G y , G z ) is the reciprocal lattice vector. The coefficients in this expansion for the two valleys are related via the time-reversal symmetry relation C − (G) = C * + (−G). The wave vectors and their corresponding coefficients C + (G) for Si are studied in Ref. [20].
The valley-orbit coupling is given by, Note that eF z z = 0z , and given o U 0 for all practical values of F z (see Eq. (50)), the valley-orbit coupling is strongly dominated by the matrix element of the interface potential U (z). Indeed, the role of the electric field is to control and shape ψ z , and the contributions from the matrix element of −eF z z in the valley-orbit coupling can be neglected [20,29]. The valley splitting is found from the above equation by E vs = 2|∆ vo | and the valley phase can be found by A. Electrical dependence of the valley splitting for an ideal quantum dot In this section, we consider a quantum dot with an ideal interface in the absence of a magnetic field and use the results of Section II A to find the electrical and interface-potential dependence of the valley splitting and the valley phase. As depicted in Figure 3, since the electric field pushes the envelope function towards the upper SiGe barrier, the probability density of the ground state at the lower SiGe/Si interface is negligible. Therefore we take U (z) = U 0 θ(z) and find that the valley-orbit coupling at zero magnetic field for an ideal quantum dot becomes, To carry on, we note that the terms with G 1 = G 2 would lead to fast oscillations in the integrand that average to zero. We therefore only consider terms with G 1 = G 2 and define, which we find to be C 0 = −0.2607 using Ref. [20].
We now take the integration by parts and find for an ideal quantum dot, Here ∆ int is the contribution that comes from the amplitude of ψ z,0 at the Si/SiGe interface: and ∆ t is the contribution that originates from the tail of ψ z,0 inside the barrier, In order to find analytical expressions for these contributions, we use Eqs. (14) and (15) and also use the expansions given by Eqs. (A8) (note that by using these equations we assume U 0 0 which is valid for all relevant valued for the electric field F z , see Eq. (50).) We finally arrive at the result and The last term in the square bracket is a number larger than 1 (having F z = 2 to 20 MV/m for a SiGe barrier, we find k 0 z 0 / Ũ 0 5). This indicates that |∆ int | is larger than |∆ t | (by nearly a factor of 6.) We conclude that the valley-orbit coupling (and hence also the valley splitting) within the leading order scales linearly with the electric field while it is independent of the interface potential (as long as U 0 0 ). This linear dependence is experimentally observed in Ref. [3] for a SiO 2 barrier (that has a much stronger interface potential U 0 = 3 eV compared with SiGe) and it is also predicted from a theory analysis assuming that the envelope function has zero amplitude inside the barrier [30].
In addition to the linear-in-electric-field term, here we find the valley-splitting also has a small nonlinear dependence on the electric field (note that 1/Ũ 0 ∝ F 2/3 z ). This non-linear contribution originates from the penetration of the envelope function into the barrier, and can be neglected so long asŨ 0 1; this holds provided, For a SiGe barrier with U 0 = 150 meV, the right side of the above inequality becomes ∼ 280 MV/m. This essentially means for all practically relevant values for F z , the valley splitting remains a linear function of the electric field. In Figure 7, we used Eqs. (48) and (49) and show |∆ int | and |∆ t | as a function of the electric field F z . The values that we show in the figure are in agreement with Refs. [21] and [29] where |∆ vo | ∼ 200 µeV is reported for F z = 15 MV/m. Finally, from Eq. (48) and (49) we find that the valley phase of a quantum dot with perfect interface up to the leading order only depends on the interface potential, that becomes φ v −0.44π for a SiGe barrier.

B. Magnetic dependence of the valley splitting for an ideal quantum dot
We now extend the results of the last section by including an in-plane magnetic field. As we have shown in Section II B, the electron envelope function Ψ xyz,0 in the presence of an in-plane magnetic field includes excited states of the out-of-plane motion ψ z,n , see Eq. (23) and (25). Compared to the ground state ψ z,0 , the excited states can have a larger amplitude at the Si/SiGe interface, and penetrate further to the barrier. As such, we generally expect that the valley splitting should increase in an in-plane magnetic field. In addition, as one can see from Figure 3, depending on the electric field, the excited states can have a sizeable amplitude at the lower SiGe/Si interface z = −d t . Therefore, we take into account both upper and lower interfaces and consider a barrier potential of the form U (z) = U 0 θ(z) + U 0 θ(−z − d t ). We use Eqs. (23) and (25) and find for the valley-orbit coupling, Using the excited states ψ z,n from Section II A, we numerically calculate the above integrals. The coefficients α n , β n and η are introduced in Section II B and as we explained there, up to a few Tesla, the terms containing α and β are strongly dominant over the correction containing η. Therefore, to the leading order, the magnetic contribution to the valley splitting scales quadratically with the magnetic field. Moreover, at finite magnetic fields, the valley splitting becomes dependent on ratio of the lateral confinement to the electrical confinement, x 0 /z 0 and y 0 /z 0 -see Eqs. (21) and (22).
Eq. (52) indicates that for an elliptical quantum dot at a fixed magnetic field, the valley splitting reaches its maximal value when the direction of the magnetic field is perpendicular to the axis with the larger radius. In the main plot of Figure 8, we show the valley splitting as a function of magnetic field at a fixed direction. In the inset plot of the figure, we show the valley splitting as a function of the direction of the field. We observe that the valley splitting for a quantum dot with ideal interface only slightly increases with the magnetic field. This also indicates that the dominant contribution to the valley splitting remains a linear function of the electric field F z at the finite magnetic fields.

C. Valley splitting of a quantum dot with disordered interface
We now consider a realistic quantum dot with miscuts and steps at the Si/SiGe interface and aim to study the valley splitting and its electromagnetic dependence. We take the interface potential given by Eq. (29) and use the resulting envelope function Eq. (31). We then find for the valley-orbit coupling, where ∆ ideal vo is the valley-orbit coupling for an ideal interface given by Eq. (52), ∆ s is the largest contribution originating from the interface disorders and reads, and ∆ (1) is a contribution that is first order with respect to D xyz,0 and reads, where, The last term, ∆ (2) , is a small and sub-leading contribution that is second-order with respect to D xyz,0 and ψ (1) xyz,0 . More details on this term can be found in the Appendix B.
Note that in general we have |∆ (2) | |∆ (1) | and |∆ (1) | < |∆ vo |, |∆ s |. However, as we discuss below, depending on the number and location of interface steps, |∆ vo + ∆ s | can become a small number. In this case, the contribution from ∆ (1) becomes (more) important in determining the valley splitting and its phase.

A single step at the interface
Let us now study the structure and effects of ∆ s and ∆ (1) . We begin by considering the simplest case; that is, when there is only a single step at the interface. We take the width of the step to be +a/4; the step potential is then obtained from Eq. (30) by taking x s4 → +∞ and x s1 , x s2 → −∞. We then take x s3 = x s ≥ 0 to be the position of the only interface step. From Eq. (54) for a quantum dot with a single interface step we obtain, In order to analyze ∆ 1s (1) , we note that a complete assessment of this term requires numerical calculations. However, we can obtain a rough estimation by only considering the largest contribution to Eq. (55), i.e. the term corresponding to m = 1 and n = 0. For further simplicity, we also drop the second-order correction to the envelope function, D (2) xyz,0 , so that we take γ 1,0 = γ 1,0 and N 0 = 1. We then arrive at the largest contribution to ∆ (1) due to m = 1 and n = 0, in which, ∆ 1s s (x s ) and ∆ 1s (1),{1,0} (x s ) are clearly out of phase with ∆ ideal vo . This can significantly modify and suppress the valley splitting. Note that ∆ 1s s (x s ) and ∆ 1s (1),{1,0} (x s ) are monotonically decreasing if the interface step is located further away from the quantum-dot center.
In the left panel of Figure 9 we present the normalized valley splitting as a function of the step location at B = 0. By neglecting ∆ (1) in the valley-orbit coupling (shown in the figure by the solid black line), we observe, as predicted, that the valley splitting is monotonically decreasing as the interface step moves closer to the quantum dot center. At x = 0, the valley splitting is suppressed by 75%. This amount of suppression is reported in Ref. [22] where a variational ansatz is used to approximate the envelope function in the presence of a single step. By taking into account the ∆ 1s (1),{1,0} (x s ) contribution (shown by the dashed-dot red line), we observe that the valley splitting only becomes further suppressed.
Remarkably, if we numerically calculate ∆ 1s (1) from Eq. (55) and take into account not only the dominant term corresponding to m = 1 and n = 0 but also the other terms as well, we observe that the valley splitting has in fact a non-monotonic behaviour as a function of the distance between the single interface step and the quantum dot center. In Appendix D, we discuss in more detail that this behaviour is due to the out-of-plane excited states ψ z,n≥1 .

Indeed, the terms ∆ ideal
vo , ∆ 1s s and ∆ 1s (1),{1,0} are originating from the ground state of the out-of-plane motion, ψ z,0 . As it is already shown, the closer the step to the quantum dot center, the further ∆ 1s s and ∆ 1s (1),{1,0} suppress ∆ ideal vo , until at some point, x s = x 1s m , the contributions due to the terms with n ≥ 1 in Eq. (55) reaches an equal footing as ∆ ideal vo + ∆ 1s s + ∆ 1s (1),{1,0} . If the interface step is placed any closer to the quantum dot center, x s < x 1s m , the latter sum continues to decrease whereas the contributions due to the terms with n ≥ 1 in Eq. (55) increase due to the increase of the coefficients γ m,n≥1 . Therefore, the valley splitting begins to rise.
In the left panel of the Figure 9 we have shown the valley phase as a function of the step location. Note that the valley phase, given by Eq. (42), is a π-periodic function defined within [−π/2, π/2]. Therefore, the sudden jump of the valley phase that occurs at x s ∼ 0.45x 0 can be removed by subtracting π from the values above the jump.
Let us now study how the valley splitting scales with the electric and magnetic fields. As we have shown in Sec. III A, the valley-orbit coupling for an ideal quantum dot ∆ 0 vo is a linear function of the electric field. As such, the ∆ 1s s term given by Eq. (57) is also linear in electric field. However, ∆ 1s (1) is a non-linear function with respect to the electric field. Given Eq. (59), the dominant coefficient γ 1,0 scales linearly with the electric field, and this gives rise to a quadratic scaling of ∆ 1s (1),{1,0} with the electric field (more accurately, if we keep c 1 in Eq. (38), it is easy to show that ∆ 1s (1),{1,0} acquires one more term that is cubic in the electric field.) Therefore, in general, the valley splitting for a quantum dot with a single interface is a nonlinear function of the electric field due to the ∆ (1) term as well as the normalization constant N 0 [1 + m,n γ 2 m,n ] −1/2 . In the upper panel of Figure 10, we show the valley splitting as a function of the electric field for several locations for the single interface step.
Since ∆ (1) grows with the electric field faster than a linear function, for x s > x 1s m , we expect the valley splitting should scale sub-linearly by the electric field. The further the step is located away from the center, the smaller ∆ (1) becomes so that the valley splitting approaches being a linear function of the electric field. For x s < x m , the valley splitting is mainly determined from the terms with n ≥ 1 in ∆ (1) . Therefore, numerical analysis is required to find the dependency of the valley splitting on the electric field.
In order to understand how the valley splitting changes with an in-plane magnetic field, we note that the confinement length x 0 (B y ) is reduced by the magnetic field. Using Eqs. (21), the effect of the magnetic field to the evolution of ∆ 1s s (x s ) and ∆ 1s (1),{1,0} (x s ) is equivalent as if x 0 = x 0 and x s is located from the center at a larger distance, Therefore, we expect that the magnetic field should always increase the valley splitting if x s > x m . Since the magnetic field effectively increases x s , the increase of the valley splitting by the magnetic field is larger at the step location where the slope of the curve given in the left panel of Figure 9 is steeper. For x s < x m , the magnetic field decreases ∆ 1s (1),{m,n≥1} (x s ) that controls the valley splitting. However, we observe that the valley splitting still slightly increases as a function of the magnetic field due to the increase of ∆ ideal vo + ∆ 1s s with the magnetic field. We display the valley splitting as a function of the in-plane magnetic field in the lower panel of Figure 10.

Two steps at the interface
We now consider a quantum dot with two interface steps having widths −a 0 /4 and a 0 /4. The step potential is then obtained from Eq. (30) by taking x s1 → −∞ and x s4 → +∞. For further simplification, we also assume that the two steps are placed symmetrically around the center so that we can write −x s2 = x s3 = x s . In this case, finding the contribution from ∆ (1) requires numerical analysis (even only for the term with m = 1 and n = 1 in Eq. (55).) However, we can still obtain qualitative understanding of the behaviour of the valley splitting by only considering the effect of ∆ s . In order to arrive to the extension of Eq. (54) for a quantum dot with two symmetrically located steps, we integrate over z by parts, similar to Sec. III A, and neglect the small integral containing ψ z,0 . With this, we arrive at If the steps are located at the center, x s = 0, we find from the above equation (taking F z = 15 MV/m) ∆ 2ss s (x s = 0) (1.93 − 0.18i)∆ int that is larger than ∆ ideal vo in Eq. (53) (note that ∆ ideal vo is largely determined by ∆ int ; see Eqs. (45) and (52) and Fig. 7.) When the steps are located away from the center, it reduces ∆ 2ss s (x s ) so that eventually at some point |∆ ideal vo | becomes larger than |∆ 2ss s |. As such, we can expect that the valley splitting of a quantum dot with two symmetrically located steps is a also non-monotonic function of the steps' location,  53) included. We observe that the effect of ∆ (1) is to further suppress the valley splitting as well as shift the step location where the valley splitting reaches its minimum. At this step location, x s = x 2ss m ∼ 0.3x 0 , the valley splitting is suppressed by more than 90%. The right panel of the figure shows how the valley phase is changed as a function of the step location. As mentioned before, the sudden jump can be removed by using the π-periodicity of the valley phase.
We now study the electromagnetic dependence of the valley splitting for a quantum dot with two symmetrically locates interface steps. Given Fig. 11, away from x 2ss m , the dominant contribution to the valley splitting is due to ∆ ideal vo and ∆ s . As such, we expect the scaling of the valley splitting with the electric field has to be approximately linear. This behaviour is shown in the panel (a) of the Figure 12. At the step locations close to x 2ss m , the contribution due to ∆ (1) becomes important so that we can expect a strong non-linear dependency on the electric field; this is shown in the panel (b) of the Figure. In the panel (c) of the figure we show how the valley splitting is changed by the magnetic field. For x s < x 2ss m , the valley splitting is mainly controlled by ∆ 2ss s . As such, since this term reduces by the magnetic field, the valley splitting is also decreasing as well. For x s > x 2ss m , ∆ ideal vo is the dominant contribution to the valley splitting; therefore, the magnetic field always increases the valley splitting by reducing ∆ 2ss s (and ∆ (1) .) Note that this effect is stronger for the step locations where the slope of the Eq. (30). In Figure 13, we present our result for the normalized valley splitting as a function of the position of each step. Note that due to the non-symmetric nature of the step potential, the valley splitting turns out not to be symmetric with respect to the location of the steps. We observe that the valley splitting can completely vanish in some specific configuration of the interface steps; this is also predicted in Ref. [22] using a simpler model. We have also studied the valley splitting for models including three and four interface steps, and we observed qualitatively similar behavior for the valley splitting as a function of electric and magnetic fields as presented in this section.

IV. SUMMARY AND CONCLUSIONS
To summarize, the valley splitting is one of the important properties for the silicon quantum dots that directly influences the lifetime and scalability of silicon spin qubits. As such, understanding the behaviour and tunability of the valley splitting is very important. In this work, we studied how the valley splitting responds to the electromagnetic field for both ideal and disordered quantum dots. We considered a realistic potential profile for a SiGe/Si/SiGe quantum dot by taking into account both lower and upper Si/SiGe interfaces as well the interface between upper SiGe layer and the insulating layer hosting the gate electrodes; see Fig. 1 and Eq. (2). While so far the out-of-plane electron motion has been studied by variational methods in a simpler potential model including only one Si/SiGe interface, we found the exact (within effective mass theory) envelope functions of the ground state as well as the excited states for the out-ofplane electron motion. This has enabled us to find the electron envelope function in the presence of finite magnetic field as well as interface disorder. In both cases, the envelope function reflects the coupling between inplane to out-of-plane degrees of freedom, see Eqs. (23) and (31). Our analysis enables us to obtain the coupling coefficients using perturbation theory, for arbitrary configurations for the interface disorder.
We showed that in an ideal quantum dot, the valley splitting, within the leading order, always scales linearly with the out-of-plane electric field; see Fig. 7. Moreover, the valley splitting slightly increases with an applied in-plane magnetic field due to the coupling to the out-of-plane excited states; see Figure 8. The presence of interface disorder can significantly modify and suppress the valley splitting. We considered a stair-like disordered interface and studied the suppression of valley splitting due to the interface miscuts.
For a quantum dot with a disordered interface, we found that depending on the number and locations of the interface steps, the valley splitting can scale nonlinearly with the electric field; see the upper panel of Figure 10 and panels (a) and (b) of Figure 12. If there is only one miscut at the interface, the magnetic field always increases the valley splitting, see the lower panel of Figure 10. However, for multiple interface miscuts, the magnetic field can both increase or even suppress the valley splitting, depending on the configuration of the miscuts, see panel (c) of Figure 12.
In the theory of spin relaxation induced by the valley coupling, one important set of quantities are the transition dipole matrix elements between the valley states [3,18]. For an ideal quantum dot, the envelope function has an in-plane mirror symmetry (Fig. 4). This, in turn, gives rise to vanishing of the in-plane dipole matrix elements. However, the presence of the interface disorder can break the in-plane mirror symmetry (Fig. 6). Our findings for the envelope function in the presence of disorder now enable the prediction of the dipole matrix elements as a function of the electromagnetic field. To our knowledge, the dipole matrix elements have always been taken as fitting parameters. Future work will be needed to further develop the theory of spin relaxation induced by the valley coupling based on a calculation of the transition dipole matrix elements and valley splitting in a magnetic field.
Inside the silicon layer where z ≤ 0, ψ z,0 is given by the Airy function of the first kind, Let us now consider the form of the eigenstate inside the SiGe barrier; the exact solution for the envelope function, up to prefactors, reads Ai(Ũ 0 −z −˜ z,0 ). However, in order to find an analytic relation for the ground state energy, we try to approximate the envelope function in the barrier. If there was no electric field inside the barrier (i.e. the −z term inζ 0 (z > 0) was absent), the eigenstates would have been proportional to exp(− Ũ 0 −˜ z,0 z). Due to the presence of the electric field, the potential barrier is reduced along z and the wavefunction can further penetrate into the barrier. To take this into account, we introduce a parameter λ into the exponent of the exponentially decaying wavefunction that allows further penetration into SiGe provided λ < 1. We thus approximate the wavefunction inside the barrier by, Using the continuity of the first derivative of ψ z,0 at the z = 0 interface, we can write the penetration parameter λ as, We now self-consistently determine the ground state energy by noting, From the above equation and using Eqs. (A3) and (A4) we finally arrive at, In order to find the solution of the Eq. (A6), we note that for the infinite potential wellŨ 0 = ∞, the ground state energy is determined by the smallest root (in absolute value) of the Airy function, −r 0 −2.3381. This suggests that for a finite but high potential wellŨ 0 1, the solution should remain close to −r 0 . As such, we consider a solution of the form, and expand Ai and Ai functions around −r 0 . We keep terms up to quadratic order in δ˜ to find, (A8b) By substituting Eqs. (A8) into Eq. (A6) and keeping terms up to quadratic order we findŨ 0 δ˜ 2 = 1 that gives δ˜ = ±1/ Ũ 0 .
In order to determine which sign is physically acceptable, we note from Eqs. (A3b) and (A4) that the wave function decays inside the barrier provided Ai (−˜ z,0 )/Ai(−˜ z,0 ) > 0. This is satisfied only if δ˜ has a negative sign. In this case, we also find λ < 1, as expected. Note that if we have kept the expansions in Eq. (A8) up to the cubic order, we would have found a correction to the˜ z,0 of orderŨ In order to simplify this, we use Eqs. (A7) and (A8) and find M 0 1/Ai (−r 0 ).
Appendix B: Valley splitting of a disordered quantum dot in magnetic field: Higher-order terms In this Appendix, we present the form of the subleading contribution ∆ (2) in Eq. (53). Using Eqs. (29) and (31)  The perturbative coefficients α, β, η and γ are given in Sections II B and II C, and using the excited states ψ z,n from Section II A, we can numerically evaluate the integrals.
Appendix C: Second-order correction to the envelope function due to the interface disorder Here we present the complete form of the second order correction due to the interface steps. According to perturbation theory, we have where we defined, ζ m,n m ,n = ψ x,m ψ z,n |U steps |ψ x,m ψ z,n E n,z − E n ,z + (m − m ) ω x , (C2) c m,n = ψ x,0 ψ z,0 |U steps |ψ x,0 ψ z,0 E 0,z − E n,z − m ω x .
In order to arrive to Eq. (34), we only keep the dominant terms; i.e., in the set of γ m,n , we keep γ 1,0 , and in the set of ζ m,n m ,n , we keep ζ 1,0 2,0 and ζ 1,0 0,0 . Finally, in the set of c m,n , we keep c 1,0 .
Appendix D: The effect of out-of-plane excited states in ∆ (1) In this appendix, we explain in more detail the effect of the excited states ψ z,n≥1 to the valley splitting for a quantum dot with single interface step. The step potential Eq. (30) in this case becomes, U 1s steps (x, z) = −U 0 θ(x − x s )θ(z)θ z − a 0 4 (D1) In Section III C, we found the contributions from ∆ 1s s and ∆ 1s (1),{1,0} are out-of-phase with ∆ ideal vo , therefore, these terms would monotonically suppress the valley splitting.
Here we show that once the out-of-plane excited states ψ z,n are taken into account, the real part of ∆ 1s (1) gives rise to a non-monotonic behaviour of the valley splitting as a function of the step location, as shown in left panel of Fig. 9.
To see this effect, let us integrate Eq. (56) over z by parts to find, f m,n iU 0 1 2k 0 ∞ xs ψ x,m ψ x,0 dx × ψ z,n (0)ψ z,0 (0) − e −ik0a0/2 ψ z,n ( a 0 4 )ψ z,0 ( a 0 4 ) , (D2) where we have neglected the small contributions containing integral over e −2ik0z ψ z,n ψ z,0 and e −2ik0z ψ z,n ψ z,0 . As we have shown in Sec. III A, the dominant contribution in ∆ ideal vo is given by ∆ int from Eq. (48) which is an imaginary quantity and is due to the amplitude of the envelope function at the Si/SiGe interface.
Given Eq. (D2) and Eq. (55), the imaginary part of ∆ 1s (1) is in opposite phase with ∆ int (note that in Eq. (D2) we have e −ik0a0/2 = −0.891 − 0.454i.) Therefore, the closer the step is located to the quantum dot center, the more the imaginary part of total valley-orbit coupling is suppressed. On the other hand, the real part of the valley-orbit coupling is increased when the step is closer to the center, and at some point it becomes the dominant contribution to the total valley-orbit coupling. We plot the imaginary and real parts of the valley orbit coupling in Fig. 14. The dashed-dot lines are obtained by using the analytic relations we obtained in Sec. III C for D 1s s and D 1s (1),{1,0} whereas the solid lines are found from numerically evaluating the valley-orbit coupling including all terms in ∆ 1s (1),{m,n} . We observe that the nonmonotonic behaviour of the valley splitting as a function of step location can be seen only by taking into account the out-of-plane excited states ψ z,n .