Coulomb and tunneling coupled trilayer systems at zero magnetic field

The ground-state electronic configuration of three coupled bidimensional electron gases has been determined using a variational Hartree-Fock approach, at zero magnetic field. The layers are Coulomb coupled, and tunneling is present between neighboring layers. In the limit of small separation between layers, the tunneling becomes the dominant energy contribution, while for large distance between layers the physics is driven by the Hartree electrostatic energy. Transition from tunneling to hartree dominated physics is shifted towards larger layer separation values as the total bidimensional density of the trilayers decreases. The inter-layer exchange helps in stabilize a"balanced"configuration, where the three layers are approximately equally occupied; most of the experiments are performed in the vicinity of this balanced configuration. Several ground-state configurations are consequence of a delicate interplay between tunneling and inter-subband exchange.


I. INTRODUCTION
Single-layer quasi two-dimensional electron gases (2DEG) such as those formed at the interface between two dissimilar semiconductors can be routinely driven to a many-body interaction regime by application of a strong magnetic field of several Teslas perpendicular to the 2DEG layer [1]. By increasing the magnetic field, the 2DEG first enters in the Integer Quantum Hall regime, and then in the Fractional Quantum Hall regime [2].
Multilayer coupled 2DEG's like trilayer systems offer much more possibilities for the theoretical [3][4][5] and experimental [6][7][8][9][10][11][12] search of the involved physics, even at zero magnetic field. New single-particle effects like layer (site) energies, tunneling coupling between layers, and many-body effects as the inter-layer Coulomb coupling should be included, with the two later effects being strongly dependent on the separation between layers.
We provide in this work an exhaustive theoretical study of trilayer systems at zero-magnetic field, within the framework of a Hartree-Fock mean field approximation, already used for bilayer systems [13], and also for a simplified version of the trilayer system considered here [4]. The model, schematically illustrated in Fig. 1, includes intra-layer and inter-layer (hopping) kinetic energies (not considered in Ref. [4]), site (layer) energies (also not considered in Ref. [4]), intra-layer and inter-layer exchange energies, and the usually dominant Hartree electrostatic energy. An additional external parameter is the total density of the system: we will see that as the density decreases, the trilayer system is driven to a single-particle tunneling dominated regime.
By doing a numerical minimization of the total energy of the system over all the internal variables, we have * dmiravet@gmail.com † proetto@cab.cnea.gov.ar found a zero-temperature extremely rich phase-diagram for the possible ground-state configurations. These full numerical results may be qualitatively understood as resulting from the diverse scaling properties of the different energy contributions to the total energy of the system, with respect to the external variables like the total density or the distance between layers. Interestingly, we have found that at zero magnetic field and for close enough layers the physics of the trilayer system becomes dominated by single-particle effects both in the high-density limit (dominated by the intra-layer kinetic energy) and in the low-density limit (dominated by the hopping or tunneling and site energies). The interaction dominated regime is then restricted to intermediate densities (exchange) and large separations between layers (Hartree). The wide ground-state exploration in the available parameter space (distance between layers, total electronic density, tunneling coupling) provided in this work may serve as a qualitative guide for the design of experimental trilayer systems. The rest of the work is organized in the following way: in Section II we explain the model and introduce the variational Hartree-Fock method we use to obtain the expression for total energy of the system. Analytical and numerical ground-states resulting from the minimization of the system total energy are presented in Section III. Section IV is devoted to the conclusions. In the two appendices we discuss some limits and the physics beyond the inter-layer exchange energy contribution (Appendix A), and we explain how to model the hopping parameter dependence with the layer separation (Appendix B).

