Dilute gas of ultracold two-level atoms inside a cavity; generalized Dicke model

We consider a gas of ultracold two-level atoms confined in a cavity, taking into account for atomic center-of-mass motion and cavity mode variations. We use the generalized Dicke model, and analyze separately the cases of a Gaussian, and a standing wave mode shape. Owing to the interplay between external motional energies of the atoms and internal atomic and field energies, the phase-diagrams exhibit novel features not encountered in the standard Dicke model, such as the existence of first and second order phase transitions between normal and superradiant phases. Due to the quantum description of atomic motion, internal and external atomic degrees of freedom are highly correlated leading to modified normal and superradiant phases.


I. INTRODUCTION
Progress in trapping and cooling of atomic gases [1] made it possible to coherently couple a Bose-Einstein condensate to a single cavity mode [2]. These experiments pave the way to a new sub-field of AMO physics; many-body cavity quantum electrodynamics. In the ultracold regime, light induced mechanical effects on the matter waves lead to intrinsic non-linearity between the matter and the cavity field [3,4]. In particular, the nonlinearity renders novel quantum phase transitions (QPT) [5], see [6]. Such non-linearity, due to the quantized motion of the atoms, is absent in the, so called, standard Dicke model (DM). Explicitly, the DM describes a gas of N non-moving two-level atoms interacting with a single quantized cavity mode [7]. The interplay between the field intensity/energy, the free atom energy, and the interaction energy leads to a quantum phase transition (DQPT) in the DM. Motivated by the novel phenomena arising from a quantized treatment of atomic motion, it is highly interesting to extend the DM to include atomic motion on a quantum scale, and in particular to analyze how it affects the nature of the DQPT. It is clear that such generalization of the DM results in new aspects of the system properties. The atomic motion is directly affected by the shape of the field induced potentials, which in return is determined by the system parameters and the field intensity. In addition to the terms contributing to the total energy in the regular DM, in this model we have to take into account for the motional energy of the atoms. This enters in a non-trivial way since the interaction energy depends on the motional states of the atoms.
The DM was first introduced in quantum optics to describe the full collective dynamics of atoms in a high quality cavity. The DM Hamiltonian, given in the rotating wave approximation (RWA), reads Here, the boson ladder operatorsâ † andâ create and annihilate a photon of the cavity mode, the Pauliσ ioperators act on atom i, ω, Ω and λ are the mode and atomic transition frequencies and effective atom-field coupling respectively. The DM is a well defined mathematical model for any value of its parameters. Note, however, that whether this model describes faithfully some physical situation for any choice of parameters is not guaranteed. For moderate λ, the DM appropriately describes the dynamics of atoms coupled to the cavity field, and it has been thoroughly discussed in terms of collapse-revivals [8], squeezing [9], entanglement [10] and state preparation [11]. Considerable interest, however, has been devoted to the normal-superradiant phase transition [12,13,14].
In the thermodynamic limit (N, V → ∞, N/V = const.) the system exhibits a non-zero temperature phase transition between the normal phase and the superradiant phase. In the normal phase, all atoms are in their ground state and the field in vacuum, while the superradiant phase is characterized by a non-zero field and a macroscopic excitation of the matter. Its critical coupling and critical temperature are [12] λ c = √ Ωω, As was shown in [13], the critical coupling can be seen as a condition on the atomic density ρ = N/V . Interestingly, the DQPT is of second order nature without the RWA, while it is first order if the RWA has been imposed [15]. The corrections due to the RWA to various physical observables have been considered [14,16].
More detailed analysis about miscellaneous aspects of the DQPT have been presented in numerous publications. Especially various extensions [17,18] as well as approximate methods [19] concerning the DQPT has been outlined. Recently, C. Emary and T. Brandes applied the algebraic Holstein-Primakoff boson representation on the DM. The method turned out to be very powerful and have since then been applied frequently to the DM [20].
Despite the numerous publications on the DQPT, the existence of this phase transition was widely discussed. If the two-level atoms in the DM correspond to atoms in a ground and excited state, and the transition is direct, then quantum mechanics forbids the transition. This can be seen either by realizing the necessity of adding the, so called, A 2 -term to the Hamiltonian, or by employing sum rules to bound the coefficients in the Dicke model to the "trivial" thermodynamical phase [23]. The argumentation of Ref. [23] can be generalized to a quite general no-go theorem for the DQPT [24], but it does not apply if the two level atoms in the DM correspond to atoms in two excited states, such as Rydberg states, or if the transition is not direct.
One valuable step towards an experimental realization of the DQPT was taken in relation with Ref. [21], where typical experimental parameters as well as losses were included. These authors considered the two levels coupled by a non-resonant Raman transition. In such conditions the atom-field coupling λ can be tuned more or less independently of the A 2 term in an effective two-level model, and one can reach the regime of DQPT. This paper, however, considers a situation in which atomic motion can be neglected due to high temperatures, e.g. the standard DQPT. Alternative situation, in which the atomic motion could be neglected, would be to consider quantum dots interacting with a cavity mode [18,22].
In this paper we extend the DM to take into account for atomic motion in a fully quantum mechanical description. The atomic motion then introduces an additional degree of freedom to the problem, leading to novel appearances of the system phase diagrams. The gas is assumed dilute such that atom-atom scattering can be neglected, and that the motion is restricted to one dimension due to tight confinement in the remaining two directions via external trapping. Furthermore, the atoms are assumed trapped by the cavity field itself. Consequently, a normal-superradiant QPT is not possible, since for a vanishing field the atomic trapping capability is lost. However, we may assume a lowest bound of the field such that at least one bound state of the trapping "potential" is guaranteed. This can be achieved by an external pumping of the cavity, which imposes a nonvanishing cavity field. Our research is partly carried out in the adiabatic regime, motivated by the ultracold atoms considered and its justification is numerically verified. In this adiabatic regime, the problem relaxes to solving a 1-D time-independent Schrödinger equation. In particular we study the case of a Gaussian mode profile utilizing this adiabatic method. The situation with a standing wave mode profile is also considered, however using a full numerical rather than adiabatic approach. For a Gaussian profile the number of bound states is crucial for the thermodynamics, and we find great divergences from the regular DM. Among these are the existence of both first and second order QPT's and multiple super-radiant phases. The second, standing wave mode, shows slight similarities to the model of [18] but the QPT is found to be of second order, and the PT survives for zero temperature and finite ω opposite to the regular DQPT.
The paper is organized as follows. In the next section we present the generalized DM which includes the motion of the atom. The adiabatic diagonalization of the single particle Hamiltonians utilized for the Gaussian mode profile is introduced and the general expression for the partition function given. The following Sec. III considers the situation of a Gaussian mode profile. We thoroughly discuss the importance of bound states. Section IV instead considers a standing wave mode profile in a fully numerical fashion. In the appendix, however, we derive analytical expressions in the regimes of tight binding which enables us with various asymptotic properties. Last we conclude with a summery in Sec. V.

