Quantum phases of correlated electrons in artificial molecules under magnetic fields

We investigate the stability of few-electron quantum phases in vertically coupled quantum dots under a magnetic field of arbitrary strength and direction. The orbital and spin stability diagrams of realistic devices containing up to five electrons, from strong to weak interdot coupling, is determined. Correlation effects and realistic sample geometries are fully taken into account within the full configuration interaction method. In general, the magnetic field drives the system into a strongly correlated regime by modulating the single-particle gaps. In coupled quantum dots different components of the field, either parallel or perpendicular to the tunneling direction, affect single-dot orbitals and tunneling energy, respectively. Therefore the stability of the quantum phases is related to different correlation mechanisms, depending on the field direction. Comparison of exact diagonalization results with simple models allows one to identify the specific role of correlations.


I. INTRODUCTION
Due to three-dimensional confinement on length scales comparable to the De Broglie wavelength, the electronic properties of semiconductor quantum dots ͑QDs͒ show several similarities with those of atoms, the most significant ones being their discrete energy spectrum and the resulting shell structure. 1,2 Fine structure due to the exchange interaction ͑Hund's rule͒ [2][3][4] and Kondo physics 5,6 have also been predicted and demonstrated. QDs are therefore regarded as artificial atoms. [7][8][9][10][11] These systems have stimulated the investigation of the fundamentals of few body physics in semiconductors, since the number of electrons ͑or holes͒ in QDs can be controlled very accurately, and almost all relevant parameters influencing their strongly correlated states, such as confinement potential or the coupling with static magnetic fields, can be tailored in the experiments, driving the system between very different regimes.
The analogy between artificial and natural atoms is extended to artificial molecules ͑AMs͒ by realizing tunnelcoupled devices. [12][13][14][15][16][17] The interdot tunneling introduces a energy scale which may be comparable to other energy scales, namely the single-particle ͑SP͒ confinement energy, the carrier-carrier interaction, and the magnetic confinement energy. In AMs it is possible to tune the coupling among QDs basically at will, thus exploring different molecular bonding regimes; 18 this is an intriguing option with respect to natural molecules, since for the latter the internuclear coupling is almost fixed by the balance between nuclear repulsion and the electrostatic attraction mediated by valence electrons.
There is a number of different techniques to fabricate AMs. Laterally coupled AMs [19][20][21][22][23] are obtained by creating an electrostatic confinement in a semiconductor heterostructure, such as a doped heterojunction, by means of the photolithographic patterning of metallic gates, deposited on the surface of the heterostructure. Another way to realize AMs is by the Stransky-Krastanov mechanism leading to the formation of self-assembled QDs with nm-scale confinement with similar sizes and regular shapes. In a stack of several layers, the formation of QDs in the top layer, coupled through tunneling with the ones in the underlying layer, has been recently demonstrated. 24,25 In the following we shall investigate electronic properties of vertically coupled AMs realized starting from a triple barrier heterostructure to form a mesa pillar, possibly using a combination of dry and wet etching. 14,18,[26][27][28][29][30][31] Typically, the effective diameter for electron conduction is smaller than the geometrical diameter of the top contact. For example, with a top contact Ϸ0.5 m in diameter, it is estimated that the dot diameter is Ϸ0.1 m when the dot contains a few electrons. Still, in this class of devices the depletion potential in the lateral direction realizes a much softer potential than in the growth ͑vertical͒ direction, where electrostatic confinement is induced by the band mismatch of the layers in the heterostructure. In this vertical geometry transport occurs perpendicularly to the plane of the dots, in response to an applied voltage between source and drain contacts. The strength of the coupling can be modulated by the width of the barrier separating the QDs. In addition, a circular Schottky gate is typically placed around the body of the pillar to control the charging of the QDs.
In AMs, as well as in single QDs, electronic states can be manipulated by a vertical magnetic field B Ќ , perpendicular to the plane of the QDs ͑parallel to the growth direction͒, which drives the system from a low-correlation ͑low-field͒ regime to a strongly correlated ͑high-field͒ regime by varying the SP splittings. 9 Recently, nontrivial transitions between different quantum states induced by B Ќ and/or by varying the interdot barrier width L b have been predicted and demonstrated in AMs with different coupling regimes. 17,18,32,33 While this vertical field configuration has been widely studied both theoretically and experimentally, 15,18,[27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42]44,43,[45][46][47][48][49][50][51][52] few works have been devoted to the in-plane field configuration, 30,41,[53][54][55] which is more complicated to simulate and less intuitive, principally because of the lack of the cylindrical symmetry and of analytical solutions for the SP states. Furthermore, since in this configuration angular momentum is not a good quantum number, the computational effort increases very rapidly with the electron number ͑see Sec. II B͒.
On the other hand, in vertically coupled QDs an in-plane magnetic field B ʈ modulates the tunneling energy, and one may hope to induce transitions between few-particle states, reminiscent of those that one can obtain by varying the tunneling energy via the barrier potential. 32 In this way, different coupling regimes would be investigated within just one sample, which is of course an advantage, since it is almost impossible to grow samples completely identical except for the barrier width. It should be stressed that when the tunneling is modulated by an in-plane magnetic field in vertically coupled QDs, the other energy scales, in particular the Coulomb interaction, are only weakly affected by the field. Due to its importance to the QD-based implementation of quantum-information processing, 41 such field configuration has been already investigated in the two-electron case, 55 whereas less attention has been devoted to the case where more than two carriers are present, in particular when B has both an in-plane and a perpendicular component, which is relevant for the transport spectroscopy of such samples. In this case we expect an interesting physics: in fact both inplane ͑due to B Ќ ͒ and vertical ͑due to B ʈ ͒ correlations are involved in the transitions between quantum states.
In this paper we theoretically investigate the stability of field-induced quantum states of few electrons in AMs. The main focus is on correlation effects, largely enhanced by the field, which determine different mechanisms ͑which may interfere between each other͒ driving the transition between different ground states, depending on the field direction and interdot coupling. In our analysis we consider three realistic samples with different interdot barrier, with the aim of investigating the phenomenology occurring between the two opposite limits of strong coupling and molecular dissociation. Our numerical approach is based on a real-space description of SP states, which takes into account the complexity of real samples. Because of the competition between kinetic energy, Coulomb interaction, and Zeeman energy, it is often necessary to treat exactly ͑in contrast, e.g., to mean-field methods͒ the few-body Hamiltonian of this kind of samples. Our method of choice to include carrier-carrier Coulomb interaction is the full configuration interaction ͑FCI͒ approach, 56 which proved to be accurate and reliable, 17,18,57,58 although limited to few electrons. In order to stress correlation effects, we contrast our numerical results with SP and Hartree-Fock predictions of the stable quantum phases. We will see that correlation effects play a crucial role in determining the electronic properties of the system, and that schematic pictures neglecting these phenomena do not lead to correct predictions.
The paper is organized as follows: In Sec. II we illustrate our numerical approach to the calculation of SP and interacting states; in Sec. III we discuss the results concerning different samples in a magnetic field of arbitrary direction; and in Sec. IV we summarize our findings.

