We develop a variational wave function for the ground state of a one-dimensional bosonic lattice gas. The variational theory is initally developed for the quantum rotor model and later on extended to the Bose-Hubbard model. This theory is compared with quasi-exact numerical results obtained by Density Matrix Renormalization Group (DMRG) studies and with results from other analytical approximations. Our approach accurately gives local properties for strong and weak interactions, and it also describes the crossover from the superfluid phase to the Mott-insulator phase.
©2004 Optical Society of America
During the last years a spectacular development in the storage and manipulation of cold atoms in optical lattices [1, 2] has taken place. Greiner et al. , to name one important example, succeeded in experimentally driving a quantum phase transition between a superfluid and a Mott-insulating phase in bosonic systems. This experimental progress has revived the interest in the Bose-Hubbard model [Eq. (1)] as a generic Hamiltonian for strongly correlated bosons, by which the quantum phase transition can be described [3, 4]. The Bose-Hubbard Hamiltonian has been used previously in condensed matter physics to study the adsorption of noble gases in nanotubes , or Cooper pairs in superconducting films with strong charging effects [6, 7]. In this context a lot of work has already been done to characterize the quantum phase transition, the statistics, and the low-energy excitations of the Bose-Hubbbard model [8, 9]. However, new interesting questions arise now due to the good tunability of the experiments with optical lattices. In particular, it becomes possible to study time-dependent processes such as driven quantum phase transitions . A theoretical understanding of such phenomena is challenging, since the characteristics of the superfluid phase—where the atoms tend to delocalize throughout the lattice and large fluctuations in the local density exist —, and the Mott-insulating phase — where the number fluctuations decrease, and a gap in the excitation spectrum opens —, must be covered at the same time. Both regions are separated by a non-analyticity of the spectrum, which implies that a perturbative study [10, 11] works best in strong coupling limit, while a Hartree-Fock-Bogoliubov mean field works best in the superfluid regime. In addition it is possible to develop a mean field theory [12, 13] based on a Gutzwiller ansatz [3, 4]: this reproduces the mean field theory in the superfluid limit, as well as the limit of infinite interaction, which raises the hope that the theory also interpolates properly between these limits.
In this paper we develop a variational description of the ground state of an ensemble of cold atoms in an one-dimensional optical lattice. Our trial wavefunction treats the connections between neighboring sites as entities which decouple in the limit of infinitely large lattices. We apply this technique first to the quantum rotor model —which describes the lattice for large and commensurate occupation per site , and has also been used to describe an array of Josephson junctions —, and next to the Bose-Hubbard Hamiltonian. The accuracy of the variational theory in both models is confirmed by several comparisons. In the quantum rotor case, we use a spin wave approximation in the limit of weak interaction and a first order perturbation theory in the limit of weak tunneling. For the Bose-Hubbard model we compare against results obtained applying the quasi-exact, numerical DMRG method to one-dimensional lattices with up to 128 sites, and also with calculations based on the Gutzwiller ansatz. Our conclusion is that the variational picture of self-regulated connections between sites provides a rather cheap and simple, but very good description of the local properties of the system in the superfluid and insulator regimes, and a fairly good interpolation across the quantum phase transition. Our variational ansatz, however, fails to describe long range properties of the superfluid phase, such as the algebraic decay with distance of the off-diagonal elements of the one-particle density matrix. The simplicity of the method suggests a possible generalization to higher dimensionalities and other physical models.
The outline of the paper is as follows: In Section 2 we introduce the quantum rotor model as a possible limit of the Bose-Hubbard Hamiltonian. Next, information about the ground state of the quantum rotor model is obtained variationally as the solution of a Mathieu equation. We can estimate energies, correlation functions and length, and the variance of the density as a function of the only free parameter. A comparison with perturbative estimates demonstrates the accuracy of the method when computing local properties. Since the quantum rotor model is only an approximate description of the optical lattice, in Section 3 we develop a similar variational theory for the Bose-Hubbard Hamiltonian. After bringing the Bose-Hubbard Hamiltonian to an appropriate form, we can estimate the local properties of its ground state. The variational solutions are compared in Sec. 3.3 with the results of DMRG studies of the Bose-Hubbard model. We confirm that the variational method describes very well the local properties of both the Mott-insulator and the superfluid regime, and provides a fairly good interpolation across the phase transition. Finally, in Section 4 we summarize our results.
2. Quantum phase model
2.1. Relation to the Bose-Hubbard model
In this section we show the equivalence of the Bose-Hubbard model
and the quantum rotor model for large and integer occupation n̄. In our notation, M is the number of lattice sites and N=n̄M the number of atoms. Both the Bose-Hubbard model and the quantum rotor model show a phase transition due to the interplay between the kinetic term proportional to J and the interaction term proportional to U. For convenience we have subtracted the ground state energy in the perfect insulator limit U/J→∞.
If we expand a configuration of the lattice using Fock states
and the number of particles per lattice site is large, nk >1, we may approximate the hopping terms as follows
Here, A±|n〉=|n±1〉 are ladder operators and 𝓟 projects on states with non-negative occupation numbers, nk ≥0. To lowest order the error |Δlj 〉 is
and its norm is bound by
where =〈(nl -n̄)2〉 is the variance in the number of particles per lattice site. For the approximation (3) to be valid, the uncertainty in the number of atoms must be small compared to the mean value, n̄≫σ, and the interaction energy must exceed the neglected terms, Ūn(n̄-1)≫Jσ.
Following the previous procedure the Bose-Hubbard model becomes
Here is approximately the density, =a†jaj-n̄ is essentially the number operator and we have used that ∑j |ψ〉=0 when we work with states that have a fixed, commensurate number of particles. In the following we will define the energy per lattice site as
Since the physically interesting states will be concentrated around large occupations, nk =n̄, the usual step now is to drop the projector, 𝓟, and move to the basis of phase states, defined by
In doing so, we obtain the identification and →-i∂ /∂ϕ j , which produces the usual representation of the quantum rotor model
with the associated state writen as
A similar derivation is possible using path integrals .
2.2. Variational ansatz
In this section we estimate the properties of the ground state of HQR variationally. Due to the previous splitting (10), any wavefunction Ψ(ϕ⃗) can only depend on the phase difference between neighboring wells, ξj =ϕ j+1-ϕj. These new quantum variables describe the connections between neigboring sites. In the limit of large lattices it seems reasonable to assume that these connections become independent from each other adopting the product state
This representation becomes exact in the Mott-insulating limit, U/J→∞, where hmott (ξ)=1, and in the superfluid limit, U/J→0, where hsf (ξ)=∑n∈ℤδ(ξ-2πn), as can be verified by direct substitution in Eq. (10).
Even though the phase representation is the best one to find a trial wavefunction, it is not the optimal one for performing computations. It is instead more convenient to work with the variables which are conjugate to the phase differences ξj=ϕ j+1-ϕj . These are the new quantum numbers, wk , given by
In terms of these numbers, the ansatz (11) reveals itself as a simple product wavefunction
with coefficients given by the Fourier transform
As sketched in Fig. 1, the wk play the roles of chemical potentials which are established between different wells: the difference between the potentials on the extremes of a site gives the fluctuations over the mean and commesurate occupation n̄. In this picture
where ∑±|w〉=|w±1〉 are new ladder operators and ∑ z |w〉=w|w〉.
By minimizing the energy associated with HQR over all states within a given ansatz we can both obtain an upper bound to the energy of the ground state and approximate its wave function. A simple computation with our product ansatz leads to the result [compare (7)]
where the expected values are computed over a single connection, 〈∑ z 〉=∑ww |hw |2, and the wavefunctions are assumed to be normalized, ∑ w |hw |2=1. Since the stationary states have a well defined parity, h̃(-w)=(-1) Ph̃w , the optimal variational state must satisfy the linear equation
which is nothing but the Fourier transform of a Mathieu equation
The estimated ground state energy per site is given by the lowest eigenvalue of either equation.
Using the product ansatz and the same approximations required to derive HQR , we can also compute other properties of the ground state. For instance, the variance of the number of atoms per lattice site
and the correlation functions
which decay exponentially with the distance. This implies that the ansatz (13) only describes properly the decay of the correlations in the Mott-insulating regime, since the correlations in the superfluid regime follow a power law decay. However, as we will see below, local observables (σ, γ1, ε …) are properly estimated even if long–range ones are not.
We have solved Eq. (17) numerically in a truncated space. The results are summarized in Fig. 2, where we also plot reference estimates arising from two other analytical methods. In the limit U≫J we compare with a first order perturbative calculation around the solution Ψ(ϕ⃗)=1, which is possible thanks to an energy gap of order 𝒪(U) in the excitation spectrum. In the limit J≫U we rather use a spin wave or harmonic approximation in which the cosineterms of the Hamiltonian HQR are expanded up to second order in the phase difference between neighboring sites [see Sec. 2.3]. This approximation is valid in the superfluid regime, where the phase does not vary much between neighboring wells. From the graphical comparison we see that the variational wavefunction provides a fairly accurate description of the ground state of the quantum rotor model in both the superfluid and insulating limits. As a side note, we must remark that this ground state has a divergent fluctuation of the number of particles per site as J→0. This differs from the expected behavior of the ground state of the original Bose-Hubard Hamiltonian, and it reminds us that HQR can only model the atomic ensemble when the variance, σ, is small compared to the mean occupation number, n̄.
2.3. Harmonic approximations to the quantum phase model
In the limit J≫U it is possible to estimate the ground state of the rotor model (6) analytically. Since we are deep in the superfluid regime, the wavefunction will be concentrated around the line ϕ1=ϕ2=…=ϕM, and we can approximate
If we remove the periodic boundary conditions on ϕi and change variables, the preceding Hamiltonian may be diagonalized, , with frequencies given by ωk =(8ρJU)1/2 |sin(πk/M)|, where k is an integer in the range -M+1<2k<M and labels the different values of the momentum in the lattice. The ground state energy [Fig. 4(a)] may be estimated as the zero-point energy of the harmonic oscillator. For large M, the sum over k may be replaced with an integral, giving
The variance of the number of particles is related to the expectation value of the momentum, and using the same procedure as above one obtains
3. The Bose-Hubbard model
In this section we apply to the Bose-Hubbard Hamiltonian the techniques that were developed in Sec. 2. We will do it in three steps: First we will develop a phase representation which is valid for all occupation numbers. Next we will prove that this representation is equivalent to a similarity transformation of the Hamiltonian which brings it to a form similar to (6), at the price of losing Hermiticity. Finally we will show how to implement the ansatz of independent connections (13) to produce estimates for the usual set of observables (ε, γ1, σ), which are to be validated with DMRG calculations.
3.1. Coherent states
The phase coherent states |ϕ〉 are defined by
Unlike the phase states defined in Sec. 2, they are not orthogonal to each other, 〈ϕ|θ〉=exp[ei(θ-ϕ) ], but they form a complete basis and an expansion like (10) is still possible. A nice property of the coherent states is that we can rewrite the operators a, a †, a † a, etc, in terms of the phases very easily. For instance, , and .
Using this representation, we obtain an effective Hamiltonian for the wavefunction Ψ(ϕ⃗), i.e. , which has the form
Here stands for the tranpose of Hcoh . This operator was already used in Ref.  to study the Bose-Hubbard model with only two sites. On the one hand, it is a non-Hermitian operator1 and we cannot do a simple variational study. On the other hand the Hamiltonian still depends on the phase differences, and it is reasonable to look for approximate eigenstates which have the form (11). This will be done in the following section.
3.2. Variational procedure for non-Hermitian operators
In this section we will find the best variational function which has the product form of Eq. (11). However, as it happened in Section 2, instead of working with phase variables it will be more convenient to develop a representation of the Bose-Hubbard Hamiltonian in terms of connections. This is once more a two-steps process. First we use a similarity transformation suggested by the definition of the coherent states
Since OajO -1=A-j and Oa†jO-1=A+j , we find
The Hamiltonians (26) and (28) are equivalent: while one is defined in terms of phase variables, the other one is defined using occupations numbers, and both are related by a Fourier transform. The second and final step is to rewrite everything in terms of connections, with the quantum numbers from Eq. (12) and the relations ∑ x =∑++∑-, ∑y=i(∑--∑+). The result is a decomposition of the Hamiltonian
into terms which are local, H 1, and terms which involve pairs of connections, H 2.
For the quantum rotor model we proved that the optimal product wavefunction was an eigenstate of a Hamiltonian which did not couple connections, like H 1. The difference now is that, since the operator Hcoh is not Hermitian, we cannot establish a variational principle and that argument is no longer valid. Nevertheless, we will again propose a variational ansatz which is an eigenstate of the local operator H 1|h̃〉⊗ M =Mεest |h̃〉⊗ N . Using the following equality
and the product ansatz |χ〉=|h̃〉⊗ M , we arrive to an upper bound for the lowest eigenvalue of the Bose-Hubbard Hamiltonian, expressed in terms of the non-Hermitian one
1The hermiticity of HBH is maintained due to an implicit projection that takes place when we reconstruct the state HBH|ψ〉 from HtcohΨ(ϕ⃗) (See Ref. ).
The way to use this variational principle is as follows. First, for a given J and U we compute the lowest eigenstate of H 1 and this way obtain h̃. After the equivalence (14), finding the ground state of the local Hamiltonian H 1 becomes equivalent to solving a modified Mathieu equation
which describes exactly the static properties of a pair of sites with open boundary conditions . Once we have εest we must still compute the correction Δεest using a rather straightforward expansion which is shown in Sec. 3.4. Surprisingly, Δεest happens to be negative, so that it is actually an improvement over the simple estimate given by εest [See Fig. 5].
From the optimal variational state, |ψ〉=O-1|h̃〉⊗ M , and the estimate for the energy, εvar =εest +Δεest , we may compute other observables. For the density fluctuations and nearest neighbor correlations we use the virial theorem
whereas for other properties one has to evaluate numerically the matrix products shown in Sec. 3.4. This allows us to prove that for the product states |h̃〉⊗M correlations decay exponentially, opposite to what is expected in the superfluid phase, whose correlations should decay algebraically. Nevertheless, as we will see next, this family of states does estimate accurately the local properties of the optical lattice.
3.3. Comparison to DMRG results
We will now compare the results for the ground state energy, the correlation functions, and the variance of the particle number provided by the two variational ansatz, (18) and (32), and the Gutzwiller ansatz  with those obtained by DMRG studies of the Bose-Hubbard model. The DMRG, developed 1992 by White [17, 18] in the area of condensed matter theory, is a very powerful numerical tool to investigate static and dynamic properties of strongly correlated quasi-one-dimensional spin, fermionic or bosonic quantum systems. The DMRG is an essentially quasi-exact numerical method. The fundamental ideas stem from real space renormalization methods: the system size is grown iteratively while the (exponentially diverging) size of the Hilbert space is kept constant by decimation. Hereby one tries to retain only that subset of states that is essential to describe the physical quantity under consideration. In DMRG these are expectation values with respect to low-lying states (“target states”), and in particular with respect to the ground state wave function.
DMRG builds up the system linearly: at each growth step, suitable density matrices for the target states are derived that yield information on the relevance of Hilbert space states. Building on this information, the states and operators are projected onto Hilbert subspaces of fixed dimension M containing the most relevant states. M is chosen to be small enough to be handled numerically, but large enough to obtain the desired accuracy; numerical results can be extrapolated in M to the exact limit of infinite M in the thermodynamic limit. However, results presented here have converged for the largest M considered and no further extrapolation was necessary.
Details on the DMRG method may for example be found in . In the case of the Bose-Hubbard model the DMRG has been used to study properties of the system [20, 21, 22]. The results of DMRG agree very well with exact diagonalization results for small systems, with quantum Monte-Carlo simulations e.g. [23, 24, 25, 26], and with 13th order perturbation theory .
We have used the DMRG to study the properties of the ground state of the Bose-Hubbard model on one-dimensional lattices with 128 sites, and commensurate fillings n̄=1,2, and 3. In Figs. 3(a–c) we show the mean energy per site ε, the nearest neighbor correlations, c 1=〈aj 〉, and the variance σ of the density, calculated both with the DMRG and with the variational estimates developed above. As expected, there are no indications of the phase transitions in these quantities, neither in the variational results nor in the numerical solutions. Rather, an inflexion of the nearest neighbor correlation points out the location of the superfluid-insulator transition which lies roughly between 3n̄ and 4n̄ (see  and ref. therein). The agreement of the two methods is fairly good above the phase transition and below it.
A more detailed comparison is provided in Fig. 4 for the case n̄=1. In this figure we plot together results from the DMRG, the variational ansatz derived above, the quantum rotor model and the well-known Gutzwiller ansatz. The Gutzwiller ansatz [28, 3, 4] is a variational ansatz made of a product of single-site wave functions, |Ψ G 〉=ΠMj=1 |Φj〉, where |Φj〈=∑∞m=0 f (j)|mj 〉 and are constants. Such a wavefunction cannot be used in the one-dimensional Mott insulator regime, because a perturbative study shows that the corrections of order 𝒪(J/U) are lost and all correlations become zero. However this ansatz gives good results in the superfluid regime, where the long-range order is well described by |Ψ G 〉, and we can use these results and those of the DMRG to assert the accuracy of our variational estimates. As Fig. 4 shows, in the Mott insulator regime, the DMRG results agree perfectly with our variational theory for the Bose-Hubbard model and for the quantum rotor model. Close to the phase transition is the point at which the quantum rotor model no longer describes well the atoms in the optical lattice due to the growth of density fluctuations. At this point we also observe a small disagreement between the DMRG and the coherent states, which is due to the growth of long range correlations and vanishes as we get deeper into the superfluid regime.
3.4. Numerical evaluation of the upper bound
In this section we will show how to compute the corrections to the local energy, Δεest , from Eq. (31). We basically need a method to compute expectation values of the operator O -2 around product states which have the form
in which at most two contiguous vectors are affected by single-connection operators. For instance, this is the case of the optimal variational state, |ψ〉=O -1|h̃〉⊗ M , where A and B are just the identity, but also of H 2|ψ〉, where A and B are ∑+, ∑- or ∑z. After some manipulations it is possible to write
with the real matrices and vectors
We have used this technique to compute numerically the correction Δεest for different lattice sizes and found small or no differences for more than 30 sites. Intuitively, this is because in the limit of large powers the matrices (HO) k become projectors on the eigenvector with the largest eigenvalue. Using the same type of expansion we may compute other correlations
For large values of Δ and large lattices, this quotient will decay exponentially as γΔ1 , where γ 1 is the largest eigenvalue of the matrix H_O.
In this work we have studied analytically and numerically the properties of the ground state of an ensemble of bosonic atoms in an 1D optical lattice. For the study of the atomic ensemble we have used both the quantum rotor model and the Bose-Hubbard model. Exploiting the fact that in these models there exists only nearest neighbor hopping and local interactions, we have developed a variational wavefunction that may be used to easily estimate local properties, such as the energy per well, the nearest neigbor correlations and the fluctuations of the density. In the case of the quantum rotor model we have verified our results with perturbative calculations around the strongly interacting regime, and with a spin wave approximation around the superfluid regime. In the case of the Bose-Hubbard model we have compared the variational estimates with numerical results obtained using the DMRG technique for a maximum density of three atoms per well. We have concluded that this procedure leads to fairly good estimates of local ground state properties of both Hamiltonians, in both the superfluid and the insulator regime, the largest disagreement being localized around the phase transition. On the other hand, we have also shown that our variational ansatz fails to describe long range properties of the superfluid phase, such as the algebraic decay of the first order correlation function with respect to distance.
We thank W. Zwerger for fruitful discussions. J. J. García-Ripoll and J. I. Cirac thank the Deutsche Forschungsgemeinschaft (SFB 631) and the Kompetenznetzwerk Quanteninformationsverarbeitung der Bayerischen Staatsregierung. C. K. and U. S. thank the Hess-Preis of the DFG and the Studienstiftung des deutschen Volkes for financial support. Work at the University of Innsbruck is supported by the Austrian Science Foundation, EU Networks and the Institute for Quantum Information. P.Z. thanks the Max Planck Institut für Quantenoptik for hospitality during his stay in Garching, and thanks the Humboldt Foundation for support during this stay.
References and links
1. M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, “Quantum phase transition from a superfluid to a Mott insulator in a gas of ultracold atoms,” Nature 415, 39 (2002). [CrossRef] [PubMed]
3. D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, “Cold Bosonic Atoms in Optical Lattices,” Phys. Rev. Lett. 81, 3108–311 (1998); [CrossRef]
4. D. Jaksch, V. Venturi, J. I. Cirac, C. J. Williams, and P. Zoller, “Creation of a Molecular Condensate by Dynamically Melting a Mott Insulator,” Phys. Rev. Lett. 89, 040402 (2002) [CrossRef] [PubMed]
6. Min-Chul Cha, M. P. A. Fischer, S. M. Girvin, M. Wallin, and A. P. Young, “Universal conductivity of two-dimensional films at the superconductor-insulator transition,” Phys. Rev. B 44, 6883–6902 (1991). [CrossRef]
7. M. P. A. Fischer, G. Grinstein, and S. M. Girvin “Presence of quantum diffusion in two dimensions: Universal resistance at the superconductor-insulator transition,” Phys. Rev. Lett.64, 587–590 (1990);M. P. A. Fischer, “Quantum phase transitions in disordered two-dimensional superconductors,” Phys. Rev. Lett.65, 923–926 (1990). [CrossRef]
8. S. Sachdev, Quantum Phase Transitions, (Cambridge Univ. Press, Cambridge, 2001).
9. M. P. A. Fischer, P. B. Weichman, G. Grinstein, and D. S. Fisher, “Boson localization and the superfluid-insulator transition,” Phys. Rev. B 40, 546–570 (1989). [CrossRef]
10. J. K. Freericks and H. Monien, “Phase diagram of the Bose Hubbard model,” Europhys. Lett. 26, 545–550 (1994). [CrossRef]
11. J. K. Freericks and H. Monien, “Strong-coupling expansions for the pure and disordered Bose-Hubbard model,” Phys. Rev. B 53, 2691–2700 (1996). [CrossRef]
12. D. van Oosten, P. van der Straten, and H. T. C. Stoof, “Quantum phases in an optical lattice,” Phys. Rev. A 63, 053601 (2001). [CrossRef]
13. A. M. Rey, K. Burnett, R. Roth, M. Edwards, C. J. Williams, and C. W. Clarck, “Bogoliubov Approach to Super-fluidity of Atoms in an Optical Lattices,” J. Phys. B , 36, 825–841 (2003). [CrossRef]
14. R. Fazio and H. van der Zant, ‘Quantum phase transitions and vortex dynamics in superconducting networks,” Phys. Rep. 355, 235 (2001) and ref. therein. [CrossRef]
15. A. van Otterlo, K.-H. Wagenblast, R. Baltin, C. Bruder, R. Fazio, and G. Schön, “Quantum phase transitions of interacting bosons and the supersolid phase,” Phys. Rev. B 52, 16176–16186 (1995). [CrossRef]
18. S. R. White, “Density-matrix algorithms for quantum renormalization groups,” Phys. Rev. B. 48, 10345–10356 (1993). [CrossRef]
19. I. Peschel, X. Wang, M. Kaulke, and K. Hallberg, Density-matrix Renormalization, (Springer-Verlag, Berlin, 1998).
20. T. D. Kühner, S. R. White, and H. Monien, “One-dimensional Bose-Hubbard model with nearest-neighbor interaction,” Phys. Rev. B 61, 12474–12489 (2000). [CrossRef]
21. T. D. Kühner, Diploma work (1997), University of Bonn.
22. S. Rapsch, U. Schollwöck, and W. Zwerger, “Density matrix renormalization group for disordered bosons in one dimension,” Europhys. Lett. 46, 559–564 (1999). [CrossRef]
24. G. G. Batrouni and R. T. Scalettar, “World-line quantum Monte Carlo algorithm for a one-dimensional Bose model,” Phys. Rev. B 46, 9051 (1992). [CrossRef]
25. P. Niyaz, R. T. Scalettar, C. Y. Fong, and G. G. Batrouni, “Phase transitions in an interacting boson model with near-neighbor repulsion,” Phys. Rev. B 50, 362–373 (1994). [CrossRef]
26. N. V. Prokof’ev, B. V. Svistunov, and I. S. Tupitsyn, Phys. Lett. A 238, 253 (1998). [CrossRef]
27. N. Elstner and H. Monien, “Dynamics and thermodynamics of the Bose-Hubbard model,” Phys. Rev. B 59, 12184–12187 (1999) and ref. therein. [CrossRef]
28. D. S. Rokhsar and B. G. Kotliar, “Gutzwiller projection for bosons,” Phys. Rev. B 44, 10328–10332 (1991). [CrossRef]
29. J. R. Anglin, P. Drummond, and A. Smerzi, “Exact quantum phase model for mesoscopic Josephson junctions,” Phys. Rev. A 64, 063605 (2001). [CrossRef]