II. GENERALIZED DICKE MODEL AND ITS PARTITION FUNCTION
We consider an gas of N ultracold identical two-level atoms, with mass m and energy level separationhΩ, interacting with a single cavity mode with frequency ω. For a low temperature gas we include atomic center-of-mass motion and mode variation. In the dipole and rotating wave approximation, the extended DM becomes Here,p i andx i the scaled center-of-mass momentum and position of atom i respectively, g(x) the effective positiondependent atom-field coupling and V is the mode volume. Throughout the paper we will use scaled variables such thath = m = 1. The case of a single atom is given by the generalized Jaynes-Cummings Hamiltonian studied by numerous authors [25,27].
In the thermodynamic limit we let V → ∞ and N → ∞ such that the atomic density is fixed; ρ = N/V ≡ ρ 0 . The partition function reads where β −1 = T and T is the scaled temperature and the trace is over the field and atomic degrees of freedom. It is convenient to perform the trace of the field in terms of Glauber's coherent states,â|α = α|α . In the thermodynamic limit one may replaceâ → α andâ † → α * in the evaluation of the partition function [12]. In other words; in the large atom number limit the photon ladder operators, or more preciselyâ/ √ N andâ † / √ N , can be treated as commuting operators. Using the fact that atomic operators mutually commute between themselves, for example [x i ,p j ] = iδ ij , the partition function can be written where the integration is over the whole complex α-plane and In the sense of ultracold atoms as considered here, the kinetic energy of the atoms is assumed smaller than the effective potential energy. Provided that the adiabatic potentials do not cross, it is then legitimized to perform an adiabatic diagonalization of the internal states [26]. In this regime, the single particle Hamiltonian relaxes to two decoupled adiabatic ones This approximation will be imposed in the next section considering a Gaussian mode profile. However, in the proceeding section dealing with the standing wave mode such an approach is not advocate, since then the curve crossings between the adiabatic potentials break adiabaticity [26]. The justification of the adiabatic approximation applied to the Gaussian mode profile will be discussed in the end of next Section. Within this regime, the problem has become one of solving for the eigenvalues of two time-independent decoupled Schrödinger equations. The adiabatic Hamiltonians depend solely on the norm |α| and in polar coordinates the angle part can therefore be integrate out. By denoting the eigenvalues E ± n (r) respectively, where r = |α| and n is a quantum number/numbers that can be either discrete and/or continuous, we get the adiabatic partition function (8) Without loss of generality we can choose ρ 0 = 1 as it only scales the effective atom-field coupling. It is worth mentioning that the numerics deal with exponentially large numbers, especially for small temperatures, which restrict the analysis to certain ranges.

