Suppression of quantum-mechanical collapse in bosonic gases with intrinsic repulsion: A brief review

It is well known that attractive potential which is inversely proportional to the squared distance from the origin gives rise to the critical quantum collapse in the framework of the three-dimensional (3D) linear Schroedinger equation. This article summarizes theoretical analysis, chiefly published in several original papers, which demonstrates suppression of the collapse caused by this potential, and the creation of the otherwise missing ground state in a 3D gas of bosonic dipoles pulled by the same potential to the central charge, with repulsive contact interactions between them, represented by the cubic term in the respective Gross-Pitaevskii equation (GPE). In two dimensions (2D), quintic self-repulsion is necessary for the suppression of the collapse; alternatively, this may be provided by the effective quartic repulsion, produced by the Lee-Huang-Yang correction to the GPE. 3D states carrying angular momentum are constructed in the model with the symmetry reduced from spherical to cylindrical by an external polarizing field. Interplay of the collapse suppression and miscibility-immiscibility transition is considered in a binary condensate. The consideration of the 3D setting in the form of the many-body quantum system, with the help of the Monte Carlo method, demonstrates that, although the quantum collapse cannot be fully suppressed, the self-trapped states, predicted by the GPE, exist in the many-body setting as metastable modes protected against the collapse by a tall potential barrier.


