Angular-Momentum Modes in a Bosonic Condensate Trapped in the Inverse-Square

1. Introduction and the Model

It is well known that attractive potential $U\left(r\right)=-U_0/\left(2r^2\right)$, with $U_0>0$, plays a special role in quantum mechanics [1]. Indeed, unlike the classical equation of motion in the same potential, which gives rise to equivalent solutions for all values of $U_0$, due to the scaling invariance of the equation, $U_0$ cannot be scaled out in the corresponding Schrödinger equation, written in the normalized form,

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

(1)

where $U_0$ is an irreducible parameter. It is known that the three-dimensional (3D) equation (1) gives rise to the spherically symmetric ground state (GS) if the strength of the pull to the center does not exceed a critical value,

$$U_0<{\left(U_0\right)}_{\mathrm{cr}}=1/4,$$

(2)

while at $U_0>1/4$ Eq. (1) does not maintain the GS, leading, instead, to the quantum collapse of wave function $\psi$ (called “fall onto the center" of the quantum particle in the book by Landau and Lifshitz [1]). Exact solutions of Eq. (1) for the developing collapse, as well as the respective solutions for the classical equation of motion, were recently analyzed in Ref. [2]. Such a difference between the classical and quantum behavior under the action of the same potential is considered as the quantum anomaly, alias “dimensional transmutation" [3,4] (quantum anomalies for particles trapped in an external potential are also known in the 2D geometry [5]).

A physical realization of the Schrödinger equation (1) is possible for a particle with a permanent electric dipole moment (such as a small polar molecule) pulled to the electric charged fixed at the center. If the particle minimizes its energy by aligning the polarization of the dipole moment with the local electric field, the attraction potential takes precisely the form adopted in Eq. (1) [7]. In this case, for the particle with the mass $\sim 100$ proton masses and the elementary electric charge set at the center, the critical value $U_0=1/4$ corresponds to a very small dipole moment, $\sim {10}^{-6}$ Debye. Therefore, the overcritical case of $U_0>1/4$ is relevant in the actual physical context.

A possible solution of the quantum-collapse problem, i.e., restoration of the missing GS in the case of $U_0>1/4$, was proposed in terms the secondary quantization, replacing the Schrödinger equation (1) by the corresponding linear quantum-field theory [3,4]. However, the solution does not predict the size of the so introduced GS. Instead, the renormalization-group technique, on which the field-theory formulation is based, postulates the existence of a GS with an arbitrary spatial size, in terms of which all other spatial scales are measured in the framework of the theory.

A completely different, many-body, solution of the collapse problem was proposed in Ref. [7]. Instead of addressing the single particle, it deals with an ultracold gas of particles pulled to the center, in the form of a Bose-Einstein condensate (BEC). In particular, ultracold gases of polar molecules such as LiCs [8] and KRb [9] are available to the experiment, and it was demonstrated too that a fixed electric charge (ion) can be embedded in BEC and held at a fixed position by means of an optical trapping scheme [10].

The suppression of the quantum collapse in this setting is secured by repulsive collisional interactions between the particles. The analysis was performed in the framework of the mean-field (MF) approximation [11], using the Gross-Pitaevskii equation (GPE) including the same potential as in Eq. (1):

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

(3)

In the framework of the MF approximation, the cubic term in Eq. (3) represents the self-repulsion in the gas. The coefficient in front of this term is set equal to 1 by means of normalization. Dipole-dipole interactions between the particles were also taken into account, applying another version of the MF approximation. It considered the interaction of a given dipole with the collective electric potential created by all other dipoles, and amounted to renormalization of coefficient $U_0$.

For a mixture of two species of particles, a two-component extension of Eq. (3) was introduced in Ref. [12]. The respective system of nonlinearly coupled GPEs also provides the suppression of the quantum collapse and creation of the GS [12].

Solutions of Eq. (3) for stationary states with chemical potential $\mu <0$, in the form of

$$\psi =e^{-i\mu t}u\left(x,y,z\right),$$

(4)

with steady-state wave function u satisfying the equation

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

(5)

These states are characterized by their norm,

$$N=\int \int \int {\left|u\left(x,y,z\right)\right|}^{2}dxdydz,$$

(6)

and angular momentum,

$$M=i\int \int \int {\psi }^{*}\left(x,y,z\right)\left(\mathbf{r}\times \nabla \right)\psi \left(x,y,z\right)dxdydz,$$

(7)

where * stands for the complex conjugate. N and M are dynamical invariants of GPE (3). It also conserves the Hamiltonian,

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

(8)

The linearized version of Eq. (3) produces exact, although unnormalizable (with diverging norm (6)), pairs of eigenmodes with $\mu =0$ and orbital and magnetic quantum numbers $l,m$ [7],

$$u=u_0{r}^{-{\sigma }_{\pm }}{\mathrm{Y}}_{lm}\left(\theta ,\phi \right),{\sigma }_{\pm }=\frac{1}{2}\pm \sqrt{\frac{1}{4}-U_l},$$

(9)

which are written in spherical coordinates $\left(r,\theta ,\phi \right)$. Here $u_0$ is an arbitrary amplitude, ${\mathrm{Y}}_{lm}$ are the spherical harmonics, and

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

(10)

Solution (9) implies that these eigenstates exists at $U_l<1/4$, cf. Eq. (2), while the collapse occurs at $U_l>1/4$.

The suppression of the collapse and recreation of the GS (with a normalizable wave function, unlike the unnormalizable one (9)) in the framework of the full nonlinear GPE (3) was demonstrated for the spherically isotropic state, corresponding to $l=m=0$ in expression (9) (i.e., the s-wave orbital, in terms of atomic physics [1]). To this end, Eq. (5) is converted, by substitution

$$u\left(r\right)=r^{-1}v\left(r\right)$$

(11)

with real function $v\left(r\right)$, into a radial equation

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

(12)

An asymptotic solution to Eq. (12) at $r\to 0$ is

$$v\left(r\right)=\sqrt{U_0/2}+\mathrm{const}·r^{s/2},s\equiv 1+\sqrt{1+8U_0}$$

(13)

(here, $\mathrm{const}$ is not determined by the asymptotic expansion). At $r\to \infty$, the bound-state wave function with $\mu <0$ decays exponentially, $v\sim exp\left(-\sqrt{-2\mu }r\right)$. A global approximation for the GS can be constructed as an interpolation, juxtaposing the asymptotic forms which are valid at $r\to 0$ and $r\to \infty$:

$${\left({\psi }_{\mathrm{GS}}(r,t)\right)}_{\mathrm{interpol}}\approx \sqrt{\frac{U_0}{2}}{e}^{-i\mu t}{r}^{-1}{e}^{-\sqrt{-2\mu }r}.$$

(14)

Singularity $\sim r^{-1}$ of wave function (14) at $r\to 0$ is acceptable, as the corresponding 3D norm integral (6) converges at $r\to 0$. In particular, the approximate wave function (14) gives rise to the following relation between the chemical potential and norm,

$${\mu }_{\mathrm{GS}}\approx -\frac{1}{2}{\left(\frac{\pi U_0}{N_{\mathrm{interpol}}}\right)}^2.$$

(15)

Actually, scaling $\mu \sim -1/N^2$, which is produced by Eq. (15), is an exact property of solutions to Eq. (5), irrespective of validity of the approximation. In the limit of $\mu \to -0$, Eq. (14) produces an exact solution of Eq. (3), ${\psi }_{\mu =0}\left(r\right)=\sqrt{U_0/2}{r}^{-1}$, with a divergent norm.

The numerical solution of Eq. (5) corroborates the conclusion following from Eqs. (13) – (15): GPE (3) does not give rise to the collapse at $U_0>1/4$, and maintains the well-defined GS at all values of $U_0$. For the setting with the elementary electric charge fixed at the center, and particles with the mass $\sim 100$ proton masses and electric dipole moment $\sim 1$ Debye, the radius of the newly predicted GS, which replaces the collapsing regime in the gas of $\sim {10}^5$ particles, is $r_{\mathrm{GS}}\simeq 3$$\mu$m [7]. Note that the MF approximation is relevant for states characterized by such a length scale. Furthermore, beyond the framework of the MF approximation, the many-body theory demonstrates that the predicted state persists as a meta-stable one, separated by a tall potential barrier from the collapsing state (the many-body analysis does not allow full suppression of the collapse) [13].

Relation (15) satisfies the anti-Vakhitov-Kolokolov criterion, $d\mu /dN>0$, which is a necessary condition for stability of families of bound states supported by a self-repulsive nonlinearity [14] (the Vakhitov-Kolokolov criterion per se, $d\mu /dN<0$, pertains to the case of the self-attraction [15,16]), Full stability of the GS created by Eq. (3) for all considered values of $U_0$ was corroborated by systematic simulations of the perturbed evolution of the bound states [7].

An obviously interesting extension of the analysis outlined above, which is the subject of the present work, is to develop it for bound states with reduced symmetry, which carry angular momentum. Previously, states with embedded vorticity, which is directly related to the angular momentum, were considered in a modified setting, with the polarization of dipole moments fixed by an external uniform electric field, directed along the z axis [17]. In that case, the spherically symmetric pulling potential in Eq. (3) is replaced by an axially symmetric one, $-\left(U_0/2r^2\right)cos\theta$. Critical values of $U_0$ above which the linear axially-symmetric model gives rise to the quantum collapse of states with magnetic quantum numbers (vorticities) $m=0$ (i.e., GS), 1 and 2 are, respectively, ${\left(U_0\right)}_{\mathrm{cr}}^{\left(m\right)}=1.28$, $7.58$, and $19.06$. For comparison, Eqs. (2) and (33) yield the set of smaller critical values, viz. ${\left(U_0\right)}_{\mathrm{cr}}^{\left(l\right)}=0.25$, $2.25$, and $6.25$ for $l=0$, 1, and 2, respectively. The nonlinear model with the axisymmetric potential also suppresses the collapse and creates stable bound states with vorticity m at $U_0>{\left(U_0\right)}_{\mathrm{cr}}^{\left(m\right)}$.

Vortex states were also addressed in the 2D version of the setting, whose linear version gives rise to the collapse at all $U_0>0$. In that case, the cubic self-repulsion in GPE is not sufficient to suppress the collapse [7]. A physically relevant alternative is provided by the quartic term, which is produced by the Lee-Huang-Yang effect, i.e., a correction to the MF theory induced by quantum fluctuations [18]. To this end, one may consider a binary BEC in which the MF self-repulsion in each component is (almost) precisely balanced by the inter-component attraction [19,20]. The resulting amended GPE predicts self-trapped states in the form of quantum droplets, which have been observed experimentally with quasi-2D [21,22] and full 3D [23,24] shapes in binary BECs with the contact interactions, as well as in dipolar BECs with long-range interactions between atomic magnetic dipole moments [25,26].

In the present context, the corresponding 2D GPE takes the form of Eq. (3), with term ${\left|\psi \right|}^2\psi$ replaced by ${\left|\psi \right|}^3\psi$ [27], leading to the solution for the wave function with the following asymptotic form at $r\to 0$, which replaces the collapsing solution of the 2D linear equation:

$$\psi \approx {\left[\frac{1}{2}\left(U_m+\frac{4}{9}\right)\right]}^{1/3}{e}^{i\mu t-im\phi }{r}^{-2/3},U_m\equiv U_0-m^2$$

(16)

(cf. expression (14) at $r\to 0$ and Eq. (10)), where $\left(r,\phi \right)$ are the 2D polar coordinates, and integer m is the vorticity. Note that the 2D norm of this wave function with the acceptable singularity $r^{-2/3}$ converges at $r\to 0$. A nontrivial problem is stability of the 2D vortex bound states corresponding to the asymptotic form (16), which exist under the condition of $U_0>m^2-4/9\equiv {\left(U_0\right)}_{min}$. Surprisingly, this problem admits an exact solution: the vortex states with $m\ge 1$ are stable provided that $U_0>(7/9)\left(3m^2-1\right)\equiv {\left(U_0\right)}_{\mathrm{stability}}$, the GS with $m=0$ being stable too [27]. Note that, for all $m\ge 1$, ${\left(U_0\right)}_{\mathrm{stability}}$ exceeds ${\left(U_0\right)}_{min}$, hence the 2D vortex states are unstable in the interval of $m^2-4/9<U_0<$$(7/9)\left(3m^2-1\right)$. The dominant instability mode is one which leads to slow drift of the vortex’ pivot off the central point, along a spiral trajectory [27].

The objective of the present work is to construct 3D states, carrying angular momentum, as solutions of GPE (3 with nonzero orbital quantum number, $l\ge 1$, and the magnetic quantum number taking values $0\le m\le l$. Essential results for this problem can be obtained in an analytical form, as shown below in Section 2. In particular, an exact existence threshold is found for all nonlinear states with $l\ge 1$, viz.,

$$U_0\ge {\left(U_0^{\left(l\right)}\right)}_{\mathrm{thr}}=l(l+1)$$

(17)

(it does not depend on m). Numerical results, which essentially corroborate the analytical predictions, are reported in Section 3, and the paper is concluded by Section 4. To the best of our knowledge, considerations of GPE-based models in similar 3D settings have not been reported before.

2. Analytical Considerations

2.1. The Framework for the Analysis of the Gross-Pitaevskii Equation

It is relevant to explicitly write the underlying 3D equation (3) in the spherical coordinates:

$$\begin{array}{c}i\frac{\partial \psi }{\partial t}=-\frac{1}{2}\left[\frac{{\partial }^2}{\partial r^2}+\frac{2}{r}\frac{\partial }{\partial r}\right.\\ \left.+\frac{1}{r^2}\left(\frac{{\partial }^2}{\partial {\theta }^2}+cot\theta ·\frac{\partial }{\partial \theta }+U_0+\frac{1}{{sin}^2\theta }\right)\frac{{\partial }^2}{\partial {\phi }^2}\right]\psi +{\left|\psi \right|}^2\psi .\end{array}$$

(18)

Stationary solutions to Eq. (18), with chemical potential $\mu <0$ and vorticity which is represented by the magnetic quantum number m, are looked for as

$$\psi =exp\left(-i\mu t+im\phi \right)u^{\left(m\right)}\left(r,\theta \right),$$

(19)

cf. Eq. (4), with real function $u^{\left(m\right)}\left(r,\theta \right)$ satisfying the stationary equation,

$$\begin{array}{c}\mu u^{\left(m\right)}=-\frac{1}{2}\left[\frac{{\partial }^2}{\partial r^2}+\frac{2}{r}\frac{\partial }{\partial r}\right.\\ \left.+\frac{1}{r^2}\left(\frac{{\partial }^2}{\partial {\theta }^2}+cot\theta ·\frac{\partial }{\partial \theta }+\left(U_0-\frac{m^2}{{sin}^2\theta }\right)\right)\right]u^{\left(m\right)}+{\left(u^{\left(m\right)}\right)}^3.\end{array}$$

(20)

To eliminate the singularity, $\sim 1/r$, of the solution, we perform substitution (11), i.e.,

$$u^{\left(m\right)}\left(r,\theta \right)\equiv r^{-1}{V}^{\left(m\right)}\left(r,\theta \right),$$

(21)

in Eq. (20), which leads to an equation for $V^{\left(m\right)}\left(r,\theta \right)$ [cf. Eq. (12)],

$$2\mu r^2{V}^{\left(m\right)}=-r^2\frac{{\partial }^2{V}^{\left(m\right)}}{\partial r^2}-\left(\frac{{\partial }^2}{\partial {\theta }^2}+cot\theta ·\frac{\partial }{\partial \theta }+\left(U_0-\frac{m^2}{{sin}^2\theta }\right)\right)V^{\left(m\right)}+2{\left(V^{\left(m\right)}\right)}^3.$$

(22)

To analyze the form of eigenmodes in the limit of $r\to 0$, one may define, straightforwardly,

$$V^{\left(m\right)}\left(r=0,\theta \right)\equiv v^{\left(m\right)}\left(\theta \right),$$

(23)

and directly set $r=0$ in Eq. (22), which leads, without any approximation, to the angular ordinary differential equation:

$$\frac{d^2{v}^{\left(m\right)}}{d{\theta }^2}+cot\theta ·\frac{d{v}^{\left(m\right)}}{d\theta }+\left(U_0-\frac{m^2}{{sin}^2\theta }\right)v^{\left(m\right)}=2{\left(v^{\left(m\right)}\right)}^3.$$

(24)

Equation (24) should be solved in the natural region of $0\le \theta \le \pi$. In fact, Eq. (23) implies that relevant solutions to Eq. (24) provide the boundary condition for solutions to Eq. (22) at $r\to 0$. The boundary condition at $r\to \infty$ amounts to the exponential decay of the wave function:

$$V^{\left(m\right)}\left(r\to \infty ,\theta \right)\sim exp\left(-\sqrt{-2\mu }r\right).$$

(25)

The artificial singularity of Eq. (24) at $\theta =0$ and $\pi$ suggests a possibility of the existence of singular solutions with the asymptotic form

$$v_{\sin \mathrm{g}}^{\left(m\right)}\left(\theta \right)\approx v_0^{\left(m\right)}{\theta }^{-1},\mathrm{or}{v}_{\sin \mathrm{g}}^{\left(m\right)}\left(\theta \right)\approx v_0^{\left(m\right)}{(\pi -\theta )}^{-1}.$$

(26)

The substitution of this ansatz in Eq. (24) yields

$${\left(v_0^{\left(m\right)}\right)}^2=\frac{1}{2}\left(1-m^2\right).$$

(27)

Obviously, it follows from Eq. (27) that the singular solution may exist only at $m=0$ (in particular, for $m=0$ and $U_0=0$ Eq. (24) has an exact singular solution, $v^{\left(0\right)}\left(U_0=0\right)={\left(\sqrt{2}sin\theta \right)}^{-1}$). Actually, the singular solution (26) is irrelevant, as the respective normalization integral includes a divergent factor, ${\int }_0^{\pi }{\left[v_{\sin \mathrm{g}}^{(m=0)}\left(\theta \right)\right]}^{2}sin\theta d\theta =\infty$.

2.2. Identifying Vortex (Sectoral) States with Equal Magnetic and Orbital Quantum Numbers, $m=l$

Vortex-state solutions of Eq. (22) with $m=l$, also known as sectoral ones, are determined by Eq. (24) with boundary conditions at $\theta =0$ and $\pi$

$$v^{\left(m\right)}\sim {\theta }^m,v\sim {(\pi -\theta )}^m$$

(28)

(by definition, m does not take negative integer values). Equation (24) with these boundary conditions admits an exact solution, suggested by the form of spherical harmonic ${\mathrm{Y}}_{ll}$, in the form of

$$v^{(m=l)}\left(\theta \right)=v_0^{(m=l)}{\left(sin\theta \right)}^m,$$

(29)

with an infinitesimal amplitude $v_0^{\left(m\right)}$, at the threshold value of the strength of the pulling potential,

$${\left(U_0^{(m=l)}\right)}_{\mathrm{thr}}=m(m+1)\equiv l(l+1),$$

(30)

which is an example of the general threshold value (17). Above the threshold, i.e., at

$$U_0>m(m+1),$$

(31)

amplitude $v_0^{\left(m\right)}$ in Eq. (29) grows from infinitesimal to finite values. This fact has a simple explanation: because the orbital kinetic energy of the state with quantum numbers $l=m$ is

$$E_{\mathrm{kin}}=m(m+1)/\left(2r^2\right),$$

(32)

condition (31) implies that the energy of pulling to the center exceeds the kinetic (centrifugal) energy.

In the simplest case, $m=l=1$ (the p-wave sectoral state, in terms of atomic physics), one can substitute expression (29),

$$v^{(m=l=1)}\left(\theta \right)=v_0^{(m=l=1)}sin\theta ,$$

(33)

as an approximation, in the nonlinear term of Eq. (24), which yields

$${\left(sin\theta \right)}^3=\frac{1}{4}\left[3sin\theta -sin\left(3\theta \right)\right].$$

(34)

Neglecting, as usual, the third harmonic in Eq. (34) and keeping the term $\sim sin\theta$, for small values of the distance from the threshold, viz.,

$$U_0-{\left(U_0^{(m=l=1)}\right)}_{\mathrm{thr}}\equiv U_0-2\ll 1,$$

(35)

the so approximated Eq. (24) predicts the squared amplitude of the approximate solution (29) as

$${\left(v_0^{(m=l=1)}\right)}^2\approx \frac{2}{3}\left(U_0-2\right).$$

(36)

Similarly, for the vortex (sectoral) state with $m=l=2$ (the d-wave orbital, in terms of atomic physics), the exact solution (29) with an infinitesimal amplitude $v_0^{(m=l=2)}$,

$$v^{(m=l=2)}\left(\theta \right)=v_0^{(m=l=2)}{sin}^2\theta ,$$

(37)

exists at $U_0=6$, as given by Eq. (30) with $l=2$. To predict the amplitude of this state above the threshold, one can use the following approximation:

$${sin}^6\theta \equiv \frac{15}{16}{sin}^2\theta -\frac{3}{8}{sin}^2\left(2\theta \right)+\frac{1}{16}{sin}^2\left(3\theta \right)\approx \frac{15}{16}{sin}^2\theta ,$$

(38)

where the second and third harmonics are omitted. The substitution of ansatz (37) and approximation (38) in the nonlinear term in Eq. (24) for $m=2$ produces the squared amplitude of the solution,

$${\left(v_0^{(m=l=2)}\right)}^2\approx \frac{8}{15}\left(U_0-6\right),$$

(39)

in the case of

$$U_0-U_0^{(m=l=2)}\equiv U_0-6\ll 1.$$

(40)

2.3. Identifying the States with $m<l$

For the case of $l=1$ (the p state, in terms of atomic orbitals) and $m=0$, Eq. (24) produces a zonal mode, which is shaped as a dipole along the z direction. In the limit case of the infinitesimal amplitude, an exact solution of this type, suggested by the form of the ${\mathrm{Y}}_{10}$ spherical harmonic,

$$v^{(l=1,m=0)}\left(\theta \right)=v_0^{(l=1,m=0)}cos\theta ,$$

(41)

exists at the same threshold value, $U_0^{(l=1,m=0)}=2$, which is given by Eq. (30) for $l=1$. Above the threshold, viz., at small $U_0-2$, the squared amplitude of the mode (41) is

$${\left(v_0^{(l=1,m=0)}\right)}^2\approx \frac{2}{3}\left(U_0-2\right).$$

(42)

It coincides with the one given Eq. (36) for $l=m=1$, which also represents the p-wave orbital.

Further, in the case of $l=2$ (the d-wave state, in terms of atomic orbitals) and $m=1$, an appropriate ansatz for the respective dipole-vortex state (also known as a tesseral mode) is suggested by the form of spherical harmonic ${\mathrm{Y}}_{21}$:

$$v^{(l=2,m=1)}\left(\theta \right)=v_0^{(l=2,m=1)}sin\left(2\theta \right).$$

(43)

For $l=2$, the threshold value of the strength of the pull to the center is $U_0^{(l=2,m=0)}=6$, as given, once again, by Eq. (30) with $l=2$. Above the threshold, i.e., at $0<U_0-6\ll 1$ (cf. Eq. (40)), the squared amplitude of solution (43) is predicted as

$${\left(v_0^{(l=2,m=1)}\right)}^2\approx \frac{2}{3}\left(U_0-6\right),$$

(44)

cf. Eq. (39).

For the set of quantum numbers $l=2$ and $m=0$, the “sandwich-shaped" ansatz (with structure $\left(-+-\right)$ along the z axis, alias the zonal mode)), as suggested by the corresponding spherical harmonic ${\mathrm{Y}}_{20}$, is

$$v^{(l=2,m=0)}\left(\theta \right)=v_0^{(l=2,m=0)}\left(1-3{cos}^2\theta \right).$$

(45)

This expression with an infinitesimal amplitude gives an exact solution of Eq. (24) at the same threshold, $U_0=6$, which corresponds to $l=2$ and $m=2$ or $m=1$, see above. Finally, above the threshold, i.e., at $0<U_0-6\ll 1$, the approximation similar to that elaborated above predicts the squared amplitude of solution (45) as

$${\left(v_0^{(l=2,m=0)}\right)}^2\approx \frac{4}{11}\left(U_0-6\right),$$

(46)

cf. Eqs. (39) and (44).

Similarly, for $l>2$, the sectoral, tesseral, and zonal states, with $m=l$, $0<m<l$, and $m=0$, respectively, exist at $U_0>l(l+1)$. Above the existence threshold, amplitudes of these states scale as $\sqrt{U_0-l(l+1)}$, cf. Eqs. (36), and (42) for $l=1$, and Eqs. (39), (44), and (46) for $l=2$.

2.4. The Global Approximation for the States with $l\ge 1$

Similar to the interpolation approximation for the GS, represented by Eqs. (14) and (15), an approximation for the states carrying the angular momentum can be constructed by juxtaposing expressions (29), or (37), or (41), or (43), or (45) for $r\to 0$, and Eq. (25) for $r\to \infty$ in Eqs. (19) and (21):

$$\begin{array}{c}{\psi }_{\mathrm{interpol}}=exp\left(-i\mu t+im\phi \right)r^{-1}{V}^{\left(l,m\right)}\left(r,\theta \right),\\ V^{\left(l,m\right)}\left(r,\theta \right)\approx v^{\left(l,m\right)}\left(\theta \right)exp\left(-\sqrt{-2\mu }r\right)\end{array}$$

(47)

The family of the bound states is characterized by their norm, considered as a function of $\mu$, cf. Eq. (15). In terms of the interpolation approximation (47), the norm is written as

$$N_{\mathrm{interpol}}=\frac{I^{(l,m)}}{\sqrt{-2\mu }},I^{(l,m)}\equiv \pi {\int }_0^{\pi }{\left[v^{\left(l,m\right)}\left(\theta \right)\right]}^{2}sin\theta d\theta .$$

(48)

As well as its GS counterpart (15), this relation satisfies the anti-Vakhitov-Kolokolov criterion, which is a necessary stability condition for these bound states.

Comparison of this approximation with numerical results demonstrates that it provides reasonable accuracy for the vortex (sectoral) states with $m=l$, while for the tesseral and zonal states, with $m<l$, a discrepancy with the numerical solutions is too large. For the sectoral states, the substitution of the angular wave function (29) in Eq. (48) yields the following result:

$$N_{\mathrm{interpol}}^{(m=l)}\left(\mu \right)=\frac{2\pi \left(2m\right)!!}{(2m+1)!!\sqrt{-2\mu }}{\left(v_0^{(m=l)}\right)}^2.$$

(49)

In this expression, coefficients ${\left(v_0^{(m=l)}\right)}^2$ should be substituted as per Eqs. (36) and (39). Note that the scaling

$$N\sim 1/\sqrt{-\mu },$$

(50)

featured by Eq. 49,is an exact property of solutions generated by Eq. (20), irrespective of the applicability of any approximation.

2.5. The Thomas-Fermi (TF) Approximation

In the case of a large strength of the pull-to-the-center potential, $U_0\gg 1$, the angular equation (24) admits the application of the Thomas-Fermi (TF) approximation, which neglects the derivatives in that equation, yielding

$${\left(v_{\mathrm{TF}}^{\left(m\right)}\left(\theta \right)\right)}^2=\left\{\begin{array}{c}\frac{1}{2}\left(U_0-\frac{m^2}{{sin}^2\theta }\right),\mathrm{at}arcsin\left(\frac{m}{\sqrt{U_0}}\right)<\theta <\pi -arcsin\left(\frac{m}{\sqrt{U_0}}\right),\\ 0,\mathrm{at}0\le \theta <arcsin\left(\frac{m}{\sqrt{U_0}}\right)\mathrm{and}\pi -arcsin\left(\frac{m}{\sqrt{U_0}}\right)<\theta \le \pi .\end{array}\right.$$

(51)

The TF approximation is relevant in the case of $U_0>m^2$, when $arcsin\left(m/\sqrt{U_0}\right)$ exists. As any TF approximation, expression (51) is continuous, featuring discontinuities of the first derivative, $dv/d\theta$, at points $\theta =arcsin\left(m/\sqrt{U_0}\right)$ and $\theta =\pi -arcsin\left(m/\sqrt{U_0}\right)$. The largest squared value of the approximate solution (51) is

$${\left(v_{\mathrm{TF}}^{\left(m\right)}\left(\theta =\frac{\pi }{2}\right)\right)}^2=\frac{1}{2}\left(U_0-m^2\right).$$

(52)

For large $U_0$, the TF approximation applies to the full equation (22) as well: neglecting both the radial and angular derivatives in it, one obtains

$${\left(V_{\mathrm{TF}}^{\left(m\right)}(r,\theta )\right)}^2=\left\{\begin{array}{c}\frac{1}{2}\left(U_0-\frac{m^2}{{sin}^2\theta }\right)-\left|\mu \right|r^2,\mathrm{at}\frac{m^2}{2{sin}^2\theta }+\left|\mu \right|r^2<\frac{U_0}{2},\\ 0,\mathrm{at}\frac{m^2}{2{sin}^2\theta }+\left|\mu \right|r^2>\frac{U_0}{2}.\end{array}\right.$$

(53)

In this solution, it is taken into account that $\mu$ is negative for bound states. As well as the TF approximation (51), the result (53) is relevant for $U_0>m^2$.

Further, if the TF approximation is used at $r\to 0$, in the form of Eq. (51), but not globally, the corresponding interpolation-type approximation, similar to Eq. (47), is

$${\left(u_{\mathrm{interpom}}^{\left(m\right)}\left(r,\theta \right)\right)}_{\mathrm{TF}}=\frac{v_{\mathrm{TF}}^{\left(m\right)}\left(\theta \right)}{r}exp\left(-\sqrt{-2\mu }r\right).$$

(54)

3. Numerical Results

To produce stationary wave functions of the bound states in the numerical form, it is relevant, first, to substitute ansatz (19) into the underlying GPE (18), where $u^{\left(m\right)}$ is written as per Eq. (21), but with time-dependent $V^{\left(m\right)}\left(r,\theta ,t\right)$. The substitution leads to a time-dependent version of Eq. (22), viz.,

$$2i{r}^2\frac{\partial V^{\left(m\right)}}{\partial t}=-r^2\frac{{\partial }^2{V}^{\left(m\right)}}{\partial r^2}-\left[\frac{{\partial }^2}{\partial {\theta }^2}+cot\theta ·\frac{\partial }{\partial \theta }+\left(U_0-\frac{m^2}{{sin}^2\theta }\right)\right]V^{\left(m\right)}+2{\left|V^{\left(m\right)}\right|}^2{V}^{\left(m\right)}.$$

(55)

The use of the reduced wave function $V^{\left(m\right)}\left(r,\theta \right)$ makes it possible to present the results free from the singular factor $r^{-1}$ (see Eq. (19)).

Stationary solutions of Eq. (22) were found by means of the imaginary-time simulations [28] of Eq. (55), performed with the help of the split-step Fourier algorithm, using a set of $1024\times 128$ spatial modes in the $\left(r,\theta \right)$ domain. The above-mention scaling invariance of Eq. (20) suggest that it is sufficient to generate the results for a fixed norm, which we set here to be $N\equiv 1$, while the chemical potential $\mu$ is adjusted to keep this value of the norm, as per Eq. (50).

3.1. The Stationary Vortex (Sectoral) Modes with $m=l$

First, Figure 1 and Figure 2 display essential results for the vortex (sectoral) bound state with quantum numbers $m=l=1$, and compares them to the analytical predictions reported in the previous section. For the strength of the pulling potential $U_0=2.5$, the chemical potential of the bound state corresponding to the fixed norm, $N=1$, which is presented in Figure 1, is $\mu =-1.64$. In this case, the r-dependence of the numerically found value of $V^{(m=l=1)}\left(r,\theta =\pi /2\right)$, which is the maximum of wave function $V^{(m=l=1)}\left(r,\theta \right)$ for fixed r, is presented by the solid line in Figure 1(a). The dashed line in the same panel displays an exponential fit to the numerical result, viz.,

$$V^{(m=l=1)}\left(r,\theta =\pi /2\right)\approx 0.549exp\left(-\sqrt{-2\mu }r\right).$$

(56)

The respective interpolation approximation, produced by Eqs. (33), (36), and (47), is

$$V^{(m=l=1)}\left(r,\theta =\pi /2\right)\approx \sqrt{\frac{2}{3}\left(U_0-2\right)}exp\left(-\sqrt{-2\mu }r\right)\approx 0.577exp\left(-\sqrt{-2\mu }r\right),$$

(57)

being quite close to its numerically found counterpart. Further, the numerically obtained profile of $v^{(1,1)}\left(\theta \right)\equiv V^{(1,1)}\left(r=0,\theta \right)$ for the same $U_0=2.5$ is displayed by the solid curve in Figure 1(b), while the dashed one shows the corresponding analytical approximation suggested by Eq. (33), with the same coefficient as used for the fitting curve in panel Figure 1(a) (see Eq. (56)), viz., $v^{\left(1.1\right)}\left(\theta \right)\approx 0.549sin\theta$. It is seen that the fit profile is virtually identical to the numerical one. If the fit coefficient $0.549$, taken from Eq. (56), is replaced by its analytically predicted counterpart, $\sqrt{(2/3)\left(U_0-2\right)}\approx 0.577$ [see Eq. (57)], the agreement between the numerical findings and analytical predictions remains very accurate.

The family of the bound states with $m=l=1$ is characterized in Figure 2(a) by the numerically found dependence of the largest value of their wave function, $V^{(1,1)}\left(r=0,\theta =\pi /2\right)\equiv v^{(m=l=1)}(\theta =\pi /2)$, on the potential strength $U_0$, plotted along with the analytical prediction (36). Although the prediction was derived under the condition of $0<U_0-2\ll 1$, it is seen that, even for $U_0=4$, the analytical result remains a highly accurate one.

The adequacy of the TF approximation for large values of the potential’s strength $U_0$ is corroborated by Figure 2(b), which displays a numerically generated profile of $v^{\left(1\right)}\left(\theta \right)$ for $U_0=20$, along with its TF-predicted counterpart, given by Eq. (51). It is observed that the approximation is reasonable. In particular, the largest value of the wave function, given by expression (52), is practically identical to its numerically obtained counterpart.

For the higher-order vortex (sectoral) mode with $m=l=2$, an example of the numerically found profile, $V^{(2,2)}\left(r=0,\theta \right)\equiv v^{(2,2)}\left(\theta \right)$, is displayed in Figure 3(a) for $U_0=6.5$, along with the analytically predicted profile, $0.564{sin}^2\theta$ (see Eq. (39)), with the coefficient produced by the fit procedure. The chemical potential of this state is $\mu =-1.44$. The analytical prediction is reproduced by the numerical solution practically exactly.

The family of the vortex modes with $m=l=2$ is represented in Figure 3(b) by the dependence of $V^{(m=l=2)}\left(r=0,\theta =\pi /2\right)\equiv v^{(m=l=2)}(\theta =\pi /2)$ on the potential’s strength $U_0$, as found from the numerical solution and predicted by Eq. (39). In comparison with the similar dependences shown in Figure 2(a) for $m=l=1$, the discrepancy between the numerical and analytical results is somewhat larger, but the analytical approximation is still quite reasonable even for $U_0=8$, when the formal applicability condition (40) does not hold.

3.2. The Bound States with $m<l$

An example of the numerically found angular profile of the dipole (zonal) mode, with $l=1$, $m=0$ and $U_0=4$, is plotted by the solid curve in Fig. Figure 4(a). The dashed line is its analytically predicted counterpart (41), with the fit coefficient $v_0^{(1,0)}=1.209$, while the respective analytical expression (42) yields a close value, $v_0^{(1,0)}\approx 1.155$. The chemical potential of this bound state is $\mu =-36$. The radial dependence $V(r,\theta =0)$ is shown in Figure 4(b) (the solid line). The family of the dipole-shaped (zonal) states is characterized in Figure 4(c) by dependences of $V^{(l=1,m=0)}\left(r=0,\theta =\pi /2\right)\equiv v^{l=1,m=0)}(\theta =\pi /2)$ on strength $U_0$ of the pulling potential, as found from the numerical solution, and as predicted analytically by Eq. (39). It is seen from these figures that the accuracy of the analytical approximation is reasonable, even for values of $U_0-2$ which are not small.

Further, the solid and dashed curves in Figure 5(a) show, respectively, an example of the numerically obtained profile of the vortex-dipole (tesseral) mode, with $l=2$, $m=1$ and $U_0=8$, and its analytical counterpart, produced by Eq. (43) with the fit coefficient $v_0^{(2,1)}=1.13$ (the chemical potential of this state is $\mu =-51$). Note that Eq. (44) gives the analytical approximation for the latter coefficient as $v_0^{(2,1)}=1.155$, which is very close to the fit value which makes the numerical and analytical profiles virtually identical in Figure 5(a). Figure 5(b) shows the numerically found radial dependence $V(r,\theta =\pi /4)$ at $U_0=8$. The family of the tesseral states is characterized in Figure 5(c) by dependences of $V^{(l=2,m=1)}\left(r=0,\theta =\pi /2\right)\equiv v^{(l=2,m=1)}(\theta =\pi /2)$ on the potential’s strength $U_0$, as given by numerical solution, and as predicted by Eq. (44). It is seen from these figures that the accuracy of the analytical approximation is reasonable, even for values of $U_0-2$ which are not small.

Finally, numerical results for the “sandwich-shaped" (zonal) state with $l=2$, $m=0$ are presented in Figure 6. An example of the respective profile of $v^{(l=2,m=0)}\left(\theta \right)$ is plotted in Figure 6(a) for $U_0=8$ and chemical potential $\mu =-30.97$, along with its analytically-predicted counterpart (45), with coefficient $v_0^{(l=2,m=0)}=0.852$, which is predicted by Eq. (46). Figure 6(b) shows the numerically found radial dependence $V(r,\theta =\pi /2)$ at $U_0=8$. The family of the zonal states is characterized by the dependence of $V^{(l=2,m=0)}\left(r=0,\theta =\pi /2\right)\equiv v^{(l=2,m=0)}(\theta =\pi /2)$ on the potential’s strength $U_0$, which is plotted in Figure 6(c), along with the respective analytical prediction given by Eq. (46). In this case too, the accuracy of the analytical approximation is reasonable.

4. Conclusion

A solution for the problem of the quantum collapse in the gas of particles pulled to the center by the well-known potential, $U\left(r\right)=-U_0/\left(2r^2\right)$, was reported in Ref. [7]. Repulsive collisions between the particles, represented by the cubic term in the respective GPE (Gross-Pitaevskii equation), suppress the collapse and create the otherwise missing spherically symmetric GS (ground state), alias the s-wave orbital, in terms of atomic physics. However, the previous analysis did not address bound states with reduced symmetry, which carry angular momentum. The present work reports results of the systematic analysis of such states, which are labeled, as usual, by the orbital and magnetic quantum numbers, l and m. These states exist when the strength of the potential exceeds a threshold value, $U_0>l(l+1)$, which does not depend on m. At a relatively small distance from the threshold, the sectoral, tesseral, and zonal states, with $m=l$, $0<m<l$, and $m=0$, respectively, are found in the approximate analytical form for $l=1$ and 2 (alias the p- and d-wave orbitals). In the general case, these states are found in the numerical form, which corroborates the accuracy of the analytical approximation.

As an extension of this work, it may also be interesting to analyze dynamical excitation of the states with angular momentum from the GS by pulses of an incident vortex field, and to consider Rabi oscillations between different bound states driven by an ac field. Also relevant is a possibility to consider the p- and d-wave orbitals in the framework of the many-body theory, instead of the mean-field approximation, following the lines of Ref. [13].