A. Derivation of the partition function for transversal motion
A Fabry-Perot cavity has eigenmodes that are, to a good approximation, Gaussian in the transverse and harmonic in the longitudinal direction. Assuming an external deep trap in the longitudinal direction and one transverse direction, we may consider the one dimensional problem in which the atom field coupling has a spatial Gaussian shape. As is well known [25], and seen from Eq. (7), only atoms in the "adiabatic" internal state corresponding to the Hamiltonian h − ad (r) will feel an attractive potential, while the others will be scattered away from the cavity field. We therefore consider only a sub "quasi" Hilbert space containing the bound states E − n (r) of where ∆ x is the transverse mode width. In order to proceed in an analytic way, we make the following approximate ansatz, . (10) The unknown constants are determined from the conditions: (i) The two functions have the same asymptotic values for x → ±∞, (ii) their maximum are the same and (iii) they share the same FWHM. Explicitly this yields The bound eigenvalues of the Hamiltonian are known to be [28] (13) Let us introduce the number of bound states, for a given set of parameters, asÑ and define the function With this, the partition function (8), considering only bound states, becomes By the variable substitution y = r 2 /N we get In the thermodynamic limit, this integral is solved by the saddle point method [29] where C 1 is a constant. Note that y has the meaning of scaled field intensity. One obstacle of the above model already mentioned in the introduction, is the fact that for a shallow potential well the number of bound states will vanish. In this limit, the cavity field can no longer serve as a trap for the atoms. Consequently the ground state is the one of zero atoms, and we cannot have a proper thermodynamic limit N → ∞. We therefore add the constrain of a minimum of one bound state is assumed. This can be met experimentally by including an external driving of the cavity mode, so that the field is non-zero throughout. Thus, the "normal" phase will contain a non-zero cavity field which is sustained by the external pumping. We have numerically checked that this does not introduce any significant changes of the analysis.
For the numberÑ of bound states, we have Naturally,Ñ depends on the system parameters. AsÑ is an integer, a change in the system parameters may bring about jumps between integer numbers ofÑ . This will cause discontinuities in the function g 1 . LettingÑ = 1 we find U 0 (y) > q 2 y . For small fields, y → 0, the potential amplitude U 0 vanishes and the bound states cease to exist. However, for small but non-zero fields, the above inequality may be met for large couplings λ and widths ∆ x .