II. MODEL
Typically, a semiconductor trilayer system is realized experimentally by confining three GaAs quantum wells among AlAs or Al x Ga 1−x As barriers of variable width and height, which gives a control on the tunneling coupling t among layers. Usually, the central GaAs quantum well is designed with a larger width than the two side-wells with the purpose of populate the central well, which tends to be less populated than the two sidewells [14,15]. In our model of Fig. 1, we can simulate this feature by imposing, for instance, that ε 2 < ε 1 = ε 3 .
The model for our Coulomb and tunneling coupled trilayer system is represented schematically in Fig. 1, and defined through the corresponding Hamiltonian in Eq. (1). It consists of three strictly two-dimensional metallic layers, at a distance d between them. Two other layers, located at z = ± h along the z-axis (h d) provide the compensating positive charge densities. Neighboring metallic layers are coupled through the tunneling term t (hopping), while all layers are Coulomb coupled through the inter-layer Coulomb interaction. Schematic view of the three layer system. Each layer of area A has a two-dimensional electronic number density given by n1, n2, and n3. Charge neutrality is provided by two positively charged layers located at z = ± h (not shown), with h d. The layers are coupled by the hopping t and the Coulomb interaction. ε1, ε2, and ε3 are the site or layer energies.
In detail, the Hamiltonian of our model is as follows, with A denoting the area of each layer. Here, c † jkσ (c jkσ ) is a creation (annihilation) operator for an electron in layer j (j = 1, 2, 3), with two-dimensional momentum k and spin σ (↑ or ↓). Each one of the five terms in the Hamiltonian represents a different physical contribution to the total energy. The first corresponds to the sum of the layer and kinetic energy of electrons in each layer, with ε jk = ε j + 2 k 2 /2m * (m * being the electron effective mass of the well-acting semiconductor). The second term represents the quantum tunneling of electrons among different layers, with a hopping amplitude parametrized by t (> 0); j 1 = j 2 , and the sum is restricted to neighboring layers. The last three Coulombrelated terms represent the attractive ion (positive layer) -electron interaction, the ion (positive layer) -ion (positive layer) interaction, and the repulsive electron -electron interaction, respectively. p m (m = L, R) denotes the uniform density of the two layers located at left and right of the three central layers. The system fulfills a global neutrality condition, N = n 1 + n 2 + n 3 = p L + p R , where N represents the total two-dimensional number density, and n j (j = 1, 2, 3) are the number densities of each metallic layer. Finally, V ij (q) = 2πe 2 q e −qdij , with d ij = 0, d, h; is the dielectric constant of the wellacting semiconductor (∼ 12.5 for GaAs). For t = 0 and ε 1 = ε 2 = ε 3 , the model reduces to the one studied previously in Ref. [4]. The layer energies are changed at will in real samples by the application of back and front gates, which in our case are represented by the two positively charged layers at ± h. In this work we consider only the symmetric configurations, i.e. p L = p R , and The presence of the last electron-electron interaction term in Eq. (1) implies that no exact solution is available for the model, and forces us to attempt its approximate solution. For this, we will employ a Hartree-Fock variational approximation, widely used for the bilayer case, either at zero [13] or with magnetic field [16][17][18], and also used in the previous trilayer zero-tunneling and zero-site energy study [4]. But before that, we find convenient to perform an exact transformation of the Hamiltonian, by defining the following operators: with 0 ≤ θ ≤ π, and 0 ≤ φ < 2π. The new operators satisfy fermion anticommutation relations {α kσ , α † k s } = δ k,k δ σ,σ , while all others anticommutators vanishes. The transformation to the "subband" basis a † kσ , b † kσ , c † kσ from the "layer" basis c † 1kσ , c † 2kσ , c † 3kσ is defined by the two angles θ and φ. For θ = φ = 0, both basis are the same. For θ = π/2, φ = 0, a † kσ , b † kσ , and c † kσ are creation operators for electrons in the ground-(symmetric, no nodes), first-excited (antisymmetric, one node), and second-excited (symmetric, two nodes) subband states of the trilayer system. The canonical transformation of Eq. (2) may be considered as a general rotation in a 3-dimensional pseudo-spin Hilbert space of a pseudo-spin spinor pointing in then = (sin θ cos φ, sin θ sin φ, cos φ) direction, with the layer index j playing the role of the pseudo-spin components. In the ground-state, the direction ofn is determined by the total energy minimum. The coefficients in the 3 × 3 transformation matrix in Eq. (2) are obtained from the eigenstates of the 3 × 3 matrixn(φ, θ) · S = n x (φ, θ)S x + n y (φ, θ)S y + n z (φ, θ)S z , with S x , S y and S z being the components of the 3 × 3 angular momentum matrices corresponding to unit pseudo-spin (or angular momentum) [19]. Write in the a, b, c basis, the hopping term in Eq. (1) becomes diagonal for the choice φ = 0, θ = π/2, but our variational ansatz given by Eq. (3) below is more flexible and allows φ and θ to take any value within their permissible range. On the other side, the transformation of Eq. (2) is not the more general one, considering that the angles φ and θ may be, in principle, σ and k dependent. For simplicity, we have not included these dependences in our calculations.
After expressing the Hamiltonian in Eq. (1) in term of the subband operators, we have taken the expectation value of the transformed Hamiltonian with the following Hartree-Fock variational ansatz for the ground statevector [20], Here, k aσ , k bσ , and k cσ are the Fermi wavectors for electrons with spin σ in subbands a, b, and c, respectively. If any of the six k ασ = 0, this means that the corresponding subband is empty. We obtain, after a lengthy calculation where and Here, r s = 1/(a * 0 √ πN ) is the dimensionless twodimensional density parameter, a * 0 = 2 /e 2 m * is the effective Bohr radius of the well-acting semiconductor, and Ry * = m * e 4 /(2 2 2 ) is the effective Rydberg. For GaAs as well-acting material, m * 0.067m 0 with m 0 denoting the bare electron mass, and 12.5, resulting in Ry * 5.83 meV, and a * 0 98.7Å. Also, η α = σ η ασ , η ασ = k α † kσ α kσ /AN are the total and spindiscriminated subband occupation factors, d * = d/a * 0 is the distance from the central layer to the two lateral layers (in units of the effective Bohr radius a * 0 ), t * = t/Ry * , and ε * i = ε i /Ry * . The expression for the quantities I αβσ (q) is given in the Appendix A. These exchange integrals correspond to the overlap area between two Fermi circles, with radius k α,σ and k β,σ , with the circles centers separated by a distance q. The overlap area is maximum for q = 0, and decreases as q increases; when q ≥ k α,σ + k α,σ the superposition area (and the associated exchange integral) becomes zero. All energies in Eq. (4) are given in units of N Ry * , i.e., in units of energy per unit area. As defined above, all 0 ≤ η ασ ≤ 1, and ασ η ασ = 1.
E K 0 corresponds to the sum of the intra-layer kinetic and layer energy contributions. For the case ε * 1 = ε * 2 = ε * 3 , the layer energy contribution reduces to an uninteresting constant term, which only depends on the total electron density [4]. E T 0 is the tunneling or interlayer kinetic energy, and is the only term where the angle φ appears. E H 0 is the Hartree electrostatic energy, and E X-intra 0 and E X-inter 0 are the intra-subband and intersubband exchange-energy contributions, respectively. For a given set of external parameters (r s , , the ground-state energy E 0 depends on eight variational parameters: η a↑ , η a↓ , η b↑ , η b↓ , η c↑ , η c↓ , θ, φ. However, the neutrality condition η a + η b + η c = 1 allows us to eliminate one of the six subband occupation factors. Regarding the hopping term, since sin θ ≥ 0 for 0 ≤ θ ≤ π, for having E T 0 ≤ 0, the condition is that (η a − η c ) cos φ ≥ 0. As noted above, the total energy in Eq. (4) is invariant under the exchange of the subband labels "a" and "c". And as the angle φ only enters through the hopping term E T 0 , its value is just determined by the sign of η a − η c . For example, by assuming that we restrict ourselves to the configurations with η a ≥ η c , the optimum value for φ is φ = 0. Having assumed instead that η a ≤ η c for all possible configurations, the optimum value for φ will be φ = π. Both choices are equivalent, and we have adopted here the first one: η a ≥ η c , φ = 0. Under these constraints, E 0 has been minimized numerically with respect to the remaining six variational parameters: five occupation factors and the layer mixing angle θ. The corresponding results are presented in Figs. 2-7. The scaling with r s and d * of the different energy contributions deserves some discussion. In the lowdensity limit, r s grows since it is inversely proportional to the square root of the density and the physics is dominated by the tunneling and layer contributions, which are single-particle effects. On the other side, in the high-density limit r s is small, and the main contributions to the total energy comes from the intra-layer kinetic and Hartree terms. For large enough d * , on the other side, the system always becomes dominated by the Hartree term, which favors an effective bilayer configuration with the central well empty, and with electrons distributed equally in the two side layers [14,15]. The intra-layer exchange interaction of Eq. (8) scales linearly with r s , is always negative, and is important in favoring spin-polarized ground-state configurations. Regarding the inter-layer exchange contribution of Eq. (9), its scaling with r s and d * is less trivial, mainly due to the presence of the ratio d * /r s in the arguments of the exponentials. As discussed in the Appendix A, in the limit d * /r s 1, and by expanding the exponentials this term acquires a leading linear dependence in d * , like the Hartree electrostatic contribution. In the opposite limit d * /r s 1, the exponential terms are small and the interlayer exchange scales linearly with r s . According to our numerical evidence, this term is always positive and, in consequence, its optimum configuration is either θ = 0 or the spin-balanced case η aσ = η bσ = η cσ , as in both cases E X-inter 0 = 0.

III. RESULTS
The total energy in the 6-parameter space is plagued by local minima, that makes the task of finding the global minimum a difficult numerical challenge. We also have to minimize fulfilling the constraint 0 ≤ η i ≤ 1, with the minimum being just at the boundary in some cases. To carry out the minimization we have partitioned the 6-parameter space in (typically) 10 6 regions. Starting from a central point for each region, we find a local minimum using a Simplex algorithm [21]. Then we found the global energy minimum as the minimum among all regions. That procedure is easy to parallelize, in particular we have implemented it using MPI (Message Passing Interface) facilities [22].
Before discussing the full numerical results, we find convenient to analyze some particular limits of the trilayer system, which admits either analytical or semianalytical solutions. Most of the ground-state configurations will appear already in these simple limits, helping the understanding of the general cases discussed later.
A. d * = 0 limit In the limit of small distance between layers, the expression for the ground-state energy simplifies greatly, having made the choice ε * 1 = ε * 3 = 0, ε * 2 = 0. We display in Fig. 2 the ground-state phase diagram which results from the numerical minimization of Eq. (10), for the case t * = 0 [23]. The more prominent feature here is the growing of configuration B (η b↑ = η b↓ = 1/2) with respect to all configurations at lower values of r s , as − ε * 2 increases. This is clear physically: as θ = 0, the configuration B is equivalent to η 2↑ = η 2↓ = 1/2, and a negative value of ε * 2 favors the central well filling. The boundary between configurations A and B does not depend on ε * 2 , since the energy of the term proportional to ε * 2 is the same in both configurations (see Eq. (10)). The value r s = 2.011, which is the boundary between A and B is determined then by the balance between the intra-layer kinetic energy and the intra-layer exchange  Table I.   TABLE I. Ground-state configurations for d * = t * = 0. In all cases θ = 0. Note that for φ = θ = 0, the subband occupation factors are the same that the layer occupation factors. The symbol "sp" is used to identify spin-polarized states, while "ml", "bl", and "tl" represent monolayer, bilayer and trilayer configurations respectively.
Configuration η a↑ η a↓ η b↑ η b↓ η c↑ η c↓ A: sp, ml 0 0 1 0 0 0 B: ml 0 0 term. Indeed, for ε * 2 = 0 we can apply the analytical considerations from Ref. [4], according to which the critical density at which a transition occurs from a configuration with p equally occupied components to another configuration with p+1 equally occupied components is given by r   s (5, 6) 1.008. These analytical results are exactly the transition points at ε * 2 = 0 in Fig. 2, obtained numerically. For t * = 0 (the case analyzed in Fig. 2), exists the symmetry η a ↔ η c . This means, for example, that configuration C, corresponding to η b↑ = η b↓ = x, and η a↑ = 1 − 2x, is degenerate with the configurations η a↓ = 1 − 2x, or η c↑ = 1 − 2x,  Table II.   TABLE II. Ground-state configurations for d * = ε * 2 = 0. In all cases θ = π/2. or η c↓ = 1 − 2x. For finite t * , we will see how some of these degeneracies are broken. Configurations A and B are the only cases that we have found in the present work where the trilayer is actually a single-layer from the point of view of the electronic distribution, with all electrons located in the central layer. The fully spinpolarized configuration A is preferred over the unpolarized configuration B at densities such that r s ≥ 2.011 by the action of the intra-subband exchange term in Eq. (10), which always favors spin-polarized configurations. As most of the transitions in this work, at the boundary between A and B configurations, the occupation factors change abruptly.
As a way to understand the effect of the hopping parameter t * , we display the ground-state configurations in the parameter space r s − t * , for d * = ε * 2 = 0 in Fig. 3. From Eq. (10) with ε * 2 = 0, the tunneling energy attains its minimum value with the choice θ = π/2, η a > η c . In the low-density limit r s 1, the tunneling term is the dominant contribution and all the electrons are in the atype subband (either in the α or β configurations). In the high-density limit r s 1, the physics is dominated by the intra-subband kinetic energy term, whose energy is optimized (lower) by spreading electrons among the three subbands a, b, and c. The intra-subband exchange energy contribution, scaling linearly with r s and favoring spin and pseudo-spin polarized states, plays an important role at intermediate densities, by favoring spinpolarized configurations over spin-unpolarized configurations as r s increases (i.e., β → α, δ → γ, ζ → ). For increasing t * , the configuration β becomes more stable and gains area in the parameter space at the expense of the remaining configurations at higher densities. In other words, at constant r s , the occupancy of the bonding a-type subband increases with t * , for all r s < r (0) s (2, 3) = 1.513, until the system falls into β configuration.
At first sight it is not clear why the ground-state configurations for d * = t * = ε * 2 = 0 displayed in Figs. 2 and 3 are not the same. The point here is that in Fig. 2 (3) we display the configurations which are the ground-states for finite values of ε * 2 (t * ). The case where the three external parameters d * , t * , and ε * 2 strictly vanish is somehow ill-defined, as the total energy then reduces to just the first two terms in Eq. (10), that left the angle θ undetermined. Since any value of θ is allowed, the ground-state configuration is not unique. The same consideration applies to the results displayed in Fig. 2 of Ref. [4].
The competition between the tunneling strength parameter t * and the central layer energy ε * 2 is shown in Fig. 4. For 0 < t * 0.2, the ε * 2 -driven configurations of Fig. 2 are the ground-state ones; for larger values of the tunneling parameter, the t * -driven configurations of Fig. 3 are the ones with lower energies. The boundary between configurations A and α may be easily determined from Eq. (10), and found to be t * = − ε * 2 /(2 √ 2). For ε * 2 = − 0.5, this gives t * 0.18, in agreement with the numerically calculated boundary in Fig. 4. The boundary does not depend on r s , due to the equivalent occupation factors in configurations A (one fully occupied spin-polarized layer) and α (one fully occupied spin-polarized subband), and the identical quadratic scaling with r s of both tunneling and ε * s energy contributions. For the same reason, the boundary between the B and β configurations does not change with r s .
Considering configurations C, D, and E, they are similar to the ones found in Fig. 2 regarding the occupancies, but with θ = 0. The minimizing angle θ may be found from Eq. (10), as follows. Calling s = sin θ, and given the simplicity of the s dependence in the last terms, E 0 (d * = 0) may be optimized analytically with respect to s, . (11) This is quite reasonable: for t * = 0, and given that in  Fig. 3. For 0 < t * 0.2, the configurations A, B, and F are the same as in Fig. 2. The three remaining configurations C, D, and E are similar to the corresponding ones in Fig. 2 regarding the subband occupancies, but with θ = 0.
configurations C, D, and E, η a = η c , the system gains some tunneling energy by allowing that θ = 0, although it is found numerically in all cases that the value of the minimizing angle is small.
It is worth of emphasize that at the (A,α) boundary crossing, the trilayer system suffers an abrupt reaccommodation of the electronic charge. This is easily seen from the following set of equations relating the occupation densities in the layer and subband basis: Since in the A configuration η b↑ = 1 and θ = 0, this means in real space that η 2↑ = 1. In the α configuration, since instead η a↑ = 1 and θ = π/2, this translates in real space to the distribution η 1↑ = η 3↑ = 1/4, η 2↑ = 1/2.
In words, at the (A,α) transition the system passes from a monolayer to a trilayer configuration, as tunneling increases. It is interesting to note that a "balanced" configuration in the {1, 2, 3} layer space is also a "balanced" configuration in the {a, b, c} subband space.  Table III. The straight line corresponds to the approximation given in Eq. (13) for the boundary between α and β configurations.
TABLE III. Ground-state configurations for ε * 2 = 0, with the layer separation dependent hopping t * (d * ), and for large values of d * . Configuration P3 is present between P4 and P5, but not visible in the figure because it is very narrow. The symbol "sp*" is used to indicate that this configuration is degenerated with a non-polarized configuration obtained flipping the spin of a component.

B. Numerical results for the full model
In Fig. 5 we display the numerical results for the ground-state configurations, in the r s − d * plane, with the layer distance dependent hopping t * (d * ), as given by Eq. (B2) in Appendix B. For d * 1, the hopping is sizable (t * (0) 1), and the ground-state configurations are the same as in Fig. 3. For d * 1, the Hartree term instead becomes predominant, and the zero-tunneling and zero-site energy ground-state configurations of Ref. [4] are obtained. The quadratic scaling of E T 0 with r s , however, stabilizes the tunnelingdriven configurations as r s increases. In other words, for enough large values of r s , the trilayer system always enter in a tunneling dominated regime. The configuration labeled P2', which is different from the spontaneous inter-layer coherent state P1 configuration found previously (η b↑ = 1, θ = π/2) [4], results from a competition between the hopping and inter-subband exchange contributions. The P2' configuration may be understood as a compromise between the tunneling stabilized β configuration, and the Hartree induced P2 configuration. In passing from β to P2, θ vanishes and the system looses tunneling energy. In P2', on the other side, the system keeps the gain in tunneling energy by allowing θ = 0, but at the same time minimizes to some extent the intersubband exchange energy contribution by inducing a more "balanced" subband space configuration (one subband is occupied in β, two subbands are occupied in P2').
Another interesting feature of the phase diagram in the Fig. 5 is the behavior of the frontier between the α (spin-polarized) and β (non spin-polarized) configurations. For d * = 0, the limit is in the expected value of r s = 2.011 (as in Fig. 3). However, when d * grows and the associated d * -dependent hopping t * (d * ) diminishes, the frontier between polarized and non-polarized configurations is moved to lower densities. Even for the lowest density considered (at r s = 3.0), it is possible to find a situation with the trilayer in a paramagnetic configuration. This is due to the d * dependence of the inter-subband exchange contribution. As explained in Appendix A, in the limit d * /r s → 0, one gets a linear dependence with d * for this term, as shown in Eq. (A5). Considering that the tunneling and Hartree energies are the same in α and β configurations, the boundary between both can be found by imposing the condition that the sum of the intra-subband kinetic energy, and the intra-and inter-subband exchange energies be the same in both cases. In the limit d * /r s → 0, one can use the expansion in Eq. (A5) for the inter-subband exchange energy and obtain the equation for the boundary valid to linear order in d * . The important point here is that the gain in energy of the β configuration comes from the inter-subband exchange interaction, which in turn is a consequence that the tunneling between layers is finite, allowing a subband mixing angle θ = 0. For even higher d * 's, the β configuration turns into P2'. When the layers are more separated, the angle θ flips from π/2 to 0, the system lost tunneling energy and configuration P2 becomes the one with the lowest energy.
The only difference between Fig. 5 and Fig. 6 is that in the former ε * 2 = 0, while in the latter ε * 2 = − 0.8. The main consequences are the disappearance of configurations P2 and P4 of Fig. 5 in Fig. 6, being replaced by configurations P3' and P3", and the growing stability of the tunneling-driven configurations on the lefthand side of the diagram. This last feature is easy to understand: by doing ε * 2 < 0 the occupation of the central layer increases, and this favors the stability of the α and β configurations particularly, that when translated to layer space have a central layer with twice the occupancy of the two side layers. This effect is also reflected in configurations γ, δ, , ζ, although somehow in a less strength, due to the decreasing influence of tunneling and layer energies for decreasing r s .
Regarding the disappearance of the P2 and P4 configurations which are present in Fig. 5 but not in Fig. 6, this is due to the fact that both configurations have zero occupancy of the central layer, and this is in direct conflict with the fact that ε * 2 < 0. The new stable configurations P3' and P3", on the other side, have a finite occupancy of the central layer.
As in the ε * 2 = 0 case, the non spin-polarized configuration β is still present even at low densities (r s > 2.0) and low and medium inter-layer distances. However, when d * increases (and consequently, the hopping decreases) the configuration A is stabilized. This configuration implies that all the electrons are in the central layer with the same spin projection. It is stabilized by the exponential decay of the tunneling parameter with d * , that somehow recreates the situation in Fig. 2, in the limit of low densities. For larger d * , the Hartree energy becomes important and a situation with the three layers occupied, like P3', is preferred.  Finally, we display in Fig. 7 the layer occupancies as function of r s , for fixed d * ( 2), t * ( 0.05) and ε * 2 ( − 0.8) [8]. For r s → 0, the layer occupancies in the figure are well reproduced by the analytical estimates given in Ref. [4], according to which the layer occupation factors in the P6 configuration are given by