I. INTRODUCTION
One of standard exercises given to students taking the course of quantum mechanics is solving the three-dimensional (3D) Schrödinger equation with an isotropic attractive potential, [1]. This exercise offers a unique example of critical phenomena in the nonrelativistic quantum theory. Indeed, the corresponding classical (Newton's) equation of motion for the particle's coordinates, r = {x, y, z}, admits obvious rescaling, t ≡t/ √ U 0 , which eliminates U 0 from Eq. (2), thus making the solution invariant with respect to the choice of a positive value of the potential strength, U 0 . However, the invariance is lost in the corresponding 3D Schrödinger equation for wave function ψ (r, t), in which U 0 cannot be removed by rescaling. This drastic difference between the classical mechanical system and its quantum-mechanical counterpart is known as the quantum anomaly , alias "dimensional transmutation" [2,3]. The consequence of the anomaly is well known: if an external trapping potential, is added to Eq. (3), to make the integral norm, convergent at r → ∞, the Schrödinger equation gives rise to the normal set of trapped modes, starting from the ground state (GS), at On the other hand, above this critical point, i.e., at U 0 > 1/4, the GS does not exist (or, formally, speaking, it has an infinitely small size, corresponding to energy E → −∞, which is known as "fall onto the center" [1], the other name for which is "quantum collapse" [2,3]).
In the 2D space, the quantum collapse, driven by the same potential (1), is more violent, taking place at any value U 0 > 0 (in other words, the respective critical value is (U 0 ) (2D) cr = 0). Finally, in the 1D case the same potential (1) gives rise to a still stronger superselection effect, which means splitting the 1D space into two non-communicating subspaces, x ≷ 0 [4].
A solution to the quantum-collapse problem in the 3D case was proposed in terms of a linear quantum-field-theory, replacing the usual quantum-mechanical wave function by the secondary-quantized field [2,3]. This approach makes it possible to introduce the GS, which is missing at U 0 > 1/4 in the framework of standard quantum mechanics. However, the solution does not predict a definite value of the size of the newly created GS. Instead, the field-theory formulation, based on the renormalization-group technique, introduces a GS with an arbitrary spatial scale, in terms of which all other spatial sizes are measured in that setting.
The present mini-review aims to summarize results produced by works which elaborated another possibility to resolve the problem of the quantum collapse. This possibility was proposed in Ref. [5], and then developed, for more general situations, in works [6] and [7]. The solution was based on the consideration of an ultracold gas of bosonic particles pulled to the center by potential (1). The gas was assumed to be in the state of the Bose-Einstein condensate (BEC) [8], and the suppression of the single-particle quantum collapse in this coherent many-body setting was provided by repulsive contact interactions between colliding particles in the gas. The solution was elaborated in the framework of the mean-field approach [8], i.e., treating the single-particle wave function, which represents all particles in the gas, as a classical field governed by the corresponding Gross-Pitaevskii equation (GPE).
The same work [5] offered a physical realization of potential (1) in the 3D space, which was previously considered as a formal exercise [1]. The realization is provided by assuming that the bosonic particles are small molecules carrying a permanent electric dipole moment, d, pulled by the electrostatic force to a point-like charge, Q, placed at the origin, which creates electric field E = Qr/r 3 . In this connection, it is relevant to mention that it has been demonstrated experimentally that a free charge (ion), immersed in an ultracold gas, may be kept at a fixed position by means of a laser-trapping technique [9]. Assuming that the orientation of the dipole carried by each particle is locked to the local field, i.e., d/d = sgn(Q) (r/r), so as to minimize the interaction energy, the respective interaction potential is U (r) = −d · E, which is tantamount to potential (1) with strength As for the dipolar molecules which may be used to build the BEC under the consideration, experimental results suggest that they may be, e.g., LiCs [10] or KRb [11]. The gas of ultracold dipolar molecules, trapped in a pancake-shaped configuration shaped by an appropriate external potential [12], with the central electric charge immersed in the gas as outlined above, provides for the realization of the 2D version of the setting. An alternative realization of the 2D setting is offered by a gas of polarizable atoms without a permanent dielectric moment, while an effective moment is induced in them by the electric field of a uniformly charged wire set perpendicular to the pancake's plane [13], or with an effective magnetic moment induced by a current filament (e.g., an electron beam) piercing the pancake perpendicularly.
In the context of 2D settings, it is relevant to mention that a quantum anomaly was also predicted in a model described by the GPE in the 2D space, for a gas of bosons with the repulsive contact interaction, trapped in the harmonic-oscillator potential (4) [14]. The anomaly breaks the specific scaling invariance of this gas, which holds in the mean-field approximation.
In terms of the GPE, the contact repulsive interaction in the bosonic gas is represented by the cubic term [8]. With the addition of this term, and taking into regard the external trapping potential (4), which is present in any experiment with ultracold atoms, the linear Schrödinger equation (3) is replaced by the GPE, which is written here in the scaled form: It is relevant to mention that the 3D GPE with the self-attractive interaction, which corresponds to the opposite sign in front of the cubic term in Eq. (8), gives rise, in the absence of the attractive potential (U 0 = 0) to the well-known supercritical wave collapse [17]. A relation of this setting to Eq. (8) is that the inclusion of the trapping potential ∼ Ω 2 gives rise to stable bound states in the form of spherically symmetric bound states and ones with vorticity m = 1 (cf. Eq. (16) below), provided that norm N does not exceed a certain critical value [18]- [21].
The energy (Hamiltonian) corresponding to Eq. (8) is The scaled variables and constants, in terms of which Eq. (8) is written, are related to their counterparts measured in physical units: where m and a s are the bosonic mass and s-scattering length, which accounts for the repulsive interactions between the particles [8], and r 0 is an arbitrary spatial scale. The total number of bosons in the gas is given by where the norm of the scaled wave function is given by Eq. (5).
It is relevant to note that, as it follows from Eq. (7) and rescaling (10), the above-mentioned critical value, U 0 = 1/4, of the strength of the attractive potential (see Eq. (6)) corresponds to a very small dipole moment, d ∼ 10 −6 Debye, if central charge Q is taken as the elementary charge, and the mass of the particle is ∼ 100 proton masses. Therefore, the overcritical case of U 0 > 1/4, the consideration of which is the main objective of the present article, is relevant in the actual physical context. Taken as Eq. (8), the GPE neglects dipole-dipole interactions between the particles. These interactions can be taken into account in the framework of another application of the mean-field approach. Indeed, the local density of the dipole moment in the gas (i.e., the dielectric polarization of the medium) is P = d |ψ(r)| 2 , hence the electrostatic field generated by the polarization, E d , is determined by the Poisson equation, ∇ · (E d + 4πP) = 0, which can be solved immediately: Then, the extra term in the GPE, accounting for the interaction of the local dipole with the collective field (12), created by all the other dipoles, is This term, if added to Eq. (8), may be absorbed into a redefinition of the scattering length accounting for the repulsion between the particles. In the underlying physical units, this amounts to where m is the mass of the dipolar molecule. For the typical value of a s ∼ 10 nm and the above-mentioned mass of the particle, ∼ 100 proton masses, Eq. (14) demonstrates that the additional term is essential for dipole moments d 0.3 Debye. The rest of the article is organized as follows. In Section II, results are reported for the basic model, outlined above, as per Ref. [5]. Particular subsections of Section II first recapitulate the description of the 3D and 2D collapse in the framework of Schrödinger equation (3), which includes the trapping potential (4), and then present main results obtained in the 3D case on the basis of Eq. (8) (with Ω = 0, as the trapping potential is not a necessary ingredient of the nonlinear model, on the contrary to the linear one). The results explicitly demonstrate the creation of the originally missing GS by the self-repulsive cubic nonlinearity at U 0 > 1/4. In addition, a subsection of Section II reports a new result, viz., a quantum phase transition in the GS of the model which includes the Lee-Huang-Yang (LHY) correction [22] to the mean-field GPE. The summary of results for the 2D nonlinear model are also presented in Section II. It is demonstrated that the cubic self-repulsive term is insufficient for the suppression of the 2D quantum collapse and restoration of the missing GS. This is possible if a quintic repulsive term is included in the GPE, which may account for three-body collisions, or if the quartic LHY correction is included in the effective two-dimensional GPE . A short subsection concluding Section II formulates a challenging problem of the consideration of the quantum collapse in the gas of fermions.
Section III addresses, along the lines of Ref. [6], the collapse suppression and creation of the GS in the 3D model with the symmetry of the effective attractive potential reduced from spherical to cylindrical by an external field which polarizes dipole moments of the particles. In this version of the model, states carrying the angular momentum are constructed too, in addition to the GS.
Section IV deals with a two-component model in 3D, which makes it possible to consider the interplay of the collapse suppression and the transition between miscibility and immiscibility in the binary system. A weak quantum phase transition, which occurs in that setting, is also briefly considered in Section IV.
Section V presents results for the basic 3D model, considered in terms of the many-body quantum theory, as per Ref. [51], with the help of variational approximation for the many-body wave functions and numerically implemented Monte Carlo method. The main result is that, strictly speaking, the quantum collapse is not fully suppressed in the many-body theory, but, nevertheless, the noncollapsing self-trapped state, predicted by the mean-field theory, exists as a metastable one, insulated from the collapse by a tall potential barrier.
The paper is concluded by Section VI, which also suggests directions for further work on this general topic.

II. THE BASIC THREE-AND TWO-DIMENSIONAL MODELS
This section summarizes results produced in Ref. [5]. The quantum phase transition driven by the LHY correction to the mean-field theory, briefly outlined in subsection 3, is a new finding.

A. The quantum collapse in the linear Schrödinger equation
First, it is relevant to recapitulate the analysis of linear Schrödinger equation (3), to which the trapping potential (4) is added: Stationary solutions of Eq. (15) in 3D spherical coordinates, (r, θ, ϕ), are looked for as where µ is the energy eigenvalue (or chemical potential, in terms of the GPE), Y lm (θ, ϕ) is the spherical harmonic with quantum numbers (l, m), and radial wave function Φ(r) is real. Substituting ansatz (16) in Eq. (15), two exact solutions for Φ(r) can be found: which exist under condition The smaller value of µ (in the case of l = 0, it defines the GS of the system under the consideration) corresponds to σ + , i.e.,, the top sign in Eq. (18). The wave function is characterized by its norm (5), where Γ is the Gamma-function. Equation (20) shows why the trapping potential ∼ Ω 2 is necessary for the existence of physically relevant (normalizable) eigenmodes of the linear Schrödinger equation (15), as norm (20) diverges (due to its weak localization at r → ∞) in the limit of Ω → 0. These solutions for the stationary wave functions do not exist at U l > 1/4 (note that the presence of the angular momentum, l ≥ 1, secures the existence of the bound states at essentially larger values of U 0 , as per Eq. (19)). The nonexistence of stationary wave functions implies that the system suffers the onset of the quantum collapse, as confirmed by simulations of time-dependent equation (15), see an example in Fig. 1. Indeed, a set of instantaneous profiles of √ r|ψ(r, t)|, shown in Fig. 1 for the weakly overcritical case, U 0 = 0.27, with l = 0, confirm the development of the self-compression (finally, collapse) of the wave function towards r = 0 (in the simulations, the collapse is eventually arrested due to a finite mesh size of the numerical scheme).
In 2D, the GS solution to Eq. (15) exists only for U 0 < 0. In the exact form, the GS wave function is given by Eqs. (16) and (17), but with (18) replaced by Direct simulations of the 2D equation (15) at U 0 > 0 also demonstrate the onset of the collapse dynamics.
B. The three-dimensional ground state (GS) created by the cubic self-repulsive nonlinearity The most essential results may be produced by GPE (8) without an external trapping potential, hence the equation simplifies to The substitution of ψ = e −iµt Φ(r) for isotropic stationary states of Eq. (22) (here, solely l = 0 is considered, cf. Eq. (16), with the intention to construct the GS, which always has l = 0) yields equation Scaling invariance of Eq. (23) at r → 0 suggests that the respective asymptotic form of the solution should be Φ ∼ 1/r, therefore solutions are looked for as with function χ(r) obeying equation Asymptotic forms of solutions to Eq. (25) can be readily constructed for r → 0 and r → ∞. First, the expansion at r → 0 yields where χ 1 is a free constant, in terms of this expansion. At r → ∞ the asymptotic form of the bound-state solution with µ < 0 is where χ 0 is an arbitrary constant, in terms of the expansion for r → ∞. A global analytical approximation can be constructed as an interpolation, stitching together the asymptotic forms (26) (where the correction term ∼ χ 1 is neglected, in the present approximation) and (27): Note that the singularity of wave function (28) at r → 0 is acceptable, as the respective integral (5) converges at small r. It is also relevant to mention that, following the substitution of the asymptotic form (28) in the effective pseudopotential in Eq. (22), which includes the nonlinear term, U pseudo (r) ≡ −(1/2)U 0 r −2 + |ψ(r)| 2 , the singularity ∼ r −2 at r → 0 cancels out in it. Note also that a more singular attractive potential, U (r) = −U 0 /r b , with U 0 > 0 and b > 2, gives rise to asymptotic form |ψ| 2 ≈ U 0 /r b of the solution at r → 0, hence the corresponding norm still converges at b < 3.
Due to the nonlinearity of Eq. (22), the chemical potential of the GS depends on its norm. Using approximation (28), it is easy to calculate µ as a function of N : In fact, scaling µ ∼ N −2 is an exact property of solutions to Eq. (22), which follows from a straightforward analysis of this equation. Note also that, in the limit of µ → −0, Eq. (28) gives a particular exact solution of Eq. (22), although its norm diverges. Equation (25) can be easily solved in a numerical form. A typical example of the numerical GS solution, along with approximation (28), is displayed in Fig. 2(a) for U 0 = 0.8, which is essentially larger than the critical value of the attraction strength, (U 0 ) (3D) cr = 1/4 (see Eq. (6)), beyond which linear Schrödinger equation (15) has no GS. Further, Figs. 2(b) and 2(c) represent the family of the GS states, by means of dependences µ(N ), for two values, U 0 = 0.8 and 0.1, which are, respectively, larger and smaller than 1/4. Thus, on the contrary to the linear Schrödinger equation, GPE (22) maintains the GS at all values of U 0 and N . In other words, the inclusion of the repulsive cubic term in Eq. (22) completely suppresses the quantum collapse in the 3D space, creating the GS where it does not exist in the linear Schrödinger equation.
The analytical approximation (28) suggests an estimate for the radial size of the GS created by the repulsive nonlinearity: It is relevant to rewrite this estimate in terms of physical units, as per Eqs. (10), (11), and (14): Note that arbitrary spatial scale r 0 , which was used in rescalings (10) and (11), cancels out in Eq. (11). Thus, GPE (22) uniquely predicts the radius of the restored GS, in terms of physical parameters of the model. It is natural that r GS , given by Eq. (32), shrinks to zero in the limit of vanishing nonlinearity, which is tantamount to N ph → 0, implying the onset of the collapse in the framework of the linear Schrödinger equation. Note also too  (6)) for the linear Schrödinger equation (15). Here, solid and dashed curves depict, respectively, the numerical results and analytical approximation given by Eqs. (28) and (29) (in panels (b) and (c), the curves follow scaling µ ∼ N −2 , which is an exact property of Eq. (22)). In particular, the analytical approximation predicts N (µ = −0.225) = 5.30 for U0 = 0.8 (the solution shown in (a)), while the numerically found counterpart of this value is Nnum(µ = −0.225) = 6.26. The convergence of the numerical and analytical curves for N (µ) at µ → 0 corresponds to the fact that Eq. (28) gives exact solution (17) in this limit.
that, if the contribution from by the dipole-dipole interactions (∼ d 2 ) dominates over the contact interactions in Eq. (32) (md 2 2 a s ), the latter result strongly simplifies, taking into regard Eq. (7): (r GS ) ph = √ 2d/|Q| N ph . Then, for Q equal to the elementary charge, d ∼ 1 Debye, and N ph ∼ 10 5 , the latter estimate predicts the GS with radius ∼ 3 µm. This result upholds the self-consistency of the model, as the mean-filed approximation (and the respective GPE) are definitely applicable for scales 1 µm [8].
It is worthy to stress that Eq. (8), which does not include the trapping potential, predicts the GS with the finite norm at U 0 < 1/4, as the norm of the corresponding stationary solutions to the linear equation (15), see Eqs. (16) and (17), diverges at Ω = 0. Lastly, simulations of Eq. (8) with random perturbations added to the stationary solutions demonstrate that the GS is always dynamically stable [5]. The stability agrees too with the anti-Vakhitov-Kolokolov criterion, dµ/dN > 0, which is a necessary stability condition for localized states supported by self-repulsive nonlinearities [15] (the original Vakhitov-Kolokolov criterion, dµ/dN < 0, is the necessary stability condition in the case of self-attraction [16,17]) .
C. The quantum phase transition induced by the Lee-Huang-Yang (LHY) correction to the mean-field theory The singularity ∼ r −1 of the stationary wave function at r → 0, as seen in Eqs. (24) and (26), suggests that, although the singularity is integrable, as the respective 3D integral for the total norm converges, the LHY correction [22] to the mean-field theory, which is relevant for higher values of the condensate's density, must be taken into regard. As shown in Ref. [23,24], the scaled GPE with this correction, represented by coefficient g LHY > 0, is Then, the asymptotic form of the stationary wave function at r → 0, which was found above in the form determined by the cubic term in Eq.
under the condition of U 0 > 2/9, while at 0 < U 0 < 2/9 the asymptotic form is determined by the linearization of Eq. (33), leading to the same result as given by Eqs. (17) and (18) with σ = σ − (and U l replaced by U 0 ): where Φ 0 is an arbitrary constant in terms of the expansion at r → 0. Note that the wave function with asymptotic form r −σ+ , which corresponds to the GS in the linear Schrödinger equation, is incompatible with the presence of the LHY term in Eq. (33), although power −2/3 in expression (33) coincides, at the critical point, U 0 = 2/9, with σ + , rather than σ − . Thus, the jump from the asymptotic form (35), produced by the linear Schrödinger equation, to one (34) generated by the LHY term at U 0 = 2/9 (in particular, the jump between σ − and σ + ), takes place at U 0 = 2/9, which is a signature of a quantum phase transition. Examples of such phase transitions were studied in many-body settings [25] and in many other systems [26]- [30].
Lastly, the LHY term may stabilize the bosonic gas pulled to the center by potential (1) even in the case of the effective attractive interaction, which is possible in a binary condensate, with intrinsic self-repulsion in each components, and dominating attraction between them, as proposed in Refs. [23] and [24], and realized experimentally in the form of "quantum droplets" (in the binary condensate of 39 K) in Refs. [31]- [33]. For the symmetric mixture, with equal wave functions of the two components, the effective GPE takes the form of Eq. (33) with the opposite sign in front of the cubic term. This model may be a subject for special consideration.
D. The two-dimensional ground state created by the quintic self-repulsive nonlinearity As mentioned above, the GPE in the form of Eq. (8) may be relevant, as a physical model, in 2D too. However, the 2D version of norm (5) of the wave function with asymptotic form ∼ r −1 at r → 0, which follows from this equation (see Eq. (28)), logarithmically diverges at small r. This means that, the cubic self-repulsion is not strong enough to suppress the collapse in the 2D geometry. On the other hand, the GPE may also include the quintic repulsive term accounting for three-body collisions, provided that the collisions do not give rise to conspicuous losses [34,35].
The 2D GPE can be derived from the underlying 3D version if tight confinement, with the respective harmonicoscillator length, a ⊥ , is imposed in direction z by the trapping harmonic-oscillator potential, reducing the effective dimension to that of the plane with remaining coordinates (x, y) [36]- [38]. In particular, if the dominating quintic terms appears in the 3D GPE with coefficient g 5 , the reduction to the 2D equation replaces it by In the scaled form, the 2D equation it written in the polar coordinates, (r, θ), as Stationary solutions to Eq. (36) (not only the GS, but also for states carrying the angular momentum) are looked for as where integer l is the azimuthal quantum number, cf. Eq. (24). The substitution of this ansatz in Eq. (36) yields an equation for real χ 2D (r): with U (2D) l ≡ U 0 − l 2 , cf. Eq. (19). Note that, unlike the 3D case, in 2D nonlinear model the analysis is possible equally well for l = 0 (the GS) and l ≥ 1.
The expansion of the solution to Eq. (38) at r → 0 yields where s = (1/2) 1 + 5 + 16U > 0, the suppression of the collapse and creation of the GS (l = 0), or the state with l ≥ 1, by the quintic self-repulsive term.
Combining the 2D asymptotic form (39), valid at r → 0, and the obvious approximation valid at r → ∞, χ 2D ≈ χ 0 exp − √ −2µr , one derives an interpolation formula for the GS, and the dependence µ(N ) following from it, cf. Eqs. (28) and (29) in the 3D case: Similar to the situation in the 3D case, Eq. (40) gives an exact wave function with a divergent norm in the limit of µ → −0, cf. Eq. (30). The approximation (40) makes it possible to define the rms radial size of the two-dimensional GS, cf. Eq. (31) in the 3D case: Note that the quintic term supports the GS in 2D even at 0 < −U (2D) l < 1/4, when the central potential is repulsive. The correctness of this counter-intuitive conclusion is corroborated by the above-mentioned fact that the analytical approximation (40) gives exact solution (41)  Generally, the results for the 2D model are more formal than those summarized above for 3D, as the realization of the dominant quintic nonlinearity in BEC is problematic, in experimentally relevant settings. On the other hand, the LHY correction to the GPE is also sufficient to provide the suppression of the quantum collapse and restoration of the GS in the 2D setting. The same dimension-reduction procedure as outlined above, will replace the original LHY coefficient in 3D equation (33) by 2/5π −3/4 a −3/2 ⊥ g LHY . Finally, the quartic LHY term determines the asymptotic form of the wave functions at r → 0 as ∼ r −2/3 , which provides for the convergence of the 2D norm, i.e., it secures the existence of the GS in the 2D model including the LHY term. E. A challenging issue: the Fermi gas pulled to the center An interesting possibility is to elaborate the 3D model for the gas of fermions pulled to the center by potential (1). In a rigorous form, this is a challenging problem, as the dynamical theory for Fermi gases cannot be reduced to a simple mean-field equation [39]. Nevertheless, there is a relatively simple approach to the description of stationary states in a sufficiently dense gas, which is based on a time-independent equation for the real fermionic wave function, Φ (r) [40][41][42], with a nonlinear term of power 7/3 generated by the density-functional approximation, even in the absence of direct interaction between the fermions, which is forbidden by the Pauli principle. In the scaled form, this equation, including potential (1) and chemical potential µ, is The asymptotic form of the solution to Eq. (43) at r → 0 is This result demonstrates a problem similar to the one stressed above in the case of the 2D model with the cubic nonlinearity: the substitution of expression (44) in the 3D integral (5) for the norm of the wave function leads to the logarithmic divergence at r → 0, hence the relatively weak nonlinearity in Eq. (43) is insufficient for the suppression of the 3D quantum collapse of the Fermi gas pulled to the center by potential (1), and a more sophisticated analysis is necessary in this case.