B. Numerical results
To study Eq. (17), we analyze the parameter dependence of the function Note that f 1 (y) is the free energy per particle, given a scaled field intensity y. Let us briefly discuss characteristics of f 1 (y) before approaching it numerically. The first part arises from the bare field, and it is energetically favorable to have a vanishing field. The second part contains the atom-field interactions plus kinetic and potential atomic energies. The interaction energy enters indirectly into the atomic potential part. Increasing the field amplitude deepens the potential well and therefore lowers the energy, and it is therefore more beneficial to have a large field. The two terms therefore compete, and in particular, the location of the maximum of f 1 (y) depends on the particular system parameters used. Thus, atomic motion, directly related to the shape and depth of the adiabatic potential, is a crucial ingredient for the QPT. If the smallest possible y maximizes the function, the system is said to be in a "normal" phase (in quotes because the field is still non-zero to guarantee at least one bound state), while if a non-minimal y is optimal the system is in a superradiant phase. In the limit of large y, the second term diverges as ln [g 1 (y)] ∼ √ y while the first term goes as ∼ −y, and we conclude that a maximum of f 1 (y) is only obtained for a finite y. These reflections are numerically verified in Fig. 1 showing f 1 (y) for four different couplings λ. We see that there is a critical coupling λ c for which λ < λ c the system is in a "normal" phase and for λ > λ c it is in a superradiant phase.
In Fig. 2 we display the critical coupling λ c as function of ω while keeping the other parameters fixed. In (a) we present two examples for different T and in (b) two examples for different ∆ x . For the plots, the minimum y is taken so that there is at least one bound state in the well. To the left of the curves the phase is superradiant, while to the right it is "normal". The structure of the phase diagram is clearly different from the one of the regular DM in which, at zero temperature, λ c ∼ √ ω. For certain couplings λ s , the system is always in a superradiant phase independent of ω. The location of these resonances are insensitive to the temperature but not to trap width ∆ x . The "sharpness" of these points makes it possible to have a "normal"-superadiant-"normal" QPT by fixing all parameters but the coupling λ which is varied around λ s . This novel feature comes about due to the varying number of bound statesÑ in the well. For say ω ≈ 1 and a weak coupling, the system is in the "normal" phase with only one bound state. As λ is increased the system goes through a QPT into a superradiant phase. A closer numerical study shows that this QPT is caused by the sudden change of going from one to two bound states in the trap. As the coupling is further increased, the discontinuity that arose from the appearance of a second bound state is of less importance and the system reenters the "normal" state. Hence, the presence of a second bound state "forced" the system out of the "normal" phase by inducing a sudden kink/maximum in the function f 1 (y). When the coupling is increased even further, the same may happen again when a third bound state is introduced in the trapping potential. Eventually, however, the system enters a superradiant phase without reentering the "normal" phase. In order to discuss the character of the QPTs we introduce the parameter that maximizes f 1 (y). Thus, I is the scaled field intensity where a discontinuity in it indicates a first order QPT and a discontinuity in its first derivative signals a second order QPT. Here y 0 is the lower bound which assures at least one bound state. The parameter I corresponding to the T = 0.4 curve of Fig. 2 (a) is presented in Fig. 3. We note that the "normal"-superradiant QPT by increasing λ is of first order, while the superradiant-"normal" QPT for growing λ is of second order. From the figure it is not clear if the second order QPT is really a QPT or a cross-over. Figure 4, showing a slice (ω = 0.8) of Fig. 3 around one singularity λ s ≈ 1.57, confirms that it is indeed a second order QPT. Thus, contrary to the regular DQPT of normal-superradiance, the QPTs can be either of first or second order nature in this extended DM. Figure 3 does not only reveal the type of PTs of Fig. 2, but also that there is a number of first order QPTs between various superradiant phases.  One example of the T − ω phase diagram is presented in Fig.5 (a). As for the λ−ω diagram, both first and second order QPT exist. Figure 5 (b), displaying the number of bound states, confirms that the sudden changes (first order QPTs) are due to additional bound states. The T −ω phase diagrams again display several different superradiant phases. In this case, however, there are both first and second order PTs between the superradiant phases. Interestingly, our numeric analysis indicates that the T − ω phase diagrams seem to be fairly independent of Ω. Contrary to the regular DM, in the ultracold regime atomic motion plays an important role for the system characteristics. The confining potential is determined for a given field intensity y, and consequently, the atomic density ρ at (x) per particle depends as well on y. In other words, apart from changes in the field intensity and in the atomic inversion, the phase transition will as wll be manifested in the atomic density. A natural consequence of the coupled system considered here is that the motional state of the atoms are entangled with the internal state. Thus, the phases (normal and superradiant) are intrinsically different from the ones of the regular DM. The regular DQPT derives from a competition between the free field energy and the interaction atom-field energy. While the free field is minimized by vacuum, the interaction energy decreases with an increasing field in-tensity. In the present model does the kinetic energy contribute to the total energy. It is thus an interplay between three terms; free field, atom-field interaction (containing the bare internal atomic energies), and motional energies. This is not so evident from the free energy per particle (19), where the kinetic energy is hidden in the second term. To illuminate the importance of the atomic motion we study the atomic inversion W defined as the probability for a single atom to be in its excited state |e minus the probability for it to be in the ground state |g . Once the adiabatic approximation has been imposed we are left with a single internal state; the lower adiabatic state |− ad = sin θ|e + cos θ|g . In particular, the angle tan 2θ = 2g(x) √ ρ 0 α/Ω √ N depends on the spatial coordinate and the inversion becomes (21) The above equation clarifies that the atomic motion enters the problem in a non-trivial way. It also shows how the internal atomic properties are taken care of even though in the adiabatic approximation the system properties derives from a single internal state |− ad . We have numerically verified the appearance of the PT in terms sudden changes in ρ at (x) (first order PT) or ∂ρ at (x)/∂x (second order PT).
In deriving the phase diagrams, tight confinement of the atoms in two directions has been assumed. If only the longitudinal motion is frozen out, one regains an effective two dimensional problem whose eigenvalues are obtained from the Schrödinger equation with potential V (x, y) = −λ exp − x 2 +y 2 ∆ 2 x . As in the one dimensional situation studied in this section, V (x, y) possesses a finite number of bound states and one would expect very similar phase diagrams for this two dimensional case as for the one dimensional model.

