Table Of Content

Quantum Criticality in a Bosonic Josephson Junction P. Buonsante,1,2 R. Burioni,1 E. Vescovi,1 and A. Vezzani3,1 1Dipartimento di Fisica, Universita` degli Studi di Parma, V.le G.P. Usberti n.7/A, 43100 Parma, Italy 2Istituto Nazionale di Ottica-CNR (INO-CNR) and European Laboratory for Non-Linear Spectroscopy (LENS) Via N. Carrara 1, I-50019 Sesto Fiorentino, Italy 3Centro S3, CNR Istituto di Nanoscienze , via Campi 213/a, 41100 Modena, Italy (Dated: January 18, 2012) 2 In this paper we consider a bosonic Josephson junction described by a two-mode Bose-Hubbard 1 model,andwethoroughlyanalyzeaquantumphasetransitionoccurringinthesysteminthelimitof 0 infinitebosonicpopulation. Wediscusstherelation betweenthisquantumphasetransition andthe 2 dynamicalbifurcation occurringin thespectrumof theDiscreteSelfTrappingequationsdescribing n thesystematthesemiclassicallevel. Inparticular,weidentifyfiveregimesdependingonthestrength a oftheeffectiveinteractionamongbosons,andstudythefinite-sizeeffectsarisingfromthefiniteness J ofthebosonicpopulation. Wedevoteaspecialattentiontothecriticalregimewhichreducestothe 7 dynamicalbifurcationpointinthethermodynamiclimitofinfinitebosonicpopulation. Specifically, 1 wehighlightananomalousscalinginthepopulationimbalancebetweenthetwowellsofthetrapping potential,aswellasintwoquantitiesborrowedfromQuantumInformationTheory,i.e. theentropy ] of entanglement and the ground-state fidelity. Our analysis is not limited to the zero temperature h case, but considers thermal effects as well. c e m I. INTRODUCTION cal level. Actually, this regime shrinks with increasing - bosonicpopulation,andreducestothesemiclassicalcrit- t a ical point in the limit of infinite population which, in Several years ago it was suggested that two noninter- t s acting condensates trapped in a double well trap at zero the present setting, has the role of a thermodynamic t. temperature would give rise to the bosonic analog of a limit [24]. Focusing on the “standard” order parame- a ter, i.e. the population imbalance of the two wells, as m Josephsonjunction[1]. Furtherstudiesofsuchabosonic well as on quantities borrowed from quantum informa- Josephsonjunctionfollowed,wheretheeffectofbosonin- - tion theory, such as the entropy of entanglement and d teractionwastakenintoaccountinasemiclassicalframe- the ground-state fidelity , we calculate the critical expo- n work based on the so-called discretized Gross-Pitaevskii o equations[2–5]. Differentinterestingregimesandbehav- nentscharacterizingthequantumphasetransitioninthe c thermodynamic limit. Moreover,we study the finite-size iorswerepredictedforthesystemdependingonthevalue [ scalingeffects,evidencing,withinthe criticalregime,the of the effective interaction among the bosons in each of anomalousscalingofthe consideredquantitiesasa func- 2 the two wells of the trapping potential. These include v the so-calledmacroscopic self-trapping phenomenon and tion of the total boson population. Our analysis is not 6 restricted to the ground-state of the system, i.e. to the a dynamical bifurcation corresponding to a localization 1 zero-temperature case. We also consider finite temper- transition in the ground-state of the system [6]. Just a 8 atures, and evidence a different scaling of the quantum few years after the above mentioned theoretical work, 3 and thermal fluctuations of the population imbalance in . brilliant experiments were carried out with ultra-cold 2 the critical regime. atomic gases where the predicted Josephson oscillations, 1 1 macroscopicnonlinearself-trappinganddynamicalbifur- The plan of the paper is as follows. In Section II we 1 cation were observed [7–9]. describe the quantum Hamiltonian employed in the de- : The quantum version of the discrete Gross-Pitaevskii scriptionofthetwo-modesystem,anddiscussitsrelation v equations—knownasdiscreteself-trapping(DST)equa- withsimilarmodels,aswellaswiththesemiclassicalDST i X tions [10] among the nonlinear community — was also equations; In Section III we illustrate how the low-lying r widely investigated [11–13]. In particular, a significant spectrum of the quantum Hamiltonian can be equiva- a amount of attention [14–23] was directed to the quan- lently investigated by means of a Schrödinger-like equa- tum counterpartof the abovementioneddynamicaltran- tionforafictitiousparticleonaone-dimensionaldomain, sition from a localized to a delocalized regime, which, in and observe that the semiclassical DST equations natu- thispaperwillbeclearlyrecognized,asaquantumphase rallyemergefromasmall-oscillationapproacharoundthe transition. minima of the potential energy. Section IV discusses the In particular we carry out a detailed analysis of the so-called Gaussian regimes, where the Schrödinger-like ground state of the system which, at the quantum level, equation can be safely approximated by a Schrödinger is described by a two-mode Bose-Hubbard model. We equation proper, featuring a harmonic potential part. identify five regimes depending on the effective bosonic The limits of validity of such an approximation are ana- interaction. The central, critical regime corresponds to lyzed in Section V, where further non-Gaussian regimes the dynamical bifurcation occurring in the solutions of arediscussed. Somemoretechnicalissuesconcerningthe the time-independent DST equations at the semiclassi- strong-couplingregimesaredeferredtoAppendixA.The 2 criticalregimeis detailedly analyzedin Sec. VI, whereas is a su(2) coherent state, as discussed in Ref. [18]. The Sec. VII is devoted to a discussion of the behavior of dynamics of the above variables is determined by the two interesting indicators borrowed from Quantum In- expectation value of the trial state on Eq. (1), which formation Theory, the entropy of entanglement and the plays the role of a semiclassical Hamiltonian, ground-statefidelity. The thermalfluctuations ofthe or- derparameterareinvestigatedinSec. VIII,whileSecIX ΩN γ = ΦH Φ = z2 εz 1 z2cosϕ (6) contains our conclusions and perspectives. H h | | i 2 2 − − − (cid:16) p (cid:17) where γ = UN and ε = v are the effective param- Ω Ω II. THE MODEL eters of the model [2]. Specifically, the dynamics is dictated by the so-called Discrete Self-Trapping equa- The Hamiltonian of the system reads tions [10] ensuing from Eq. (6), and the corresponding fixed-point equations are known to exhibit bifurcations H = UH ΩH vH (1) at γ 2/3 1 = ε2/3 [19, 28]. The manifestation of the 2 int− 2 kin− 2 bias sem|ic|lass−icalbif|ur|cations in the features of the quantum system has been repeatedly reported in the literature. where Forinstancetheemergenceofstructuresinthespectrum Hint = a†1a†1a1a1+a†2a†2a2a2 of the quantum problem (1) at the same energies as the semiclassical fixed-point solutions has been investigated =(cid:16)[n1(n1 1)+n2(n2 (cid:17)1)], (2) in Refs. [14, 26, 28–30] − − In the following we will be mainly concerned with the Hkin = a†1a2+a†2a1 (3) mirror-symmetric case, ε = 0, for which the bifurca- Hbias =(cid:16)(n1 n2) (cid:17) (4) tion corresponds to a dynamical transition at γ = 1 − − in the ground-state of Hamiltonian (6). For γ 1 ≥ − and the lattice operators a†j, aj and nj = a†jaj respec- this is symmetric, z = 0, while for γ < −1, due to the tively create, destroy and count bosons at lattice site j. nonlinear nature of the dynamical equations, the mirror The operators in Eqs. (2) and (3) account for interac- symmetry of the Hamiltonian is spontaneously broken, tionsamongbosonsatthesamesiteandfortheirkinetic z = 1 γ 2. energy, respectively. The relative importance of these Al±thoug−h th−e ground-state of the quantum Hamilto- p termsiscontrolledbythe interactionstrengthU andthe nian (1) is always symmetric, n = n , a crossover 1 2 hopping amplitude Ω>0. The operatorHbias allows for in its internal structure can behreciognihzediin the vicin- a local energy offset between the two sites, breaking the ity of the mean-field critical value. Several indicators mirror symmetry of the system. have been employed to highlight this crossover: energy Notice that Hamiltonian (1) commutes with the total gaps[14,31,32],fluctuationsinthepopulationimbalance number operator, N = n1 +n2, which therefore can be [15, 16, 19, 20, 23], structure of the occupation number consideredasaconstantforallpracticalpurposes. Owing distribution[15,16,21,23,33]fidelity[31],entanglement to this fact, the inclusion of a further term of the form entropy[23, 32,34–36], coherence visibility [16, 23], gen- ΓHdip = Γn1n2, accounting for an interaction between eralized purity [37]. Refs. [17, 18, 22, 24, 31] consider bosons atdifferent lattice produces a Hamiltonian ofthe similarissuesonsystemscomprisingmorethantwosites. same class as H [25]. Note indeed that H = 1(N2 dip 2 − WementionthatHamiltonian(1)canbemappedonto N H ). − int a Lipkin-Meshkov-Glick (LMG) model with vanishing Also notice that a simple mapping exists be- anisotropyparameter. Someoftheresultsweillustratein tween the spectral features of Hamiltonian (1) the followinghavebeen obtainedin that contextmaking for attractive U < 0 and repulsive U > 0 in- use of numerical or analitical techniques different from teractions. This can be appreciated by ob- ours[38–42]. serving that D UH + vH ΩH D 1 = 2 int 2 bias− 2 kin − − −U2Hint− v2Hbia(cid:2)s− Ω2Hkin , where D =e(cid:3)iπ(n1−n2) is a unitary transform, [26]. (cid:2)Hamiltonian (1) can be stu(cid:3)died by assuming that the III. THE SCHRO¨DINGER-LIKE EQUATION state of the system is wellapproximatedby a trialwave- FOR THE LOW-LYING EIGENSTATES function of the form Under suitable conditions, the eigenvalue equation for N |Φi= √1N! r1+2 zeiϕ2a†1+r1−2 ze−iϕ2a†2! |0i (5) Hpraomxiimltoanteiadnby(1a),SHch|rΨöidi=ngEer|-Ψliikeiseqvueraytiosantfisofraactfiocritliytioaups- particle on a one-dimensional domain. Several slightly where 1 z 1andϕaretworealdynamicalvariables different versions of this approach have been employed − ≤ ≤ describingthepopulationimbalanceandthemacroscopic the literature [5, 19, 20, 25, 26, 36, 43–46]. All of them phase difference between the two wells [3, 27]. Eq. (5) arebasedonthe expansionofthe eigenstateonthe Fock 3 basis of the site occupation numbers andretainonlythelowestordertermsintheexpansionof the interaction and kinetic term in the eigenvalue equa- ν N a†j tion, we obtain the following Schrödinger-like equation Ψ = c N ν ν , ν = 0 , (7) [20]: | i ν| − i1⊗| i2 | ij (cid:16)√ν(cid:17)! | i ν=0 X ΩN 2 d d 1 z2 +V (z) ψ(z)=E¯ψ(z) (10) where 0 is the vacuum state, aj 0 = 0. Introducing 2 −N2dz − dz ε the nor|mialized population imbalan|cie of the Fock state (cid:20) p (cid:21) relevant to expansion coefficient cν wherez [ 1, 1]isthecontinuumlimitofthemeshgrid ∈ − and the effective potential is [19] (N ν) ν 2ν z = − − =1 (8) ν N − N V (z)= γz2 1 z2+εz (11) ε 2 − − the eigenvalue equation reads p Owing to the N 2 coefficient in frontofthe derivative − UN2 ΩN 1+z 1 z 1 part in Eq. (10), for large populations the low-lying so- z2c ν − ν c 4 ν ν − 2 "s 2 (cid:18) 2 − N(cid:19) ν−1 lutions of such equation will be strongly localized in the vicinity of the minima of the potential energy, Eq. (11). 1 z 1+z 1 vN This suggest that the Schrödinger-likeequation (10) can − ν ν c + z c =E¯c (9) ν+1 ν ν ν be further simplified by introducing the (small) devia- s 2 2 − N # 2 (cid:18) (cid:19) tion from the minima of Eq. (11). Note that the mini- where the ”rescaledenergy” E¯ =E+U 2N N2 /4 mizationofEq.(11)isequivalenttothedeterminationof − − thelowest-energystationarysolutionoftheDiscreteSelf- Ω/2=E+Ω[γ(N 2) 1]differsfromtheoriginaleigen- 2 2 − − (cid:0) (cid:1) Trappingequationsarisingfromthe mean-fieldHamilto- value by a unimportant population-dependent constant. nian(6)[24]. Thisiseasilyverifiedinthemirrorsymmet- Now, for large boson populations N, both the square riccase,ε=0,where,uptothefourthorderinthesmall rootsappearingininEq.(9)andthecoefficientsc can ν 1 displacement from the minimum point(s), u = z z , beexpandedinpowersofN 1. Noteindeedthatt±helat- − γ − equation (10) is equivalent to ter can be seen as the values that a continuous function ψ takes on at the points of the grid mesh in the interval d2 [ 1,1]definedbyz withν =0,1,2, ,N. Specifically, +Au2+Bu3+Cu4 ϕ(u)= ϕ(u) (12) − ν ··· −du2 E if we set c = 2 ψ(z ) = 2 ψ(z 2N 1) [47] with(cid:20)ψ(z)=ϕ(z z ) and (cid:21) ν±1 N ν±1 N ν ∓ − − γ q q z A B C γ E N2(γ+1) N2 E¯N N2 γ 1 0 0 + ≥− 4 16 Ω 2 1 N2γ2(γ2 1) N2γ5 γ2 1 N2γ6(5γ2 4) N2(1+γ2) E¯N γ < 1 1 − − − γ − ± − γ2 4 ∓ 4 16 4 − Ω r p (13) IV. HARMONIC REGIME where 1 We note that the leading order in Eq. (12) is almost γ > 1 1 N√γ+1 − everywhere the second, except at the mean-field critical σ2 = = (15) fpoouinrtthγ. T=he−re1fowrhe,efroertthheefilorwstesntoenn-evragnyislhevinegls,tearwmayisfrtohme γ 2√A  N γ 1γ2 1 γ <−1 | | − criticality and for sufficiently large populations, we can  p let B = C = 0 and Eq. (12) reduces to the Schrödinger Actually, for γ < 1 Eq. (12) describes the small − equation for a harmonic oscillator, whose discrete spec- oscillations about either of the two equivalent minima trum is of the effective potential, z = 1+γ 2. Therefore γ ± − a generic solution of Eq. (10) is well approximated by =√4A n+ 1 = 1 n+ 1 (14) a superposition of one harmonic spolution for each well. En 2 σ2 2 The symmetry of the problem dictates that these super- (cid:18) (cid:19) γ (cid:18) (cid:19) 4 positions have a definite parity, 2 ϕ (z z ) ϕ (z+ z ) ψn±(z)= n −| γ| √±2 n | γ| (16) 1.5 Ej 1 ∆ Also, since in this approximationthe two localized func- 0.5 tions ϕ (z z ) have a vanishing overlap, both ψ (z) have annene±rgyγvery close to that in Eq. (14). Inn±par- 0 2 ticular, the ground-state of the system ψ (z) is well ap- 0 proximated by a symmetric superposition ψ+(z) of two 1.5 0 Gaussian functions of the form Ej 1 ∆ ϕ0(u)= 1 e−4uσ2γ2 (17) 0.5 (σ √2π)1/2 γ 0 −2 −1.5 −1 −0.5 0 γ with σ2 given in (the lower of) Eq. (15). γ For γ > 1 the effective potential has a single mini- − mum atz =0,andthe problemcanbe mapped exactly γ FIG.1: (coloronline)LowlyingenergygapsforN =50(up- onto a quantum harmonic oscillator, as far as the width perpanel)andN =5000(lowerpanel). Graythicklinesrep- of the eigenfunctions of the latter problem do not ex- resenttheanalyticresultsinEqs. (22)and(23). Specifically, ceed the interval in which the harmonic term dominates the dashed line is Ω γ2 1 and the solid line is Ω√γ+1. over higher-order terms. The energy spectrum is again Symbols represent npumer−ical data. Specifically, circles, up- given by Eq. (14) and, in particular, the ground-state is wardtriangles,squaresanddownwardtrianglesare∆E with n ψ (z)=ϕ (z), Eq. (17), with σ2 givenin (the upper of) n=1,2,3,4, respectively. 0 0 γ Eq. (15). The most natural physical quantity which is in prin- ciple measurable, is given by the population imbalance γ1 γ γ2 each energy level is (almost) twofold de- ≪ ≪ corresponding to the operator [20] in the original two- generate, so that we have site problem: 0, γ γ γ zˆ= a†1a1−a†2a2, (18) ∆E2n+1 ≈(cid:26)Ω√γ+1, γ31 ≪≪γ ≪≪γ42 (22) N and Its value and fluctuations are easily evaluated on the Ω γ2 1, γ γ γ ground state of the system ∆E2n ≈ Ω√γ+−1, γ13 ≪γ ≪γ24 (23) (cid:26) p ≪ ≪ zˆ = dsz ψ (z)2 =0 (19) In Fig. 1 we compare Eqs. (22) and (23) with numerical 0 h i | | data for ∆E . The plots clearly show that the analytic Z n approximationsbecomemoreaccuratewithwithincreas- and [48] ing boson population. zˆ2 = dsz2 ψ (z)2 σ2 +z2 (20) h i | 0 | ≈ γ γ Z V. LIMITS OF VALIDITY AND REGIMES where denotes the expectation value on the ground state ofh·tihe system. Moreover one can also consider the In the present section we describe the limits of valid- expectationvalue onthenthexcitedstateobtaining: ity of the harmonic solutions described above and the n h·i different regimes taking over at different values of the effective parameter γ. A crucial assumption in passing zˆ2 = dsz2 ψ (z)2 σ2(2n+1)+z2 (21) h in | n | ≈ γ γ fromEq.(9)toEq.(10)isthatthe coefficients cν canbe Z regarded as continuous functions of the population im- The fluctuations of the population imbalance have balance of the Fock state defined in Eq. (8). This may been investigated by several authors, [15, 16, 19, 23], be problematic for eigenstates in the mid-spectrum of andanalyticalexpressionsforthisquantityhavebeenre- Hamiltonian (1), but is very reasonable for states close portedinRefs. [15,19,23]. InbothGaussianregimesthe to the extrema ofthe same spectrum. However,the con- (low lying) energy spectrum of Eq. (12) is well approxi- tinuousapproachcanproveillposedevenforthe ground matedbyEq. (14). SinceEn =E¯n−Ω2[γ2(N2−2N)−1], state of the system. It is clear that, for the continuous and E¯ depends on as described in table (13), it is limitofEq.(10)tobewelldefined,theGaussinwidthσ n n γ E easy to calculate the lowest energy gaps in the spectrum must be much larger than distance between two neigh- of H, ∆E = E E . As we mentioned earlier, for boringpointsinthediscreteEq.(9),namely2N 1. This n n n 1 − − − 5 means γ γ √N for γ < 1 and γ γ N2 VI. CRITICAL REGIME 1 4 ≫ ∼ − − | | ≪ ∼ for γ > 1 [25]. − According to our previous analysis, as the effective We observe that the same estimates can be obtained parameter γ approaches the mean-field critical value by imposing that the expansion parameter in the strong γ = 1, the system exhibits the hallmarks of a phase c − coupling regimes be actually small, as discussed in some transition, where the minimum z of the effective po- γ detail in Sec. A. This ensures that for γ γ and γ tential, table (13), plays the role of the order parame- 1 ≫ ≪ γ the first-order perturbative description of the ground ter. Note thatthis quantityis inprinciplemeasurableas 4 state applies. the absolute value of the population imbalance operator z = zˆ . The fluctuations of this order parameter γ h| |i The parameter range γ1 γ γ4 where the continu- zˆ2 zˆ 2 =σ2 (24) ous approach applies can b≪e div≪ided into three regimes. h| | i−h| |i γ For γ γ γ the ground state of the system is well coincide with the Gaussian width of Eq. (15). 1 2 approxi≪mate≪dbyasymmetricsuperpositionoftwoGaus- Thebehaviorforzγ andσγ2 describedintable(13)and sians. On the other hand, for γ γ γ the ground Eq. (24),respectively,holdstrueinthemirror-symmetric 3 4 ≪ ≪ state is well approximated by a single Gaussian. case, v = ε = 0. It proves useful to study behaviour of theorderparameterinthepresenceofavanishingenergy offset between the sites of the dimer, v 1. According These two regimes are separated by an intermediate | |≪ to ourdiscussionsaboutthe Gaussianregimes,the value critical regime γ < γ < γ , surrounding the mean- 2 3 of the order parameter will coincide with the absolute field critical point γ = 1, where the higher order c − minimum of the effective potential which, owing to the terms in Eq. (12) cannot be neglected and the prob- presence of a non-vanishing ε, is now non degenerate for lem is not Gaussian any more. The boundaries of this any value of γ. The minimization of Eq. (11) then gives regime can be obtained from a comparison between the leading and next-to-leading term in the perturbative ex- zˆ ε γ < 1 pmaantseiofnorotfhethleoweeffrecbtoivuendpoγten<tiaγl. cAanrbeeasoobntaabilneedesetii-- h|zˆ|iε =(cid:26) hh||zˆ||iiεε==00+− γγ(+εγ12−1) γ >−−1 (25) 2 c ther by imposing Aσ2 B σ3 or by requiring that the whereonlylineartermsinεhavebeentakenintoaccount. γ ≫| | γ distance between the two equivalent harmonic minima, Eq. (25) allows us to introduce a sort of susceptibility 2 1 γ 2, is much larger than the width of the rele- vγapcn−tγ−G2a∼u−ssNia−n2/g3r.ouSnimdislatartlye,stσhγe.uIpnpbeortbhocuansdesγo3n>e fiγncdiss χz = ddǫh|zˆ|iε =( |γγ+1|(1γ12−1) γγ ><−−11 (26) obtained by requiring Aσ2 Cσ4, which results once γ ≫ γ Equations (13), (15), (24), (26) provide some critical again into γ γ N 2/3 [19, 45]. As discussed in 3 − c ∼ − exponents for the system. As γ approaches its critical Sec. VI,theenergygapsinthecriticalregimedependon value γ = 1, the order parameter vanishes (from be- wthheatsyhsatepmpenpsopinultahtieonG,a∆usEsia∼n rNeg−im1/e3s,.aTthviasriparnocveidwesitha low) asch|zˆ|−i ∼ (γc − γ)αz, its fluctuations diverge as pNr−eddicetpeedndinenEtqcsu.t(o2ff2)toantdhe(2v3a)n.ishing of ∆E at γ = γc ha|szˆ|χ2zi−∼h||γzˆ|i−2γ∼c||−γα−χ,γwc|h−eαrσeanditssusceptibilitydiverges 1 α =α = , α =1 (27) z σ χ In summary, one can recognize five different regimes 2 depending on the value of the effective parameter γ. For The above discussion rigorously applies in the limit of γ γ andγ γ thecontinuousapproximationofSec- infinite population, which plays the role of an effective 1 4 ≪ ≫ tionIIIdoesnotapply,butthesystemcanbeanalyzedby thermodynamic limit. Indeed, as we observe in Section resorting to first-order perturbative theory, as discussed V, only in this limit does the critical regime shrink to a in Appendix A. For γ γ γ the continuous ap- single point corresponding to the critical value. For any 1 2 ≪ ≪ proximation applies, and the ground-state of the system finitepopulationsthereexistsafiniteregionsurrounding is a symmetric superposition of two Gaussians centered thecriticalpointwheretheGaussianresultsdonotapply. at the equivalent harmonic minima of the effective po- In this region the quartic term dominates and Eq. (10) tential. For γ γ γ the system is critical, in a can be approximated as 2 3 ≪ ≪ sensethatwillbe clarifiedshortly. Note thatinthe limit d2 of large populations the critical regime reduces to the +Cz4 ψ (z)= ψ (z) (28) mean-field critical point. For γ γ γ the ground- −dz2 n En n 3 4 (cid:20) (cid:21) ≪ ≪ state is a Gaussian centered at the unique minimum of whereC islistedintheupperrowoftable(13). Itiseasy the effective potential, z = 0. Note that in estimating γ to provethat in this criticalregime the eigenfunctions of the location of the boundaries γ we only considered 1 4 Eq. (28) have form the dependence on the (large) bo−sonic population, and neglected other numerical factors. ψ (z)=N1/6φ (zN1/3), (29) n n 6 2.2 0.7 N = 400 N = 800 2 )0.6 N = 1600 2 i N = 3200 | 1.8 h|z0.5 N = 6400 − E n1.6 2i0.4 ∆ z| 3 h| 1/N 1.4 3/(0.3 2 N 1.2 0.2 1 0.1 −6 −4 −2 0 2 4 0.8 N2/3(γ − γ c) 102 103 104 N FIG. 3: (color online) Scaling of the fluctuations of the (absolute valueof) thepopulation imbalance operator. FIG. 2: (color online) Dependence on the population of the first few energy gaps ∆E at γ = 1 (from bottom to top, n − n=1,2 ,5). The data points are numerically determined ··· values, thesolid lines are just guides to the eye. ThecomparisonwithEqs. (15),(24)and(31)thenallows us to conclude that where φn(ζ) is the solution of the population- 2 3 independent problem τ = , ξ = (34) −3 2 d2 ζ4 −dζ2 + 16 φn(ζ)=E¯nφn(ζ) (30) The correctness of the scaling assumption in Eqs. (33)- (cid:20) (cid:21) (34)isdemonstratedinFig. 3,whereanicedatacollapse and the eigenvalue corresponding to ψ (z) is = is observed for a wide range of boson populations. We n n N2/3¯ . This result allows us to estimate the behaEviour mention that – in the context of the LMG model – the n E of the (low-lying) energy gaps in the spectrum of Hamil- criticalexponentsin(34)weredeterminednumericallyin tonian (1) in the critical regime. We get ∆E N 1/3, Ref.[38]andmakinguseofananalyticapproachdifferent n − ∼ where we used the definition of in table (13) and the from ours in Refs. [39–41]. E fact that the spectrum of problem (30) does not depend We conclude by remarking that the rescaled energies onN. This estimate,whose correctnessdemonstratedin ¯ and the expectation values for the rescaled problem n E Fig. 2, providesanalternateroute forthe determination (30) can be estimated by applying the Bohr-Sommerfeld oftheboundariesofthecriticalregion. Indeed,thesecan quantization rule. This gives be assessedbyrequiringthatthe predictioninthe Gaus- saisanthraetgiinmeths,eEcqrist.ic(a2l2r)egainmde,(2i.3e).,Nar1e/3o.fTthheissagmiveesoorndceer 2En1/4 ¯ ζ4dζ =π(n+ 1), n N, (35) n again γc−γ2 ∼N−2/3 and γ3−γc ∼N−2/3. Z−2En1/4rE − 16 2 ∈ As to the fluctuations of the order parameter, partic- whichprovidesquite satisfactoryestimates ofthe energy ularizing Eq. (29) to the ground-state of the system, levels of the system, n=0, we get h|zˆ|2i−h|zˆ|i2 =N−2/3 dζdζζ2φ|φ(0ζ()ζ)2|2 (31) En =ΩC1N−1/3(n+ 21)2/3 (36) R | 0 | Within a standard scaling approRach it is reasonable to where C is a suitable constant. Within the same ap- 1 assume that there exists a correlation length, dictating proximation the mean square displacement of the eigen- the effective size of the system, which diverges at the function φ (ζ) can be evaluated as well. In particular, n critical point for the simclassical orbit of energy ¯ , one gets n E γ γ ξ. (32) γ c − N ∼| − | ζ4 Any physical quantity should depend on the boson pop- φn(ζ)2ζ2dζ ζ2δ(p2+ ¯n)dpdζ ≈ 16 −E ulationN onlythroughtheratioN/ ,sothatitshould Z Z be Nγ =C (n+ 1)2/3 (37) 2 2 N zˆ2 zˆ 2 =Nτf (33) h| | i−h| |i (cid:18)Nγ(cid:19) where C2 is a suitable constant. 7 VII. QIT INDICATORS 104 104 As we mention above,quantities borrowedfromquan- tum information theory have been employed to investi- 103 gate the the incipient quantumphase transitionbetween the localizedanddelocalizedground-stateofthe system. 102 102 Owingtotheparticularlysimplestructureofthesystem, F someofthesequantitiesturnouttobeentirelyequivalent χ 101 to usual indicators. This is the case of the generalized 100 −1.1 −1.05 −1 purity ofRef. [37],which basicallycoincideswith the coherence visibility consideredin Refs. [16, 23, 36], i.e. the off-diagonal term of the one-body density matrix. Simi- larly, the Fisher information of Ref. [23] is equivalent to 10−2 the population imbalance fluctuations considered in sev- eralearlierworks[15,16,19,20]aswellasinthe present −2 −1.5 −1 −0.5 0 0.5 1 γ paper. The fidelity and entanglement entropy considered in Refs. [31] and [23], respectively, appear to be less trivial FIG. 4: (color online) The analytic expression obtained for thefidelitysusceptibilityintheGaussianregimesiscompared quantities [49]. The former is nothing but the overlap withthenumericalresultsforthesamequantity. Thecircles, between two ground-states of the system corresponding upwardtriangles, squares, downward triangles and diamonds to slightly different values of the parameter driving the are obtained from exact diagonalization of Eq. (1) for sys- transition Ψ Ψ [50]. Since this quantity is almost γ γ+κ tems containing 500, 1000, 2000, 4000 and 8000 bosons, re- h | i everywhere extremely close to unity, a more convenient spectively. In all cases we used Eq. (38) with κ=10−4. The indicatorisprovidedbytheso-calledfidelitysusceptibility thick gray lines represent Eq. (40) for thesame boson popu- [51] lations. Notice how χF(γ) depends on the boson population only for γ < γc. The inset contains a magnification of the 1 Ψ Ψ 1 d2 region enclosed by the dashed black frame, where Eq. (40) γ γ+κ χ = lim −h | i = Ψ Ψ (38) f κ 0 κ2 2dκ2h γ| γ+κi κ=0 does not apply. The thin solid lines are guides to the eye. → (cid:12) (cid:12) Inabipartitesystem,theentropyofentangle(cid:12)mentofa pure state Ψ is givenby Tr(ρ logρ ), where ρ is the ofthe entanglemententropyis quite differentin the spin j j j reducedden|siitymatrixofs−ubsystemj,obtainedfromthe and boson context, and the results are not readily com- fulldensitymatrix Ψ Ψ aftertracingoverthestatesof parable. As to the ground-state fidelity, a comparison the othersubsystem| [5ih2].|Inthecaseunderinvestigation will be made with the results in Ref. [42]. the two subsystem are of course the two sites composing Based on the results sketched in the previous sections the Bose-Hubbard dimer and, owing to number conser- and in the appendices, we are able to give an analytic vation, the entropy of entanglement of an eigenstate of description of the behaviour of the fidelity susceptibility Hamiltonian (1) is and entropy of entanglement, Eqs. (38) and (39). Let us start with the former quantity. As we discuss N in Section III, in the Gaussian regimes the ground-state 1 S = c 2log c 2 (39) of the system is strictly related to the ground state of e −log (N +1) | ν| 2| ν| 2 Xν=0 a harmonic oscillator. Specifically, for γ1 ≪ γ ≪ γ2 the ground state is well approximated by a symmet- where the cν’s are the expansion coefficients in Eq. (7) ric superposition , ψ+(z), of two Gaussian functions 0 and we introduced a normalization factor. This quan- 1 tity has been considered in a few earlier works on the of square width σγ2 = N|γ| γ2−1 − centered at Bacotsioe-nHsuibtbwarads duismeder.toInintvheestpigraesteentcheeofprreecpuurlssoivreoifnttehre- Lzγik=ewi±se, 1fo−r γγ−l2l,γas deγs(cid:16)crtibheedpgbroyunEdqs(cid:17)s.ta(t1e6)isanwdell(1d7e)-. 3 4 Mott insulator–superfluid quantum phase transition oc- scribedpby a single G≪aussianϕ (z) of square width σ2 = 0 γ curring on infinite lattices [53]. In attractive systems 1 it was employed for the characterizationof Schro¨dinger- N√1+γ − centered at zγ = 0, as described by Eq. (17). Straightforwardcalculations result in cat-like states in the presence of a small energy offset v (cid:0) (cid:1) bofoferttwshmeeeaqnlultafihnnetiuttwemobcosoistuoennstopefroptphauerltastyoifsotntehsme[2[m336]e.]a,na-nfidelidnttrhaensstituidony χF(γ)= 8|γ|31√Nγ2−1 + 16γ(22γ(γ−21−)21)2 γγ1≪γγ≪γγ2 (40)  64(γ+1)2 3≪ ≪ 4 These two indicators have been considered for the LMG model [42, 54] which, as we mention, can be whose correctness is demonstrated in Fig. 4. It is mapped onto Hamiltonian (1). However,the natural de- worth noticing that Eq. (40) predicts a dependence on composition of the system necessary for the calculation the boson population only in the lower Gaussian regime 8 γ γ γ , whereas for γ γ γ the fidelity su1sc≪eptibi≪lity2turns out to be N3 ≪indep≪ende4nt. This can 104 − be explained by the fact that in the latter regime the small variation κ in the effective parameter affects only thewidthofthesingleGaussianrepresentingtheground- state of the system. Conversely, for γ γ γ , the 103 1 2 ≪ ≪ ) change affects both the centers and the widths of the c γ two Gaussians comprising the ground state of the sys- ( F tem. Also notice that Eq. (40) predicts a divergence in χ the fidelity susceptibility atthe mean-fieldcriticalpoint, 102 γ =γ . However,foranyfiniteN,thereisasmallregion c surrounding such point where the Gaussian approxima- tions do not apply. This is demonstrated in the inset of Fig. 4, which clearly shows that for finite populatons χF(γ) features a finite peak in the vicinity of γc. The 101 103 104 larger the population, the sharper and the closer to the N critical point is the peak. In order to study the behavior of χ (γ) in the critical FIG. 5: Scaling of the fidelity susceptibility at the mean- F region it proves useful to exploit Eq. (28). As it is dis- field critical value. The data points are obtained from exact diagonalization of system containing N = 500 2n bosons, cussed in Ref. [55], the limit in Eq. (38) can be carried · with n = 0, ,5. The thin solid line corresponds to 7.1 out exactly through a perturbative expansion involving 10−3N4/3. ··· · the entire spectrum of the problem at a given value of thedrivingparameterγ. Consideringtheharmonicterm in Eq. (28) as a perturbation over the quartic term, we where get log πe γ γ γ N4/3 dζζ2φ (ζ)φ (ζ) 2 C(γ)= 2(cid:18)2|γ|√γ2−1(cid:19) 1≪ ≪ 2 (44) χF(γc)= 16 ¯ 0 ¯ n (41) log2 2√πγe+1 γ3≪γ≪γ4 nX6=0(cid:20)R E0−En (cid:21) Note that,despite(cid:16)Eq. (39(cid:17)) forbids that Se(γ) exceed where the eigenfunctions φ (ζ) and eigenvalues ¯ are unity, the quantity in Eq. (43) diverges in the vicinity of n n dfined by (29) and (30). Rigorously speaking theEφ (ζ) the mean-fieldcriticalpoint. However,this is not alarm- n provides a good approximation for an eigenstate of the ing, since at any finite boson population there is a small original problem, Eq. (1), only for sufficiently small val- intervalaroundsuchpoint where the ground-state is not ues of the quantum number n. However, only the first GaussianandEq. (43) does not apply. It could be easily few termsinthe sumappearinginEq. (41)areexpected proven that the requirement that Eq. (43) be less than to provide a significant contribution, owing to the lo- unityprovidesanestimatefortheboundariesofthenon- calization of the ground state φ (ζ) and the energy gap Gaussian regime that is equivalent to the one derived in 0 appearing in the denominator. Section V. An interesting feature revealed by Eqs. (43)- (44) is that in the N limit the entanglement en- The N4/3 behavior predicted in Eq. (41) and con- → ∞ tropytendstoaconstant,independentoftheinteraction firmed by the numerical results in Fig. 5 can be rec- to hopping ratio, S (γ) = 1 γ = γ . Note indeed that ognized as the superextensive divergence of the fidelity e 2 ∀ 6 c quantity in Eq. (44) does not depend on N, so that the susceptibility which, according to Ref. [51], is the hall- secondterminEq.(43)vanisheswithincreasingpopula- markofaquantumphasetransition. Thesameexponent tion. wasestimatednumericallyinRef.[42]. Also,ourEq.(40) The behavior of the entropy of entanglement in the agrees with a similar result in the same paper, although critical regime can be obtained by plugging ψ(z) = we obtain a different behaviorof the sub-leading terms. N1/6φ (N1/3z) into Eq. (42), where φ (ζ) is once again Let us now turn to the entropy of entanglement. Re- 0 0 the ground state of Eq. (30). Straightforward calcula- calling that the continuous version of Eq. (39) is tions give Se ≈−log12N Z−11dz|ψ(z)|2log2(cid:18)N2 |ψ(z)|2(cid:19), (42) Se(γ)= 23 − log12N Z dζ|φ0(ζ)|2log2(cid:2)2|φ0(ζ)|2(cid:3) (45) Since the integral in the second term of Eq. (45) does it is easy to check that in the Gaussian regimes notdependonN,wegetthatlim S (γ )= 2. Thus the mean-field critical point is siNgn→a∞lledebycthe3entropy 1 C(γ) of entanglement as a peak which, in the large popula- S (γ)= + (43) e 2 2log N tion limit, turns into a discontinuity. Figure 6 shows 2 9 18 1 N =50 N =500 16 0.9 N =5000 1 14 N =∞ − 0.8 ] ) N =∞ 2312 γ − (0.7 e ) S γc 10 ( 0.6 e S 8 [ 0.5 6 0.4 4 −6 −4 −2 0 2 4 102 104 106 γ N FIG. 6: (color online) The analytic expression for the entan- FIG. 7: Scaling of the entanglement entropy at the mean- glemententropy,Eq.(43),iscomparedwithnumericalresults field critical value. The data points are obtained from exact atdifferentbosonpopulations. Numericaldata(symbols)and diagonalizationofsystemcontainingN =50 5n bosons,with · analyticresults(solidlines)correspondingtothesameboson n=0, ,5. The thin solid line corresponds to 3.4log N 2 ··· − populationhavethesamecolor. Theblackdashedhorizontal 0.11. line is the limit of the Gaussian result, Eq. (43). The black star is thethermodynamiclimit of theentanglement entropy atthemean-fieldcriticalpoint,Eq.(45). Theverticaldashed where the thermal fluctuations are line signals the mean-field critical point. coth(cid:18)Ω√2kγB2T−1(cid:19), γ γ γ σ2 = N γ√γ2 1 1 ≪ ≪ 2 (48) asucltosmopbatraiisnoendbveiatweexeanctEdqisa.g(o4n3a)l-i(z4a4t)ioann.dTnhuemagerreiceamlernet- T,γ  cothN|(cid:16)√Ω|2γ√k+Bγ1+T−1(cid:17), γ3 ≪γ ≪γ4 is quite satisfactory,and improveswith increasing boson and zγ is givenin Table (13). Therefore for kBT <∆E2 population. Note that even at relatively large popula- (i.e. the lowest energy gap of the system) σ2 is es- T,γ tions the entropy of entanglement is still rather different sentially givenby the fluctations of the quantum ground from its limit. This is because the finite size corrections states (20), whereas for large temperatures the thermal decrease logarithmically with N. The different behavior fluctuation dominates and we have at criticality is demonstrated in Fig 7, which shows the 2kBT , γ γ γ correctness of Eq. (45). σ2 ΩN|γ|(γ2−1) 1 ≪ ≪ 2 (49) T,γ ≈ 2kBT , γ γ γ  ΩN(γ+1) 3 ≪ ≪ 4 Eq. (49) highlights the typical linear dependence on VIII. THERMAL FLUCTUATIONS temperatureofthethermalfluctuationsinharmonicsys- tems. We remark that both thermal and quantum fluc- tuationsdecreasewiththetotalpopulationas1/N. This We now turn to the thermal fluctuations of the order is demonstrated by the data collapse in Fig. 8. parameter zˆ Alsointhethermalregimeweexpecttheharmonicap- proximationto apply only whenthe higher ordertherms zˆ2 zˆ 2 = ne−βEn hzˆ2in−hzî2n (46) in Eq. (12) are negligible. This means that for γ1 < h iT −h iT P n(cid:0)e−βEn (cid:1) γ < γ2, σT2,γ should be smaller than A/B while for γ < γ < γ , σ2 should be smaller than A/C. In both 3 4 T,γ P cases for γ 1 we get kT . N(γ +1)2. Note indeed where and denotes the thermal averageand the T n ≈ − h·i h·i that, as evidentfromFig.8,the higher the temperature, quantum expectation value on the state of energy E , n the larger must the population be for the harmonic ap- respectivley. proximation (dash-dotted line) to apply. Evaluating E , zˆ2 and zˆ using standard proper- ties of quantumnhahmoinnic oschilliantors (see Section IV) we Let us now consider the regime γ2 < γ < γ3 where the quartic term dominates. The values of E , zˆ and get zˆ2 can be estimated within the semiclasnsichaliTBohr- n h i Sommerfeld approximation discussed in Sec. VI. Plug- zˆ2 zˆ 2 =z2+σ2 (47) ging Eqs. (36) and (37) into Eq. (46) we get, for h iT −h iT γ T,γ 10 N = 2000 102 2T,γ101 NN == 1200000 σ 3 N = 100 / 2 N 2NσT,γ101 100 A. 101 100 γ 2T, σ 2 / 1N 100 10−1 10−2 10−1 100 101 102 B. kT/Ω 10−2 10−1 100 101 102 FIG. 8: (color online) Scaling of imbalance fluctuations as kBT/Ω a function of temperature for γ = 0.2. Black, red, blue − andmagentasolidlinecorrespondto100,200,1000and2000 FIG. 9: (color online) Scaling of imbalance fluctuations as a bosons, respectively. Dashed and dash-dotted lineis thethe- function of temperature for γ = γc = 1 and different bo- − oretical predictions (15) and (49), respectively. son populations. Thesamedataare multiplied bya different function of the boson population in the two panels, to highlight the different scaling laws in the quantum and thermal γ =γ = 1, regime. ThethickgraylineinpanelBisthetypicalbehaviour c − oftheclassical thermalfluctuationsforapurequarticpoten- zˆ2 zˆ 2 σ ne−βNC11/Ω3(n+21)34C2(n+ 12)23. tial, Eq.(51). h iT −h iT ≈ γ,T ≈ P N23 ne−βNC11/Ω3(n+21)34 1 (50) It is easy to check that also in tPhis critical regime σ IV KT/Ω∼N(γ+1)2 γ,T reducestothepurelyquantumcontributionwhenkBT is 0.8 KT/Ω∼|γ+1|1/2 smallerthantheenergygapbetweenthelowesteigenpair, II C Ω/N1/3. At higher temperatures we have KT/Ω∼N−1/3 1 0.6 Ω / γ+1∼±N−1/3 e−βNC11/Ω3n43C2n23dn kBT 21 KT0.4 σγ,T ≈ RN32 e−βNC11/Ω3n43dn ≈(cid:18)NΩ(cid:19) , (51) II I 0.2 I which is the typicRal behaviour of classical thermal fluc- III tuations when a pure quartic potential is present. We 0 remark that within the critical region the fluctuations −1.1 −1 −0.9 −0.8 γ of the population imbalance due quantum effects scale as N 2/3 with increasing population, while the thermal − fluctuations scale as N 1/2. This is shown in the upper FIG.10: (coloronline)Inthisfigurewerepresentthedifferent − thermalandquantumregimesintheregionaroudthecritical and lower panel of Fig. 9, respectively. pointγ = 1thedataareobtainedforN =100. InregionsI- The different quantum and classical regimes around − IVrepresentquantumharmonic,thermalharmonic,quantum γ = 1 can be summarized as in Figure 10. Note that quarticand thermal non-harmonicbehaviours respectively. − theboundariesbetweendifferentregions,signalledbythe coloured lines, are actually crossovers between different behviours. The only true critical point is γ = 1, in the tuations should decay as N 1/2 for not too large tem- − − limit of infinite N. peratures,andthey shouldbe independent ofN forvery In regions I and III σ2 coincides with the result ap- γ,T large kBT/Ω. plying at zero-temperature for purely quartic and har- monicpotentials,respectively. RegionII ischaracterized by thermal fluctuations in a harmonic potential, while IX. CONCLUSIONS AND PERSPECTIVES in region IV we have a classical behaviour in an anhar- monic potential. InregionsI and II σ2 decaysas N 1, γ,T − while in region III quantum critical effects produce fluc- In this paper we addressed the quantum counterpart tuations decaying as N 2/3. Finally, in region IV fluc- ofthetransitionfromalocalizedtoadelocalizedground- −