Fabien Trousselet

Pamela Rueda-Fonseca

Arnaud Ralko

Competing supersolids of Bose-Bose mixtures in a triangular lattice

We study the ground state properties of a frustrated two-species mixture of hard-core bosons on a triangular lattice, as a function of tunable amplitudes for tunnelling and interactions. By combining three different methods, a self-consistent cluster mean-field, exact diagonalizations and effective theories, we unravel a very rich and complex phase diagram. More specifically, we discuss the existence of three original mixture supersolids: (i) a commensurate with frozen densities and supersolidity in spin degrees of freedom, in a regime of strong interspecies interactions; and (ii) when this interaction is weaker, two mutually competing incommensurate supersolids. Finally, we show how these phases can be stabilized by a quantum fluctuation enhancement of peculiar insulating parent states.

I. INTRODUCTION

In nowadays condensed matter physics, common fascinating collective behaviors and novel quantum phases are reported in various domains of physics, more specifically in bosonic systems encountered in quantum magnetism, ultracold atoms on optical lattices and strongly correlated materials. Thanks to their versatility, an important amount of exotic phases has been reported in the literature these last years, both experimental and theoretical, such as different types of superfluids 1 , insulators 2,3 , Bose metals [4][5][6] or supersolid phases [7][8][9][10][11][12] . In the latter, the system enters a phase combining crystalline order and superfluidity, and typically arising from the quantum melting of a Mott insulator; this phenomenon is at the origin of intense scientific activities and debates; Experimentally, indications for supersolidity were supposed to be found in 4 He 8 . However, necessary conditions for continuous symmetry breaking questioned this interpretation 9 and recent experiments have clearly ruled out this scenario 14 .

Meanwhile, supersolids are easier achieved on lattice systems; whereas on a square geometry with nearest neighbor interactions a soft-core description is required 15 , they can also be found in hard-core bosonic models when frustration is induced by either further neighbor interactions 10 and/or lattice geometry, e.g. triangular [11][12][13] . An other interesting direction to stabilize supersolidity is to increase the number of degrees of freedom in frustrated systems. In that respect, bosonic mixtures with several species of bosons, usually encountered in optical lattices 16,17 , either heteronuclear 18 or homonuclear 19 , as well as bilayer systems with interactions but no hopping between layers 20 are very promising. Such systems allow for an even broader variety of quantum phases than their single-species counterparts, depending on the intra-and inter-species interactions, and on the dimension and lattice connectivity.

Theoretically, compared to an already rich literature on mixtures in 1D 21 only few works have focused on 2D cases for instance on square 22 and triangular [START_REF] Mathey | Fermi mixtures have been considered on the triangular lattice[END_REF][24][25] lattices. In all those systems, when contact and dipolar interactions 20,26,27 are taken independently, interactioninduced insulators, e.g. density-homogeneous or density wave 28,29 are favored. Hence, their competition can, along with quantum fluctuations, trigger various unconventional phases 11,20,27,30 , especially in systems with kinetic or interaction frustration.

FIG. 1: Triangular lattice hosting bosons of two species labelled a (blue) and b (brown). the bosons can hope from two neighboring sites with the energy -t, and interact via an intra-species interaction V and a point contact inter-species interaction U .

In this paper, we provide a theoretical study of a twospecies bosonic mixture on a triangular lattice as depicted in Fig. 1. We aim to focus on the most simple model with competing interactions in frustrated geometry in order to study possible mechanisms for stabilizing exotic phases. The richness of the phase diagram of a related one-species model with hard core constraint 12 allows us to expect even more exotic physics for a mixture of two mutually interacting species. Even though such interactions can be encountered in cold atom systems 11,12,20 with polarized dipoles 26,31,32 , the present work describes the maximally frustrating cases for which all interactions are repulsive. We address the following issues: (i) how on-site and nearest neighbor interactions compete and (ii) how lattice frustration impacts on the stabilization of non-conventional phases of such a mixture.

To achieve these goals, we use mainly a cluster mean field theory (CMFT) and exact diagonalizations (ED) on periodic clusters; both methods are detailed, along with the model we consider, in Section II. Note that preliminary Quantum Monte-Carlo (QMC) simulations have been performed to support our findings (see text). We present in Section III an overview of the phase diagram obtained this way, before focusing on the most interesting phases. First, we describe in Section IV some commensurate spin-like phases. These are characterized with means of perturbative approaches; they have frozen densities and, for one of them, spin-like supersolidity. Next, we analyze in Section V incommensurate phases found in this study: these include two original two-species supersolid (SS) phases, which belong to our main findings. Eventually we address in Section VI the nature of phase transitions involving these peculiar phases, before some concluding remarks in Section VII.

II. MODEL AND METHOD

We study a two-species (spin) extended Bose-Hubbard model on a triangular lattice:

H α = i,j -t α (b † iα b jα + h.c.) + V α n iα n jα -µ α i n iα H ab = α H α + U i n ia n ib . (1)
H α is the one-species Hamiltonian (α = a, b), with t α and V α respectively the nearest neighbor hopping and interaction amplitudes; µ α is the chemical potential for bosons of species α created by operators b † iα at site i. In this work, we focus on the limit of hard core bosons (U α,α ≫ |t α |, V α , µ α ) with an implicit onsite intra-species repulsion U αα i n iα (n iα -1)/2. Together with the repulsive interspecies coupling U , this allows us to maximize the effects of frustration. We have studied the more general case in function of independent µ a and µ b and found that the richest physics was found in the symmetric case µ a = µ b , including the novel supersolid regimes. Hence, when the system is (a, b)-symmetric (t α = t, V α = V and µ α = µ), as considered in this work since the most original two-spin phases arise there, the particlehole transformation b ′ iα = b † iα allows for a mapping between µ * > 0 and µ * < 0, with µ * = µ -3V -U/2 a rescaled chemical potential. Note that, such a transformation, if restricted to a single species and for µ * = 0 amounts to changing the sign of U . This U < 0 case is of interest since the corresponding terms mimic qualitatively the effects of interlayer short-range interactions in dipolar cold atom bilayers 20 .

We compute the ground states of H ab on periodic clusters with up to N = 12 sites (see Fig. 2), using either CMFT (if not precised) or ED methods. Note that the former method has been employed successfully in many one-species bosonic systems on various lattices such as triangular 12,28,35 , pentagonal 33 and hexagonal 2 . It is thus expected to be also very efficient in the present model.

N 3 , Λ 1 3 N 6 , Λ 1 2 N 12, Λ 2 3

FIG. 2:

Clusters considered in our CMFT analysis with N = 3, 6 and 12 sites. Internal and external bonds correspond respectively to continuous and dashed lines. For each cluster, the cluster scaling parameter λ is given.

In the former case, the N i (N e ) internal (external) bonds are treated exactly (at the mean field level) 12,34,35 and correlations are better taken into account as N increases. As in Ref. 35 , we thus define the scaling parameter λ = N i /(3N ) which quantifies finite boundary effects; the Thermodynamic Limit (TL) is achieved for λ → 1 (infinite lattice). The mean-field parameters determined self-consistently are the densities niα = n iα and the superfluid fractions (SF) φ iα = b iα . In addition, to evidence 3-fold symmetry breaking, we define the order parameter

M α = |n α (k)|, the diagonal S d ab = | (n a -n b )(-k)(n a -n b )(k) | and off-diagonal S od ab = | i,j b ia b † ib b jb b † ja |/N 2 correlation
functions. We used the Fourier transform n α (k) = 1 N s n sα e ik•rs at point k = (4π/3, 0), corner of the Brillouin zone.

III. OVERVIEW OF THE PHASE DIAGRAM

Let us first focus on the (µ * , U ) phase diagram in which two-species phases -detailed along the paper -emerge, illustrated here for t/V = 0.15. For these parameters, the one-species Hamiltonian H α is known to present a rich phase diagram with either empty/full, homogeneous superfluid, √ 3× √ 3 solid at density 1/3 or 2/3, or supersolid phases 11 . When U is switched on, the correlations between the two species increase, and the resulting GS can be either a product of two one-species phases or a speciesentangled state as depicted in the phase diagram shown in density, various supersolids with 3-sublattice structures are stabilized; they will be discussed in detail in Section V. Finally, an insulating regime with n = 1 is stabilized for U > 1.2V in the vicinity of µ * = 0; it will be the object of Section IV. For µ * > 0, we obtain an equivalent phase diagram thanks to the particle-hole symmetry.

IV. COMMENSURATE SPIN-LIKE PHASES

For strong inter-species repulsion U ≫ V, t, |µ * | (triangular uppermost domain on Fig. 3-up) the system is Mott insulating, with φ α = 0 and n = 1. Indeed, U imposes the local constraint of single occupancy defined as nia + nib = 1. This is reflected by n = 1 plateaus in both ED (Fig. 4-a) and CMFT (Figs. 4-b and7) results. In order to better describe this regime, we define spin 1/2 operators σ z i = (n ia -n ib )/2 and σ + i = b † ia b ib . At second order of the perturbation theory, we obtain an effective XXZ model

H XXZ = - J ⊥ 2 i,j (σ + i σ - j + σ - i σ + j ) + J z i,j σ z i σ z j , (2)
where J ⊥ = 4t 2 /U and J z = 2V + 4t 2 /U . Interestingly, this model predicts a (2m z , -m z , -m z ) SS with diagonal σ z = 0 and off-diagonal σ + = 0 order parameters when J z /J ⊥ > 4.6(1); otherwise, a SF with only off-diagonal order 36,37 . In the present context, these phases correspond respectively to a spin supersolid (SSS) and a spin superfluid (SSF). While the density is uniform, the SSS has a (2m z , -m z , -m z ) structure where

V. INCOMMENSURATE SUPERSOLIDS