III. THE THREE-DIMENSIONAL MODEL WITH CYLINDRICAL SYMMETRY
The presentation in this section follows the original analysis reported in Ref. [6].

A. Formulation of the model
The previous section addressed the most fundamental spherically symmetric configuration in the 3D space. Because the geometry plays a crucially important role in determining properties of the bound states produced by the model, it is interesting to consider physically relevant settings in 3D with the spatial symmetry reduced from spherical to a lower one. In particular, it is possible to consider the model in which a strong uniform external field is applied to the quantum gas, so that all the dipole moments, carried by the particles, are polarized not towards the center, but in a fixed direction (z), so that d = de z . This configuration gives rise the cylindrically symmetric potential of the interaction of the dipolar particle with the fixed attractive center: where cos θ ≡ z/r. If the polarizing external field is electric, it acts on the central charge too. For this reason, a more relevant situation corresponds to the case when the ultracold gas is composed of Hund A type of molecules, with mutually locked electric and magnetic dipoles. Then, external uniform magnetic field may be employed to align the dipoles in the fixed direction [44].
The 3D GPE with potential (45) is Along with the consideration of the GS, it is also relevant to construct eigenmodes carrying the orbital angular momentum, which corresponds to the azimuthal quantum number, m: where the spherical coordinates are used again, cf. Eq. (16), and real eigenmode φ should be found as a solution of the equation following from the substitution of ansatz (47) in Eq. (46):