II. THEORETICAL MODEL
We consider N interacting electrons in an AM structure. Since we are interested in weakly confined devices, the energy region of our concern is relatively close to the semicon-ductor band gap and we may use the single-band approximation for the conduction band. In this framework the carriers are described by the effective-mass Hamiltonian Here m * , , and g * are the effective mass, dielectric constant, and g-factor, respectively, B is the Bohr magneton, S is the total spin, A͑r͒ is the vector potential at position r, and B is the magnetic field, with B = ٌ ϫ A. Equation ͑1͒ neglects nonparabolicity effects, but it includes otherwise in the confinement potential V͑r͒ the full three-dimensional nature of quantum states in realistic samples, such as layer width, tunneling, and finite band offsets. In the following the potential V͑r͒ describes two identical vertically coupled disk-shaped QDs. Since in this kind of samples the confinement is much tighter in the growth direction z than in the xy plane, we separate the potential as V͑r͒ = V͑x , y͒ + V͑z͒, where V͑z͒ represents two identical square quantum wells of width L W , separated by a barrier of width L b and conduction band mismatch V 0 . We perform the usual choice of a parabolic in-plane confinement of natural frequency 0 : Adequacy of Eq. ͑2͒ has been demonstrated by both theoretical calculations 59,60 and far infrared spectroscopy experiments. 61,62 Note, however, that our numerical approach does not assume any symmetry and may treat different confinement potentials; in particular, the vector potential A͑r͒ is not limited to describe the z-directed field. We suppose that the active region of the AM is sufficiently well-separated from the top and substrate contacts that we can reasonably neglect interactions with the leads in our model. 18 The energy splitting ប 0 between SP eigenvalues of the in-plane confinement, 10 given by Eq. ͑2͒, is an input parameter of the model. Such quantity mimics the effects of electrode screening and interaction with the environment, and in real charging experiments it is modified by the Schottky gate. 18 Here, for the sake of simplicity, ប 0 is kept fixed as N is varied. A qualitative discussion of this issue is postponed to the end of Sec. III. Schematically, our algorithm is the following. First, we numerically calculate the electron SP states, mapping the SP Hamiltonian on a real-space grid. Then, the SP orbitals thus obtained are used to evaluate Coulomb matrix elements and finally to represent the interacting Hamiltonian on the basis of Slater determinants ͑SDs͒, according to the FCI approach. These steps are described below in detail.