In contrast to the above mentioned SSS in which supersolidity comes from the species (spin-like) degrees of freedom, in SS phases discussed here n α may be incommensurate and have finite φ α . From the species-resolved observables shown in Fig. 6 and7, we identify three such SS phases. (i) For large U/V = 2 (right column) and µ * /V about ±1.0(2), the finite M α indicates a 3-fold order (3FO for both species, while only one species has a non-zero φ α . This phase, dubbed 1SF/2CO in Fig. 3, is a rearrangement of two one-species phases of H α , the √ 3 × √ 3 solid and the SS obtained from it by enhancing the quantum fluctuations 11,12 . Two sublattices are (almost) filled by a and b bosons respectively, while on the remaining, an incommensurate density for bosons of one species (e.g. a) accounts for superfluidity (φ a > 0) and a population imbalance (n a > n b ). Within this structure,

n Α U V 2 M Α Φ Α Μ V BPB 3SS 2 3
FIG. 6: Species-resolved observables (nα, Mα, φα). Each species corresponds to either continuous or dashed line. All data come from a N = 6 CMFT with fixed t/V = 0.15 and U = 0.5V The arrow points out a tiny region which disappears in the TL (see text). The shaded regions correspond to the different phases of Fig. 3.

the repulsion energy ∝ U is minimized thanks to the localization of bosons of a single species. (ii) A distinct phase is found for small U/V = 0.5 (left column) and µ * /V about ±0.5 (3), with a finite φ α for both species as well as a 3-fold symmetry breaking.

As shown in the typical snapshot obtained by CMFT in Fig. 8(c), it can be seen as a superposition of two onespecies SS (these would be obtained at U = 0); both species contribute symmetrically to superfluidity. This phase, called 3-fold SS (3SS) is the first example of a collective two-species supersolid and is obtained from the parent n = 2/3 solid by a defect condensation upon doping as µ * increases. As in this parent state, the weak inter-species coupling merely forces the localized a and b bosons to occupy distinct sublattices. (iii) The most remarkable incommensurate SS phase is achieved, upon increasing t/U , when quantum fluctuations become too strong for the density-uniform n = 1 SSS phase. This original two-species SS is called the bosonic pinball (BPB) due to a structure very similar to its fermionic coun-terpart with similar interactions, the pinball liquid, and is depicted in Fig. 8(b,d). The latter has almost one localized electron per site on one sublattice (pins) and a metallic behavior on the remaining hexagonal lattice (balls) 28,29,38 . Here, the particles are bosonic and the species play the role of the spins. The structure is depicted on Fig. 1(b). One sublattice, forming a triangular super-lattice with nia + nib close to 1, is filled in majority by one type of bosons. The two-spin superfluid character is carried by the remaining bosons on the complementary hexagonal lattice. This is shown in the BPB region on Fig. 6 and Fig. 9, by a coexistence of solid (M α = 0) and superfluid (φ α = 0) orders in both species. This lattice symmetry breaking is reminiscent of the parent SSS phase, the partial quantum melting of which involves a condensation of two types of defects coming from doublyoccupied and empty sites. Finally, unlike the 3SS and 1SF/2CO, the BPB can be stabilized at the n = 1 commensurability when µ * = 0.

VI. THERMODYNAMIC LIMITS AND PHASE TRANSITIONS

All the phases described above exist in the TL. Indeed, in Fig. 9 affected by size effects thus Fig. 3 is representative of the TL. In the strong correlation limit U, V ≫ t [39], perturbation theory allows to locate phase transitions where the defect energy vanishes (condensation), e.g. for the 2/3 → 3SS transition. Consider a defect consisting of an extra boson (say of species a) inserted in the empty sublattice of the 2/3 crystal with the energy cost 3V . Via second order processes, such a defect can hop to second neighbors with an -t eff = -t 2 /U -t 2 /V . The terms -t 2 /U and -t 2 /V account respectively for processes where (i) the extra a boson hops via the b-filled sublattice, or (ii) a vacancy on the a-filled sublattice is created and then deleted. In this limit, the defects condense and lead to supersolidity for µ c = 3V -6t eff , e.g. µ c /V ≃ -0.65 for U = 0.5V in good agreement with the results of Fig. 6 and7. This defect condensation mechanism implies a gauge symmetry breaking, additionally to the lattice symmetries already broken in the 2/3 phase. This indicates that the transition is of second order, as confirmed by the absence of discontinuities in φ and M as function of µ (Fig. 9). This corresponds to the standard picture of defect condensation accounting for an incommensurate SS; the SF density (here inhomogeneous) is carried by defects on top of density modulations 11 . In contrast, the transitions between (i) homogeneous SF and 2/3 solid, and (ii) the BPB and 3SS are found of first order, characterized by hysteresis in the CMFT. We find distinct transition points if µ is, from one CMF calculation to another, either stepwise increased or stepwise decreased (see [40]). In case (ii), n a -n b (Fig. 9(b)) has distinct size-dependence, but a vanishing value in the TL cannot be ruled out for the BPB, despite a non-linear behavior. n a -n b is strictly 0 for the 3SS. For the BPB, the non-linear behavior makes difficult to extract n a -n b as λ → 1, and various fits can give either zero (linear) or finite (e.g. weighted a + λ b ) value. However, BPB and 3SS have different symmetries as obtained by CMFT (Fig. 8(b,c)) and confirmed by preliminary QMC calculations. Indeed, to check this point, we have performed Stochastic Series Exchange (SSE) QMC simulations on clusters up to 36 × 36 sites for various U in function of µ * , for which real-space correlations n sa n iα (for site i far away from the reference site s to get rid of shortdistance effects) confirm the existence and symmetries of both the 3SS and the BPB (Fig. 8(d,e)). The complete analysis of this model by SSE-QMC, being beyond the scope of this paper, will constitute a separate work.

VII. CONCLUSION

We study an interacting two-species bosonic mixture on a triangular lattice by combining CMFT, ED and perturbative methods, supported by SSE-QMC results. We focus on a region of parameter space in which peculiar phases arise due to the competition between frustration and quantum fluctuations. Within a very rich and complex phase diagram, we evidence three original mixture supersolids, a commensurate spin supersolid (SSS) and two mutually competing incommensurate phases arising from the partial quantum melting of parent states. The most interesting phase, dubbed bosonic pinball (BPB) due to an inner structure very reminiscent of its fermionic counterpart 29 , results from strong inter-species effects. It is worth mentioning that this rich physics is found for attractive U (for specific parameters), a situation more directly connected to dipolar cold atom bilayer experiments. We hope this work will stimulate further investigations in this direction.

Fig. 3 .

 3 (i) As µ * < 0 increases, the first non-trivial phase encountered is a homogeneous two-spin superfluid with φ a = φ b = 0 (φ α = 1 N i φ iα ). When n a and n b are large enough, V and U terms drive the system into a 3fold ordered insulator characterized by φ a = φ b = 0 and a total density n = 2/3 (n a = n b = 1/3), the two-species counterpart of the √ 3 × √ 3 phase. (ii) At even larger

FIG. 3 :

 3 FIG.3: Phase diagrams as obtained from CMFT calculations on a N = 3 cluster at t = 0.15 and V = 1. Top: wide range of µ * showing two-spin superfluids, 2/3 and 4/3 plateaus, XXZ physics domain and original collective two-spin physics (white rectangle). Bottom: zoom of the white rectangle and all new two-spin phases (see text for acronyms). Solid (dashed) lines are phase boundaries subsisting (vanishing) in the Thermodynamic Limit (TL).

FIG. 4 :

 4 FIG. 4: Strong U regime characterized using either N = 12 (ED: a) or N = 6 (CMFT: b) clusters. (a) n, S d ab and S od ab as a function of µ * /V for U = 2V and t = 0.15V . The n = 1 plateau and both finite S d ab and S od ab evidence the SSS; (b) n and φ = φa + φ b as a function of µ * /U distinguishing the SSF (V = 0) and the SSS (V = 0.2U ) regimes.

FIG. 5 :

 5 FIG. 5: Strong U regime characterized using either N = 6 (CMFT: a) or N = 12 (ED: b) clusters. (a) Species-averaged order parameter M and superfluid fraction Φ, evidencing the quantum fluctuation enhancement from the SSS to the BPB phase at V = 1 and t = 0.15; (b) comparison of the kinetic energy gain ǫ * ≃ -kt 2 /U for H ab (filled) and HXXZ (empty) for three (t, V ) sets.

FIG. 7 :

 7 FIG.7: Same as Fig.6, but for U = 2V .

FIG. 8 :

 8 FIG.8: Examples of two-spin supersolids with 3-sublattice structure obtained by CMFT density maps on the 12-site cluster (a,b,c) and by QMC correlations for 36 × 36 sites cluster (d,e). (a) the n = 1 spin supersolid (SSS) at µ * = 0 and U = 2V ; (b) and (c) respectively the bosonic pinball (BPB) at µ * = -0.15V and the 3-supersolid (3SS) at µ * = -0.25V , for t = 0.15V and U = 0.5V . na -n b is 0 for the 3SS and finite for the BPB for all CMFT clusters (Fig.9for the TL study). (d) and (e) are respectively QMC results for the BPB at µ * = -0.12V and the 3SS at -0.72V . Note that the 3SS and BPB have distinct symmetries (under π/6 rotations).

FIG. 9 :

 9 FIG. 9: M and φ as a function of µ * for t = 0.15V and U = 0.5V obtained by CMFT on clusters of N = 3 (continuous), 6 (dashed) and 12 (symbols) site. See Fig. 2 for the cluster shapes. Cluster scalings as function of λ (see text) are shown for b) na -n b with different possible fits giving either zero (continuous) or finite (dashed) TL values (see text) and c) M (filled symbols) and φ/2 (empty symbols) within both the BPB at |µ * | = t (squares) and the 3SS at |µ * | = 3t (circles).

Acknowledgments

A.R. and F.T. acknowledge financial support by the Agence Nationale de la Recherche under grant No. ANR 2010 BLANC 0406-0. P.R.-F. would like to thank Néel Institute for kind hospitality.


Thomas Grange

Electron transport in quantum wire superlattices

Electronic transport is theoretically investigated in laterally confined semiconductor superlattices using the formalism of non-equilibrium Green's functions. Velocity-field characteristics are calculated for nanowire superlattices of varying diameters, from the quantum dot superlattice regime to the quantum well superlattice regime. Scattering processes due to electron-phonon couplings, phonon anharmonicity, charged impurities, surface and interface roughness and alloy disorder are included on a microscopic basis. Elastic scattering mechanisms are treated in a partial coherent way beyond the self-consistent Born approximation. The nature of transport along the superlattice is shown to depend dramatically on the lateral dimensionality. In the quantum wire regime, the electron velocity-field characteristics are predicted to deviate strongly from the standard Esaki-Tsu form. The standard peak of negative differential velocity is shifted to lower electric fields, while additional current peaks appear due to integer and fractional resonances with optical phonons.

I. INTRODUCTION

Electron transport in superlattices (SLs) has been widely studied in the case of 1D periodic arrangement of 2D semiconductor layers 1-3 (Fig. 1a). In contrast, the nature of charge transport in superlattices of lower dimensionality remains an open question. The nature and efficiency of dissipative processes are well known to strongly depend on the dimensionality, with longest lifetimes and coherence times in 3D quantum confined structures 4-7 . Quantum dot (QD) SLs are hence expected to display radically different non-equilibrium transport properties than quantum well (QW) SLs. Potential applications for high performance thermoelectric converters 8 , photovoltaics 9 , and quantum cascade lasers [START_REF] Wingreen | Quantum Electronics[END_REF]11 further motivates the fundamental understanding of the nonequilibrium transport properties of SLs with 3D quantum confinement.

Though ballistic transport in a QW SL is ideally a pure 1D problem, scattering processes couples all spatial degrees of freedom (Fig. 1a). In QW SLs, the in-plane free electron motion possesses a continuous dispersion. During the scattering processes, energy transfer occurs between the axial and lateral motions. Hence the effective energy conservation laws in the motion along the SL axis can greatly differ from the 3D one of the involved scattering mechanisms. For example, in QW heterostructures, 3D elastic scattering processes due to static defects act effectively as 1D inelastic scattering mechanisms with respect to the motion along the SL. In contrast, in the limit of purely 1D SL with 0D lateral motion (Fig. 1d), the lateral state is frozen so that the 3D energy conservation laws of the scattering processes apply directly to the 1D motion along the SL.

The effect of dissipation on 1D quantum electron transport has been investigated in finite quantum region such as atomic wires 12 and stacked QDs 13 . Quantum dissipative transport in QD SLs has been studied in Ref. 14 for high electric fields. However full velocity-field characteristics of 1D SLs have not been calculated so far. In particular it is still unkown how the Esaki-Tsu negative differential velocity (NDV) 1 evolves when the lateral dimensionality is reduced towards QD SLs.