B. The linear Schrödinger equation with the cylindrical symmetry
The analysis of the present model should start with identifying conditions for the existence of the GS in the respective linear Schrödinger equation, obtained by dropping the cubic term in Eq. (48). At r → 0, an asymptotic solution to the linear equation is looked for as Φ (r, θ) = r −σ χ lin (θ). The substitution of ansatz (49) in the linearized version of Eq. (48) and dropping term µΦ, which is negligible for the asymptotic analysis at r → 0, leads to an equation that can be written in terms of ξ ≡ cos θ: For U 0 = 0, equation (50) with integer values may be solved in terms of the associated Legendre functions, l ≥ m being the orbital quantum number. Note that the singular wave function (49) is 3D-normalizable, at r → 0, for σ < 3/2, i.e., as it follows from Eq. (51), solely for the GS, with m = l = 0 and σ = 1. The onset of the quantum collapse is signalled by a transition in Eq. (50) from real eigenvalues σ to complex ones. Because the effective eigenvalue in the equation is ǫ ≡ σ 2 − σ, i.e., σ = (1 + √ 1 + 4ǫ)/2, the transition to complex σ happens at point ǫ = −1/4 (i.e., σ = 1/2). At U 0 = 0, Eq. (50) cannot be solved in terms of standard special functions. A result of a numerical solution is displayed in Fig. 4. It demonstrates that, for given azimuthal quantum number m, eigenvalue σ decreases, with the increase of U 0 from zero to some critical value (U 0 ) cr , from σ (U 0 = 0) = m + 1 to σ (U 0 = (U 0 ) cr ) = 1/2. For lowest values of m, the numerically found critical values of the potential strength, at which σ = 1/2 is attained are (U 0 ) cr (m = 0, 1, 2) = 1.28, 7.58, 19.06.
In the framework of the linear Schrödinger equation, the quantum collapse takes place at U 0 > (U 0 ) cr (m). It is relevant to compare critical values (52) of the strength of the axisymmetric potential with those given by Eq. (19) for the spherically isotropic one: The comparison naturally shows that the critical strengths are much lower for the spherical potential, which provides stronger pull to the center.