A. Single-particle states
The SP energies ␣ and wave functions ␣ ͑r͒ are obtained from the numerical solution of the eigenvalue problem associated to the Hamiltonian which is the SP term appearing in the first line of Eq. ͑1͒. We include in Eq. ͑3͒ a magnetic field of arbitrary direction via the vector potential A = 1 2 ͑B ϫ r͒ ͑symmetric gauge͒. Specifically, we find the eigenstates of Eq. ͑3͒ by mapping it on a real-space grid of N grid = ͟ i=1 3 N i points, identified by the grid vectors .. ,N k and ê 1,2,3 = x , ŷ , ẑ. The resulting finite-difference equation can be rewritten in terms of a discrete eigenvalue problem, as described in detail elsewhere. 63 This results in the diagonalization of a large, sparse matrix which is performed by the Lanczos method. 56

B. Few-particle states
Once SP states are obtained, we proceed including the Coulomb interactions between charge carriers. The full many-body Hamiltonian H can be rephrased, by means of the second-quantization formalism, 64 as where Ĥ 0 is the single particle Hamiltonian and Ĥ Z is the Zeeman Hamiltonian Here, ĉ ␣ † ͑ĉ ␣ ͒ creates ͑destroys͒ an electron in the spinorbital ␣ ͑r , s͒ = ␣ ͑r͒ ͑s͒, where ͑s͒ is the spinor wave function, and Ŝ is the total spin vector. In the above equations, N SP is the number of SP states that are taken into account.
Since typical SP energy spacings computed from the Schrödinger equation associated with Eq. ͑3͒ are in the meV range, while characteristic values of the Zeeman term Ĥ Z are typically two orders of magnitude smaller ͑e.g., B B = 5.79 ϫ 10 −2 meV at 1 T͒, the effect of Ĥ Z is first neglected in the calculation of the electronic ground states ͑GSs͒. Then, the effect of Ĥ Z , which just lifts the ͑2S +1͒-fold degeneracy of the total spin S, is included as a perturbation to the first order in the field.
Then, our algorithm proceeds forming a basis of SDs, filling N SP spin-orbitals ␣ ͑r , s͒, calculated numerically as described in Sec. II A, with N electrons in all possible ways. We assume that the Fock space generated by the SDs is approximately complete, namely a generic eigenstate of the many-body Hamiltonian Ĥ can be expressed as a linear combination of SDs. One may reduce the computational effort by exploiting the symmetries of the system. In fact, by using a suitable combination of the SDs which diagonalize symmetry-related operators, the Hamiltonian can be written in a block-diagonal form. 56 The larger the symmetry of the system, the smaller the dimensions of the blocks and the time necessary to complete the full diagonalization. In the present case, an in-plane component of the magnetic field removes the cylindrical symmetry C ϱh , which is often exploited in the vertical field arrangement of cylindrically symmetric QDs. However, the Hamiltonian still commutes with the square total spin Ŝ 2 and z-component Ŝ z operators, respectively. Therefore the subspaces are appropriately labeled by the values of S and the minimum positive value of S z consistent with S. Once the effective subspace is selected, our code performs a unitary transformation in order to rewrite Ĥ in a block-diagonal form. Finally, matrices obtained in this way are handled by a Lanczos routine included in the package DONRODRIGO 65 run on a SP3 IBM system.

