Spin and rotational symmetries in unrestricted Hartree–Fock states of quantum dots

Ground state energies are obtained using the unrestricted Hartree–Fock (HF) method for up to four interacting electrons parabolically confined in a quantum dot subject to a magnetic field. Restoring spin and rotational symmetries we recover Hund's first rule. With increasing magnetic field, crossovers between ground states with different quantum numbers are found for fixed electron number that are not reproduced by the unrestricted HF approximation. These are consistent with the ones obtained with more refined techniques. We confirm the presence of a spin blockade due to a spin mismatch in the ground states of three and four electrons.


Introduction
In recent years, semiconductor quantum dots [1] have been the subject of many experimental and theoretical investigations. Their electronic properties can be controlled with a high accuracy by applying external gate voltages and magnetic fields. Circular dots are especially interesting. At low magnetic field B, the addition energy of a two-dimensional circular dot has been found to exhibit pronounced peaks for N = 2, 6, 12 electrons [2]. This suggests that levels group themselves into shells, in analogy with atoms and nuclei and that closed-shell configurations are particularly stable. The behaviour of open-shell states has been found to be consistent with Hund's first rule [3]. At nonzero magnetic field, transitions between ground states involving total spins and angular momenta have been observed [4]- [6]. These have a profound effect on the transport properties of a quantum dot. For instance, if the total spins of two ground states with N and N + 1 electrons differ by more thanh/2, spin blockade of electron transport is predicted [7]. This has been recently observed in the conductance of quantum dots [8].
Theoretically, the electronic structure of quantum dots has been studied using many different techniques (see [9] and references therein). Exact diagonalization (ED) [10]- [18], configuration interaction (CI) [19,20], stochastic variational method [21] and the pocket state method [22] allow ground and excited state energies and their quantum numbers to be calculated with very good accuracy. Also quantum Monte Carlo (QMC) methods [23]- [27] have been employed: they provide accurate estimates for ground and excited state energies, although total spin symmetry is not always preserved [25,26]. By means of all these techniques, shell structure and Hund's rule have been analysed in detail. In addition, 'magic' values for the total angular momentum have been predicted to occur for dot states with a given total spin. They occur if electrons in a quantum dot strongly interact and arrange themselves in a rotating Wigner molecule [13,26] [28]- [31]. All the above methods are computationally very expensive and can be used for relatively low electron numbers N 13: only with Monte Carlo methods have electron numbers up to N = 24 been reached [24] at zero magnetic field.
Other methods are used for treating systems with larger electron numbers, like the Hartree-Fock (HF) approach [32]- [38] and density functional theory [39]- [42]. They generally provide less accurate estimates of the ground state energy and the resulting wavefunctions can have unphysically broken symmetries. The HF methods are paradigmatic, in that the variational ground 3 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT state wavefunction is a single Slater determinant which does not properly include correlations and in general is not an eigenfunction of the total spin [43]. In contrast with space restricted HF (RHF) methods, space unrestricted HF (UHF) methods [32,34] systematically allow for symmetry breaking allowing, as a starting point for the calculations, wavefunctions without rotational invariance even in the case of circularly symmetric dots. While UHF yields a lower estimate for the ground state energy, it can have severe drawbacks when considering physical properties of the ground state like the total spin. For instance, UHF calculations sometimes fail to predict Hund's rule, in contradiction to experiments, and in contrast to results obtained with other methods. In addition, wavefunctions with broken symmetries do not allow the total spin and the angular momentum to be determined. Methods to restore the correct symmetries were pioneered in the 1950s [44]- [46] and were used in the context of quantum dots, mostly for the rotational symmetry [47]- [51]. Restoration of the spin symmetry in quantum dots has received much less attention [47].
In this paper, we apply a systematic projection procedure to restore both the spin and the rotational symmetries of ground state variational wavefunctions obtained by UHF calculations for quantum dots with up to four electrons, including a magnetic field. It is our aim to show that, after all the symmetries are restored, the wavefunctions are considerably improved and show several physical features which are not reproduced by the straightforwardly applied UHF method. Restoring symmetries introduces correlations, absent within the single UHF Slater determinant, leading to better energy estimates.
We demonstrate the efficiency and accuracy of the projection procedure by comparing our results with those of the methods mentioned above. The main findings are: (i) at zero magnetic field, Hund's first rule is recovered for four electrons which has been claimed earlier to be violated by UHF [32]; (ii) for nonzero magnetic field, many crossovers between ground states with different total spins and angular momenta of up to four electrons are found that are completely missed by using the HF method alone.
The paper is organized as follows. In section 2, the UHF procedure is briefly sketched, emphasizing the role of broken symmetries. In section 3, the projection techniques employed in the symmetry restoration are discussed. Results for zero and nonzero magnetic field are discussed in section 4. The paper is concluded by pointing out the perspectives for obtaining results for higher numbers of electrons that are not accessible by other methods.

