Confinement-induced resonance from the generalized Gross-Pitaevskii equations
Alexander Yu. Cherny





cherny@theor.jinr.ru
Bogoliubov Laboratory of Theoretical Physics, Joint Institute for Nuclear Research, 141980, Dubna, Moscow region, Russia





The confinement-induced resonances for trapped bosons in the cigar-shaped and pancake geometries are studied within the generalized Gross-Pitaevskii equations, which are a simplified version of the Hartree-Fock-Bogoliubov approximation. Although the Hartree-Fock-Bogoliubov method is considered applicable only for small interparticle  interactions, the resonance denominators for the chemical potential are obtained in both quasi-one and quasi-two dimensions. A useful integral representation of the one-particle Green's function are found for the cylindrical confinement. We find the position of a smoothed resonance for the chemical potential in the pancake geometry at positive scattering length.


Confinement-induced resonance from the generalized Gross-Pitaevskii equations
    Alexander Yu. Cherny
    April 29, 2024
=============================================================================







§ INTRODUCTION



The confinement-induced resonance (CIR) is an efficient tool to tune effective interparticle interactions in low-dimensional systems of cold gases <cit.>. The resonance was predicted in the cigar-shaped  <cit.> and pancake <cit.> geometries
and confirmed experimentally <cit.>. The cases of very tight and anisotropic confinements were studied in Refs. <cit.> and <cit.>, respectively.

The tight cylindrical or flat atomic waveguide creates, in effect, a low-dimensional system of cold atoms, and the problem is how to describe the corresponding interactions in low dimensions. The standard approach <cit.> is to obtain interparticle interactions in low dimensions from the two-body scattering problem, thus reducing the effective dimension of the many-body system.

For identical bosons, it seems more attractive to develop a scheme that directly describes the wavefunction of the Bose-Einstein condensate  (BEC) for trapping potentials of arbitrary shape. This is all the more interesting because various optical traps are currently being produced experimentally, including an optical-box trap <cit.>. The well-known Gross-Pitaevskii (GP) equation <cit.> is not suitabe for this purpose, since it does not reproduce the resonance denominators for effective interactions in low dimensions. For instance, in the cigar-shaped geometry, the wavefunction of BEC coincides up to a factor with the ground state of the trapping harmonic oscillator in the x-y plane provided its frequency is sufficiently high. Then one can integrate out the x and y coordinates in the interaction term of the GP functional and obtain the effective coupling constant in one dimension <cit.>: g_1D= 2ħ^2a/(m l^2). This result is valid only away from the resonance when a≪ l. Here a, m, and l are the 3D scattering length in free space, boson mass, and length of the two-dimensional oscillator, respectively. On the other hand, the correct coupling constant contains the resonance denominator <cit.>, see Eq. (<ref>) below.

In the paper <cit.>, a system of non-linear equations in both stationary and time-dependent cases was proposed and called the generalized Gross-Pitaevskii (GGP) equations. Apart from the standard one-body wavefunction of BEC, the equations also contain the two-body wavefunction of BEC.  This approach is nothing else but a simplified version of the full Hartree-Fock-Bogoliubov (HFB) approximation (see Appendix <ref> below). The GGP equations are applicable for arbitrary shape of the trapping potential. In particular, the resonance denominators were obtained <cit.> in both cases of cigar-shaped and pancake geometries. However, the resulting denominators of Ref. <cit.>, which determine the position of the resonances, contained incorrect coefficients, since the solutions were derived with an additional quite rough approximation. As is shown in this paper, the GGP equations indeed yield the exact coefficients, which are consistent with the two-body scattering problem <cit.>.

The GGP equations are nonlinear, but the effective coupling constants in low dimensions are determined by the short-range correlations, where the nonlinearity is not important. The problem can be reduced to the two-body problem and the corresponding Green's functions (GFs). We find and examine the integral representations for the one-particle GFs at low energies for the flat and cylindrical harmonic confinements. In the latter case, the equation for the GF differs from the previously found integral representation <cit.>, which is actually the Fourier transform of the Feynman propagator for a harmonic oscillator with contour rotation in the complex plane. We believe that the found representation is more convenient for obtaining the low-energy expansions of the GF.

This paper is organized as following. In the next section, the GGP equations are written and discussed in the case of the tight confinements. In Sec. <ref>, the integral representations and expansions of the GFs are obtained for the cylindrical and flat confinements. The CIRs for the effective coupling constant and chemical potential are found within the suggested scheme in Secs. <ref> and <ref>. The main results and prospects are discussed in Conclusions. In Appendix <ref>, the GGP are derived through the HFB approximation.




§ GENERALIZED GROSS-PITAEVSKII EQUATIONS



