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 \u2192 \u221e$, 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 \u226b 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 $ \u223c 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.

## I. INTRODUCTION

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 optics^{1} and recently observed in cold-atom experiments.^{2,3} This state exists only for a specific value of $ g = g c = \u2212 4 \pi \u210f 2 / ( m N \u2009 ln \u2009 r )$, with *r* = 8.567; a stronger (weaker) attraction yields a collapsing (expanding) state.^{1} For *g* = *g _{c}*, 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 Son^{8} 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 *B _{N}* scales as $ B N / B 2 \u221d r N$. Here, $ B 2 = \u2212 4 \u210f 2 e \u2212 2 \gamma / ( 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 \u2264 4$.

^{9–14}Bazak and Petrov

^{15}have calculated

*B*for

_{N}*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 *g _{c}*, 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 = \u2212 2 \u210f 2 \pi / ( 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 \u2192 \u221e$ 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 \u226a 1$ expansion, the problem becomes essentially quantum as soon as one passes from $ N = \u221e$ 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 $ \u223c 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 \u223c 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 \u2272 N \u2272 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}

## II. MEAN-FIELD ANALYSIS IN THE LARGE-*N* LIMIT

*g*is the coupling constant, and $ L \u0302 2 = \u2212 [ \u2202 \theta 2 + cot \u2009 \theta \u2009 \u2202 \theta + ( 1 / \u2009 sin 2 \theta ) \u2202 \phi 2 ]$. To regularize the local interaction term in Eq. (1), we use a cutoff for the interaction at an angular momentum $ l c \u226b 1$ (equivalent to the momentum cutoff at $ \kappa = l c$). The coupling constant is then taken as

^{23,24}The uniform solution is thus stable with respect to small amplitude modulations for

*g*larger than $ g d = \u2212 2 \pi / N$. Just below this point, the lowest-lying dipole mode (

*l*= 1) becomes unstable. On the other hand, since

*g*<

_{d}*g*, 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 $ \delta \theta $ and show that its variational energy diverges as $ \u223c ( g \u2212 g c ) / \delta \theta 2$ in the asymptotic limit $ \delta \theta \u2192 0$.

_{c}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 $ \Delta E max \u2248 0.0336 N$ is attained at *g* = *g _{c}*. The gas–soliton transition at

*g*is thus discontinuous. The barrier decreases as

_{c}*g*approaches

*g*and disappears completely at

_{d}*g*, 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.

_{d}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 $ \u223c 1 / | \mu | \u223c 1 / ( | g | | \Phi | 2 )$ and, therefore, it contains $ \u223c 1 / | g |$ atoms. The MF theory is thus valid when $ | g | \u226a 1$. In our case, this inequality is equivalent to $ N \u226b 1$, since we are interested in the region close to $ g c = \u2212 4 \pi / ( N \u2009 ln \u2009 r ) = \u2212 1.862 \pi / N$.^{26} A conceptual problem when going from $ N = \u221e$ 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 Son^{8} who argue that the soliton has a finite size proportional to *a* (although with a prefactor exponentially small for large *N*) and finite energy $ \u221d 1 / a 2$. Compared to the energy of the uniform state, proportional to *g*, the soliton energy varies exponentially fast, $ \u221d 1 / a 2 \u221d e 4 \pi / 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 \u226b 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*.

## III. RENORMALIZATION OF THE INTERACTION AND THE TWO-BODY PROBLEM

*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 \pi )$ can be identified with the energy

*E*

_{2}of a pair of atoms, which can be calculated directly from the two-body Schrödinger equation for the relative motion,

*g*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

_{r}*g*or in powers of

*g*, the corresponding series deviate from each other at the beyond-MF level since $ g \u2212 g r \u223c g 2 \u2009 ln \u2009 \kappa \u226a | 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).

_{r}*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)

## IV. NEARLY FLAT SOLITON AND THE TRANSITION POINT

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 *g _{r}*. 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 $ \u27e8 r i j 2 \u27e9$. In the spherical geometry, $ \u27e8 r i j 2 \u27e9$ is calculated from the distribution of the three-dimensional (chord) distances $ r i j = | r i \u2212 r j |$, where **r**_{i} is the three-dimensional coordinate corresponding to the point Ω_{i} on the sphere.

^{8}

^{,}

*E*and on a plane by

*B*)

*g*Eq. (6). One can show that

_{r}^{28}We emphasize that $ B N \u223c 1 / \u27e8 r i j 2 \u27e9$ is a beyond-MF energy scale. It is indeed much smaller than the MF kinetic or interaction energies, which both scale as $ \u223c N / \u27e8 r i j 2 \u27e9$, 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.

*B*. 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.

_{N}*N*and $ 1 / \u27e8 r i j 2 \u27e9$),

*B*given by Eq. (10) cannot be exponentially large in

_{N}*N*as there are no other exponentially large terms in Eq. (13). We, thus, recover the MF result that the transition point corresponds to

*g*=

_{r}*g*. Substituting this equality into the left-hand side of Eq. (13) tells us that at the crossing $ B N = \u2212 N / ( 2 \u2009 ln \u2009 r ) + 0.199 N$ and, using Eq. (11), we obtain the following estimate for the critical soliton size:

