Quantitative modeling of spin relaxation in quantum dots

We use numerically exact diagonalization to calculate the spin-orbit and phonon-induced triplet-singlet relaxation rate in a two-electron quantum dot exposed to a tilted magnetic field. Our scheme includes a three-dimensional description of the quantum dot, the Rashba and the linear and cubic Dresselhaus spin-orbit coupling, the ellipticity of the quantum dot, and the full angular description of the magnetic field. We are able to find reasonable agreement with the experimental results of Meunier et al. [Phys. Rev. Lett. 98, 126601 (2007)] in terms of the singlet-triplet energy splitting and the spin relaxation rate, respectively. We analyze in detail the effects of the spin-orbit factors, magnetic-field angles, and the dimensionality, and discuss the origins of the remaining deviations from the experimental data.


I. INTRODUCTION
Quantum dots (QDs) are well-known examples of confined and quantized systems in semiconductor heterostructures. Among other potential applications in, e.g., quantum optics, QDs have been experimentally realized as controllable quantum bits. 1 On the road to quantum computing applications it is essential to understand the electronic and spin-dependent relaxation processes in QDs. 2,3 The coupling between quantized states and lattice vibrations in semiconductor structures has been examined already for decades. [4][5][6] In QDs the spin relaxation has attracted attention both theoretically 7-14 and experimentally. [15][16][17][18][19][20] Recently, Meunier et al. 15 measured the triplet-singlet relaxation rate in a two-electron (quasi) two-dimensional (2D) QD. They also derived a semiempirical model reproducing the general trends of the experimental data. However, a good agreement was found using a spin-orbit (SO) coupling that depended on the magnetic field, and that was significantly smaller than the SO strength reported in earlier studies. 21 In this work we describe a two-electron QD using both 2D and three-dimensional (3D) models, respectively, and include (i) the description of a tilted external magnetic field, (ii) the ellipticity of the QD, and (iii) both the Rashba and the linear and cubic Dresselhaus spin-orbit (SO) interaction. We point out that none of the previous theoretical approaches have aimed at such a complete description, which is shown to be important when comparing with experiment at the most detailed level. The two-electron states are calculated numerically exactly by diagonalizing the Hamiltonian, and the phonon-induced relaxation rate is obtained from Fermi's golden rule. We find a reasonable agreement with the experimental results, 15 provided that the 3D description and a moderate ellipticity of the QD are applied. The remaining discrepancy between experiment and theory can be caused by a more complex dot geometry of the experiment, and/or it indicates the limit of the widely used effective mass approximation for electrons confined by a harmonic potential.