We consider a dilute many-body system of identical bosons, interacting with a pairwise short-range[A potential is of the short-range type if it decreases at infinity as 1/r^α with the exponent α > 3 or faster (see, e.g., Ref. <cit.>). For such a potential, only the s-wave scattering amplitude survives in the limit of zero scattering energy. It is called the scattering length [see Eq. (<ref>) below].]  potential V(r), in an external field V_ext(r). For simplicity, the bosons are assumed to be spinless. The equilibrium properties of such system at zero temperature are described by Eqs. (16) and (17) of Ref. <cit.>, which can be written in the form

    [Ĥ(r_1) -μ]ϕ(r_1)+∫ dr_2ϕ(r_2)V(|r_1-r_2|)φ(r_1,r_2)=0,   
       [Ĥ(r_1) +Ĥ(r_2)-2μ]ψ(r_1,r_2)+V(|r_1-r_2|)φ(r_1,r_2)=0,

where μ is the chemical potential. Here we set by definition

    ψ(r_1,r_2)= φ(r_1,r_2) - ϕ(r_1)ϕ(r_2),

and the one-body Hamiltonian is given by

    Ĥ(r)=-ℏ^2∇_r^2/2m + V_ext(r)-E_0

with E_0 being its ground-state energy. By definition, ϕ(r) ≡⟨Ψ̂(r)⟩ and φ(r_1,r_2) ≡⟨Ψ̂(r_1)Ψ̂(r_2)⟩ are the anomalous averages of the field Bose operators. The derivation of these equations and their relation to the HFB approximation are given in Appendix <ref>.

The function ϕ(r) is the one-body function of inhomogeneous BEC and can be considered as the order parameter. The function φ(r_1,r_2) is called pair wavefunction. It describes two bosons in the Bose-Einstein condensate, and this function is not reduced in general to the product of the one-body functions. The exact definitions can be done through the one-body and two-body matrices <cit.>. The property of decay of correlations dictates the boundary condition

    φ(r_1,r_2) ≃ϕ(r_1)ϕ(r_2)  when |r_1-r_2|≫ξ_h,

where ξ_h is the healing (or coherence) length defined through the chemical potential: ξ_h=ħ/√(μ m). The system of equations (<ref>) and (<ref>) is nonlinear due to the presence of the condensate wavefunction ϕ(r_2) under the integral in Eq. (<ref>). The GGP equations (<ref>) and (<ref>), being a particular case of HFB approximation, are applicable to weak interactions when the mean distance between bosons is much smaller than the healing length.

Various solutions for different external potentials were considered in Ref. <cit.>. In particular, the 1D and 2D harmonic CIRs were obtained [see Eqs. (48) and (57) of Ref. <cit.>]. However, the used method did not enable us to obtain the precise coefficients in the resonance denominators, which determine the CIR positions. In this paper, we show that the correct coefficients can be derived from Eqs. (<ref>) and (<ref>).

We assume that the radius r_0 of the short-range interaction potential V(r) is much smaller than the characteristic length of the external potential V_ext(r) and the mean distance between bosons. This implies that at low energies, physical processes are governed by the s-wave scattering length a of the interaction potential, which is obtained from the symmetric solution of the two-body Scrödinger equation at zero energy in free space.

Thus, our task is to get the chemical potential of N strongly confined bosons at zero temperature as a function of the 3D scattering length a and oscillator frequency ω in the case of cylindrical and flat harmonic confinements.




§ GREEN'S FUNCTIONS OF THE ONE-BODY SCHRÖDINGER EQUATION AT LOW ENERGIES FOR CYLINDRICAL AND FLAT HARMONIC CONFINEMENTS



To obtain the chemical potentials of the Bose-Einstein condensates for the cigar-shaped and pancake geometries, we need to calculate the corresponding GFs for the one-body Hamiltonian (<ref>), which is defined as

    Ĝ_E=-(Ĥ-E)^-1

with the scattering energy E, which is assumed to be sufficiently small (see the discussions below).

In this section, we find useful integral representations for the GFs in the cases of cylindrical and flat harmonic confinements. These representations are very convenient for obtaining the low-energy expansions of the GFs.



 §.§ Cylindrical confinement



The cylinder harmonic confinement implies that a particle is confined in the x-y plane but can move freely in z-direction with a finite momentum, whose absolute value is denoted by ħ p. The one-particle Hamiltonian (<ref>) is given by

    Ĥ=-ℏ^2∇^2/2m_0 + m_0ω^2(x^2+y^2)/2-E_0,

where the shift E_0=ħω is introduced to get zero ground-state energy.