_{c}*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 $ \u223c 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 \u2273 16$, as follows from Eq. (14). In Sec. V, we discuss what happens with the gas–soliton transition for finite *N*.

## V. NUMERICAL RESULTS

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 \u2212 1 ) / ( 8 \pi E N )$ as a function of $ x = 1 / ( g r N ) = ln \u2009 ( 2 e \u2212 1 / 2 / a ) / ( 2 \pi N )$. With this choice of coordinates, in the weakly interacting limit ( $ x \u2192 \u2212 \u221e$), 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 \u2212 1 ) / ( 8 \pi )$. 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.

*x*and

*y*, explicitly reads

*N*, not only for $ N \u226b 1$. However, we find that $ \Delta E sphere ( N )$ converges very fast to the large-

*N*asymptote Eq. (12). For

*N*= 3, we numerically find $ \Delta E sphere ( 3 ) = 0.60 ( 2 )$, and the distinction is considerable only for

*N*= 2, where $ \Delta E sphere ( 2 ) = \u2212 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 *g _{r}* <

*g*and

_{c}*y*= 0 (collapsed state) for

*g*>

_{r}*g*, consistent with the MF description of Sec. II. The transition in the limit $ N \u2192 \u221e$ is marked by the vertical line at $ x = \u2212 ln \u2009 r / ( 4 \pi ) = \u2212 0.171$ (equivalent to

_{c}*g*=

_{r}*g*).

_{c}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 \u2272 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 \u223c 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).

## VI. MONTE CARLO FOR BOSONS ON A SPHERE

*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

*r*regularizes the wave function and automatically takes care of the boundary condition when two atoms are on the opposite poles of the sphere ( $ \theta i j = \pi $). For the Jastrow factor, we choose $ \chi ( r ) = K 0 ( 2 e \u2212 \gamma r / a ) e \u2212 r / b$. With this choice $ \Psi 0$ satisfies the Bethe–Peierls condition at short interparticle distances ( $ \theta i j \u2248 r i j \u2192 0$) as in the two-body case (see Sec. III), and we control the system size by tuning the variational parameter

_{ij}*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 ^{4}He 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 \u2273 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 \u2192 \u2212 0.171$ as the number of particles is increased.

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.

**r**

_{i}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 $ \Psi 0$ with respect to

**r**

_{i}projected to the sphere surface. For the Jastrow product (16), this task reduces to calculating the projected gradient of

*χ*,

*χ*by using

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}

## VII. CONCLUSIONS

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* = *g _{c}*, the barrier equals $ \Delta E max = 0.0336 N$, and we can compare it with the dipole mode frequency in the uniform phase $ \epsilon 1 , 1 = 1 \u2212 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 $ \Delta E max \u2272 \epsilon 1 , 1$, the role of the barrier is not very important, and the crossover is smooth. Accordingly, the condition $ \Delta E max = \epsilon 1 , 1$ leads to a characteristic atom number $ N \u2248 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.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX A: DERIVATION OF EQ. (8)

*N*-body energy in powers of $ | g | \u226a 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

*g*

^{3}in the form (see Ref. 27)

*l*and

*m*correspond to the angular momentum and its projection for spherical harmonics. The quantities

*μ*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 $ \nu \u2032$. The summation is facilitated by the fact that $ V 0 \nu 0 \nu \xaf = V \nu \xaf 0 0 \nu \xaf = V \nu 0 \nu \xaf 0 = g / ( 4 \pi )$. The more general matrix element $ V \nu \u2032 \nu \nu \xaf \u2032 \nu \xaf$ 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 \u2032$ 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

*g*by using the relation $ 1 / g = 1 / g r + ln \u2009 ( e \u2212 \gamma + 1 / 2 / l c ) / ( 2 \pi )$, which follows from Eqs. (2) and (6). Note that the cutoff dependence drops out.

_{r}### APPENDIX B: DERIVATION OF EQ. (11)

From Eqs. (9) and (10) obtained by Hammer and Son,^{8} one can conclude that the product $ | B N | \u27e8 r i j 2 \u27e9 = 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 | \u27e8 r i j 2 \u27e9$ tends to a universal constant.

^{1}The corresponding wave functions form a family of self-similar states parametrized by the size

*σ*,

*ρ*asymptote.

^{1,8}The normalization integral, the kinetic energy, and the interaction energy corresponding to the solution $ \Psi \sigma $ can be found from the identities

*σ*by

One can check that if *g* = *g _{c}*, the total energy $ B MF ( \Psi \sigma ) = 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 / \sigma 2 \u223c N / \sigma 2$, defining the MF energy scale. The leading beyond-MF contribution, which we denote $ B BMF$, is smaller by the factor $ | g | \u223c 1 / N \u226a 1$. The shape of $\Psi $ 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 $ \sigma = \sigma N$ is determined by minimizing the energy functional including the beyond-MF term $ B BMF ( \Psi \sigma )$. Calculating this quantity is a complicated task, which requires diagonalizing the Bogoliubov Hamiltonian for fluctuations around the stationary condensate solution $ \Psi \sigma $. In this manner, one should, in principle, obtain the size

