We investigate the ground-state properties of N bosons with attractive zero-range interactions characterized by the scattering length a > 0 and confined to the surface of a sphere of radius R. We present the analytic solution of the problem for N = 2, mean-field analysis for N , and exact diffusion Monte Carlo results for intermediate N. For finite N, we observe a smooth crossover from the uniform state in the limit a / R 1 (weak attraction) to a localized state at small a/R (strong attraction). With increasing N, this crossover narrows down to a discontinuous transition from the uniform state to a soliton of size R / N. The two states are separated by an energy barrier, tunneling under which is exponentially suppressed at large N. The system behavior is marked by a peculiar competition between space-curvature effects and beyond-mean-field terms, both breaking the scaling invariance of a two-dimensional mean-field theory.

In two dimensions, N bosonic atoms with attractive point-like interactions form bound states with very peculiar properties. Addressing this problem by using the mean-field (MF) density functional theory characterized by the coupling constant g and normalization N leads to a solution called the Townes soliton, first studied in nonlinear optics1 and recently observed in cold-atom experiments.2,3 This state exists only for a specific value of g = g c = 4 π 2 / ( m N ln r ), with r = 8.567; a stronger (weaker) attraction yields a collapsing (expanding) state.1 For g = gc, the shape of the soliton is defined only up to an arbitrary scaling factor, and the total MF energy vanishes independent of the soliton size.4–7 

Hammer and Son8 have argued that the two-dimensional coupling constant depends logarithmically on the typical length, which breaks the MF degeneracy and fixes a specific soliton size. They predict that in the limit of large N, the soliton shrinks by r = 2.927 every time one adds another atom, while the energy BN scales as B N / B 2 r N. Here, B 2 = 4 2 e 2 γ / ( m a 2 ) is the dimer energy, and γ is the Euler constant. These predictions complement the exact few-body numerics performed over a few decades for N 4.9–14 Bazak and Petrov15 have calculated BN for N up to 26.

In this paper, we consider the problem of N bosons confined on a spherical surface and interacting via zero-range attraction. In the MF approximation on a flat surface, there is a critical coupling constant gc, which separates the collapsed ( g < g c) and uniform ( g > g c) ground states. However, the curved geometry breaks the MF scaling invariance. As a result, the uniform state remains metastable in a finite interval g d < g < g c, where the uniform and collapsed states are separated by an energy barrier. These are characteristics of a first-order phase transition, the system size taken as the order parameter and the point g = g d = 2 2 π / ( m N ) being the spinodal on the uniform phase side. Here, the barrier and the frequency of the lowest (dipole) excitation of the uniform condensate vanish. These results are obtained in the N limit, for which the classical Gross–Pitaevskii description is exact. However, due to the singular behavior of the collapsed state and the fact that it does not constitute a good starting point for a 1 / N 1 expansion, the problem becomes essentially quantum as soon as one passes from N = to finite N. We argue that in this case the discontinuous gas-collapse transition turns into an avoided crossing between the uniform state and the Townes soliton of size R / N, the tunneling between these states being exponentially suppressed at large N. To confirm these predictions, we solve the two-body problem analytically, perturbatively analyze the gas and the soliton solutions, and perform diffusion Monte Carlo (DMC) and variational Monte Carlo (VMC) analysis up to N = 128.

Our results indicate that the change from the smooth crossover to discontinuous transition (for all practical purposes) happens approximately at N 10, which is quantitatively large. We attribute this effect to the relatively weak breaking of the scaling invariance, which results in a quantitatively low energy barrier between the phases. This means that linear superpositions of the uniform and soliton states and macroscopic tunneling dynamics could be observable in this system for rather large atom numbers 10 N 20.

Note that the uniform-to-soliton transition for bosons on a sphere is very different from what happens with their one-dimensional counterpart. Indeed, for attractive one-dimensional bosons on a ring, this transition is continuous in both infinite-N and finite-N cases.16,17

Our study is inspired by the growing experimental interest in shell-shaped bosonic gases18–20 and by the goal of producing low-dimensional atomic devices whose curved geometry acts as a tunable degree of freedom for quantum simulation.21,22

We start our discussion with the MF analysis. Setting = m = R = 1, the Hamiltonian of the system reads
H ̂ = 1 2 d Ω ( ψ ̂ Ω L ̂ 2 ψ ̂ Ω + g ψ ̂ Ω ψ ̂ Ω ψ ̂ Ω ψ ̂ Ω ) ,
(1)
where ψ ̂ Ω is the operator creating a boson at position Ω = ( θ , φ ) on the sphere, g is the coupling constant, and L ̂ 2 = [ θ 2 + cot θ θ + ( 1 / sin 2 θ ) φ 2 ]. To regularize the local interaction term in Eq. (1), we use a cutoff for the interaction at an angular momentum l c 1 (equivalent to the momentum cutoff at κ = l c). The coupling constant is then taken as
g = 4 π / ln ( | B 2 | / κ 2 ) = 2 π / ln [ 2 e γ / ( κ a ) ] .
(2)
On the MF level, we approximate ψ ̂ Ω = ψ ̂ Ω = Φ ( Ω ) and minimize Eq. (1) with the normalization constraint d Ω | Φ ( Ω ) | 2 = N. In this manner, if we also assume the cylindrical symmetry (no dependence on the azimuthal angle φ), we arrive at the Gross–Pitaevskii equation
[ 1 2 ( θ 2 + cot θ θ ) + g | Φ ( θ ) | 2 μ ] Φ ( θ ) = 0.
(3)
The first obvious solution of Eq. (3) is the uniform condensate Φ = N / ( 4 π ) corresponding to the energy E N = g N 2 / ( 8 π ). The spectrum of the Bogoliubov modes in this case reads ε l , m = l ( l + 1 ) [ l ( l + 1 ) + g N / π ] / 2.23,24 The uniform solution is thus stable with respect to small amplitude modulations for g larger than g d = 2 π / N. Just below this point, the lowest-lying dipole mode (l = 1) becomes unstable. On the other hand, since gd < gc, in the interval g d < g < g c, the uniform state can only be metastable. Indeed, one can consider the Townes soliton solution with the angular width δ θ and show that its variational energy diverges as ( g g c ) / δ θ 2 in the asymptotic limit δ θ 0.