In this article, we present a detailed theory of electron transport in semiconductor SLs laterally confined in nanowires using the framework of non-equilibrium Green's functions (NEGF). The formalism presented here allows to describe the transport in SLs for any lateral dimensionality from the regime of coupled 0D QDs (Fig. 1d) to the one of coupled 2D QWs (Fig. 1a). For this purpose, we include in our modeling the scattering mechanisms that are known to be important in planar QW heterostructures as well as the ones specific to quasi-0D systems. Specifically, we consider the mechanism of anharmonic polaron decay which has been shown to be dominant in InGaAs self-assembled QDs 6,15 . In addition, as the nature of elastic scattering mechanisms evolves from incoherent to coherent as the dimensionality decreases, we include a coherent treatment of elastic scattering processes for lateral state-conserving processes. The electronic transport properties are shown to depend dramatically on the SL lateral dimensionality. With increasing lateral confinement, the drift velocityvoltage characteristics evolve from the Esaki-Tsu form for usual QW SL, dominated by a single NDV effect, towards richer characteristics where (i) the linear conductivity increases, (ii) the standard Esaki-Tsu NDV occurs at much lower voltages, and (iii) large peaks due to integer and fractional resonances with optical phonons appear.

The paper is organized as follows.

The theory of electron transport in laterally confined SLs is described in section II. The application to the transport in GaAs/AlGaAs nanowire superlattices of varying diameters is presented and discussed in section III.

II. THEORY

Our aim is to model the non-equilibrium electron transport in a nanowire SL heterostructure. We first introduce the general form of the Hamiltonian. We then introduce an electronic basis set, and present the NEGF formalism. We then treat the relevant scattering terms arising from electron-phonon interactions and static disorder effects. Finally we present the method used for solving the self-consistent NEGF equations.

A. Hamiltonian of the nanowire superlattice

We consider charge carriers in a nanowire SL in interaction with phonon modes. The following Hamiltonian is used to model the system:

Ĥ = Ĥ0 + Ve + Ve-ph + Ĥvib + V m-f e-e , (1)
where Ĥ0 is the electronic Hamiltonian of a single charge carrier in the idealized heterostructure under a homogeneous electric field (the length gauge will be used throughout this paper); Ve is the static electronic potential due to disorder effects such as interface and surface roughness, charged impurities and alloy disorder; Ve-ph represents the interaction terms between electrons and phonons; Ĥvib is the vibrational Hamiltonian of the crystal, including its anharmonic part; V m-f e-e represents the mean-field potential arising from the Coulomb interaction with other charge carriers. The standard caret notation is used to refer to operators.

B. Electronic structure and basis states

We consider cylindrical nanowires of diameter D = 2R having SL heterostructures along their longitudinal axis. We use cylindrical coordinates (z, ρ, θ) where z is the nanowire axis (see Fig. 1b). The material composition and modeling parameters are assumed to depend only on the z coordinate. The electronic structure is modeled considering a single band within the envelope function approximation. We consider a basis set of the form:

Ψ α,n (z, ρ, θ) = ζ α (z)φ n (ρ, θ), (2)
where the ζ i (z) are defined below and φ n (ρ, θ) are eigenstates of the Schrödinger equation on a homogeneous disk (i.e. the assumed nanowire cross-section):

φ n (ρ, θ) = 1 √ πRJ (mn+1)) (χ ln mn ) J mn χ ln mn ρ R e imnθ (3)
where J m is the Bessel function of order m and χ l m is its lth zero. The indexes m n and l n are chosen such that they are associated with the nth eigenvalue of the 2D lateral motion at the z coordinate:

E n (z) = 2 χ ln mn 2 2m * (z)R 2 (4)
where m * (z) is the material-dependent effective mass.

For each lateral mode we are then left with a z-dependent electronic Hamiltonian:

n| Ĥ0 |p = δ n,p ĥ0 -eF ẑ + E n (z) , ( 5
)
where F is the electric field applied along the z axis and

ĥ0 = -2 2 ∂ ∂z 1 m * (z) ∂ ∂z + V b (z), (6)
with V b (z) being the material-dependent band-offset potential. The eigenstates of the periodic ĥ0 Hamiltonian satisfy the Bloch theorem and can be classified by wavevector κ z ∈ [-π/l p , π/l p ] and miniband index ν

ĥ0 ϕ ν κz (z) = ε ν κz ϕ ν κz (z) (7)
In order to discretize the basis set in the z direction, various choices are possible. A full real-space discretization 16,17 by e.g. finite differences requires a large basis size in order to quantitatively describe the miniband structure, rendering the NEGF numerical implementation challenging. On the other hand, modespace approaches are also possible, and allow the description of the transport properties with a much smaller basis 14,18,19 . Here, the Hilbert space is reduced by keeping only a finite number of minibands, i.e. the low lying states of the unbiased structure. An energy cut-off E c is fixed and only the minibands separated from the ground state by less than E c are kept. Formally we define the following projection operator Pz =

ε ν 0 -ε 0 0 <Ec ν,κz |ϕ ν κz ϕ ν κz | (8)
In the following, the non-equilibrium dynamics will be computed within the subspace P generated by Pz . The various operators Ô defined on the full Hilbert space need to be transformed in operators Ôr acting within the P subspace. The restrictions of the ĥ0 unperturbed Hamiltonian and ẑ position operator are simply defined by: ĥ0r = ĥ0 Pz (9) ẑr = Pz ẑ Pz (10) Within the P subspace, we consider the localized basis which is constructed by diagonalizing the ẑr position operator 14,19 :

ẑr |ζ α = ζ α |ζ α (11)
where the |ζ α (α = 1, 2, ...) are sorted by increasing ζ α eigenvalues. In addition, their phase is chosen in order to fullfill periodicity:

ζ α+pM (z) = ζ α (z -pL) ( 12
)
where p is an integer, M is the number of minibands and L is the period length.

All the terms V(j) in the Hamiltonian other than ĥ0 are functions of the sole ẑ operator. In order to conserve the commutation relation [ V(j) , ẑ] = [ V(j)r , ẑr ] = 0, i.e. to conserve their local nature, we use the following transformation:

V(j)r = α |ζ α ζ α | V(j) |ζ α ζ α | (13)
This transformation allows us to keep exactly the same form for the current operator 34 :

ĵzr = i [ Ĥr , ẑr ] = i [ ĥ0r , ẑr ] (14)
We also keep only a finite number of lateral states satisfying E i -E 0 < E c . In the following NEGF simulations, it will be checked that the E c cut-off value is chosen high enough so that adding higher minibands and/or higher lateral states does not change significantly the calculated current. In the following, we will use exclusively operators acting in the P subspace and all the r subscripts indicating restricted operators will be omitted for simplicity.

C. Non-equilibrium Green's functions

Electron Green's functions

The NEGF formalism is generally introduced in terms of four Green's functions (GFs) [START_REF] Datta | Electronic transport in mesoscopic systems[END_REF][START_REF] Haug | Quantum kinetics in transport and optics of semiconductors[END_REF][START_REF] Caroli | [END_REF][23] . Here the lesser, upper, retarded and advanced electron GFs are defined respectively in terms of the following quantities:

G < α,β,n (t 1 , t 2 ) = i ĉ+ β,n (t 2 )ĉ α,n (t 1 ) (15a) G > α,β,n (t 1 , t 2 ) = -i ĉα,n (t 1 )ĉ + β,n (t 2 ) (15b) G R α,β,n (t 1 , t 2 ) = -iΘ(t 1 -t 2 ) {ĉ α,n (t 1 ), ĉ+ β,n (t 2 )} (15c) G A α,β,n (t 1 , t 2 ) = iΘ(t 2 -t 1 ) {ĉ α,n (t 1 ), ĉ+ β,n (t 2 )} (15d)
where ĉ+ α,n and ĉα,n are respectively the creation and annihilation operators in the electronic state |Ψ α,n . It is convenient to introduce the GF operators Ĝ(t 1 , t 2 ) defined here by:

α, n| Ĝ(t 1 , t 2 )|β, n = δ n,n G α,β,n (t 1 , t 2 ), (16)
where G refers to any of the four GFs. Since we do not consider spin-dependent effects, the spin index is not indicated here for simplicity and all quantities are implicitly spin diagonal. The spectral function reads in operator form:

 = i( ĜR -ĜA ) = i( Ĝ> -Ĝ< ) (17)
The retarded GF can be expressed from the spectral function:

ĜR (t 1 , t 2 ) = -iΘ(t 1 -t 2 ) Â(t 1 , t 2 ) ( 18
)
The lesser GF (15a) is an extension of the concept of density matrix for two different times. Its value at equal times corresponds to the density matrix:

ρ(t) = -i Ĝ< (t, t). (19)
Here, off-diagonal terms of the GFs are considered in the z direction but not in the lateral directions where only the diagonal terms of the |φ n eigenmodes are considered. The calculations are thus independent of the basis choice along the z-axis, allowing full coherent description of transport along the nanowire. In contrast, the use of only diagonal GF terms in the lateral directions does not allow us to describe coherent transport effects in the lateral directions. This seems a reasonable approximation for thin nanowires where the lateral scattering potentials are weak compared to the lateral quantization energies. On the other hand, in the limit of large diameter nanowire and planar heterostructures, this is equivalent to considering only GF terms which are diagonal in the in-plane momentum, similarly to previous implementations of NEGF in planar quantum well heterostructures [16][17][18] . In this later case, lateral localization effects might become relevant at low temperature and/or in presence of strong disorder effects, which remains beyond the scope of this work.

In steady state, the GFs depend only on the time difference (t 2 -t 1 ). In this case we will use one-time-dependent GFs Ĝ(t 1 -t 2 ) = Ĝ(t 1 , t 2 ). Energy-dependent GFs are then defined by Fourier-transforming these quantities Ĝ

(E) = 1 dt Ĝ(t)e iEt/ (20)
Moreover, in steady state, the advanced and retarded GFs are linked by

ĜA (E) = ĜR (E) + , (21)
so that the GFs are forming only two independent quantities.

Phonon Green's functions

The Hamiltonian of the vibration of the crystal reads

Ĥvib = Ĥph + Va (22) Ĥph = ω q,p â+ q,p âq,p (23)
where Ĥph and Va are respectively the harmonic and anharmonic parts, and â+ q,p and âq,p are the creation and annihilation operators for a phonon mode of branch p with wavevector q. The phonon GFs are defined by 24 :

D < q,p (t 1 , t 2 ) = -i Â+ q,p (t 2 ) Âq,p (t 1 ) (24)
D > q,p (t 1 , t 2 ) = -i Âq,p (t 1 ) Â+ q,p (t 2 ) ( 25
) D R q,p (t 1 , t 2 ) = -iΘ(t 1 -t 2 ) [ Âq (t 1 ), Â+ q,p (t 2 )] , (26)
with Âq,p = âq,p + â+ -q,p .

In steady-state, we will use the same one-time notation D(t) and energy-dependent GFs D(E) as for electronic GFs (Eq. 20). When only the harmonic part of the Hamiltonian is taken into account, the lesser and retarded equilibrium phonon GFs for a mode of frequency ω read respectively:

D <(0) ω (t) = -i[(N ω + 1)e iωt + N ω e -iωt ] (28) D R(0) ω (t) = iΘ(t)[e iωt -e -iωt ], (29)
where N ω is the Bose factor at the energy ω. In the energy domain, it reads:

D <(0) ω (E) = -2πi[(N ω + 1)δ(E + ω) + N ω δ(E -ω)] (30)
D R(0) ω (E) = 1 E + i0 + -ω - 1 E + i0 + + ω . ( 31
)
Harmonic phonon GFs are usually directly used in electron transport calculations. However, in this work, we will use instead anharmonic GFs for optical modes which will be derived below (Sec. II F). This will allow us to account for the mechanism of anharmonic polaron decay.

D. Equations of motion

The equations of motion in the NEGF formalism can be expressed in terms of the so-called Dyson and Keldysh relations; in steady-state, they read respectively [START_REF] Datta | Electronic transport in mesoscopic systems[END_REF][START_REF] Haug | Quantum kinetics in transport and optics of semiconductors[END_REF][START_REF] Caroli | [END_REF][23]

