Suppression of Quantum-Mechanical Collapse in Bosonic Gases with Intrinsic Repulsion: A Brief Review

Abstract

It is known that attractive potential ~$-1/r^2$ gives rise to the critical quantum collapse in the framework of the three-dimensional (3D) linear Schrödinger 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.

1. Introduction

One of standard exercises given to students taking a course in quantum mechanics is solving the three-dimensional (3D) Schrödinger equation with an isotropic attractive potential [1],

$$U\left(r\right)=-\frac{U_0}{2r^2},U_0>0.$$

(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, $\mathbf{r}=\left\{x,y,z\right\}$,

$$\frac{d^2\mathbf{r}}{d{t}^2}=-\frac{\partial U}{\partial \mathbf{r}}\equiv U_0\frac{\mathbf{r}}{r^4},$$

(2)

admits obvious rescaling, $t\equiv \tilde{t}/\sqrt{U_0}$, which eliminates $U_0$ from Equation (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 $\psi \left(\mathbf{r},t\right)$,

$$i{\psi }_t=-\frac{1}{2}{\nabla }^2\psi -\frac{U_0}{2r^2}\psi ,$$

(3)

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,

$$U_{\mathrm{trap}}=\frac{1}{2}{\Omega }^2{r}^2,$$

(4)

is added to Equation (3), to make the integral norm

$$N=\int {\left|\psi \left(\mathbf{r}\right)\right|}^{2}d\mathbf{r}$$

(5)

convergent at $r\to \infty$, the Schrödinger equation gives rise to the normal set of trapped modes, starting from the ground state (GS), at

$$U_0<{\left(U_0\right)}_{\mathrm{cr}}^{\left(3\mathrm{D}\right)}=1/4.$$

(6)

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\to -\infty$, 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 ${\left(U_0\right)}_{\mathrm{cr}}^{\left(2\mathrm{D}\right)}=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\gtrless 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, based on the renormalization-group technique, the field-theory formulation 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,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]; that is, 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 $\mathbf{E}=Q\mathbf{r}/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., $\mathbf{d}/d=\mathrm{sgn}\left(Q\right)\left(\mathbf{r}/r\right)$), in order to minimize the interaction energy, the respective interaction potential is $U\left(r\right)=-\mathbf{d}·\mathbf{E}$, which is tantamount to potential (1) with strength

$$U_0=2|Q|d.$$

(7)

As for the dipolar molecules which may be used to build the BEC under the consideration, experimental results suggest that they may be, for example, 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 considering 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:

$$i{\psi }_t=-\frac{1}{2}\left({\nabla }^2+\frac{U_0}{r^2}-{\Omega }^2{r}^2\right)\psi +{\left|\psi \right|}^2\psi .$$

(8)

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 Equation (8), gives rise—in the absence of the attractive potential ($U_0=0$)—to the well-known supercritical wave collapse [15]. A relation of this setting to Equation (8) is that the inclusion of the trapping potential ~${\Omega }^2$ gives rise to stable bound states in the form of spherically symmetric bound states and ones with vorticity $m=1$ (cf. Equation (16) below), provided that norm N does not exceed a certain critical value [16,17,18,19].

The energy (Hamiltonian) corresponding to Equation (8) is

$$E=\frac{1}{2}\int \left[{\left|\nabla \psi \right|}^2-\left(\frac{U_0}{r^2}-{\Omega }^2{r}^2\right){\left|\psi \right|}^2+\left|{\psi }^4\right|\right]d\mathbf{r}.$$

(9)

The scaled variables and constants, in terms of which Equation (8) is written, are related to their counterparts measured in physical units:

$$\mathbf{r}=\frac{{\mathbf{r}}_{\mathrm{ph}}}{r_0},t=\frac{\hslash }{m{r}_0^2}{t}_{\mathrm{ph}},\psi =2\sqrt{\pi a_s}{r}_0{\psi }_{\mathrm{ph}},U_0=\frac{m}{{\hslash }^2}{\left(U_0\right)}_{\mathrm{ph}},\Omega =\frac{m{r}_0^2}{\hslash }{\Omega }_{\mathrm{ph}},$$

(10)

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

$$N_{\mathrm{ph}}=\int {\left|{\psi }_{\mathrm{ph}}\left({\mathbf{r}}_{\mathrm{ph}}\right)\right|}^{2}d{\mathbf{r}}_{\mathrm{ph}}\equiv \frac{r_{0}N}{4\pi a_s},$$

(11)

where the norm of the scaled wave function is given by Equation (5).