C. Suppression of the quantum collapse by the repulsive nonlinearity under the cylindrical symmetry
As well as in the isotropic setting, cf. Eq. (24), the repulsive cubic term in Eq. (48) may balance the attractive potential ∼ −r −2 if, at r → 0, the wave function contains the singular factor r −1 (rather than generic r −σ in Eq. (49)). Then, the substitution of Φ(r, θ) = r −1 χ(r, θ), transforms Eq. (48) into the following equation: (recall ξ ≡ cos θ). Note that Eq. (54) makes it possible to write the norm of the 3D wave function as To analyze solutions to Eq. (55) at r → 0, one may expand them as assuming s > 0, which leads to the following equation for χ 0 (ξ), that does not admit an exact solution: cf. Eq. (50). Bound states produced by Eq. (55) were found by means of numerical methods in Ref. [6]. Typical profiles of solutions for function χ (r, ξ), generated by Eq. (55), are displayed in Fig. 5, for m = 0, 1, 2 and fixed norm, N = 2π.
A crude analytical form of the solutions is provided by the Thomas-Fermi approximation (TFA), which neglects all derivatives in Eq. (55) [8]: Actually, this approximation for m ≥ 1 exists only for U 0 > 3 √ 3/2 m 2 (otherwise, Eq. (59) yields χ 2 TF ≡ 0). Families of the bound states with different quantum numbers m are presented in Fig. 6 by a set of curves showing the chemical potential, µ, versus nonlinearity strength, U 0 , for a fixed norm (N = N 0 ≡ 2π; producing the results for a fixed norm is sufficient, as the scaling invariance of Eq. (55) implies an exact property, µ (U 0 , N ) = (N/N 0 ) −2 µ (U 0 , N 0 ), the same as mentioned above for the isotropic configuration).
As seen in Fig. 6(a), this simple approximation is reasonably close to its numerically found counterpart. Lastly, the stability of the bound states against perturbations was verified in Ref. [6] by direct simulations of the underlying GPE (46), demonstrating complete stability that the families of the states for m = 0, 1, and 2.  (45), with the the reduced (cylindrical) symmetry, and for the fixed norm, N = 2π, as obtained in Ref. [6]. The dashed curve in (a) additionally shows dependence (60) predicted by the TF approximation.