: ĜR (E) = E Î -Ĥ0 -Vc -ΣR (E) -1 (32) Ĝ< (E) = ĜR (E) Σ< (E) ĜA (E) ( 33
)
where Vc is the part of the interaction which is treated coherently. Here it reads Vc = V m-f e-e + V rand e

, where V m-f e-e is the mean-field Coulomb potential and V rand e is the part of the electronic disordered potential which is generated randomly and treated coherently (see below). Similarly to the spectral GF, the spectral self-energy is defined as

Γ = i( ΣR -ΣA ) = i( Σ< -Σ> ), (34)
and the retarded self-energy can be expressed as:

ΣR (t) = -iΘ(t) Γ(t) = Θ(t)( Σ< (t) -Σ> (t)). ( 35
)
In the energy domain, this relation reads:

ΣR (E) = - i 2 Γ(E) + 1 2 H[ Γ](E) (36)
where H denotes the Hilbert tranform

H[ Γ](E) = P dE Γ(E )/π(E -E ).

E. Self-energies due to electron-phonon interaction

We consider the interactions of electrons with longitudinal optical (LO), surface optical (SO) and longitudinal acoustic (LA) phonon modes laterally confined in the nanowire. Each interaction is of the form:

Ve-ph = q,p fq,p Âq,p .

(37)

The form factors fq,p are given in appendix A for the different modes. In the self-consistent Born approximation (SCBA), the self-energies due to electron-phonon interaction are of the form:

Σ< (t) = i q,p fq,p Ĝ< (t)D < q,p (t) fq,p . (38)
Within the {|ζ α } basis considered here, the nonvanishing coupling terms being diagonal with respect to the axial wavefunctions, the self-energy reads

Σ < αβn (t) = i q,p n f (n,n ) q,p (α)f (n ,n) q,p (β)G < αβn (t)D < q,p (t), (39)
where

f (n,n ) q,p (α) = Ψ α,n | fq,p |Ψ α,n . (40)

Optical phonons

For longitudinal optical (LO) and surface optical (SO) phonons, we assume wavevector-independent optical phonon GFs. Indeed only long-wavelength optical phonons with negligible dispersion are effectively coupled to the relevant electronic states. We can then write

Σ <(iO) αβn (t) = iD < LO (t) n W (n,n ) iO (α, β)G < αβn (t), (41)
where

W (n,n ) iO (α, β) = q f (n,n ) q,iO (α)f (n ,n) q,iO (β). ( 42
)
In the energy domain, it reads

Σ <(iO) αβn (E) = i n W (n,n ) iO (α, β) × dE 2π D < LO (E )G < αβn (E -E ). ( 43
)
If the phonon GF D < LO is taken harmonic, only quantized energy exchanges are allowed. In QDs, it has been shown that energy exchanges strongly differing from the optical phonon energy take place due to simultaneous electron-phonon interaction and anharmonic couplings among phonons 6,15 . In these previous works, the electron-phonon interaction was diagonalized exactly, and the Fermi golden rule was used to calculate scattering rates among polaronic states induced by anharmonic couplings. Here, instead, the SCBA is used in order to treat the electron-optical-phonon interaction, and anharmonicity is included in the optical-phonon GFs. As already checked in Ref. 14, the SCBA very well reproduces the exact polaron formation, especially for temperatures T verifying k b T < E LO where E LO is the optical phonon energy. The calculation of anharmonic GFs is presented below.

Acoustic phonons

In contrast to optical phonons, the acoustic phonons are assumed to be harmonic and to have a linear dispersion. The corresponding self-energy can be expressed as:

Σ < αβn (t) = i n K (n,n ) α,β (t)G < αβn (t), (44)
where

K (n,n ) α,β (t) = q,LA D <(0) ωq (t)f (n,n ) q,LA (α)f (n ,n) q,LA (β) (45)
F. Green's functions of anharmonic phonons

We now calculate the anharmonic phonon GFs involved in the electron-optical-phonon self-energy. We retain only the cubic couplings in the anharmonic terms, and thermal equilibrium population of phonon modes is assumed. As we consider only diagonal terms in the phonon's GFs, the Dyson equation reads

D R q,p (E) = 1 D R(0) q,p (E) -1 -Π R q,p (E) , (46)
where Π R is the phonon retarded self-energy. The lesser GF is then calculated using the Keldysh relation:

D < q,p (E) = D R q,p (E)Π < q,p (E)D A q,p (E). ( 47
)
The phonon self-energy is calculated in the Born approximation, taking into account cubic anharmonic terms. The phonon lesser self-energy reads

Π < q,p (t) = i q1,p1,q2,p2 |V a (q, q 1 , q 2 )| 2 D <(0) q1,p1 (t)D <(0) q2,p2 (t),
(48) in which we have taken the harmonic phonon GFs in the right hand side. The cubic anharmonic coupling V a is expressed in terms of the Grüneisen constant as reported in Ref. 15. The phonon spectral and retarded self-energies are then given respectively by

Γ j (E) = i(Π < j (-E) -Π < j (E)) (49) Π R j (E) = 1 2 H(Γ j )(E) - iΓ j (E) 2 (50)
where H denotes the Hilbert transform. Applying the Dyson equation, we obtain the phonon lesser GF:

D < j (E) = 4E 2 j Π < (E) E 2 -E 2 j + E j H(Γ j )(E) 2 + E 2 j Γ 2 j (E) . (51)

G. Elastic scattering

Disordered potentials due to randomly distributed charged impurities, alloy disorder, rough interfaces and rough surfaces break the ideal symmetry of the structure, and couple the dynamics in the various directions. In particular, they couple the 1D motion along the z superlattice axis to the the lateral motion. In planar QW heterostructures, as the lateral dispersion forms a continuum, the static disorder produces an incoherent and irreversible evolution within the z-electron basis. In contrast, in the limit of purely 1D heterostructures, where only one lateral state is involved, the evolution is completely reversible, as elastic couplings induce a unitary evolution within the z-electron basis. A treatment of elastic scattering beyond the SCBA is thus required in these low-dimensional structures.

In previous studies, the surface roughness in nanowires has been treated coherently by random generation of surfaces [START_REF] Wang | [END_REF]26 . Here, we consider in addition interface roughness, impurity scattering and alloy scattering. Instead of generating microscopic random realizations of all surfaces, interfaces, charge positions, and alloy atom positions, we develop a simpler approach based on their correlations. We believe that it is sufficient for low or moderate disorder effects, which are anyway often not known precisely from the microscopic point of view. The disordered potentials are split into diagonal and off-diagonal parts with respect to the lateral eigenstates:

V diag e = n φ n | Ve |φ n , (52)
V off-diag e = n =m φ n | Ve |φ m . ( 53
)
The diagonal component is treated coherently, while the off-diagonal part can only be treated incoherently since we consider only diagonal elements of the GFs with respect to the lateral states.

1.

Elastic scattering self-energies

The off-diagonal elastic coupling component V off-diag e , which does not conserve the lateral state, is treated within the SCBA. The corresponding self-energy reads:

Σ < αβn (t) = i n =n V nn (α)V n n (β) G < αβn (t), (54)
where

V nn (α) = Ψ α,n | Ve |Ψ α,n . (55)

Elastic scattering coherent terms

In order to treat coherently the diagonal component V diag e , we first calculate the following covariance matrix

M n (α, β) = V nn (α)V nn (β) . ( 56
)
We then generate random potentials with correlations that obey this calculated covariance. To this purpose, we need to rotate the covariance matrix in its principal component basis. More precisely, the M n symmetric matrix can be diagonalized in

M n = R n D n R + n , (57)
where D n is a diagonal matrix and R n is a unitary matrix. We then generate random diagonal matrices d rand n in the following way: for each diagonal element of index i, we generate a d n (i, i) random number obeying a Gaussian distribution with a variance D n (i, i). The random potential is then obtained by a back transformation into the original basis:

V rand e = n |φ n V rand e,n φ n |, (58)
V rand e,n = R n d rand n R + n . (59)
The V rand potential is then included in V c in the Dyson equation (32) in order to be treated coherently. The all NEGF calculation is then made for different generated random potentials until the distribution of the calculated observables reached the desired accuracy.

H. Field-periodic boundary conditions

We assume the nanowire SL to contain an infinite number of periods of length L. The possible formation of electrical field domains is not considered in this work 3,27 . Instead, the electric field is assumed to have the same periodicity as the SL. The field-periodic boundary condition for the mean-field component of the electrostatic potential V m-f reads

V m-f (z + pF L) = V m-f (z), ( 60
)
where p is an arbitrary integer. Note that the externally applied electric field F is already included in H 0 . The Poisson equation reads

∂ 2 V m-f ∂z 2 = - e ε 0 s (ρ e (z) + ρ d (z)) , (61)
where ρ e (z) and ρ d (z) are the charge densities of carriers and dopants respectively. From Eq. 61, this field-periodic boundary assumption implies that the electron probability distribution ρ e (z) is periodic. In the length gauge used here, the periodicity condition for the GFs reads

G α+pM,β+pM,n (E) = G α,β,n (E + pF L) ( 62
)
where M is the number of minibands. Note that in nanowires, the SL periodicity can be broken by the inclusion of the coherent disordered potentials (II G 2). In order to investigate such effects, the simulated period can be increased to several SL periods. In this case, the period length L in Eqs. 60 and 62 is replaced by L p = m p L where m p is the number of SL periods included in one computational period. In order to know the needed number of periods, the integer m p is increased until the calculated observables tend towards constant values.

I. Electron density and current

The expectation values of the various observables can be derived using the relation between the density matrix and the lesser GF at equal times (Eq. 19). The electron probability distribution is given by:

ρ e (z) = -2ie α,n dE 2π G < ααn (E)|ζ α (z)| 2 , ( 63
)
where the factor 2 arises from the spin index. Using the expression of the current operator in terms of the restricted operators (Eq. 13), the current reads

J z = - 2e V α,β β|[ ĥ0 , ẑ]|α n dE 2π G < αβn (E). ( 64
)

J. Numerical details

The self-consistent problem formed by the Keldysh equation, the Dyson relation, the expression of the various self-energies and the electrostatic Poisson equation is solved iteratively. Starting from an initial guess of the GFs, we calculate iteratively (i) the lesser Σ< and upper Σ> self-energies; (ii) the retarded self-energy (Eq. 36); (iii) the coherent part of the interaction, comprising the mean-field electrostatic potential and the random potential; (iv) the retarded GF from the Dyson equation; (v) the lesser GFs from the Keldysh relation. This procedure is repeated until a self-consistent solution is reached, with two convergence criteria based on both the lesser GF and the calculated current. Before each step (i), the GFs are renormalized in order to fulfill the charge neutrality for each period. It is checked that this renormalization factor converges accurately towards unity as the convergence is achieved.

The iterative steps (i) and (ii) are performed numerically in the time-domain. It is advantageous since the self-energies within the SCBA involves the product in the time domain. In contrast, energy domain calculations would involve numerically expensive convolutions, as a broad continuum of anharmonic optical phonons is considered. On the contrary, the iterative steps (iv) and (v) are more easily computed in the energy domain, so that Fourier transforms of the GFs and the self-energies are used respectively in between steps (v)-(i) and (ii)-(iv). Though adaptive energy grids have been shown to be useful in order to resolve GF peaks with a reduced number of energy points 17 , a homogeneous energy grid is used here in order to make use of a fast Fourier transform algorithm.

