{'fraction_non_alphanumeric': 0.043226948393163336, 'fraction_numerical': 0.018059627191769456, 'mean_word_length': 4.3481733471379975, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum Criticality in a Bosonic Josephson Junction

In this paper we consider a bosonic Josephson junction described by a two-mode Bose-Hubbard model, and we thoroughly analyze a quantum phase transition occurring in the system in the limit of infinite bosonic population. We discuss the relation between this quantum phase transition and the dynamical bifurcation occurring in the spectrum of the Discrete Self Trapping equations describing the system at the semiclassical level. In particular, we identify five regimes depending on the strength of the effective interaction among bosons, and study the finite-size effects arising from the finiteness of the bosonic population. We devote a special attention to the critical regime which reduces to the dynamical bifurcation point in the thermodynamic limit of infinite bosonic population. Specifically, we highlight an anomalous scaling in the population imbalance between the two wells of the trapping potential, as well as in two quantities borrowed from Quantum Information Theory, i.e. the entropy of entanglement and the ground-state fidelity. Our analysis is not limited to the zero temperature case, but considers thermal effects as well.


I. INTRODUCTION
Several years ago it was suggested that two noninteracting condensates trapped in a double well trap at zero temperature would give rise to the bosonic analog of a Josephson junction [1]. Further studies of such a bosonic Josephson junction followed, where the effect of boson interaction was taken into account in a semiclassical framework based on the so-called discretized Gross-Pitaevskii equations [2][3][4][5]. Different interesting regimes and behaviors were predicted for the system depending on the value of the effective interaction among the bosons in each of the two wells of the trapping potential. These include the so-called macroscopic self-trapping phenomenon and a dynamical bifurcation corresponding to a localization transition in the ground-state of the system [6]. Just a few years after the above mentioned theoretical work, brilliant experiments were carried out with ultra-cold atomic gases where the predicted Josephson oscillations, macroscopic nonlinear self-trapping and dynamical bifurcation were observed [7][8][9].
In particular we carry out a detailed analysis of the ground state of the system which, at the quantum level, is described by a two-mode Bose-Hubbard model. We identify five regimes depending on the effective bosonic interaction. The central, critical regime corresponds to the dynamical bifurcation occurring in the solutions of the time-independent DST equations at the semiclassi-cal level. Actually, this regime shrinks with increasing bosonic population, and reduces to the semiclassical critical point in the limit of infinite population which, in the present setting, has the role of a thermodynamic limit [24]. Focusing on the "standard" order parameter, i.e. the population imbalance of the two wells, as well as on quantities borrowed from quantum information theory, such as the entropy of entanglement and the ground-state fidelity , we calculate the critical exponents characterizing the quantum phase transition in the thermodynamic limit. Moreover, we study the finite-size scaling effects, evidencing, within the critical regime, the anomalous scaling of the considered quantities as a function of the total boson population. Our analysis is not restricted to the ground-state of the system, i.e. to the zero-temperature case. We also consider finite temperatures, and evidence a different scaling of the quantum and thermal fluctuations of the population imbalance in the critical regime.
The plan of the paper is as follows. In Section II we describe the quantum Hamiltonian employed in the description of the two-mode system, and discuss its relation with similar models, as well as with the semiclassical DST equations; In Section III we illustrate how the low-lying spectrum of the quantum Hamiltonian can be equivalently investigated by means of a Schrödinger-like equation for a fictitious particle on a one-dimensional domain, and observe that the semiclassical DST equations naturally emerge from a small-oscillation approach around the minima of the potential energy. Section IV discusses the so-called Gaussian regimes, where the Schrödinger-like equation can be safely approximated by a Schrödinger equation proper, featuring a harmonic potential part. The limits of validity of such an approximation are analyzed in Section V, where further non-Gaussian regimes are discussed. Some more technical issues concerning the strong-coupling regimes are deferred to Appendix A. The critical regime is detailedly analyzed in Sec. VI, whereas Sec. VII is devoted to a discussion of the behavior of two interesting indicators borrowed from Quantum Information Theory, the entropy of entanglement and the ground-state fidelity. The thermal fluctuations of the order parameter are investigated in Sec. VIII, while Sec IX contains our conclusions and perspectives.

II. THE MODEL
The Hamiltonian of the system reads where and the lattice operators a † j , a j and n j = a † j a j respectively create, destroy and count bosons at lattice site j. The operators in Eqs. (2) and (3) account for interactions among bosons at the same site and for their kinetic energy, respectively. The relative importance of these terms is controlled by the interaction strength U and the hopping amplitude Ω > 0. The operator H bias allows for a local energy offset between the two sites, breaking the mirror symmetry of the system.
Notice that Hamiltonian (1) commutes with the total number operator, N = n 1 + n 2 , which therefore can be considered as a constant for all practical purposes. Owing to this fact, the inclusion of a further term of the form ΓH dip = Γn 1 n 2 , accounting for an interaction between bosons at different lattice produces a Hamiltonian of the same class as H [25]. Note indeed that Also notice that a simple mapping exists between the spectral features of Hamiltonian (1) for attractive U < 0 and repulsive U > 0 interactions.
This can be appreciated by ob- is a unitary transform, [26].
Hamiltonian (1) can be studied by assuming that the state of the system is well approximated by a trial wavefunction of the form where −1 ≤ z ≤ 1 and ϕ are two real dynamical variables describing the population imbalance and the macroscopic phase difference between the two wells [3,27]. Eq. (5) is a su(2) coherent state, as discussed in Ref. [18]. The dynamics of the above variables is determined by the expectation value of the trial state on Eq. (1), which plays the role of a semiclassical Hamiltonian, where γ = UN Ω and ε = v Ω are the effective parameters of the model [2]. Specifically, the dynamics is dictated by the so-called Discrete Self-Trapping equations [10] ensuing from Eq. (6), and the corresponding fixed-point equations are known to exhibit bifurcations at |γ| 2/3 − 1 = |ε| 2/3 [19,28]. The manifestation of the semiclassical bifurcations in the features of the quantum system has been repeatedly reported in the literature. For instance the emergence of structures in the spectrum of the quantum problem (1) at the same energies as the semiclassical fixed-point solutions has been investigated in Refs. [14,26,[28][29][30] In the following we will be mainly concerned with the mirror-symmetric case, ε = 0, for which the bifurcation corresponds to a dynamical transition at γ = −1 in the ground-state of Hamiltonian (6). For γ ≥ −1 this is symmetric, z = 0, while for γ < −1, due to the nonlinear nature of the dynamical equations, the mirror symmetry of the Hamiltonian is spontaneously broken, Although the ground-state of the quantum Hamiltonian (1) is always symmetric, n 1 = n 2 , a crossover in its internal structure can be recognized in the vicinity of the mean-field critical value. Several indicators have been employed to highlight this crossover: energy gaps [14,31,32], fluctuations in the population imbalance [15,16,19,20,23], structure of the occupation number distribution [15,16,21,23,33] fidelity [31], entanglement entropy [23,32,[34][35][36], coherence visibility [16,23], generalized purity [37]. Refs. [17,18,22,24,31] consider similar issues on systems comprising more than two sites.
We mention that Hamiltonian (1) can be mapped onto a Lipkin-Meshkov-Glick (LMG) model with vanishing anisotropy parameter. Some of the results we illustrate in the following have been obtained in that context making use of numerical or analitical techniques different from ours [38][39][40][41][42].

III. THE SCHRÖDINGER-LIKE EQUATION FOR THE LOW-LYING EIGENSTATES
Under suitable conditions, the eigenvalue equation for Hamiltonian (1), H|Ψ = E|Ψ is very satisfactorily approximated by a Schrödinger-like equation for a fictitious particle on a one-dimensional domain. Several slightly different versions of this approach have been employed the literature [5,19,20,25,26,36,[43][44][45][46]. All of them are based on the expansion of the eigenstate on the Fock basis of the site occupation numbers where |0 is the vacuum state, a j |0 = 0. Introducing the normalized population imbalance of the Fock state relevant to expansion coefficient c ν the eigenvalue equation reads where the "rescaled energy" differs from the original eigenvalue by a unimportant population-dependent constant. Now, for large boson populations N , both the square roots appearing in in Eq. (9) and the coefficients c ν±1 can be expanded in powers of N −1 . Note indeed that the latter can be seen as the values that a continuous function ψ takes on at the points of the grid mesh in the interval [−1, 1] defined by z ν with ν = 0, 1, 2, · · · , N . Specifically, and retain only the lowest order terms in the expansion of the interaction and kinetic term in the eigenvalue equation, we obtain the following Schrödinger-like equation [20]: where z ∈ [−1, 1] is the continuum limit of the mesh grid and the effective potential is [19] V ε (z) = γ 2 Owing to the N −2 coefficient in front of the derivative part in Eq. (10), for large populations the low-lying solutions of such equation will be strongly localized in the vicinity of the minima of the potential energy, Eq. (11). This suggest that the Schrödinger-like equation (10) can be further simplified by introducing the (small) deviation from the minima of Eq. (11). Note that the minimization of Eq. (11) is equivalent to the determination of the lowest-energy stationary solution of the Discrete Self-Trapping equations arising from the mean-field Hamiltonian (6) [24]. This is easily verified in the mirror symmetric case, ε = 0, where, up to the fourth order in the small displacement from the minimum point(s), u = z − z γ , equation (10) is equivalent to with ψ(z) = ϕ(z − z γ ) and

IV. HARMONIC REGIME
We note that the leading order in Eq. (12) is almost everywhere the second, except at the mean-field critical point γ = −1 where the first non-vanishing term is the fourth. Therefore, for the lowest energy levels, away from criticality and for sufficiently large populations, we can let B = C = 0 and Eq. (12) reduces to the Schrödinger equation for a harmonic oscillator, whose discrete spectrum is where Actually, for γ < −1 Eq. (12) describes the small oscillations about either of the two equivalent minima of the effective potential, z γ = ± 1 + γ − 2. Therefore a generic solution of Eq. (10) is well approximated by a superposition of one harmonic solution for each well. The symmetry of the problem dictates that these super-positions have a definite parity, Also, since in this approximation the two localized functions ϕ n (z ± z γ ) have a vanishing overlap, both ψ ± n (z) have an energy very close to that in Eq. (14). In particular, the ground-state of the system ψ 0 (z) is well approximated by a symmetric superposition ψ + 0 (z) of two Gaussian functions of the form with σ 2 γ given in (the lower of) Eq. (15). For γ > −1 the effective potential has a single minimum at z γ = 0, and the problem can be mapped exactly onto a quantum harmonic oscillator, as far as the width of the eigenfunctions of the latter problem do not exceed the interval in which the harmonic term dominates over higher-order terms. The energy spectrum is again given by Eq. (14) and, in particular, the ground-state is ψ 0 (z) = ϕ 0 (z), Eq. (17), with σ 2 γ given in (the upper of) Eq. (15).
The most natural physical quantity which is in principle measurable, is given by the population imbalance corresponding to the operator [20] in the original two- Its value and fluctuations are easily evaluated on the ground state of the system ẑ = ds z|ψ 0 (z)| 2 = 0 (19) and [48] where · denotes the expectation value on the ground state of the system. Moreover one can also consider the expectation value · n on the nth excited state obtaining: The fluctuations of the population imbalance have been investigated by several authors, [15,16,19,23], and analytical expressions for this quantity have been reported in Refs. [15,19,23]. In both Gaussian regimes the (low lying) energy spectrum of Eq. (12) is well approximated by Eq. (14).
, andĒ n depends on E n as described in table (13), it is easy to calculate the lowest energy gaps in the spectrum of H, ∆E n = E n − E n−1 . As we mentioned earlier, for  22) and (23). Specifically, the dashed line is Ω γ 2 − 1 and the solid line is Ω √ γ + 1. Symbols represent numerical data. Specifically, circles, upward triangles, squares and downward triangles are ∆En with n = 1, 2, 3, 4, respectively. and In Fig. 1 we compare Eqs. (22) and (23) with numerical data for ∆E n . The plots clearly show that the analytic approximations become more accurate with with increasing boson population.

V. LIMITS OF VALIDITY AND REGIMES
In the present section we describe the limits of validity of the harmonic solutions described above and the different regimes taking over at different values of the effective parameter γ. A crucial assumption in passing from Eq. (9) to Eq. (10) is that the coefficients c ν can be regarded as continuous functions of the population imbalance of the Fock state defined in Eq. (8). This may be problematic for eigenstates in the mid-spectrum of Hamiltonian (1), but is very reasonable for states close to the extrema of the same spectrum. However, the continuous approach can prove ill posed even for the ground state of the system. It is clear that, for the continuous limit of Eq. (10) to be well defined, the Gaussin width σ γ must be much larger than distance between two neighboring points in the discrete Eq. (9), namely 2 N −1 . This means γ ≫ γ 1 ∼ − √ N for γ < −1 and |γ| ≪ γ 4 ∼ N 2 for γ > −1 [25].
We observe that the same estimates can be obtained by imposing that the expansion parameter in the strong coupling regimes be actually small, as discussed in some detail in Sec. A. This ensures that for γ ≫ γ 1 and γ ≪ γ 4 the first-order perturbative description of the ground state applies.
The parameter range γ 1 ≪ γ ≪ γ 4 where the continuous approach applies can be divided into three regimes. For γ 1 ≪ γ ≪ γ 2 the ground state of the system is well approximated by a symmetric superposition of two Gaussians. On the other hand, for γ 3 ≪ γ ≪ γ 4 the ground state is well approximated by a single Gaussian.
These two regimes are separated by an intermediate critical regime γ 2 < γ < γ 3 , surrounding the meanfield critical point γ c = −1, where the higher order terms in Eq. (12) cannot be neglected and the problem is not Gaussian any more. The boundaries of this regime can be obtained from a comparison between the leading and next-to-leading term in the perturbative expansion of the effective potential. A reasonable estimate for the lower bound γ 2 < γ c can be obtained either by imposing A σ 2 γ ≫ |B|σ 3 γ or by requiring that the distance between the two equivalent harmonic minima, 2 1 − γ −2 , is much larger than the width of the relevant Gaussian ground states σ γ . In both cases one finds γ c − γ 2 ∼ N −2/3 . Similarly, the upper bound γ 3 > γ c is obtained by requiring A σ 2 γ ≫ C σ 4 γ , which results once again into γ 3 − γ c ∼ N −2/3 [19,45]. As discussed in Sec. VI, the energy gaps in the critical regime depend on the system population, ∆E ∼ N −1/3 , at variance with what happens in the Gaussian regimes. This provides a N −dependent cutoff to the vanishing of ∆E at γ = γ c predicted in Eqs. (22) and (23).
In summary, one can recognize five different regimes depending on the value of the effective parameter γ. For γ ≪ γ 1 and γ ≫ γ 4 the continuous approximation of Section III does not apply, but the system can be analyzed by resorting to first-order perturbative theory, as discussed in Appendix A. For γ 1 ≪ γ ≪ γ 2 the continuous approximation applies, and the ground-state of the system is a symmetric superposition of two Gaussians centered at the equivalent harmonic minima of the effective potential. For γ 2 ≪ γ ≪ γ 3 the system is critical, in a sense that will be clarified shortly. Note that in the limit of large populations the critical regime reduces to the mean-field critical point. For γ 3 ≪ γ ≪ γ 4 the groundstate is a Gaussian centered at the unique minimum of the effective potential, z γ = 0. Note that in estimating the location of the boundaries γ 1−4 we only considered the dependence on the (large) bosonic population, and neglected other numerical factors.
VI. CRITICAL REGIME According to our previous analysis, as the effective parameter γ approaches the mean-field critical value γ c = −1, the system exhibits the hallmarks of a phase transition, where the minimum z γ of the effective potential, table (13), plays the role of the order parameter. Note that this quantity is in principle measurable as the absolute value of the population imbalance operator z γ = |ẑ| . The fluctuations of this order parameter coincide with the Gaussian width of Eq. (15).
The behavior for z γ and σ 2 γ described in table (13) and Eq. (24), respectively, holds true in the mirror-symmetric case, v = ε = 0. It proves useful to study behaviour of the order parameter in the presence of a vanishing energy offset between the sites of the dimer, |v| ≪ 1. According to our discussions about the Gaussian regimes, the value of the order parameter will coincide with the absolute minimum of the effective potential which, owing to the presence of a non-vanishing ε, is now non degenerate for any value of γ. The minimization of Eq. (11) then gives where only linear terms in ε have been taken into account. Eq. (25) allows us to introduce a sort of susceptibility Equations (13), (15), (24), (26) provide some critical exponents for the system. As γ approaches its critical value γ c = −1, the order parameter vanishes (from below) as |ẑ| ∼ (γ c − γ) αz , its fluctuations diverge as |ẑ| 2 − |ẑ| 2 ∼ |γ − γ c | −ασ and its susceptibility diverges as χ z ∼ |γ − γ c | −αχ , where The above discussion rigorously applies in the limit of infinite population, which plays the role of an effective thermodynamic limit. Indeed, as we observe in Section V, only in this limit does the critical regime shrink to a single point corresponding to the critical value. For any finite populations there exists a finite region surrounding the critical point where the Gaussian results do not apply. In this region the quartic term dominates and Eq. (10) can be approximated as where C is listed in the upper row of table (13). It is easy to prove that in this critical regime the eigenfunctions of Eq. (28) have form where φ n (ζ) is the solution of the populationindependent problem and the eigenvalue corresponding to ψ n (z) is E n = N 2/3Ē n . This result allows us to estimate the behaviour of the (low-lying) energy gaps in the spectrum of Hamiltonian (1) in the critical regime. We get ∆E n ∼ N −1/3 , where we used the definition of E in table (13) and the fact that the spectrum of problem (30) does not depend on N . This estimate, whose correctness demonstrated in Fig. 2, provides an alternate route for the determination of the boundaries of the critical region. Indeed, these can be assessed by requiring that the prediction in the Gaussian regimes, Eqs. (22) and (23), are of the same order as that in the critical regime, i.e. N 1/3 . This gives once again γ c − γ 2 ∼ N −2/3 and γ 3 − γ c ∼ N −2/3 . As to the fluctuations of the order parameter, particularizing Eq. (29) to the ground-state of the system, n = 0, we get Within a standard scaling approach it is reasonable to assume that there exists a correlation length, dictating the effective size of the system, which diverges at the critical point Any physical quantity should depend on the boson population N only through the ratio N/N γ , so that it should be The comparison with Eqs. (15), (24) and (31) then allows us to conclude that The correctness of the scaling assumption in Eqs. (33)-(34) is demonstrated in Fig. 3, where a nice data collapse is observed for a wide range of boson populations. We mention that -in the context of the LMG model -the critical exponents in (34) were determined numerically in Ref. [38] and making use of an analytic approach different from ours in Refs. [39][40][41].
We conclude by remarking that the rescaled energies E n and the expectation values for the rescaled problem (30) can be estimated by applying the Bohr-Sommerfeld quantization rule. This gives which provides quite satisfactory estimates of the energy levels of the system, where C 1 is a suitable constant. Within the same approximation the mean square displacement of the eigenfunction φ n (ζ) can be evaluated as well. In particular, for the simclassical orbit of energyĒ n , one gets φ n (ζ) 2 ζ 2 dζ ≈ ζ 2 δ(p 2 + ζ 4 16 −Ē n ) dp dζ where C 2 is a suitable constant.

VII. QIT INDICATORS
As we mention above, quantities borrowed from quantum information theory have been employed to investigate the the incipient quantum phase transition between the localized and delocalized ground-state of the system. Owing to the particularly simple structure of the system, some of these quantities turn out to be entirely equivalent to usual indicators. This is the case of the generalized purity of Ref. [37], which basically coincides with the coherence visibility considered in Refs. [16,23,36], i.e. the off-diagonal term of the one-body density matrix. Similarly, the Fisher information of Ref. [23] is equivalent to the population imbalance fluctuations considered in several earlier works [15,16,19,20] as well as in the present paper.
The fidelity and entanglement entropy considered in Refs. [31] and [23], respectively, appear to be less trivial quantities [49]. The former is nothing but the overlap between two ground-states of the system corresponding to slightly different values of the parameter driving the transition Ψ γ |Ψ γ+κ [50]. Since this quantity is almost everywhere extremely close to unity, a more convenient indicator is provided by the so-called fidelity susceptibility [51] In a bipartite system, the entropy of entanglement of a pure state |Ψ is given by −Tr(ρ j log ρ j ), where ρ j is the reduced density matrix of subsystem j, obtained from the full density matrix |Ψ Ψ| after tracing over the states of the other subsystem [52]. In the case under investigation the two subsystem are of course the two sites composing the Bose-Hubbard dimer and, owing to number conservation, the entropy of entanglement of an eigenstate of Hamiltonian (1) is where the c ν 's are the expansion coefficients in Eq. (7) and we introduced a normalization factor. This quantity has been considered in a few earlier works on the Bose-Hubbard dimer. In the presence of repulsive interactions it was used to investigate the precursor of the Mott insulator-superfluid quantum phase transition occurring on infinite lattices [53]. In attractive systems it was employed for the characterization of Schrödingercat-like states in the presence of a small energy offset v between the two sites of the system [36], and in the study of the quantum counterpart of the mean-field transition for small finite boson populations [23]. These two indicators have been considered for the LMG model [42,54] which, as we mention, can be mapped onto Hamiltonian (1). However, the natural decomposition of the system necessary for the calculation of the entanglement entropy is quite different in the spin and boson context, and the results are not readily comparable. As to the ground-state fidelity, a comparison will be made with the results in Ref. [42].
Based on the results sketched in the previous sections and in the appendices, we are able to give an analytic description of the behaviour of the fidelity susceptibility and entropy of entanglement, Eqs. (38) and (39).
Let us start with the former quantity. As we discuss in Section III, in the Gaussian regimes the ground-state of the system is strictly related to the ground state of a harmonic oscillator. Specifically, for γ 1 ≪ γ ≪ γ 2 the ground state is well approximated by a symmetric superposition , ψ + 0 (z), of two Gaussian functions of square width σ 2 γ = N |γ| γ 2 − 1 −1 centered at z γ = ± 1 − γ −2 , as described by Eqs. (16) and (17). Likewise, for γ 3 llγ ≪ γ 4 the ground state is well described by a single Gaussian ϕ 0 (z) of square width σ 2 γ = N √ 1 + γ −1 centered at z γ = 0, as described by Eq.

(17). Straightforward calculations result in
whose correctness is demonstrated in Fig. 4. It is worth noticing that Eq. (40) predicts a dependence on the boson population only in the lower Gaussian regime γ 1 ≪ γ ≪ γ 2 , whereas for γ 3 ≪ γ ≪ γ 4 the fidelity susceptibility turns out to be N −independent. This can be explained by the fact that in the latter regime the small variation κ in the effective parameter affects only the width of the single Gaussian representing the groundstate of the system. Conversely, for γ 1 ≪ γ ≪ γ 2 , the change affects both the centers and the widths of the two Gaussians comprising the ground state of the system. Also notice that Eq. (40) predicts a divergence in the fidelity susceptibility at the mean-field critical point, γ = γ c . However, for any finite N , there is a small region surrounding such point where the Gaussian approximations 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 larger the population, the sharper and the closer to the critical point is the peak. In order to study the behavior of χ F (γ) in the critical region it proves useful to exploit Eq. (28). As it is discussed in Ref. [55], the limit in Eq. (38) can be carried out exactly through a perturbative expansion involving the entire spectrum of the problem at a given value of the driving parameter γ. Considering the harmonic term in Eq. (28) as a perturbation over the quartic term, we get where the eigenfunctions φ n (ζ) and eigenvaluesĒ n are dfined by (29) and (30). Rigorously speaking the φ n (ζ) provides a good approximation for an eigenstate of the original problem, Eq. (1), only for sufficiently small values of the quantum number n. However, only the first few terms in the sum appearing in Eq. (41) are expected to provide a significant contribution, owing to the localization of the ground state φ 0 (ζ) and the energy gap appearing in the denominator. The N 4/3 behavior predicted in Eq. (41) and confirmed by the numerical results in Fig. 5 can be recognized as the superextensive divergence of the fidelity susceptibility which, according to Ref. [51], is the hallmark of a quantum phase transition. The same exponent was estimated numerically in Ref. [42]. Also, our Eq. (40) agrees with a similar result in the same paper, although we obtain a different behavior of the sub-leading terms.
Let us now turn to the entropy of entanglement. Recalling that the continuous version of Eq. (39) is it is easy to check that in the Gaussian regimes where Note that, despite Eq. (39) forbids that S e (γ) exceed unity, the quantity in Eq. (43) diverges in the vicinity of the mean-field critical point. However, this is not alarming, since at any finite boson population there is a small interval around such point where the ground-state is not Gaussian and Eq. (43) does not apply. It could be easily proven that the requirement that Eq. (43) be less than unity provides an estimate for the boundaries of the non-Gaussian regime that is equivalent to the one derived in Section V. An interesting feature revealed by Eqs. (43)- (44) is that in the N → ∞ limit the entanglement entropy tends to a constant, independent of the interaction to hopping ratio, S e (γ) = 1 2 ∀ γ = γ c . Note indeed that quantity in Eq. (44) does not depend on N , so that the second term in Eq. (43) vanishes with increasing population.
The behavior of the entropy of entanglement in the critical regime can be obtained by plugging ψ(z) = N 1/6 φ 0 (N 1/3 z) into Eq. (42), where φ 0 (ζ) is once again the ground state of Eq. (30). Straightforward calculations give Since the integral in the second term of Eq. (45) does not depend on N , we get that lim N →∞ S e (γ c ) = 2 3 . Thus the mean-field critical point is signalled by the entropy of entanglement as a peak which, in the large population limit, turns into a discontinuity. a comparison between Eqs. (43)- (44) and numerical results obtained via exact diagonalization. The agreement is quite satisfactory, and improves with increasing boson population. Note that even at relatively large populations the entropy of entanglement is still rather different from its limit. This is because the finite size corrections decrease logarithmically with N . The different behavior at criticality is demonstrated in Fig 7, which shows the correctness of Eq. (45).

VIII. THERMAL FLUCTUATIONS
We now turn to the thermal fluctuations of the order parameterẑ where · T and · n denotes the thermal average and the quantum expectation value on the state of energy E n , respectivley.
Evaluating E n , ẑ 2 n and ẑ n using standard properties of quantum hamonic oscillators (see Section IV) we get where the thermal fluctuations are and z γ is given in Table ( 13). Therefore for k B T < ∆E 2 (i.e. the lowest energy gap of the system) σ 2 T,γ is essentially given by the fluctations of the quantum ground states (20), whereas for large temperatures the thermal fluctuation dominates and we have Eq. (49) highlights the typical linear dependence on temperature of the thermal fluctuations in harmonic systems. We remark that both thermal and quantum fluctuations decrease with the total population as 1/N . This is demonstrated by the data collapse in Fig. 8. Also in the thermal regime we expect the harmonic approximation to apply only when the higher order therms in Eq. (12) are negligible. This means that for γ 1 < γ < γ 2 , σ 2 T,γ should be smaller than A/B while for γ 3 < γ < γ 4 , σ 2 T,γ should be smaller than A/C. In both cases for γ ≈ −1 we get kT N (γ + 1) 2 . Note indeed that, as evident from Fig. 8, the higher the temperature, the larger must the population be for the harmonic approximation (dash-dotted line) to apply.
(50) It is easy to check that also in this critical regime σ γ,T reduces to the purely quantum contribution when k B T is smaller than the energy gap between the lowest eigenpair, C 1 Ω/N 1/3 . At higher temperatures we have which is the typical behaviour of classical thermal fluctuations when a pure quartic potential is present. We remark that within the critical region the fluctuations 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 and lower panel of Fig. 9, respectively. The different quantum and classical regimes around γ = −1 can be summarized as in Figure 10. Note that the boundaries between different regions, signalled by the coloured lines, are actually crossovers between different behviours. The only true critical point is γ = −1, in the limit of infinite N .
In regions I and III σ 2 γ,T coincides with the result applying at zero-temperature for purely quartic and harmonic potentials, respectively. Region II is characterized by thermal fluctuations in a harmonic potential, while in region IV we have a classical behaviour in an anharmonic potential. In regions I and II σ 2 γ,T decays as N −1 , while in region III quantum critical effects produce fluctuations decaying as N −2/3 . Finally, in region IV fluc-  tuations should decay as N −1/2 for not too large temperatures, and they should be independent of N for very large k B T /Ω.

IX. CONCLUSIONS AND PERSPECTIVES
In this paper we addressed the quantum counterpart of the transition from a localized to a delocalized ground-state occurring in an interacting bosonic gas trapped in a double well potential, i.e. a bosonic Josephson junction. We identified and discussed five regimes depending on the value of the effective interaction among bosons. For sufficiently small effective interactions we analyzed the low-lying spectrum of the system Hamiltonian (a twomode Bose-Hubbard model) by recasting the problem in terms of a Schrödinger-like equation for a single particle in a one-dimensional domain. In particular we studied how the results of the analysis on the quantum system reduce to those of the semiclassical treatment based on the Discrete Self-Trapping equations as the bosonic population is increased. We devoted a special attention to the critical regime which, in the theromodynamic limit of infinite bosonic population, collapses to a bifurcation point in the semiclassical picture. We extended our analysis to quantities borrowed from Quantum Information Theory, namely the entropy of entanglement and the ground-state fidelity, and to finite temperature effects. In the latter respect, we evidenced that, in the critical regime, the quantum and thermal fluctuations of the population imbalance of the two wells exhibit different scaling behaviors, which may be amenable to quantitative measurement.
As we mention in Section V, when interaction dominates over kinetic energy the continuous approach of Section III does not apply. However, in this situation an approximation of the ground-state of the system can be given by resorting to perturbative theory. For large repulsive interactions, first order perturbative theory gives where, assuming without loss of generality that N is even, Equation (A1) applies when the factor in front of the first order perturbative term is small, λ ≪ 1. It is easy to check that, for this to be true, it should be γ ≫ γ 4 , where γ 4 ∼ N 2 is the upper bound of the "single Gaussian" regime discussed in Section V. Note that in passing from the "single Gaussian" regime, γ 3 ≪ γ ≪ γ 4 , to the present "strong coupling" regime, γ ≫ γ 4 , there is a change in the behaviour of the fluctuations of the population imbalance operator, Eq. (18), on the ground-state of the system, as discussed in Ref. [25]. Indeed in the former regime this is clearly ẑ 2 − ẑ 2 = σ 2 γ , where σ 2 γ ∼ N −1 is given in Eq. (15). Conversely, in the strong coupling regime of perturbative state (A1), we get where the last approximate equality applies in the limit of small perturbations and large populations. Since for γ ∼ γ 4 the fluctuations in Eqs. (15) and (A4) are of the same order of magnitude, there is no intermediate regime between the "single Gaussian" and the repulsive "strong-coupling" regimes Similar arguments apply in the attractive "strongcoupling" regime. At first order in perturbation theory the ground state is where |Ψ 0 = 1 We observe that, unlike the repulsive case, there is no cross-over of the fluctuations of operator (18) in passing from the strong-coupling regime γ ≪ γ 1 to the "double Gaussian" regime described in Section V. Note indeed that Eq. (A6) coincides with Eq. (20) at large negative γ. The same results apply when one considers a single peak, which gives But this is exactly the same behavior as σ 2 γ in Eq. (15) for large populations and large (negative) γ.