Singular mean-field states: A brief review of recent results

This article provides a focused review of recent findings which demonstrate, in some cases quite counter-intuitively, the existence of bound states with a singularity of the density pattern at the center, while the states are physically meaningful because their total norm converges. One model of this type is based on the 2D Gross-Pitaevskii equation (GPE) which combines the attractive potential ~ 1/r^2 and the quartic self-repulsive nonlinearity, induced by the Lee-Huang-Yang effect (quantum fluctuations around the mean-field state). The GPE demonstrates suppression of the 2D quantum collapse, driven by the attractive potential, and emergence of a stable ground state (GS), whose density features an integrable singularity ~1/r^{4/3} at r -->0. Modes with embedded angular momentum exist too, and they have their stability regions. A counter-intuitive peculiarity of the model is that the GS exists even if the sign of the potential is reversed from attraction to repulsion, provided that its strength is small enough. This peculiarity finds a relevant explanation. The other model outlined in the review includes 1D, 2D, and 3D GPEs, with the septimal (seventh-order), quintic, and cubic self-repulsive terms, respectively. These equations give rise to stable singular solitons, which represent the GS for each dimension D, with the density singularity ~1/r^{2/(4-D). Such states may be considered as a result of screening of a"bare"delta-functional attractive potential by the respective nonlinearity.


I. Introduction
A. Singular states pulled to the center by attractive fields 1. Outline of the topic A commonly adopted condition which distinguishes physically relevant states produced by various two-and threedimensional (2D and 3D) models originating in quantum mechanics, studies of Bose-Einstein condensates (BECs), nonlinear optics, plasma physics, etc., is that the respective fields, such as wave functions in quantum mechanics and mean-field (MF) description of BEC, or local amplitudes of optical fields, must avoid singularities at r → 0, where r is the radial coordinate. A usual example of the relevance of this condition is provided by the quantum-mechanical Schrödinger equation with the attractive Coulomb potential, U(r) ∼ −r −1 : the stationary real wave function of states arXiv:2007.02320v1 [cond-mat.quant-gas] 29 Jun 2020 with integer azimuthal quantum number (alias vorticity) l ≥ 0 has expansion u(r) = u 0 + u 1 r + O r 2 r l (1) at r → 0, thus avoiding a singularity, despite the fact that the trapping potential is singular [1]. In quantum mechanics, a critical role is played by a more singular attractive potential, viz., with U 0 > 0, which gives rise to the quantum collapse, alias "fall onto the center " [1]. This well-known phenomenon means nonexistence of the ground state (GS) in the 3D and 2D Schrödinger equations with potential (2). In 3D, the collapse occurs when U 0 exceeds a finite critical value, U (3D) 0 cr (in the notation adopted below in Eq. (8), = 0, i.e., the 2D collapse happens at any U 0 > 0.
In both 3D and 2D cases, potential (2) may be realized as electrostatic pull of a particle (small molecule), carrying a permanent electric dipole moment, to a charge placed at the center, assuming that the local orientation of the dipole is fixed by the minimization of its energy in the central field [2]. In addition to that, in the 2D case the same potential (47) may be realized as attraction of a magnetically polarizable atom to a thread carrying electric current (e.g., an electron beam) transversely to the system's plane, or the attraction of an electrically polarizable atom to a uniformly charged transverse thread (other 2D settings in Bose-Einstein condensates (BECs) under the action of similar fields were considered in Refs. [3] and [4]).
A fundamental issue is stabilization of the 3D and 2D quantum-mechanical settings with pulling potential (2) against the collapse, with the aim to create a missing GS. A solution was proposed in Refs. [5]- [7], which replaced the original quantum-mechanical problem by one based on a linear quantum-field theory. While such a model produces a GS, it does not answer a natural question: what an effective radius of the GS is for given parameters of the setting, such as U 0 in (2) and the mass of the quantum particle, m. Actually, the field-theory solution defines the GS size as an arbitrary spatial scale, which varies as a parameter of the respective field-theory renormalization group. A related problem is the definition of self-adjoint Hamiltonians in the linear quantum theory including the interaction of a particle with singular potentials [8,9].
Another solution was proposed in Ref. [2], which replaced the 3D linear Schrödinger equation by a nonlinear Gross-Pitaevskii equation (GPE) [10] for a gas of particles attracted to the center by potential (47), with repulsive inter-particle collisions. In the framework of the mean-field (MF) approximation, the 3D GPE gives rise to the missing GS at all values of U 0 > U (3D) 0 cr . The radius of the GS is fully determined by model's constants, i.e., U 0 , m, the scattering length of the inter-particle collisions, and the number of particles, N . For typical values of the physical parameters, an estimate for the radius is a few microns. Beyond the framework of the MF, it was demonstrated that the many-body quantum theory, applied to the same setting, does not, strictly speaking, create GS, but the interplay of the attraction to the center and inter-particle repulsion gives rise to a metastable state. For sufficiently large N , the metastable state is nearly tantamount to GS, being separated from the collapse regime by a very tall potential barrier [11]. Further, the mean-field GS was also constructed in the 3D gas embedded in a strong uniform electric field, which reduces the symmetry of the effective pulling potential from spherical to cylindrical [12], as well as in the two-component 3D gas [13].
In 2D, the problem is more difficult, as the repulsive cubic term in GPE, which represents inter-atomic collisions in the MF approximation [10], is not strong enough to suppress the 2D quantum collapse and create the GS. The main issue is that, in 3D and 2D settings alike, the MF wave function, ψ(r), produced by the respective GPE, has density |ψ(r)| 2 diverging ∼ r −2 at r → 0. In terms of the integral norm, (D = 3 or 2 is the dimension), the density singularity r −2 is integrable in 3D, but not in 2D, where it gives rise to a logarithmic divergence of the norm: where r cutoff is a cutoff (smallest) radius, which may be determined by the size of particles in the condensate. Actually, the cubic self-repulsion is critical in 2D, as any stronger nonlinear term is sufficient to stabilize the 2D setting. In 3D, the critical value of the repulsive-nonlinearity power, which also leads to the logarithmic divergence of N , is 7/3; it is relevant to mention that the respective nonlinear term, |ψ| 4/3 ψ, represents the effective repulsion in the density-functional model of the Fermi gas [14]- [17].
A solution for the 2D setting is offered by the quintic defocusing nonlinearity [2], that may account for three-body repulsive collisions in the bosonic gas [18,19]. However, a difficulty in the physical realization of the quintic term is the fact that three-body collisions give rise to effective losses, kicking out particles from BEC to the thermal component of the gas [20]- [22].

New results included in the review
Recently, much interest was drawn to quasi-2D and 3D self-trapped states in BEC in the form of "quantum droplets", filled by a nearly incompressible two-component condensate, which is considered as an ultradilute quantum fluid. This possibility was predicted in the framework of the 3D [24] and 2D [25][26][27] GPEs which include the Lee-Huang-Yang (LHY) corrections to the MF approximation [28]. They represent effects of quantum fluctuations around the MF states. The two-component structure of the condensate makes it possible to provide nearly complete cancellation between the inter-component MF attraction and intra-component repulsion (which, in turn, can be adjusted by means of the Feshbach-resonance technique [29]) and thus create stable droplets through the balance of the relatively weak residual MF attraction and LHY-induced quartic self-repulsion. The quantum droplets of oblate (quasi-2D) [30,31] and fully 3D (isotropic) [32,33] shapes were created in a mixture of two different spin states of 39 K atoms, as well as in a mixture of 41 K and 87 Rb atoms [34]. Further, it was predicted that 2D [35][36][37] and 3D [38] droplets with embedded vorticity have their stability regions too. The LHY effect has also opened the way to the creation of stable 3D droplets in single-component BEC with long-range interactions between atoms carrying magnetic dipole moments [39]- [43], although dipolar-condensate droplets with embedded vorticity are unstable [44].
The LHY effect in the 3D model of BEC pulled to the center by potential (2) was considered in Ref. [23], with a conclusion that the LHY term gives rise to a quantum phase transition at U 0 = 2/9 (note that it is close to but smaller than the above-mentioned critical value for the linear Schrödinger equation, U (3D) 0 cr = 1/4). The phase transition manifests itself in the change of the asymptotic form of the GS stationary wave function at r → 0: This phase transition may be categorized as one of the first kind, as power α in terms r −α in Eq. (5) undergoes a finite jump at U 0 = 2/9, from 1/3 to 2/3. In Section 2 of the present review we summarize recent results which demonstrate that the stabilization of the GS in the quasi-2D bosonic gas pulled to the center by potential (47) may be provided by the LHY correction to the GPE [45]. This possibility is essential because, as mentioned above, the alternative, in the form of the quintic self-repulsion, is problematic in the BEC setting. The underlying three-dimensional GPE, which includes the LHY quartic defocusing term, is where Ψ stands for equal wave functions of two components of the BEC, W (r) is the general trapping potential, a > 0 is the scattering length of inter-particle collisions, δa ≷ 0, with |δa| a, represents the above-mentioned small disbalance of the inter-component attraction and intra-component repulsion, and the last term in Eq. (6) is the LHY correction to the MF equation [24].
The reduction of Eq. (6) to the 2D form, with coordinates (x, y), under the action of tight confinement applied in the z direction, was derived in Ref. [25], producing the GPE with a cubic term multiplied by an additional logarithmic factor, However, this limit implies extremely strong confinement in the z direction, with the transverse size a ⊥ ξ, where the healing length is ξ = 32 √ 2/3π (a/ |δa|) 3/2 a ≈ 5 (a/ |δa|) 3/2 a [24]. For experimentally relevant parameters [30]- [33], an estimate is ξ 30 nm. On the other hand, a realistic size of the confinement length in the experiment is few µm, implying relation a ⊥ ξ, opposite to the above-mentioned one necessary for the derivation of Eq. (7). Therefore, it is relevant to reduce Eq. (6) to the 2D form, keeping the same nonlinearity as in Eq. (6).
Results for both GS and vortex states, produced by the analysis of the 2D equation (9) [45], are summarized, as a part of the present review, in Section 2. Essential conclusions are that all the GS solutions (with zero vorticity) are stable, while vortex modes are stable if U 0 exceeds a certain critical value, which depends on the vorticity.
B. Singular solitons: previously known results

Singular solitons in free space
The condition of the convergence of the integral norm is equally relevant for self-trapped states, which are predicted, as localized solutions, by models such as the nonlinear Schrödinger equation (NLSE) for a complex wave function, ψ (x, y, z, t): where σ = 1 and −1 correspond, respectively, to the self-repulsive (defocusing) and attractive (focusing) signs of the nonlinearity. Typical physical realizations of the NLSE feature cubic (ν = 1) and quintic (ν = 2) nonlinearities. The commonly known bright-soliton solutions of the self-focusing (σ = −1) NLSE in 1D, with arbitrary real chemical potential µ < 0, are free of singularities: At ν ≥ 2, these 1D solutions are subject to instability against the wave collapse, i.e., catastrophic shrinkage of the self-focusing field [49]. In the case of σ = +1 (self-defocusing), Eq. (10 gives rise to an exact solution in the form of a singular soliton, For ν ≤ 2, this singular solution is physically irrelevant because it gives rise to a divergent total norm, N = +∞ −∞ |ψ(x)| 2 dx. In particular, in terms of the usual GPE, the divergence of the norm implies that the creation of BEC states in the form of singular solitons (12) would require an infinite number of atoms. In optics, where the 1D version of Eq. (10), with t replaced by the propagation distance, z, governs paraxial propagation of a stationary light beam in a planar nonlinear waveguide with transverse coordinate x, the divergence implies an infinite power of the beam. On the other hand, at ν > 2 the norm of the singular solution (12) converges, the exact result being where Γ is the Euler's Gamma-function. Note that, in the limit of µ → −0, solution (12) takes the form of whose norm diverges at all values of ν, in agreement with Eq. (13).
The singularity produced by Eq. (15) is physically admissible, producing a convergent norm, for ν > 1 at D = 2, and for ν > 2/3 at D = 3. Thus, in the 3D case only the interval of is a physically relevant one for NLSE (10) with the self-defocusing nonlinearity. If the norm converges, Eq. (10) implies that it is related to the chemical potential by the following scaling formula: It is worthy to note that, at ν > 2/D, relation (18) satisfies the anti-Vakhitov-Kolokolov (anti-VK) criterion, dN/dµ > 0, which is a necessary stability condition for a family of bound states maintained by any self-defocusing nonlinearity [54] (the VK criterion proper, dN/dµ < 0, is necessary for stability of soliton families created by selffocusing nonlinearity [49,55]). In particular, the entire 3D existence interval (17) is compatible with the anti-VK condition, which suggests, in particular, that 3D solutions generated by the cubic self-repulsion may be stable, which is indeed true (see Section 3 below). The 3D equation (10) with ν > 1 and self-focusing sign of the nonlinearity produces the same scaling relation (18). However, it is expected that the respective solution families are unstable, as the relation does not satisfy the VK criterion.

Solitons pinned to a singular potential
In the 1D setting, singular solitons may be regularized, leading to localized states with a finite norm, if a deltafunctional attractive potential with strength ε > 0 is added to Eq. (10): An exact solution to Eq. (19) is with shift ξ, which removes the singularity in solution (20), defined by relation The finite amplitude of the regularized solution (20) is As it follows from Eq. (20), such modes, pinned to the delta-functional attractive potential, exist in a finite interval of negative values of the chemical potential: Exact solutions for solitons pinned to the same potential, but produced by Eq. (19) with the opposite (self-focusing) sign in front of the nonlinear term, were recently considered in Ref. [56].
For the exact solution given by Eqs. (20) and (21), it is easy to calculate the norm in the case of the cubic nonlinearity, ν = 1: Note that the N (µ) dependence given by Eq. (24) obviously satisfies the anti-VK criterion.
It is relevant to mention that the attractive delta-functional potential (generally, a complex one, which includes a local gain) may also produce stable dissipative solitons (in particular, exact ones) in the framework of the 1D complex Ginzburg-Landau equation, i.e., NLSE with complex coefficients in front of the second-derivative and nonlinear terms [57][58][59][60].

New results included in the review
Recent work [53] has reported results for families of stable 1D, 2D, and 3D singular solitons produced by Eq. (10) with, respectively, septimal (seventh-order, corresponding to ν = 3), quintic, and cubic self-repulsive nonlinearities. The results, which are directly relevant to the topic of the present review, are summarized in Section 3. In the same section we explain that all the relevant nonlinear terms, including the seemingly "exotic" septimal one, naturally occur in physical media (chiefly, in optics). An essential conclusion is similar to that presented in Section 2 for the singular states pulled to the attractive center: the repulsive nonlinearity helps to create stable singular GSs with finite norms. On the other hand, 2D states with embedded vorticity can be found in the model with the attractive sign of the quintic nonlinearity, but they are completely unstable.
with real radial function obeying the equation At r → 0, an asymptotic expansion of a relevant solution of Eq. (26) is cf. Eq. (5). Obviously, the density singularity corresponding to this asymptotic solution, u 2 ∼ r −4/3 , is weak enough to make the 2D integral norm (3) convergent at r → 0. Asymptotic expression (28) suggests substitution in Eq. (26), from which one can derive an equation for the singularity-free radial functions χ(r): Accordingly, the asymptotic form (30) of the solution at r → 0 is replaced by a singularity-free expansion, Usually, the presence of integer vorticity l ≥ 1 implies that the amplitude vanishes at r → 0 as r l , which is necessary because the phase of the vortex field is not defined at r = 0. However, the indefiniteness of the phase is also compatible with the amplitude diverging at r → 0. In the linear equation, this is the divergence of the standard Neumann's cylindrical function, Y l (r) ∼ r −l , which makes the respective 2D state unnormalizable for all values l ≥ 1. In the present case, Eq. (28) demonstrates that the interplay of the attractive potential and quartic self-repulsion curtails the divergence to the level of u ∼ r −2/3 , for all values of l, thus securing the normalizability of all the states under the consideration. This conclusion may be compared to what was found in Ref. [2], where the quintic repulsive term produced another integrable singularity of the 2D density, with u(r) ∼ r −1/2 .
The asymptotic form of the solution, given by Eq. (31) is meaningful if it yields χ(r) > 0 [otherwise, Eq. (26) cannot be derived from Eq. (9)], i.e., for U l > 0, as well as for −4/9 < U l < 0. The latter interval implies, according to Eq. (27), For the vortex states, with l ≥ 1, condition (32) means, in any case, U 0 > 0. However for the GS with l = 0, Eq. (32) admits an interval of negative values of U 0 , namely, While the existence of the bound state under the combined action of the repulsive potential, with U 0 < 0, and the defocusing quartic nonlinearity is a counter-intuitive finding, it is closely related to the above-mentioned fact, first reported in Ref. [50], that the 2D equation with the self-defocusing nonlinearity acting in the free space (without any potential) admits singular solutions with asymptotic form (15). If potential (2) is added to the 2D version of Eq. (10), expression (15) is replaced by which is tantamount to Eq. (28) in the case of ν = 3/2. In the limit of r → ∞, the asymptotic form of the solution to Eq. (30) is where χ 0 is an arbitrary constant, and µ must be negative. Then, a coarse approximation for the global solution can be obtained as an interpolation bridging asymptotic expressions (31) and (35): As mentioned above, families of localized states are usually characterized by dependences N (µ). In particular, calculating norm (3) for the approximate solution given by Eq. (36), one obtains (U l + 4/9) 2/3 where Γ(2/3) ≈ 1.354 is the value of the Gamma-function. Note that this dependence satisfies the anti-VK criterion.
Note that TF radius r 0 keeps the same value, as given by Eq. (38), in the presence of the MF defocusing cubic term with σ = 1 in Eq. (26), although the shape of the GS is more complex than one given by Eq. (38) for σ = 0. In this case, the asymptotic limit of the respective N (σ=1) TF (µ) dependence at µ → −∞ is the same as given by Eq. (39), while in the limit of µ → −0 it features a weaker singularity: Even in the case of the focusing sign of the MF term, corresponding to σ = −1 in Eq. (26), the LHY-induced quartic nonlinearity is able to stabilize the condensate against the combined action of the MF self-attraction and pull to the center. In this case, the TF approximation, applied to Eq. (26), cannot be easily resolved to predict u TF (r), but it readily produces an inverse dependence, for r as a function of u: Then, looking for a maximum of expression (41), which is attained at u max = 2/3, it is easy to find the corresponding size of the TF state: Equation (42) suggests that the GS exists, in the case of σ = −1, for |µ| exceeding a threshold value, |µ| > (|µ|) thr = 4/27.
According to Eq. (42), the norm diverges at µ → −4/27 as The analytical predictions reported in this subsection are compared to their numerical counterparts in the following one.  added to the stationary states ( * stands for the complex conjugation), the stability condition being, as usual [61], that all eigenvalues Λ must be pure-imaginary. The so predicted (in)stability was verified by direct simulations of Eq. (9). The results were produced for the model including the cubic term in Eq. (9), with σ = ±1, as well as for the most fundamental case of σ = 0, when the nonlinearity is represented solely by the LHY term.
A typical stable GS, obtained as a numerical solution of Eq. (30) with σ = 0, U 0 = 10, and µ = −1, is displayed in Fig. 1, together with its counterpart produced by the TF approximation (the interpolation-based approximation, given by Eq. (36), is not displayed here, as it is not relevant for large values of U 0 ). It is seen that the TF approximation is very close to its numerical counterpart in the inner zone, r < r 0 (see Eq. (38)), while in the outer one, r > r 0 , the TF approximation yields zero, being inaccurate, in this sense. Actually, the overall impact of the difference between the numerical and TF solutions in the outer zone is less significant due to the presence of factor r −2/3 in expression (29) for the full solution, u(r). As a result, in the case presented in Fig. 1  U 0 belonging to interval (33), was also confirmed by the numerical results. Figure 3 displays numerically found GSs for U 0 = −0.40, taken close to the limit value, U 0 = −4/9 ≈ −0.44, for all the three values of the coefficient in front of the MF cubic term in Eq. (30), viz., σ = 0 and ±1. The same figure shows that the interpolating approximation for these solutions, provided by Eq. (36), is quite accurate in this case (the TF approximation (38) does not apply to U 0 < 0). It is seen too that the MF cubic term does not strongly affect the solution. Dependences N (µ) for families of the GS solutions, obtained from Eq. (30) without and with the repulsive or attractive MF cubic term (σ = 0 and σ = ±1, respectively), are displayed in Figs. 4(a,b), for different values of strength U 0 of the potential, both U 0 > 0 and U 0 < 0. In panel (a), the N (µ) curves produced by the TF approximation as per Eq. (39) are compared to their numerical counterparts. The same panel demonstrates that the interpolating approximation is very accurate for U 0 = −0.4 (while its accuracy is poor for U 0 > 0). Panel (c) in Fig. 4 confirms that, in the presence of the attractive cubic term (σ = −1), the GS exists at |µ| > (|µ|) thr , as predicted by the TF approximation in Eq. (43). Up to the accuracy of the numerical results, the threshold value is indeed (|µ|) thr = 4/27, as given by Eq. (43). This finding is explained by the fact that the width of the GS diverges in the limit of |µ| → (|µ|) thr , as seen in Eq. (42), hence in this limit the derivatives become negligible in Eq. (26), making the TF approximation asymptotically exact. In principle, N (µ) must steeply diverge at |µ| → 4/27, according to Eq. (44), but it is difficult to collect numerical data very close to the threshold, as the GS is extremely broad in this limit.
The computation of eigenvalues for small perturbations, as well as direct simulations, demonstrate full stability of the GS solutions for both positive and negative values of U 0 (at which the GS exists), at all µ < 0. In particular, Fig.  3(b) demonstrates the stability in the counter-intuitive case of the repulsive potential, with U 0 = −0.4. The stability is not affected either by the TF cubic term, holding for σ = 0 and ±1.
Lastly, an analytical consideration, verified by numerical data [45], has demonstrated that vortex modes with l = 1 (see Eq. (25)), are stable at U 0 > 14/9, and unstable at U 0 < 14/9. For the vortices with l = 2, the critical value of the pulling-potential strength, below which they are unstable, is much larger, viz., U 0 = 77/9. Simulations of the evolution of unstable vortex modes demonstrate that the vortical pivot drifts in the outside direction, along an unwinding spiral trajectory. Eventually, the pivot is ousted to periphery, and the vortex mode transforms into a stable GS with zero vorticity. In this case, the norm of the residual GS is essentially smaller than the original value, due to intensive emission of small-amplitude waves in the course of the transient evolution.

III. Singular solitons in one, two, and three dimensions
This section summarizes original results which were recently reported in Ref. [53]. In the 2D case, the results are related to those presented in the previous section for the model with the external potential.
As mentioned in the Introduction, relevant self-defocusing nonlinearities which are necessary to support singular selftrapped states in 3D, 2D, and 1D settings, are represented by the cubic, quintic, or septimal terms, that correspond, respectively, to ν = 1, 2, and 3 in Eq. (10). Realizations of the cubic and quintic nonlinearities of either sign (focusing or defocusing) in diverse physical media are well known [62][63][64][65]. In particular, such nonlinear terms with controllable strengths (including the "exotic" septimal term), can be realized in optical waveguides filled by suspensions of metallic nanoparticles, control parameters being the density and size of the particles [66][67][68].
Results outlined below provide not only solutions in analytical and numerical forms, but also an interpretation of the physical meaning of the singular solitons.
A. Analytical results

The one-dimensional model with the septimal nonlinearity
Singular solitons created by the 1D version of Eq. (10) with the seventh-order defocusing (σ = 1) nonlinearity, ν = 3, are looked for as with real function U(x) satisfying equation In the application to the planar optical waveguide, t is not time, but the propagation distance (often denoted z), x is the transverse coordinate, and −µ is the propagation constant. The exact solution of Eq. (47) is given by Eq. (12): The asymptotic form of solutions to Eq. (47) at x → 0 does not depend on µ, Note that expression (49) is an exact solution of Eq. (47) with µ = 0, but its integral norm diverges at |x| → ∞. For µ < 0, the exponentially decaying asymptotic form of Eq. (48) at |x| → ∞ is Lastly, the norm of the 1D soliton family is given by Eq. (13), with the numerical coefficient which is accidentally close to π. This N (µ) dependence satisfies the anti-VK criterion, which, as mentioned above, is necessary for the stability of localized modes supported by repulsive nonlinearities [54].
2. Physical interpretation of the 1D singular soliton: screening of a "bare" δ-functional potential Although the existence of the stable singular solitons under the action of the septimal self-repulsive nonlinearity is firmly established by the above analysis, this result may seem counter-intuitive, as it is commonly believed that localized modes may only be supported by self-attraction. The purport of the result may be understood by comparing Eq. (47) to a modified equation, which includes an attractive delta-functional potential with strength ε > 0, cf. Eq. (19). An exact solution to Eq. (52) is produced by Eqs. (20), (21), and (22): (recall A ≡ |U(x = 0)| is the amplitude of the mode pinned to the attractive delta-functional potential). The approximate value of the offset in Eq. (54) corresponds to the limit of ε √ −2µ. This solution may be considered as a regularized version of its singular counterpart (49), which obviously converges to the singular state for ε → ∞. It is also relevant to calculate the value of the Hamiltonian, corresponding to Eq. (52), In the limit of ε √ −2µ, it is the negative sign suggesting that solution may be the GS.
A physical interpretation of the ability of the self-repulsive 1D model to create singular solitons is offered by the fact that, according to Eq. (54), the strength ("bare charge") ε of the delta-functional attractive potential diverges in the limit of ξ → 0, which brings one back to the underlying singular solution given by Eqs. (49)-(50). This observation implies that an infinitely large "bare attractive charge", embedded in the self-defocusing septimal medium, is completely screened by the nonlinearity, which builds the singular soliton with the convergent norm, for that purpose. This mechanism roughly resembles the renormalization procedure in quantum electrodynamics, where an infinite bare charge of the electron cancels with other diverging factors, making it possible to produce finite observable predictions.

The two-dimensional model with the quintic nonlinearity
In the 2D setting, it is relevant to consider Eq. (10) with the quintic self-defocusing term: In terms of optics, Eq. (58) models the paraxial propagation of light in a bulk waveguide, with t being the propagation distance (usually denoted z).
In polar coordinates (r, θ), solutions of Eq. (58) with integer vorticity, l = 0, 1, 2, ..., is looked for as with real amplitude function U(r) satisfying the radial equation, cf. Eq. (47): For 2D singular solitons with l = 0, Eq. (60) produces the following expansion at r → 0: with the respective 2D integral norm converging at r → 0. For µ = 0, is an exact solution of Eq. (60), but its norm diverges at r → ∞. The asymptotic form of the solution at r → ∞ is found from the linearized version of Eq. (60), [cf. Eq. (50)] where C is a constant, and the second term in the parenthesis its a correction to the lowest approximation.
For vortex states given by Eq. (59) with l ≥ 1, a singular solution with convergent norm is obtained with the opposite (self-focusing) sign in front of the quintic term in Eq. (58), the asymptotic approximation at r → 0 being cf. Eq. (61). For µ = 0, Eq. (64) represents an exact solution of Eq. (60) with the opposite (self-focusing) sign in front of the quintic term, but its integral norm diverges at r → ∞. An exact scaling relation for the norm of the solutions, with l = 0 and l ≥ 1 alike, follows from Eq. (58) with either sign of the quintic term: (for l = 0, a numerically found value is const ≈ 21.2). This dependence satisfies the anti-VK criterion, dN/dµ > 0 hence the respective GS solutions with l = 0, maintained by the defocusing quintic nonlinearity, may be (and indeed are) stable, but it contradicts the VK condition per se, dN/dµ < 0, hence the vortex states, which exist in the case of self-focusing, are definitely unstable, as corroborated by numerical simulations [53].

Interpretation of the 2D singular soliton: screening of a ring-shaped attractive potential
Similar to what is outlined above for the 1D model, it is possible to introduce a version of Eq. (58) with a deltafunctional potential; however, it is concentrated on a ring of a small radius, ρ, instead of the single point. The so modified equation (60) with l = 0 is 1 2 cf. Eq. (52). At r > ρ the solution of Eq. (66) is sought for in a form similar to that given by Eq. (61), while inside the ring it is taken as U = const ≈ 2 −3/4 ρ −1/2 . Eventually, straightforward manipulations, which take into regard the jump of dU/dr at r = ρ, yield a relation between the strength of the delta-functional potential concentrated on the ring and the ring's radius: cf. Eq. (54). Then, in the limit of ρ → 0, an effective "charge" of the 2D attractive potential is Thus, the 2D singular soliton may be construed as a solution providing the screening of the finite "charge" by the defocusing quintic nonlinearity.

Effects of additional nonlinear terms on 1D and 2D singular solitons
Even if the septimal or quintic term represents the dominant nonlinearity in the underlying NLSE, the lower-order terms, i.e., cubic and/or quintic ones (with respective coefficients g 3 and g 5 ), should be included in the realistic model of the light propagation [66][67][68]. In particular, the accordingly amended septimal 1D equation (47) is replaced by The additional terms produce negligible corrections to the asymptotic singular form (49) at x → 0 [53]: Further, if, in the 2D setting, the cubic term (the same as in Eq. (69)) is added to Eq. (60), it also produces a negligible correction to the singular asymptotic form given by Eq. (61) [53]: δU 2D (r) ≈ −2 −5/4 g 3 r 1/2 , cf. Eq. (70).
6. The 3D model with the cubic nonlinearity In the 3D case, the relevant equation is the NLSE with the usual cubic term: It has a plethora of physical realizations for both the defocusing (σ = 1) and focusing (σ = −1) signs of the nonlinearity [64], including the GPE for BEC with, respectively, repulsive or attractive interactions between atoms [10]. For isotropic stationary states, ψ = exp(−iµt)U(r), where r is the 3D radial coordinate, the real radial functions obeys the commonly known equation, For the asymptotic consideration of singular solutions at r → 0, the term on the right-hand side of Eq. (72) is negligible. Dropping this term and looking for solutions in the form of where r 0 is an arbitrary radial scale, one can transform Eq. (72) into the following form (which is an exact transformation for µ = 0): For σ = 1 (defocusing), the appropriate solution to Eq. (74) is also determined by the balance of the potential force and friction, the relevant region being τ → −∞, i.e., r r 0 : Both solutions (75) and (76) eventually translate into the asymptotic expression for the 3D density given above by Eq. (16).
In the opposite limit of r → ∞, a straightforward consideration yields an asymptotic expression for the solution in the form of with constant C, cf. Eq. (63). An exact scaling relation between the 3D norm and chemical potential, as it follows from the cubic NLSE (71), is the same as its counterpart (65) in the 2D model with the quintic nonlinearity: This relation satisfies the anti-VK criterion, hence the family of the 3D singular solitons may be (and indeed is) stable in the case of the defocusing cubic nonlinearity, σ = 1.
7. Interpretation of the 3D singular solitons: screening of an attractive spherical potential Similar to what is presented above for the 2D model, one can augment the 3D model by the attractive deltafunctional potential concentrated on a sphere of radius ρ, the respective stationary equation being (cf. Eq. (66)) 1 2 where it is assumed that, although ρ is small, it must be larger than r 0 in Eq. (73). A stationary solution to Eq. (79) is sought for in the form given by Eqs. (73) and (76) at r > ρ, and as U = const ≈ 2ρ ln (ρ/r 0 ) −1 at r < ρ. Then, calculations similar to those outlined above in the 2D setting lead to a relation between ρ and the strength of the attractive potential, ε 3D = 1/ (2ρ), cf. Eq. (67). In the limit of ρ → 0 (which is imposed along with r 0 → 0, so as to keep condition ρ > r 0 valid), the respective "charge" of the 3D attractive potential is Q 3D = 4πρ 2 · ε 3D = 2πρ → 0, cf. Eq. (68) Thus, one may realize the 3D singular soliton as a state which provides the screening of the vanishingly small "charge".
IV. Numerical results for the 1D, 2D, and 3D singular solitons The numerical scheme for producing singular solitons as solutions of the NLSEs with the self-repulsive nonlinearity must be adjusted to the fact that, in the analytical form, the solutions take infinite values at the origin. In Ref. [53], a finite-difference scheme was used for this purpose. It was defined on a grid with spacing ∆, constructed so that closest to the origin were points with coordinates (x, y, z) = (±∆/2, ±∆/2, ±∆/2) in the 3D case, and similarly in 1D and 2D. At these points, boundary conditions with large but finite values of |ψ| were fixed according to the asymptotically exact analytical expressions (49), (61), and ( 16). A typical shape of the 1D singular soliton, produced by exact solution (48) with µ = −1, is plotted in Fig. 5 with stepsize ∆x = 10 −5 . In agreement with the fact that the entire family of the 1D singular solitons satisfies the anti-VK criterion, numerical results corroborate the full stability of the family. In the 2D setting, numerical solution of the quintic equation (58) has corroborated the existence and full stability of the singular solitons with zero vorticity (l = 0). As an illustration, Fig. 6 displays a global view of the stable 2D soliton produced by direct simulations of Eq. (58), starting from the initial conditions taken as per exact solution (62) corresponding to µ = 0 (the formal divergence of its integral norm at r → ∞ is restricted by the finite size of the integration domain). In particular, he numerical simulations confirm the stability of the 2D singular solitons against azimuthal perturbations which attempt to break the axial symmetry of the solitons.
In the 3D setting, numerical solution of Eq. (71) produces a family of 3D singular solitons. Comparison of these solutions with the analytical prediction, given by Eqs. (73) and (76) for small r, is not straightforward, as the analytical expression contains indefinite parameter r 0 . Therefore, an example, displayed in Fig. 7(a) presents the comparison of the numerical solution to analytical profile const · r −1 at small r, with the constant selected as the best-fit parameter (const = 0.001 for the case of µ = −1, displayed in Fig. 7(a)).
Further, one can try to use the asymptotic form of the 3D solution at large r, given by Eq. (77), with C set equal to the best-fit value of const selected at r → 0, as a global analytical approximation. Figure 7(a) demonstrates that such an approximation is extremely close to its numerical counterpart at all values of r.
Finally, direct simulations of the evolution of the 3D singular solitons confirm stability of the soliton family, see an example in Fig. 7

V. Conclusion
The aim of this paper is to produce a brief review of findings recently reported in works [53] and [45], that demonstrate the existence of several types of singular bound states in 3D, 2D, and 1D models with self-repulsive nonlinearities. These states feature a singular structure of the density at r → 0, while the total norm converges, making the states physically relevant solutions. One model combines the attractive potential ∼ −1/r 2 in the 2D space and the LHY (quartic) self-repulsive term. The respective GPE (Gross-Pitaevskii equation) is a model of BEC composed of particles carrying permanent electric dipole moments, which are pulled to the central charge. Previously, it was found, in the framework of the MF and many-body quantum settings alike [2,11], that the quantum collapse, driven by this attractive potential, can be effectively suppressed in the 3D space by the usual cubic nonlinearity (which represents, as usual, repulsive interactions of colliding particles), while the same cubic nonlinearity cannot stabilize the 2D condensate. A brief review of those results was presented in Ref. [23].
The new result, reported in Ref. [45] and summarized in Section 2 of the present review, is that the quartic selfrepulsive term is sufficient to suppress the 2D quantum collapse and restore the missing GS (ground state), which is a completely stable one. States with embedded angular momentum (vorticity) are constructed too. They are stable if the strength of the attractive potential exceeds a certain threshold value, which depends on the vorticity. An essential peculiarity is that stable 2D GS modes exist, counter-intuitively, even when the attractive central potential is replaced by the repulsive one, with a sufficiently small strength, as defined by Eq. (33). The latter feature is closely related to another topic relevant to the consideration of singular states, which was elaborated in Ref. [53] and is summarized in Section 3 of the review. It deals with NLSEs (nonlinear Schrödinger equations) in 1D, 2D, and 3D free space, which contain the septimal, quintic, and cubic repulsive terms, respectively. These equations demonstrate the existence of stable singular solitons, that realize the model's GS in all the cases, with the 1D solitons found in an exact analytical form. The physical interpretation of these counter-intuitive findings is provided by the possibility to construe the singular solitons as a result of screening, by the respective self-repulsive nonlinearity, of a delta-functional attractive potential, whose integral strength is divergent in 1D, finite in 2D, and vanishingly small in 3D.
A common feature of both models considered in this paper is that, on the contrary to the stable zero-vorticity GSs, two-dimensional vortex states are completely unstable. Therefore, a challenging problem for further work is search for physically relevant modifications of the 2D models which may allow the existence of stable vortices. A still more challenging objective is to construct 3D states with embedded vorticity and investigate their stability. Another subject for further work may be singular dissipative solitons in complex Ginzburg-Landau equations [53]. Quite challenging may also be extension of the analysis for singular modes (if any) in models of Fermi gases based on density-functional equations.