Phonon-assisted relaxation between triplet and singlet states in a self-assembled double quantum dot

We study theoretically phonon-induced spin dynamics of two electrons confined in a self-assembled double quantum dot. We calculate the transition rates and time evolution of occupations for the spin-triplet and spin-singlet states. We characterize the relative importance of various relaxation channels, including two-phonon processes, as a function of the electric and magnetic fields. The simulations are based on a model combining the eight-band \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{k}\!\cdot \!\varvec{p}$$\end{document}k·p method and configuration-interaction approach. We show that the electron g-factor mismatch between the Zeeman doublets localized on different dots opens a relatively fast triplet-singlet phonon-assisted relaxation channel. We also demonstrate, that the relaxation near the triplet-singlet anticrossing is slowed down up to several orders of magnitude due to vanishing of some relaxation channels.


Model
We consider two self-assembled, vertically stacked InAs/GaAs QDs 16 . The composition gradient in the dots is modeled as a trumpet shape 33 with the maximum In content of 0.7 and the minimum of 0.4 (see Fig. 1). The mathematical models of the dot geometry and composition are described in Ref. 34 . The heights of the lower ("l") and the upper ("u") dot are taken h l = 8.5 a G and h u = 9.5 a G respectively, where a G = 0.565325 nm is the GaAs lattice constant. The dots are placed on wetting layers of 0.4 In content and thickness of a single a G . The parameters defining approximate QD base radii are chosen as r l = 28 a G and r u = 30 a G . The material intermixing is accounted for via a Gaussian blur (where we took 0.6 nm of the standard deviation). The simulations are performed for the axially oriented magnetic (B) and electric (F) fields.
We calculate single-electron states |φ n � using the eight-band k·p model 35,36 . We take into account strain distribution 37 and piezoelectric field 38,39 in the system. We used computational box 200a G × 200a G × 200a G for strain simulations and 100a G × 100a G × 70a G for finding the single-particle states. The model and the details of its implementation are given in Ref. 40 , here we use an extended version, described in Ref. 41 . The electric field is incorporated in the Hamiltonian via the standard band-diagonal term where z 0 is taken in the middle of the computational box. To preclude escaping the electron through the barrier in high electric fields, we take constant values: V (F) l = V (F) (z l ) for z < z l and V (F) u = V (F) (z u ) for z > z u , where z l , z u is the bottom of the lower dot (counted with the wetting layer) and the top of the upper one, respectively. Finally, we calculate the two-electron states where |vac.� is the crystal ground state, the coefficients c   www.nature.com/scientificreports/ single-electron states (i.e. two spin-dependent s-type states in each dot). In an idealized case (no spin-orbit interaction and no tunnel coupling), one can obtain the well known singlet/triplet configurations where (N l , N u ) describes the nominal occupations in the dots, a † (l/u)(↑/↓) is the creation operator for the electron state in the lower/upper dot and ↑ / ↓ denotes the spin orientation. Because of the SO coupling, spin is no longer a good quantum number. Furthermore, spatial configurations are mixed by the tunnel coupling. However, spinmixing effects are generally small, and the above notation is useful to classify well localized states (far from the tunnel resonances).
The spin-related properties of many-particle states can be characterized by the D j irreducible representations of the full rotational group, where j is related to the total angular momentum 42,43 . In the case of a two electron system, the direct product of representations gives D 1/2 ⊗ D 1/2 = D 0 + D 1 . In consequence, states can belong to the one-dimensional trivial representation D 0 (singlet states) or to the three-dimensional D 1 (threefold degenerated triplet states). In the presence of the SO coupling, the geometrical symmetry breaking affects also the spin degree of freedom. For the QDM considered here, at B = 0 T the system is described by the C 2v symmetry point group. In this group, D 1 splits into one-and two-dimensional representations 42,43 , which lifts the degeneracy of the T ± and T 0 triplet states.
Spin-orbit interaction creates various channels for phonon-assisted spin-flip processes. One class of such effects is driven by spin admixture mechanisms, where the single-particle state with some nominal spin orientation gets a contribution of the opposite spin 30,[44][45][46] . The other class contains direct spin-phonon coupling mechanisms [46][47][48] , which in QDs are typically weaker compared to the channels due to the spin admixture 30,44,45 .
The phonon-induced transition rates between the two-electron states can be calculated using the Fermi golden rule 49 where ω mn = (E n − E m )/ is related to the energy difference between the initial and final state, n B (ω) is the Bose-Einstein distribution, denotes the phonon branch, and ω k, = c k with branch-dependent speed of sound c . Finally, for the two-electron states where where H int (k, ) is the carrier-phonon interaction Hamiltonian for a single phonon mode. We take into account the deformation potential and piezoelectric electron-phonon couplings 50,51 . The detailed description of the model and carrier-phonon Hamiltonian is given in Ref. 41 .