The conclusion that the ground state is, respectively, the collapsed state for g < g c and the uniform (gas) state for g > g c, is confirmed by our numerical analysis of Eq. (3) based on imaginary-time propagation. In practice, we see the collapsed state as a solution localized on a few sites of our spatial grid. This solution is nonuniversal since its size is determined by the finite grid spacing, and it shrinks to a delta-function distribution when increasing the density of the mesh. The imaginary-time propagation method does not predict other nonuniform solutions. On the other hand, by using a shooting method, which can in principle also describe maxima or saddle points of the energy functional, we do observe a universal (independent of the grid) nonuniform nodeless solution in the interval g d < g < g c. We identify it as a saddle point separating the uniform state from the collapsed one. Its energy is always larger than the energy of the uniform state. The maximal energy barrier Δ E max 0.0336 N is attained at g = gc. The gas–soliton transition at gc is thus discontinuous. The barrier decreases as g approaches gd and disappears completely at gd, which we identify as the spinodal point, consistent with the Bogoliubov analysis. We should note that the gas–soliton transition breaks the rotational symmetry on a sphere in the same manner as it breaks the translational symmetry in the flat case, and Eq. (3) explicitly takes this symmetry breaking into account.

The MF (or classical) Gross–Pitaevskii equation is the leading-order description of a weakly interacting Bose condensate. It holds when the number of atoms per healing volume is large (see, for instance, Ref. 25). In two dimensions, the healing volume scales as 1 / | μ | 1 / ( | g | | Φ | 2 ) and, therefore, it contains 1 / | g | atoms. The MF theory is thus valid when | g | 1. In our case, this inequality is equivalent to N 1, since we are interested in the region close to g c = 4 π / ( N ln r ) = 1.862 π / N.26 A conceptual problem when going from N = to finite N is caused by the singular MF description of the soliton state. Its size vanishes, and the binding energy diverges for g < g c, which is not a good starting point for a perturbation theory. This difficulty is resolved by the beyond-MF analysis of Hammer and Son8 who argue that the soliton has a finite size proportional to a (although with a prefactor exponentially small for large N) and finite energy 1 / a 2. Compared to the energy of the uniform state, proportional to g, the soliton energy varies exponentially fast, 1 / a 2 e 4 π / g. We thus conjecture that for large but finite N the gas–soliton transition is a crossing of these two states. Although not directly applicable to describe the crossing, the MF theory allows us to make two important statements. First, the MF description suggests that the barrier separating the two phases grows linearly with N, the crossing narrows down exponentially with N. Second, the MF results also suggest that for sufficiently large N the soliton at the crossing is small and is only weakly influenced by the local curvature. In the following Secs. III and IV, we will use the inequality N 1 to perturbatively describe, respectively, the uniform solution and the soliton solution. The transition point can then be located by comparing the corresponding energies. In Sec. VI, we will present our numerical results obtained for finite N.

In contrast to the soliton state, the uniform solution profits from a well-behaved MF description. Therefore, the MF formula E N = g N 2 / ( 8 π ) is sufficient for the asymptotic large-N analysis. Nevertheless, we would like to renormalize g and express it in a cutoff independent form. One way of performing this task is to note that the energy of a Bose gas in the weakly interacting limit is the interaction energy shift for a single pair multiplied by the number of pairs. In this approach, the renormalized g / ( 4 π ) can be identified with the energy E2 of a pair of atoms, which can be calculated directly from the two-body Schrödinger equation for the relative motion,
L ̂ 2 ψ ( θ ) = E 2 ψ ( θ ) .
(4)
The relative wave function ψ ( θ ) is regular everywhere except θ 0, where it satisfies the Bethe–Peierls boundary condition ψ ( θ ) ln ( θ / a ). The suitable solution is the associated Legendre function ψ ( θ ) P 1 / 2 + s [ cos ( π θ ) ], and the parameter s = E 2 + 1 / 4 satisfies
ln ( 1 / a ) = [ Ψ 0 ( 1 / 2 + s ) + Ψ 0 ( 1 / 2 s ) ] / 2 + ln ( e γ / 2 ) ,
(5)
where Ψ 0 is the digamma function. Equation (5) implicitly defines the function E 2 ( a ), which we find to smoothly connect the limits of strong ( a 1) and weak ( a 1) interactions (see Fig. 1). In the strongly interacting limit, we obtain the expansion E 2 = 4 exp ( 2 γ ) / a 2 1 / 3 + o ( 1 ), and, therefore, reproduce the flat-surface asymptotics E 2 B 2. For discussing the weakly interacting limit a 1, let us introduce
g r = 2 π ln ( 2 e 1 / 2 / a )
(6)
and expand Eq. (5) for small | E 2 |. In this manner, we obtain the series
E 2 = g r 4 π ( g r 4 π ) 3 + o ( g r 3 ) .
(7)
Thus, gr can be identified as the renormalized coupling constant suitable for describing the two-dimensional scattering in our particular geometry. Note that the energy can be equivalently expanded either in powers of g or in powers of gr, the corresponding series deviate from each other at the beyond-MF level since g g r g 2 ln κ | g |. We also note that the second-order term in Eq. (7) vanishes because of our choice of the constant under the logarithm in Eq. (6).
Fig. 1.

Inverse of the energy along the crossover from the weakly attractive (large negative x) to strongly attractive (large positive x) regimes. The cases of N = 2 and N = are shown as the black dashed-dotted curve and the piecewise linear black solid curve, respectively. The transition in the thermodynamic limit takes place at x = 0.171. The red, orange, pink, green, blue, and purple crosses are the Monte Carlo results, respectively, for N = 3 , 4 , 8 , 16 , 32, and 128 bosons. The dashed curves represent the large-x asymptote Eq. (15), which takes into account the leading-order curvature-induced energy shift for a nearly flat soliton.

Fig. 1.

Inverse of the energy along the crossover from the weakly attractive (large negative x) to strongly attractive (large positive x) regimes. The cases of N = 2 and N = are shown as the black dashed-dotted curve and the piecewise linear black solid curve, respectively. The transition in the thermodynamic limit takes place at x = 0.171. The red, orange, pink, green, blue, and purple crosses are the Monte Carlo results, respectively, for N = 3 , 4 , 8 , 16 , 32, and 128 bosons. The dashed curves represent the large-x asymptote Eq. (15), which takes into account the leading-order curvature-induced energy shift for a nearly flat soliton.