C. Validity of the adiabatic approximation
We conclude this section by analyzing the adiabatic approximation. By a simple rotation, the amplitudes α appearing in the single atom Hamiltonian (6) can be taken real. For real α, the two last terms of h(α) is readily diagonalized by the unitary transformation [26] where the angle θ was given right above Eq. (21). Momentum transforms as UpU † =p − (σ + +σ − ) ∂θ, ∂θ ≡ ∂θ/∂x. Due to the spatial dependence of θ = θ(x), the transformed Hamiltonian is non-diagonal;h(α) ≡ U h(α)U † = h ad (|α|) + h cor (α). Here, h ad (|α|) is the adiabatic Hamiltonian (7) and h cor (α) contains the nonadiabatic corrections. Explicitly one finds [26] h cor (α) = 1 2 The kinetic energy is smaller or of the same order as max x V + ad (x, |α| 2 ), which provides a measure of h cor (α) in comparison to the adiabatic Hamiltonian h ad (|α|). For typical parameters, ω = Ω = 1, ∆ x = 3, and y = 2, the terms of h cor (α) are at least one order of magnitude smaller than the terms of h ad (|α|), which justifies the use of the adiabatic approximation.

IV. LONGITUDINAL THERMODYNAMICS
In the previous section we studied how the DQPT was modified due to motion of the atoms in a finite potential well, assuming the atomic motion to be frozen out in the longitudinal and one transverse direction. Here we instead assume the atoms to move freely along the center axis of the Fabry-Perot cavity while tightly bound in the transverse directions.
The corresponding single atom Hamiltonian (6) reads where µ is the scaled photon wave number which will be set to unity hereafter, µ = 1. This Hamiltonian, with the field still quantized, has been considered in several papers, see for example [27]. Normally L ≫ 2π, where L is the cavity length, so neglecting boundary effects is not a crude approximation [27]. The Hamiltonian is of the form of a generalized Mathiue equation, and hence the spectrum E ν (k) is described by a band index ν and a quasi momentum k extending over the first Brillouin zone. Due to the internal two-level structure of the atom, the Brillouin zone is twice the size of what is imposed by the periodicity of the mode [27]. Clearly, E ν (k) depends on the field amplitude α. The corresponding eigenfunctions are written Φ k,ν (x) = φ (e) k,ν (x)|e + φ (g) k,ν |g . For a constant coupling g(x) = λ 0 , these Bloch functions are simple plane waves giving a constant energy shift independent of system parameters such as λ, ω, and Ω. For a standing wave mode coupling on the other hand, the Bloch functions cannot be decoupled from the internal atomic states, and consequently the atomic motion will affect the structure of the phase diagrams as will be demonstrated below.
In the previous section we assumed the adiabatic regime and diagonalized the Hamiltonian in its internal degrees of freedom. In this case, the adiabatic potentials V ± ad (x) cross and adiabaticity breaks down in the range where the QPTs occur. Fortunately, the Hamiltonian is easily diagonalized numerically by truncation the dimension of the Hamiltonian matrix. We present, however, asymptotic analytical results in the Appendix which relies on the adiabatic approximation. These analytical results enable us to extract the limiting situation of large field amplitudes. Furthermore, as a numerical diagonalization directly renders several of the Bloch bands we do not restrict the analysis to just the lowest one. However, it turns out that for most of the presented examples only the lowest band contribute to the dynamics due to the low temperatures considered. Exceptions are in the plots of the critical temperature where we indeed go to rather high temperatures and the excited bands become important.
The partition function is written like in the previous section as where C 2 is a constant and the free energy per particle with Here, as above, y = |α| 2 /N represent the scaled field intensity. Shown in the Appendix, the second part of f 2 (y) scales asymptotically as ∼ √ y for large y. Thus, a maximum of the free energy can only be obtained for finite or zero field intensities y.
As in Sec. III, we derive the critical atom-field coupling λ c and temperature T c . In Fig. 6 we show the results of how the critical coupling depends on ω (a) and Ω (b) for different temperatures. The critical temperature as function of ω and Ω is displayed in Fig. 7, where the inserted numbers indicate the values of the coupling λ. In both cases, the critical quantities show clear differences compared to the ones of the regular DM, (2). The critical coupling scales as λ c ∼ √ ω for fixed Ω just like in the regular DM. However, the Ω-dependence is not possessing the same structure as (2), and especially for small values on Ω the critical coupling λ c is non-zero. This was also found in a nearest-neighbor coupling model studied in [18].
In the regular DM, for a fixed λ and Ω the critical temperature diverges for small ω and goes to zero for large ω. Here we note that for high temperatures the regular behavior is regained, while for low temperatures, and in particular zero temperature, a QPT takes place for finite ω. Fixing ω instead and vary Ω we get even more surprising results. In the regular DM, there is an upper temperature for which the QPT is lost and at zero temperature a QPT occurs for finite Ω. In our model, a similar phase-diagram is obtained for a range of parameters ω and λ, but there also exist parameter regimes where no QPT occurs for zero temperature.
Like in the previous section the nature of the PTs is studied by introducing the scaled field intensity I maximizing f 2 (y). It is found that the QPT's are of second order character in all cases. The numbers to each curve display the dimensionless temperature T . Note in particular that for Ω → 0, the critical coupling λc = 0.

