Quantum phase transitions of atom-molecule Bose mixtures in a double-well potential

arXiv:1406.6833v2 [quant-ph] 20 Oct 2014 Quantum phase transitions of atom-molecule Bose mixtures in a double-well potential A. Rela˜no,1J. Dukelsky,2P. P´erez-Fern´andez,3and J.M. Arias4 1Departamento de F´ısica Aplicada I and GISC, Universidad Complutense de Madrid 2Instituto de Estructura de la Materia, CSIC, Serrano 123, 28006 Madrid 3Departamento de F´ısica Aplicada III, Universidad de Sevilla 4Departamento de F´ısica At´omica, Molecular y Nuclear, Universidad de Sevilla and Unidad Asociada to IEM (CSIC) Madrid (Dated: October 21, 2014) Abstract The ground state and spectral properties of Bose gases in double-well potentials are studied in two different scenarios: i) an interacting atomic Bose gas, and ii) a mixture of an atomic gas interacting with diatomic molecules. A ground state second-order quantum phase transition (QPT) is observed in both scenarios. For large attractive values of the atom-atom interaction, the ground-state is degenerate. For repulsive and small attractive interaction, the ground-state is not degenerate and is well approximated by a boson coherent state. Both systems depict an excited state quantum phase transition (ESQPT). In both cases, a critical energy separates a region in which all the energy levels are degenerate in pairs, from another region in which there are no degeneracies. For the atomic system, the critical point displays a singularity in the density of states, whereas this behavior is largely smoothed for the mixed atom-molecule system. 1 I. INTRODUCTION Bose-Einstein condensation in dilute atomic gases provides a valuable tool to study fundamental problems of quantum systems at a macroscopic scale. Among such phenomena, the quantum transition from Josephson oscillations to self-trapping is one of the most interesting. It comprises the transition from a regime in which two spatially separated Bose-Einstein condensates (BEC) oscillate, to another one in which the macroscopic wave function remains trapped in one of the two trap sides. This transition has been experimentally observed in [1] by means of a condensate of 87Rb atoms in an optical trap. It has also been theoretically studied in Refs. [2–4]. Using a two-site Bose-Hubbard Hamiltonian as a description of two BECs that can tunnel between the two wells, it has been shown that the quantum phase transition (QPT) develops at a critical value of the atom-atom interaction. The imbalance, that is, the difference of population between the two wells, acts as an order parameter for the transition, changing from zero in the Josephson oscillation regime to a non-zero value in the self-trapping regime. This is an example of QPT which has been studied in other systems. In [5], the same kind of transition was reported in a non-linear version of the same system, and it was related to a classical bifurcation in the Hamiltonian. In [6], a similar bifurcation was linked to a transition from zero to non-zero values in the imbalance. In this work, we expand these studies by focusing not only in the ground state QPT, but also looking at quantum critical phenomena in the energy spectrum. In particular, we will describe the so called excited-state quantum phase transitions (ESQPTs). These occur when some kind of non-analytic behavior [7] arise in an excited region of the energy spectrum, normally a singularity in the density of states or in the flow of the energy spectrum through the critical energy line. In contrast to usual quantum phase transitions [8], they are not characterized by a critical value of the control parameter λc, but by a critical energy Ec above the ground state, at which the non-analyticity takes place. Usually, ESQPTs are linked to QPTs; below the critical parameter λcno transition appears, while above λcthe non-analytic behavior characteristic of the QPT propagates to a certain critical excitation energy Ec. However, in some cases an ESQPT occurs without the corresponding QPT. As a paradigmatic example, we mention the realization of the Lipkin-Meshkov-Glick model studied in [6]. 2 ESQPTs entail dramatic dynamic consequences, like the enhancement of the decoherence if the ESQPT takes place in the environment with which a quantum system interacts [9], the onset of chaos [10] or the emergence of symmetry-breaking equilibrium states after a quench [11]. They appear in models pertaining to different branches of Physics, like the interacting boson model (IBM) [12, 13], Lipkin-Meshkov-Glick model [14], vibron model [15], Dicke and Jaynes-Cummings models [10, 16, 17], kicked top [18] and microwave Dirac billiards [19]. It is worth mentioning that, despite the fact that the excitation energy is linked to the temperature of an isolated system, ESQPTs are qualitatively different from thermal phase transitions. The existence of a critical energy does not necessarily entail a corresponding critical temperature. For example, the Dicke model restricted to the maximum angular momentum sector displays an ESQPT, but has no thermal phase transition [20]. In order to investigate these issues, we treat two related systems. The first consists of an atomic Bose gas trapped in a two-well potential, and the second incorporates the interaction of the atoms with a diatomic molecule. Mixed condensates of atoms and molecules have been experimentally realized by means of the manipulation of a Feshbach resonance [21–23]. Mixtures of atoms and molecules in double-well potentials have been studied previously in [24]. In [25] a similar system was used to model the atom-atom interaction by means of a Feshbach resonance. The organization of the paper is as follows. We introduce the two models in Sect. II. In Sect. III and Sect. IV the ground state quantum phase transitions are investigated numerically as well as analytically for both models; first, within a mean field approach, for which we later incorporate fluctuations. Sect. V is devoted to analyze ESQPTs in both models. Finally, in Sect. VI the main conclusions of this work are summarized. II. THE MODELS We shall study a bosonic atom gas confined in a double-well potential as a model of a bosonic Josephson junction [26]. It has been recently shown that a restriction to the lowest two modes is a valid approximation for a wide range of the bosonic interaction [4]. In this limit, the model reduces to the two-site Bose-Hubbard Hamiltonian. Each site represents the lowest mode of the right and left wells, respectively, H1=−Ja† LaR+a† RaL+U Na† La† LaLaL+a† Ra† RaRaR,(1) 3 where the a† kand ak(k=L, R for left and right wells) operators are the atom creation and annihilation operators in each well and obey the usual boson commutation rules, Nis the total number of atoms, Jis the hopping strength, and Uis atomic interaction. Note that the scaling factor 1/N in front of the atomic interaction is needed to ensure that the energy per particle E/N is well-defined in the thermodynamic limit. The static and dynamic properties of the Hamiltonian (1) have been extensively studied in recent years [27–30]. Moreover, it has been shown that (1) is a quantum integrable model and a particular limit of Richardson- Gaudin bosonic families [31–33]. In addition to (1), we will consider here an extension of the two-site Hubbard model by the explicit inclusion of an interacting diatomic molecule. The resulting two-channel Hamiltonian is: H2=−Ja† LaR+a† RaL+U Na† La† LaLaL+a† Ra† RaRaR +ωb†b−g √2Nhb†(aLaL+aRaR) + a† La† L+a† Ra† Rbi,(2) where Nis the total number of atoms plus twice the number of molecules, N= 2b†b+ a† LaL+a† RaR, and b†(b) creates (annihilates) a molecule. The operators band b†follow the usual bosonic commutation rules. Because of the shape of the atom-molecule interaction (two atoms are destroyed to create a molecule, and viceversa) they represent a dimeric entity, giving rise to the previously written conservation of the total number of particles N. The frequency ωrepresents the self-energy of the molecule. When it goes to zero, the Hamiltonian (2) describes a Feshbach resonance interacting with the atomic cloud. In what follows we will use the hopping term Jas the energy unit. This Hamiltonian is in general non-integrable. However, as a function of the interaction parameters ω,U, and git has two integrable limits which will allow us to explore the transition from regularity to chaos. As in the previous case, the scaling factors 1/N and 1/√Nin front of the atom-atom and the atom-molecule interaction ensure that the energy per particle is well-defined in the thermodynamic limit. For both Hamiltonians, (1) and (2), Nis a conserved quantity. They also conserve a kind of parity, related to the interchange between atoms in left and right wells. In connection with this symmetry, it is worth mentioning that the population imbalance operator that counts the difference between the number of atoms in the left and right wells ˆ I=a† RaR−a† LaL,(3) 4 does not commute with the parity operator. Therefore, its expectation value is zero in any common eigenstate of the Hamiltonian and the parity operator. III. PHASE DIAGRAM In this section we analyze and compare the ground-state phase diagram of both models, the double well with and without the interaction with the diatomic molecule. In order to carry out this comparison, we select the atomic interaction Uthat is common to both models as the control parameter, and keep fixed ωand g. As it will be seen below, the ground state of one of the phases is doubly degenerate due to spontaneous symmetry-breaking. Thus, the population imbalance is a good candidate of order parameter for this transition. In fact, it has been used in [4] to characterize the macroscopic self-trapping in the two-site Bose- Hubbard model. However, it is not possible to rely directly on the population imbalance to study numerically this transition, because both Hamiltonians, (1) and (2), have a well defined parity and, therefore, the imbalance expectation value is zero in any of their eigenstates. One possible solution is to introduce a small symmetry-breaking term to break the degeneracies, as it has been done in [4]. In this work we do not follow this way; we will use, instead, a new order parameter derived from a quantum version of the Fisher information, defined as [34] FQF I =(bnR−bnL)2−(hbnR−bnLi)2.(4) For eigenstates with well-defined parities, this expression reduces to: FQF I =(bnR−bnL)2.(5) Therefore, to obtain a quantity proportional to the number of left or right-well atoms, we will use pFQF I as the order parameter of the transition. As it will be seen below, this magnitude behaves in the same way as the population imbalance in the exact solution of both Hamiltonians. A. Mean-field approximation The mean-field or semi-classical approximation is expected to provide accurate results for this kind of two-level systems in the large Nlimit. In fact, it has been shown that the 5 ground state energy in this approximation is exact to leading order in N[35]. The great advantage of this approach is that it can be solved analytically. We use boson coherent states as variational wavefunctions to model the ground state of both Hamiltonians. The corresponding wavefunctions are: |Ψ1i=e√N(γRa† R+γLa† L)|0i,(6) for the case without molecule, and |Ψ2i=e√N(β 2b†+γRa† R+γLa† L)|0i,(7) for the case with molecule. They are eigenstates of the annihilation operators aL,aR(both) and b(only the last), with eigenvalues √NγR,√NγL, and √Nβ/2, respectively; hence, they fulfill the mathematical requirement for coherent states. The former could be obtained from the second one as the limit in which the atom-molecule coupling parameters gand the corresponding amplitude in the coherent state (β) go to zero, and the self-energy of the molecule ωgoes to infinity: g→0, ω→ ∞and β→0. Therefore we will mostly concentrate on the second model, and treat the other one as a particular case in the appropriate limit. It is important to note that both coherent states break the parity symmetry and will give a finite population imbalance different from zero in the symmetry-broken phase. The energy per particle using the coherent state |Ψ2iis E N=hΨ2|H2|Ψ2i NhΨ2|Ψ2i=−2JγRγL+ω 2β2−gβ(γ2 R+γ2 L) + U(γ4 R+γ4 L),(8) where the variational parameters are assumed to be real, provided g > 0. By minimizing this energy surface for the corresponding set of parameters J,ω,gand U, optimal variational parameters β, γRand γLare obtained, providing an approximation to the exact ground state wavefunction and energy. This method is only exact, as mentioned above, for the leading order in N[35]. Because of the conservation of the total number of particles in the system, there is a constraint for the variational parameters: β2+γ2 R+γ2 L= 1. This condition is included in the minimization by using a Lagrange multiplier which introduces a new parameter, λ. When the energy surface, with the constraint of the particle number conservation, is minimized, the optimal values of the variational parameters as a function of the control 6 parameters are determined from the set of equations, β2+γ2 R+γ2 L= 1,(9) −g(γ2 R+γ2 L) + ωβ −2λβ = 0,(10) gβγR−2Uγ3 R+JγL+λγR= 0,(11) JγR+gβγL−2Uγ3 L+λγL= 0.(12) This set of four non linear equations has two classes of solutions: a) the symmetric, with γR=γL=γ; and the non-symmetric, b) γR6=γL. The former corresponds to the symmetric phase; the latter, to the non-symmetric phase. The solutions in the symmetric phase are β=                                −g 2U+(1−i√3)A 63 √4U(B+√B2+4A3)1/3−(1+i√3)(B+√B2+4A3)1/3 12 3 √2Uif U∈(Uc,0) −2J−ω+√12g2+(2J+ω)2 6gif U= 0 −g 2U−A 33 √4U(B+√B2+4A3)1/3+(B+√B2+4A3)1/3 63 √2Uif U∈(0,∞) (13) γ=r1−β2 2,(14) with A=−9g2+ 6U(2J−2U+ω),(15) B= 54g(2JU −g2+Uω).(16) The range of variation of the Hamiltonian parameters for this solution is: g∈R− {0}, ω and J∈R.Ucis the critical value for which the system undergoes the QPT. For the non-symmetric phase, U∈(−∞, Uc), the solutions are 7 β=−g 4U+(1 −i√3)A′ 12 3 √2UB′+p(B′)2+ 4(A′)31/3− (1 + i√3) B′+p(B′)2+ 4(A′)31/3 24 3 √2U,(17) γR=s1−β2 2+p−J2+U2(β2−1)2 2U,(18) γL=−J 2UγR .(19) with A′= 12U(ω−4U)−9g2,(20) B′= 54g(2Uω −g2).(21) For this solution the Hamiltonian parameters are restricted to g∈R−{0}, ω ∈Rand J∈ R−{0}. The critical value of U(Uc) can be computed using equations (13) and (17). The variational parameter βhas to be a continuous function of Uto ensure the continuity of the ground state energy. By equating both expressions of βin the symmetric, (13), and nonsymmetric, (17), phases we can determine numerically the critical value Uc. For instance, for a system defined by J= 1, ω= 5 and g= 5, the value obtained for Ucis Uc=−1.14018. A similar analysis can be performed for the limit without coupling to the molecule obtaining the following expressions for the variational parameters γR=                q1 2+ Sign(U)√U2−1 2Uif U < Uc 1 √2if U≥Uc (22) γL=                  −J 2U 1 r1 2+Sign(U)√U2−1 2U if U < Uc 1 √2if U≥Uc (23) 8 0.2 0.4 0.6 0.8 1.0 Variational parameters γR γL -4 -2 0 2 4 Control Parameter U 0.0 0.2 0.4 0.6 0.8 <nR-nL>/N √ FQFI N=200 √ FQFI N=100 (a) (b) 0.0 0.2 0.4 0.6 0.8 1.0 Variational parameters γR γL -4 -2 0 2 4 Control Parameter U 0.0 0.2 0.4 0.6 0.8 <nR-nL>/N √ FQFI N=3200 (c) (d) FIG. 1: (Color online) Panels (a) and (c) show the variational parameters γRand γL, as a function of the control parameter U. Panels (b) and (d) show the normalised expectation value of the invalance operator b Iin solid line, and the square root of the Fisher information, obtained by numerical diagonalization, in open symbols. Panels (a) and (b) are for the mixed atom-molecule system (two numerical calculations for different system sizes are shown in (b), (black) circles correspond to N= 100, and (orange) squares, to N= 200) whereas panels (c) and (d) present the system without molecule. In this case, the critical value for the coupling constant is Uc=−1. Note that Eqs. (9)-(12) are invariant under the transformation γL→ −γL,γR→ −γR, which implies that there always exist two variational solutions. In the symmetric phase, both of them give rise to b I= 0, corresponding to the same physical state. On the contrary, in the non-symmetric phase both have non-zero imbalance, b I1=−b I2, implying that the ground state is degenerate. So, in this case any linear combination of these two solutions is also a ground-state of the system, and thus it is possible to build a parity-projected solution, with b I= 0. As a consequence, the critical coupling Ucseparates two different regions: the symmetric one, in which the populations of the two wells are always equal, and the non-symmetric one, in which one can find symmetry-breaking ground states. As the coherent states used in this section break the symmetry for γL6=γR, they provide an accurate description of the transition in terms of the imbalance. 9 -4 -2 0 2 4 U 0 0.05 0.1 Excitation energy (arb. units) (a) -4 -2 0 2 4 U 0 0.05 0.1 Excitation energy (arb. units) (b) FIG. 5: (Color online) First few excited state energy gaps as a function of the control parameter U. Symbols are the exact results for N= 200, while lines provide the beyond mean field approximation. The rest of the parameters are fixed to J= 1, ω= 5 and g= 5 in panel (a), and J= 1 in panel (b). The vertical dot line marks the separation between the two phases. V. EXCITED-STATE QUANTUM PHASE TRANSITIONS In some collective many-body quantum systems, QPTs are accompanied by ESQPTs, giving rise to a critical energy Ecfor certain values of the system control parameter [9–17]. Traditionally, ESQPTs have been linked to singularities in the density of states or in one of its derivatives, depending on the number of system’s degrees of freedom in the semiclassical limit [38]. It has been shown that for a system with just one effective degree of freedom, as the Tavis-Cummings model, a λ-singularity in the density of states characterizes the ESQPT. For the Dicke model, which has two effective degrees of freedom, the singularity making the ESQPT appears in the first derivative of the density of states [17]. This fact entails an increasing difficulty for detecting such a critical behavior in complex systems, with many effective degrees of freedom. In this work, we rely on a different way to identify ESQPTs. In general, the critical energy separates two different regions of the spectrum which share some of the properties of the ground state at one and at either side of the QPT. For instance, in Dicke and Lipkin models the critical energy separates one region in which all eigenstates are doubly-degenerate from another in which there are no degeneracies and all the levels have a well defined parity [11]. So, in order to characterize the ESQPTs in our system, we study the distance between 16 -10 -5 0 5 10 -10 -5 0 5 10 E/N U (a) -10 -5 0 5 10 -10 -5 0 5 10 E/N U (b) FIG. 6: (Color online) Complete phase diagram. Panel (a) depicts the case with molecule, obtained with N= 320 particles, and J= 1, ω= 5 and g= 5. Panel (b) depicts the case without molecule, with N= 1000 particles and J= 1. In both cases, dark filled region (blue online) indicates the symmetry-breaking phase, in which the eigenstates are doubly-degenerate, and light filled region (yellow online) the normal phase, in which there are no degeneracies and the parity is a good quantum number for all the eigenstates. Vertical dotted lines enclose the region in which there is no ESQPT. consecutive energy levels with opposite parity. We follow the same criterion as in Fig. 2. Figure 6 shows the energy per particle versus the strength of the Uparameter for the two systems. Panel (a) corresponds to the case with a diatomic molecule and a total of N= 320 particles, with J= 1, ω= 5 and g= 5; while panel (b) corresponds to the case without molecule, with N= 1000 particles and J= 1. The lowest line is the ground-state energy, while the uppermost line gives the limit of the calculated spectra for the system size used (highest calculated excited state). There are many excited states in between both lines. In both panels two vertical dotted lines are plotted: the left line gives the ground state QPT, Uc, while the right one marks an equivalent singularity for the highest excited state, e Uc. The energy region with no degeneracies is called normal (light filled region, yellow online), while the region with double degeneracies is a reflection of the appearance of a broken-symmetry (dark filled region, blue online). Our calculations clearly identify an ESQPT, that is marked by the line in between light and dark regions. It separates normal (non-degenerate) from broken-symmetry (degenerate) phases. Figure 6 shows that both systems behave in the same qualitative way. For U < Uc, the lower part of the spectrum is degenerate, whereas the upper part is not degenerate. For large values of the control parameter, U > e Uc, the 17 spectrum behaves in the opposite way: the upper part is degenerate, whereas the lower part is not. In between, Uc≤U≤e Uc, all the energy levels are not degenerate. This entails that, for attractive interactions U < 0 there is a QPT at a critical coupling Uc, and a critical energy Ecfor U < Uc. On the contrary, if the interaction is repulsive U > 0, an analogue of a QPT takes place on the highest-excited level at a critical coupling e Uc, and also a critical energy appears for U > e Uc. No ESQPT is observed in the region Uc≤U≤e Uc. The critical line can be easily estimated from the numerics in the case without the molecule and gives Ec=     U 2−1, U < Uc=−1; U 2+ 1, U > e Uc= 1. (35) This result entails that the spectrum is turned upside down when changing the interaction from attractive to repulsive, being the critical points symmetric with respect to the free case, U= 0. The same analysis is not that clean for the system with a molecule. The precise values for e Ucand Ec(U) depend on the molecular parameters ωand g. For ω6= 0 and g6= 0 the symmetry around U= 0 is broken. Furthermore, due to the smaller values of the total number of particles Naccessible to the our current computational power, it has not been possible to obtain a precise estimate of the behavior in the thermodynamic limit. For g=ω= 5 a numerical fit to the case with N= 320 particles gives rise to Ec=     0.46U−1.89, U < Uc∼ −1.3 0.42U+ 2.55, U > e Uc∼1.7 (36) where both critical values, Ucand e Uc, correspond to the finite-size precursors, which are still far from the true critical points in the thermodynamic limit. It is possible to study ESQPT for fixed values of the system parameters by looking at the density of states ρ(E) as a function of the excitation energy. A singular behavior in ρ(E) signals an ESQPT at the critical energy Ec. It has been recently shown that ρ(E) presents a logarithmic singularity at E=Ecfor models with just one semiclassical degree of freedom. On the contrary, the density of states for models with more than one semiclassical degree of freedom is not singular. In these cases, the logarithmic singularity occurs in its derivative ρ′(E) [38]. This feature has been tested in the integrable Tavis-Cummings model with one semiclassical degree of freedom, and in the non-integrable Dicke model with two 18 semiclassical degrees of freedom [17]. Following a similar procedure, we fix the value of the interaction U= 9 and calculate the density of states as a function of the excitation energy for the two systems. From the results of Eqs. (35) and (36) we estimate the critical energy to be around E∼6.3 in the case with molecule, and E∼5.5 in the case without molecule. 0 0.05 0.1 0.15 0.2 0.25 0 2 4 6 8 10 ρ(E) E/N (a) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 4 5 6 7 8 9 ρ(E) E/N (b) FIG. 7: Density of states ρ(E): Panel (a) for N= 320, U= 9 and a molecule at resonance. Panel (b) for N= 2000 without molecule. Figure 7 shows the density of states ρ(E) versus the excitation energy. Panel (a) corresponds to the case of a molecule at resonance and N= 320 particles. It shows no neat trace of a singular behavior around the estimated critical energy Ec∼6.3. This is however compatible with a logarithmic singularity in the derivative of the density of states, which is very difficult to see in systems of this size. Panel (b) corresponds to N= 2000 without molecule. It clearly shows a precursor of the logarithmic singularity at Ec= 5.5. The results obtained with these two systems are therefore, compatible with those derived in Ref. [17] for the Tavis-Cummings and Dicke models. In order to get a deeper insight into the properties of the ESQPTs, we plot in Fig. 8 the ratio of degenerate energy levels as a function of the scaled energy E/N. The calculation has been done in the following way: a) we choose a set of 100 energy levels, starting from the ground state; b) we count the number of pairs whose relative energy difference is less than 10−6; c) we normalize this value to 100; d) we associate this value to an energy equal to the mean energy of the interval; e) we repeat the calculation starting from the 11th energy level; f) we proceed in the same way up to covering the whole spectrum. In the same figure we plot the results for the system with (panel (a)) and without the molecule (panel (b)). 19 0 20 40 60 80 100 0 2 4 6 8 10 % Degenerated E/N (a) 0 20 40 60 80 100 6.4 6.6 6.8 7 0 20 40 60 80 100 3 4 5 6 7 8 9 % Degenerated E/N (b) 0 20 40 60 80 100 5.3 5.4 5.5 5.6 FIG. 8: Ratio of degenerate energy levels as a function of the energy: panel (a) for N= 320, U= 9 and a molecule at resonance. Panel (b) for N= 2000 without molecule. The case without molecule shows an abrupt transition from non-degenerate to degenerate energy levels at the expected value for the critical energy Ec∼5.5. The case with molecule is not as sharp, though the critical behavior is still clearly observed. For this case, a small number of degeneracies appear distributed randomly in the lower part of the spectrum. We attribute this behavior to the spectral fluctuations of the energy levels, and relate it to the crossings and avoided crossings shown in Fig. 5 in the symmetric phase. In spite of this behavior at low energies, we observe an abrupt change from non-degenerate to degenerate in the spectrum at an energy Ec∼6.6. The 5% discrepancy with our estimate is attributed to finite-size effects. Figure 9 shows the behavior of the order parameter pFQF I for both systems, with molecule (panel (a)) and without molecule (panel (b)). It has been obtained by calculating the square root of the Fisher information for all the eigenstates of the systems. The results are similar to those obtained for the density of states. For the case without the molecule the quantum Fisher information displays a sharp singularity at the estimated critical energy Ec∼5.5. On the contrary, there is no such a signature in the case with molecule. These results are compatible with the fact that ESQPTs are softer in systems with more than one semiclassical degree of freedom, as it is pointed out in [16, 38]. In the case with the molecule, large fluctuations between different eigenstates blur the expected singular behavior of the order parameter at the estimated critical energy. On the other hand, in the limit without molecule, which has just one semiclassical degree of freedom, a cusp singularity at 20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 1 2 3 4 5 6 7 8 9 10 √FQ F I E/N (a) 6.1 6.3 6.5 6.7 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 3 4 5 6 7 8 9 10 √FQ F I E/N (b) 0.1 0.14 0.18 5.3 5.4 5.5 5.6 FIG. 9: Fisher information as a function of the energy: Panel (a) for N= 320, U= 9 and a molecule at resonance. Panel (b) for N= 2000 without molecule. -0.1 0 0.1 0.2 0.3 0.4 0.5 0 1 2 3 4 5 6 7 8 9 η E/N 0 0.1 0.2 0.3 6 6.3 6.6 6.9 FIG. 10: Spectral statistics for U= 9 and N= 360, even-parity states, as a function of the scaled energy E/N. the critical energy is clearly seen. Finally, it is interesting to analyze the degree of chaos in the spectra as a function of the energy. In Ref. [10] it was conjectured a connection between the ESQPT and the onset of chaos in the Dicke model. Here we analyze the nearest-neighbor spacing distribution as the simplest measure of quantum chaos. While integrable systems follow a Poisson distribution PP(s) = exp(−s), non-integrable and fully chaotic systems are well described by the Wigner surmise PW(s) = π 2sexp −π 4s2. Given a set of spacings si= (Ei+1 −Ei)/hsi, a natural 21 way to measure the degree of chaos is by means of the ratio η=vars−varW varP−varW ,(37) where “var” means the variance of the distribution. The parameter ηvaries from η= 0 for fully chaotic systems, to η= 1, for integrable systems. In Fig. 10 we show the parameter ηas a function of the scaled energy E/N. The calculation has been done as follows: a) a set of 250 levels, starting from the ground state has been chosen, b) the set has been unfolded, discarding a 20% of higher and lower levels, c) the parameter ηhas been calculated, d) the value of ηhas been linked to the mean energy value of the set, e) the procedure has been repeated increasing the energy. It is important to note that all the eigenstates have a definite value of parity. As it can be seen in Fig. 10, a transition from a fully chaotic to an intermediate regime takes place at E∼Ec. Below this critical point, the system is more or less chaotic, with the exception of the region very close to the ground state, which is closer to integrability. This is fully compatible with the complex behavior shown in Fig. 5 for the lowest excited states in the symmetric phase. On the contrary, above the critical energy the system is far from chaos. The interplay between chaos and regularity as a function of energy is qualitatively similar to the Dicke model [10, 37]. For Dicke and the present model, chaos characterizes the normal phase, in which there are no degeneracies, whereas the symmetry-breaking phase is closer to integrability. Despite this fact, there seems to be no clear connection between the ESQPT and the onset of chaos [37]. Dynamics of systems with two semiclassical degrees of freedom suffering ESQPTs seems to be more chaotic in the phase in which there are no degeneracies in the spectrum. VI. CONCLUSIONS By means of analytic and numerical results, we have studied the QPT signaled by the population imbalance in two kinds of double-well BECs: one comprising bosonic atoms tunneling between the two wells, and the second including the interaction of the atoms with a diatomic molecule. Our main results are summarized as follows: Both systems display a second-order QPT at a certain critical value Ucof the atom-atom interaction. In both cases, for U < Ucthe ground state is degenerate due to the breaking of the left-right symmetry in the atomic system, giving rise to a finite population imbalance. 22 On the contrary, for U > Ucthe ground state is not degenerate and the imbalance is zero. This critical value of the control parameter was obtained by means of a mean-field calculation, and it was numerically tested by means of exact diagonalizations and finite-size scaling using the gap between the ground state and the first excited state as a signature of the transition. We have also related the population imbalance with the square root of the Fisher information, and tested the mean-field with exact diagonalizations. Beyond mean-field techniques that take into account the quantum fluctuations have been used to describe the low energy spectrum. A comparison with exact diagonalizations shows the goodness of these theoretical results. In particular, it is worth mentioning that the theory successfully accounts for the crossings between levels corresponding to opposite parities, which take place when the degree of freedom of the diatomic molecule is added. This formalism greatly improves the description of the critical point, showing that the secondorder QPT is strongly influenced by correlations and fluctuations that are beyond the meanfield approximation. Both systems also display ESQPTs at certain critical energies depending on the value of the interaction U. Contrary to the ground state QPT, the nature of these transitions differs in both models. The two-site Bose-Hubbard Hamiltonian is characterized by a λdivergence in the density of states. This behavior in the energy spectrum is comparable to the Lipkin and the Tavis-Cummings models, both having a single semiclassical degree of freedom. On the contrary, the density of states for the case in which the atoms interact with a diatomic molecule, shows no divergence. This behavior is qualitatively similar to the Dicke model, in which the divergence arises in the derivative of the density of states. Finally, we can infer a correlation between the ESQPT and the onset of chaos in the systems with atoms and molecules. A similar conclusion has been obtained for the Dicke model [10], though it seems that the transition between the integrable and the chaotic regime does not coincide with the critical energy at which the ESQPT takes place [37]. Taking into account the strong numerical limitations for the atom-molecule system, we conjecture that quantum chaos only appears in the non-degenerate region of the spectrum. Whether this is a generic mechanism, linking chaos and critical phenomena in the energy spectrum of systems with more than one semiclassical degree of freedom remains still an open question. 23 ACKNOWLEDGMENTS This work is has been partially supported by the Spanish MINECO grants FIS2012- 35316, FIS2012-34479, and FIS2011-28738-c02-01, by Junta de Andaluc´ıa under group number FQM-160 and Project P11-FQM-7632, and by the Consolider-Ingenio 2010 Programme CPAN (CSD2007-00042). [1] M. Albiez, R. Gati, J. F¨olling, S. Hunsmann, M. Cristiani, and M. K. Oberthaler, Phys. Rev. Lett. 95, 010402 (2005). [2] G.S. Paraoanu, S. Kohler, F. Sols, A. J. Leggett, J. Phys. B 34, 4689 (2001). [3] A. P. Tonel, J. Links, A. Foerster, J. Phys. A 38, 1235 (2005). [4] B. Juli´a-D´ıaz, D. Dagnino, M. Lewenstein, J. Martorell, and A. Polls, Phys. Rev. A 81, 023615 (2010). [5] V. S. Shchesnovich and V. V. Konotop, Phys. Rev. Lett. 102, 055702 (2009). [6] T. Zibold, E. Nicklas, C. Gross, and M. K. Oberthaler, Phys. Rev. Lett. 105, 204101 (2010). [7] M. A. Caprio, P. Cejnar, and F. Iachello, Ann. Phys. (NY) 323, 1106 (2008). [8] S. Sachdev, Quantum Phase Transitions (Cambridge University Press, Cambridge, 1999). [9] A. Rela˜no, J. M. Arias, J. Dukelsky, J. E. Garc´ıa-Ramos, and P. P´erez-Fern´andez, Phys. Rev. A78, 060102 (2008); P. P´erez-Fern´andez, A. Rela˜no, J. M. Arias, J. Dukelsky, and J. E. Garc´ıa-Ramos, Phys. Rev. A 80, 032111 (2009). [10] P. P´erez-Fern´andez, A. Rela˜no, J. M. Arias, P. Cejnar, J. Dukelsky, and J. E. Garc´ıa-Ramos , Phys. Rev. E 83, 046208 (2011). [11] R. Puebla, A. Rela˜no, and J. Retamosa, Phys. Rev. A 87, 023819 (2013); R. Puebla and A. Rela˜no, Europhys. Lett. 105, 50007 (2013). [12] J. M. Arias, J. Dukelsky, and J. E. Garc´ıa-Ramos, Phys. Rev. Lett. 91, 162502 (2003). [13] P. Cejnar, J. Jolie, and R. F. Casten, Rev. Mod. Phys. 82, 2155 (2010). [14] P. Ribeiro, J. Vidal, and R. Mosseri, Phys. Rev. E 78, 021106 (2008). [15] F. P´erez-Bernal and O. ´ Alvarez-Bajo, Phys. Rev. A 81, 050101(R) (2010). [16] P. P´erez-Fern´andez, P. Cejnar, J. M. Arias, J. Dukelsky, J. E. Garc´ıa-Ramos, and A. Rela˜no, Phys Rev. A 83, 033802 (2011). 24 [17] T. Brandes, Phys. Rev. E 88, 032133 (2013). [18] V. M. Bastidas, P. P´erez-Fern´andez, M. Vogl, and T. Brandes, Phys. Rev. Lett. 112, 140408 (2014). [19] B. Dietz, F. Iachello, M. Miski-Oglu, N. Pietralla, A. Richter, L. von Smekal, and J. Wambach, Phys. Rev. B 88, 104101 (2013). [20] M. A. Alcalde, M. Bucher, C. Emary, and T. Brandes, Phys. Rev. E 86, 012101 (2012). [21] E. Timmermans, P. Tommasini, R. Cˆot´e, M. Hussein, and A. Kerman, Phys. Rev. Lett. 83, 2691 (1999). [22] G.-R. Jin, C. K. Kim, and K. Nahm, Phys. Rev. A 72, 045602 (2005). [23] G. Santos, A. Tonel, A. Foerster, and J. Links, Phys. Rev. A 73, 023609 (2006). [24] A. Motohashi and T. Nikuni, Phys. Rev. A 82, 033631 (2010). [25] J. C. Sanders, O. Odong, J. Javanainen, and M. Mackie, Phys. Rev. A 83, 031607(R) (2011). [26] G. J. Milburn, J. Corney, E. M. Wright, and D. F. Walls, Phys. Rev. A 55, 4318 (1997). [27] A. J. Leggett, Rev. Mod. Phys. 73, 307 (2001). [28] M. Holthaus and S. Stenholm, Eur. Phys. J. B 20, 451, (2001). [29] M. Chuchem, K. Smith-Mannschott, M. Hiller, T. Kottos, A. Vardi, and D. Cohen, Phys. Rev. A 82, 053617 (2010). [30] M. Mel-Messeguer, B. Juli´a-D´ıaz, A. Polls, J. Low Temp. Phys. 165, 180 (2011). [31] J. Dukelsky, C. Esebbag, and P. Schuck, Phys. Rev. Lett. 87, 066403 (2001). [32] H.-Q. Zhou, J. Links, R. H. McKenzie and X. -W. Guan, J. Phys. A 36, L113 (2003). [33] G. Ortiz, R. Somma, J. Dukelsky, S. Rombouts, Nucl.Phys. B 707, 421 (2005). [34] G. Mazzarella, L. Salasnich, A. Parola, and F. Toigo, Phys. Rev. A 83, 053607 (2011). [35] R. Gilmore and D.H. Feng, Nucl. Phys. A 25, 189 (1978). [36] J. P. Blaizot and G. Ripka, ”Quantum Theory of Finite Systems”, The MIT Press (1986). [37] M. A. Bastarrachea-Magnani, S. Lerma-Hern´andez, and J. G. Hirsch, Phys. Rev. A 89, 032102 (2014). [38] P. Str´ansk´y, M. Macek, and P. Cejnar, Annals of Physics 345, 73 (2014) 25