Close modal
Equation (7) is the energy of two atoms calculated up to terms g 3 g r 3. One can also generalize it to arbitrary N by using the standard perturbation expansion of the Hamiltonian (1) at small g. Up to the third order, this procedure gives (see the derivation in  Appendix A, based on Ref. 27)
E N = [ g r 4 π ( g r 4 π ) 3 ] N ( N 1 ) 2 + ( g r 4 π ) 3 N ( N 1 ) ( N 2 ) .
(8)
As expected, Eq. (8) reproduces the perturbative two-body result [Eq. (7)] and also predicts the leading nonpairwise contribution to the energy in the weakly interacting limit. Note that the three-body and, more generally, other beyond-MF terms, are suppressed by at least one additional power of g 1 / N compared to the leading two-body term E N = g r N 2 / ( 8 π ).

Let us assume (and a posteriori verify) that for large N there is a localized soliton solution of Eq. (1), degenerate with the uniform solution at a certain critical gr. By localized, we mean that the soliton size is much smaller than the sphere radius. In this nearly flat limit, we can address the problem perturbatively.

As a measure of the cloud size, it is convenient to introduce the rms separation between pairs of atoms r i j 2 . In the spherical geometry, r i j 2 is calculated from the distribution of the three-dimensional (chord) distances r i j = | r i r j |, where ri is the three-dimensional coordinate corresponding to the point Ωi on the sphere.

As we mention in the introduction, in the limit N , the size of the soliton on a plane equals8,
r i j 2 = a e N ( ln r ) / 2 + o ( N )
(9)
and its energy (to avoid confusion, we always denote energies on a sphere by E and on a plane by B)
B N = a 2 e N ln r + o ( N ) = e N ln r + 4 π / g r + o ( N ) ,
(10)
where the second equality in Eq. (10) follows from the definition of gr Eq. (6). One can show that
r i j 2 | B N | = 0.553 ,
(11)
which is a stronger statement than what can be obtained by simply multiplying Eqs. (9) and (10).28 We emphasize that B N 1 / r i j 2 is a beyond-MF energy scale. It is indeed much smaller than the MF kinetic or interaction energies, which both scale as N / r i j 2 , but have different signs and almost cancel each other. A detailed derivation of Eq. (11) is presented in  Appendix B, where we also review some properties of flat solitons.
Passing now to the curved geometry, one can show that the leading-order effect of the curved spherical surface on a nearly flat soliton is to shift its energy by Δ E sphere according to
E N = B N + Δ E sphere = B N 0.199 N .
(12)
To obtain Eq. (12), we use the flat Townes soliton as the unperturbed solution and the difference between the kinetic energy operator on a sphere and on a plane as the perturbation (refer to  Appendix C for more details). It is important to emphasize that for the validity of Eq. (12), the curvature shift Δ E sphere should be small compared to the MF energy scale N / r i j 2 and not to the beyond-MF energy scale BN. As a result, without invalidating Eq. (12), both terms on its right-hand side can be comparable to each other. As we will now see, at the transition point, the curvature effect indeed competes with the beyond-MF effect.
Combining the results of Sec. III on the energy of the uniform phase and Eq. (12) for the soliton, we obtain the condition for their crossing (valid for large N and 1 / r i j 2 ),
g r N 2 / ( 8 π ) + o ( N ) = B N + Δ E sphere + o ( N , 1 / r i j 2 ) .
(13)
Substituting Eq. (10) into Eq. (13), we find that the critical coupling constant equals g r = 4 π / ( N ln r ) + o ( N 1 ) g c. This conclusion just follows from the fact that at the crossing, BN given by Eq. (10) cannot be exponentially large in N as there are no other exponentially large terms in Eq. (13). We, thus, recover the MF result that the transition point corresponds to gr = gc. Substituting this equality into the left-hand side of Eq. (13) tells us that at the crossing B N = N / ( 2 ln r ) + 0.199 N and, using Eq. (11), we obtain the following estimate for the critical soliton size:
r i j 2 | crit = 4.06 / N .
(14)
In Eq. (13), we keep track of possible corrections; we include the beyond-MF correction o(N) for the uniform phase on the left-hand side and an estimate of higher-order beyond-MF and finite-curvature corrections for the soliton on the right-hand side. Using Eq. (14), one can a posteriori verify that these corrections are small.

We have thus confirmed that for large N the gas–soliton transition is a crossing of the uniform state, which occupies the whole sphere surface and a soliton of size 1 / N. In passing, we note that in order to observe a significant size change at the transition point the number of atoms should be N 16, as follows from Eq. (14). In Sec. V, we discuss what happens with the gas–soliton transition for finite N.

In Fig. 1, we show the energy obtained by solving the N-boson problem on a sphere by means of the diffusion and variational Monte Carlo methods, details of which will be discussed in Sec. VI. The dots in Fig. 1 stand for the DMC results for N = 3 (red), N = 4 (orange), N = 8 (pink), N = 16 (green), N = 32 (blue), and N = 128 (purple). The N = 2 result [Eq. (5)] is plotted as the black dashed-dotted curve.

For presenting the data, we choose to plot the quantity y = ( N 1 ) / ( 8 π E N ) as a function of x = 1 / ( g r N ) = ln ( 2 e 1 / 2 / a ) / ( 2 π N ). With this choice of coordinates, in the weakly interacting limit ( x ), the curves y(x) for different N should all tend to the straight black solid y = x line, equivalent to E N = g r N ( N 1 ) / ( 8 π ). We note that the convergence of the finite-N results to this line with increasing N is not monotonic and is relatively slow, which is very well explained by the cubic terms in Eq. (8). We do not show the corresponding curves to avoid cluttering.

In the opposite strongly interacting limit ( x + ), our numerical results agree with Eq. (12), which, in variables x and y, explicitly reads
y = N 1 8 π 1 ( B N / B 2 ) exp ( 4 π N x + 1 2 γ ) Δ E sphere ( N ) .
(15)
In fact, we calculate the curvature-induced offset Δ E sphere ( N ) for arbitrary N, not only for N 1. However, we find that Δ E sphere ( N ) converges very fast to the large-N asymptote Eq. (12). For N = 3, we numerically find Δ E sphere ( 3 ) = 0.60 ( 2 ), and the distinction is considerable only for N = 2, where Δ E sphere ( 2 ) = 1 / 3 (see Sec. III). The dashed curves in Fig. 1 show Eq. (15) for N = 3, 4, 8, and 16 using the ratios B N / B 2 calculated in Ref. 15.

We observe that with increasing N the energy curves in Fig. 1 tend to the piecewise linear function, y = x (gas phase) for gr < gc and y = 0 (collapsed state) for gr > gc, consistent with the MF description of Sec. II. The transition in the limit N is marked by the vertical line at x = ln r / ( 4 π ) = 0.171 (equivalent to gr = gc).