V. CONCLUSIONS
In this work we have studied a new regime in the DM. The atoms are assumed trapped by the cavity field itself and ultracold such that their center-of-mass kinetic energy is of the order of the atom-field interaction. This calls for a full quantum mechanical treatment of the atomic motion and at the same time take into account for spatial mode variations. The analysis is motivated from our earlier findings, where we demonstrated that atomic motion greatly affects the system dynamics in many-body cavity QED systems [4]. Expectedly, we have shown that this is also true for the DM. The analysis is restricted to considering one dimensional problems, and both the case of a Gaussian and a standing wave mode profile were treated. However, we motivated that similar phenomena are expected also for higher dimensional situations. In particular, we made evident that the varying number of bound states in the potential formed by The dimensionless critical temperature Tc as function of the system parameters. The inserted numbers give the value of λ. We note that a QPT is possible at zero temperature and finite ω, and moreover, that for zero Ω the QPT may vanish for small temperatures. In (a) Ω = 1 while in (b) ω = 1.
a Gaussian mode profile induces novel first order QPTs. Additionally, we found that for certain couplings λ s , the system is superradiant for any field frequency ω. Moreover, great differences with the regular DM was also encountered for a standing wave mode profile.

APPENDIX A: TIGHT BINDING APPROACH
Here, we analytically consider the large field asymptotic expressions for the energy per particle f 2 (y) in the case of a standing wave mode profile. In this regime one may utilize the tight binding approximation [30] to derive the spectrum. We further assume the adiabatic approximation to be valid and that we can restrict the analysis to the lowest energy band. The adiabatic potentials are given by V ± ad (x) = ±h Ω 2 4 + λ 2 cos 2 (x)y, where y = |α| 2 /N is as before the scaled field intensity. A convenient base for writing down the periodic Hamiltonian in matrix form is to use the Wannier states w ± j (x) = x|j w,± [30]. The function w ± j (x) is localized in the jth "well" of the potentials V ± ad (x). By the tight binding approximation we assume w,± i|h ± ad |j w,± = 0 unless i = j or i = j ± 1. Within the validity regime of this approximation we may as well replace the Wannier functions by Gaussian functions [4]. The widths of the Gaussians are given by approximating the potential wells by harmonic oscillators, giving where and x j is the position of the jth potential well. To avoid un-physical contributions from the non-orthogonality of the Gaussians we impose dx w G (x−x j )w G (x−x i ) = δ ij . We further introduce the matrix elements where we only consider i = 0, 1. We note that the Wannier functions are directly related to the depth of the corresponding potential and therefore their width σ will also depend on y. This explains the field intensity dependence of E i (y). Another important observation is that w ± G (x−x j ) are localized where |V ± ad (x)| are close either to its maximum or its minimum, resulting in different coupling elements J ± i (y), indicated by the ±-superscript. In this notation we get the lowest band tight binding energy E ± 1 (k) = E 0 (y) + J ± 0 (y) + E 1 (y) + J ± 1 (y) 2 cos(k).
The part in front of the cosine function is strictly negative resulting in that the ground state energy is given by k = 0. The kinetic energy integrals of (A4) are readily solvable, and one finds E 0 (y) = 1 4σ 2 , E 1 (y) = − 1 8σ 4 exp − π 2 4σ 2 2σ 2 + π 2 . (A6) The potential integrals of (A4) are not analytically solvable for the given potentials (A1). Instead we make the same kind of approximation as in section III V ± ad (x) ≈ ±A ± B cos 2 (x), (A7) and identify Within this regime we find We emphasize that the width σ 2 depends on the field intensity y; The applied approximations are only reliable for z < 1 [4], and it turns out that the QPTs occur beyond these approximations. Nonetheless, we may find the asymptotics for the free energy f 2 (y). In the large y limit we find that ln [g 2 (y)] ∼ √ y. Consequently, the field intensity I will always be finite, regardless of parameter choices. We have verified numerically the y square-root dependence of ln [g 2 (y)] for large intensities.