It is convenient to use the following notation for radius vector: r=(ρ,z) with ρ=(x,y). We calculate the GF (<ref>) in coordinate representation with the radial coordinate ρ=√(x^2+y^2) being equal to zero and E=ħ^2p^2/2m_0:

    G_p(z-z')≡⟨ z,ρ=0|Ĝ_E|ρ'=0,z'⟩
       =-1/2π^2l^2∑_n=0^∞p.v.∫_-∞^∞dqe^iq(z-z')/2ħω n+ħ^2/2m_0(q^2-p^2).

Here l_0=√(ħ/(m_0ω)) is the harmonic oscillator length, and p.v. denotes the Cauchy principal value. This kind of regularization ensures the reality of the GF in coordinate representation <cit.>.

Here we use the expression for the GF in Dirac notation ⟨r |Ĝ_E|r' ⟩ =-∑_k⟨r |k⟩⟨ k|r' ⟩ /(E_k-E) with |k⟩ and E_k being the eigenfunctions  and eigenvalues of a Hamiltonian, respectively. For the Hamiltonian (<ref>), k=(n,m,q) is a multi-index, where n=0,1,2,… is the main quantum number, m=0,±1,±2,… is the orbital momentum in 2D, and q is the z-component of the wavevector. The eigenvalues and normalized eigenfunctions are given by E_k=ħω(2n+|m|)+q^2/2m_0 and φ_k(ρ,α,z)=e^iqz/√(2π)ϕ_n,m(ρ,α), respectively. Here ϕ_n,m(ρ,α) are the eigenvalues for the 2D oscillator in the polar coordinates x=ρcosα, y=ρsinα. One can show that for an arbitrary main quantum number, ϕ_n,0(ρ=0,α)=1/(√(π)l_0) and ϕ_n,m(ρ=0,α)=0 for m≠0. Since we are interested at the GF of the Hamiltonian (<ref>) only at ρ=ρ'=0, we are left with one sum over the main quantum number and thus arrive at Eq. (<ref>).

When the scattering momentum is sufficiently small, p<2/l_0, only the first term in the sum in Eq. (<ref>) is singular. Separating this term and using relations ∫_-∞^∞dq/2πe^iq z/(q^2+ϰ^2) =e^-ϰ|z|/2ϰ and p.v.∫_-∞^∞dq/2πe^iq z/(p^2-q^2) =sin p|z|/2p yield

    ħ^2/2m_0G_p(z)=   1/π l_0^2[sin p|z|/2p-l_0/4Λ(2|z|/l_0,-p^2l_0^2/4)], 
    Λ(ξ,ϵ)=   ∑_n=1^∞e^-ξ√(n+ϵ)/√(n+ϵ),

where dimensionless variables ξ =2|z|/l_0 and ϵ =-p^2l_0^2/4 are introduced.

The first term of the GF (<ref>) is proportional to the one-dimensional GF of a free particle, which obeys the relation (∂^2/∂ z^2+p^2)sin p|z|/2p =δ(z).  At large distances |z|≫ l_0, the second term of the GF falls off exponentially as it follows directly from the expansion (<ref>). Then the asymptotics of the scattering part of the wavefunction, proportional to the GF, is essentially one-dimensional, as expected.

Let us calculate the GF at low momenta |ϵ|≪ 1 and small coordinates ξ≪ 1. To expand Λ(ξ,ϵ) in ξ and ϵ around zero, we use the identity for exponential e^-x=1/√(π)∫_0^∞dt/t^3/2e^-1/te^-x^2t/4 and find the integral representation

    Λ(ξ,ϵ)=1/√(π)∫_0^∞dt/t^3/2e^-1/te^-ξ^2(1+ϵ)t/4Φ(e^-ξ^2 t/4,1/2,1+ϵ),

where Φ(x,s,α)=∑_n=0^∞x^n/(n+α)^s is the Lerch transcendent, or Lerch zeta function <cit.>. Upon expanding the integrand

    e^-ξ^2(1+ϵ)t/4Φ(e^-ξ^2t/4,1/2,1+ϵ)=   2 √(π)/√(t)ξ +ζ(1/2) -ϵ/2ζ(3/2) 
        + O(ξ^2) + O(ϵ^2)

and substituting it into Eq. (<ref>), we obtain

    Λ(ξ,ϵ)=   2/ξ+ζ(1/2)+ξ/2-ϵ/2ζ(3/2)+O(ξ^2)+ O(ϵ^2)
       +O(ξϵ)

with ζ(x) being the Riemann zeta function. The symbol O(x) is the big O notation, which designates terms of the order of x or smaller. When |z|≪ l_0 and p≪ 1/l_0, the expansion of the GF (<ref>) takes the form

    ħ^2/2m_0G_p(z)=   -1/4π|z|-1/4π l_0ζ(1/2)+|z|/4π l_0^2-p^2l_0/32πζ(3/2)
       +O(z^2/l_0^2)+O(|z|  p^2 l_0)+O(p^4 l_0^4).





 §.§ Flat confinement




The flat harmonic confinement along the z-axis is considered by close analogy with the previous subsection. A particle, being confined in z-direction, can move freely in the x-y plane. The Hamiltonian reads

    Ĥ=-ℏ^2∇^2/2m_0 + m_0ω^2z^2/2-E_0

with the shift E_0=ħω/2. A complete set of its eigenfunctions is given by φ_q,n(ρ,z) =ϕ_n(z)e^iq·ρ/(2π) with ϕ_n(z) being the normalized eigenfunctions of the linear harmonic oscillator. We need only their absolute values at z=0: when n is odd then ϕ_n(0)=0, while for even n

    |ϕ_2k(0)|=1/π^1/4 l_0^1/2√((2k)!)/2^kk!,   k=0,1,2,…,

see, e.g., the textbook <cit.>.

Using the same method as in the previous subsection, we find the expression for the GF at z=z'=0 and E=ħ^2p^2/2m_0

    G_p(|ρ-ρ'|)≡⟨ρ,z=0|Ĝ_E|z'=0,ρ'⟩
       =-∑_k=0^∞|ϕ_2k(0)|^2p.v.∫d^2q/(2π)^2e^iq·(ρ-ρ')/2ħω k+ħ^2/2m_0(q^2-p^2),

where |ϕ_2k(0)| is given by Eq. (<ref>). When the momentum is small, p<2/l_0, only the first term in the sum is singular. Separating it and using the identities p.v.∫d^2q/(2π)^2e^iq·ρ/p^2-q^2 =Y_0(pρ)/4, ∫d^2q/(2π)^2e^iq·ρ/ϰ^2+q^2 =K_0(ϰρ)/2π for positive p and ϰ yield

    ħ^2/2m_0G_p(ρ)=   1/√(π)l_0[Y_0(pρ)/4-1/2πS(2ρ/l_0,-p^2l_0^2/4)], 
    
      S(ξ,ϵ)=   ∑_k=1^∞K_0(ξ√(k+ϵ))(2k)!/4^k(k!)^2

with the dimensionless variables ξ=2ρ/l_0 and ϵ=-p^2l_0^2/4. Here Y_0(x) and K_0(x) are the Bessel and modified Bessel functions of the second kind, respectively (see, e.g. <cit.>, Ch. 9).

Similar to the one-dimensional case considered in Sec. <ref>, the first term is proportional to the GF of a free particle in two dimensions: (∂^2/∂ρ^2+1/ρ∂/∂ρ+p^2)Y_0(pρ)/4 =δ(ρ). When ρ≫ l_0, the second term in Eq. (<ref>) is exponentially small, which actually implies two-dimensional scattering.

The relation K_0(x)=∫_0^∞dt/2te^-te^-x^2/(4t) enables us to sum up the series (<ref>) under the integral and find its integral representation in elementary functions

    S(ξ,ϵ)=∫_0^∞dt/2texp(-t-ϵ ξ^2/4t)[1/√(1-e^-ξ^2/(4t))-1].

This relation is a particular case of a more general equation, which was obtained by another method in Refs. <cit.>. This general equation is actually the Fourier transform of the well-known Feynman propagator for a harmonic oscillator.

It is easy to obtain the expansion

    S(ξ,ϵ)=   √(π)/ξ+lnξ +γ-√(π)-ln 2+C-ξ√(π)/8-ϵln 2
       +O(ξ^2)+ O(ϵ^2)+O(ξϵ),

where γ=0.5772… is Euler's constant and

    C=∫_0^∞dy  e^-y/2y[1/√(1-e^-y))(1+√(1-e^-y)))+1-1/√(y)].

Its numerical value is given by C=0.8035589….

To find the GF at ρ≪ l_0 and p≪ 1/l_0, we substitute the expansions (<ref>) and Y_0(x)=2ln(xe^γ/2)/π + O(x^2ln x) into Eq. (<ref>) and finally arrive at

    ħ^2/2m_0G_p(ρ)=   -1/4πρ+1/4π^3/2l_0ln(π p^2l_0^2/2B)
       +1/8π l_0[ ρ/l_0-ln2/√(π)p^2l_0^2]+O(ρ^2 p^2ln pρ)
       +O(ρ^2/l_0^2)+O(ρ p^2l_0)+O(p^4l_0^4),

where we put by definition

    B≡2πexp[2(C-√(π))]=0.904916….

The first two terms of GF (<ref>) were first calculated in <cit.> with the constant B being equal to one. In the subsequent paper <cit.>, its numerical value was estimated to be 0.915. Further, the value of B being equal to 0.9049 was obtained in Ref. <cit.>, which is consistent with Eq. (<ref>).
We believe that this is the most accurate estimation.







§ CONFINEMENT-INDUCED RESONANCE IN ONE DIMENSION: CYLINDRICAL GEOMETRY





 §.§ The tight confinement



When the bosons are trapped in the x-y transverse direction with the axially symmetric harmonic potential V_ext(ρ) =mω^2 ρ^2/2, the order parameter depends only on the transverse coordinates, but it is homogeneous along the z-axis.

The tight confinement implies that the chemical potential is small in scale of ħω, that is, l≪ξ_h, in terms of the oscillator length l =√(ħ/(mω)). Then with a good accuracy one can put for the wavefuncton of BEC <cit.> ϕ =ϕ(ρ_1) =√(n_1D)ϕ_0(ρ_1) with ϕ_0(ρ_1) =exp[-ρ^2_1/(2l^2)]/l√(π) being the ground state wavefunction of the one-particle Schrödinger equation with Hamiltonian (<ref>), and n_1D is the linear density of the Bose-Einstein condensate. Thus, the nonlinear term in Eq. (<ref>) is treated as a perturbation, and the chemical potential is obtained by multiplying Eq. (<ref>) by ϕ_0(ρ_1) and integrating out over r_1:

    μ=1/L∬ dr_1dr_2ϕ_0(ρ_1)ϕ_0(ρ_2)V(|r_1-r_2|)φ(r_1,r_2),

where L is the length of the system along z-axis.

To take the integral (<ref>), we need to examine the properties of the pair wavefunction. In the case of transverse confinement, Eqs. (<ref>) and (<ref>) are translationally invariant under a spatial translation z_1→ z_1+b, z_2→ z_2+b for any b, which yields φ =φ(ρ_1,ρ_2,|z_1-z_2|). Moreover, for the harmonic confinement, the spatial variables in Eq. (<ref>) are separated into the center-of-mass position and relative coordinates by analogy with the two-body problem. We can look for a solution in the form

    φ(r_1,r_2)=φ(ρ,z)exp[-(ρ_1+ρ_2)^2/(4l^2)]

and the same relation for ψ(ρ_1,ρ_2,z), see Eq. (<ref>). Here ρ=ρ_1-ρ_2 and z=z_1-z_2 are the relative coordinates, and ρ =|ρ|. The boundary conditions (<ref>) read φ(ρ,z) ≃ n_ 1Dexp[-ρ^2/(4l^2)]/(π l^2) when z≫ξ_h.

In the range r≪ l, one can use the approximation

    φ(ρ,z)=ηφ^(0)_ 3D(r)n_ 1D/π l^2,

since r_0≪ l. Here r=√(ρ^2+z^2), and φ^(0)_ 3D(r) is the 3D solution of the two-body Schrödinger equation in the center-of-mass system at zero energy

    -(ħ^2/m)∇^2φ^(0)_ 3D(r)+V(r)φ^(0)_ 3D(r)=0.

The solution is assumed to be regular in the vicinity of the origin and obeys the boundary condition when r≫ r_0

    φ^(0)_ 3D(r)≃     1-a/r, 
    
    a=    m/4πħ^2∫ drV(r)φ^(0)_ 3D(r).

Substituting Eqs. (<ref>) and (<ref>) into Eq. (<ref>) and using the difference of the  scales r_0≪ l and Eq. (<ref>), we are left with

    μ=η2ħ^2a/m l^2n_ 1D.


Here the point is that the constant η≠1, which is crucial for the appearance of CIR. This constant can be recovered from Eq. (<ref>).



 §.§ The application of the Green's function



Substituting Eq. (<ref>) into Eq. (<ref>), we arrive at the equation for the wavefunction φ(ρ,z):

    [-ħ^2∇^2/2m_0+m_0ω^2ρ^2/2-E_0-2μ]ψ(ρ,z)=-V(r)φ(ρ,z)

with m_0=m/2 being the reduced mass. It is sufficient to find the wavefunction φ(ρ,z) at ρ=0. According to Eqs. (<ref>) and (<ref>), it has the asymptotics

    φ(ρ=0,z)≃ηn_ 1D/π l^2(1-a/|z|)

in the short-range regime r_0≪ |z|≪ l.

The solution of Eq. (<ref>) can be found with the Green's function (GF) at E=2μ=ħ^2p^2/2m_0 as ψ(ρ,z) =∫ dr' ⟨ z,ρ|Ĝ_E|ρ',z'⟩ V(r')φ(ρ',z'). In the short-range regime |z|≪ l and at low energies p≪ 1/l, the GF at ρ=ρ'=0 was obtained in Sec. <ref> and given by Eq. (<ref>). At m_0=m/2 and l_0=√(ħ/(m_0 ω))=√(2)l, we get

    φ(ρ=0,z)=n_ 1D/π l^2+G_p(z)∫ dr'V(r')φ(ρ',z'),

where Eq. (<ref>) is used: φ(ρ=0,z) =n_ 1D/(π l^2) +ψ(ρ=0,z).

The potential V(r) is localized within r≲ r_0, where the wavefunction φ(ρ,z) takes the form (<ref>). Using the relation (<ref>) and the first two terms of Eq. (<ref>) for the Green's function, we find

    φ(ρ=0,z)≃n_ 1D/π l^2[1-η  a(1/|z|+ζ(1/2)/√(2)l)].

Comparing this asymptotics with Eq. (<ref>) finally yields

    η=(1+a ζ(1/2)/√(2)l)^-1,

where ζ(1/2)=-1.46035….



 §.§ The chemical potential



The chemical potential is given by Eqs. (<ref>) and (<ref>)

    μ   =g_1Dn_ 1D,
    
    g_1D   = 2ħ^2a/m l^21/1+a ζ(1/2)/√(2)l.

Taking into account the relation g_1D=-2ħ^2/m  a_1D for the coupling constant in 1D, we arrive at the expression for the effective 1D scattering length a_1D in the two-body problem (see, e.g., Ref. <cit.>), first obtained in Ref. <cit.>

    a_1D=-l^2/a(1+a ζ(1/2)/√(2)l).







§ CONFINEMENT-INDUCED RESONANCE IN TWO DIMENSIONS: FLAT GEOMETRY





 §.§ The tight confinement



The flat confinement of identical bosons is provided by harmonic forces along the z-axis: V_ext(z) =mω^2 z^2/2. The order parameter is homogeneous within the x-y plane, and in the case of the tight confinement l≪ξ_h, one can put with a good accuracy ϕ =ϕ(z_1) =√(n_2D)ϕ_0(z_1) with the oscillator ground-state wavefunction ϕ_0(z_1) =exp[-z_1^2/(2l^2)] /(l^1/2π^1/4). Here n_2D is the 2D density of the Bose-Einstein condensate. The chemical potential can be obtained from Eq. (<ref>) by analogy with the previous section

    μ=1/S∬ dr_1dr_2ϕ_0(z_1)ϕ_0(z_2)V(|r_1-r_2|)φ(r_1,r_2),

where S is the area of the system occupied in the x-y-plane.

The pair wavefunction is invariant under translations along the x-y plane: φ(r_1,r_2)=φ(z_1,z_2,ρ). It can be written for the harmonic confinement as φ(r_1,r_2)=φ(z,ρ)exp[-(z_1+z_2)^2/(4l^2)].

The approximation for the pair wavefunction for r≪ l

    φ(z,ρ)=ηφ^(0)_ 3D(r)n_ 2D/l√(π),

is analogous to Eq. (<ref>). In the same manner, we find from Eqs. (<ref>), (<ref>) and (<ref>)

    μ=4πħ^2/mn_ 2D u,

where we introduce a more convenient dimensionless parameter

    u=η a/√(2π)l.


The prefactor η can also be found by analogy with Sec. <ref>.



 §.§ The application of the Green's function



On the one hand, in the shot-range regime r_0≪ρ≪ l, Eqs. (<ref>) and (<ref>) yield

    φ(z=0,ρ)≃ηn_ 2D/l√(π)(1-a/ρ).

One the other hand, one can use the GF for Eq. (<ref>) (see Sec. <ref>) and obtain in the same regime

    φ(z=0,ρ)≃    n_ 2D/√(π) l+G_p(ρ)∫ dr'V(r')φ(z',ρ')
    
                            =    n_ 2D/√(π) l+G_p(ρ)ηn_ 2D/l√(π)4πħ^2a/m

with the GF (<ref>), where 2μ=ħ^2p^2/2m_0, m_0=m/2, l_0=√(2)l=√(2ħ/(m ω)).



 §.§ The chemical potential




In the quasi-1D condensate, considered in the previous section <ref>, the GF remains finite in the limit p→ 0. By contrast, in the quasi-2D condensate, the GF (<ref>) exhibits a logarithmic divergence as p→ 0. This implies that the chemical potential should be determined self-consistently. Indeed, comparing Eq. (<ref>) with Eq. (<ref>) and using the definition (<ref>) yield the algebraic equation for u

    1/u+ln u=δ, 
       δ=√(2π)l/a-ln[8π^2 n_ 2Dl^2/B],

where the constant B is defined by Eq. (<ref>). Thus, the chemical potential is given by Eqs. (<ref>), (<ref>), and (<ref>).

Equation (<ref>) has no real solutions when δ<1 and two solutions when δ>1: one is greater than 1, and the other is less than 1.  It follows from the thermodynamic stability condition ∂μ/∂ n_ 2D =4πħ^2/mu/1-u>0 that only the solution u<1 is physical. This fits well with the condition of weak coupling, which tells us that the healing length should be at least smaller than the mean distance between bosons on the plane: ħ/√(μ m)>1/√(n_ 2D), see Eq. (<ref>). When the dimensionless parameter δ is big then u≃ 1/δ and we arrive at the weak-coupling relation for the chemical potential

    μ≃4πħ^2 n_ 2D/mδ,

which allows us to interpret δ, given by Eq. (<ref>), as the resonance denominator.  Note, however, that this approximation overestimates the values of the chemical potential, see Fig. <ref>a.



The CIR arises (see Fig. <ref>b) when the 3D scattering length is negative while the 2D density is sufficiently small for δ to be positive and close to 1, which is its minimum possible value within our approach. At positive scattering length, a reminiscent of the resonance is also observed, when the dimensionless chemical potential reaches the local maximum at

    l/a =√(2/π)=0.7978…

as it follows from Eq. (<ref>). This differs from experimentally found <cit.> value l/a=1.19(3), at which the CIR is realized at positive scattering lengthes.

Relations (<ref>), (<ref>), and (<ref>) should be used with some caveats. First, the quasi 2D regime imposes the restriction on the oscillator length: the characteristic wavevector of bosons p≃√(n_ 2D) should be smaller than 2/l as discussed in Sec. <ref>. This yields l/a≲2/√(n_ 2Da^2), which for the chosen parameter n_ 2Da^2=0.05 in Fig. <ref>b is of order of 9. Second, we assume that the range of interaction potential is much smaller than the oscillator length: r_0≪ l. Thus, in the vicinity of l=0, the approximation of zero-range interactions (r_0=0) is used. We believe that the extremely tight confinement r_0≳ l also can be described by means of the GGP equations with methods similar to those considered in Ref. <cit.>.





§ CONCLSUSIONS



Using the GGP equations (<ref>) and (<ref>), we obtained the correct resonance denominators, which determine the positions of CIRs in the cigar-shaped and pancake geometries, see Eqs. (<ref>), (<ref>), and (<ref>), (<ref>).

The GGP equations, being a particular case of the HFB approximation (see Appendix <ref>), are applicable for a weakly interacting Bose gas, strictly speaking. Nevertheless, the GGP equations correctly reproduce the effective coupling in the low-dimensional geometries. The reason is that these equations describe well the short-range correlations even inside the radius of interaction potential. This confirms that the GGP equations can be considered as a universal tool for describing BECs in various regimes. For instance, the method has good prospects for studying the extremely tight confinements, crossover between 3D and 2D or 1D regimes, and the dipolar BEC in various regimes, including quantum droplets <cit.>.

In addition, the integral representation for GF was found for the cylindrical confinement [Eqs. (<ref>), (<ref>)], which is valid for small scattering wavevectors p<2/l_0 and arbitrary values of coordinates. It is convenient for obtaining the low-density expansion (<ref>) of the GF.

Within our approach, we found the position (<ref>) of the smoothed resonance at positive scattering length in the pancake geometry. The position is independent of the density. This is nearly two-third of the experimentally observed value <cit.>. The discrepancy might be attributed to condensate depletion at finite temperatures or to many-body effects that are beyond the used weak-coupling approximation.



The author is grateful to Azis Embareck for useful and fruitful discussions and Vladimir Melezhik for careful reading of the manuscript and useful comments.





§ HARTREE-FOCK-BOGOLIUBOV APPROXIMATION



Let make some remarks about the Hartree-Fock-Bogoliubov (HFB) method. We consider only the stationary HFB theory at zero temperature. The HFB approximation for fermions was suggested by Bogoliubov <cit.> as a variational scheme formulated in terms of the normal and anomalous averages of two Fermi-field operators. These scheme can be simplified <cit.> if we note that these averages are given by bilinear forms of u and v parameters of the Bogoliubov transformation. Then the variational procedure can be reformulated directly in terms of u and v parameters, and the resulting equations are called Bogoliubov-de Gennes equations.

For bosons, the relations analogous to the Bogoliubov-de Gennes equations can be written down <cit.>. They are called the Bogoliubov equations. In addition to the u and v parameters, these equations contain the wavefunction of the inhomogeneous Bose-Einstein condensate.

On the other hand, one can also formulate the HFB method for bosons as a variational scheme with three variational functions: the condensate wavefunction and normal and anomalous averages of two Bose-field operators. Below we follow this scheme <cit.> and derive the generalized GP equations by further simplifying the HFB scheme.




 §.§ Hartree-Fock-Bogoliubov approximation for bosons in the variational formulation


The ground-state energy and mean number of interacting bosons in an external field is given by (see, e.g., textbooks <cit.>)

    E =   -ℏ^2/2m∫ d x.∇^2_x⟨Ψ̂^†(x')Ψ̂(x)⟩|_x'=x
       +∫ d x  V_ ext(x)⟨Ψ̂^†(x)Ψ̂(x)⟩
       +1/2∫ d x_1 d x_2 V(x_1,x_2) ⟨Ψ̂^†(x_1)Ψ̂^†(x_2)
    Ψ̂(x_2)Ψ̂(x_1)⟩,
    
    
    N=   ∫ d x ⟨Ψ̂^†(x)Ψ̂(x)⟩,

respectively. Here x=( r,σ) are the coordinate and spin of a particle, respectively, and ∫ d x⋯ = ∑_σ∫ d r⋯. The brackets ⟨⋯⟩ stand for the ground-state average. The Bose field operators are denoted as Ψ̂(x) and Ψ̂^†(x).

The Bose-Einstein condensation implies that the global gauge symmetry is broken. Then one can separate the condensate contribution to the field operators: Ψ̂(x)=ϕ(x) +ϑ̂(x) and Ψ̂^†(x)=ϕ^*(x) +ϑ̂^†(x). Here ϕ(x) and ϑ̂(x) are the numerical and operator parts, respectively, for which we have ϕ(x)=⟨Ψ̂(x)⟩ and ⟨ϑ̂(x)⟩=0. Then the two-boson correlators read

    ⟨Ψ̂^†(x)Ψ̂(x')⟩ =    ϕ^*(x)ϕ(x') +⟨ϑ̂^†(x) ϑ̂(x')⟩,
    ⟨Ψ̂(x)Ψ̂(x')⟩ =    ϕ(x)ϕ(x') +⟨ϑ̂(x) ϑ̂(x')⟩.

The normal two-boson correlator (<ref>) is nothing else but the one-body density matrix, and the anomalous correlator (<ref>) can be treated as the two-body wavefunction of the Bose-Einsten condensate <cit.>.

The diagonal part of the two-body density matrix, which determines the interaction term, can also be found by separating the condensate part and using Wick's theorem for the correlators of ϑ̂ and ϑ̂^† operators

    ⟨Ψ̂^†(x_1)Ψ̂^†(x_2)Ψ̂(x_2)Ψ̂(x_1)⟩
    =ϕ^*(x_1)ϕ(x_2)^*ϕ(x_2)ϕ(x_1)
       +ϕ(x_1)ϕ(x_2)⟨ϑ̂^†(x_1) ϑ̂^†(x_2)⟩
    +ϕ^*(x_1)ϕ^*(x_2)⟨ϑ̂(x_2) ϑ̂(x_1)⟩
       +ϕ(x_1)ϕ^*(x_2)⟨ϑ̂^†(x_1) ϑ̂(x_2)⟩
    +ϕ^*(x_1)ϕ(x_2)⟨ϑ̂^†(x_2) ϑ̂(x_1)⟩
       +⟨ϑ̂^†(x_1) ϑ̂^†(x_2)⟩⟨ϑ̂(x_2) ϑ̂(x_1)⟩
       +⟨ϑ̂^†(x_1) ϑ̂(x_1)⟩⟨ϑ̂^†(x_2) ϑ̂(x_2)⟩
       +⟨ϑ̂^†(x_1) ϑ̂(x_2)⟩⟨ϑ̂^†(x_2) ϑ̂(x_1)⟩.

Wick's theorem is applicable, since the effective HFB Hamiltonian contains only a bilinear form of creation and annihilation operators <cit.>.

The normal ⟨ϑ̂^†(x) ϑ̂(x')⟩ and anomalous ⟨ϑ̂(x) ϑ̂(x')⟩ averages are not independent quantities. The relations between them at zero temperature can be obtained <cit.> using the fact that the HFB ground state is the quasiparticle vacuum. The relations do not depend explicitly on parameters of the Hamiltonian and take the form

    ∫ d x ⟨ϑ̂^†(x_1)ϑ̂^†(x)⟩⟨ϑ̂(x)
    ϑ̂(x_2)⟩ =⟨ϑ̂^†(x_1)ϑ̂(x_2)⟩
       ∫ d x ⟨ϑ̂^†
    +∫ d x ⟨ϑ̂^†(x_1)ϑ̂(x)⟩⟨ϑ̂^†(x)ϑ̂(x_2)⟩,


    ∫ d x ⟨ϑ̂^†(x)ϑ̂(x_1)⟩⟨ϑ̂(x)ϑ̂(x_2)⟩
    =∫ d x ⟨ϑ̂(x_1)ϑ̂(x)⟩
       ∫ d x ⟨ϑ̂^†(x_1)ϑ̂(x)⟩⟨ϑ̂(x)×⟨ϑ̂^†(x)ϑ̂(x_2)⟩.


Finally, the stationary HFB equations can be obtained by means of variation of E-μ N, given by Eqs. (<ref>)-(<ref>), under the conditions (<ref>) and (<ref>). The variational functions are ϕ(x), ϕ^*(x), ⟨ϑ̂(x)ϑ̂(x')⟩, ⟨ϑ̂^†(x)ϑ̂^†(x')⟩, and ⟨ϑ̂^†(x)ϑ̂(x')⟩. We do not write explicitly the resulting equations.

The full HBF scheme is rather cumbersome, and additional simplifications of Eq. (<ref>) are often used <cit.>. The simplest approximation is to neglect all the quantum correlators and leave only the wavefunction of the condensate. Then we arrive at the well-known GP functional with the help of the pseudopotential.



 §.§ Generalized Gross-Pitaevskii equations as a simplified version of the Hartee-Fock-Bogoliubov approximation



The GGP equations can be obtained <cit.> from the full HFB theory when the normal two-boson terms are completely neglected in Eq. (<ref>)

    ⟨Ψ̂^†(x_1)Ψ̂^†(x_2)Ψ̂(x_2)Ψ̂(x_1)⟩≃ϕ^*(x_1)ϕ(x_2)^*ϕ(x_2)ϕ(x_1)
       +ϕ(x_1)ϕ(x_2)⟨ϑ̂^†(x_1) ϑ̂^†(x_2)⟩
    +ϕ^*(x_1)ϕ^*(x_2)⟨ϑ̂(x_2) ϑ̂(x_1)⟩
       +⟨ϑ̂^†(x_1) ϑ̂^†(x_2)⟩⟨ϑ̂(x_2) ϑ̂(x_1)⟩.

Besides, if the condensate depletion is small, one can neglect the second term in the r.h.s of Eq. (<ref>), which is of the next order of the condensate depletion. Thus, we are left with the relation

    ⟨ϑ̂^†(x_1)ϑ̂(x_2)⟩≃∫ d x ⟨ϑ̂^†(x_1)ϑ̂^†(x)⟩⟨ϑ̂(x)ϑ̂(x_2)⟩.

In the approximation (<ref>), Eq. (<ref>) turns into identity.

Then the variational procedure yields <cit.> the GGP equations (<ref>)-(<ref>).



