Angular-momentum modes in a bosonic condensate trapped in the inverse-square potential

In the mean-field approximation, the well-known effect of the critical quantum collapse in a 3D gas of particles pulled to the center by potential U(r) = -U_0/r^2 is suppressed by repulsive interparticle interactions, which create the otherwise non-existing s-wave ground state. Here, we address excited bound states carrying angular momentum, with the orbital and magnetic quantum numbers, l and m. They exist above a threshold value of the potential's strength, U_0>l(l+1). The sectoral, tesseral, and zonal modes, which correspond to m = l, 0<m<l, and m = 0, respectively, are found in an approximate analytical form for relatively small values of U_0 - l(l+1). Explicit results are produced for the p- and d-wave states, with l = 1 and 2, respectively. In the general form, the bound states are obtained numerically, confirming the accuracy of the analytical approximation.


I. INTRODUCTION AND THE MODEL
It is well known that attractive potential U (r) = −U 0 / 2r 2 , 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, 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, while at U 0 > 1/4 Eq. (1) does not maintain the GS, leading, instead, to the quantum collapse of wave function ψ (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 charge 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 ∼ 100 proton masses and the elementary electric charge set at the center, the critical value (2) corresponds to a very small dipole moment, d cr ∼ 10 −6 Debye.Therefore, the overcritical case of U 0 > 1/4 is relevant in the actual physical context.
On the other hand, recently new techniques have been elaborated for the creation of efficient optical trapping potentials for ultracold atoms, the commonly known ones being red-and blue-shifted optical lattices [8][9][10].More sophisticated potentials can be "painted" in space by fast moving laser beams [11].Recently, a technique based on the use of digital micromirror devices was elaborated for the creation of various traps for 2D Bose-Einstein condensates (BECs) [12].In particular, it may be used for inducing anharmonic potentials ∼ r k , with k ≥ 2, that were explored as a setting for the control of 2D vortex dynamics in BEC [13].In this connection, it is relevant to mention that the pull-to-the-center potential ∼ −1/r 2 for an atom in the 3D space can be induced by red-shifted laser illumination provided by beams isotropically converging towards r = 0, resembling that igniting the inertial nuclear fusion (with the power larger by many orders of magnitude) [14,15].Indeed, the energy of the detuned dipole interaction of the atom with the optical field of local intensity I(r) is where c 0 is the light speed in vacuum, ω 0 the frequency of the atomic optical transition, ∆ its detuning from the driving frequency, and Γ the linewidth [16].Then, in the geometric-optics approximation, the substitution of where P is the integral illumination power, in Eq. ( 3), produces the potential sought for.Then, an estimate for ω, corresponding for the wavelength of the laser beams λ ∼ 1 µm and the lightest atom (hydrogen), demonstrates that the critical value (2) corresponds to extremely small powers, P cr ∼ 10 −7 W, hence the overcritical case is fully relevant for the experimental implementation.
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 of 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 BEC.In particular, ultracold gases of polar molecules such as LiCs [17] and KRb [18] 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 [19].
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 [20], using the Gross-Pitaevskii equation (GPE) with the same potential as in Eq. ( 1), In the properly scaled form, the GPE id written as where the cubic term in Eq. ( 5) represents the self-repulsion in the gas, in the framework of the MF approximation.
The coefficient in front of this term is set equal to 1 by means of normalization.Dipole-dipole interactions between the particles in the gas were also taken into account, applying another version of the MF approximation.It considers the interaction of a given dipole with the collective electric potential created by all other dipoles, and amounts to renormalization of coefficient U 0 .
In addition to the irreducible coefficient U 0 , another control parameter of the MF theory is the norm, Alternatively, it is possible to fix that the norm is 1, while the cubic term in Eq. ( 5) will acquire coefficient N as a free parameter.For a mixture of two species of particles, a two-component extension of Eq. ( 5) was introduced in Ref. [21].The respective system of nonlinearly coupled GPEs also provides the suppression of the quantum collapse and creation of the GS [21].
Solutions of Eq. ( 5) for stationary states with chemical potential µ < 0, in the form of with steady-state wave function u satisfying the equation These states are characterized by their norm (6) and angular momentum, where * stands for the complex conjugate.N and M are dynamical invariants of GPE (5).It also conserves the Hamiltonian, The linearized version of Eq. ( 5) produces exact, although unnormalizable (with diverging norm (6)), pairs of eigenmodes with µ = 0 and orbital and magnetic quantum numbers l, m [7], which are written in spherical coordinates (r, θ, φ).Here u 0 is an arbitrary amplitude, Y lm are the spherical harmonics, and Solution (11) implies that these eigenstates exists at U l < 1/4, cf.Eq. ( 2), while the quantum collapse occurs at The suppression of the collapse and recreation of the GS (with a normalizable wave function, unlike the unnormalizable one (11)) in the framework of the full nonlinear GPE (5) was demonstrated for the spherically isotropic state, corresponding to l = m = 0 in expression (11) (i.e., the s-wave orbital, in terms of atomic physics [1]).To this end, Eq. ( 8) is converted, by substitution with real function v(r), into a radial equation An asymptotic solution to Eq. ( 14) at r → 0 is (here, const is not determined by the asymptotic expansion).At r → ∞, the bound-state wave function with µ < 0 decays exponentially, v ∼ exp − √ −2µr .A global approximation for the GS can be constructed as an interpolation, juxtaposing the asymptotic forms which are valid at r → 0 and r → ∞: Solutions with µ > 0 correspond to delocalized states with the following asymptotic form at r → ∞: These states are not considered here in detail.Singularity ∼ r −1 of wave function (16) at r → 0 is acceptable, as the corresponding 3D norm integral (6) converges at r → 0. In particular, the approximate wave function (16) gives rise to the following relation between the chemical potential and norm, Actually, scaling µ ∼ −1/N 2 , which is produced by Eq. ( 18), is an exact property of solutions to Eq. ( 8), irrespective of validity of the approximation.In the limit of µ → −0, Eq. ( 16) produces an exact solution of Eq. ( 5), ψ µ=0 (r) = U 0 /2r −1 , with a divergent norm.
The numerical solution of Eq. ( 8) corroborates the conclusion following from Eqs. ( 15) -( 18): GPE ( 5) 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 ∼ 100 proton masses and electric dipole moment ∼ 1 Debye, the radius of the newly predicted GS, which replaces the collapsing regime in the gas of ∼ 10 5 particles, is r GS ≃ 3 µm [7].Note that the MF approximation is relevant for states characterized by such a length scale.
As concerns the implementation of the model with the pull-to-the center optical potential provided by Eqs. ( 3) and ( 4), an estimate demonstrates that the same order of magnitude of the GS radius is expected for the gas of ∼ 10 6 atoms under the action of the isotropically converging illumination with power P ∼ 10 −3 W.
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) [22].For the number of particles ≳ 100, the barrier is practically impenetrable, hence the MF-predicted solution becomes an effective GS.
Relation (18) satisfies the anti-Vakhitov-Kolokolov criterion, dµ/dN > 0, which is a necessary condition for stability of families of bound states supported by a self-repulsive nonlinearity [23] (the Vakhitov-Kolokolov criterion per se, dµ/dN < 0, pertains to the case of the self-attraction [24,25]).The stability of the GS created by Eq. ( 5) 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 [26].In that case, the spherically symmetric pulling potential in Eq. ( 5) is replaced by an axially symmetric one, − U 0 /2r 2 cos θ.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, (U 0 ) (m) cr = 1.28, 7.58, and 19.06.For comparison, Eqs. ( 2) and ( 36) yield the set of smaller critical values, viz.(U 0 ) (l) cr = 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 > (U 0 ) (m) cr .Vortex states were also addressed in the 2D version of the setting, whose linearized form gives rise to the collapse at all positive values of U 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 [27].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 [28,29].The resulting amended GPE predicts self-trapped states in the form of quantum droplets, which have been observed experimentally with quasi-2D [30,31] and full 3D [32,33] shapes in binary BECs with the contact interactions, as well as in dipolar BECs with long-range interactions between atomic magnetic dipole moments [34,35].
In the present context, the corresponding 2D GPE takes the form of Eq. ( 5), with term |ψ| 2 ψ replaced by |ψ| 3 ψ [36], leading to the solution for the wave function with the following asymptotic form at r → 0, which replaces the collapsing solution of the 2D linear equation: (cf. expression ( 16) at r → 0 and Eq. ( 12)), where (r, φ) 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 → 0. A nontrivial problem is stability of the 2D vortex bound states corresponding to the asymptotic form (19), which exist under the condition of U 0 > m 2 − 4/9 ≡ (U 0 ) min .Surprisingly, this problem admits an exact solution: the vortex states with m ≥ 1 are stable provided that U 0 > (7/9) 3m 2 − 1 ≡ (U 0 ) stability , the GS with m = 0 being stable too [36].
Note that, for all m ≥ 1, (U 0 ) stability exceeds (U 0 ) min , hence the 2D vortex states are unstable in the interval of m 2 − 4/9 < U 0 < (7/9) 3m 2 − 1 .The dominant instability mode is one which leads to slow drift of the vortex' pivot off the central point, along a spiral trajectory [36].
The objective of the present work is to construct 3D states, carrying angular momentum, as solutions of GPE (5 with nonzero orbital quantum number, l ≥ 1, and the magnetic quantum number taking values 0 ≤ m ≤ 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 ≥ 1, viz., (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.

II. ANALYTICAL CONSIDERATIONS
A. The framework for the analysis of the Gross-Pitaevskii equation It is relevant to explicitly write the underlying 3D equation ( 5) in the spherical coordinates: Stationary solutions to Eq. ( 21), with chemical potential µ < 0 and vorticity which is represented by the magnetic quantum number m, are looked for as cf. Eq. ( 7), with real function u (m) (r, θ) satisfying the stationary equation, To eliminate the singularity, ∼ 1/r, of the solution, we perform substitution (13), i.e., in Eq. (23), which leads to an equation for V (m) (r, θ) [cf.Eq. ( 14)], For the consideration of the form of eigenmodes in the limit of r → 0, one may define, straightforwardly, and directly set r = 0 in Eq. ( 25), which leads, without any approximation, to the ordinary differential equation with respect to the angular coordinate: Equation ( 27) should be solved in the natural region of 0 ≤ θ ≤ π.In fact, Eq. ( 26) implies that relevant solutions to Eq. ( 27) provide the boundary condition for solutions to Eq. ( 25) at r → 0. The boundary condition at r → ∞ amounts to the exponential decay of the wave function: (recall the bound states may exist solely with µ < 0).The presence of the artificial singularity in Eq. ( 27) at θ = 0 and π suggests a possibility of the existence of singular solutions with the asymptotic form The substitution of this ansatz in Eq. ( 27) yields Obviously, it follows from Eq. ( 30) that the singular solution may exist only for m = 0 (in particular, for m = 0 and U 0 = 0 Eq. ( 27) has an exact singular solution, v (0) (U 0 = 0) = √ 2 sin θ −1 ).Actually, the singular solution ( 29) is irrelevant, as the respective normalization integral includes a divergent factor, B. Identifying vortex (sectoral) states with equal magnetic and orbital quantum numbers, m = l Vortex-state solutions of Eq. ( 25) with m = l, also known as sectoral ones, are determined by Eq. ( 27) with boundary conditions at θ = 0 and π (by definition, m does not take negative integer values).Equation ( 27) with these boundary conditions admits an exact solution, suggested by the form of spherical harmonic Y ll , in the form of with an infinitesimal amplitude v (m) 0 , at the threshold value of the strength of the pulling potential, which is an example of the general threshold value (20).Above the threshold, i.e., at (m) 0 in Eq. ( 32) 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 condition (34) 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 (32), as an approximation, in the nonlinear term of Eq. ( 27), and make use of identity Neglecting, as usual, the third harmonic in Eq. ( 37) and keeping the term ∼ sin θ, for small values of the distance from the threshold, viz., the so approximated Eq. ( 27) predicts the squared amplitude of the approximate solution (32) as Similarly, for the vortex (sectoral) state with m = l = 2 (the d -wave orbital, in terms of atomic physics), the exact solution (32) with an infinitesimal amplitude v exists at U 0 = 6, as given by Eq. ( 33) with l = 2.To predict the amplitude of this state above the threshold, one can use the following approximation: sin 6 omitting the second and third harmonics.The substitution of ansatz (40) and approximation (41) in the nonlinear term of Eq. ( 27) with m = 2 produces the squared amplitude of the solution, in the case of C. Identifying the states with m < l For the case of l = 1 (the p state, in terms of atomic orbitals) and m = 0, Eq. ( 27) 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 Y 10 spherical harmonic, exists at the same threshold value, U (l=1,m=0) 0 = 2, which is given by Eq. ( 33) for l = 1.Above the threshold, viz., at small U 0 − 2, the squared amplitude of the mode (44) is It coincides with the one given Eq. ( 39) 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 Y 21 : For l = 2, the respective threshold value of the strength of the pull to the center is U (l=2,m=0) 0 = 6, as given, once again, by Eq. ( 33) with l = 2. Above the threshold, i.e., at 0 < U 0 − 6 ≪ 1 (cf.Eq. ( 43)), the squared amplitude of solution ( 46) is predicted as cf. Eq. (42).
D. The global approximation for the states with l ≥ 1 Similar to the interpolation approximation for the GS, represented by Eqs. ( 16) and ( 18), an approximation for the states carrying the angular momentum can be constructed by juxtaposing expressions (32), or (40), or (44), or (46), or (48) for r → 0, and Eq. ( 28) for r → ∞ in Eqs. ( 22) and ( 24): The family of the bound states is characterized by their norm, considered as a function of µ, cf.Eq. ( 18).In terms of the interpolation approximation (50), the norm is written as As well as its GS counterpart (18), 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 ones, with m < l, a discrepancy with the numerical solutions is relatively large.For the sectoral states, the substitution of the angular wave function (32) in Eq. ( 51) yields the following result: In this expression, coefficients v (m=l) 0 2 should be substituted as per Eqs.( 39) and ( 42).Note that the scaling featured by Eq. 52,is an exact property of solutions generated by Eq. ( 23), irrespective of the applicability of any approximation.

E. The Thomas-Fermi (TF) approximation
In the case of a large strength of the pull-to-the-center potential, U 0 ≫ 1, the angular equation ( 27) admits the application of the Thomas-Fermi (TF) approximation, which neglects the derivatives in that equation, yielding The TF approximation is relevant in the case of U 0 > m 2 , when arcsin m/ √ U 0 exists.As any TF approximation, expression (54) is continuous, featuring discontinuities of the first derivative, dv/dθ, at points θ = arcsin m/ √ U 0 and θ = π − arcsin m/ √ U 0 .The largest squared value of the approximate solution (54) is For large U 0 , the TF approximation applies to the full equation ( 25) as well: neglecting both the radial and angular derivatives in it, one obtains In this solution, it is again taken into account that µ is negative for bound states.As well as the TF approximation (54), the result (56) is relevant for U 0 > m 2 .Further, if the TF approximation is used at r → 0, in the form of Eq. ( 54), but not globally, the corresponding interpolation-type approximation, similar to Eq. ( 50), is  22) and ( 24)) at θ = π/2, as obtained from the numerical solution of Eq. ( 25) (the solid line), and as produced by the exponential fit to it, see Eq. ( 59) (the dashed line).(b) The solid line: the numerically produced profile of v (1,1) (θ) ≡ V (1,1) (r = 0, θ) (see Eq. ( 26)), for U0 = 2.5.The chemical potential of this bound state is µ = −1.64.The dashed line shows the analytical prediction for the same profile, produced by Eq. ( 36), with coefficient v (m=l=1) 0 replaced by the same fit value, 0.549, as in panel (a) (see Eq. ( 59)).