A maximum coherence length l c is set in the calculations, so that we consider only the GF terms G α,β,n for |ζ α -ζ β | ≤ l c . Due to the field-periodic boundary conditions, the GFs G α,β,n are calculated for α belonging to a given period, while β varies away from α within the set coherence length l c . In the simulation of the nanowire SL, we check that l c is taken large enough by increasing it until convergence of the calculated current. In the regime of QW SL, we have verified that there is no difference in the calculated observables considering one or several periods, provided the doping densities remain low enough to prevent the formation of field domain instabilities.

III. RESULTS AND DISCUSSIONS

A. Transport in nanowire SLs: transition from QW to QD SLs

As an application of the theory presented above, we present calculations of the electronic transport in a GaAs/Al 0.3 Ga 0.7 As nanowire SL. One period consists in 5 nm GaAs (well) and 5 nm Al 0.3 Ga 0.7 As (barrier). The SL is assumed to be homogeneously n-doped with a concentration of 10 16 cm -3 . Assumed values for surface and interface roughnesses are given in appendix B. The fundamental subband has a width of 9.7 meV. It is separated from the first excited subband by a minigap of 122 meV. In the following calculations only the fundamental subband is considered (the cut-off value E c is taken around 100 meV). Fig. 2 shows the drift velocity as a function of the nanowire diameter for two different electric fields along the SL. For large nanowire diameters, the current density tends toward a constant value. This is interpreted as the 2D regime for the lateral motion. In contrast, for smaller diameters, the quantum wire regime is reached, and the current density varies with the lateral confine-

1 0 0 1 0 1 E 1 -E 0 ( m e V )
FIG. 2: (Color online). Drift velocity in a nanowire superlattice as a function of the nanowire diameter for various electric fields. The superlattice period consists of 5 nm GaAs followed by 5 nm AlGaAs. The temperature is set to 300 K. ment. In order to confirm this interpretation, we plot in Fig. 3 the quantities

A lat (E) = 1 s n A α,α,n (E) (65) n lat (E) = - 1 s Im n G < α,α,n (E) , (66)
which represent respectively the spectral function and the electron density summed over the lateral states; s = πR 2 is the nanowire cross-section and α is the single basis state per period considered. For a nanowire diameter of 150 nm, the smooth shape of the spectral function indicates that the lateral quantization effects are overcome by broadening effects. Its flat shape is characteristic of a 2D subband. On the other hand, for nanowire diameters below 17 nm, the current displays relatively small variations in Fig. 2. This corresponds to a quasi-0D regime for the lateral motion, where only the ground lateral state is occupied, as shown in Fig. 3 for a nanowire diameter of 15 nm. In between these two limit regimes, large variations of the current densities are calculated in Fig. 2, corresponding to an intermediate dimensional regime where several but still distinct lateral states are involved, as depicted in Fig. 3 for nanowire diameters of 30 and 50 nm.

Fig. 4 shows the calculated velocity-voltage characteristics for various nanowire diameters. For a nanowire diameter of 150 nm, where the 2D lateral regime is reached in good approximation, the standard Esaki-Tsu behavior is observed: negative differential velocity (NDV) is expected to occur when the carriers are excited beyond an inflection point in the miniband dispersion 1 . For smaller diameters, the characteristics deviate from this simple standard behavior. The most remarkable difference is the appearance of a large peak when the voltage drop per period matches the LO-phonon energy ( ω LO =37 meV in GaAs). In addition, smaller resonances become visible at ω LO /2 and ω LO /3. The interpretation is the following: in the 2D regime, the lateral motion acts as a continuous energy reservoir, which renders the 3D energy conservation laws of the scattering mechanisms not directly visible on the SL transport properties. In particular, though LO phonons provide the most efficient inelastic scattering processes, their quasi-monochromatic nature is not evidenced on the velocity-voltage characteristics. In contrast, in the limit of 1D SLs with 0D lateral motion, the energy is exchanged directly between the quantized SL levels (i.e. the Wannier-Stark states) and the phonons. Hence resonances in the transport occur when two consecutive Wannier-Stark levels are separated by the LO-phonon energy ω LO . The smaller resonances around ω LO /2 and ω LO /3 are attributed to LO-phonon assisted transport between states distant of 2 and 3 periods respectively.

In the inset of Fig. 4, we observe an increase in the electron mobility with decreasing NW thickness. This is attributed to the decrease of the efficiency of the scattering processes with decreasing dimensionality. In addition, the position of the first NDV peak shifts to smaller voltage values with decreasing NW thickness. This can be understood within the Esaki and Tsu semi-classical model 1 , where the NDV peak is predicted to occur for edF τ / = 1 (τ being the scattering time). Beyond this critical point, the Esaki-Tsu model predicts that electrons are accelerated beyond an inflection point in the energy-momentum dispersion relation and hence experience a negative effective mass. This edF = /τ critical point can also be interpreted as the transition between the domain of validity of miniband transport and Wannier-Stark hopping models 3 . In both pictures, the overall reduction of scattering processes with decreasing dimensionality explains qualitatively the observed shift of the NDV peak towards lower voltages.

B. Role and nature of elastic scattering processes

In order to analyze the role and the nature of elastic scattering processes, simulations with different treatments of elastic scattering are shown in Fig. 5. First, we discuss the comparison between full calculations with the one neglecting elastic scattering. These two calculations notably differ for thick NWs. In contrast, for thin NWs, the relative difference in the calculated transport is very small, indicating that elastic scattering only plays a minor role. This strong reduction of elastic scattering processes with decreasing dimensionality is consistent with the intuitive expectation that elastic mechanisms are suppressed by the discretization of energy levels. Now we compare the calculations using the standard SCBA, the "selective SCBA" defined above and the full calculation. In the QW limit the standard SCBA, the "selective SCBA" and the full calculation are found to give the same results. In contrast, in the quantum wire regime, the standard SCBA totally differs from the two other models. The reason is that the standard SCBA fails when the fluctuations of the energy levels become larger than their linewidths. In contrast, in the present full calculation, the main coupling elements, which are diagonal with respect to the lateral eigenstates, are treated coherently. If we now compare the full calculation to the selective SCBA calculation, we can see that the selective SCBA provides an excellent approximation to the full calculation, down to 15 nm thin nanowires. Our interpretation is that the typical energy fluctuations of the levels stay smaller than the miniband width and are too weak to induce localization effects.

C. Role of phonon anharmonicity

In Fig. 6 we show a comparison of the calculations with and without the inclusion of phonon anharmonicity for various NW diameters. While the inclusion of anhar-monicity is found to be almost irrelevant for thick NWs in the QW SL limit, striking differences are observed for thin NWs. Indeed the direct LO-phonon emission is gradually suppressed with decreasing dimensionality towards 0D systems, while the mechanism of anharmonic polaron decay becomes dominant. For very thin NW diameters corresponding to the pure 1D SL limit, note that the simulation without phonon anharmonicity does not converge due to the strong reduction of the linewidths below our numerical possibilities in energy resolution. Hence our calculations show a clear transition in the dominant energy dissipation pathway as the lateral dimensionality is reduced, from direct LO-phonon emission to polaron anharmonicity.

IV. CONCLUSION

Within the NEGF formalism, we have developed a theory of electronic transport in nanowire SLs. This model allows the calculation of transport in a wide range of NW SL thickness, for lateral dimensionality ranging from 0D to 2D, i.e. from QD SL to QW SL. We have studied how the transport evolves with this change of lateral dimensionality. While velocity-voltage characteristics are dominated by the standard Esaki-Tsu NDV in QW SLs, electron-phonon resonances appear when the lateral quantum regime is reached. In addition, the electron mobility increases and the NDV peaks are strongly shifted to lower electric fields. This is accompanied by important changes in the dominant scattering mechanisms. In QW SLs, the transport is mainly controlled by LO-phonon scattering and elastic scattering processes. In contrast, when the lateral dimensionality is reduced, these two mechanisms are progressively suppressed, and anharmonic polaron decay becomes dominant.

Interaction with longitudinal optical modes

The polar interaction between electrons and optical phonons confined in a cylindrical nanowire is treated within the dielectric continuum model 28 . The zonecenter LO phonons involved in the couplings with electrons are assumed to be monochromatic with frequency ω LO . The interaction with the longitudinal optical (LO) modes of the nanowire reads

Ĥe-LO = m,l,qz f LO m,l,qz (ρ, θ, z) ÂLO m,l,qz , (A1)
where

ÂLO m,l,qz = â+ m,l,qz + â-m,l,-qz (A2)
is the sum of phonon creation and annihilation operators with opposite momentum. The form factor reads

f LO m,l,qz (ρ, θ, z) = C LO m,l,qz J m χ l m ρ R e imφ e iqzz , (A3)
where

|C LO m,l,qz | 2 = e 2 ω LO 2πε 0 LJ 2 m+1 (χ l m )(χ l2 m + R 2 q 2 z ) 1 ε ∞ - 1 ε s .
(A4) ε 0 is the vacuum permittivity; ε s and ε ∞ are respectively the static and high frequency relative permittivities.

Interaction with surface optical modes

The interaction with the surface optical modes (SO) reads 28 Ĥe-SO = m,l,qz C SO m,l,qz I m (q z ρ)e imθ e iqzz ÂSO m,l,qz , (A5)

where