Layer occupation factors
with ξ(d * ) = (1 + 3/4d * ) −1 . Evaluating these expressions at d * = 2, we obtain that η a = η c = 5/11 0.45, η b = 1/11 0.09, in good agreement with the numerical values at r s = 0 in Fig. 7. To obtain these analytical estimates, only the intra-layer kinetic and Hartree energies contributions to the total energy are considered, as dictated by the scaling behavior of the total energy in the high-density limit. This abrupt redistribution of charge among the layers is reminiscent of the abrupt charge transfer between the ground-and first-excited subbands of a single semiconductor quantum well, observed experimentally in Ref. [24] and discussed theoretically in Ref. [25].
For increasing r s , the sequence of transitions is as follows: P6 → P5 → P3' → P3" → P3' → A. The transitions P6 → P5, P5 → P3', and P3' → A are clearly discontinuous, as at each one of them a given layer occupation factor abruptly passes from a finite value to zero. The interesting "reentrant" sequence P3' → P3" → P3', on the other side, involves only smooth occupation changes at the boundaries. The full sequence can be understood as a consequence of the progressive filling of the central well as r s increases. The quadratic scaling with r s of the term proportional to ε * 2 in Eq. (5) makes this term to become dominant for large r s , culminating with the stabilization of the A configuration, where all electrons are in the central layer, and fully spin polarized. The "reentrant" sequence P3' → P3" → P3' is interesting, since it reflects once more the importance of the inter-layer exchange energy contribution E X-inter 0 , and its competition with the tunneling energy E T 0 . In the P3' configuration, θ = 0, as this is one of the two possible ways to cancel the positive contribution of E X-inter 0 , although that leads to a vanishing of the gain in tunneling energy too. But for 1.68 r s 2.05, the P3" "balanced" configuration becomes the one with the lowest energy. Here, since is minimized using the second option available for making it as small as possible: equal occupancy of all the subbands. This allows that θ = 0, and leads to a gain of the tunneling energy. As r s increases further, the equal occupancy constraint cannot be maintained any more, and the system returns to the unbalanced P3' ground-state configuration, but now with a preferential occupancy of the central layer.
It should be noted that the three parameters d * , t * , and ε * 2 used in Fig. 7 were directly obtained from the triple semiconductor quantum well system studied experimentally in Ref. [8], as explained in Appendix B. After adjustment of these parameters, our model reproduces qualitatively the main features of the more elaborated calculations reported in Ref. [8], using densityfunctional-theory in the Local Density Approximation. Due to the associated computational cost, these types of calculations are usually restricted to a given particular set of parameters. In particular, it is reported in Ref. [8] that η 1 η 2 η 3 1/3 at N 10.8 × 10 10 cm −2 (r s 1.74, corresponding to the upper empty white circle in Fig. 6), and that η 1 /η 2 η 3 /η 2 2 at N 16.7 × 10 10 cm −2 (r s 1.4, corresponding to the lower empty white circle in Fig. 6) [26]. These LDA determined layer occupation factors are in good qualitative agreement with the ones from our model, as shown in Fig. 7. This gives us some confidence that our simple model is able to reproduce the results of more elaborate calculations, and that it may be useful for the design of real samples and for understanding the physics of trilayer semiconductor systems.
One may wonder how much of the findings of the present work can be extrapolated to other apparently similar 2DEG's, as for example those formed at the interface between two band insulators, such as SrT iO 3 and LaAlO 3 [27]. These are fascinating systems, showing evidence of superconducting and ferromagnetic prop-erties, simultaneously [28]. Unfortunately, our simple theoretical model is not suitable for capturing the main physics needed for describing these 2DEG's at oxide interfaces. In the first place, the experimental results suggest the presence of two type of carriers at the interface [29]: one type is mobile, and presumably responsible of the superconducting features, and the other type of carriers seems to be localized, and responsible of the magnetic (ferromagnetic) features. Our model has only mobile carriers, which are the ones that gives all the possible ground states presented above. Other important difference is that we have assumed from the very outset that our semiconductor-based 2DEG has translational invariance in the x − y plane, an assumption that seems to be not justified in the case of the 2DEG at oxide interfaces [28]. And lastly, all the evidence in this latter case points to the importance of correlation effects in describing their physics [30], beyond the reach of our variational HF approximation, which on the other side is reasonable for the treatment of the weakly-correlated semiconductor systems discussed here.