B. 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. 4(a).The dashed line is its analytically predicted counterpart (44), with the fit coefficient v  4(c) by dependences of V (l=1,m=0) (r = 0, θ = π/2) ≡ v l=1,m=0) (θ = π/2) on strength U 0 of the pulling potential, as found from the numerical solution, and as predicted analytically by Eq. (42).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.
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.
In physical units, characteristics of the predicted stationary states are similar to the above-mentioned ones for the GS.In particular, in the case of the gas containing ∼ 10 5 dipolar particles with the dipole moment ∼ 1 Debye, pulled to the central electric charge, the eigenstates carrying the angular momentum will have the size of several microns.Essentially the same size of the sectoral, tesseral, and zonal eigenstates is expected in the gas of ∼ 10 6 light atoms pulled to the center by the converging optical potential defined as per Eqs.( 3) and (4) with illumination power ∼ 10 −3 W, The next objective of the analysis is to accurately investigate stability of the stationary modes, found above in the numerical and approximate analytical forms.A direct approach to the stability analysis may be based on real-time simulations of Eq. (58).Technically, this is a challenging problem due to the presence of the artificial singularities in Eq. ( 58) at θ = 0 and π.Preliminary simulations demonstrate that, in some cases, apparent artifacts appear in the form of growing singularities near θ = 0 and π, resembling the above-mentioned irrelevant stationary singular solutions given by Eqs. ( 29) and (30).Accurate dynamical simulations is a subject for a separate work, which will be reported elsewhere.In any case, it is relevant to mention that the well-known splitting instability of vortex states in nonlinear systems with self-attractive nonlinearity [38,39], which is driven by azimuthal perturbations, has no reason to exist in the present case of the self-repulsion.
As a further 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. [22].