*σ*(and, therefore, $ \u27e8 r i j 2 \u27e9$) and the energy

_{N}*B*up to the preexponential factors [we are speaking about the terms

_{N}*o*(

*N*) in Eqs. (9) and (10)], eventually arriving at Eq. (11). We will now show that determining the product $ B N \sigma N 2$ does not require these complicated calculations; it is sufficient to know the scaling properties of $ B BMF ( \Psi \sigma )$.

*κ*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 / \kappa $ 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 / \Lambda 2$. This means that the beyond-MF term $ B BMF ( \Psi \sigma , \kappa )$ calculated with the cutoff

*κ*for the condensate wave function $ \Psi \sigma $ and the analogous quantity $ B BMF ( \Psi \sigma \u2032 , \kappa \u2032 )$ calculated with the cutoff $ \kappa \u2032 = \kappa \sigma / \sigma \u2032$ for $ \Psi \sigma \u2032$ are related by

*σ*for a given fixed

*κ*. In view of Eq. (B9), it is sufficient to calculate $ B BMF ( \Psi \sigma \u2032 , \kappa \u2032 )$, which we write as $ B BMF ( \Psi \sigma \u2032 , \kappa \u2032 ) = B BMF ( \Psi \sigma \u2032 , \kappa ) + \Delta B$. Here, $ \Delta B = B BMF ( \Psi \sigma \u2032 , \kappa \u2032 ) \u2212 B BMF ( \Psi \sigma \u2032 , \kappa )$ is the difference between the beyond-MF terms calculated for the same condensate wave function $ \Psi \sigma \u2032$ but different cutoff momenta. This is a local term since it takes into account excitations with de Broglie wave lengths between $ 1 / \kappa $ and $ 1 / \kappa \u2032$, both assumed to be much smaller than the soliton size $ \sigma \u2032$. When calculating $ \Delta B$, one can thus use the Bogoliubov theory for homogeneous condensates (and then average over the density) or just simply note that changing $ \kappa \u2192 \kappa \u2032$ corresponds to a renormalization of the coupling constant given by the second-order Born integral. The result is [by

*o*(1), we mean terms $ \u226a 1$ in the asymptotic limit $ | g | \u223c 1 / N \u2192 0$]

*σ*up to the leading beyond-MF correction as

*σ*leads to the result

*B*and

_{N}*σ*separately require calculating $ B BMF$, whereas the product Eq. (B12) does not. Our theory implies that $ | g \u2212 g c | \u226a | 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

_{N}*g*by

*g*in the last term in Eq. (B11) as well as in Eq. (B12). Equation (11) then follows from Eqs. (B12) and (B6).

_{c}*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 \u2212 g c | \u226a | g |$, discussed above, and $ \kappa \u226b 1 / \sigma 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 / \sigma N$ by at least one power of

*N*. The obvious choice

*g*=

*g*satisfies the first constraint. It corresponds to $ \kappa c = 2 a \u2212 1 \u2009 exp [ N ( ln \u2009 r ) / 2 \u2212 \gamma ] = \sigma N \u2212 1 \u2009 exp [ o ( N ) ]$ as follows from Eqs. (2), (9), and (B6). Since

_{c}*o*(

*N*) is not necessarily asymptotically large,

*κ*may not be larger than $ \sigma N \u2212 1$, and we cannot ensure that the interaction is sufficiently short ranged. Let us then choose $ \kappa = \sigma N \u2212 1 N$ and show that the corresponding

_{c}*g*satisfies $ | g \u2212 g c | \u226a | g |$. To this end, we write

### APPENDIX C: DERIVATION OF EQ. (12)

*θ*. In this case, we have $ d \Omega = 2 \pi \u2009 sin \u2009 \theta d \theta = 2 \pi \rho d \rho $, and the operator $ L \u0302 2$ defined after Eq. (1) can be written as

*ρ*is restricted to be between 0 and 2 with a proper boundary condition for $\Psi $ at

*ρ*= 2. However, the corresponding finite-size corrections are exponentially small for $ \sigma \u226a 1$, and we can extend the domain up to $ \rho = \u221e$. The leading-order curvature-induced energy shift is thus triggered by the integral in Eq. (C2) and can easily be calculated for $ \Psi = \Psi \sigma $. Using Eq. (B8), we obtain

*σ*, 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.

## REFERENCES

*a*. This peculiarity of the two-dimensional scattering is usually not a problem since most interaction-related physical observables depend on $\nln\n\u2009\na$. Exponentially large

*a*is thus not unusual. For instance, for low-energy quasi-two-dimensional collisions, the effective two-dimensional scattering length $\na\n\u221d\n\nl\n0\n\u2009\nexp\n(\n\u2212\n\n\n\pi \n/\n2\n\nl\n0\n/\n\na\n\n3\nD\n)$ is exponentially large when the ratio of the three-dimensional scattering length $\n\na\n\n3\nD\n<\n0$ to the confinement oscillator length

*l*

_{0}is (moderately) small, see