III. RESULTS
In this section we present the calculated field-dependent GSs of AMs, covering the regimes from strong to weak coupling. The application of the FCI method ensures that our results have comparable accuracy in different regimes. To this end, we have selected three realistic samples based on a AlGaAs/ GaAs heterostructure, labeled B1, B2, and B3, with interdot barriers L b = 2.5, 3.2, and 4.0 nm, respectively. The sample with intermediate barrier shows a full molecular character, while B1 and B3 represent the single QD and the molecular dissociation limits, respectively.
Except for the barrier width, parameters are common to all samples. The lateral confinement is ប 0 = 3.5 meV, the band-offset between barrier and well materials of the heterostructure is V 0 = 300 meV, each well width of the double quantum well potential is L W = 12 nm. Material parameters appropriate to GaAs are m * / m e = 0.067, = 12.9, and g * = −0.44. All parameters refer to a set of samples 18,27 described in literature and processed at the University of Tokyo. A. Single-particle states SP states in an arbitrary field have been obtained according to Sec. II A with N 1 = N 2 = 80, N 3 = 128; the resulting 819 200ϫ 819 200 sparse Hamiltonian matrix is diagonal-ized by means of the Lanczos algorithm. The grid parameters have been chosen in order to reproduce the lowest analytical zero field ͑Fock-Darwin͒ energies nm with an accuracy better than 0.5 meV.
Before presenting our results for specific samples, let us summarize the properties of SP states in a single QD with parabolic in-plane confinement and vertical magnetic field B Ќ . In this case SP states are analytical and given by the so-called Fock-Darwin ͑FD͒ states, with energies nm = ប⍀͑2n + ͉m͉ +1͒ − ͑ប c /2͒m, n =0,1,2,... and m =0, ±1, ±2,... being the principal and azimuthal quantum numbers, respectively. The effective oscillator frequency ⍀ is defined by ⍀ = ͱ 0 2 + c 2 / 4, and c = eB / m * c is the cyclotron frequency. A perpendicular magnetic field splits the degeneracies of the states and reduces the energy gap between energy levels with increasing angular momentum ͑see solid lines in Fig. 1͒; this favors transitions in the few-electron states in order to increase in-plane correlations, as explained in the next section. In symmetric AMs, FD levels are replicated, rigidly shifted by ⌬ SAS , the energy gap between the lowest symmetric ͑S͒ and antisymmetric ͑AS͒ states arising from the double-well potential along the growth direction ͑highest confined states along the growth direction need not be considered if L W is sufficiently small͒. 18 Figure 1 shows the S ͑solid lines͒ and AS ͑dashed lines͒ replicas of FD states for the three samples as a function of B Ќ . Note that the S/AS labeling is valid only as far as the field is entirely along the z axis, namely the motions in the xy plane and along z are completely uncoupled, and the splitting ⌬ SAS does not depend on B Ќ . The symmetry group, in this highly symmetric case, is C ϱh . Figure 2 shows calculated SP levels for B1, B2, and B3 as a function of the in-plane field B ʈ . Contrary to the verticalfield configuration, the parallel field B ʈ , which we take to lie along the x axis, couples the motion along the in-plane y and the vertical z directions, and it reduces the symmetry group of the system from C ϱh to C 2h . Therefore SP wave functions lose their well-defined component of the orbital angular momentum as well as the S/AS character along z. Besides, the states labeled S and AS when B ʈ = 0, respectively, get closer as B ʈ increases, and tunneling is progressively suppressed. Clearly, the tunneling suppression occurs at different values of the field for different samples, since the in-plane field may significantly affect tunneling only when c ʈ = eB ʈ / m * c ϳ ⌬ SAS / ប. 55 Although the application of B ʈ strongly reduces the symmetry group of the system, it still preserves two residual symmetries: ͑i͒ the reflection x → −x with respect to the yz plane perpendicular to the field direction and ͑ii͒ the rotation of 180°around the x axis. The four irreducible representations of the residual group C 2h , namely A u , A g , B u , and B g , 66 are associated to the appropriate SP states in Fig. 2, together with the alternative labeling at zero field. It is possible to understand qualitatively the behavior of the SP states with the parallel field B ʈ by taking into account the energy gap between the orbitals with the same symmetry; the smaller the gap, the larger the "repulsion" between levels. For example, we expect strong repulsion between the second and the fourth state, and a weaker effect for the first and the fifth state because of the larger energy gap between them.
The effect of the in-plane field on SP states is demonstrated in Fig. 3, which shows the SP charge densities ͉͑x , y͉͒ 2 at a fixed value of z inside one well for the three lowest levels ͑their quantum numbers are indicated in the figure͒, for three selected values of B ʈ . The field increases the confinement, squeezing the electronic wave functions. Moreover, FD orbitals with different angular momentum and belonging to different energy shells are mixed by B ʈ , leading to charge density modulation and nodal surfaces.
In presence of a magnetic field of arbitrary direction ͑tilted magnetic field͒ also the C 2h symmetry is broken, and we observe simultaneously the reduction of both tunneling and intradot energy gaps between FD states as the field is increased. 55 An example of the SP states' behavior in a magnetic field which is tilted at = 45°with respect to the growth direction is reported in Fig. 4.