The unrestricted HF approximation
The Hamiltonian for a system of N interacting electrons, parabolically confined in the x-y-plane and subject to a perpendicular magnetic field B = Be z (e z is the z-axis unit vector, here and in the followingh = c = 1) is where H 0 (r, p, s z ) = p + eA(r)

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
with Coulomb interaction potential V(r) = e 2 /4πε 0 ε|r|, B = rotA, effective electron mass m * , confinement frequency ω 0 , effective g-factor g * and the Bohr magneton µ B . The z component of the ith spin is s zi = ±1/2. Furthermore, −e is the electron charge and ε 0 (ε) the vacuum (relative) dielectric constant. The Hamiltonian (1) commutes with the total angular momentum L = l 1 + · · · + l N (z component in two dimensions), the z-component of the total spin S z = s z1 + · · · + s zN and the total spin S = s 1 + · · · + s N . The single-particle term H 0 is exactly diagonalizable. It yields the Fock-Darwin (FD) spectrum n,m,s z = (2n + |m| + 1) + with the corresponding harmonic eigenfunctions φ n,m,s z (r) [52]. Here, n and m are principal and angular momentum quantum numbers. Further, we introduced the cyclotron frequency ω c = eB/m * , and the effective confinement frequency = (ω 2 0 + ω 2 c /4) 1/2 . When including interactions, the problem is in general not solvable. Here, we deal with the electron interactions by using as a starting point the HF approximation. The interacting, manybody ground state wavefunction is written as a single Slater determinant consisting of N orbitals of the form ψ s z (r), assumed to be eigenfunctions of s z ψ s z with a, b denoting the spatial parts, and α, β the spinors corresponding to s z = ±1/2, respectively. Denoting N + and N − the number of spin up and spin down orbitals, respectively, the Slater determinant is the eigenfunction of S z . Slater determinants are not in general eigenfunctions of S 2 , unless |S z | = N/2 when the HF solution has total spin quantum number S = N/2. Thus, HF solutions with S z < N/2 are not consistent with the spin symmetry of the Hamiltonian. This is known as spin contamination [43]: a single Slater determinant is in general a superposition of many eigenfunctions with different total spins but with a given S z . The determinant is variationally optimized in order to minimize the energy giving rise to the well known N coupled integro-differential equations for the orbitals ψ s z i (r), where N s z = N + δ s z ,1/2 + N − δ s z ,−1/2 and the electron density is

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
Within a given sector with fixed S z , equation (7) is solved self-consistently, starting from an initial guess for the orbitals. Spatially unrestricted initial guesses are used: orbitals are assumed such that they are not eigenfunctions of the angular momentum [32,34,47] (UHF method). With this, solutions of the HF equations corresponding to non-rotationally invariant electron densities, are obtained. These symmetry-broken solutions enhance the number of variational degrees of freedom, thus generally leading to better energy estimates. Therefore, in addition to the lack of total spin symmetry already present in RHF calculations, UHF solutions do not possess the spatial symmetry of (1). Their predictive power for total spin and angular momentum of the ground state is therefore limited. As a consequence, the symmetries of UHF wavefunctions must be suitably adjusted in order to address these properties.
It is important to note that for given N and S z , many local minima of the energy surface can be found depending on the initial guess used for the UHF calculation. In general, one finds a sequence of states | . For a given S z , one has to perform extensive scans over the space of initial conditions in order to achieve confidence that the first state in the above sequence is the best UHF state, with the lowest attainable energy. The UHF ground state is then defined as the state with the lowest among the energies E S z 1 . This also determines the value of S z for the UHF ground state.