Let us now discuss how the cloud size changes as we vary the interaction strength. In Fig. 2, we show the rms separation between atoms as a function of the scattering length a (rescaled by an N-dependent coefficient). In the strongly interacting limit, we are dealing with a localized and almost flat soliton, the size of which is proportional to a [see Eq. (9)]. Accordingly, we fit this linear dependence and rescale the horizontal axis to make the curves for different N collapse to a single line. In the opposite weakly interacting limit, the distribution of atoms on our unit sphere is uniform, and the rms size tends to 2. The way these two asymptotes are approached depends on N. We find that for N 10 the derivatives of the curves are monotonic (within our precision), and for larger N, we start seeing a nonmonotonic structure, which eventually transforms into an abrupt change at a certain critical a. Our calculations are consistent with the scenario that there is a (narrow) crossing region where the exact ground state is a linear superposition of the soliton state and the uniform gas state. However, beginning with N 20, our DMC scheme does not account for these effects. For N = 32 and 128, the DMC scheme predicts the energies and rms sizes of the two phases, but due to the importance sampling (necessary for these large N), the system gets stuck in one of the phases and cannot tunnel to the other. Neglecting this macroscopic tunneling, one can say that the system undergoes a first-order phase transition (see the discussion in Sec. VI). We find that for N = 32 and, particularly, for N = 128 the rms size of the soliton at the transition point is in quantitative agreement with Eq. (14).

Fig. 2.

The rms separation versus a for different values of N as in Fig. 1 (we use the same color code). The interatomic distance is defined as the three-dimensional chord length. The horizontal axis is proportional to the scattering length a rescaled by an N-dependent coefficient to make all curves tend to a unique line at small a. The nonmonotonic behavior of the slope of the curves for N 16 signals the abruptly changing size in the vicinity of the gas–soliton transition.

Fig. 2.

The rms separation versus a for different values of N as in Fig. 1 (we use the same color code). The interatomic distance is defined as the three-dimensional chord length. The horizontal axis is proportional to the scattering length a rescaled by an N-dependent coefficient to make all curves tend to a unique line at small a. The nonmonotonic behavior of the slope of the curves for N 16 signals the abruptly changing size in the vicinity of the gas–soliton transition.

Close modal
To find the ground state of the N-boson problem numerically, we employ the variational (VMC) and diffusion (DMC) Monte Carlo methods, general details of which can be found, for instance, in Ref. 29. The standard method has to be adapted for solving the Schrödinger equation on a sphere. As the variational (for VMC) and guiding (for DMC) wave function, we take the Jastrow product
Ψ 0 ( Ω 1 , , Ω N ) = i < j χ ( r i j ) .
(16)
Expressing Ψ in terms of chord distances rij regularizes the wave function and automatically takes care of the boundary condition when two atoms are on the opposite poles of the sphere ( θ i j = π). For the Jastrow factor, we choose χ ( r ) = K 0 ( 2 e γ r / a ) e r / b. With this choice Ψ 0 satisfies the Bethe–Peierls condition at short interparticle distances ( θ i j r i j 0) as in the two-body case (see Sec. III), and we control the system size by tuning the variational parameter b.

In projection Monte Carlo methods, a first-order phase transition can be captured by imposing appropriate symmetry on the guiding wave functions. Probably, the most famous example of a zero-temperature quantum phase transition is that of the liquid–solid phase transition in 4He for which Monte Carlo calculations provide very good agreement with experiments.30 The distinguishing feature between the two phases is the order parameter, with the liquid phase having translational invariance, which is instead broken in the solid phase, where particles get localized near the lattice sites. On the variational level, close to the phase transition point, the variational energy has two minima as a function of the localization size. One corresponds to a solid (localization size smaller than the mean interparticle distance), while the other to a liquid (infinite localization size). The variational energy increases by going away from the minimum, so that each phase is stable against a weak perturbation, the phase with the lower energy being the ground state and the one with the higher energy being metastable. In contrast, in a second-order phase transition, only a single variational minimum exists, so that only one phase is stable at a time.

In our case, the variational parameter b describes the degree of localization of the state and allows one to discriminate between the two phases. For a sufficiently large number of particles ( N 20), indeed, we observe the double-minimum structure, characteristic of a first-order phase transition, see Fig. 3. The two minima are characterized by different values of b and correspond, respectively, to the gas phase and to the localized soliton phase. The double minimum is observed in the variation energy in a narrow region of x and the critical x (when the two minima are degenerate) approaches the MF prediction x 0.171 as the number of particles is increased.

Fig. 3.

The variational energy per particle E VMC / N as a function of the localization parameter b for N = 32 and three close values of x near the transition: for weak attraction (upper curve), there is a single minimum at large b 1 physically corresponding to a uniform gas state, for strong attraction (bottom curve), the minimum is located at significantly smaller b and describes a localized state, the middle curve shows the double minimum structure, typical for a first-order phase transition. The left and right insets show snapshots of particles' coordinates in VMC simulation in the minimum of the lower curve (marked by triangle, b = 14) and in the minimum of the upper curve (square, b = 44), respectively.

Fig. 3.

The variational energy per particle E VMC / N as a function of the localization parameter b for N = 32 and three close values of x near the transition: for weak attraction (upper curve), there is a single minimum at large b 1 physically corresponding to a uniform gas state, for strong attraction (bottom curve), the minimum is located at significantly smaller b and describes a localized state, the middle curve shows the double minimum structure, typical for a first-order phase transition. The left and right insets show snapshots of particles' coordinates in VMC simulation in the minimum of the lower curve (marked by triangle, b = 14) and in the minimum of the upper curve (square, b = 44), respectively.

Close modal

DMC method is based on solving the Schrödinger equation in imaginary time, and it allows one to find the ground-state properties in the limit of long projection time. The convergence is significantly improved by using a guiding wave function (16) with values of optimal b obtained by minimizing the variational energy. The procedure is in principle exact and is independent of the choice of the guiding function, as long as all relevant configurations are allowed. We find that for N = 32 and 128, [Eq. (16)], with the variationally optimal values of b, well describes both the gas and soliton solutions, and for these particle numbers, we do profit from importance sampling. However, since in the many-body configurational space, the “gas” guiding function to a large extent excludes the “soliton” region and vice versa, in practice, the diffusion process cannot account for the tunneling between the two phases. Quantitative description of the tunneling process in this case is an interesting future project going beyond the scope of this paper.