II. MODEL AND METHOD
Our two-electron QD exposed to an external magnetic field B is described by the Hamiltonian where r i = (x i , y i , z i ) with i = 1, 2 are the coordinates of the two electrons. The single-electron Hamiltonian is given by a sum of the kinetic, external potential, and Zeeman terms as We consider magnetic fields tilted from the xy plane with an angle θ, and tilted azimuthally from the y axis with an angle φ. Thus the magnetic field has an expression The vector potential in the symmetric gauge, A = (B × r) /2, then becomes which approaches A = B 0 (−y, x)/2 in the 2D limit (z → 0). The harmonic confinement potential is asymmetric both laterally and vertically: the ellipticity in the xy plane can be tuned with the δ-parameter, and the inplane and off-plane confinement strengths are fixed to ω 0 = 2.5 meV and ω z = 11.9 meV, respectively. The material parameters are chosen to be those of GaAs, i.e., effective mass m * = 0.067m e , relative permittivity ǫ r = 12.4 and the gyromagnetic ratio g * = −0.44. The constant µ B = e/(2m e ) in Eq. (3) is the Bohr magneton.
A harmonic form has been shown to describe well the electron confinement in both lateral and vertical semiconductor (GaAs-type) quantum dots. 33 The model was established soon after the first Coulomb-blockade experiments, when the measured and calculated addition energies were compared. 34 More recently, the harmonic model has been shown to explicitly yield the measured singleeletron spectrum, 35 and also the spin-blockade oscillations for both lateral and vertical QD devices up to about 50 electrons. 36,37 The SO interaction in Eq. (1) is given by with the Rashba term and the Dresselhaus term 25 containing both linear and cubic contributions, The coupling strengths of the Rashba and Dresselhaus terms are set by parameters α and γ, respectively, and σ x and σ y are the Pauli matrices. We choose α = 0 and γ = 11.15 eVÅ 3 as our default values, but consider also different magnitudes to assess the sensitivity of the spin relaxation to the SO coupling. The measurements of Meunier et al. 15 and 2D calculations 8,23 suggest relatively small coupling constants. This is in agreement with the experiment of Zumbühl et al. 24 reporting γ = 9 eVÅ 3 which is notably smaller than the commonly used value of 27.5 eVÅ 3 . The total wave function is a product of the xy (planar) component ψ (r 1 , r 2 ), obtained from exact diagonalization of the Hamiltonian in Eq. (1) in a (large) basis of spin-symmetrized products of Cartesian singleelectron harmonic-oscillator eigenstates, 22,26 and a simple z (perpendicular) component ψ z (z 1 , z 2 ) incorporating the material thickness. The z component of the wave function is conveniently modeled in terms of the ground state wave function of a harmonic potential, with confinement strength ω z defining an effective confinement length, L z = 4 /(m * ω z ), in the z direction.
The spectrum is obtained by diagonalizing the Hamiltonian in a basis of spin symmetrized two-electron harmonic-oscillator states. Basis sizes |n x , n y , n z up to n max = 6 for the (x, y) coordinates and n z = 0 for the vertical coordinate are included. Analytical matrix elements are then obtained for all terms except for the electron-electron interaction in Eq. (1).
It was shown by Popsueva et al. 22 that, in the 2D case with ω x = ω y = ω 0 , all matrix elements can be evaluated analytically. In particular, for any double pairs of basis functions (i, j) the four-dimensional integrals over the electron-electron interaction can be expressed as where The variable s is here defined from the Bethe integral, When this expression is introduced in calculating the matrix element between any pair of basis functions (i, j), the expression above is derived with s = s 2 x + s 2 y . When the integrals have been analytically evaluated once, the matrix elements for any confinement ω 0 are obtained by a simple multiplication.
In 3D the electron-electron repulsion amounts to more complicated six-dimensional integrals. When extending with harmonic confinement also in the z direction a similar scaling relation can be obtained, i.e., Here erf(x) = 2/ √ π x 0 exp(−t 2 )dt is the error function. The matrix elements for the electron-electron interaction can thus be obtained for any ω 0 by scaling as in 2D when they have been calculated for a given value of the vertical confinement. We note also that the formula above provides a continuous route from 3D to 2D, since For 2D calculations the finite thickness can be consistently taken into account by replacing Eq. (8) by Eq. (11) when calculating the interaction terms. The triplet-singlet relaxation rate is calculated using Fermi's golden rule, which -in the case of phonon emission -may be written as, 27 where T and S refer to the initial singlet and final triplet states, respectively, V is the normalization volume, ∆E is the singlet-triplet energy splitting, q = (q , q z ) is the momentum of the released phonon, and H ph = i=1,2 e −i(q ·ri+qzzi) is the phonon coupling. Linear dispersion relations have been adopted, i.e., we set ω S,T ≡ ∆E/ = c l,t q. Here c l = 4720 m/s (c t = 3340 m/s) denotes the longitudinal (transverse) speed of sound. 15 The electrons couple to longitudinal acoustic phonons through the deformation potential coupling 28,29 as well as the piezoelectric coupling, whereas coupling to transverse phonons only takes place through the latter one, (16) In these expressions ρ = 5300 kg/m 3 is the GaAs mass density, Ξ d = 6.7 eV is the deformation potential constant, and h 14 = 1.4 × 10 9 V/m is the piezoelectric constant. There are two transverse phonon modes, and hence M 3 (q) is considered twice when computing the rate.
The matrix element S| H ph |T in Eq. (13) separates into a product of a planar (xy) and a perpendicular (z) component, i.e., Here |ψ (S) and |ψ (T ) are the planar two-electron singlet and triplet states, respectively. We have assumed that the z component of the wave function is frozen, i.e., no excitations are considered in that direction. The expansion coefficients a k,l are obtained from the diagonalization of the Hamiltonian in Eq. (1). n k and m k are the single-particle harmonic quantum numbers of one of the electrons in the x and y direction, respectively. The factor of two in Eq. (17) arises from the fact that the total wave function of each state is antisymmetric. The singleparticle matrix elements Φ ni,n i ′ = φ ni |e −iqxx |φ n i ′ are essentially associated Laguerre polynomials in q 2 multiplied by a damping exponential function, Hereq x = q x /(2m * ω 0 ) 1/2 , N = min(n, n ′ ), M = (|n − n ′ | − η)/2, and η = 1 if n + n ′ is odd and 0 otherwise (analogously in y) 32 . The phonon induced triplet-singlet relaxation rate is now calculated by inserting Eqs. (14)- (17) into Eq. (13) and performing the integration over the phonon momentum q.