Restoring symmetries
In order to reflect the symmetries of the Hamiltonian in the UHF solutions, projection operators [44,45] can be applied to the wavefunctions. To avoid confusion, in this section operators are denoted by a hat. ConsiderP L andP where we used the commutation rules presented above and the idempotency of projection operators. The spin projectorP S z S can be written as [44] Its action on the UHF wavefunction is [44,53] (11) where N < = min{N + , N − } and C q (S, S z , N, N + ) 6 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT are the Sanibel coefficients [53,54]. The term |T q = |T (1) q + · · · + |T (n q ) q is the sum of all of the distinct Slater determinants obtained by interchanging in the initial determinant | S z all q spinor pairs with opposite spins. The projection operator on the total angular momentum L is [45] Acting on |T q with exp (iLγ) results in |T q (γ) , a determinant where all the orbitals are rotated by γ around the z-axis. Combining (11) and (14) we have The projected state (15) is a sum of many Slater determinants. This indicates that a high degree of correlation has been introduced by applying the projection technique to the UHF scheme. As is clear from (15), spin projection implies calculations which involve many Slater determinants. Therefore the second form of (9) is especially useful from the computational point of view. A detailed description of the procedure is developed in appendix A for N = 3. We have implemented numerically the evaluation of (9), using (15) and standard theorems for the evaluation of the Hamiltonian matrix elements [46]. In this work, we will determine the projected ground state, defined as the projected state with the lowest energy. This procedure, labelled 'projected HF' (PHF) in the following, is far from being trivial: in particular, it is not sufficient to project the UHF ground state alone to determine the PHF ground state. If several UHF minima E S z i are almost degenerate, all the corresponding | For a given UHF solution (corresponding to fixed B and S z ), the spin projection involves the evaluation of all the overlapping matrix elements T 0 |Ĥ|T (i) q (γ) (see appendix A for details). This constitutes a computational overhead, especially for high particle numbers. For given S z and S, the number of matrix elements to be evaluated is given by n S z = N < q=0 n q . The worst case scenario occurs for the minimal values of S z . As an example, for even N, the number of matrix elements is (a similar formula with the same asymptotic behaviour holds for odd N). Numerically performing the angular momentum projection requires the discretization of γ ∈ [0, 2π] in n L points, followed by a fast Fourier transform (FFT) procedure. The lower bound for the latter is determined by the maximum angular momentum |L| L max desired, independent of N. We point out that, once the sampling has been performed, the FFT simultaneously produces all the angular momenta 7 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT projections for |L| L max . We have checked that, in order to achieve good convergence with L max = 20 (much higher than the highest |L| involved in the results discussed below), one needs n L = 256. Thus, for the case N = 4, S z = 0, fixed S and |L| 20, we need to evaluate n S z n L = 1536 matrix elements. This compares favourably with respect to more refined methods: for instance, ED for N = 4, S z = 0, S = 2, L = 14, needs up to 19 774 Slater determinants [15]. Also for increasing numbers of electrons PHF performs well: in the case N = 6, S z = 0, S = 0, L = 0 (not discussed in the present work), CI needs 661 300 configurational state functions (linear combination of Slater determinants) [20] while our procedure requires 5120 matrix elements to evaluate not only L = 0 but all |L| 20. We expect similar performance gains for even larger numbers of electrons.
In the next section, we will show that the PHF ground state has lower energy than the UHF one. Most importantly, all quantum numbers of the projected ground state can be determined, in contrast to UHF.