|C SO m,l,qz | 2 = e 2 ω m,l,qz 2πε 0 LI m (q z R)(I m-1 (q z R) + I m+1 (q z R)) × 1 q z R 1 ∞ - - 1 s - (A6) and = - I m (q z R)((K m-1 (q z R) + K m+1 (q z R)) K m (q z R)((I m-1 (q z R) + I m+1 (q z R)) ext (A7)
where ext is the relative permittivity outside the nanowire.

ω 2 m,l,qz = 1 + 0 -∞ ∞ - ω 2 TO (A8)

Interaction with acoustic phonons

We consider the deformation potential interaction with acoustic phonons confined laterally in the nanowire. In principle, the boundary condition couple longitudinal acoustic (LA) and transverse acoustic (TA) modes. Coupled LA-TA modes interacting with electrons have been considered by Yu et al 29 . However, there, only phonon modes with axial symmetry have been considered. We expect this approximation to be valid in the thin nanowire limit but to completely fail for large nanowires where LA phonons do not necessarily propagate parallel to the nanowire axis. Here, instead, we neglect the coupling between LA and TA modes (i.e. between compressive and shear modes) at the surface of the nanowire, and consider the interaction of electrons with all the LA modes. This treatment is exact in the bulk limit and is expected to provide a reasonable approximation for nanowires of moderate size, as far as the NW diameter is large compared to the lattice constant. We believe such assumption for the boundary conditions does not change significantly the results for the electron transport properties, as (i) the nanowire thicknesses studied here are relatively large compared to the atomic scale and (ii) the efficiency of the acoustic phonon scattering mechanism is always much weaker than other scattering mechanisms. The displacement fields for purely compressive phonon modes can be written as

u = ∇φ (A9)
We consider the free boundary condition, so that ∇u = 0 at the surface and φ can be written as

φ = m,l,qz C m,l,qz J m χ l m ρ R e imφ e iqzz (
with the normalization condition:

d 3 r|u m,l,qz (r)| 2 = 2µω m,l,qz , (A12)
µ is the material density and ω m,l,qz the phonon frequency.

The integration using the form of Eq. A10 gives:

|C m,l,qz | 2 = 2πLµω(χ l2 m + R 2 q 2 z )J 2 m+1 (χ l m ) (A13)
The Hamiltonian describing the electron-LA-phonon interaction is taken of the form Ĥe-LA = D ∇.û, (A14)

where D is the deformation potential of the band of interest. Here it gives:

Ĥe-LA = D(χ l2 m /R 2 + q 2 z ) φ, (A15) Imperfections and intrinsic atomic disorder break the structure symmetry and the crystal periodicity, thus responsible for elastic scattering processes. We consider below the effects of surface roughness of the nanowires, interface roughness of the superlattice heterostructures, randomly located charged impurities and alloy disorder.

Interface roughness

We note z 0 i the mean z-coordinate of the ith interface. The deviation from this average position at the in-plane coordinate r is denoted by δz i (r). The interface roughness can be statistically characterized by (i) the distribution of the interface position along the heterostructure growth axis, and (ii) its in-plane autocorrelation. All interfaces are assumed to obey the same statistics. We note P (z) the distribution of the δz i interface deviation and σ its root mean square. We note C(r) the in-plane autocorrelation form factor so that:

δz i (r 1 )δz i (r 2 ) = σ 2 C(|r 2 -r 1 |) (B1)
We note ∆V the heterostructure band offset. The coupling correlations involved in the calculation of elastic scattering processes read:

V (ir) nn (α)V (ir) n n (β) = ∆V 2 Y nn i W (i) αβ , (B2)
where the lateral form factor reads Y nn = dr 1 dr 2 φ n (r 1 )φ n (r 2 )C(|r 2 -r 1 |), (B3) and the axial form factor for the ith interface reads

W (i) αβ = z 0 i+1 z 0 i-1 dz i P (z i -z 0 i )w (i) α (z i )w (i) β (z i ), (B4)
with

w (i) α (z i ) = [z 0 i ,zi] dz|ζ α (z)| 2 . (B5)
In the numerical study, P (z) is taken as a Gaussian distribution with a standard deviation σ = 0.15 nm and the autocorrelation C(r) decreases as an exponential e -r/λ with a correlation length λ = 8 nm.

Surface roughness

We consider small fluctuations of the surface around its ideal cylindrical position. The surface of the nanowires are considered as infinite barriers for the electrons. Motivated by the fact that a homogeneous change δR in the nanowire diameter would induce a change ∂En ∂R δR in the n-th lateral state energy, we make the following assumption for the surface roughness couplings 30 :

n|V (sr) (θ, z)|n = ∂E n ∂R ∂E n ∂R δR θ,z (B6)
We assume an exponential autocorrelation of the form In the numerical study, the standard deviation is taken as δR = 0.15 nm and the surface correlation length is taken as L s c = 8 nm.

  FIG. 1: (Color online). Schematic of superlattices with different lateral geometries. (a) shows a standard quantum well superlattice consisting of 2D (x,y) infinite layers of various materials stacked periodically along the z direction. (b),(c) and (d) represent nanowire superlattices of different lateral sizes. Dotted arrows are sketched of classical trajectories of charge carriers through the various superlattices under an electric field applied along the z direction. Scattering processes occur due to various effects (interaction with phonons, structural disorder...). It illustrates the transition from transport involving 3D scattering processes towards a pure 1D motion with increasing quantum lateral confinement. Cylindrical coordinates (ρ,θ,z) used in the calculations are indicated for the (b) nanowire.

  o w i r e d i a m e t e r ( n m )

  FIG. 3: (Color online).Spectral function A lat (E) (upper panel) and electron density n lat (E) (lower panel) for various nanowire diameters (see text for definitions). The electric field is 20 kV/cm and the temperature is 300 K.

FIG. 4 :

 4 FIG. 4: (Color online). Drift velocity-voltage characteristics at 300K for various diameters of the nanowire superlattice.

  FIG. 5: (Color online). Drift velocity-voltage characteristics at 300K for different treatments of the scattering mechanisms.

  FIG.6:(Color online). Drift velocity-voltage characteristics at 300K with and without the inclusion of phonon anharmonicity.

  u m,l,qz ÂLA m,l,qz ,

  ,qz (ρ, θ, z) ÂLA m,l,qz , (A16)where the form factor reads:f LA m,l,qz (ρ, θ, z) = C LA m,l,qz J m χ l m ρ R e imφ e iqzz , (A17)|C LA m,l,qz | 2 =where c s the sound celerity. The dispersion relation readsω m,l,qz = c s χ l2 m R 2 + q 2 z . (A19)Appendix B: Elastic scattering mechanisms

δR(θ 1

 1 , z 1 )δR(θ 2 , z 2 ) = δR 2 e -|z2-z1|/L s c × e -R| arg[e i(θ 2 -θ 1 ) ]|/L s c (B7) in which we have factorized the axial and azimuthal correlation terms in order to simplify the calculations. Finally we find that the correlations of the coupling terms are given byV e -πR/L s c if m n = m n R[1-e -πR/L s c (-1) mn -m n ] πL s c [(mn-m n ) 2 +(R/L s c ) 2 ] if m n = m n ,(B9) and the axial form factor reads W (sr) αβ = dz dz |ζ α (z)| 2 |ζ β (z )| 2 e -|z-z |/L s c . (B10)

V. ACKNOWLEDGMENTS

This work has been supported by the Alexander von Humboldt foundation and the Austrian Science Fund FWF through Project No. F25-P14 (SFB IR-ON). S. Birner, M. Branstetter, C. Deutsch, P. Greck, G. Koblmüller, M. Krall, T. Kubis, S. Rotter, K. Unterrainer, and P. Vogl are gratefully acknowledged for fruitful discussions.

Appendix A: Electron-phonon interactions in cylindrical nanowires

In this appendix we give the form factor of the couplings of electrons with longitudinal-optical (LO), surface-optical (SO), and longitudinal-acoustical (LA) phonon modes confined in a cylindrical nanowire. We consider the full confinement of phonons by the nanowire surfaces but neglect the effect of the SL interfaces.

Charged impurities

The random location of charged impurities is also a source of elastic scattering. For simplicity of the calculation and for clear comparison with the 2D limit, we assume that the filling factor of the nanowires is unity. The Coulomb potential created by ionized impurities at the position r d reads

where L s is the screening length. Following Nelander and Wacker 31 , L s is calculated within the Debye screening model for a 3D homogeneous electron gas with the same average electron density and at the lattice temperature. After performing an in-plane Fourier transform this potential reads 32 :

The coupling correlation terms are obtained by averaging over the different possible position of the charged impurities:

After some algebra we find:

B14) where F k are Hankel transform coefficients:

with m n,n = |m n -m n |, and the axial form factor reads:

(B17)

Alloy disorder

We consider a ternary alloy of the form A x B 1-x C. We assume a random and uncorrelated distribution of the atoms A and B. We assume that, at the atomic scale, the carriers experience local band offset potentials V AC or V BC with probability x and 1 -x, respectively 16,[START_REF] Bastard | Wave mechanics applied to semiconductor heterostructures[END_REF] . Here these potentials are taken as the one of the binary compounds. The mean potential reads

while the alloy scattering potential is defined by:

The variance of the local band offset is then:

We denote V the volume occupied by a binary pair of atoms (V = a 3 /4 for zinc-blende structures, a being the lattice constant). The alloy disorder coupling terms read

where x(z) is the z-dependent alloy composition, and where the axial and lateral form factors read respectively:


Karim Ferhat

Arnaud Ralko

Phase diagram of the 1/3-filled extended Hubbard model on the Kagome lattice

We study the phase diagram of the extended Hubbard model on the kagome lattice at 1/3 filling. By combining a configuration interaction approach to an unrestricted Hartree-Fock, we construct an effective Hamiltonian which takes the correlations back on top of the mean-field solution. We obtain a rich phase diagram with, in particular, the presence of two original phases. The first one consists of local polarized droplets of metal standing on the hexagons of the lattice, and an enlarged kagome charge order, inversely polarized, on the remaining sites. The second, obeying a local ice-rule type constraint on the triangles of the kagome lattice, is driven by an antiferromagnetically coupling of spins and is constituted of disconnected six-spin singlet rings. The nature and stability of these phases at large interactions are studied via trial wave functions and perturbation theory.

I. INTRODUCTION

Lattice systems with geometrical frustration are at the center of intense theoretical and experimental research activity of current condensed matter, thanks to the extremely rich exotic phases encountered. One of the most important features of frustrated systems is the presence, in the classical limit, of highly degenerate ground states. These states, under the presence of quantum fluctuations, yield unconventional phases generic either in fermionic or bosonic systems such as spin liquids 1,2 , valence bond physics 3,4 , pinball phases 5,6 , fractionalized defects 7 , and even supersolids 8,9 . Hence, frustration is at the crossroads of many disciplines, ranging from quantum antiferromagnets to optical lattices, via itinerant electrons and bosonic systems.

Intensively studied in quantum magnetism for exotic phases and spin liquid ground states, the role of the frustration is more and more considered in fermionic systems [5][6][7]10,11 . Among the typical two-dimensional lattices possessing a geometrical frustration, one can mention the triangular, checkerboard 8,12 , Cairo-pentagonal 13,14 , and kagome [15][16][17][18] . The latter is the most frustrated two-dimensional lattice and, as a function of the electron filling, possesses several induced local constraint phases. At half filling, under strong Hubbard repulsion, several possible ground states have been proposed for the Heisenberg model [START_REF] Mendels | Introduction to Frustrated Magnetism[END_REF] , ranging from valence bond solids [20][21][22][23][24] to various spin liquids such as a U(1) algebraic 25 , a triplet-gapped 26 and a Z 2 topological state with anyonic excitations [27][28][29] . At intermediate interaction and for various fillings, original electronic features have recently been pointed out, such as peculiar Mott transition 30,31 , anomalous quantum Hall effects 32,33 , Fermi surface instabilities 34 , competing orders 35 , and even anyonic states 36 , to name a few. Surprisingly, at lower filling, only a few works have been done so far. In particular, the 1/3-filled system takes a special role since, in the non-interacting limit, the Fermi sea is at the Dirac cones at the corner of the Brillouin zone (the K point) and the ground state (GS) is a semi-metal (SM). With interactions, let us cite studies on topological phases 10 , Mottinsulator transitions for spinless fermions 37 or fractionalization of defects with an extra bow-tie term 7 . Under the presence of the Hubbard interaction, spin charge density waves with antiferromagnetic spin stripes are stabilized 10 , but this phase was enforced to preserve the translations. With nearestneighbor interaction, a local constraint emerges, often called the ice-rule, very reminiscent of those observed in spin-ices. Thus, the ground state manifold is macroscopically degenerate with a finite entropy at zero temperature. This has to be connected with the physics of hard-core bosons which arises in corner-sharing frustrating unit lattices such as the checkerboard 8,18 , the pyrochlore and the Cairo-pentagonal lattices 13,14 for example. On such a manifold, the quantum fluctuations select a resonating plaquette phase by an order-to-disorder mechanism through a perturbation theory of the third order similar to that of Resonating Valence Bond (RVB) physics, and the local icerule constraint enforces the system to have exactly one doubly occupied site per triangle 7 . Despite the presence of all the ingredients -frustration, local constraint, ground state degeneracy, etc-needed to observe a very rich and exotic physics, the phase diagram of this model is still lacking in the literature, even at the mean-field level.

In this paper, we provide a theoretical study of the phase diagram of the extended Hubbard model at 1/3 filling at low and intermediate on-site and nearest neighbor interactions. In Sec. II are presented the model and methods, as well as the physical quantities used to determine the phases. In particular, we briefly introduce the configuration interaction (CI) method [START_REF]Relativistic and Electron Correlation Effects in Molecules and Solids[END_REF][START_REF] Friedman | [END_REF][40][41][42][43] that we combine to an Unrestricted Hartree-Fock in order to partially bring back the correlations lost in the mean-field decoupling. Section III gives an overview of the thermodynamic limit phase diagram and a brief description of the encountered phases, emphasizing two original ones, the pinned metal droplet phase (PMD) and a spin charge density wave (SCDW). The identification of the phases, the extrapolation to the thermodynamic limit, and the stability of the phases upon the correlations by the CI are discussed in Sec. IV. Finally, the structure of the phases beyond the UHF and CI is discussed in Sec. V. We present a variational theory that helps in the understanding of the mechanism stabilizing the two peculiar phases, one being the presence of local polarized droplets of metal isolated by pins, and the other consisting of disconnected hexagons of six-spin singlets.

II. MODEL AND METHODS

The model considered in this paper is the extended Hubbard model on the kagome lattice:

H = -t i,j (ξ † i ξ j + h.c) + U i n i,↑ n i,↓ + V i,j n i n j (1)
where

ξ † i = (c † i,↑ , c † i,↓
) is a spinor of fermion (electron) creation operators at site i, n i,σ = c † i,σ c i,σ is the number of electrons with spin σ and n i = n i,↑ + n i,↓ is the total density. U is the on-site Hubbard term, V is the nearest neighbor Coulomb repulsion and t > 0 is the hopping term. Throughout the paper, the lowest band is completely filled, corresponding to the 1/3 filling. As reported in the literature, at small interactions U and V , the system is a semi-metal (SM), while at large V (keeping U small) a resonating plaquette charge density wave (CDW) is stabilized 7,10,11,37,44 . Indeed, in this large V limit, the system is governed by a local constraint enforcing each triangle of the lattice 45 to have exactly only one doubly occupied site and the others empty. This ensures the ground state, at the electrostatic limit (t = 0), to be macroscopically degenerate. Then, a ring exchange of particles on a hexagon of order t 3 /V 2 lifts this degeneracy and selects the CDW. How those two phases behave upon intermediate values of U is one of the questions addressed in this work.

The kagome lattice with periodic boundary conditions is defined by its linear size L, and has three sites per unit cell, called A, B and C in Fig. 1. The number of sites is given by N = 3L 2 . Most of the computations of the ground states performed in this paper have been done by (i) using an Unrestricted Hartree-Fock (UHF) and (ii) partially bringing back on the top of this mean-field solution the correlations by using a Configuration Interaction on single and double excitations (SDCI) of the UHF Slater determinant [START_REF]Relativistic and Electron Correlation Effects in Molecules and Solids[END_REF][START_REF] Friedman | [END_REF][40][41][42][43] . As a reminder, in the standard UHF approach, the interaction terms n i,σ n j,σ ′ are decoupled as a Hartree term n i,σ n j,σ ′ + n i,σ n j,σ ′ -n i,σ n j,σ ′ and a Fock term

c † i,σ c j,σ ′ c † j,σ ′ c i,σ + c † i,σ c j,σ ′ c † j,σ ′ c i,σ - c † i,σ c j,σ ′ c † j,σ ′ c i,σ
. The sets n i,σ and c † i,σ c j,σ ′ , for i and j ∈ [1, N ], are computed self-consistently, as well as the single electron basis wave functions. The configuration inter-action (CI) method [START_REF]Relativistic and Electron Correlation Effects in Molecules and Solids[END_REF][START_REF] Friedman | [END_REF][40][41][42][43] , then consists in constructing an effective Hamiltonian H eff ij by projecting the microscopic Hamiltonian H in a subspace made of Slater determinants {φ i }, constructed in this single electron basis obtained by the UHF. If all Slater determinants are kept in the basis, then the full Hilbert space is spanned and H eff ij is nothing else but H. This is usually referred to as the full CI (FCI) in the literature [START_REF]Relativistic and Electron Correlation Effects in Molecules and Solids[END_REF] . Obviously, in this case, the limitations are the same as those for the exact diagonalizations (EDs), the Hilbert space is increasing exponentially, and only small cluster sizes are available. An advantage to working with the Slater basis however, is that a simple energy criterion can be exploited. Indeed, at small and intermediate values of V and U , the correlations are expected to be well taken into account by single and double excitations from the Fermi sea. Hence, it is reasonable to truncate the basis to those processes and to restrict our effective Hamiltonian to this subspace. This is the so called single-double CI (SDCI). Note that this has been intensively studied in quantum chemistry where it has been shown that, for certain molecules, up to 80% of the correlation energy was restored. Here, in addition to our results, we show in Sec. IV-B an improvement of about 60% for certain phases.

In order to characterize the phase diagram, we define a set of order parameters. The structure factor M (q) = ψ 0 |n(-q)n(q)|ψ 0 /N , where n(q) = n ↑ (q) + n ↓ (q) indicates spontaneous translational symmetry breaking of the ground state. Since at very large Coulomb repulsion V the system is governed by the local ice-rule constraints, we define the local triangle operators

ξ ± = 6 N t w 0 ρ ± A + w 1 ρ ± B + w 2 ρ ±
C with ρ ± i = n i↑ ± n i↓ and w = e i2π/3 . ξ + = 1 for a phase obeying an ice-rule of two particles per triangle with exactly one doubly occupied site, and ξ -= 1 if the ice-rule of the type of two particles with one empty site is realized. Both order parameters are zero in the homogeneous phase, and finite otherwise. To study the internal spin ordered phases, we define the operator ρ - A ρ - B ρ - C , which is finite if there are no vacant sites.

III. OVERVIEW OF THE PHASE DIAGRAM

The phase diagram of the 1/3-filled extended Hubbard model on the kagome lattice extrapolated to the thermodynamic limit from the UHF results up to 3 × 36 × 36 site clusters is depicted in Fig. 2. Four phases have been obtained in the (V, U ) plane, the semi-metal (SM), a charge density wave (CDW), a spin charge density wave (SCDW), and a pinned metal droplet phase (PMD). As already mentioned, the first two phases have already been reported elsewhere 37,44,45 . For the CDW, a perturbative effective model of order t 3 /V 2 maps the present model onto a quantum dimer model on the hexagonal lattice 11,37,44,45 . Indeed, in the classical limit t = 0 and V = 0, the ground state is macroscopically degenerate, and the configurations that minimize the energy obey a rule of two electrons per triangle with one doubly occupied site. A finite t lifts this degeneracy and in the limit t ≪ V the system is effectively described by a hardcore quantum dimer model Phase diagram in the (V, U ) plane. All points are extrapolations at the thermodynamic limit obtained by size scalings. The four phases are the semi-metal (SM), the resonating charge density wave (CDW), the spin charge density wave (SCDW), and the pinned metal droplets (PMDs). Dashed (continous) lines and empty (filled) circles correspond to first-order (second-order) order transitions. Results for the scan line (purple) are detailed in Fig. 3.

(QDM) on the honeycomb lattice 11,37,44 whose ground state consists of resonating plaquettes. The mean-field treatment cannot capture the resonating plaquettes of the QDM, but its solution can be thought of as its electrostatic counterpart 10 . This can be seen in the UHF local densities in Fig. 4. A second order transition between the SM and the CDW is obtained. To the best of our knowledge, the two other phases reported in this work, namely the PMD and SCDW phases, have not been reported so far 46 . Starting from the CDW phase in the large V regime, by increasing U with U > 2V , it is more favorable for the system to forbid double occupancies. This is happening through a first order phase transition. In contrast to the recent paper of Wen et al. [10], since we do not restrict our phases to preserve translational symmetries, we obtain a √ 3 × √ 3 unit-cell phase with six-spin antiferromagnetic hexagonal rings. All spin-singlet rings are disconnected by empty sites, as depicted in Fig. 4-a (see the text for more details). We call this phase the spin charge density wave (SCDW). Its phase transition with the SM is of second order. Eventually, at lower V interaction, the antiferromagnetic rings are not favored anymore, and the system prefers to separate up and down spins via another first order transition. One spin sector crystallizes on one third of the sites in order to form an enlarged kagome structure, the other sector forming disconnected half-filled hexagonal metal rings on the remaining sites. This phase, called the pinned metal droplets (PMDs) phase, since the hexagons (metal droplets) are disconnected by localized electrons (pins), is stabilized thanks to the kinetic energy of those rings, as discussed in Sec. V. Except for the semi-metal, all the phases reported in this phase diagram break the translational symmetry, as a finite M (K) indicates (see Fig. 3). We provide in Sec. IV the complete analysis about the nature of these phases by UHF, and then we describe the SDCI method and apply it to confirm their stability upon the correlations. We have followed the following strategy. First, we identify the phases by the UHF. Then, we check their stability deep in their domains and close to the phase boundaries by constructing a configuration interaction effective Hamiltonian in the single and double excitations (SDCI) Slater basis (see the next section for more details). We also verify that no new phases arise. Finally, in Sec. V, we unravel the mechanism of our phases by constructing a trial wave function and a perturbative theory based on our UFH and SDCI results at large U and V limit. This validates the present phase diagram. V t order parameters FIG. 3: (Color online). Behavior at the thermodynamic limit of the order parameters defined in the text along the scan line depicted in Fig. 2 corresponding to U = 10 -5V : (blue) M (K), (gray) ρ - A ρ - B ρ - C , (green) ξ -, and (red) ξ + . Four distinct domains are observed, separated by continuous (second order) and dashed (first order) lines. Inset: size scaling of the CDW order parameter close to the transition with the SM; a clear change of behavior indicates the location of the phase transition.

A. UHF thermodynamic limit

In this section, we first present the thermodynamic limit (TL) phase diagram, displayed in Fig. 2, extracted from finitesize scalings of the UHF results. Then, we confirm the existence and the stability of all the phases beyond the mean-field level thanks to the single and double configuration interaction (SDCI). We have performed finite size scalings on the four order parameters, M (q), ξ ± , and ρ - A ρ - B ρ - C on clusters up to 3×36×36 sites. To illustrate our results, we consider the scan line displayed in Fig. 2, which passes through all the phases. The results are displayed in Fig. 3. As one can see, size effects are important (inset), as expected for the kagome lattice, but thanks to the large size clusters available, precise size scalings are performed. Our analysis reveals two second order phase transitions, between the SCDW and the SM and between the SM and the CDW. The finite order parameter M (K) (blue) at the TL indicates the spontaneous symmetry breaking of translations at the K point and thus a √ 3 × √ 3 unit cell of nine sites. The quantity ρ - A ρ - B ρ - C (gray) is zero as soon as on each triangle at least one zero magnetization site is present or finite otherwise. We observe that the PMD phase is the only one with a finite TL value, and all sites of each triangle are occupied. Finally, ξ +(-) (red, green) provide insights about the internal charge (magnetization) structure of each triangle. In Table I is reported which order parameter is compatible with which phase. All of them have been identified unambiguously

Phases M (K) ξ + ξ -ρ - A ρ - B ρ - C SM 0 0 0 0 CDW 0 0 PMD SCDW 0
TABLE I: Classification of the possible phases that may occur in the extended Hubbard model on the kagome lattice. Such phases can be distinguished from the density-density correlation M (K), the internal triangular structures in charge (ξ + ) and magnetization (ξ -), and the presence or not of a s z = 0 site per triangle (ρ

- A ρ - B ρ - C ).
within the UHF approach as seen in Fig. 3. The corresponding density maps of the PMD and the SCDW, for the charge and the magnetization, are very close to the ones obtained by the SDCI (see below) depicted in Fig. 4.

B. Configuration interaction: beyond the Hartree-Fock

In order to get further insights about the nature of the phases and their stability beyond the UHF, we have developed a configuration interaction method which aims to bring back the correlations on top of a mean-field single particle solution, here obtained by the UHF (see Sec.II). In Table II are reported the ground state energies obtained by ED (FCI), SDCI and UHF on a 2 × 2 cluster, the largest available by ED at the 1/3 filling. Keeping all single and double excitations above the Fermi sea corresponds to a total of 3489 Slater determinants. The parameters are chosen in the three non trivial phases, the CDW, the SCDW and the PMD. For each of them, the SDCI clearly improves the GS energies, hence showing that correlations have been incorporated.

In our model, if the phases obtained by UHF are not stable under the correlations (for example the UHF is known to overestimate ferromagnetism), we expect the SDCI to stabilize other solutions. We were able to distinguish the four phases identified at the mean-field level. The snapshots of the local densities n i,σ of the PMD and the SCDW phases are depicted in Fig. 4. By turning on the correlations on top of the UHF solution, we found that all the phases are stable up to the largest cluster available by the CI, in our case the 6 × 6 unit cells (108 sites). Note that for such a large cluster, the corresponding basis contains 263226945 determinants, which is not reachable by our computers. In order to reduce this size, it is possible to keep among the single and double excitations only those for which the excitation energy does not exceed a chosen cutoff Λ CI . By adjusting correctly this cutoff, one can find a compromise between the improvement of the GS wave function and the size of the basis we want to deal with. For the 3 × 6 × 6 site cluster, we have kept 1850220 Slater determinants for energy improvement of the order of 10% on the GS energies. On the 3 × 3 × 3 site cluster, for which no cutoff is required (97039 Slater determinants), the UHF energy is 14.4466(8) and 11.4323(0) for the SDCI in the SCDW phase at U = 8 and V = 2 [an improvement of 20.97(5)%]. In the PMD phase at U = 8 and V = 0, the UHF gives -24.9667(2) against -28.7069(6) for the SDCI [14.98(0)%]. In Fig. 4 are depicted the local charge and spin densities obtained by SDCI on the 108 site cluster, for the PMD and the SCDW phases. For information, in our computations, the change in the density values between the SDCI and the UHF never exceeded 3-4%, even close to the phase boundaries. This indicates that (i) these phases are stable upon the correlations, and (ii) the UHF phase diagram at the TL of Fig. 2 captures the essential physics. Let us conclude this section by a remark concerning the phase boundaries of Fig. 2. In this work, we have employed the SDCI to validate the presence of the main phases beyond the UHF and to verify whether or not new phases could be stabilized. Hence, we have focused deep in their domains, and made checking close to the main transitions. The question of the exact location of the boundaries within this approach requires an optimized methodology in selecting the active determinants, in order to reduce the size of the basis while keeping the correlations strong [41][42][43] . This is left for another work.

V. TRIAL WAVE FUNCTION AT LARGE U AND V

Let us now try to understand the characteristics of the PMD and the SCDW phases by considering the large U and V limit.

At large U ≫ V , we find a phase which is the separation between pins (one electron) on each site of an enlarged kagome lattice and balls of three polarized electrons on each hexagon, that we expect to be in a local metal state (delocalized electrons on six sites). Indeed, as a first approximation, one could think that these hexagons are isolated by the pins, as suggested by the real space density maps of the PMD depicted in Figs. 4(c) and 4(f). To validate this picture, we have calculated the energy of the trial wavefunction (TWF)

|ψ = [⊗ i |i ] ⊗ [⊗ h | h ]
where the first product is over single polarized electron states {|i } standing on the sites of the enlarged kagome lattice and the second is over all the remaining hexagons, | h being the ground state of an extended Hubbard model of three inversely polarized electrons on a ring:

H = -t i∈ (c † i c i+1 + h.c) + V n i n i+1 . (2)
By computing the ground state energy and the kinetic operator P = (c + 0 c 1 c + 2 c 3 c + 4 c 5 + h.c.) (the sites are labeled clockwise from zero to five), an inflexion point appears about V ≃ 3t. For V < 3t, the ground state is delocalized (metal). At V = 0, the ground state on each hexagon is a simple metal of three electrons, and in terms of local momenta k on the six-site ring, the wave function is trivially given by the product state of the single electron Bloch states

| h = |k 0 ⊗ |k + ⊗ |k -,
where k 0 = 0, and k ± = ±π/3, with energies E 0 = -2t and E ± = -t respectively. By applying the perturbation theory on |ψ at large U , we find a correction to the energy as

E k 0 = -2t -4 t 2 U , E k ± = -t -3 t 2 U
, and for the pins

E i = -24 t 2
U , for a total trial energy E TWF = -4t -34 t 2 U . In Fig. 5 are depicted the GS energies at V = 0 obtained by UHF and the trial energy obtained from the TWF with metallic hexagons; their agreement is excellent.

We turn now tn the SCDW. In the large U limit, at t = 0, there are no double occupancies and the system, for any finite V , enters an ice-rule constraint of two electrons per triangle. Its ground state manifold is macroscopically degenerate and consists of closed loops of various lengths. When the quantum fluctuations t are turned on, the UHF and the SDCI select the smallest possible loops as observed in Fig. 4-d. From a perturbation theory, there are two types of processes which can lift this degeneracy: an antiferromagnetic coupling of spins of order ∼ t 2 /(U -V ) at the second order, and resonating plaquette terms (ring exchange of three electrons on a hexagon) in ∼ t 3 /V 2 at the third order, similar to those of the QDM 37 5: Agreement between the UHF energy and the trial wave function based on metallic hexagons for the PMD phase. We show several scans from t = 0.5 (darker points) to t = 1 (lighter points). Each scan is made from U = 10 to 50. The dashed line is the energy expected in the PMD phase E = -4t -34 t 2 U .

at U = 0 and V ≫ t. Unfortunately, neither the UHF (meanfield) nor the SDCI (up to double excitations) allow us to distinguish which process is dominant in functions of U and V . All we can say is that, at large enough V and U ≫ V , the Heisenberg mechanism will be obviously favored. But let us go back to Eq. ( 2) and |ψ . The inflexion point in P appearing around V ≃ 3 indicates that ∼ t 3 /V 2 ring exchange terms become dominant for higher V and thus the PMD is destroyed at such large values of V . At the mean-field level, the UHF snapshot would be similar to that of Fig. 4-d, but now the resonating hexagons are no longer the same as the Heisenberg loop 6 discussed above. Hence, two scenarios are possible at the transition of the PMD: (i) the system enters a phase with resonating plaquettes separated by pins based on our trial wave function, and then a further transition to a spin state driven by Heisenberg couplings is stabilized (with no pins), or (ii) the Heisenberg state is selected right away. Our UHF and SDCI results present a phase transition about V ≃ 0.25 up to U = 12t, which is much smaller than V ≃ 3, at which the internal structure of the hexagon is changing. This seems to indicate that scenario (ii) is favored, or in this case, based on the results of Fig. 4-d, a ground state made of six-spin singlet hexagons of energy per site E 6 ≃ -0.467(1) is selected among the degenerate configurations. This is what we call the spin charge density wave (SCDW). Finally, let us compare our results with previous work. In Wen et al. 10, since no translational symmetry breaking is allowed, a stripe phase with disconnected antiferromagnetic chains is obtained. In this case, the quantum fluctuations of the perturbation theory lowers the classical energy by E c = 1 4 -log 2 ≃ -0.443(1). The hexagonal spin chain is then more likely when the symmetry constraint is removed, as confirmed by the snapshots of the SCDW depicted in Fig. 4(c) and 4(d).

VI. CONCLUDING REMARKS

In this paper, we have studied an extended Hubbard model of spinful electrons on the kagome lattice at 1/3 filling. In addition to two already known phases, the semi-metal (SM) and the charge density wave (CDW), we have established the presence of two original phases, the pinned metal droplets (PMDs) and the spin charge density wave (SCDW) stabilized at intermediate U and V .

We have first used an unrestricted Hartree-Fock approach in order to extract the phase diagram at the TL. We have then constructed an effective configuration interaction Hamiltonian by projecting the microscopic model onto a basis of Slater determinants containing the UHF wave function, and all single and double excitations (SDCI). This allows for bringing back partial correlations on top of the UHF solution. For various sets of parameters, the ground-state energy (and related GS) is strongly improved (up to 60% of the correlation energy). We have confirmed that (i) no new phases are stabilized by the correlations and (ii) the two original phases, namely the PMD and the SCDW, are very robust against the correlations.

Finally, to get more insights about the nature of two new phases, we have performed a perturbation theory on a trial wave function. We have shown that the properties of the pinned metal droplets (PMDs) are perfectly captured by the UHF and explained by our trial wave function. It consists of polarized pins (electrons) standing on the sites of an enlarged kagome lattice, separating droplets of metal living on the remaining hexagons. The second important phase, the spin charge density wave (SCDW), is driven by Heisenberg inter-actions and leads to loop-6 ring isolated spin states. Note that we have discussed the possibility of a third phase that could arise in the system, close to the PMD but with resonating plaquettes instead of metal droplets on the remaining hexagons. In some sense, this would correspond to an internal crossover of the hexagon states, the pins remaining unchanged. However, within our approaches, this phase does not seem to be stabilized at the intermediate parameters considered in this work. Instead, the PMD and the six-spin singlet SCDW are the more stable.

We recently became aware of related work by F. Pollmann et. al (Ref. [46]), where the authors have considered an effective model of Eq. ( 1) at large U and V limit. We thank F. Pollmann for stimulating discussions about their results corroborating the existence of the phases reported in this manuscript.

FIG. 1 :

 1 FIG. 1: (Color online). a) Example of a 3 × 3 kagome lattice with three sites per unit cell dubbed A, B and C. b) The Brillouin zone with k1 and k2 being the reciprocal vectors of u1 and u2 and the elementary translational vectors in real space defined in (a). The highest symmetry points Γ; M and K are depicted.