III. RESULTS
A. Singlet-triplet energy splitting The tilting angle θ has a relatively large effect on the energies, so that the closing point of the singlet-triplet energy gap ∆E moves to larger magnetic fields when the angle is increased. Figure 1(b) shows the gap as a function of B in comparison with the experimental values (circles). First, the downward "bending" of ∆E at small fields is due to the ellipticity; in a circular QD the gap increases linearly as B approaches zero. In this respect, we find a good agreement with the experiment, and may conclude that the ellipticity of the real QD device is close to our chosen value δ = 1.3. Secondly, the best overall agreement is found with a tilting angle θ = 50 • which is smaller than θ exp = 68 • ± 5 • reported in Ref. 15.
Another angle affecting the energy gap is the azimuthal direction of the magnetic field with respect to the axis of the ellipticity. In Figure 1(c) we plot the variation of ∆E(B) as a function of the azimuthal angle φ for three magnetic field strengths. We observe that the variation increases with increasing magnetic field, so that close to the crossing (B ∼ 3T ) the variation is ∼ 20%. Combining this observation with Fig. 1(b) leads to a conclusion that the best agreement with the experimental result is obtained with φ = 90 • , i.e., when the shortest axis of the ellipse is parallel to the tilting direction of the magnetic field. In the following we set φ = 90 • unless stated otherwise.
With a fixed δ obtained from the bending as B → 0 [see Fig. 1(b)], the singlet-triplet energy gap essentially depends on ω 0 and ω z . SO parameters are so small that they play only a minor role in the energy gap. The parameter values used here, ω 0 = 2.5 meV and ω z = 11.9 meV, give a reasonable agreement with the experiment in the energy separation and in the spin-relaxation rate (see below). An increase in both parameters and thus smaller vertical extension can also lead to good agreement in the energy gap, but a poor agreement in the spin relaxation rate. This effect is further illuminated in Sec. III.C. We point out that the results are not sensitive to using other types of frozen ground-state vertical basis functions and/or a larger number of vertical basis functions. The latter, as well as representing the vertical dimension by an infinite barrier at z = 0 and a harmonic confinement for z > 0 does not lead to significant changes. We conclude that good agreement with the experiment is obtained for a narrow range of harmonic confinement parameters. However, the theoretical magnetic tilting angle needs to be reduced to θ = 50 • to arrive to this level of agreement. The origin of this discrepancy pose a challenge for future experiments and theory. In the remaining part of this work we will stick to θ = 50 • .

B. Spin relaxation rate
In Fig. 2 we show the calculated total relaxation rate (solid line) in comparison with the experimental result (circles with error bars) in Ref. 15. Dashed and dashdotted lines correspond to the contributions from the piezoelectric coupling and the deformation potential coupling, respectively. We have used γ = 11.15 eVÅ 3 and α = 0 for the Dresselhaus and Rashba SO coupling strengths, respectively. Apart from the smallest and largest energy gaps, where the relaxation rate is very small, we obtain a very good agreement with the experimental measurements. We find that interchanging the values between the two couplings has only a minor effect on the transition rate; in fact, this is true for any combination with α 2 + (γ p 2 z ) 2 = 0.58 meVÅ. The experiment of Meunier et al. 15 gives a measure of this quantity but cannot resolve the relative magnitude of the two terms. Measurements of angular resolved phonons 30 or external fields controlling the Rashba term are needed to address the comparison between theory and experiment at a more detailed level.
In general, the coupling to the phonon bath is strongest when the QD size matches with the phonon wavelength. 2,7,15,27 The piezoelectric coupling (dashed line) is found to dominate the relaxation rate in Fig. 2 at small energy splittings (∆E 0.5 meV), whereas the deformation potential coupling dominates at larger ∆E. This is expected due to the √ ∆E and 1/∆E dependence of these couplings, respectively. In the former case, the coupling occurs through slower transverse phonons, which yields a peak at smaller ∆E.
Next we examine the effect of the SO coupling strength on the relaxation rate in more detail. Figure 3 shows the scaled relaxation rates for four different sets of α and γ in the units of meVÅ and eVÅ 3 , respectively. The black solid line corresponds to Fig. 2. The other curves are divided by the actual α 2 + (γ p 2 z ) 2 and multiplied by the square of the reference value, γ = 11.15 eVÅ 3 . Overall, no significant qualitative changes in the rate are obtained within a realistic parameter range. In particular, the peak position is insensitive to both γ and α, and thus we cannot achieve a better agreement with the experimental  15 (circles with error bars). The other parameters are the same as in Fig. 1. The solid line shows the total rate. The dashed and dashdotted lines show the contribution from the piezoelectric and deformation potential coupling, respectively. The spin-orbit coupling parameters are γ = 11.15 eVÅ 3 and α = 0 for the Dresselhaus and Rashba terms, respectively. rate in Fig. 2 (circles) by tuning either the Dresselhaus or Rashba SO coupling, or both. Furthermore, the fact that all curves are very close to each other indicates that the transition rate scales to a good agreement as The scaling formula can be derived from an interaction Hamiltonian V I ∼ αV 1 I + γV 2 I in the first-order perturbation limit when the interfering paths can be neglected, i.e., From Fig. 3 we see that the scaling for fixed θ, φ is valid. This suggests that the spin relaxation through the Rashba or Dresselhaus coupling in fact follows noninterfering paths within our basis states. Generally, however, the SO coupling is anisotropic and depends on the orientation of the magnetic field with respect to the crystal axis. 2 Consequently, the relaxation rate is anisotropic as well. 31 Hence, in view of a tilted magnetic field it is difficult to precisely assort the Rashba and Dresselhaus contributions in our QD device.
In Fig. 4 we show how the relaxation rate for fixed magnetic field strength (here B = 2 T) depends on the field direction with any combination of θ and φ. As pointed out above, the relaxation rate is strongly dependent on the energy splitting and qualitatively follows a form ∆E n exp (−c∆E 2 ) with positive constants n and c. At B = 2 T the rate is small for both θ ∼ 0 and θ ∼ 90 • , corresponding to small-∆E and large-∆E limits in the qualitative formula, respectively. At intermediate values of θ we find an area of higher relaxation rates that increases as a function of φ. In this region of increased rates the energy splitting gives phonon wavelengths comparable with QD size. Detailed analysis of the angular dependence is not the scope of this work, but we finally point out that these results could be used as a guideline in the analysis of forthcoming experiments where the direction of the magnetic field can be varied; in principle the angular properties could be used indirectly to obtain precise structural information of the QD device.