To mention a few technical details of our DMC scheme specific to the spherical geometry, we find it convenient to work with particles' three-dimensional coordinates ri rather than with their polar angles and azimuths. In particular, the chord distance evaluation in this case is much more straightforward. The drift for particle i is governed by the gradient of Ψ 0 with respect to ri projected to the sphere surface. For the Jastrow product (16), this task reduces to calculating the projected gradient of χ,
¯ r i χ ( | r i r j | ) = χ ( r i j ) ( r i r j ) r i r j r i j .
(17)
Neglecting finite time step corrections, we realize the diffusion and drift in a three-dimensional manner (as if the atoms were not confined), then projecting back to the sphere surface. The local energy can be expressed in terms of derivatives of χ by using
L ̂ r i 2 χ ( | r i r j | ) = 1 + r i r j 2 χ ( r i j ) + ( r i j 4 r i r j r i j ) χ ( r i j ) .
(18)

Finally, we mention that another Monte Carlo method, path integral Monte Carlo, has recently been employed to study supersolidity of bosons on a sphere with soft-core and dipolar interactions at finite temperature.31 

In this paper, we address the problem of N bosons on a spherical surface interacting via attractive zero-range potential and investigate how this system crosses over from the uniform gas state to the localized soliton state as one increases the attraction. Our main finding is that for sufficiently large N, we are essentially dealing with a first-order transition or, more precisely, a narrow avoided crossing between two states, which are significantly different in size and separated from each other by an energy barrier. On the other hand, for low values of N, the energy and rms size are smooth functions of the scattering length.

Our exact Monte Carlo calculations show that the change between these two regimes occurs when N is roughly between 10 and 20. These relatively high values may be explained by a quantitatively weak geometry-induced breaking of the two-dimensional MF scaling invariance and, therefore, by a numerically small barrier. Indeed, in the MF approximation at g = gc, the barrier equals Δ E max = 0.0336 N, and we can compare it with the dipole mode frequency in the uniform phase ε 1 , 1 = 1 g c / g d = 0.262. Here, we have in mind that the dipole mode marks the initial “direction” for the system to tunnel toward the soliton side.32 We see that for large N, the barrier is indeed much higher than the dipole frequency, the tunneling is suppressed, and the gas–soliton transition looks discontinuous (from the practical viewpoint). By contrast, when Δ E max ε 1 , 1, the role of the barrier is not very important, and the crossover is smooth. Accordingly, the condition Δ E max = ε 1 , 1 leads to a characteristic atom number N 8 where the crossover physics changes. This estimate is consistent with our Monte Carlo results.

We mention that the crossover scenario on a sphere is different from the previously solved model of attractive one-dimensional bosons on a ring where the gas–soliton transition is always continuous (no barrier), independent of N.16,17 In perspective, it would therefore be interesting to study other geometries where the confinement can be used as a knob for controlling the crossover type and for elucidating the relative role of the space curvature, finite size, manifold topology, MF, and beyond-MF terms. Being able to control the barrier and N gives a way to prepare and probe superpositions of different macroscopic or mesoscopic quantum states.

We acknowledge support from ANR grant “Droplets” No. ANR-19-CE30-0003-02 and from the EU Quantum Flagship (PASQuanS2.1, 101113690). G.E.A. acknowledges support by the Spanish Ministerio de Ciencia e Innovación (MCIN/AEI/10.13039/501100011033, Grant No. PID2020-113565GB-C21) and by the Generalitat de Catalunya (Grant No. 2021 SGR 01411). We also acknowledge Santander Supercomputacion support group at the University of Cantabria who provided access to the supercomputer Altamira Supercomputer at the Institute of Physics of Cantabria (IFCA-CSIC), member of the Spanish Supercomputing Network, for performing simulations/analyses. ICFO group acknowledges support from: European Research Council AdG NOQIA; MCIN/AEI [PGC2018-0910.13039/501100011033, CEX2019-000910-S/10.13039/501100011033, Plan National FIDEUA PID2019-106901GB-I00, Plan National STAMEENA PID2022-139099NB-I00, project funded by MCIN/AEI/10.13039/501100011033 and the “European Union NextGenerationEU/PRTR” (PRTR-C17.I1), FPI]; QUANTERA MAQS PCI2019-111828-2; QUANTERA DYNAMITE PCI2022-132919, QuantERA II Programme co-funded by European Union's Horizon 2020 program under Grant Agreement No. 101017733; Ministry for Digital Transformation and of Civil Service of the Spanish Government through the QUANTUM ENIA project call—Quantum Spain project, and by the European Union through the Recovery, Transformation and Resilience Plan—NextGenerationEU within the framework of the Digital Spain 2026 Agenda; Fundació Cellex; Fundació Mir-Puig; Generalitat de Catalunya (European Social Fund FEDER and CERCA program, AGAUR Grant No. 2021 SGR 01452, QuantumCAT U16-011424, co-funded by ERDF Operational Program of Catalonia 2014-2020); Barcelona Supercomputing Center MareNostrum (FI-2023-1-0013); funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union, European Commission, European Climate, Infrastructure and Environment Executive Agency (CINEA), or any other granting authority. Neither the European Union nor any granting authority can be held responsible for them (EU Quantum Flagship PASQuanS2.1, 101113690, EU Horizon 2020 FET-OPEN OPTOlogic, Grant No. 899794), EU Horizon Europe Program (This project has received funding from the European Union's Horizon Europe research and innovation program under Grant Agreement No. 101080086 NeQSTGrant Agreement No. 101080086—NeQST); ICFO Internal “QuantumGaudi” project; European Union's Horizon 2020 program under the Marie Sklodowska-Curie Grant Agreement No. 847648; “La Caixa” Junior Leaders fellowships, “La Caixa” Foundation (ID 100010434): CF/BQ/PR23/11980043.

The authors have no conflicts to disclose.