FIG. 2 :

 2 FIG.2:(Color online). Phase diagram in the (V, U ) plane. All points are extrapolations at the thermodynamic limit obtained by size scalings. The four phases are the semi-metal (SM), the resonating charge density wave (CDW), the spin charge density wave (SCDW), and the pinned metal droplets (PMDs). Dashed (continous) lines and empty (filled) circles correspond to first-order (second-order) order transitions. Results for the scan line (purple) are detailed in Fig.3.

FIG. 4 :

 4 FIG. 4: (Color online). Snapshots of local charges (a,c) and spins (b,d) for the PMD phase (a,b) at U = 8 and V = 0 and SCDW (c,d) at U = 8 and V = 2 obtained by SDCI on a cluster of 108 sites. The larger the circle, the higher the density (charge and spin). The colors green and purple correspond respectively to spin sectors up and down. The larger and lighter disks in (a) and (b) illustrate the metal droplets standing on the hexagons of the PMD phase, pinned by localized electrons (green circles in (b)).

TABLE II :

 II Ground state energies of the different phases obtained by the three methods used in this work, here on a 2 × 2 cluster (12 sites) for various parameter sets.

Acknowledgements -A.R. is grateful to F. Mila for insightful remarks about the stability of the phases. We thank S. Fratini for enlightening discussions and his careful rereading of the manuscript. A.R. thanks G. Bouzerar and M.-B. Lepetit for their valuable comments about the methods used in this paper. This work is supported by the French National Research Agency through Grants No. ANR-2010-BLANC-0406-0 NQPTP and No. ANR-12-JS04-0003-01 SUBRIS-SYME.