Results
In this Section we discuss the calculated phonon-assisted relaxation rates. We also present the simulations of quantum kinetics for the two-electron singlet and triplet states.
Single-phonon spin relaxation rates. The energy spectrum for the lowest two-electron states is shown in Fig. 2a. At F = 0 , the ground state has (approximately) a configuration of S(0, 2), where both electrons are localized in the upper dot (which is assumed to be the larger one). The next three energy levels are nominally triplet states T ± ≡ T ± (1, 1) , and T 0 ≡ T 0 (1, 1) . At B = 0 , the energy splitting between them is up to several neV, hence can be neglected. The next state forms nominally the S(1, 1) configuration. Finally, the energetically highest state is S(2, 0).
For the nonzero electric field, we obtained the well known structure of energy branches with avoided crossings. Since the positive electric field decreases (increases) the energy in the lower (upper) dot, the states can be tuned into resonance. In consequence, at F ≈ 5.82 kV/cm and F ≈ 23.84 kV/cm, there are pronounced avoided crossings, corresponding to the tunneling of a single electron between the dots. The anticrossing at F ≈ 14.86 kV/ cm is very narrow, which is due to two-particle character of the involved tunneling process. Furthermore, at |S(0, 2)� = a † u↑ a † u↓ |vac.�, |T + (1, 1)� = a † u↑ a † l↑ |vac.�, www.nature.com/scientificreports/ F ≈ 14.87 kV/cm there is a "sweet spot", which corresponds to minimal difference (approximately 0.19 meV) between the lowest singlet and triplet configurations. First, we calculate the phonon-assisted relaxation rates between the two lowest singlet states (red lines in Fig. 2b,c). Depending on the magnitude of the electric field, the rates can describe different tunnel transitions: . The two pronounced maxima corresponding to the energy resonances from Fig. 2a are related to the increased overlap between the wavefunctions (which is due to spatial delocalization). The emission of longitudinal acoustic phonons from quantum dots preferentially takes place in the strongest confinement direction 52 , which in our case is the z-th axis. A phonon-assisted transition can be thought of as taking place at the one QD or the other. Constructive interference of quantum amplitudes for these two processes happens when the phonon wave is out of phase at these two points, hence, the phonon emission is expected to be enhanced for the frequencies ω = πc l (2n + 1)/D , where n has integer values and c l is the (longitudinal) sound velocity 32,53 . This is manifested by the oscillations of the transition rates 53,54 .
Since the transitions between the singlet states are spin conserving, they do not significantly depend on magnetic field. The rates on the order of tens of ns −1 for the dots separated by about 10 nm are consistent with former predictions based on the single-band k·p approximation 32 .
Next, we calculate phonon-assisted transition rates from the triplet states to the lowest singlet state (Fig. 2b,c). Depending on the electric field, such processes can either involve tunneling [ T(1, 1) → S(0, 2) and For the interdot tunnel transitions the rates oscillate, while no oscillations appear for the charge conserving processes. In order to check the accuracy of our CI approach, we calculated the transition rates in an extended basis of 24 electron states. The results for F = 14.87 kV/cm, F = 20 kV/cm, and F = 26 kV/cm are marked as points in Fig. 2b,c. At F = 14.87 kV/cm (the sweet spot) the rates are up to 65% larger compared to the minimal basis of 4 states. This can be partially attributed to changes in the energy differences (in the extended basis the T(1, 1) − S(1, 1) splittings are up to 13% larger). Also in the region of the inter-dot tunneling the rates are shifted, which is further enhanced by larger T(1, 1) − S(2, 0) energy splittings (at F = 26 kV/cm these energy differences are up to 36% larger in the extended basis) and by possible changes in the localization of the electrons. This is, however, a merely quantitative correction. Moreover, in our calculations, we cover rates spanning over many orders of magnitude. As extending the basis is demanding in the computational cost (in particular for the two-phonon processes), we perform further calculations in the original basis of the 4 states, keeping in mind that the actual values near the sweet spot may be larger.
The transitions T ± → S(0, 2) and T ± → S(2, 0) can be considered (approximately) as a single-electron tunneling accompanied by the spin-flip, while the other electron plays the role of a passive spectator. Such transitions are driven (primarily) by the spin admixture mechanism related to the Dresselhaus spin-orbit coupling 55 . Also the rates Ŵ(T ± → S(1, 1)) corresponding to the charge conserving processes, at low magnetic field crucially depend on the Dresselhaus coupling.
The relaxation processes T 0 → S(0, 2) and T 0 → S(2, 0) exhibit different behavior compared to those involving the T ± . We have verified that, in this case, the spin admixture mechanisms described above play minor role. Instead, the dominant transition channel is related to the difference of the single-particle g-factors for the states localized on different dots. This is a real-space spin-orbital effect connecting the position to the spin degree of freedom. It can be interpreted with the effective Zeeman Hamiltonian where g 1/2 and σ (1/2) z are g-factors and the (z-th) Pauli matrices defined for a given single-particle orbital (1 or 2). While such a Hamiltonian conserves the projection of the total angular momentum J z , it does not commute with the J 2 operator, hence mixes the states belonging to different representations of the rotation group. This mixing effect was shown to be an important factor determining transport properties of a gate-defined double quantum dot 29 . Such an effect takes place only for many-body configurations. In the case of two electrons, the T ± states are unaffected because �S|H Z |T ± � = 0 . On the other hand, �S(1, 1)|H Z |T 0 � = 1 2 µ B (g 1 − g 2 )B z which mixes the states. This gives rise to a spin-mixing that leads to a phonon-assisted relaxation 30,44,45 . For the QDM system considered here, g u ≈ −1.05 vs. g l ≈ −0.88 for the upper and the lower dot respectively, opening a significant T 0 → S relaxation channel. In consequence, the transition rate Ŵ(T 0 → S(1, 1)) shows a minimum at F ≈ 14.88 kV/cm, which is near the point where the single-particle g-factors are equal [ g ≈ (g u + g l )/2 ] due to the delocalization of the electron states. This is very close to, although not exactly coinciding with, the S-T "sweet spot" at F ≈ 14.87 kV/cm.
The magnetic-field dependence of the lowest two-electron states at F = 26 kV/cm is shown in Fig. 3a. The energies of all the states contain a diamagnetic contribution ∝ B 2 . In addition, the energies of the spin-polarized states T ± show a Zeeman shift. Since the single-electron g-factors are negative, the T + state has a lower energy than T − . The calculated phonon-assisted relaxation rates between the states are shown in Fig. 3b. For the considered magnetic-field range, one can approximate the dependence Ŵ(T ± → S(2, 0)) ≈ ∓aB + c . The dominant component c depends (mainly) on the Dresselhaus spin-admixture mechanism 55 . The relaxation from the T 0 state follows Ŵ(T 0 → S(2, 0)) ≈ a 0 B 2 , where the parameter a 0 mainly depends on the difference between the singleparticle g-factors of the Zeeman doublets a 0 ∝ (g u − g l ) 2 , consistent with the discussion in terms of the effective Zeeman Hamiltonian, presented above. Finally, the relaxation T 0 → T + can be viewed as a single-electron spin-flip in the individual dots. The rate Ŵ(T 0 → T + ) depends on magnetic field as ∝ B 5 , which is consistent with the results for one electron in a single QD 30,44,45 .
Next, we calculate the energy branches (Fig. 4a) and the relaxation rates (Fig. 4b) for the electric field corresponding to the sweet spot ( F = 14.87 kV/cm). The observed dependencies of the rates result from the interplay of various SO mechanisms (which can interfere constructively or destructively), and Zeeman energy splitting (important for the considered energy scale). We have checked that if the influence of the Zeeman splitting were (artificially) neglected, the rates could be very well fitted by the expression Ŵ(T ± → S(1, 1)) = aB 2 ∓ bB + c . Hence, the pronounced minimum of Ŵ(T + → S(1, 1)) visible at B ≈ 0.65 T results (mainly) from a combination of the negative linear and the positive quadratic contributions. The triplet-singlet transitions are several orders of magnitude slower compared to the rates shown in Fig. 3b. This is partly related to the fact, that the phonon spectral density decreases at low frequencies, and the energy differences between the S and T states are now relatively small. In addition, a strong delocalization of the electron states decreases some of the spinorbit coupling mechanisms, that give rise to T ± → S(1, 1) and T 0 → S(1, 1) relaxation. On the other hand, the magnitude of Ŵ(T 0 → T + ) is very similar (with the difference of about 10%) to the case of F = 26 kV/cm. This is due to the fact, that the T 0 → T + transition involves spin-flips in the individual dots, which are not sensitive to the electric field. Two-phonon processes. As temperature grows, the role of two-phonon processes in spin relaxation increases, in particular between levels separated by a small energy difference [56][57][58][59] . In this section we calculate the two-phonon correction to the spin relaxation rates as a function of temperature and compare it to our results www.nature.com/scientificreports/ for single-phonon processes, presented in the preceding section. In addition, we calculate the singlet-triplet dephasing rate due to elastic phonon scattering involving virtual transitions to singlet states (2,0) and (0,2) 60,61 . The time evolution of the reduced density matrix ρ that describes the electron part of the system is determined (for a pure initial state) by the time-convolutionless (TCL) equation, which in the interaction picture can be written as 62 where K(t) is a generator. Within the TCL projection operator technique 62 , the generator can be expanded up to a desired order in the coupling (here represented by the H int part of the Hamiltonian). To simulate two-phonon processes, one needs to take into account terms up to the fourth order where the explicit definitions of K n (t) are given e.g. in Refs. 60,62 . Note, that K 4 (t) describes not only two-phonon effects, but contains also perturbative corrections to the single-phonon processes. The relaxation rates from a state i to j are calculated within the Markov limit [ K n (t → ∞) ≡ K (∞) n ] by investigating the terms in the equation of motion for �j|ρ|j� that are proportional to �i|ρ|i� . The rate includes the second-order contribution Ŵ (2) (i → j) , corresponding to the single-phonon transition described in the previous Section. Performing calculations similar to those described in Refs. 60,61 , we further obtain the fourth-order contribution, accounting for the two-phonon processes d dt ρ(t) = K(t)ρ(t), www.nature.com/scientificreports/ where P denotes the Cauchy principal value, ω ij = E j − E i / , and the spectral densities are given by In Eq. (1) the indices i, j refer to the singlet-triplet doublet of single-occupation states around the "sweet spot", while α, β run through all the states except to i and j. In the derivation of this result, we have selected only terms corresponding to energy-conserving two-phonon transitions between the two states, neglecting contributions that do not have a resonant two-phonon structure.
To assess the importance of the contributions coming from the two-phonon processes, we calculated the transition rates for two regimes of parameters. Figure 5a,b presents the temperature dependence of Ŵ (2) (T 0 → S(1, 1)) and Ŵ (4) (T 0 → S(1, 1)) for the electric field corresponding to the "sweet spot" F = 14.87 kV/ cm. While for low temperatures the single-phonon processes are much faster (by several orders of magnitude), the Ŵ (4) (T 0 → S(1, 1)) starts to dominate from T ≈ 25 K (at B = 0.1 T) or T ≈ 14 K (at B = 1.0 T). On the other hand, for the transitions T ± → S(1, 1) (not shown here) the single-phonon processes dominate in the whole temperature range. The relative importance of the two-phonon processes in the case of low energy splitting between the involved states is due to the super-ohmic character of the phonon spectral density that makes single-phonon transitions unfavorable at low phonon frequencies, while inelastic phonon scattering gains importance as the range of occupied phonon states enhances with temperature. We also performed calculations for the regime of tunneling, at F = 26 kV/cm and B = 0.1 T (not shown here). In that case, the two-phonon contribution is not significant (at T = 60 K it is two orders of magnitude smaller than the single-phonon rate).
(1) www.nature.com/scientificreports/ While singlet-triplet transitions rely on spin-orbit coupling and are therefore relatively slow, dephasing of singlet-triplet superpositions can also result from a two-phonon process that can be understood as elastic scattering of a phonon, during which the system undergoes a virtual transition to one of the doubly occupied singlet states and back 60,61 . Within the singlet-triplet doublet around the "sweet spot", this process is strong only for the singlet, when it is spin-conserving, while for the triplet it is very inefficient. This leads to distinguishability between the two states and, in consequence, to pure dephasing of singlet-triplet superpositions. The decay of coherence between the S(1, 1) and T 0 states appears in the second order as a result of transitions between these two states, as well as thermally activated spin-conserving transitions from S(1, 1) to S(2, 0) and S(0, 2), which are exponentially suppressed at low temperatures. The corresponding dephasing rate γ (2) , extracted from the K 2 term in the TCL expansion in the Markov limit, has the form (using short-hand indexing '1' and ' 3' for the states S(1, 1) and T 0 , respectively) which is consistent with the usual Lindblad equation obtained via Born-Markov approximation in the weak coupling limit 62 . To account for the fourth-order contribution, we use the approach from Ref. 61 and extract the dephasing rate γ (4) from the TCL generator K 4 , where indices α , β run through all the states except to S(1, 1) and T 0 . As shown in Fig. 5c,d, the dephasing is indeed much faster (a few orders of magnitude at the typical experimental temperatures about 4 K) than the relaxation between the two states. The two-phonon pure dephasing mechanism dominates also over the singlephonon contributions to dephasing, related to real transitions. This dephasing rate grows even more at electric fields closer to the resonances between the state S(1, 1) and S(2, 0) or S(0, 2) 63 .

Relaxation kinetics.
To investigate the kinetics of the system, we numerically solve the Master equation (within the Markov and secular approximations) assuming the initial occupation of the T 0 state. We focus on the regime of the tunnel transitions, where we took F = 26 kV/cm, and we consider magnetic fields of B = 0.1 T and B = 1 T. At T = 0 K (Fig. 6a,b) the evolution is exponential, where the only significant process is the direct relaxation T 0 → S(2, 0) . However, for a non-zero temperature, the picture becomes more complicated. The transitions to upper states become possible, and the final occupations form a temperature-dependent equilibrium (according to the Gibbs distribution). Furthermore, the Bose-Einstein distribution of phonons enhances the rates between states separated by a small energy difference. As shown in Fig. 6c,d, for T = 60 K the evolution   www.nature.com/scientificreports/ of occupations is no longer exponential and involves more states. The relaxation to the lowest singlet state can occur directly T 0 → S(2, 0) , but also through T 0 → S(1, 1) → S(2, 0) . Although the phonon-assisted transitions T 0 ↔ T ± are negligible, the states T ± get occupied through T 0 → S(1, 1) → T ± , and T 0 → S(2, 0) → T ± transitions. For a weak magnetic field (Fig. 6c), the occupation dynamics of all the triplet states exhibit comparable timescales. On the other hand, the higher magnetic field (Fig. 6d) leads to two distinct regimes: the fast transitions T 0 → S(2, 0) and T 0 → S(1, 1) in the first stage of the evolution, then slow transitions to the T ± . Such behavior results from different magnetic-field dependencies of the involved transition rates (see Fig. 3b). In order to quantitatively assess the importance of distinct transition channels for the kinetics presented in Fig. 6, we calculate the transition rates as a function of temperature (Fig. 7). For low temperatures, the relaxation is clearly dominated by the direct transition T 0 → S(2, 0) . However, with increasing temperature and at low magnetic fields, the transitions involving the T ± start to play an important role in the system dynamics. On the other hand, at higher magnetic fields ( B = 1 T in Fig. 7b), the transitions from/to the T ± states are very slow compared to the T 0 → S(2, 0) in the whole range of temperatures.

Conclusions
We have studied quantum kinetics of two electrons in a quantum dot molecule. With a realistic model of the QD system geometry, 8-band k·p method, and configuration-interaction approach, we have calculated two-electron states. We have investigated the phonon-assisted transitions between the triplet and singlet states in the presence of external magnetic and electric fields. We have considered the triplet-singlet transitions accompanied by tunneling as well as the case of occupation conserving relaxation. We have identified channels of the triplet-singlet relaxation that become important in different parameter regimes. While for weak magnetic fields the tunnel transitions related to spin-admixture mechanisms are dominating, the regime of moderate and strong magnetic fields favors another mechanism related to the difference between the electron g-factors in the dots. We have also shown, that near the "sweet spot", lifting of this mechanism leads to a considerably longer lifetime of the triplet T 0 state. We have also demonstrated a non-exponential quantum kinetics resulting from the interplay of various direct and sequential processes contributing to the relaxation. Finally, we studied the influence of two-phonon processes on the transition rates. We have shown that they can significantly contribute to the relaxation process at the "sweet spot", at temperatures of several Kelvins and above. We have also confirmed that singlet-triplet coherence near the "sweet spot" is limited by the pure dephasing process due to elastic phonon scattering, which is much faster than relaxation.