Andrea Tononi: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Grigori E. Astrakharchik: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Supervision (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Dmitry S. Petrov: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the authors upon reasonable request.

A complementary manner of renormalizing the interaction is to expand the N-body energy in powers of | g | 1 starting from the kinetic energy term in Eq. (1) as the unperturbed Hamiltonian. This standard perturbation theory works for finite N and also predicts the nonpairwise energy contribution. For N bosons on a unit sphere, we obtain the ground-state energy up to terms of order g3 in the form (see Ref. 27)
E N = [ g ( 1 ) + g ( 2 ) + g ( 3 ) ] N ( N 1 ) 2 + g 3 N ( N 1 ) ( N 2 ) 6 ,
(A1)
where
g ( 1 ) = V 00 00 = g 4 π ,
(A2)
g ( 2 ) = ν | V 0 ν 0 ν ¯ | 2 2 ξ ν = ( g 4 π ) 2 [ ln e 1 2 γ l c 2 + o ( 1 ) ] ,
(A3)
g ( 3 ) = ν ν V 0 ν 0 ν ¯ V ν ν ν ¯ ν ¯ V ν 0 ν ¯ 0 4 ξ ν ξ ν V 00 00 ν | V 0 ν 0 ν ¯ | 2 4 ξ ν 2 = ( g 4 π ) 3 [ ln 2 e 1 2 γ l c 2 1 + o ( 1 ) ] ,
(A4)
and
g 3 = 6 ν V 0 ν 0 ν ¯ V ν ¯ 0 0 ν ¯ V ν 0 ν ¯ 0 4 ξ ν 2 = 6 ( g 4 π ) 3 [ 1 + o ( 1 ) ] .
(A5)
In Eqs. (A3)–(A5), ν = { l , m } , ν = { l , m } , ν ¯ = { l , m } , ν ¯ = { l , m } , ν = l = 1 l c m = l l, and ξ ν = l ( l + 1 ) / 2, where l and m correspond to the angular momentum and its projection for spherical harmonics. The quantities
V μ ν ζ η = g d Ω ψ ζ * ( Ω ) ψ μ * ( Ω ) ψ η ( Ω ) ψ ν ( Ω )
(A6)
are the interaction matrix elements for two-body transitions from single-particle states μ and ζ to ν and η, and the corresponding wave functions are spherical harmonics. The last terms in Eqs. (A3)–(A5) are the final results obtained by summing over ν and ν . The summation is facilitated by the fact that V 0 ν 0 ν ¯ = V ν ¯ 0 0 ν ¯ = V ν 0 ν ¯ 0 = g / ( 4 π ). The more general matrix element V ν ν ν ¯ ν ¯ in the first sum in Eq. (A4) is not necessary to calculate explicitly. The result in this case is obtained by summing over the projections m and m and using the addition theorem for spherical harmonics. Finally, to arrive at Eq. (8) of the main text, we rewrite the expansion (A1) in powers of gr by using the relation 1 / g = 1 / g r + ln ( e γ + 1 / 2 / l c ) / ( 2 π ), which follows from Eqs. (2) and (6). Note that the cutoff dependence drops out.

From Eqs. (9) and (10) obtained by Hammer and Son,8 one can conclude that the product | B N | r i j 2 = exp [ o ( N ) ]. This does not exclude that this product may be a power of N since the term o(N) can, in principle, be logarithmic. Here, by developing the beyond-MF description of the flat soliton, we show that its size and energy are related by Eq. (11), i.e., the product | B N | r i j 2 tends to a universal constant.

The planar version of the Hamiltonian Eq. (1) reads
H ̂ plane = 1 2 d 2 ρ ( ψ ̂ ρ ρ 2 ψ ̂ ρ + g ψ ̂ ρ ψ ̂ ρ ψ ̂ ρ ψ ̂ ρ ) ,
(B1)
and the corresponding mean-field Gross–Pitaevskii energy functional is
B MF ( Ψ ) = 1 2 d 2 ρ [ | ρ Ψ ( ρ ) | 2 + g | Ψ ( ρ ) | 4 ] .
(B2)
Minimizing Eq. (B2) with the normalization constraint d 2 ρ | Ψ ( ρ ) | 2 = N leads to stationary solutions only for g = g c = 4 π / ( N ln r ).1 The corresponding wave functions form a family of self-similar states parametrized by the size σ,
Ψ σ ( ρ ) = N ln r 8 π σ 2 f ( ρ / σ ) ,
(B3)
where f ( ρ ) is a nodeless solution of
f ( ρ ) + f ( ρ ) / ρ f ( ρ ) + f 3 ( ρ ) = 0.
(B4)
The function f ( ρ ) is found numerically, it has a bell shape with exponentially decaying large-ρ asymptote.1,8 The normalization integral, the kinetic energy, and the interaction energy corresponding to the solution Ψ σ can be found from the identities
d ρ ρ f 2 ( ρ ) = d ρ ρ [ f ( ρ ) ] 2 = 1 2 d ρ ρ f 4 ( ρ ) = 4 ln r .
(B5)
The rms separation of atoms in the Townes soliton is related to the size parameter σ by
r i j 2 = σ 2 d ρ ρ 3 f 2 ( ρ ) / d ρ ρ f 2 ( ρ ) = 1.541 σ ,
(B6)
where we use
d ρ ρ 3 f 2 ( ρ ) = 2.211 .
(B7)
For future reference, we also provide the integral
d ρ ρ 3 [ f ( ρ ) ] 2 = 2.968 .
(B8)

One can check that if g = gc, the total energy B MF ( Ψ σ ) = 0 for any σ. This result, in a more general form, follows from the stationarity condition.4 Although the kinetic energy and the interaction energy in Eq. (B2) cancel each other, they are separately of order g N 2 / σ 2 N / σ 2, defining the MF energy scale. The leading beyond-MF contribution, which we denote B BMF, is smaller by the factor | g | 1 / N 1. The shape of Ψ is fixed by the dominant MF forces. The beyond-MF term is too weak to induce significant changes in the shape, but it breaks the degeneracy associated with the parameter σ. The optimal σ = σ N is determined by minimizing the energy functional including the beyond-MF term B BMF ( Ψ σ ). Calculating this quantity is a complicated task, which requires diagonalizing the Bogoliubov Hamiltonian for fluctuations around the stationary condensate solution Ψ σ. In this manner, one should, in principle, obtain the size σN (and, therefore, r i j 2 ) and the energy BN up to the preexponential factors [we are speaking about the terms o(N) in Eqs. (9) and (10)], eventually arriving at Eq. (11). We will now show that determining the product B N σ N 2 does not require these complicated calculations; it is sufficient to know the scaling properties of B BMF ( Ψ σ ).

In contrast to Eq. (B2), the Hamiltonian Eq. (B1) is not scaling invariant because of the (implicit) cutoff κ in the interaction. We can nevertheless use the following scaling property: a scaling transformation (change of coordinates by a factor Λ) of a many-body eigenstate of the Hamiltonian Eq. (B1) corresponding to a cutoff κ gives an eigenstate of the same Hamiltonian but with 1 / κ rescaled by Λ. This rescaling of the cutoff momentum is equivalent to rescaling the scattering length a by Λ [see Eq. (2)]. The corresponding energy gets rescaled by 1 / Λ 2. This means that the beyond-MF term B BMF ( Ψ σ , κ ) calculated with the cutoff κ for the condensate wave function Ψ σ and the analogous quantity B BMF ( Ψ σ , κ ) calculated with the cutoff κ = κ σ / σ for Ψ σ are related by
B BMF ( Ψ σ , κ ) = B BMF ( Ψ σ , κ ) ( σ / σ ) 2 .
(B9)
Let us now assume that we know B BMF ( Ψ σ , κ ) for a certain reference size σ , and we want to calculate B BMF ( Ψ σ , κ ), i.e., we want to know how the beyond-MF term scales with σ for a given fixed κ. In view of Eq. (B9), it is sufficient to calculate B BMF ( Ψ σ , κ ), which we write as B BMF ( Ψ σ , κ ) = B BMF ( Ψ σ , κ ) + Δ B. Here, Δ B = B BMF ( Ψ σ , κ ) B BMF ( Ψ σ , κ ) is the difference between the beyond-MF terms calculated for the same condensate wave function Ψ σ but different cutoff momenta. This is a local term since it takes into account excitations with de Broglie wave lengths between 1 / κ and 1 / κ , both assumed to be much smaller than the soliton size σ . When calculating Δ B, one can thus use the Bogoliubov theory for homogeneous condensates (and then average over the density) or just simply note that changing κ κ corresponds to a renormalization of the coupling constant given by the second-order Born integral. The result is [by o(1), we mean terms 1 in the asymptotic limit | g | 1 / N 0]
Δ B = d 2 ρ | Ψ σ ( ρ ) | 4 2 κ κ 2 π kdk ( 2 π ) 2 g 2 k 2 + o ( 1 ) = N 2 g 2 ln r ( 4 π σ ) 2 ln κ κ + o ( 1 ) = N 2 g 2 ln r ( 4 π σ ) 2 ln σ σ + o ( 1 ) .
(B10)
Using the fact that B MF ( Ψ σ ) = ( g g c ) d 2 ρ | Ψ σ ( ρ ) | 4 / 2 and Eqs. (B9) and (B10), we can write the energy of a soliton of size σ up to the leading beyond-MF correction as
B ( Ψ σ ) = N 2 ( g g c ) ln r 8 π σ 2 + B BMF ( Ψ σ , κ ) ( σ σ ) 2 N 2 g 2 ln r ( 4 π σ ) 2 ln σ σ .
(B11)
Minimizing Eq. (B11) with respect to σ leads to the result
B N σ N 2 = N 2 g 2 ln r 32 π 2 .
(B12)
The point is that BN and σN separately require calculating B BMF, whereas the product Eq. (B12) does not. Our theory implies that | g g c | | g |, so that the first term on the right-hand side of Eq. (B11) is consistent with the beyond-MF energy scale. Not to exceed the accuracy, we should then also replace g by gc in the last term in Eq. (B11) as well as in Eq. (B12). Equation (11) then follows from Eqs. (B12) and (B6).
To complete the proof of Eq. (11), we still have one technical point to address. Even when a is fixed, we can choose different g to model the interaction, provided that we tune κ in accordance with Eq. (2). This freedom is limited by two constraints: | g g c | | g |, discussed above, and κ 1 / σ N. The second constraint allows us to neglect finite-range effects and guarantees that our theory can be used to describe zero-range interactions. In fact, for the accuracy claimed in Eq. (B10), we need κ to exceed 1 / σ N by at least one power of N. The obvious choice g = gc satisfies the first constraint. It corresponds to κ c = 2 a 1 exp [ N ( ln r ) / 2 γ ] = σ N 1 exp [ o ( N ) ] as follows from Eqs. (2), (9), and (B6). Since o(N) is not necessarily asymptotically large, κc may not be larger than σ N 1, and we cannot ensure that the interaction is sufficiently short ranged. Let us then choose κ = σ N 1 N and show that the corresponding g satisfies | g g c | | g |. To this end, we write
g = 2 π ln ( 2 e γ / κ c a ) + ln ( κ c / κ ) g c g c 2 2 π ln κ c κ
(B13)
and note that κ c / κ = exp [ o ( N ) ], which means | g g c | = g 2 o ( N ) | g | as we need. This point completes the proof of Eq. (11).
Let us find how the energy of a nearly flat soliton is shifted under the influence of the sphere. The idea is to perturbatively calculate the curvature-induced MF energy shift assuming that σ 1 and using the flat-soliton wave function Eq. (B3) as the unperturbed solution. Consider a cylindrically symmetric function Ψ ( ρ ) where ρ = 2 sin ( θ / 2 ) is the chord distance between the soliton center, assumed to be at the north pole, and the point on the sphere with polar angle θ. In this case, we have d Ω = 2 π sin θ d θ = 2 π ρ d ρ, and the operator L ̂ 2 defined after Eq. (1) can be written as
L ̂ 2 = 2 ρ 2 1 ρ ρ + ρ 2 4 ( 2 ρ 2 + 3 ρ ρ ) .
(C1)
The MF energy functional on the sphere can thus formally be written as
E MF ( Ψ ) = B MF ( Ψ ) + d 2 ρ ρ 2 8 Ψ * ( ρ ) ( 2 ρ 2 + 3 ρ ρ ) Ψ ( ρ ) ,
(C2)
where B MF ( Ψ ) is given by Eq. (B2), and the domain of ρ is restricted to be between 0 and 2 with a proper boundary condition for Ψ at ρ = 2. However, the corresponding finite-size corrections are exponentially small for σ 1, and we can extend the domain up to ρ = . The leading-order curvature-induced energy shift is thus triggered by the integral in Eq. (C2) and can easily be calculated for Ψ = Ψ σ. Using Eq. (B8), we obtain
Δ E sphere ( Ψ σ ) = N ln r 32 [ f ( ρ ) ] 2 ρ 3 d ρ = 0.199 N ,
(C3)
which is Eq. (12). Although Eq. (C3) is formally a MF correction, it is by a factor of σ 2 weaker than the MF energy scale N / σ 2. It becomes comparable to the beyond-MF scale if σ 2 1 / N, which is exactly what happens at the transition point, as we mention in the main text. Note that Eq. (C3) is independent of σ, and thus it only shifts the energy, the optimal σ still being determined by the beyond-MF energy term Eq. (B11) derived in the flat-surface approximation.
1.
R. Y.
Chiao
,
E.
Garmire
, and
C. H.
Townes
,
Phys. Rev. Lett.
13
,
479
(
1964
).
2.
B.
Bakkali-Hassani
,
C.
Maury
,
Y.-Q.
Zou
,
É.
Le Cerf
,
R.
Saint-Jalm
,
P. C. M.
Castilho
,
S.
Nascimbene
,
J.
Dalibard
, and
J.
Beugnon
,
Phys. Rev. Lett.
127
,
023603
(
2021
).
3.
C.-A.
Chen
and
C.-L.
Hung
,
Phys. Rev. Lett.
127
,
023604
(
2021
).
4.
S. N.
Vlasov
,
V. A.
Petrishchev
, and
V. I.
Talanov
,
Izv. Vyssh. Uchebn. Zaved., Radiofiz.
14
,
1353
(
1971
).
5.
L. P.
Pitaevskii
,
Phys. Lett. A
221
,
14
(
1996
).
6.
L. P.
Pitaevskii
and
A.
Rosch
,
Phys. Rev. A
55
,
R853
(
1997
).
7.
B.
Bakkali-Hassani
and
J.
Dalibard
, “
Townes soliton and beyond: Non-miscible Bose mixtures in 2D
,” in
Proceedings of the International School of Physics “Enrico Fermi,” Course 211 - Quantum Mixtures with Ultra-Cold Atoms
,
2022
.
8.
H.-W.
Hammer
and
D. T.
Son
,
Phys. Rev. Lett.
93
,
250408
(
2004
).
9.
L. W.
Bruch
and
J. A.
Tjon
,
Phys. Rev. A
19
,
425
(
1979
).
10.
S. K.
Adhikari
,
A.
Delfino
,
T.
Frederico
,
I. D.
Goldman
, and
L.
Tomio
,
Phys. Rev. A
37
,
3666
(
1988
).
11.
E.
Nielsen
,
D. V.
Fedorov
, and
A. S.
Jensen
,
Phys. Rev. A
56
,
3287
(
1997
).
12.
E.
Nielsen
,
D. V.
Fedorov
, and
A. S.
Jensen
,
Few-Body Syst.
27
,
15
(
1999
).
13.
I. V.
Brodsky
,
M. Y.
Kagan
,
A. V.
Klaptsov
,
R.
Combescot
, and
X.
Leyronas
,
Phys. Rev. A
73
,
032724
(
2006
).
14.
O. I.
Kartavtsev
and
A. V.
Malykh
,
Phys. Rev. A
74
,
042506
(
2006
).
15.
B.
Bazak
and
D. S.
Petrov
,
New J. Phys.
20
,
023045
(
2018
).
16.
L. D.
Carr
,
C. W.
Clark
, and
W. P.
Reinhardt
,
Phys. Rev. A
62
,
063611
(
2000
).
17.
R.
Kanamoto
,
H.
Saito
, and
M.
Ueda
,
Phys. Rev. A
67
,
013608
(
2003
).
18.
R. A.
Carollo
,
D. C.
Aveline
,
B.
Rhyno
,
S.
Vishveshwara
,
C.
Lannert
,
J. D.
Murphree
,
E. R.
Elliott
,
J. R.
Williams
,
R. J.
Thompson
et al,
Nature
606
,
281
(
2022
).
19.
F.
Jia
,
Z.
Huang
,
L.
Qiu
,
R.
Zhou
,
Y.
Yan
, and
D.
Wang
,
Phys. Rev. Lett.
129
,
243402
(
2022
).
20.
Y.
Guo
,
E.
Mercado Gutierrez
,
D.
Rey
,
T.
Badr
,
A.
Perrin
,
L.
Longchambon
,
V. S.
Bagnato
,
H.
Perrin
, and
R.
Dubessy
,
New J. Phys.
24
,
093040
(
2022
).
21.
A.
Tononi
and
L.
Salasnich
,
Nat. Rev. Phys.
5
,
398
(
2023
).
22.
L.
Amico
,
D.
Anderson
,
M.
Boshier
,
J.-P.
Brantut
,
L.-C.
Kwek
,
A.
Minguzzi
, and
W.
von Klitzing
,
Rev. Mod. Phys.
94
,
041001
(
2022
).
23.
S.
Prestipino
and
P. V.
Giaquinta
,
Phys. Rev. A
99
,
063619
(
2019
).
24.
A.
Tononi
and
L.
Salasnich
,
Phys. Rev. Lett.
123
,
160403
(
2019
).
25.
C.
Mora
and
Y.
Castin
,
Phys. Rev. A
67
,
053615
(
2003
).
26.
We note that small g 1 / N implies exponentially large a. This peculiarity of the two-dimensional scattering is usually not a problem since most interaction-related physical observables depend on ln a. Exponentially large a is thus not unusual. For instance, for low-energy quasi-two-dimensional collisions, the effective two-dimensional scattering length a l 0 exp ( π / 2 l 0 / a 3 D ) is exponentially large when the ratio of the three-dimensional scattering length a 3 D < 0 to the confinement oscillator length l0 is (moderately) small, see
D. S.
Petrov
and
G. V.
Shlyapnikov
,
Phys. Rev. A
64
,
012706
(
2001
).
27.
A.
Pricoupenko
and
D. S.
Petrov
,
Phys. Rev. A
103
,
033326
(
2021
).
28.
References 8 and 15 implicitly assume without verification that the terms o(N) in Eqs. (9) and (10) are constants. However, the behavior o ( N ) = ln N cannot be excluded, making Eq. (11) nontrivial. We thank one of the reviewers for raising the concern.
29.
J.
Boronat
and
J.
Casulleras
,
Phys. Rev. B
49
,
8920
(
1994
).
30.
P. A.
Whitlock
,
D. M.
Ceperley
,
G. V.
Chester
, and
M. H.
Kalos
,
Phys. Rev. B
19
,
5598
(
1979
).
31.
M.
Ciardi
,
F.
Cinti
,
G.
Pellicane
, and
S.
Prestipino
,
Phys. Rev. Lett.
132
,
026001
(
2024
).
32.
A possible way of solving this tunneling problem is by using the instanton approach, see
J. A.
Freire
and
D. P.
Arovas
,
Phys. Rev. A
59
,
1461
(
1999
).