Note that as it follows from Equation (7) and rescaling (10), the above-mentioned critical value, $U_0=1/4$, of the strength of the attractive potential (see Equation (6)) corresponds to a very small dipole moment, $d\sim {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 Equation (8), the GPE neglects dipole–dipole interactions between the particles. These interactions can be considered 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 $\mathbf{P}=\mathbf{d}{\left|\psi \left(\mathbf{r}\right)\right|}^2$, hence the electrostatic field generated by the polarization, ${\mathbf{E}}_d$, is determined by the Poisson equation, $\nabla ·\left({\mathbf{E}}_d+4\pi \mathbf{P}\right)=0$, which can be solved immediately:

$${\mathbf{E}}_d=-4\pi \mathbf{P}\equiv -4\pi \mathbf{d}{\left|\psi \left(\mathbf{r}\right)\right|}^2.$$

(12)

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

$$-\left(\mathbf{d}·{\mathbf{E}}_d\right)\psi \equiv 4\pi d^2{\left|\psi \right|}^2\psi .$$

(13)

If added to Equation (8), this term 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

$$a_s\to {\left(a_s\right)}_{\mathrm{eff}}\equiv a_s+m{d}^2/{\hslash }^2,$$

(14)

where m is the mass of the dipolar molecule. For the typical value of $a_s\sim 10$ nm and the above-mentioned mass of the particle (∼100 proton masses), Equation (14) demonstrates that the additional term is essential for dipole moments $d\gtrsim 0.3$ Debye.

The rest of the article is organized as follows. In Section 2, results are reported for the basic model outlined above, as per Ref. [5]. Particular subsections of Section 2 first recapitulate the description of the 3D and 2D collapse in the framework of the Schrödinger Equation (3), which includes the trapping potential (4), and then present main results obtained in the 3D case on the basis of Equation (8) (with $\Omega =0$, as the trapping potential is not a necessary ingredient of the nonlinear model, 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 2 reports a new result, viz. a quantum phase transition in the GS of the model which includes the Lee–Huang–Yang (LHY) correction [20] to the mean-field GPE. The summary of results for the 2D nonlinear model are also presented in Section 2. 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 2 formulates a challenging problem of the consideration of the quantum collapse in the gas of fermions.

Along the lines of Ref. [6], Section 3 addresses 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 in addition to the GS.

Section 4 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 4.

Section 5 presents results for the basic 3D model, considered in terms of the many-body quantum theory, as per Ref. [21], 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. Nevertheless, the non-collapsing 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 6, which also suggests directions for further work on this general topic.

2. 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 Section 2.3, is a new finding.

2.1. The Quantum Collapse in the Linear Schrödinger Equation

First, it is relevant to recapitulate the analysis of the linear Schrödinger Equation (3), to which the trapping potential (4) is added:

$$i{\psi }_t=-\frac{1}{2}\left({\nabla }^2+\frac{U_0}{r^2}-{\Omega }^2{r}^2\right)\psi .$$

(15)

Stationary solutions of Equation (15) in 3D spherical coordinates, $\left(r,\theta ,\phi \right)$, are looked for as

$${\psi }_{3\mathrm{D}}=exp(-i\mu t)Y_{lm}\left(\theta ,\phi \right)\Phi \left(r\right),$$

(16)

where $\mu$ is the energy eigenvalue (or chemical potential, in terms of the GPE), $Y_{lm}\left(\theta ,\phi \right)$ is the spherical harmonic with quantum numbers $\left(l,m\right)$, and radial wave function $\Phi \left(r\right)$ is real. Substituting ansatz (16) in Equation (15), two exact solutions for $\Phi \left(r\right)$ can be found:

$$\begin{array}{c}\Phi \left(r\right)={\Phi }_0{r}^{-{\sigma }_{\pm }}exp\left(-\Omega r^2/2\right),\end{array}$$

(17)

$$\begin{array}{c}\mu =\Omega \left(\frac{3}{2}-{\sigma }_{\pm }\right),{\sigma }_{\pm }\equiv \frac{1}{2}\pm \sqrt{\frac{1}{4}-U_l},\end{array}$$

(18)

which exist under condition

$$U_l\equiv U_0-l\left(l+1\right)<1/4.$$

(19)

The smaller value of $\mu$ (in the case of $l=0$, it defines the GS of the system under the consideration) corresponds to ${\sigma }_{+}$; that is, the top sign in Equation (18). The wave function is characterized by its norm (5),

$$N=4\pi {\int }_0^{\infty }{\Phi }^2\left(r\right)r^{2}dr=2\pi {\Phi }_0^2{\Omega }^{-\left(1\mp \sqrt{\frac{1}{4}-U_l}\right)}\Gamma \left(1\mp \sqrt{\frac{1}{4}-U_l}\right),$$

(20)

where $\Gamma$ is the Gamma-function. Equation (20) shows why the trapping potential $\sim {\Omega }^2$ is necessary for the existence of physically relevant (normalizable) eigenmodes of the linear Schrödinger Equation (15), as norm (20) diverges in the limit of $\Omega \to 0$ due to its weak localization at $r\to \infty$.

These solutions for the stationary wave functions do not exist at $U_l>1/4$ (note that the presence of the angular momentum, $l\ge 1$, secures the existence of the bound states at essentially larger values of $U_0$, as per Equation (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 Figure 1. Indeed, a set of instantaneous profiles of $\sqrt{r}|\psi (r,t)|$, shown in Figure 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$. Note that in the simulations, the collapse is eventually arrested due to a finite mesh size of the numerical scheme.

In 2D, the GS solution to Equation (15) exists only for $U_0<0$. In the exact form, the GS wave function is given by Equations (16) and (17), but with (18) replaced by

$$\mu =\Omega \left(1-{\sigma }_{\pm }\right),{\sigma }_{\pm }=\pm \sqrt{-U_0}.$$

(21)

Direct simulations of the 2D Equation (15) at $U_0>0$ also demonstrate the onset of the collapse dynamics.

2.2. 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

$$i{\psi }_t=-\frac{1}{2}\left({\nabla }^2+\frac{U_0}{r^2}\right)\psi +{\left|\psi \right|}^2\psi .$$

(22)

The substitution of $\psi =e^{-i\mu t}\Phi \left(r\right)$ with real $\Phi \left(r\right)$ for isotropic stationary states of Equation (22) (here only $l=0$ is considered, cf. Equation (16), with the intention to construct the GS, which always has $l=0$) yields equation

$$\mu \Phi =-\frac{1}{2}\left(\frac{d^2\Phi }{d{r}^2}+\frac{2}{r}\frac{d\Phi }{dr}+\frac{U_0}{r^2}\right)\Phi +{\Phi }^3.$$

(23)

The scaling invariance of Equation (23) at $r\to 0$ suggests that the respective asymptotic form of the solution should be $\Phi \sim 1/r$. Therefore, solutions are looked for as

$$\Phi \left(r\right)=\frac{\chi \left(r\right)}{r},$$

(24)

with function $\chi \left(r\right)$ obeying equation

$$\mu \chi =-\frac{1}{2}\left[{\chi }^{″}+\left(\frac{U_0}{r^2}-{\Omega }^2{r}^2\right)\chi \right]+\frac{{\chi }^3}{r^2}.$$

(25)

Asymptotic forms of solutions to Equation (25) can be readily constructed for $r\to 0$ and $r\to \infty$. First, the expansion at $r\to 0$ yields

$$\chi \left(r\right)=\sqrt{U_0/2}+{\chi }_1{r}^{s/2},s\equiv 1+\sqrt{1+8U_0},$$

(26)

where ${\chi }_1$ is a free constant, in terms of this expansion. At $r\to \infty$, the asymptotic form of the bound-state solution with $\mu <0$ is

$$\chi ={\chi }_{0}exp\left(-\sqrt{-2\mu }r\right),$$

(27)

where ${\chi }_0$ is an arbitrary constant, in terms of the expansion for $r\to \infty$. A global analytical approximation can be constructed as an interpolation, stitching together the asymptotic forms (26) (where the correction term ~${\chi }_1$ is neglected in the present approximation) and (27):

$$\psi (r,t)=\sqrt{\frac{U_0}{2}}{e}^{-i\mu t}{r}^{-1}{e}^{-\sqrt{-2\mu }r}.$$

(28)

Note that the singularity of wave function (28) at $r\to 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 Equation (22) (which includes the nonlinear term, $U_{\mathrm{pseudo}}\left(r\right)\equiv -(1/2)U_0{r}^{-2}+{\left|\psi \left(r\right)\right|}^2$), the singularity $\sim r^{-2}$ at $r\to 0$ cancels out in it. Note also that a more singular attractive potential, $U\left(r\right)=-U_0/r^b$, with $U_0>0$ and $b>2$, gives rise to asymptotic form ${|\psi |}^2\approx U_0/r^b$ of the solution at $r\to 0$. Hence, the corresponding norm still converges at $b<3$.

Due to the nonlinearity of Equation (22), the chemical potential of the GS depends on its norm. Using approximation (28), it is easy to calculate $\mu$ as a function of N:

$$\mu =-\frac{1}{2}{\left(\frac{\pi U_0}{N}\right)}^2.$$

(29)

In fact, scaling $\mu \sim N^{-2}$ is an exact property of solutions to Equation (22), which follows from a straightforward analysis of this equation. Note also that in the limit of $\mu \to -0$, Equation (28) gives a particular exact solution of Equation (22),

$${\psi }_{\mu =0}\left(r\right)=\sqrt{U_0/2}{r}^{-1},$$

(30)

although its norm diverges at $r\to \infty$b.

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 Figure 2a for $U_0=0.8$, which is essentially larger than the critical value of the attraction strength, ${\left(U_0\right)}_{\mathrm{cr}}^{\left(3\mathrm{D}\right)}=1/4$ (see Equation (6)), beyond which linear Schrödinger Equation (15) has no GS. Further, Figure 2b,c represent the family of the GS states by means of dependences $\mu \left(N\right)$ for two values, $U_0=0.8$ and $0.1$, which are, respectively, larger and smaller than $1/4$. Thus, in contrast 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 Equation (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:

$$r_{\mathrm{GS}}^2\equiv \frac{4\pi }{N}{\int }_0^{\infty }{\left|\psi \left(r\right)\right|}^2{r}^{4}dr=\frac{N^2}{2{\pi }^2{U}_0^2}.$$

(31)

It is relevant to rewrite this estimate in terms of physical units, as per Equations (10), (11), and (14):

$${\left(r_{\mathrm{GS}}\right)}_{\mathrm{ph}}\equiv r_0{r}_{\mathrm{GS}}=\frac{2\sqrt{2}\left({\hslash }^2{a}_s+m{d}^2\right)N_{\mathrm{ph}}}{m{\left(U_0\right)}_{\mathrm{ph}}}.$$

(32)

Note that the arbitrary spatial scale $r_0$ which was used in rescalings (10) and (11) cancels out in Equation (11). Thus, GPE (22) uniquely predicts the radius of the restored GS in terms of the physical parameters of the model.

It is natural that $r_{\mathrm{GS}}$, given by Equation (32), shrinks to zero in the limit of vanishing nonlinearity, which is tantamount to $N_{\mathrm{ph}}\to 0$, implying the onset of the collapse in the framework of the linear Schrödinger equation. Note also that, if the contribution from the dipole–dipole interactions (∼$d^2$) dominates over the contact interactions in Equation (32) ($m{d}^2\gtrsim {\hslash }^2{a}_s$), the latter result strongly simplifies, taking into regard Equation (7): ${\left(r_{\mathrm{GS}}\right)}_{\mathrm{ph}}=\left(\sqrt{2}d/|Q|\right)N_{\mathrm{ph}}$. Then, for Q equal to the elementary charge, $d\sim 1$ Debye, and $N_{\mathrm{ph}}\sim {10}^5$, the latter estimate predicts the GS with radius $\sim 3$ $\mathsf{\mu }$m. This result upholds the self-consistency of the model, as the mean-field approximation (and the respective GPE) are definitely applicable for scales ≳1 $\mathsf{\mu }$m [8].

It is worth stressing that Equation (8), which does not include the trapping potential ($\Omega =0$), 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 Equations (16) and (17)) diverges at $\Omega =0$. Lastly, simulations of Equation (8) with random perturbations added to the stationary solutions demonstrate that the GS is always dynamically stable [5]. The stability also agrees with the anti-Vakhitov–Kolokolov criterion, $d\mu /dN>0$, which is a necessary stability condition for localized states supported by self-repulsive nonlinearities [22] (the original Vakhitov–Kolokolov criterion, $d\mu /dN<0$, is the necessary stability condition in the case of self-attraction [15,23]) .

2.3. The Quantum Phase Transition Induced by the Lee–Huang–Yang (LHY) Correction to the Mean-Field Theory

As seen in Equations (24) and (26), the singularity ~$r^{-1}$ of the stationary wave function at $r\to 0$ suggests that, although the singularity is integrable (as the respective 3D integral for the total norm converges), the LHY correction [20] to the mean-field theory, which is relevant for higher values of the condensate’s density, should be taken into regard. As shown in Refs. [24,25], the scaled GPE with this correction, represented by coefficient $g_{\mathrm{LHY}}>0$, is

$$i{\psi }_t=-\frac{1}{2}\left({\nabla }^2+\frac{U_0}{r^2}\right)\psi +{\left|\psi \right|}^2\psi +g_{\mathrm{LHY}}{\left|\psi \right|}^3\psi .$$

(33)

Then, the asymptotic form of the stationary wave function at $r\to 0$, which was found above in the form determined by the cubic term in Equation (22), $\Phi \left(r\right)\approx \sqrt{U_0/2}{r}^{-1}$, is replaced by

$${\Phi }_{\mathrm{LHY}}\left(r\right)\approx {\left(U_0/2-1/9\right)}^{1/3}{r}^{-2/3},$$

(34)

under the condition of $U_0>2/9$. Meanwhile, at $0<U_0<2/9$, the asymptotic form is determined by the linearization of Equation (33), leading to the same result as given by Equations (17) and (18) with $\sigma ={\sigma }_{-}$ (and $U_l$ replaced by $U_0$):

$${\Phi }_{\mathrm{LHY}}\left(r\right)\approx {\Phi }_0{r}^{-\left(1/2-\sqrt{1/4-U_0}\right)},$$

(35)

where ${\Phi }_0$ is an arbitrary constant in terms of the expansion at $r\to 0$. Note that the wave function with asymptotic form $r^{-{\sigma }_{+}}$, which corresponds to the GS in the linear Schrödinger equation, is incompatible with the presence of the LHY term in Equation (33), although power $-2/3$ in expression (33) coincides at the critical point, $U_0=2/9$, with ${\sigma }_{+}$, rather than ${\sigma }_{-}$.

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 ${\sigma }_{-}$ and ${\sigma }_{+}$), 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 [26] and in many other systems [27,28,29,30,31].

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. This is possible in a binary condensate, with intrinsic self-repulsion in each component, and dominating attraction between them, as proposed in Refs. [24,25], and realized experimentally in the form of “quantum droplets” (in the binary condensate of ${}^{39}$K) in Refs. [32,33,34]. For the symmetric case, with equal wave functions of the two components, the effective GPE takes the form of Equation (33) with the opposite sign in front of the cubic term. This model may be a subject for special consideration.

2.4. The Two-Dimensional Ground State Created by the Quintic Self-Repulsive Nonlinearity

As mentioned above, the GPE in the form of Equation (8) may also be relevant as a physical model in 2D. However, the 2D version of norm (5) of the wave function with asymptotic form ~$r^{-1}$ at $r\to 0$, which follows from this equation (see Equation (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 [35,36].

The 2D GPE can be derived from the underlying 3D version if tight confinement, with the respective harmonic-oscillator length, $a_{\perp }$, is imposed in the z direction by the trapping harmonic-oscillator potential, reducing the effective dimension to that of the plane with remaining coordinates $\left(x,y\right)$ [37,38,39]. In particular, if the dominating quintic terms appear in the 3D GPE with coefficient $g_5$, the reduction to the 2D equation replaces it by ${\left(\sqrt{3}\pi a_{\perp }^2\right)}^{-1}{g}_5$.

In the scaled form, the 2D equation is written in the polar coordinates, $\left(r,\theta \right)$, as

$$i{\psi }_t=-\frac{1}{2}\left({\psi }_{rr}+\frac{1}{r}{\psi }_r+r+\frac{1}{r^2}{\psi }_{\theta \theta }+\frac{U_0}{r^2}\right)\psi +{\left|\psi \right|}^4\psi .$$

(36)

Stationary solutions to Equation (36) (not only the GS, but also for states carrying the angular momentum) are looked for as

$${\psi }_{2\mathrm{D}}\left(r,t\right)=e^{-i\mu t+il\theta }{r}^{-1/2}\chi \left(r\right),$$

(37)

where integer l is the azimuthal quantum number, cf. Equation (24). The substitution of this ansatz in Equation (36) yields an equation for real ${\chi }_{2\mathrm{D}}\left(r\right)$:

$$\mu \chi =-\frac{1}{2}\left[{\chi }^{″}+\left(U_l^{\left(2\mathrm{D}\right)}+\frac{1}{4}\right)r^{-2}\chi \right]+r^{-2}{\chi }^5,$$

(38)

with $U_l^{\left(2\mathrm{D}\right)}\equiv U_0-l^2$, cf. Equation (19). Note that, unlike the 3D case, in the 2D nonlinear model the analysis is possible both for $l=0$ (the GS) and $l\ge 1$.

The expansion of the solution to Equation (38) at $r\to 0$ yields

$$\chi ={\left[\frac{1}{2}\left(U_l^{\left(2\mathrm{D}\right)}+\frac{1}{4}\right)\right]}^{1/4}+{\chi }_1{r}^s,$$

(39)

where $s=\left(1/2\right)\left(1+\sqrt{5+16U_l^{\left(2\mathrm{D}\right)}}\right)$, and ${\chi }_1$ is an arbitrary constant in terms of this expansion, cf. Equation (26) in the 3D case. The solution with a finite norm exists at $U_l^{\left(2\mathrm{D}\right)}>-1/4$, representing at $U_l^{\left(2\mathrm{D}\right)}>0$ the suppression of the collapse and creation of the GS ($l=0$), or the state with $l\ge 1$ by the quintic self-repulsive term.

Combining the 2D asymptotic form (39), valid at $r\to 0$, and the obvious approximation valid at $r\to \infty$, ${\chi }_{2\mathrm{D}}\approx {\chi }_{0}exp\left(-\sqrt{-2\mu }r\right)$, one can derive an interpolation formula for the GS and the dependence $\mu \left(N\right)$ following from it, cf. Equations (28) and (29) in the 3D case:

$$\begin{array}{ccc}\hfill {\psi }_{2\mathrm{D}}\left(r,t\right)& =& {\left[\frac{1}{2}\left(U_l^{\left(2\mathrm{D}\right)}+\frac{1}{4}\right)\right]}^{1/4}{e}^{-i\mu t+il\theta }{r}^{-1/2}{e}^{-\sqrt{-2\mu }r},\hfill \\ \hfill \mu & =& -\left(U_l^{\left(2\mathrm{D}\right)}+\frac{1}{4}\right){\left(\frac{\pi }{2N}\right)}^2.\hfill \end{array}$$

(40)

Similar to the situation in the 3D case, Equation (40) gives an exact wave function with a divergent norm in the limit of $\mu \to -0$,

$${\psi }_{2\mathrm{D}}^{(\mu =0)}\left(r\right)={\left[\frac{1}{2}\left(U_l^{\left(2\mathrm{D}\right)}+\frac{1}{4}\right)\right]}^{1/4}{e}^{il\theta }{r}^{-1/2},$$

(41)

cf. Equation (30). The approximation (40) makes it possible to define the rms (root-mean square) radial size of the two-dimensional GS, cf. Equation (31) in the 3D case:

$$r_{\mathrm{GS}}^{\left(2\mathrm{D}\right)}\equiv \sqrt{\frac{2\pi }{N}{\int }_0^{\infty }{\left|{\psi }_{2\mathrm{D}}\left(r\right)\right|}^2{r}^{3}dr}=\frac{N}{\pi \sqrt{\left(U_l^{\left(2\mathrm{D}\right)}+1/4\right)}}.$$

(42)

Note that the quintic term supports the GS in 2D even at $0<-U_l^{\left(2\mathrm{D}\right)}<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) for $\mu \to 0$, including the case of $0<-U_l^{\left(2\mathrm{D}\right)}<1/4$.

An example of the stable GS, and curves $\mu \left(N\right)$ for the GS in 2D are displayed along with the analytical approximation (40) in Figure 3 (referring to $l=0$, although the replacement of $U_0$ by $U_l$ actually makes no difference in the plots). The $\mu \left(N\right)$ curves are shown for both signs of the central potential, $U_0=-0.18$ and $U_0=0.05$. Simulations of the perturbed evolution in the framework of Equation (36) confirm the stability of the GS families.

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 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 the 3D Equation (33) by $\sqrt{2/5}{\pi }^{-3/4}{a}_{\perp }^{-3/2}{g}_{\mathrm{LHY}}$. Finally, the quartic LHY term determines the asymptotic form of the wave functions at $r\to 0$ as $\sim r^{-2/3}$, which provides for the convergence of the 2D norm. That is, it secures the existence of the GS in the 2D model including the LHY term.

2.5. 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 for Fermi gases the dynamical theory cannot be reduced to a simple mean-field equation [40]. Nevertheless, there is a relatively simple approach to the description of stationary states in a sufficiently dense gas, based on a time-independent equation for the real fermionic wave function, $\Phi \left(\mathbf{r}\right)$ [41,42,43,44], 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 $\mu$—is

$$\mu \Phi =-\frac{1}{3}{\nabla }^2\Phi +{\Phi }^{7/3}-\frac{U_0}{2r^2}\Phi .$$

(43)

The asymptotic form of the solution to Equation (43) at $r\to 0$ is

$$\Phi \left(r\right)=\sqrt{\left(3+4U_0\right)/8}{r}^{-3/2}+O\left(r^{1/2}\right).$$

(44)

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\to 0$. Hence, the relatively weak nonlinearity in Equation (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.

3. The Three-Dimensional Model with Cylindrical Symmetry

The presentation in this section follows the original analysis reported in Ref. [6].

3.1. Formulation of the Model

The previous section addressed the most fundamental spherically symmetric configuration in 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 $\mathbf{d}=d{\mathbf{e}}_z$. This configuration gives rise to the cylindrically symmetric potential of the interaction of the dipolar particle with the fixed attractive center:

$$U\left(\mathbf{r}\right)=-\mathbf{d}·{\mathbf{E}}_Q=-\frac{1}{2}{U}_0{r}^{-2}cos\theta ,$$

(45)

where $cos\theta \equiv z/r$.

If the polarizing external field is electric, it also acts on the central charge. For this reason, a more relevant situation corresponds to the case when the ultracold gas is composed of Hund A-type molecules, with mutually locked electric and magnetic dipoles. Then, an external uniform magnetic field may be employed to align the dipoles in the fixed direction [45].

The 3D GPE with potential (45) is

$$i\frac{\partial \psi }{\partial t}=-\frac{1}{2}\left({\nabla }^2\psi +\frac{U_0}{r^2}cos\theta \right)\psi +{|\psi |}^2\psi .$$

(46)

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:

$$\psi =e^{-i\mu t}{e}^{im\phi }\Phi (r,\theta ),$$

(47)

where the spherical coordinates are used again, cf. Equation (16), and real eigenmode $\Phi$ should be found as a solution of the equation following from the substitution of ansatz (47) in Equation (46):

$$\mu \Phi =-\frac{1}{2}\left[\frac{{\partial }^2}{\partial r^2}+\frac{2}{r}\frac{\partial }{\partial r}-\frac{m^2}{r^2{sin}^2\theta }+\frac{1}{r^{2}sin\theta }\frac{\partial }{\partial \theta }\left(sin\theta \frac{\partial }{\partial \theta }\right)+\frac{U_0}{r^2}cos\theta \right]\Phi +{\Phi }^3.$$

(48)

3.2. The Linear Schrödinger Equation with 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 Equation (48). At $r\to 0$, an asymptotic solution to the linear equation is sought as

$$\Phi \left(r,\theta \right)=r^{-\sigma }{\chi }_{\mathrm{lin}}\left(\theta \right).$$

(49)

The substitution of ansatz (49) in the linearized version of Equation (48) and dropping the term $\mu \Phi$ (which is negligible for the asymptotic analysis at $r\to 0$) leads to an equation that can be written in terms of $\xi \equiv cos\theta$:

$$\frac{d}{d\xi }\left(\left(1-{\xi }^2\right)\frac{d{\chi }_{\mathrm{lin}}}{d\xi }\right)+\left({\sigma }^2-\sigma -\frac{m^2}{1-{\xi }^2}+U_0\xi \right){\chi }_{\mathrm{lin}}\left(\xi \right)=0.$$

(50)

For $U_0=0$, Equation (50) with integer values

$$\sigma =l+1$$

(51)

may be solved in terms of the associated Legendre functions, $l\ge m$ being the orbital quantum number. Note that the singular wave function (49) is 3D-normalizable at $r\to 0$ for $\sigma <3/2$, that is, as it follows from Equation (51), solely for the GS, with $m=l=0$ and $\sigma =1$.

The onset of the quantum collapse is signalled by a transition in Equation (50) from real eigenvalues $\sigma$ to complex ones. Because the effective eigenvalue in the equation is $ϵ\equiv {\sigma }^2-\sigma$ (i.e., $\sigma =(1+\sqrt{1+4ϵ})/2$), the transition to complex $\sigma$ happens at point $ϵ=-1/4$ (i.e., $\sigma =1/2$). At $U_0\ne 0$, Equation (50) cannot be solved in terms of standard special functions. The result of a numerical solution is displayed in Figure 4. It demonstrates that, for given azimuthal quantum number m, with the increase of $U_0$ from zero to some critical value ${\left(U_0\right)}_{\mathrm{cr}}$, eigenvalue $\sigma$ decreases from $\sigma \left(U_0=0\right)=m+1$ to $\sigma \left(U_0={\left(U_0\right)}_{\mathrm{cr}}\right)=1/2$. For the lowest values of m, the numerically found critical values of the potential strength at which $\sigma =1/2$ is attained are

$${\left(U_0\right)}_{\mathrm{cr}}(m=0,1,2)=1.28,7.58,19.06.$$

(52)

Thus, in the framework of the linear Schrödinger equation, the quantum collapse takes place at $U_0>{\left(U_0\right)}_{\mathrm{cr}}\left(m\right)$.

It is relevant to compare critical values (52) of the strength of the axisymmetric potential with those given by Equation (19) for the spherically isotropic one:

$${\left(U_0\right)}_{\mathrm{cr}}^{\left(\mathrm{iso}\right)}(m=0,1,2)=\frac{1}{4}+m(m+1)\equiv 0.25,2.25,6.25.$$

(53)

The comparison naturally shows that the critical strengths are much lower for the spherical potential, which provides a stronger pull to the center.

3.3. Suppression of the Quantum Collapse by the Repulsive Nonlinearity under the Cylindrical Symmetry

As in the isotropic setting, cf. Equation (24), the repulsive cubic term in Equation (48) may balance the attractive potential ~$-r^{-2}$ if, at $r\to 0$, the wave function contains the singular factor $r^{-1}$ (rather than generic $r^{-\sigma }$ in the linear Equation (49)). Then, the substitution of

$$\Phi (r,\theta )=r^{-1}\chi (r,\theta )$$

(54)

transforms Equation (48) into the following equation:

$$\mu \chi =-\frac{1}{2}\left[\frac{{\partial }^2\chi }{\partial r^2}+\frac{1-{\xi }^2}{r^2}\frac{{\partial }^2\chi }{\partial {\xi }^2}-\frac{2\xi }{r^2}\frac{\partial \chi }{\partial \xi }+\left(U_0\xi -\frac{m^2}{1-{\xi }^2}\right)\frac{\chi }{r^2}\right]+\frac{{\chi }^3}{r^2},$$

(55)

(recall $\xi \equiv cos\theta$). Note that Equation (54) makes it possible to write the norm of the 3D wave function as

$$\frac{N}{2\pi }={\int }_0^{\infty }{r}^{2}dr{\int }_0^{\pi }sin\theta d\theta {\left|\psi \left(r,\theta \right)\right|}^2={\int }_0^{\infty }dr{\int }_{-1}^{+1}d\xi {\chi }^2\left(r,\xi \right).$$

(56)

To analyze solutions to Equation (55) at $r\to 0$, one may expand them as

$$\chi (r,\xi )={\chi }_0\left(\xi \right)+{\chi }_1\left(\xi \right)r^{s/2},$$

(57)

assuming $s>0$, which leads to the following equation for ${\chi }_0\left(\xi \right)$, which does not admit an exact solution:

$$\left(1-{\xi }^2\right)\frac{d^2{\chi }_0}{d{\xi }^2}-2\xi \frac{d{\chi }_0}{d\xi }+\left(U_0\xi -\frac{m^2}{1-{\xi }^2}\right){\chi }_0-{\chi }_0^3=0,$$

(58)

cf. Equation (50).

Bound states produced by Equation (55) were found by means of numerical methods in Ref. [6]. Typical profiles of solutions for function $\chi \left(r,\xi \right)$, generated by Equation (55), are displayed in Figure 5 for $m=0,1,2$ and fixed norm $N=2\pi$.

A crude analytical form of the solutions is provided by the Thomas–Fermi approximation (TFA), which neglects all derivatives in Equation (55) [8]:

$${\chi }_{\mathrm{TFA}}^2\left(r,\xi \right)=\left\{\begin{array}{c}\frac{1}{2}{U}_0\xi -\frac{m^2}{2\left(1-{\xi }^2\right)}-|\mu |r^2,\mathrm{at}{r}^2<\frac{1}{2|\mu |}\left(U_0\xi -\frac{m^2}{1-{\xi }^2}\right),\\ 0,\mathrm{at}{r}^2\ge \frac{1}{2|\mu |}\left(U_0\xi -\frac{m^2}{1-{\xi }^2}\right).\end{array}\right.$$

(59)

Actually, this approximation for $m\ge 1$ exists only for $U_0>\left(3\sqrt{3}/2\right)m^2$ (otherwise, Equation (59) yields ${\chi }_{\mathrm{TF}}^2\equiv 0$).

Families of the bound states with different quantum numbers m are presented in Figure 6 by a set of curves showing the chemical potential, $\mu$, versus nonlinearity strength, $U_0$, for a fixed norm ($N=N_0\equiv 2\pi$; producing the results for a fixed norm is sufficient, as the scaling invariance of Equation (55) implies an exact property, $\mu \left(U_0,N\right)={\left(N/N_0\right)}^{-2}\mu \left(U_0,N_0\right)$, the same as mentioned above for the isotropic configuration). Figure 6b,c display the $\mu (U_0$) dependences in relatively narrow intervals of values of $U_0$, stressing that the dependences are obtained above the critical values for the linear Schrödinger equation given by Equation (52), where the linear equation fails to produce any bound state.

TFA based on Equation (59) makes it possible to predict the $\mu \left(U_0\right)$ dependence for the GS ($m=0$) in an analytical form:

$${\mu }_{\mathrm{TFA}}^{\left(\mathrm{GS}\right)}=-\left(2/225\right)U_0^3.$$

(60)

As seen in Figure 6a, 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 of the families of the bound states for $m=0,1,$ and 2.

4. The Two-Component System in Three Dimensions: The Suppression of Quantum Collapse in Miscible and Immiscible Settings

This section summarizes results of the analysis reported in Ref. [7].

4.1. The Formulation of the Model and Analytical Considerations

The generalization of basic model (22) for a binary bosonic gas, with component wave functions ${\psi }_1$ and ${\psi }_2$, is provided by the system of nonlinearly coupled GPEs:

$$\begin{array}{ccc}\hfill i\frac{\partial {\psi }_1}{\partial t}& =& -\frac{1}{2}{\nabla }^2{\psi }_1+(|{\psi }_1{|}^2+\gamma |{\psi }_2{|}^2){\psi }_1-\frac{V_0}{r^2}{\psi }_1,\hfill \\ \hfill i\frac{\partial {\psi }_2}{\partial t}& =& -\frac{1}{2}{\nabla }^2{\psi }_2+(\gamma |{\psi }_1{|}^2+|{\psi }_2{|}^2){\psi }_2-\frac{V_0}{r^2}{\varphi }_2,\hfill \end{array}$$

(61)

where $\gamma$ is the relative strength of the inter-component repulsion, while the coefficients of the self-repulsion are scaled to be 1. To remain consistent with Ref. ([6]), parameter $V_0\equiv U_0/2$ is now used as the strength of the potential pulling particles to the center.

Spherically symmetric bound states with chemical potentials ${\mu }_n<0$, $n=1,2$, of the two components are looked for as

$${\psi }_n\left(r,t\right)=\frac{{\chi }_n\left(r\right)}{r}exp\left(-i{\mu }_{n}t\right),$$

(62)

with real radial functions ${\chi }_n\left(r\right)$ obeying the coupled equations

$$\begin{array}{ccc}\hfill {\mu }_1{\chi }_1& =& -\frac{1}{2}{\chi }_1^{″}-\frac{V_0}{r^2}{\chi }_1+\left({\chi }_1^2+\gamma {\chi }_2^2\right)\frac{{\chi }_1}{r^2},\hfill \\ \hfill {\mu }_2{\chi }_2& =& -\frac{1}{2}{\chi }_2^{″}-\frac{V_0}{r^2}{\chi }_2+\left({\chi }_2^2+\gamma {\chi }_1^2\right)\frac{{\chi }_2}{r^2},\hfill \end{array}$$

(63)

cf. Equations (24) and (25). In terms of these functions, the norms of the components are

$$N_n\equiv \int \left|{\varphi }_n\left(\mathbf{r}\right)\right|d\mathbf{r}=4\pi {\int }_0^{\infty }{\left[{\chi }_n\left(r\right)\right]}^{2}dr,$$

(64)

and the rms radial size of the trapped mode in each component is defined as

$$〈r_n^2〉=\frac{{\int }_0^{\infty }{\left[{\chi }_n\left(r\right)\right]}^2{r}^{2}dr}{{\int }_0^{\infty }{\left[{\chi }_n\left(r\right)\right]}^{2}dr},$$

(65)

cf. Equation (31).

An expansion of solutions to Equations (63) at $r\to 0$ is looked for as

$${\chi }_n\left(r\right)={\chi }_n^{\left(0\right)}\left[1-c_n^{\left(1\right)}{r}^{s/2}-c_n^{\left(2\right)}{r}^{s/2+2}+\cdots -d_n^{\left(1\right)}{r}^2-d_n^{\left(2\right)}{r}^4+\cdots \right],$$

(66)

with $s>0$, cf. Equation (26) (here, $c_1\ne c_2$ is possible, but power s must be the same for ${\chi }_1$ and ${\chi }_2$), which leads to a system of algebraic equations for leading-order coefficients ${\chi }_n^{\left(0\right)}$:

$$\begin{array}{ccc}\hfill {\chi }_1^{\left(0\right)}\left[{\left({\chi }_1^{\left(0\right)}\right)}^2+\gamma {\left({\chi }_2^{\left(0\right)}\right)}^2\right]& =& V_0{\chi }_1^{\left(0\right)},\hfill \\ \hfill {\chi }_2^{\left(0\right)}\left[{\left({\chi }_2^{\left(0\right)}\right)}^2+\gamma {\left({\chi }_1^{\left(0\right)}\right)}^2\right]& =& V_0{\chi }_2^{\left(0\right)}.\hfill \end{array}$$

(67)

Equations (67) give rise to solutions of two types, corresponding to mixed and demixed states in the binary gas:

$$\begin{array}{ccc}\hfill {\chi }_1^{\left(0\right)}& =& {\chi }_2^{\left(0\right)}\equiv {\chi }_{\mathrm{mix}}^{\left(0\right)}=\sqrt{V_0/\left(1+\gamma \right)};\hfill \end{array}$$

(68)

$$\begin{array}{ccc}\hfill {\chi }_1^{\left(0\right)}& \equiv & {\chi }_{\mathrm{demix}}^{\left(0\right)}=\sqrt{V_0},{\chi }_2^{\left(0\right)}=0.\hfill \end{array}$$

(69)

The numerical analysis performed in Ref. [6] demonstrates that demixed modes do not exist at $\gamma <1$, when the mutual repulsion is weaker than the self-repulsive nonlinearity, while mixed ones are completely unstable in the opposite case, $\gamma >1$. Thus, unlike other systems featuring miscibility–immiscibility transitions [46,47,48], in the present situation the transition point is not altered, under the action of the confining potential, in comparison with the commonly known free-space point, $\gamma =1$ [49].

Further analysis demonstrates a change in the structure of the r-dependent corrections in Equation (66) for the miscible system, with $\gamma <1$: at $V_0<1/2$, the dominant terms are $\sim r^{\left(1+\sqrt{1+16V_0}\right)/2}$, while at $V_0>1/2$ these are terms $\sim 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\equiv {\left(V_0\right)}_{\mathrm{phase}-\mathrm{trans}}$, expansion (66) is replaced by

$${\chi }_n\left(r\right)=\frac{1}{\sqrt{2\left(1+\gamma \right)}}\left[1+\frac{{\mu }_1+{\mu }_2}{4}{r}^{2}ln\left(\frac{r_0}{r}\right)+{(-1)}^n\frac{1+\gamma }{4\gamma }\left({\mu }_1-{\mu }_2\right)r^2\right].$$

(70)

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 Equations (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 [50].

Lastly, at $r\to \infty$, Equations (63) yield an exponential asymptotic form of the solution,

$${\chi }_n\left(r\right)\approx {\chi }_n^{(\infty )}\left(1-\frac{V_0}{\sqrt{-2{\mu }_n}r}\right)exp\left(-\sqrt{-2{\mu }_n}r\right),$$

(71)

where constants ${\chi }_n^{(\infty )}$ are indefinite in terms of the asymptotic expansion at $r\to \infty$.

4.2. Numerical and Additional Analytical Results for Trapped Binary Modes

4.2.1. Mixed Ground States

Figure 7a shows a typical profile for the mixed GS produced by a numerical solution of Equation (63) at $V_0=1$ for $\gamma =0.9$ and equal norms of the two components, $N_1=N_2=4\pi$.

The simplest global analytical approximation for the GS wave function is provided by the interpolation, similar to that introduced in the single-component setting, cf. Equation (28):

$${\chi }_n\left(r\right)\approx {\chi }_{\mathrm{mix}}^{\left(0\right)}{e}^{-\sqrt{-2{\mu }_n}r}.$$

(72)

The substitution of this interpolation in Equations (64) and (65), along with expression (68), leads to predictions for the chemical potentials and rms radius of the two components as functions of their norms (which are also valid in the case of $N_1\ne N_2$):

$$\begin{array}{ccc}\hfill {\mu }_n& =& -2{\left[\frac{\pi V_0}{(1+\gamma )N_n}\right]}^2,\hfill \end{array}$$

(73)

$$\begin{array}{ccc}\hfill 〈r_n^2〉& =& {\left[\frac{(1+\gamma )N_n}{2\pi V_0}\right]}^2.\hfill \end{array}$$

(74)

Comparison of expression (73) with numerical results is shown in Figure 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\equiv N$, which yields (for ${\chi }_1={\chi }_2\equiv \chi$):

$${\chi }_{\mathrm{TFA}}\left(r\right)=\left\{\begin{array}{c}\sqrt{\left(V_0+\mu r^2\right)/\left(1+\gamma \right)},\mathrm{at}r<R_0\equiv \sqrt{V_0/(-\mu )},\\ 0,\mathrm{at}r>R_0,\end{array}\right.$$

(75)

cf. TFA for the potential with the cylindrical symmetry, given by Equation (59). The substitution of approximation (75) in Equations (5) and (65) yields the predictions for the chemical potential and effective size of the GS:

$$\begin{array}{ccc}\hfill {\mu }_{\mathrm{TFA}}& =& -\frac{64{\pi }^2{V}_0^3}{9{(1+\gamma )}^2{N}^2},\hfill \end{array}$$

(76)

$$\begin{array}{ccc}\hfill 〈r_{\mathrm{TFA}}^2〉& =& \frac{5}{{\pi }^3}{\left[\frac{3\left(1+\gamma \right)N}{16V_0}\right]}^2\equiv \frac{5}{4\pi }{R}_0^2,\hfill \end{array}$$

(77)

(recall $R_0$ is the TFA cutoff radius defined in Equation (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 Figure 7b. 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 Figure 3a at $V_0=2$ and $\gamma =0.9$ for $N_1=4\pi$ and $N_2=2\pi$. The imbalanced mixed states with ${\mu }_1\ne {\mu }_2$ and $N_1\ne N_2$ feature equal values of ${\chi }_{1,2}(r=0)$, in agreement with Equation (68).

In the case of the strong pull to the center, $V_0\gg 1$, TFA can be generalized for imbalanced states, fixing $\left|{\mu }_1\right|\le \left|{\mu }_2\right|$ 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. [51]. In the inner layer,

$$r^2<r_0^2\equiv \frac{1-\gamma }{\gamma {\mu }_1-{\mu }_2}{V}_0,$$

(78)

both wave functions are different from zero:

$${\chi }_n^{\left(\mathrm{inner}\right)}\left(r\right)=\sqrt{\frac{V_0}{1+\gamma }-\frac{\gamma {\mu }_{3-n}-{\mu }_n}{1-{\gamma }^2}{r}^2}.$$

(79)

In the outer layer, only one component is present, in the framework of TFA: ${\chi }_2\equiv 0$,

$${\chi }_1^{\left(\mathrm{outer}\right)}\left(r\right)=\left\{\begin{array}{c}\sqrt{V_0+{\mu }_1{r}^2},\mathrm{at}{r}_0^2\le r^2\le R_0^2\equiv -V_0/{\mu }_1,\\ 0,\mathrm{at}{r}^2\ge R_0^2.\end{array}\right.$$

(80)

Both components of the TFA solution given by Equations (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 Figure 8b,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 Equation (72), produces a sufficiently accurate prediction in that case.

4.2.2. The Immiscible Ground State

As said above, in the case of $\gamma >1$, relevant states are immiscible ones. The two-layer TFA may be applied to produce an immiscible GS. In the inner layer,

$$r^2<r_0^2=\frac{(\gamma -1)V_0}{\gamma {\mu }_1-{\mu }_2},$$

the approximation yields

$${\chi }_1\left(r\right)=\sqrt{V_0+{\mu }_1{r}^2},{\chi }_2\left(r\right)=0.$$

(81)

In the outer layer, which is $r_0^2<r^2<R_0^2=V_0/(-{\mu }_2)$, the result is

$${\chi }_1\left(r\right)=0,{\chi }_2\left(r\right)=\sqrt{V_0+{\mu }_2{r}^2}.$$

(82)

That is, 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.

5. The Mean-Field Predictions versus the Many-Body Quantum Theory

5.1. Introduction to the Section

The analysis presented above was performed in Refs. [5,6,7] in the framework of the mean-field theory, that is, the respective GPEs (possibly including the beyond-mean-field LHY corrections, see Equation (33)). A relevant issue is the comparison of the basic mean-field predictions (e.g., 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\left(r\right)=-U_0/r^2$. That is, $U_0$ in Equation (1) is replaced by $2U_0$ to make the notation consistent with that in Ref. [21], which addressed the present problem. Results produced in that work are recapitulated in the present section.

The many-body Hamiltonian representing the setting under consideration is

$$\widehat{H}=-\sum _{j=1}^N\left(\frac{{\hslash }^2{\nabla }_j^2}{2m}+\frac{U_0}{r_j^2}\right)+\sum _{j<k}^N{V}_{\mathrm{int}}(|{\mathbf{r}}_j-{\mathbf{r}}_k|),$$

(83)

where ${\mathbf{r}}_j$ are coordinates of the j-th particle in the 3D space, m is the particle’s mass, and $V_{\mathrm{int}}\left(r\right)$ is the potential of the repulsive interaction between the particles. In the framework of the mean-field theory, $V_{\mathrm{int}}\left(r\right)$ 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 Equations (88) and (89).

Before introducing the many-body wave function, the single-particle one is adopted as per the following ansatz:

$$f_1\left(r\right)=r^{\beta }exp(-\alpha r^2),$$

(84)

where $\alpha \ge 0$ determines the inverse localization length, which affects the system’s size and, consequently, the density. Alternatively, $\alpha$ can be interpreted in terms of an effective external harmonic confinement with frequency $\Omega =2\alpha \hslash /m$, cf. Equation (4). At $r_j\to 0$, the shape of the wave function is controlled by parameter $\beta$ in ansatz (84).

5.2. 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 $\alpha$ and $\beta$ 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_{\mathrm{int}}(|{\mathbf{r}}_i-{\mathbf{r}}_j|)$, does not appear. The variational energy, $E_{\mathrm{var}}^{\left(1\right)}={\left[\int f_1^2\left(\mathbf{r}\right)d\mathbf{r}\right]}^{-1}\int f_1\left(\mathbf{r}\right)H{f}_1\left(\mathbf{r}\right)d\mathbf{r}$, with $f_1$ taken as per Equation (84), is

$$E^{\left(1\right)}=\alpha \left[1-\frac{8U_0-1}{2(1+2\beta )}\right].$$

(85)

For a fixed localization size, $\alpha =\mathrm{const}$, this energy is a decreasing function of $\beta$ if $U_0$ is smaller than the critical value for the onset of the collapse, $U_0=1/8$, which is tantamount to one which appears in Equation (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].

In the framework of the local-density approximation, the chemical potential of the state with uniform density n is taken as ${\mu }_{hom}=gn$, where $g=4\pi {\hslash }^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 ${\mu }_{\mathrm{loc}}=gn$, where this time n is a function of the coordinates, rather than a constant. Additionally, the external potential,

$$\mu ={\mu }_{\mathrm{loc}}-\frac{U_0}{r^2}+\frac{1}{2}m{\Omega }^2{r}^2,$$

(86)

where the harmonic-oscillator confinement with the respective length scale $a_{\mathrm{ho}}=\sqrt{\hslash /(m\Omega )}$ is added to make the size of the system finite, cf. potential (4) used above. Solving Equation (86) for the density, one obtains the following density profile:

$$n\left(r\right)=\frac{1}{g}\left\{\begin{array}{cc}\mu -\frac{1}{2}m{\Omega }^2{r}^2+U_0{r}^{-2},\hfill & \mathrm{at}r<R_{\mathrm{TFA}},\hfill \\ 0,\hfill & \mathrm{at}\ge R_{\mathrm{TFA}},\hfill \end{array}\right.$$

(87)

where the radius of the gaseous cloud is taken as per TFA, $R_{\mathrm{TFA}}=\sqrt{\mu +\sqrt{{\mu }^2+2m{U}_0{\Omega }^2}}/(\sqrt{m}\Omega )$. The density at the center features an integrable divergence in Equation (87), reflecting the presence of the attractive central potential, cf. Equation (24). Finally, the chemical potential itself is fixed by the normalization condition, $4\pi {\int }_0^{R_{\mathrm{TFA}}}n\left(r\right)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. [21], viz., the hard-sphere potential of diameter R,

$$V_{\mathrm{hard}}\left(r\right)=\left\{\begin{array}{cc}\infty ,\hfill & r<R\hfill \\ 0,\hfill & r\ge R\hfill \end{array}\right.,$$

(88)

and its soft-sphere counterpart,

$$V_{\mathrm{soft}}\left(r\right)=\left\{\begin{array}{cc}{V}_0,\hfill & r<R\hfill \\ 0,\hfill & r\ge R\hfill \end{array}\right.,$$

(89)

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

$$a_s=R[1-tanh\left(kR\right)/\left(kR\right)],$$

(90)

where the momentum corresponding to the height of the soft-sphere potential is

$$k\equiv \sqrt{m{V}_0}/\hslash .$$

(91)

For hard-sphere potential (88), the effective s-wave scattering length is identical to the diameter of the sphere, $a_s=R$.

5.3. 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. [21], the variational Monte-Carlo method was employed, which samples the probability distribution, $p={|\psi |}^2$, for a known many-body wave function, $\psi$, allowing one to calculate the variational energy as a function of trial parameters, such as $\alpha$ and $\beta$ in Equation (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\left(r\right)$, taken as per Equation (84), and a pairwise product of two-particle Jastrow terms [53], $f_2\left(r\right)$:

$$\psi ({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N)=\prod _{j=1}^N{f}_1\left(r_j\right)\prod _{j<k}^N{f}_2(|{\mathbf{r}}_j-{\mathbf{r}}_k|).$$

(92)

The Jastrow factor $f_2\left(r\right)$ in Equation (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

$$f_2^{\left(\mathrm{hard}\right)}\left(r\right)=\left\{\begin{array}{cc}0,\hfill & r<R\hfill \\ 1-R/r,\hfill & r\ge R\hfill \end{array}\right.,$$

(93)

while for the soft-sphere potential (89), it is [1]

$$f_2^{\left(\mathrm{soft}\right)}\left(r\right)=\left\{\begin{array}{cc}Asinh\left(kr\right)/r,\hfill & r<R\hfill \\ 1-a_s/r,\hfill & r\ge R\hfill \end{array}\right.,$$

(94)

where k is given by Equation (91), and constant A is determined by the condition of the continuity of $f_2\left(r\right)$ at $r=R$.

5.4. 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=10,000$. For small values of $\alpha$ (which corresponds to weak localization), the energy may be negative. This is not visible in the log-log plot of Figure 10. As the localization gets tighter, the energy becomes positive, as the two-body interaction helps the system to resist the trend of collapsing. For very tight localization, $\alpha \to \infty$, the collapse is observed for small values of N, with the energy diverging towards $-\infty$.

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, $\hslash \Omega$, 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

$$E=CN\alpha ,$$

(95)

with $C=4$. It is shown in Figure 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

$$E=N{E}^{\left(1\right)}+N(N-1)V_0.$$

(96)

The energy of the interparticle interactions, revealed by the calculations, is $E_{\mathrm{int}}=V_{0}N(N-1)/2$. The asymptotic energy (96) is shown in Figure 10 by dashed lines.

A clear conclusion is that the increase of the number of particles indeed causes a strong rise of the potential barrier, which stabilizes the metastable energy minimum corresponding to the gaseous state. This can also be concluded from Equation (96), where the contribution due to the repulsive interactions scales as $N^2$ for large N, while the term corresponding to 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_{\mathrm{barrier}}$, is estimated as the maximum value of the energy per particle (see Figure 10), and is shown in Figure 11. For a large system’s size, the barrier can be approximated by comparing the two basic energy scales given by Equations (95) and (96). The resulting asymptotic approximation for the barrier’s height is

$$E_{\mathrm{barrier}}=6N{V}_0/(8U_0-3+2C),$$

(97)

which is shown in Figure 11 by the dashed line.

6. Discussion and Conclusions

This article aims to produce a review of results reported in works [5,6,7,21] 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/\left(2r^2\right)$. 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\to 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 do not give rise to losses) is 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 also considered, 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-dependent GPEs.

In work [21], 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 also exist in the full many-body setting, 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 the possibility of the spontaneous symmetry breaking of the GS in the latter case.

As mentioned above (see Equations (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/\left(2r^2\right)$.