Hole spin relaxation in InAs/GaAs quantum dot molecules

We calculate the spin–orbit induced hole spin relaxation between Zeeman sublevels of vertically stacked InAs quantum dots. The widely used Luttinger–Kohn Hamiltonian, which considers coupling of heavy- and light-holes, reveals that hole spin lifetimes (T1) of molecular states significantly exceed those of single quantum dot states. However, this effect can be overcome when cubic Dresselhaus spin–orbit interaction is strong. Misalignment of the dots along the stacking direction is also found to be an important source of spin relaxation.


Introduction
The spin of carriers confined in nanostructures has been intensively studied in recent years due to its promising applications in spintronics and spin-based quantum information processing [1,2]. In particular, hole spins in quantum dots (QDs) have received great attention as a consequence of their long decoherence times and versatility. The confinement in heterostructures is responsible for the suppression of the main decoherence mechanisms present in bulk systems [3,4]. Additionally, the p-type symmetry of the valence band orbitals causes a weak hyperfine interaction with the lattice nuclei, thus giving rise to decoherence times potentially longer than those of electron spins [5][6][7][8][9][10][11]. This has enabled successful hole spin initialization [12] and coherent control [10,13]. Double quantum dots (DQDs) are a natural extension which should facilitate the use of independent optical transitions for spin preparation, manipulation and readout [14], as well as the scalability towards multiple qubit architectures [15]. Also, DQDs are more versatile than single QDs since the coupling between the two QDs offers an additional control mechanism, as the tunneling can be tuned by using externally applied electric fields [16][17][18][19].
Using the spin of holes in qubits requires control over the hole spin relaxation (T 1 ) and decoherence (T 2 ) times, the latter being related to the former at low temperatures [20]. In the presence of external magnetic fields, the main mechanism of spin relaxation for the valence band is usually phonon scattering mediated by spin-orbit interaction (SOI) [21,22]. Indeed, the strong SOI of holes is responsible for some of its characteristic properties, e.g. the g-factor anisotropy or the origin of antisymmetric ground states in DQDs [23][24][25][26]. The three main SOI mechanisms are the light hole-heavy hole (LH-HH) mixing, the bulk inversion asymmetry (BIA or Dresselhaus SOI) [27] and the structural inversion asymmetry (SIA or Rashba SOI) [28]. Several works have theoretically addressed the hole spin relaxation in single QDs taking into account different SOI mechanisms and have showed that one or another prevail depending on the QD traits [20,[29][30][31][32]. By comparison, the spin relaxation of holes in DQDs is still poorly understood. This is inspite of their promising prospects for the development of quantum information architectures [10,15,19]. Understanding the hole spin dynamics in DQDs is also of relevance for recently proposed exciton spin based qubits [33], since hole relaxation usually determines the excitonic spin lifetime.
In the present work, we investigate the hole spin relaxation between Zeeman split sublevels in vertically coupled DQDs. We consider InAs/GaAs DQDs with various relative positions of the individual dots while maintaining their dimensions unaltered. In particular, we study the spin relaxation for different interdot barrier thicknesses and dots alignments. The hole states are simulated using three-dimensional (3D) Hamiltonians including not only quadratic-in-k LH-HH coupling present in the Luttinger-Kohn Hamiltonian, but also cubic Dresselhaus SOI, which was found to be potentially important in single InAs QDs [31]. Calculations are carried out for varying strength of an electric field applied along the DQD stacking direction. This makes possible to study the transition from atomic-like states confined in one of the constituent QDs to fully molecular-like states, which are obtained when the electric field tunes the energy of the upper and lower dots to be the same [19].
We show that T 1 of molecular spin-orbitals is often larger than that of holes localized in single QDs, with lifetime extensions reaching 600% in some cases, which is a result of the higher symmetry of the system under resonant electric fields. Dresselhaus SOI however plays an important role in the description of the hole spin dynamics. Its inclusion in the Hamiltonian provides new channels of spin admixture, decreasing T 1 up to one order of magnitude and reducing the differences between molecular and single QD states. Nevertheless, the most drastic factor reducing T 1 is misalignment between the QDs forming the DQD. The severe symmetry breaking enhances SOI, leading to pronounced shortenings of T 1 , as well as to the appearance of spin anticrossings in the energy spectra, which constitute spin-hot spots with extremely fast relaxation.
The paper is organized as follows. In section 2 we present the details of the model we use to compute the hole states and the spin lifetimes. In section 3 we show and discuss the results of the calculations for a DQD system with strong and weak hole tunneling. Finally, conclusions are given in section 4.

Theoretical model
We investigate the spin relaxation of holes confined in vertically coupled DQDs grown along the [0 0 1] direction. The system is subject to external magnetic and electric fields applied along the growth direction. The Hamiltonian that describes the hole states reads where H L is the Luttinger Hamiltonian [34] and H B represents the terms coming from the implementation of the magnetic field. The next two terms in (1) are the 3D confining potential V QD and the externally applied electric field ( ) = F F 0, 0, z , with e standing for the hole charge and I the × 4 4 identity matrix. Finally, H BIA is the spin-orbit Hamiltonian corresponding to the Dresselhaus SOI. Rashba SOI is disregarded in this study because the system asymmetry in the growth direction is suppressed under resonant electric fields, which lead to the formation of molecular states with effective parity symmetry [23]. All the same, we have carried out a series of calculations (not shown) that confirm the negligible influence of Rashba SOI in the vicinity of the resonant field.
The Luttinger Hamiltonian H L is a four-band Hamiltonian which includes the spin-orbit coupling between light-holes (LH) and heavy-holes (HH) subbands up to second order in k [34]. The matrix representation of H L in terms of the Bloch basis functions J z = +3/2, +1/2, −1 Here m 0 is the free electron mass and γ i are the Luttinger mass parameters. For the sake of simplicity, we use constant Luttinger parameters throughout the entire nanostructure.
The uniform axial magnetic field is described by the vector potential in the symmetric gauge The implementation follows the procedure described in [35]. The resulting Hamiltonian H B has the form 2 . The last term of (7) corresponds to the Zeeman splitting with κ standing for the hole g factor, µ B the Bohr magneton and J z the matrix representation of the third component of the angular momentum with quantum number J = 3/2.
The last two elements in (1), H BIA and H SIA , are additional terms describing the Dresselhaus and Rashba SOI, respectively [36]. As stated above, we disregard H SIA due to its negligible influence under resonant electric fields. H BIA includes linear and third order in k terms and is given by: AB BA , 1 2 and cp stands for cyclic permutations of the preceding terms. The matrix form of Hamiltonian (8) can be found in [31].
Since H h is a four-band Hamiltonian, its eigenfunctions are four-component spinors of the form: and ⟩ |J z are the envelope and Bloch parts of the J z component. We study the spin relaxation of hole states mediated by the phonon bath. The transitions considered involve the Zeeman split sublevels of lowest energy, i.e. from the first excited to the hole ground state. The energy splitting of these states is small and, thus, only long-wave acoustic phonons can participate and the linear dispersion regime holds, Here, λ c is the phonon velocity of the longitudinal or two transversal phonon modes λ. The Hamiltonian that describes the coupling between holes and phonons is h ph pz dp (10) where e is the hole charge, I is the × 4 4 identity matrix, φ pz the piezoelectric potential and H dp the deformation potential term. These are the two scattering mechanisms that dominate the hole spin relaxation [30]. The deformation potential coupling H dp is given by the strain Hamiltonian [37], formally identical to (2)-(6) with k k i j replaced by the strain component ε ij and the mass coefficients by the deformation ones, and the piezoelectric interaction by the potential [38] ( ) (11) where ε r stands for the relative dielectric constant, h 14 is the piezoelectric constant and ε ij is the strain tensor component. The complete expressions and derivation of the piezoelectric potential and the deformation potential operators for the three phonon branches is presented in [31].
The transition rate between the initial hole state ⟩ |Ψ i h and the final hole state ⟩ |Ψ f h is estimated using the Fermi golden where λ − H q h ph is the hole-phonon coupling Hamiltonian, equation (10), The sum is done over all directions of q and the three phonon branches of bulk zinc-blende crystals, one longitudinal and two transversal. All calculations are carried out at zero temperature, so that only phonon emission processes are possible.
The multi-band Hamiltonian (1) is of fully 3D nature, as in DQDs the vertical and lateral dimensions can be comparable, which prevents the adiabatic separation of variables usually employed for single QDs [1]. Besides, we are interested in analyzing the effect of misalignment between the QDs forming the DQD, which implies breaking the in-plane symmetry through V QD . We then use Cartesian coordinates. It is also worth noting that H h includes third-order derivatives through the H BIA term. Since SOI is a fine effect, precise estimates of its influence require a very accurate description of such derivatives. To solve this challenging problem, we integrate H h numerically using a finite differences scheme. In order to achieve the desired accuracy while maintaining a reasonable mesh, we have explored increasing the points of the derivatives discretization. After a series of convergence tests, we found that a 5-point stencil central difference scheme offers a good description at a reasonable computational cost. In addition, we consider QDs with cuboidal shape which are perfectly adjusted to the 3D regular mesh. This simplified geometry grants comparable accuracy in the estimate of potential and kinetic energy terms, while allowing us to study the influence of the basic features of DQDs on the spin dynamics, namely the effect of interdot barrier thickness and that of misalignment. The 3D discretization of (1) yields a huge sparse matrix that is diagonalized by means of the Arnoldi iterative method. The size of these matrices ranges from × 73 644 73 644 for close aligned QDs up to × 150 326 150 326 for distant misaligned QDs.

Numerical results and discussion
The system studied is a DQD of InAs embedded in a GaAs matrix as represented in figure 1. The QDs are identical with cuboidal shape (L x = 20 nm, L y = 20 nm and L z = 3 nm) and are separated by a distance d. The parameters used in the calculations are summarized in table 1. They all correspond to the QD material InAs, except for the ones describing the phonons (c l , c t and ρ) which correspond to the matrix material GaAs as we assume bulk phonons. The confining potential V QD is zero inside the QDs and V 0 outside, where V 0 = −0.33 eV is the valence band offset between InAs and GaAs [39]. Finally, we take κ = 4/3 for the hole g factor as suggested in [40].
We investigate the dependence of the hole spin lifetime on the external electric field F z . We consider two different interdot barriers d = 3 nm and d = 9 nm as an example of a DQD system with strong and weak tunneling, respectively. In addition, for each d we study the possibility of the two QDs being perfectly aligned and also misaligned. The misalignment consists in a shift in x in opposite directions of the two QDs by an offset ∆ = 3.3 x nm as depicted in figure 1. This value corresponds to a DQD with large yet realistic lateral offset [42]. All calculations are carried out at a magnetic field B = 2 T. For this field, hole-phonon coupling via deformation potential prevails over piezoelectric coupling [31]. Figure 2 illustrates the energy spectra and hole spin lifetimes of a DQD with interdot distance d = 3 nm as a function of the external electric field F z . An analysis of the low-energy hole states reveals that they have mainly HH character 1 . Thus, the transition between the Zeeman-split sublevels, indicated in figures 2(a) and (b) by orange arrows, is essentially a transition from a HH with J z = +3/2 (⇑ in figure 2) to a HH with J z = −3/2 (⇓). Panels (a) and (c) of figure 2 show the energy spectrum and spin lifetime for a DQD with no misalignment. For a finite electric field the wave function tends to localize in one of the dots as represented in the insets of figure 2(a). The change of localization when varying F z gives rise to a charge transfer anticrossing at F z = 0, where hole states of both QDs are in resonance and the wave function forms delocalized molecular orbitals. Since the barrier thickness we consider is beyond the critical distance at which the hole ground state changes from bonding to antibonding character [23], the two states of the Zeeman doublet we consider are antibonding.

Strong tunneling regime
Calculations of the hole spin lifetime are shown in figure 2(c) for two different situations: considering only the LH-HH mixing derived from the standard Luttinger Hamiltonian H L (black line) and considering the additional influence of H BIA as well (dashed line). When only LH-HH mixing is taken into account, one can see that T 1 presents a maximum for molecular states (F z = 0) that slowly decreases as we move toward localized, atomic-like states ( F 10 z | | > kV cm −1 ). This remarkable result is reminiscent of the T 1 enhancement recently reported for single QDs at certain values of F z [32]. The inclusion of H BIA , however, reduces T 1 by one order of magnitude and moderates the longer lifetime of the molecular regime ( ≈ F 0 z ). Since spin lifetime is connected with SOI-induced spin admixture [1], in order to understand the above results we analyze the strength of the SOI, which can be qualitatively estimated from symmetry considerations. In general, a lowering in symmetry implies the activation of additional SOI mechanisms [36] and hence shorter T 1 . In particular, this is specially relevant at F z = 0, where the symmetry is highest. H L in (1) has T d symmetry, corresponding to the bulk zinc-blende point group. The confining potential that defines the aligned DQD system, V QD , reduces the symmetry to D 4h . Then, the application of an external magnetic field in the axial direction further reduces it to C 4h . We can take this as the starting point, H L at F z = 0 in figure 2(b). Next we add other factors like external electric fields or additional SOI terms, which further reduce the symmetry and hence T 1 . Thus, a finite electric field F z lifts the parity symmetry in z, reducing the system symmetry to C 4 , which explains the shorter T 1 for individual QDs as compared to the symmetric DQD in the H L curve of figure 2(c). If we include H BIA to the calculation instead, the reduction is more important, group C 2 , and thus the decrease of T 1 is more pronounced. In this case, the introduction of an external electric field no longer reduces the symmetry and has negligible effect on T 1 .
The results for a DQD with misalignment are illustrated in figures 2(b) and (d). Energetically, the main consequence of introducing a lateral offset to the QDs position is a reduction in the hole tunneling which, in turn, diminishes the magnitude of the charge anticrossing in the energy spectrum [42]. As for the hole spin lifetimes, by comparing with figure 2(c), figure 2(d) shows that misalignment of the QDs roughly preserves the shape of T 1 estimated from either H L or H BIA , but it causes an additional decrease in T 1 of one order of magnitude or more.
These results can be explained following the same reasoning as for a DQD with no misalignment. Now, the interplay of confinement potential V QD and magnetic field removes all exact symmetry elements. Therefore, the hole spin lifetime is reduced in comparison to the aligned case. The symmetry breaking becomes even more efficient in the presence of an applied electric field, resulting in shorter lifetimes.

Weak tunneling regime
In this section we investigate the same situations as above but we consider now a DQD system with an interdot barrier d = 9 nm. Figure 3(a) illustrates the energy spectrum when the QDs are aligned. A clear difference with the previous case, figure 2(a), is observed since now the J z = +3/2 antibonding and J z = −3/2 bonding states cross near the resonant field (see grey dotted boxes in figure 3). The inset of figure 3(a) contains a zoom of the crossing at | | ≈ F 0.5 z kV cm −1 . It shows that no spin anticrossing takes place in spite of including BIA spinorbit terms in the calculations.
The absence of anticrossings can be understood analyzing the symmetry of the hole states. H h belongs to the C 2 point group. The symmetry of its spinorial eigenstates Ψ n , equation (9), is then obtained from the double group C 2 . As shown in the labels of figure 3(a), the two states crossing each other have E 3/2 and E −3/2 symmetry, respectively (see appendix for more details). The different symmetry of the states justifies the absence of spin anticrossings in the spectrum.
The results for the spin lifetime are presented in figure 3(c). All calculations are also carried out considering the transition between the Zeeman split antibonding states (see orange arrows in figure 3). In general, we find similar lifetimes compared to the aligned DQD with strong tunneling. As in the strong tunneling case, H L predicts the longest T 1 values for the molecular regime (F z = 0), although the enhancement is now smaller and takes place for a narrower range of electric fields, which is a consequence of the weaker tunneling energy. Interestingly, in this case adding H BIA preserves the T 1 maximum at zero electric field. This is because the dominant term of H BIA (b 41 in (8)) scales roughly proportional to ⟨ ⟩ k z 2 [1]. In a DQD with thick interdot barrier a small electric field is enough to localize the wave function in one of the QDs, rapidly increasing ⟨ ⟩ k z 2 and reducing T 1 . In contrast, when the tunneling is important, stronger fields are needed to break the molecular character. As a consequence, in figure 2(c) ⟨ ⟩ k z 2 did not change significantly in the range considered and the T 1 maximum was less pronounced. μeV when Dresselhaus SOI is included. The presence of a spin anticrossing implies the admixture of hole states with the same symmetry. The reason is that misalignment reduces the symmetry to the C 1 point group, and all the states belong to the totally symmetric irreducible representation.
The obtained lifetime results are summarized in figure 3(d). Similarly to the strong tunneling case, misalignment of the QDs reduces the spin lifetimes by one order of magnitude. The main difference with respect to the strong tunneling case (figure 2(d)) is the appearance of two sharp dips in the hole lifetime at ≈ ± F 0.5 z kV cm −1 , where T 1 decreases by two orders of magnitude. These correspond to the spin anticrossings of figure 3(b). These anticrossings form so-called spinhot spots (see, e.g. [43]), where spin mixing is maximized. While the strong spin mixing can be benefitial for spin control protocols [42], our calculatios show that it also leads to severely reduced T 1 .
We have checked the robustness of the results in this section versus small changes in the model parameters (geometry, size and magnetic field). As discussed above, the qualitative trends are a consequence of the symmetry in the Hamiltonian rather than a specific set of parameters.

Conclusions
We have investigated the hole spin lifetime in InAs/GaAs DQDs as a function of the axial electric field strength. We have explored the effect of changing the QDs relative position (alignment and distance) and the introduction of Dresselhaus SOI.
The results reveal that severe changes in the distance d do not translate in large changes of the spin lifetime. A clear correlation between the symmetry of the system and T 1 is observed, which follows from the enhancement of SOI induced spin admixture with lowering symmetry. Thus, we show that the Luttinger Hamiltonian yields maximum T 1 for DQDs under resonant electric fields, but it decreases when the electric field pushes the wave function towards one of the QDs. This is a consequence of the higher symmetry of DQDs under resonant electric fields, when wave functions with parity symmetry are obtained.
Cubic Dresselhaus SOI lowers the symmetry to a C 2 point group, consequently reducing T 1 about one order of magnitude. In fact, a strong Dresselhaus SOI, as that found for DQDs in the strong tunneling regime, can even suppress the different behaviour between single QD and DQD states.
Misalignment of the QDs, which is often observed in DQDs, reduces the system symmetry to C 1 . When severe, it can reduce of T 1 over one order of magnitude. It is also responsible for the appearance of spin anticrossings in the energy spectra, which are absent for aligned DQDs. These are spin-hot spots, where spin mixing is maximized and a drastic decrease of T 1 is observed.
In summary, while the increase of symmetry, reached by the formation of molecular orbitals induced by resonant electric fields, favors spin purity and long T 1 values, misalignment or defects as well as strong Dresselhaus spin-orbit interaction play in opposite direction and can eventually overcome the effects of the resonant field. As for the envelope parts, f J n z , we consider that the symmetry of the matrix element operators in the 4-band Hamiltonian H h , obtained with the help of (1), is: must have A symmetry and it follows that the spinor symmetry is E −3/2 . For the other Zeeman sublevel, Ψ 2 , the dominant component is the spin up HH, so that ( ) f 3/2 2 has A symmetry and the spinor symmetry is E 3/2 . A similar procedure is followed to obtain the symmetry of the remaining hole states in figure 3.