Results
Before presenting our results, we provide some technical details of the numerical method. We rewrite equations (7) by employing the non-interacting FD-basis and obtain a nonlinear Pople-Nesbet eigenvalue problem [43]. For each value of B, the 75 lowest-lying FD states per spin direction have been assumed to represent the basis. This guarantees a fair convergence of the UHF procedure for N 4. By comparing results with those obtained with a smaller basis set (55 elements), the relative uncertainty of ground state energies has been assured to be ∼10 −6 in all of the data shown below. For a given number of electrons, many different UHF solutions for all possible values of S z have been determined. Subsequently, the projection procedure has been implemented, taking particular care of the cases with almost degenerate UHF states as discussed above.

Zero magnetic field
At zero magnetic field, expressing energies in units ω 0 and lengths in units 0 = (m * ω 0 ) −1/2 , the Hamiltonian (1) depends only upon the dimensionless parameter which represents the relative strength of the interaction. Table 1 shows ground state energies for N = 2, 3, 4 electrons with λ = 1.89, chosen in order to compare our results with those obtained with diffusion Monte Carlo (DMC) [23] (third column). The second column shows the UHF energy, the third column the PHF ground state energy while the fourth column shows the DMC data. Since at zero magnetic field the spin multiplets are degenerate with respect to S z , the latter is not mentioned in columns 3 and 4. As already discussed, projection improves the ground state energy. The quantum numbers for the PHF ground state agree with the ones obtained with DMC. The computed energies differ for at most about 4% (N = 2 case) and improve with increasing N. Our results also agree with other published data: for N = 2, a singlet ground state has also been found by means of ED [10,11,55].   [23]. In columns 3 and 4 S z is not indicated since at zero magnetic field spin multiplets degenerate. Interaction parameter λ = 1.89 (see equation (17)). Energies in units of ω 0 . For N = 3, 4 electrons, the quantum numbers predicted by PHF agree with those predicted by the CI method [19,20] and the path integral Monte Carlo (PIMC) method [25]. The ground state for N = 4 is a triplet, in agreement with Hund's first rule. Previously published calculations [32] found that the UHF ground state for four electrons has S z = 0. Our UHF calculation, although performed with a confinement energy different from that used in [32], confirms this result (table 1, column 2). A ground state with S z = 0 previously has been interpreted [32] as violating Hund's first rule. However, the latter deals with the total spin of the electrons, an undefined quantity in UHF solutions. Projecting the total spin allows us to find that the ground state has S = 1 which, due to the degeneracy of the spin multiplets, is compatible with S z = 0.
To probe the spin properties of the N = 4 case, we study it for increasing λ. For strong interactions, the UHF method is commonly assumed to favour spin polarized states [55]. Thus, a crossover of the UHF ground state to S z = 2 is expected. Table 2 shows a comparison of UHF, PHF and CI [20] results. The UHF ground state for λ 4 has S z = 2. This is not compatible with a triplet state, and indicates a violation of Hund's first rule. The violation is striking since S z = 2 is an uncontaminated S = 2 state. On the other hand, as shown in column 3, for λ 2 PHF predicts a triplet ground state, consistent with Hund's rule and in qualitative agreement with the CI results (column 4). For λ 6, the PHF ground state originates from a low lying excited state (local minimum) of UHF (see below). When increasing the strength of the interaction, it is expected that correlations, measured by the difference between the energies of the ground state obtained by the RHF approximation and the exact ground state, become increasingly important [51]. We observe that the relative difference between ground state energies calculated with PHF and CI decreases from ≈2% to ≈1% as λ increases from 2 to 8. This can be attributed to an increasing amount of correlations introduced by the PHF procedure when interactions grow. This trend was already pointed out for N = 2 when spin and rotational symmetries were restored [51].
We close this section outlining the determination of the PHF ground state in the case N = 4 with λ = 6. Table 3 shows the lower energies (column 2) of different UHF states (column 1) found for S z = 0, 1, 2. The UHF ground state is |   (17)). Column 2: UHF ground state. Column 3: PHF ground state. Column 4: CI ground state, extracted from [20]. In columns 3 and 4 S z is not indicated since for B = 0 spin multiplets are degenerate. Energies in units of ω 0 .   also checked that UHF states with higher energies (not shown in the table) do not yield a better PHF ground state energy.