IV. THE TWO-COMPONENT SYSTEM IN THREE DIMENSIONS: THE SUPPRESSION OF THE QUANTUM COLLAPSE IN MISCIBLE AND IMMISCIBLE SETTINGS
This section summarizes results of the analysis reported in Ref. [7].

A. The formulation of the model and analytical considerations
The generalization of basic model (22) for a binary bosonic gas, with component wave functions ψ 1 and ψ 2 , is provided by the system of nonlinearily coupled GPEs: (for the consistency with Ref. ([6]), parameter V 0 ≡ U 0 /2 is now used as the strength of the potential pulling particles to the center), where γ is the relative strength of the inter-component repulsion, and coefficients of the self-repulsion are scaled to be 1. Spherically symmetric bound states with chemical potentials µ n < 0, n = 1, 2, of the two components are looked for as with real radial functions χ n (r) obeying the coupled equations, cf. Eqs. (24) and (25). In terms of these functions, the norms of the components are and the rms radial size of the trapped mode in each component is defined as cf. Eq. (31). An expansion of solutions to Eqs. (63) at r → 0 is looked for as with s > 0, cf. Eq. (26) (here, c 1 = c 2 is possible, but power s must be the same for χ 1 and χ 2 ), which leads to a system of algebraic equations for leading-order coefficients χ (0) n : Equations (67) give rise to solutions of two types, corresponding to mixed and demixed states in the binary gas: The numerical analysis performed in Ref. [6] has demonstrates that demixed modes do not exist at γ < 1, when the mutual repulsion is weaker than the self-repulsive nonlinearity, while mixed ones are completely unstable in the opposite case, γ > 1. Thus, unlike other systems featuring the miscibility-immiscibility transitions [45]- [47], in the present situation the transition point is not shifted, under the action of the confining potential, from the commonly known free-space point, γ = 1 [48]. Further analysis demonstrates a change in the structure of the r-dependent corrections in Eq. (66) for the miscible system, with γ < 1: at V 0 < 1/2, the dominant terms are ∼ r (1+ √ 1+16V0)/2 , while at V 0 > 1/2 these are terms ∼ r 2 . This break of analyticity, which happens, with the increase of V 0 , at V 0 = 1/2, implies that a weak quantum phase transition happens at this value of V 0 , although the well-defined GS exists equally well at V 0 < 1/2 and V 0 > 1/2. Precisely at V 0 = 1/2 ≡ (V 0 ) phase−trans , expansion (66) is replaced by Note that the present phase transition is weak in comparison with the above-mentioned one, driven by the LHY correction to the mean-field theory, which gives rise to the jump between the different asymptotic forms of the wave function, given by Eqs. (34) and (35). For the comparison with the present setting, based on the binary BEC, especially relevant are previously investigated phase transitions in binary fluids [49]. Recall that the onset of the quantum collapse in the linear version of the model occurs at the critical value V (cr) which is a quarter of the value of the potential's strength at the phase-transition point, (V 0 ) phase−trans = 1/2. Lastly, at r → ∞ Eqs. (63) yield an exponential asymptotic form of the solution, where constants χ   7: (a) The numerically found profile of wave functions χ1(r) = χ2(r), at V0 = 1, γ = 0.9, and N1 = N2 = 4π, of the GS in the miscible binary system, as found in Ref. [7], and its comparison with the analytical approximation given by Eq.
Comparison of expression (73) with numerical results is shown in Fig. 7(b) by the dashed line. This approximation is accurate for sufficiently small V 0 , but becomes inaccurate for large V 0 . For larger V 0 , TFA can be applied to the mixed balanced mixture, with N 1 = N 2 ≡ N , which yields (for cf. TFA for the potential with the cylindrical symmetry, given by Eq. (59). The substitution of approximation (75) in Eqs. (5) and (65) yields the predictions for the chemical potential and effective size of the GS: (recall R 0 is the TFA cutoff radius defined in Eq. (75)). Analytical approximations (72) and (75) (shown by the short-and long-dashed lines, respectively) are compared to the numerically found profile of the GS in Fig. 7(b). A general conclusion (see details in Ref. [7]) is that, quite naturally, TFA works better for larger V 0 , while interpolation (72) is more accurate for smaller V 0 . Numerically generated profiles of imbalanced mixed GSs are displayed in Fig. 3(a) at V 0 = 2 and γ = 0.9 for N 1 = 4π and N 2 = 2π. The imbalanced mixed states with µ 1 = µ 2 and N 1 = N 2 feature equal values of χ 1,2 (r = 0), in agreement with Eq. (68).
In the case of the strong pull to the center, V 0 ≫ 1, TFA can be generalized for imbalanced states, fixing |µ 1 | ≤ |µ 2 | for the definiteness' sake. Then, TFA is constructed in a two-layer form, technically similar to that applied to the so-called symbiotic gap solitons in Ref. [50]. In the inner layer, both wave functions are different from zero: In the outer layer, only one component is present, in the framework of TFA: χ 2 ≡ 0, Both components of the TFA solution, given by Eqs. (78)-(80), are continuous at r = r 0 and r = R 0 . The two-layer TFA for a typical imbalanced GS is compared to its numerical counterpart in Figs. 8(b,c). The analysis reported in Ref. [7] also includes the consideration of a two-component system with attraction between the components, in the case when only one component is subject to the action of the pull-to-the-center potential, while the other one plays the role of a buffer. In particular, the interpolation, similar to that based on Eq. (72), produces a sufficiently accurate prediction in that case.

Immiscible ground states
As said above, in the case of γ > 1 relevant states are immiscible ones. The two-layer TFA may be applied to produce an immiscible GS. In the inner layer, the approximation yields In the outer layer, which is r 2 0 < r 2 < R 2 0 = V 0 /(−µ 2 ), the result is i.e., TFA predicts complete separation between the components in the immiscible state. Figure 9 compares the approximation to numerical results. Of course, the immiscible components are not completely separated in the numerical solution.

V. THE MEAN-FIELD PREDICTIONS VERSUS THE MANY-BODY QUANTUM THEORY
A. Introduction to the section The analysis presented above was performed in Refs. [5]- [7] in the framework of the mean-field theory, i.e., the respective GPEs (possibly including the beyond-mean-field LHY corrections, see Eq. (33)). A relevant issue is comparison of the basic mean-field predictions, such as the suppression of the quantum collapse and creation of the originally missing GS, with the consideration of the many-body system of repulsively interacting quantum bosons, pulled to the center by potential (1), which is taken here as U (r) = −U 0 /r 2 , i.e., U 0 in Eq. (1) is replaced by 2U 0 , to make the notation consistent with that in Ref. [51]. This problem was addressed in Ref. [51]. Results produced in that work are recapitulated in the present section. The many-body Hamiltonian representing the setting under the consideration iŝ where r j are coordinates of the j-th particle in the 3D space, m is the particle's mass, and V int (r) is the potential of the repulsive interaction between the particles. In the framework of the mean-field theory, V int (r) is characterized solely by the s-wave scattering length [8], while the many-body system should be introduced with a particular form of the interaction potential. Two basic forms of the interaction potential chosen for the analysis are specified below, see Eqs. (88) and (89). Before introducing the many-body wave function, the single-particle one is adopted as per the following ansatz: where α ≥ 0 determines the inverse localization length, which affects the system's size and, consequently, the density. Alternatively, α can be interpreted in terms of an effective external harmonic confinement with frequency Ω = 2α /m, cf. Eq. (4). At r j → 0, the shape of the wave function is controlled by parameter β in ansatz (84).
B. The single-particle solution The single-particle problem, defined by Hamiltonian (83) with N = 1, can be studied by means of the variational method, treating α and β in ansatz (84) as variational parameters. In the single-particle sector, the system is steered by the competition of the external potential and kinetic energy, while the interparticle potential, V int (|r i − r j |), does not appear. The variational energy, E For a fixed localization size, α = const, this energy is a decreasing function of β if U 0 is smaller than the critical value for the onset of the collapse, U 0 = 1/8 -the same which appears in Eq. (6). On the other hand, a metastable state may appear in the many-body system with repulsive interparticle interactions. Actually, it corresponds to the mean-field GS predicted by the solution of the GPE in Ref. [5], see further details below In the framework of the local-density approximation, the chemical potential of the state with uniform density n is take as µ hom = gn , where g = 4π 2 a s /m is the coupling constant. This choice corresponds to the short-range interaction potential determined by the s-wave scattering length as per the Born approximation. Further, the chemical potential in the presence of the external field is approximated by the sum of the local chemical potential µ loc = gn, where, this time, n is a function of the coordinates, rather than a constant, and the external potential, where the harmonic-oscillator confinement, with the respective length scale, a ho = /(mΩ), is added to make the size of the system finite, cf. potential (4) used above. Solving Eq. (86) for the density, one obtains the following density profile: where the radius of the gaseous cloud is taken as per TFA, R TFA = µ + µ 2 + 2mU 0 Ω 2 /( √ mΩ). The density at the center features an integrable divergence in Eq. (87), reflecting the presence of the attractive central potential, cf. Eq. (24). Finally, the chemical potential itself is fixed by the normalization condition, 4π RTFA 0 n(r)r 2 dr = N . To study the expected scenarios of the system's evolution, two different potentials of the inter-particle interaction were introduced in Ref. [51], viz., the hard-sphere potential of diameter R, and its soft-sphere counterpart, with finite V 0 in the latter case. By varying height V 0 of the soft-sphere potential, one can alter the respective s-wave scattering length, which is where the momentum corresponding to the height of the soft-sphere potential is For hard-sphere potential (88), the effective s-wave scattering length is identical to the diameter of the sphere, a s = R.

C. The Monte-Carlo method
An efficient way to calculate the energy of a many-body system is to use the Monte-Carlo technique. In Ref. [51], the variational Monte-Carlo method was employed, which samples the probability distribution, p = |ψ| 2 , for a known many-body wave function, ψ, allowing one to calculate the variational energy as a function of trial parameters, such as α and β in Eq. (84). The well-known Metropolis algorithm [52] was used for the implementation of the method.
The many-body trial wave function was chosen as a product of single-particle terms, f 1 (r), taken as per Eq. (84), and a pairwise product of two-particle Jastrow terms [53], f 2 (r): The Jastrow factor f 2 (r) in Eq. (92) is chosen as a solution of the linear Schrödinger equation for two-body scattering. In this way, interparticle correlations, which are important in the context of the metastability of the many-body system, are retained in the analysis. For the hard-sphere potential, the two-body solution is given by while for the soft-sphere potential (89), it is where k is given by Eq. (91), and constant A is determined by the condition of the continuity of f 2 (r) at r = R.
D. Numerical results for the many-body system Figure 10 shows the variational energy, calculated by the Monte Carlo method for a fixed radius of the soft sphere, R = 1.3a s , and a wide range of values of the number of particles, from N = 2 up to N = 10000. For small values of α, which corresponds to the weak localization, the energy may be negative (this is not visible in the log-log plot of Fig. 10). As the localization gets tighter, the energy becomes positive, as the two-body interaction helps the system to resist the trend to collapsing. For very tight localization, α → ∞, the collapse is observed for small values of N , with the energy diverging towards −∞.
For large N , the energy calculated with ansatz (84) does not immediately lead to the fully collapsed state. The localization energy, proportional to the energy scale, Ω, of the trapping potential, may become a dominating term in the energy, while the system's size is still large enough, so that the fully-collapsed state is not realized. The energy in the corresponding regime is numerically approximated as with C = 4. It is shown in Fig. 10 by the dashed-dotted line. For still tighter localization, it has been found that, in the limit of the full collapse, when all particles overlap, the energy is well approximated by formula The energy of the interparticle interactions, revealed by the calculations, is E int = V 0 N (N − 1)/2. The asymptotic energy (96) is shown in Fig. 3

by dashed lines.
A clear conclusion is that the increase of the number of particles indeed causes strong rise of the potential barrier which stabilizes the metastable energy minimum corresponding to the gaseous state. This can be also concluded from Eq. (96), where the contribution due to the repulsive interactions scales as N 2 for large N , while the term corresponding the attractive central potential scales as N .
The energy barrier between the state described by ansatz (84) and the free state with zero energy, E barrier , is estimated as the maximum value of the energy per particle (see Fig. 3), and is shown in Fig. 4. For a large system's size, the barrier can be approximated by comparing the two basic energy scales given by Eqs. (95) and (96). The resulting asymptotic approximation for the barrier's height is FIG. 11: The energy barrier between the state with α = 0 and α → ∞ in the many-body system, as a function of the number of particles, N , for the data shown in Fig. 10, as per Ref. [51]. The dashed line depicts the asymptotic approximation (97) for the large system.
which is shown in Fig. 4 by the dashed line.

VI. DISCUSSION AND CONCLUSION
This article aims to produce a brief review of results reported in works [5]- [7] and [51], that offer a solution to the known problem of the quantum collapse, alias "fall onto the center" [1], in nonrelativistic quantum mechanics. The quantum collapse occurs in the three-dimensional Schrödinger equation with 3D isotropic attractive potential −U 0 /(2r 2 ). This equation does not have a GS (ground state) if the attraction strength, U 0 , exceeds a final critical value. In that case, the Schrödinger equation gives rise to a nonstationary wave function which collapses, shrinking towards the center. The solution of the collapse problem was proposed in the above-mentioned original works in terms of the gas of bosons pulled to the center by the same potential, with repulsive contact interactions between the particles. The intrinsic repulsion is represented by the cubic term in the respective GPE (Gross-Pitaevskii equation). The setting may be realized as the 3D gas of polar molecules carrying a permanent electric dipole moment and pulled to a central electric charge. The analysis, performed in the framework of the mean-field theory, predicts suppression of the collapse in the gas, and the creation of the missing GS.
An original result, added in this article to the review of the previously published findings, is the quantum phase transition occurring in the 3D model which includes the beyond-mean-field LHY (Lee-Huang-Yang) correction in the GPE, in the form of the self-repulsive quartic term. The phase transition manifests itself by a jump of the asymptotic structure of the wave function (for r → 0) at the critical value of the strength of the attractive potential.
In the 2D version of the model, the cubic self-repulsion is not sufficient to suppress the quantum collapse. In this case, it can be suppressed if a quintic self-repulsive term, representing three-body collisions (provided that they so not give rise to losses) can be added to the underlying GPE. On the other hand, the LHY quartic term, added to the 2D GPE, is sufficient to suppress the quantum collapse and restore the respective GS.
Polarization of dipole moments in the 3D gas by an external uniform field reduces the symmetry of the central attractive potential from spherical to cylindrical. This modification of the system predicts both the GS and stabilized states carrying the angular momentum. A binary condensate, modelled by the system of nonlinearly-coupled GPEs, is considered too, making it possible to study the interplay of the suppression of the collapse in the 3D space and the miscibility-immiscibility transition in the binary BEC.
In addition to the systematic numerical analysis of these mean-field settings, the original works have produced many results by means of analytical approximations, such as combined asymptotic expansions and TFA (Thomas-Fermi approximation). All the states predicted by the mean-field theory in these settings are shown to be completely stable as solutions to the respective time-depending GPEs.
In work [51], the consideration of the same 3D setting was performed in terms of the many-body quantum theory, by means of the variational approximation for the many-body wave function, numerically handled with the help of the Monte-Carlo method. The analysis has demonstrated that, although the quantum collapse cannot be fully suppressed in terms of the many-body theory, the self-trapped states predicted by the mean-field model exist in the full many-body setting too, as metastable ones, protected against the onset of the collapse by a tall potential barrier, whose height steeply grows with the increase of the number of particles in the gas.
As an extension of the work on the topic of this article, it may be interesting to construct modes carrying the angular momentum in the isotropic 3D model, and also to consider the model with a set of two mutually symmetric attractive centers. In particular, it may be relevant to explore a possibility of the spontaneous symmetry breaking of the GS in the latter case.
As mentioned above (see Eqs. (43) and (44), a challenging issue is to develop a consistent analysis for the gas of fermions pulled to the center by the same potential, −U 0 / 2r 2 .