B. Few-particle states
In the following, we will investigate the stability of fewparticle ͑2 ഛ N ഛ 5͒ quantum phases as a function of an applied field with both in-plane and vertical components. In AMs, due to tunneling and the ensuing formation of bonding and antibonding levels, the number of SP orbitals which ensures convergence of the FCI calculation depends on the splitting, ⌬ SAS , and on the number of charge carriers, N, we take into account. Results presented in this work are obtained using N SP = 20 for 2 ഛ N ഛ 4. In order to limit the computational effort, which increases very rapidly with the number of electrons, we reduce to N SP = 15 the set of SP states for N = 5. With these values, the calculated SP states, which at arbitrary values of the field are in general nondegenerate, turn continuously into close zero field degenerate shells. The convergence of the FCI calculations at zero field with respect to the selected SP states has been carefully investigated in Ref. 56.
In order to stress the role of correlations, we will contrast the exact correlated states, obtained by the FCI method, with ͑i͒ the quantum phases predicted in a SP picture, i.e., filling SP states according to the Aufbau principle, and ͑ii͒ a Hartree-Fock-like approach.
In the SP picture the few-electron noninteracting GS energy is where ␣ is the SP energy of the ␣th orbital and n ␣ is the occupation number of the ␣th SP orbital with spin . The second model is reminiscent of the Hartree-Fock approach: In this approximation, we calculate the energy of an electronic configuration as the expectation value of the interacting Hamiltonian on the most weighted SD singled out from the FCI expansion of the correlated wave function. We call this the single SD ͑SSD͒ approximation. The SSD energy is given by Here the interaction between charge carriers is described by the integrals U ␣␤ = ͵͵ ͉ ␣ ͑r 1 ͉͒ 2 e 2 ͉r 1 − r 2 ͉ ͉ ␤ ͑r 2 ͉͒ 2 dr 1 dr 2 , ͑12͒ where U ␣␤ is the direct Coulomb integral, accounting for the repulsion between two electrons occupying orbitals ␣ and ␤, and K ␣␤ is the exchange integral, giving the exchange interaction between electrons with parallel spins. U ␣␤ and K ␣␤ correspond to V ␣␤␤␣ and V ␣␤␣␤ in Eq. ͑8͒, respectively. Note that in both SP and SSD models the GS is given by just one configuration, while, in general, the true GS is a linear combination of different electronic SDs, due to correlation. We also neglect the Zeeman coupling to the field at this level. Below we shall discuss the GS of few electrons as a function of the magnetic field intensity and direction. In several cases, we will indicate schematically ͑see the following fig-ures͒ the GS configuration, as predicted either by the FCI method or by the SSD calculation, in terms of arrows, pointing either upwards or downwards, filling in either left or right boxes ͑with respect to a diagonal line͒. The arrows stand for electron spin, while superposed boxes indicate SP orbitals of increasing energy. In the vertical-field configuration, left and right boxes indicate S and AS levels, respectively ͑right boxes are not shown if no AS states are occu-pied͒. If a finite B ʈ component is present, the boxes indicate the orbitals which evolve in a continuous manner from those at B ʈ =0 as B ʈ is switched on ͑see Fig. 2͒. For example, the lowest-energy S ͑AS͒ s and p levels at B ʈ = 0 evolve continuously into A g and A u ͑B u and A g ͒ levels, respectively, as B ʈ increases, while remaining identified pictorially by the same boxes. For FCI calculation, we show only the most weighted configurations. In some cases FCI and SSD methods predict the same stable phase, in which case we use dotted boxes. When the two methods disagree, we use solid boxes and indicate the exact ͑FCI͒ result only. Typically, for a given sample and number of carriers, the set of stable GSs predicted by the SSD method is a subset of the FCI stable phases.
In all stability diagrams shown below, calculated GSs are indicated with dots, and lines are only a guide to the eye through the calculated points. Due to numerical limitations, mainly related to the lack of cylindrical symmetry which requires numerical calculation of the SP states, we have limited our calculations to a relatively coarse grid of points in the ͑B Ќ , B ʈ ͒ plane. Although we believe that we have determined most of the GSs which are stable in the different regimes, it is possible that we miss some phases which are stable in a small range of fields.
We first consider the two-electron case. 55 In a SP picture, the lowest SP state is doubly filled in a singlet ͑S =0͒ configuration which is stable in all field regimes, since no level crossing occurs in the lowest SP state ͑see Figs. 1 and 2͒. Figure 5 shows, instead, that FCI predicts a singlet/triplet transition at finite fields, both along the B Ќ and B ʈ axes. For B2 the SSD phase boundary is also reported in Fig. 5 ͑no qualitative differences are expected for the other samples͒.
Let us focus for the moment on the SSD stability diagram of Fig. 5. At low magnetic field the ground state has a singlet character, with two electrons sitting on the same orbital. As the field is increased, we observe a singlet/triplet transition, which is different in character depending on whether B Ќ or B ʈ is varied. 55 In the former case, only S levels s and p ͑equivalent to A g , A u levels, respectively͒ are involved: the triplet configuration minimizes the Coulomb interaction by promoting an electron from a s to a p orbital, which is more spatially delocalized, i.e., U ss Ͼ U sp . In this sense, only inplane degrees of freedom play a role in the B Ќ -induced transition. The "Hartree" energy gain U sp − U ss compensates the cost in terms of SP energy due to the s → p orbital promotion. Moreover, the exchange interaction energy gain, given by the exchange integral K sp , favors the "ferromagnetic" triplet configuration ͑S =1͒ with respect to the "antiferromagnetic" singlet when the two electrons sit on the s and p orbitals, respectively. Therefore the electrostatic energy of the triplet configuration ϷU sp − K sp is further reduced with respect to the singlet one, ϷU sp .
In the in-plane configuration, instead, the tunneling energy separating the lowest-energy A g and B u orbitals is very effectively reduced by B ʈ ͑see Fig. 2͒. Both orbitals are therefore involved in the formation of singlet and triplet states. In the SSD approximation singlet and triplet states, the former made of a doubly occupied A g orbital and the latter of two parallel-spin electrons sitting on the A g and B u orbitals, respectively, tend to the same orbital energy. However, the small exchange interaction K A g B u and the Zeeman term of the S = 1 state, both included in the SSD calculation, weakly favor the triplet when the tunneling energy is reduced above a critical value of B ʈ . Although the SSD scheme is able to describe, on a qualitative level, the singlet/triplet stability regions, neglecting correlation effects in the SSD scheme strongly underestimates the stability region of the singlet state, favoring the polarized phase: in fact, the description of the singlet state in terms of two electrons sitting on a A g orbital is poor especially at high field, where an increasing contribution from the B u orbital gives rise to a spatial correlation along the growth direction, i.e., with the two electrons sitting on opposite QDs, as we have shown in a previous work. 55,67 The vertical correlation lowers the energy of the singlet state, making it stable in a larger B ʈ -region, although the Zeeman coupling favors the triplet configuration. Indeed, the stability of the singlet state in FCI may be interpreted as resulting from a "superexchange" effect. 68 The FCI singlet configuration can be better understood after a unitary transformation of the orbital basis to states localized on either one or the other dot. By schematizing the AM as a two-site lattice at half filling, the appropriate conceptual framework, when the QDs are weakly coupled, is that of the Hubbard model, 69,70 where the system has a tendency towards "antiferromagnetism," since two electrons in a singlet state can lower their energy by virtually hopping from one to the other dot, while the same process is prohibited for a triplet state by Pauli blocking. When tunneling energy does not pay off for antiparallel alignment, the triplet phase becomes favorable.
On the other hand, a SSD scheme would neglect also the correlation regarding the triplet state when both B Ќ and B ʈ are present, namely the central region in the diagrams of Fig.  5. Here the triplet state is given by the sum of mainly two components, the first one involving A g and A u orbitals, and the second one involving A g and B u orbitals, with coefficients depending on the field strength and direction. In this case the effects induced by both B ʈ and B Ќ compete in a nontrivial manner, as well as vertical and in-plane correlations. 71 Finally, we note that the stability region of the singlet state with respect to B ʈ depends on the interdot coupling, and it decreases going from B1 to B3. In fact, as we have seen before, the in-plane field affects significantly the tunneling only when c ʈ = eB ʈ / m * c is comparable with ⌬ SAS / ប. 55 Consistently with the above interpretation of the singlet/triplet transition, the critical value of B Ќ depends on tunneling strength only very weakly. We next consider the N = 3 case. In the SP noninteracting scheme ͑Fig. 6͒ only two phases are present, both with S =1/2. The field B Ќ at which the two GS energies cross is given by the intersection of the S s ͑ϵA g ͒ and p ͑ϵA u ͒ levels ͑Fig. 1͒. Figure 7 shows the stability diagram obtained from the FCI calculation; for B2 the phase boundary obtained from the SSD approach is also shown ͑dashed line͒. Although the SSD approximation seems to describe reasonably well the dynamics along the B Ќ axis, where GS configurations coincide with the FCI ones, except a slight shift of the phase boundary at lower fields with respect to FCI, this scheme fails completely at finite B ʈ . Indeed, SSD severely underestimates the stability region of the unpolarized phase ͑S =1/2͒ with respect to the polarized one ͑S =3/2͒. This fact is mainly connected to correlation effects, as we shall see in the following.
Let us stick for a moment to the vertical field configuration. As already for N = 2, also for N = 3 increasing B Ќ determines an enhancement of the in-plane correlation, leaving unaffected the motion along the growth direction. For B1 and B2, electrons in the S orbitals undergo a number of transitions where p and d levels are successively populated in order to minimize the Coulomb repulsion. Eventually, the so-called maximum density droplet ͑MDD͒, [72][73][74][75] i.e., the densest spin-polarized configuration possible, is reached. For B3 this phase is not observed at B ʈ = 0 in the range of B Ќ explored ͑Fig. 7͒. If B ʈ is small or zero, electrons are prevented from occupying AS orbitals by the kinetic energy cost ⌬ SAS . As B ʈ is increased, however, ⌬ SAS decreases. The range of fields where AS occupation takes place depends, of course, on the zero field tunneling energy, as it is easily seen comparing the GS character at low B ʈ of the B3 sample, where electrons occupy both S and AS states, with that of B1 and B2 samples, where AS states are not occupied. From this point of view, increasing B ʈ is analogous to decreasing L b . We will find a similar behavior for N =4,5.
On the other hand, in addition to the reduction of the tunneling splitting ⌬ SAS , a finite value of B ʈ determines a redistribution of the carrier density, coupling the y and z degrees of freedom, reducing the spatial symmetry, and inducing orbital hybridization, as discussed in Sec. III A. From this point of view, increasing B ʈ is not equivalent to decreasing L b , since it induces GS configurations different from those occurring at large L b . Moreover, as we show below, the high-B ʈ phases are brought about by genuine correlation effects, and in fact they are missed by the SSD approximation.
Let us consider the spin phase at high B ʈ and low B Ќ in Fig. 7, having the B u orbital doubly occupied in its mostweighted SD. At high B ʈ , the orbitals A g , B u are almost degenerate ͑see Fig. 2͒. Therefore we would expect the two configurations with two electrons either in the A g or in the B u level to occur with like probabilities. However, intraorbital Coulomb repulsion is lower for the more delocalized B u state. This is shown for B3 in Fig. 8͑b͒, where we report the behavior of U ␣␣ ͓defined in Eq. ͑12͔͒, where ␣ takes the value A g and B u , respectively. The modulation of U A g A g and U B u B u with B ʈ should not come as a surprise, since the SP orbitals are strongly affected by the field ͑Fig. 3͒, with the appearance of nodal surfaces even in the lowest-energy SP orbital A g . These integrals allow one to estimate the value of the Coulomb repulsion between two electrons both sitting in the A g or in the B u orbital, as a function of the in-plane field B ʈ . In Fig. 8͑a͒ we show the ground-and lowest excited-state energies for B3 vs B ʈ . The larger values of U A g A g at intermediate fields partly explains the occurrence of the SD with two electrons sitting on the B u orbital. On the other hand, such GS is not observed in a SSD picture, i.e., neglecting correlation, as shown in Fig. 7 ͑note that the B ʈ field range is larger in Fig. 8 than in Fig. 7͒. At high values of B ʈ , U A g A g Ӎ U B u B u and the state with two particles on the B u level becomes almost degenerate with the one with two particles on the A g level ͓Fig. 8͑a͔͒. The results for B1 and B2 samples ͑not shown͒ are similar.
The reduced Coulomb repulsion between electrons on the B u orbital is a key concept to understand the stability of the unpolarized phase S =1/2. Indeed, the SD with the A g orbital singly occupied and the B u orbital doubly occupied is always present in the GS when applying B ʈ ͑even though it is not indicated in the graphs unless it is the dominating one͒, together with other important configurations. The reason is that correlation, namely the mixing of different SDs, allows the system to lower its energy. At high B ʈ , the SD mentioned above becomes the dominating one for all considered samples, although close in energy to the SD with the A g orbital doubly occupied. As in other cases, in tilted magnetic fields the typical SDs appearing in the GS expansion, in the central region of the diagrams, represent the competition be- tween vertical correlation, driven by B ʈ , and in-plane correlation, driven by B Ќ .
For N = 4 the situation is more complicated due to the larger number of possible phases, and the strong dependence on the interdot coupling strength. SP, SSD, and FCI predictions are presented separately in Figs. 9-11, respectively. The predictions of the SP scheme in Fig. 9, if compared with those of FCI of Fig. 11, are clearly inadequate to describe the physics of the system: only two noninteracting singlet phases are stable, differing for the double occupancy of either the A u or the B u orbital. 76 As expected, the FCI calculation ͑Fig. 11͒ predicts the familiar MDD state formation through successive occupation of p and d orbitals; this is also predicted, even though approximately, by the SSD approximation ͑Fig. 10͒. On the other hand there are many differences between SSD ͑Fig. 10͒ and FCI ͑Fig. 11͒ predictions when B ʈ is also increased. One major difference is the enhanced stability of the S = 0 phase in FCI calculations. This phase involves the SD made of doubly occupied A g and B u orbitals, which is the dominant one in a large field range for all samples. Analogously to the N = 3 case, the system is allowed to lower its energy through occupation of the B u level as B ʈ is increased. The true GS, however, turns out to be a nontrivial linear combination of several other SDs, which suggests that correlations are dominant here. Indeed, by neglecting correlation effects, like in the SSD approach, the only possible energy gain is obtained by aligning the spins in a S = 1 phase, leading to the wrong SSD prediction that the S = 1 phase is always stable at small B Ќ , as shown in Fig. 10.
Such partially polarized phases, instead, are typically stable only in small ranges of the field, and their stability also depends on the sample parameters, as we show below. For N = 4 the S = 1 phase at B Ϸ 0 is stable only in a small field range, and its stability is reduced if the interdot cou-pling is reduced; in sample B3 this phase is not stable, no matter how small the field is ͑Fig. 11͒. Indeed, this phase requires occupation of two S p orbitals, and it results from the payoff between the gain in exchange energy and the kinetic energy cost. Therefore this configuration becomes rapidly unfavorable with respect to the S = 0 phase having the B u orbital doubly occupied, as B ʈ increases. In fact, the stability range of the S = 1 phase decreases as the zero field tunneling energy ⌬ SAS decreases with respect to the s − p energy splitting ប 0 . Note that, according to the SSD data ͑Fig. 10͒, the S = 1 phase is given by only one SD, instead of the FCI superposition of two main configurations with coefficients dependent on the value of the field. This fact, which indicates the progressive filling of the B u level while increasing the field, is neglected by the less accurate SSD picture, and it is described by the SSD model in terms of a transition between the two SDs mentioned above. We note in Fig. 11 that a S = 1 phase sets in also between the unpolarized ͑S =0͒ and the completely polarized ͑S =2͒ spin case. This phase, characterized by a nontrivial combination of configurations involving both S and AS states, loses its stability as the interdot barrier is increased.
Similar considerations apply to the N = 5 case ͑see Figs. 12-14 for SP, SSD, and FCI schemes, respectively͒. Typically, several transitions occur in the GS within the same spin configuration. For example, along the B ʈ axis of Fig. 14 several configurations occur, reflecting the crossing between p levels in the S and AS minibands ͑see Fig. 2͒. As the barrier L b is decreased, the splitting between s and p orbitals becomes of the same order of magnitude of ⌬ SAS , and the  number of transitions in the GS increases too. Also in this case a SSD approach is inappropriate. As in the N = 4 case, we have an intermediate phase ͑S =3/2͒ between the unpolarized ͑S =1/2͒ and the completely polarized ͑S =5/2͒ ones, which is only stable in the intermediate coupling regime; this phase is absent in the SSD results ͑Fig. 13͒.
The above results were obtained by keeping ប 0 fixed with respect to N. In single-electron tunneling experiments, however, ប 0 changes with the gate voltage, in order to inject ͑extract͒ electrons into ͑from͒ the AM. 18 A relation linking 0 with N may be derived by imposing that the electron density is approximately constant during the charging processes, 0 ϰ N −1/4 ͑Refs. 9 and 18͒. This allows us to qualitatively understand what changes the AM phases undergo if we try to mimic the dependence of 0 on N. Indeed, if ប 0 decreases as N increases, the ratios of this energy to the other relevant energy scales-namely ⌬ SAS , ប c , and Coulomb matrix elements, respectively-decrease as well. This implies that ͑i͒ the AM phase boundaries in the ͑B ʈ , B Ќ ͒ plane move in order to "squeeze" the stability regions of "ionic" phases; ͑ii͒ the transition to the MDD phase along the B Ќ axis is faster; and ͑iii͒ correlation effects are stronger.
Finally, we note that the present results for N =4,5 at B ʈ = 0 are in agreement with the ones shown by Rontani et al. dealing with similar samples, 18 except for little discrepancies due to different parameters used in the calculations, in particular the band-offset V 0 and the lateral confinement energy ប 0 .

IV. CONCLUSIONS
We have investigated the relative stability of quantum phases of interacting electrons in AMs under a magnetic field. Such phases can be identified, e.g., in transport experiments in the Coulomb blockade regime, where the addition energy is measured. 2 A very sensitive test of few-electron phases is nonlinear transport through charged QDs, where a few excited states in a small energy window are accessed. 17,18,77 In addition to predict the stability regions of the different spin phases, this work confirms that correlation effects play a fundamental role in determining the physics of the system, and a simple picture neglecting these effects could fail in interpreting the relevant experiments.
In particular, we have compared FCI predictions with ͑i͒ a SP picture, neglecting all interactions, and ͑ii͒ a SSD scheme, where only direct and exchange interactions between charge carriers are taken into account. These comparisons allowed us to identify the effects of correlations, particularly those induced in the coupled structure by an inplane magnetic field, in analogy to what happens in single QDs with an applied vertical field. In many respects, only the exact diagonalization of the few-body Hamiltonian is able to provide a correct description of correlated states.