C. Effects of dimensionality
Finally we turn our attention to the effects of the dimensionality on the singlet-triplet energy splitting and on the spin relaxation rate. We compare our 3D scheme as described in the beginning of Sec. II to a bare 2D approach, where only the in-plane confinement and 2D basis functions are applied and the 2D approximation for the electron-electron interaction is used (corresponding to ω z ≫ ω 0 . This leads to relatively stronger electronic repulsion in the ground state than in the excited states, which in turn yields a smaller energy splitting. The results of the comparison are summarized in Fig. 5. First we notice that the tilted magnetic field has a very different effect on the energy states in 3D and 2D (a). In 2D the applied magnetic field closes the energy gap already at B ∼ 2 T, whereas in 3D the closing occurs at B ∼ 3 T. Therefore the 2D model cannot yield a reasonable agreement with the experimental splitting as shown in Fig. 5(b). It should be noted that if the confinement strength ω 0 in the 2D model is slightly increased to yield the correct splitting at B = 0, the magnetic-field dependence is still wrong, and the rate becomes much worse than the 3D result in comparison with the experiment. Furthermore, as shown in Fig. 5(c) the 2D model is unable to reproduce the experimental spin relaxation rate. In contrast, the 3D model has a very good agreement with the experiment as discussed above within Fig. 2.

IV. SUMMARY AND DISCUSSION
To summarize, we have proceeded towards quantitative modeling of spin relaxation in semiconductor quantum dots containing two electrons. Our approach covers (i) a three-dimensional description of the quantum dot device, (ii) numerically exact treatment of the electronelectron interaction, (iii) a tilted magnetic field, (iv) the ellipticity of the quantum dot, and (v) both the Rashba and the Dresselhaus spin-orbit coupling, and the latter with both linear and cubic terms. We have attributed the observed nonlinear behavior of the singlet-triplet energy gap at B 0.5 T to the ellipticity of the quantum dot. Then we have found a good agreement between theory and experiment in the triplet-singlet energy splitting and in the relaxation rate, although the optimal tilting angle of the magnetic field was found to be considerably smaller than the experimental value. Finally we have explicitly shown that the three-dimensional model is essential to obtain a reasonable agreement with the experiment in both the singlet-triplet energy splitting and the spin relaxation rate.
In view of our results supplied with aspects (i)-(v) listed above, the origins behind the remaining quantitative deviations between experiment and theory -especially regarding the different angle of the magnetic field -is unclear. However, as the most probable scenario we may suggest that the actual confining potential deviates from harmonic, e.g., by being more strongly confined at the edge of the quantum dot. This can be induced by the gate geometry of the lateral device and/or the nonuniformity of the confinement in the vertical (z) direction. These conditions can increase the high sensitivity of the spectrum to the tilted magnetic field (Fig. 1) even further. In this respect, we hope that the present study motivates more efforts in this direction, both theoretically and experimentally.