Magnetic field
For nonzero magnetic field, the projection procedure provides physical features that are not correctly reproduced by the UHF approximation. In the following, we assume standard parameters for GaAs: m * = 0.067m e , ε = 12.4 and g * = −0.44. We start with the N = 2 case with ω 0 = 3.37 meV as a confinement energy, in order to compare our results with the available ED data [55]. Figure 1(a) shows the ground state energy of the PHF ground state. Correspondingly, figure 1(b) shows L (solid) and S (dashed) quantum numbers. With increasing B, an alternating sequence of singlet and triplet states is observed as L increases with unitary steps. PHF results compare well with those obtained by means of more accurate methods [10,11,55]. They predict both the same sequence of increasing L and of singlet-triplet transitions. Quantitatively, the crossover fields are not perfectly reproduced. For instance the first singlettriplet transition occurs at B ≈ 2 T [55], while we find B = 1.38 T. The sequence of transitions between states with different total spins shown by the PHF solution cannot be obtained using the UHF approximation. Figure 1 found in the whole magnetic field range, in contrast to the PHF data. Figure 1(d) shows a comparison of the UHF and PHF ground states in the 0 T B 2 T region. Clearly, at low magnetic fields the PHF singlet solution has a better energy. In addition, it is hard to relate the S z transition occurring in the UHF calculation with the genuine singlet-triplet transition exhibited by PHF. For discussing a quantum dot with four electrons, we choose a confinement energy ω 0 = 6 meV in order to be able to compare our results with the ED calculations reported earlier [16]. The PHF ground state energy (figure 2(a)), shows kinks associated with spin and angular momentum transitions ( figure 2(b)). We find a sequence of transitions which agrees with those reported earlier [16]: a low field triplet-singlet transition obtained with PHF occurs at B ≈ 0.7 T, higher than the value B ≈ 0.5 T predicted by ED [16]. Other PHF transitions are shifted to smaller values, similarly to the singlet-triplet transitions for N = 2. As a comparison, we show (figure 2(c)) the UHF ground state energy (solid line) and S z (dashed line) for N = 4 in the same magnetic field range. Apart from the higher energy estimate provided by UHF, we observe that the spin transitions scenario is very different w.r.t. the one predicted by PHF. The UHF ground state has S z = 0 up to B ≈ 3.8 T, in contrast with the above mentioned tripletsinglet transition obtained with PHF for B ≈ 0.7 T. Subsequently, a narrow transition to S z = 1 is found. Then, for B 4.2 T the UHF ground state is fully polarized. The sharp contrast with the PHF results confirms that the UHF method is not reliable for qualitatively predicting the spin properties of the system, as we already observed in the N = 2 case.
One of the important predictions of the results in [16] is a spin blockade [7] that in the ED calculation occurs at 4.96 T < B < 5.18 T, due to a violation of the total spin selection rules.