IV. CONCLUSIONS
The possible ground-state configurations of a trilayer system have been determined, within the framework of a variational Hartree-Fock approximation. The metallic layers are Coulomb coupled through the inter-layer Hartree and exchange interactions, and also due to the tunneling between the neighboring layers. At highdensity the system becomes quite simple, as the only remaining terms in the total energy in this regime are the intra-layer kinetic energy, and the Hartree classical contribution. The low-density regime is dominated by two single-particle effects introduced in this work: the tunneling between layers and the site energies. The interlayer exchange interaction is found to play an important role, helping in the stabilization of the balanced configuration, in whose vicinity most of the experimental samples are designed. It is expected that the results presented in this work for a wide range of parameters, may serve as a qualitative guide for the design of experimental samples, and also be useful for the understanding of the physics of trilayer semiconductor systems at zero magnetic field.
Here, φ loc (z − d) and φ loc (z) are the envelope normalized wavefunctions corresponding to the left and right quantum wells, respectively, and V b (z) is the potential barrier between the two wells. Consistently with our trilayer model in Fig. 1, we have approximated each quantum well by an isolated attractive deltapotential of strength α. Within this model, φ loc (z) = α * /(2a * 0 ) exp (−α * |z|/2a * 0 ), V b (z) = −α δ(z) + δ(z − d) , with α having units of energy times length, and α * = α/(a * 0 Ry * ). Replacing in Eq. (B1), and imposing the constraint that t * (0) 1 (a reasonable physical choice) we obtain The free parameter α * can be fixed now by imposing the second condition that for a given value of d * , the hopping parameter should be equal to some convenient value, usually taken from a more elaborate calculation. Solving equation above for α * , it yields α * = − 2 ln(t * ) d * . (B3) As an example of the use of these equations, using Eq. (3.6) in Ref. [3], the value of t * may be estimated from the electronic subband structure of a particular trilayer. For instance, using the Local-Density-Approximation theoretical results in Ref. [8] corresponding to a trilayer with d * 2, we obtain that t * 0.05. Replacing in Eq. (B3), α * 2.97. This gives us the dependence of the hopping parameter with the distance between layers employed in Fig. (6). The parameter ε * 2 was estimated by following the considerations of Hanna and MacDonald in Ref. [3], for the same trilayer system.