Conclusions
We have applied the UHF method to quantum dots with up to four electrons in the presence of a magnetic field, and for varying strength of the interaction. We have used a systematic projection approach for simultaneously restoring the total spin and the rotational symmetries of the UHF wavefunctions. The projected wavefunctions are superpositions of many different Slater determinants. Thus, they contain important correlations, missed by UHF solutions. We found that the energies of the ground states obtained by PHF are systematically lower than those of the UHF ground states, although still higher than the ones found with other methods such as ED, CI or QMC. The PHF provides ground states with total spins and angular momenta in qualitative agreement with the 'exact' results. The data in table 2 show that by means of PHF, important correlations are introduced for increasing λ. This is signalled also by the tendency of the PHF method to reproduce the CI results [20] with increasing accuracy. We recover Hund's first rule for four electrons at zero magnetic field. With increasing magnetic field, crossovers between ground states with different quantum numbers are found for fixed electron number that are not reproduced by the UHF approximation. These are consistent with the ones obtained with the 'exact' techniques. We have confirmed the presence of a spin blockade due to a spin mismatch in the ground states of three and four electrons.
We conclude that the PHF approach is an important technique when dealing with ground state properties like total spin and angular momentum of correlated electrons in quantum dots.
In view of what has been discussed in section 3, we expect that PHF will be very useful to obtain interesting results for quantum dots with larger electron numbers that, so far, have not been accessible by other methods, especially in the presence of a magnetic field.

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
Such quantum dots have recently been experimentally investigated [56]. In particular, we expect that the wave function obtained by PHF will be useful to obtain transition matrix elements that are important for understanding the transport behaviour.
In order to illustrate the procedure leading to the evaluation of (9), we discuss here the case of three electrons. The possible UHF solutions have S z = ±1/2 or S z = ±3/2. For the latter case no spin projection is required since it is a pure S = 3/2 state, so that only angular momentum projection has to be performed. Therefore, we concentrate here on the S z = 1/2 case. The UHF wavefunction reads from equation (5) (the case S z = −1/2 is obtained by interchanging all α and β spinors, without any conceptual difference with respect to the following discussion). This state does not have defined total spin. It is a superposition of the doublet S = 1/2 and the quadruplet S = 3/2 states. The spin projection allows the total spin part we are interested in to be selected. To do so, we generate all the possible shuffled Slater determinants |T q starting from |T 0 and weight them by the corresponding Sanibel coefficients (12), from equation (11) where For simplicity, in the following we consider the S = 1/2 case only. In order to obtain (15) we applyP L on (A.2). The action of the rotation generator exp(iLγ) on the determinants composing (A.2) affects the spatial part of the orbitals only: a i → a i (γ) and b i → b i (γ). Where a i (γ), b i (γ) are the UHF orbitals rotated by an angle γ over the z-axis. Employing (9)  The evaluation of (A.6) and (A.7) is done using the standard theorems for many body wavefunctions [46]. The overlap terms in the denominator are T 0 |T 0 (γ) = det{d (0) (γ)} and T 0 |T 1 (γ) = T 0 |T (1) 1 (γ) + T 0 |T (2) 1 (γ) = det{d (1) (γ)} + det{d (2) (γ)} respectively, with The evaluation of the numerator, requiring the calculation of one-and two-body operator matrix elements, is more involved. The single particle part is where d (i) (k|l) (γ) is the (k, l) entry of the first-order cofactor of d (i) (γ) and h (i) kl (γ) = u k |Ĥ 0 |u l (γ) M (i) kl . To save space in the above expression we have introduced the notation u i ∈ {a 1 , a 2 , b 1 } and the matrices In a similar way the two body operator matrix element is where d (i) (k 1 k 2 |l 1 l 2 ) (γ) is the (k 1 , k 2 , l 1 , l 2 ) entry of the second-order cofactor of d (i) (γ) and v (i) k 1 k 2 l 1 l 2 (γ) = u k 1 u k 2 |V |u l 1 (γ)u l 2 (γ) M (i) k 1 l 1 M (i) k 2 l 2 . The evaluation of these terms is lengthy but straightforward. Finally, the integrals in (A.6) and (A.7) are numerically evaluated by means of FFT once the spin projection has been performed.