1 Introduction↩︎

Ultralight bosonic dark matter, also known as Bose–Einstein condensate dark matter (BECDM) or fuzzy dark matter, provides a well–motivated framework in which dark matter is described by a coherent macroscopic wave function governed by the Gross Pitaevskii Poisson (GPP) system of equations. This model has been extensively explored as an alternative to cold dark matter at large and small scales, whose most important property is that wave effects can suppress small-scale structure and its collapse leads to the formation of solitonic cores [1]–[7].

Stationary solutions of the GPP system play a central role in this model, as they define important configurations. The simplest of them being the ground-state solution [8], [9], which acts as an attractor of this dark matter under very general initial conditions, from simple spherical and axisymmetric cases [10], [11], until scenarios of structure formation (e.g. [12]–[15]), as well as at local single structure collapse gravitational cooling or kinetic relaxation [16]–[18] or multi-core mergers [19]–[21]. Their properties, including scaling relations, stability criteria, and maximum mass configurations, have been widely studied both numerically and analytically [18], [22]–[25].

Beyond these spherical ground-state solutions there are others that may play an astrophysically relevant role, like the mixed state solutions called gravitational atoms [26], which have some chances to explain the anisotropy of satellite galaxies around major galaxies [27]. However, there are other intriguing solutions that may imprint some observational effects that may support or rule out BECDM, these are the vortex solutions. These states carry quantized angular momentum and exhibit a density depletion along the symmetry axis. Vortex configurations have been studied extensively in laboratory Bose-Einstein condensates [28]–[31], and have entered the astrophysical context within BECDM, where they may arise in rotating halos or during dynamical processes [32]–[36]. Recent studies have also explored their stability and potential observational signatures in galactic systems [37]–[39]. However, their structure and stability remain less well understood than those of non-rotating solitonic cores.

On the other hand, an additional ingredient of astrophysical relevance is the presence of central compact objects. Observations indicate that most galaxies host supermassive black holes at their centers [40]–[42], and their interaction with ultralight dark matter has been the subject of increasing interest. Previous works have studied accretion processes, dynamical friction, and the impact of black holes on scalar field configurations [43]–[48]. In fact, recent numerical studies have investigated the interplay between black holes and BECDM in dynamical scenarios, including mergers and binary systems [49]–[52]. A key point in these analyses is the interplay between the black-hole horizon scale and the characteristic length scale of the dark matter, namely the de Broglie wavelength. For ultralight bosons in the mass range relevant to BECDM, this wavelength is typically in the range \(\sim 10\,{\rm pc}\) to \(\sim 1\,{\rm kpc}\), depending on the boson velocity dispersion and mass. As a consequence, the black-hole horizon is generically much smaller than the coherence scale of the condensate, and local horizon-scale physics is expected to have a negligible impact on the large-scale dynamics of the dark matter wave function. In this regime, the black hole can be consistently modeled as an effective point-like Newtonian source, whose effect enters through the external gravitational potential rather than through detailed horizon-scale interactions and accretion. This justifies our treatment of the black hole as a Newtonian point mass throughout this work.

In this framework, we study the interplay between central black holes and vortex solutions of BECDM. We construct families of stationary axisymmetric solutions of the GPP system in the presence of a central black hole, focusing on both solitonic cores (\(m=0\)) as a test case and the most relevant scenario involving line-vortex configurations (\(m=1\)). Using an imaginary-time evolution method, we construct equilibrium solutions across a broad region of parameter space in terms of self-interaction and black hole mass, and analyze their structure and stability.

We make three main contributions. First, we provide a unified description of stationary configurations including both non-rotating and rotating states in the presence of a central potential sourced by a point-like mass. Second, we characterize their stability using scale-invariant parametrizations based on the turning-point criterion. Finally, we derive empirical scaling relations for the maximum mass and the onset of instability as functions of the black-hole mass.

By showing that there are stable solutions for a vortex of BECDM around a supermassive black hole, for a wide range of parameters, we motivate the study of potential formation processes and astrophysical observable signatures.

The paper is organized as follows. In Sec. 2, we introduce the Gross–Pitaevskii–Poisson system and its stationary equations, including the associated energy functional and scaling properties. In Sec. 3, we describe the numerical method based on imaginary–time evolution. The main results, including the structure and stability of the stationary solutions, are presented in Sec. 4. Finally, our conclusions are summarized in Sec. 5.

2 Model and equations↩︎

2.1 GPP Evolution Equations and Units↩︎

In our system, Ultralight Bose-Einstein Condensate Dark Matter (BECDM) evolves under its self-gravity and that of the gravitational field due to a point particle that models a central black hole. Thus, the dynamics of the condensate is governed by the coupled Gross–Pitaevskii–Poisson (GPP) system with the potential due to the black hole:

\[\begin{align} i\hbar \frac{\partial \tilde{\Psi}}{\partial \tilde{t}} &= -\frac{\hbar^2}{2m_B}\tilde{\nabla}^2 \tilde{\Psi} + m_B \left( \tilde{V} + \tilde{V}_{\rm BH} \right)\tilde{\Psi} + \tilde{g}\,|\tilde{\Psi}|^2 \tilde{\Psi}, \\ \tilde{\nabla}^2 \tilde{V} &= 4\pi G \tilde{\rho} , \end{align} \label{eq:GPP-physical}\tag{1}\]

where \(\tilde{\Psi}\) is the macroscopic wave function or order parameter of the condensate, \(m_B\) denotes the boson mass, the nonlinear coupling \(\tilde{g} = \frac{4\pi\hbar^2 \tilde{a}_s}{m_B^{2}}\) models the short-range self-interactions characterized by the \(s\)-wave scattering length \(\tilde{a}_s\), the gravitational potential \(\tilde{V}\) is sourced by the condensate’s density \(\tilde{\rho} = m_B|\tilde{\Psi}|^2\), while the external potential \(\tilde{V}_{\rm BH} = -G\tilde{M}_{\rm BH}/\tilde{r}\) represents the gravitational field of a central black hole of mass \(\tilde{M}_{\rm BH}\).

As usual, the GPP system is rewritten in dimensionless code units, for which we introduce a characteristic length scale \(x_0=\mathrm{kpc}/\lambda\) and parametrize the boson mass as \(m_B=m_{22}\times10^{-22}\,\mathrm{eV}/c^2\). Then following for example [37], we define the rescaled variables

\[\begin{align} \tilde{\vec{x}} &= x_0\,\vec{x}, & \tilde{t} &= t_0\, t, & \tilde{\rho} &= \frac{M_0}{x_0^3}\, \rho, \nonumber \\ \tilde{V} &= v_0^2\, V, & \tilde{a}_s &= a_0\, a_s, & \tilde{M}_{\rm BH} &= M_0\, M_{\rm BH}, \label{eq:scales-def} \end{align}\tag{2}\]

where the associated characteristic scales are

\[\begin{align} t_0 &= \frac{m_B x_0^{\,2}}{\hbar} \simeq 5.096\times10^{-2} \left(\frac{m_{22}}{\lambda^2}\right)\,\mathrm{Gyr}, \nonumber \\ v_0 &= \frac{\hbar}{m_B x_0} \simeq 19.20 \left(\frac{m_{22}}{\lambda}\right) \,\mathrm{km}\,\mathrm{s}^{-1}, \nonumber \\ M_0 &= \frac{\hbar^2}{4\pi G m_B^2 x_0} \simeq 6.82\times10^{6} \left(\frac{\lambda}{m_{22}^{2}}\right)\, \mathrm{M}_\odot, \nonumber \\ a_0 &= \frac{G m_B^{3} x_0^{2}}{\hbar^2} \simeq 3.228\times10^{-75} \left(\frac{m_{22}^{3}}{\lambda^2}\right)\,\mathrm{cm}. \end{align}\]

Finally, with these definitions, the dimensionless GPP equations become

\[\begin{align} i\,\frac{\partial \Psi}{\partial t} &= H(\Psi)\,\Psi, \\ \nabla^2 V &= \rho , \end{align} \label{eq:GPP-dimensionless}\tag{3}\]

where the Hamiltonian operator is given by

\[H(\Psi)= -\frac{1}{2}\nabla^2 + V + V_{\rm BH} + g\,|\Psi|^2, \label{eq:H-dimensionless}\tag{4}\]

and the black hole potential reads \(V_{\rm BH}=-M_{\rm BH}/(4\pi r)\).

In our analysis we also exploit the following scaling symmetry of the GPP system, that includes the black hole mass as follows:

\[\begin{align} && \{t,\vec{x},\Psi,V,g,M_{\rm BH}\} \;\longrightarrow\; \nonumber\\ &&\nonumber\\ &&\{\lambda^{-2}t,\lambda^{-1}\vec{x}, \lambda^{2}\Psi,\lambda^{2}V, \lambda^{-2}g,\lambda\,M_{\rm BH}\}, \label{eq:scale-symmetry} \end{align}\tag{5}\]

which leaves the system (3 ) invariant.

2.2 Stationary system, energy functional, and stability criteria↩︎

Stationary solutions of the GPP system are obtained by assuming harmonic time dependence of the condensate wave function,

\[\Psi(\vec{x},t)=e^{-i\mu t}\,\psi(\vec{x}),\]

where \(\mu\) denotes the dimensionless eigenenergy. Substituting this ansatz into Eq. (3 ) leads to the following nonlinear eigenvalue problem

\[H(\psi)\,\psi = \mu\,\psi , \label{eq:GPE-stationary-H}\tag{6}\]

constrained by the Poisson equation \(\nabla^2 V = |\psi|^2\).

Energy functional. Equation (6 ) is the result of a constrained variational principle, and stationary solutions correspond to extrema of the enthalpy functional

\[\mathcal{H}[\psi] = E[\psi] - \mu\,N[\psi],\]

where \(N\) is the particle number,

\[N[\psi]=\int |\psi|^2\,d^3x,\]

which in this case coincides with the total mass \(M\equiv N\) of the stationary solution. For our system, the total energy has different contributions

\[E[\psi]=K[\psi]+I[\psi]+W[\psi]+W_{\rm BH}[\psi],\]

where

\[\begin{align} K[\psi] &= \int \frac{1}{2}|\nabla\psi|^2\,d^3x, \\ I[\psi] &= \int \frac{g}{2}|\psi|^4\,d^3x, \\ W[\psi] &= \int \frac{1}{2}V|\psi|^2\,d^3x, \\ W_{\rm BH}[\psi] &= \int V_{\rm BH}|\psi|^2\,d^3x , \end{align}\]

corresponding to the kinetic, self-interaction, self-gravitational, and black hole gravitational energy contributions, respectively. The stationarity condition \(\delta\mathcal{H}=0\) recovers the eigenvalue equation (6 ). This decomposition is the basis of the numerical method used to construct the solutions below.

Turning points. The scaling symmetry (5 ) that leaves the GPP system invariant induces the transformation

\[\psi_\lambda(\vec{x}) = \lambda^{2}\psi(\lambda^{-1}\vec{x}), \label{eq:lambda-scaling}\tag{7}\]

which maps stationary configurations into others physically equivalent and the relevant quantities scale as

\[\begin{align} && \{M,\mu,K,I,W,W_{\rm BH},E\} ~~~\rightarrow \nonumber\\ &&\nonumber\\ &&\{\lambda M,\lambda^{2}\mu, \lambda^{3}K,\lambda^{3}I, \lambda^{3}W,\lambda^{3}W_{\rm BH}, \lambda^{3}E\},\label{eq:scale-symmetry2} \end{align}\tag{8}\]

which implies that the combination \(\sqrt{|g|}\,M\) is invariant under rescaling, and therefore provides a scale-independent characterization of families of stationary solutions.

Along a continuous family of stationary solutions parametrized by \(\mu\), stability can change at turning points, where the mapping between \(\mu\) and the conserved quantity \(M\) has a maximum value, leading to the Vakhitov-Kolokolov condition [53]

\[\frac{dM}{d\mu}=0 . \label{eq:turning-point}\tag{9}\]

Since \(\sqrt{|g|}\,M\) is scale invariant, an equivalent condition can be expressed as

\[\frac{d(\sqrt{|g|}\,M)}{d(g\,\psi_{0})}=0, \label{eq:maxgM}\tag{10}\]

where \(\psi_{0}=\max\{|\psi|\}\). This means that critical points of \(\sqrt{|g|}\,M\) as a function of \(g\,\psi_{0}\) separate stable and unstable branches of stationary solutions. This turning point criterion has been widely used in the analysis of self-gravitating condensates and boson star configurations (see, e.g., Refs. [18], [25], [53]–[55]).

2.3 Axisymmetric stationary configurations↩︎

Axisymmetric stationary configurations of the GPP system are characterized by the integer winding number \(m\), which labels the topological charge. For each value of \(m\), stationary solutions are solutions of nonlinear eigenvalue problems sharing the same differential operator but subject to different boundary conditions at the axis of symmetry. We denote the eigenfunctions and eigenvalues of such problems by \(\phi_m\) and \(\mu_m\), respectively.

In cylindrical coordinates \((r_\perp,\varphi,z)\), the stationary wave function can be written as [34], [35], [37], [38], [54], [56]:

\[\psi_m(r_\perp,\varphi,z)=\phi_m(r_\perp,z)\,e^{im\varphi},\]

where \(\phi_m(r_\perp,z)\) is real and axisymmetric. Substitution into the stationary Gross-Pitaevskii equation yields the nonlinear eigenvalue problem

\[H_m(\phi_m)\,\phi_m = \mu_m\,\phi_m , \label{eq:axisymmetric-eigenproblem}\tag{11}\]

with the Hamiltonian operator given by

\[\begin{align} H_m(\phi_m) & = & -\frac{1}{2}\left( \frac{\partial^2}{\partial r_\perp^2} + \frac{1}{r_\perp}\frac{\partial}{\partial r_\perp} + \frac{\partial^2}{\partial z^2} - \frac{m^2}{r_\perp^2} \right) \nonumber \\ & +& V + V_{\rm BH} + g\,\phi_m^2 . \label{eq:Hm-operator} \end{align}\tag{12}\]

The difference between a solitonic core solution and a line-vortex solution is determined by the boundary condition at the symmetry axis \(r_\perp=0\). For vanishing winding number (\(m=0\)), regularity allows a finite central density, \(\phi_0(0,0)\), along with the conditions

\[\begin{align} \partial_z\phi(r_\perp,z=0)&=&0,\nonumber\\ \partial_{r_\perp}\phi(r_\perp=0,z) &=&0 \label{eq:bcmeq0}, \end{align}\tag{13}\]

at the equatorial plane and the axis respectively. This condition is sufficient even with the potential \(V_{BH}\) as demonstrated in the Appendix.

Now, in the case of nonzero winding number (\(m\neq0\)), regularity of the wave function in the presence of quantized circulation requires the condition

\[\phi_m(0,z)=0 ,\label{eq:bcmne0}\tag{14}\]

which enforces a zero density along the symmetry axis and defines a line-vortex.

In both cases, localization of the solution requires \(\phi_m\to0\) as \(r_\perp^2+z^2\to\infty\). This isolation condition allows the gravitational potential \(V\) to be determined from the Poisson equation sourced by \(\rho=\phi_m^2\). These conditions suffice to guarantee regularity at the origin as shown in the Appendix, also for \(m > 0\).

This formulation shows that solitonic cores and line vortices are solutions of the same stationary GPP system, differing only in their topological charge and the corresponding boundary condition at the symmetry axis. In the following section we describe the numerical framework used to construct these stationary solutions and to analyze their stability.

3 Numerical framework↩︎

3.1 Numerical methods↩︎

Imaginary-time evolution method. Stationary axisymmetric configurations of the GPP system are constructed using an imaginary-time evolution method [37], which relaxes the system toward stationary solutions by minimizing the energy functional at fixed normalization. All configurations, including solitonic cores (\(m=0\)) and line vortices (\(m\neq0\)), are obtained by evolving the same stationary equation, just using the two mentioned boundary conditions (13 ) and (14 ) respectively.

Working in dimensionless units and defining the imaginary time \(\tau = i t\), it is possible to define an evolution equation for the stationary amplitude \(\phi_m\):

\[\frac{\partial \phi_m}{\partial \tau} = -\left[ H_m(\phi_m) - \mu_m \right]\phi_m , \label{eq:imaginary-GPE}\tag{15}\]

where \(H_m\) is the nonlinear Hamiltonian operator associated with winding number \(m\). Equation (15 ) is discretized in space using second-order finite differences and integrated in imaginary time using a first-order explicit Euler scheme,

\[\phi_m^{n+1} = \phi_m^{n} - \Delta\tau\, \left[ H_m(\phi_m^{n}) - \mu_m \right]\phi_m^{n}.\]

This evolution can be interpreted as a gradient-descent method in function space used to solve the stationary problem. In analogy with finding a root of a function via gradient methods, in our case the goal is to find a configuration \(\phi_m\) such that

\[\frac{\delta \mathcal{H}}{\delta \phi_m} = 0,\]

where \(\mathcal{H}=E-\mu N\) is the constrained functional to be minimized. The functional derivative reads

\[\frac{\delta \mathcal{H}}{\delta \phi_m} = \frac{\delta E}{\delta \phi_m}-\mu\,\frac{\delta N}{\delta \phi_m} = H_m(\phi_m)\phi_m - \mu_m \phi_m .\]

and therefore, Eq. (15 ) can be written as

\[\frac{\partial \phi_m}{\partial \tau} = - \frac{\delta \mathcal{H}}{\delta \phi_m},\]

which shows that the imaginary-time evolution drives the system toward a critical point of \(\mathcal{H}\), i.e., a solution of the nonlinear eigenvalue problem \(H_m(\phi_m)\phi_m=\mu_m\phi_m\).

At each imaginary–time step the wave function is explicitly rescaled to fix its maximum amplitude,

\[\phi_m \rightarrow \frac{\phi_m}{\max(\phi_m)} .\]

This procedure exploits the scaling symmetry of the GPP system (5 ) and allows all stationary solutions to be constructed with \(\phi_0 :=\psi_{0}=1\). Physical solutions with arbitrary amplitudes are then recovered after each time-step by applying the corresponding scaling transformation.

The imaginary-time evolution is continued until convergence is achieved, as determined by the residual of the stationary equation,

\[\| H_m(\phi_m)-\mu_m \phi_m \|_2 < 10^{-4}.\]

Fourier–space Poisson solver. At each iteration during the evolution, we obtain the gravitational potential \(V\) in the Hamiltonian \(H_m\) by solving the Poisson equation:

\[\left( \frac{\partial^2}{\partial r_\perp^2} +\frac{1}{r_\perp}\frac{\partial}{\partial r_\perp} +\frac{\partial^2}{\partial z^2} \right)V = \rho, \qquad \rho = \phi_m^2. \label{eq:poisson-axisymmetric}\tag{16}\]

for which we apply a Fourier transform along the \(z-\)direction,

\[\hat{V}(r_\perp,k_z) =\mathcal{F}_z[V] = \int V(r_\perp,z)e^{ik_z z}\,dz,\]

and similarly \(\hat{\rho}=\mathcal{F}_z[\rho]\). Specifically we compute the transform and its inverse using the fast Fourier transform (FFT) as illustrated in [21], [57]. In Fourier space, Eq. (16 ) reduces to

\[\frac{d^2\hat{V}}{dr_\perp^2} +\frac{1}{r_\perp}\frac{d\hat{V}}{dr_\perp} -k_z^2\hat{V} = \hat{\rho}. \label{eq:poisson-fourier}\tag{17}\]

that we solve on the discrete radial domain given by \(r_i=i\,\Delta r_\perp\), which yields a tridiagonal system:

\[a_i \hat{V}_{i-1}+b_i \hat{V}_i+c_i \hat{V}_{i+1} = \hat{\rho}_{i},\]

with \(a_i=\frac{1}{\Delta r_\perp^2}-\frac{1}{2r_i\Delta r_\perp}\), \(b_i=-\frac{2}{\Delta r_\perp^2}-k_z^2\), and \(c_i=\frac{1}{\Delta r_\perp^2}+\frac{1}{2r_i\Delta r_\perp}\). For each Fourier mode, the system is solved using the Thomas algorithm [58], and the solution \(V(r_\perp,z)\) is reconstructed via the inverse FFT.

Boundary conditions are imposed as follows. Regularity at the axis requires

\[\left.\frac{d\hat{V}}{dr_\perp}\right|_{r_\perp=0}=0,\]

while at \(r_\perp=r_{\max}\) we impose a monopolar boundary condition

\[V(r_\perp,z)\big|_{r_\perp=r_{\max}} = -\frac{M}{4\pi\sqrt{r_\perp^2+z^2}}.\]

which have to be fulfilled.

3.2 Numerical setup and parameter space↩︎

All solutions are constructed on a finite two-dimensional domain \((r_\perp,z)\in[0,20]\times[-20,20]\) in dimensional units, uniformly discretized with a \(128\times128\) grid and spatial resolution \(\Delta z = 2\Delta r_\perp = 0.3125\).

The imaginary-time step is chosen according to a diffusive stability condition

\[\Delta\tau = 0.25\,\min(\Delta r_\perp,\Delta z)^2 = 2.5\times10^{-2},\]

which ensures stability and convergence of the explicit Euler scheme.

Parameter space. We construct stationary configurations only for winding numbers \(m=0\) and \(m=1\), corresponding to solitonic core and singly quantized line-vortex configurations. The self-interaction coefficient is varied in the range

\[g\in[-2,-0.5), \label{eq:g-range}\tag{18}\]

while the dimensionless black-hole mass covers the range

\[M_{\rm BH}\in[0,20). \label{eq:MBH-range}\tag{19}\]

In order to provide a physical reference, we recall that the dimensionless quantities are related to their physical counterparts through the scalings introduced in Sec. 2 . For the fiducial values \(\lambda=1\) and \(m_{22}=1\), the characteristic mass scale is \(M_0 \simeq 6.82\times10^{6}\,M_\odot\). The explored interval \(M_{\rm BH}\in[0,20]\) therefore corresponds to physical black-hole masses in the range

\[\tilde{M}_{\rm BH}\in[0,1.36\times10^{8}]\,{\rm M}_\odot ,\]

which goes from the absence of a central black hole up to masses typical of supermassive black holes in galactic nuclei. In particular, values \(M_{\rm BH}\sim 1\) correspond to \(\tilde{M}_{\rm BH}\sim 7\times10^{6}\,M_\odot\), comparable to the mass of the black hole at the center of the Milky Way (\(\mathrm{Sgr\,A^\ast}\)) [40]–[42], while the upper end of the interval reaches the regime of intermediate and moderately massive supermassive black holes.

Another boundary of our parameter space is the restriction to ground states. Although excited stationary states can, in principle, be obtained within the imaginary-time approach by enforcing orthogonality with respect to previously computed solutions, we restrict our analysis to the ground state for the two cases \(m=0,1\) [35], [55]. This choice is motivated by both numerical and physical considerations. It is well known that ground-state configurations for \(m=0\) are dynamical attractors in BECDM simulations [12]–[15], as well as in local collapse simulations [16], [18] and not only for pure BECDM, but also in the presence of gas acting as baryonic matter [59], and in the presence of black holes, which promote condensation toward stationary core solutions in the fuzzy dark matter regime [48]. In contrast, excited states typically exhibit multiple stability branches [10], leading to a more intricate stability structure (e. g. [55]) whose possible formation -and therefore astrophysical relevance- has not yet been studied.

Finally, we restrict our analysis to the attractive self-interaction regime, \(g<0\), which is also motivated by stability considerations. Stationary solutions correspond to critical points of \(\mathcal{H}\), with stability determined by the sign of \(\delta^2 \mathcal{H}\). A practical criterion is given by the Vakhitov-Kolokolov condition (9 ),\(\frac{dN}{d\mu} < 0\), which characterizes stable configurations, whereas \(dN/d\mu > 0\) signals instability. While the repulsive case \(g>0\) typically satisfies this condition along the ground-state branch, the attractive regime \(g<0\) admits unstable configurations, leading to a richer dynamical structure.

Once we have defined the numerical framework and parameter space, we proceed to present the corresponding stationary solutions. The following section analyzes their structural properties and stability as functions of the winding number, self-interaction strength, and black-hole mass.

4 Results↩︎

4.1 Structure of stationary solutions↩︎

We begin by illustrating the structure of the stationary solutions on the equatorial plane (\(z=0\)). Figures 1 and 2 show the density profiles along the radial coordinate for both solitonic core (\(m=0\)) and line-vortex (\(m=1\)) configurations.

Effect of the black hole mass. Figure 1 shows solutions for the self-interaction strength fixed to \(g=-2\) and increasing black-hole mass \(M_{\rm BH}\in[0,20)\). As \(M_{\rm BH}\) increases, the gravitational potential deepens, leading to a systematic contraction of the density profiles. For \(m=0\), this produces a stronger central concentration, whereas for \(m=1\) the density peak shifts inward while preserving the characteristic zero density at the symmetry axis.

Effect of self-interaction. Figure 2 shows solutions in the absence of a central black hole (\(M_{\rm BH}=0\)) for different values of the self-interaction strength \(g\in[-2,-0.5)\). As \(|g|\) increases, the attractive interaction becomes stronger and the profiles contract progressively. For \(m=0\), the configurations become more centrally peaked, while for \(m=1\) the density maximum moves toward the axis, maintaining the vortex line at the origin.

In all cases, the qualitative behavior among the cases \(m=0\) is the same, namely a finite central density with different compactness, whereas for \(m=1\) all solutions show a ring-like distribution of matter that vanishes at the axis due to the vortex nature of the solution.

4.2 Stability and invariant parametrization↩︎

In order to analyze the stability of the solutions, we exploit the scaling symmetry of the GPP system and construct scale-invariant quantities. In particular, we consider the invariant combinations \(\sqrt{|g|}\,M\) and \(g\,\phi_{0}\), where \(\phi_{0}=\max\{\phi\}\). Notice that this is only a normalization convention and does not reduce generality. Because the Gross-Pitaevskii-Poisson system has the scaling symmetry (5 ), physical solutions with arbitrary central amplitude can be reconstructed by applying the corresponding rescaling. For this reason, results are reported using the scale-invariant combination \(g ~\phi_0\), which is independent of the parameter \(\lambda\) in Eq. (5 ), and parametrizes the corresponding scale-equivalent families of stationary solutions generated from our parameter space. Therefore our results can be mapped into the corresponding scale-equivalent families of solutions with arbitrary central field values by using any \(\lambda \ne 1\), whose physical properties are determined by the scaling relations (5 ) and (8 ).

In Figure 3 we show the invariant mass \(\sqrt{|g|}\,M\) as a function of \(g\,\phi_{0}\) for a number of families corresponding to the solitonic core (\(m=0\)) and line-vortex (\(m=1\)) configurations, for different values of the black-hole masses \(M_{\rm BH}\).

Notice that for each value of \(M_{\rm BH}\), the solutions define a continuous family of solutions in the invariant plane. Each family exhibits a turning point, characterized by an extremum of \(\sqrt{|g|}\,M\), where the condition in Eq. (10 ) is satisfied. According to the turning-point criterion, this critical point separates stable from unstable branches of the family.

For the case \(m=0\), the turning point separates a stable branch at lower values of \(g\,\phi_{0}\) from an unstable branch at higher values. A similar structure is found for \(m=1\), although the corresponding branches are shifted due to the non-zero angular momentum of the solution and the associated vortex structure.

As \(M_{\rm BH}\) increases, the branches are systematically displaced in the invariant plane, indicating that the central black hole modifies both the maximum mass and the location of the stability threshold.

4.3 Critical points and scaling relations↩︎

The turning points identified in Fig. 3 define the maximum–mass configurations along each family of solutions. These points satisfy the condition in Eq. (10 ) and determine the onset of instability.

Extracting these maxima, we construct the dependence of the invariant quantities \((\sqrt{|g|}M)_{\max}\) and \((g\,\phi_{0})_c\) on the critical parameter \((\sqrt{|g|}M_{\rm BH})_c\). The resulting relations are shown in Figs. 4 and 5.

Maximum mass. Figure 4 shows \((\sqrt{|g|}M)_{\max}\) as a function of \((\sqrt{|g|}M_{\rm BH})_c\) for both \(m=0\) and \(m=1\). The numerical data are well described by

\[\begin{align} (\sqrt{|g|}M)_{\max}^{(m=0)} &= 12.8 - 11.9\left(1-e^{-0.158\,x^{0.998}}\right), \tag{20}\\ (\sqrt{|g|}M)_{\max}^{(m=1)} &= 34.0 - 33.7\left(1-e^{-0.0325\,x^{1.07}}\right),\tag{21} \end{align}\]

where \(x=(\sqrt{|g|}M_{\rm BH})_c\). In both cases, the maximum mass decreases as the black–hole mass increases, indicating that the central potential reduces the amount of mass that can be supported in equilibrium.

However, the relative variation differs significantly between the two configurations. Over the explored range, the normalized change satisfies

\[\frac{\Delta(\sqrt{|g|}M)_{\max}}{(\sqrt{|g|}M)_{\max}} \approx\left\{ \begin{array}{cc} 0.93, & m=0\\ &\\ 0.53, & m=1 \end{array} \right.\]

showing that solitonic cores experience a bigger relative change in maximum mass across the parameter space. This indicates that \(m=0\) configurations are significantly more sensitive to the central gravitational field, whereas vortex solutions (\(m=1\)) respond more moderately, consistent with a greater structural robustness likely associated with angular momentum support.

Critical amplitude. Figure 5 shows the critical invariant \((g\,\phi_{0})_c\) as a function of \((\sqrt{|g|}M_{\rm BH})_c\). The data can be accurately fitted by the formulas

\[\begin{align} (g\,\phi_{0})_c^{(m=0)} &=& -0.712 - 0.0142\,y^{1.50}, \tag{22}\\ (g\,\phi_{0})_c^{(m=1)} &=& -0.898 - 2.16\times10^{-6}\,y^{3.55},\tag{23} \end{align}\]

where \(y=(\sqrt{|g|}M_{\rm BH})_c\). For \(m=0\), the critical amplitude decreases with increasing black-hole mass, indicating that the onset of instability occurs at progressively lower central densities. For \(m=1\), the variation is much weaker, reflecting the stabilizing role of angular momentum in vortex configurations. A comment is in order regarding the phenomenological formulas (20 )–(23 ). These expressions are intended only as empirical fits that describe the behavior of the solutions within the parameter space explored in this work. They are not meant to be used as extrapolations beyond the range of parameters considered in our simulations.

Overall, these results show that the presence of a central black hole modifies both the maximum mass and the location of the stability threshold, with a stronger effect on solitonic core solutions than on line vortices.

These results are consistent with previously reported limits in specific regimes of the parameter space. In the absence of a central black hole (\(M_{\rm BH}=0\)), the maximum mass for \(m=0\) agrees with our previous determination based on a Sturm-Liouville formulation solved using genetic algorithms [25], which is in turn consistent with the results obtained via the shooting method in [18].

In the complementary limit of vanishing self-interaction (\(g=0\)) and finite black-hole mass, the resulting profiles are consistent with the stationary solutions reported in Ref. [48], also constructed using shooting techniques.

These consistency regimes validate the imaginary-time approach used across different regions of the parameter space and confirm its ability to recover known limiting solutions.

Moreover, these results are also in qualitative agreement with analytical scalings derived in the Thomas-Fermi and self-gravitating regimes [22], [23].

5 Conclusions↩︎

We have constructed stationary solutions of self-gravitating Bose-Einstein condensate dark matter in the presence of a central black hole, within the Gross-Pitaevskii-Poisson framework for solitonic cores and line-vortex configurations. For this we implemented a novel method that uses imaginary time evolution. We also analyzed stability of the families in a wide parameter space in a scale-invariant formulation.

We have shown that both the black–hole potential and the attractive self–interaction produce a systematic compactification of the solutions. Increasing either \(M_{\rm BH}\) or \(|g|\) leads to configurations that are more centrally concentrated. However, the physical origin of this behavior differs in each case: the former deepens the external gravitational potential, whereas the latter strengthens the intrinsic self-binding of the condensate.

The stability analysis, based on scale-invariant quantities and the turning-point criterion, reveals that the maximum mass along each family of solutions is strongly affected by the presence of the central black hole. In particular, the invariant mass \((\sqrt{|g|}M)_{\max}\) decreases monotonically with increasing \((\sqrt{|g|}M_{\rm BH})_c\).

Another interesting result is the distinct response of the two cases analyzed. While vortex configurations (\(m=1\)) exhibit a larger absolute variation in their maximum mass, solitonic cores (\(m=0\)) experience a stronger relative reduction, indicating a higher sensitivity to the central gravitational field. This behavior suggests that angular momentum provides a mechanism of structural support, making vortex solutions more robust against the influence of the black hole.

We have also shown that the dependence of the maximum mass and the critical amplitude on the black–hole mass can be modeled by simple empirical relations in terms of scale-invariant variables. These relations provide a compact characterization of the stability properties of the system and may serve as a tool for connecting theoretical models of ultralight dark matter with astrophysical environments hosting central black holes.

Finally a comment on the astrophysical relevance of these configurations is in turn. It has already been shown that solitonic solutions around a black hole can be formed and in fact are attractors within the Fuzzy Dark Matter model under very general conditions[48], thereby its relevance in SMBH-FDM astrophysics (e.g [43], [47], [60]–[63]). The line-vortex solutions on the other hand, have not been studied in detail in the presence of black holes, and since we have shown now in this paper that they have a stable branch in each family of parameters, there is enough motivation to follow two natural steps including a demonstration of their stability in a full 3D environment and the search of a possible formation mechanism of these solutions, at least locally near the \(z=0\) plane.

6 Data availability↩︎

The numerical results of all the equilibrium configurations constructed in this paper, as well as the data needed to reproduce all the plots are publicly available at [64].

Acknowledgments↩︎

This research is supported by SECIHTI Grant No. CFB-2025-I-759, Laboratorio Nacional de Cómputo de Alto Desempeño Grant No. 2026-8, and CIC-UMSNH Grant No. 4.9.

7 Regularity of solutions↩︎

7.1 Case \(m=0\)↩︎

Consider the axial equation for the core configuration with \(m=0\)

\[-\frac{1}{2}\left[\partial_{r_\perp}^2\phi_0+\frac{1}{r_\perp}\partial_{r_\perp}\phi_0+\partial_z^2\phi_0\right]+V_T(r_\perp,z)\phi_0=\mu\phi_0.\nonumber\]

The regularity conditions (13 )on the axis and equatorial plane are respectively \(\partial_{r_\perp}\phi_0(0,z)=0\) and \(\partial_z\phi_0(r_\perp,0)=0\). Assuming spherical symmetry of the total potential, \(V_T(r_\perp,z)=V_T(r)\) with \(r=\sqrt{r_\perp^2+z^2}\), the ground state inherits this symmetry so that \(\phi_0(r_\perp,z)=\phi(r)\).

Using the chain rule, \(\partial_{r_\perp}\phi_0=\phi'(r)\,r_\perp/r\) and \(\partial_z\phi_0=\phi'(r)\,z/r\), the regularity conditions are automatically satisfied. The Laplacian becomes \(\phi''(r)+\frac{2}{r}\phi'(r)\), and the equation reduces to

\[-\frac{1}{2}\left[\phi''(r)+\frac{2}{r}\phi'(r)\right]+V_T(r)\phi(r)=\mu\phi(r),\nonumber\]

or equivalently

\[-\frac{1}{2r^2}\frac{d}{dr}\left(r^2\frac{d\phi}{dr}\right)+V_T(r)\phi=\mu\phi.\nonumber\]

Now, defining \(u(r)=r^2\phi'(r)\) gives \(\phi'(r)=u/r^2\) and \(u'(r)=2r^2[V_T(r)-\mu]\phi\), the system becomes

\[\begin{align} \phi'(r)=u/r^2, \nonumber \\ u'(r)=2r^2[V_T(r)-\mu]\phi. \nonumber \end{align}\]

We recall that the regularity condition requires \(\phi(0)=\phi_c<\infty\) and \(u(0)=0\). If the potential contains a Newtonian term \(V_T(r)=-\alpha/r+V(r)\) with \(V(r)\) finite at the origin, then near \(r=0\) one has \(u'(r)=-2\alpha\phi_c r+\mathcal{O}(r^2)\), which integrates to \(u(r)=-\alpha\phi_c r^2+\mathcal{O}(r^3)\). Hence \(\phi'(r)=-\alpha\phi_c+\mathcal{O}(r)\) and therefore \(\phi(r)=\phi_c-\alpha\phi_c r+\mathcal{O}(r^2)\), implying finally that \(\lim_{r\to0}\phi(r)=\phi_c<\infty\).

7.2 Regularity of the Axial \(m>0\) Vortex Solution↩︎

For a vortex configuration with winding number \(m>0\), the axial stationary equation is

\[\begin{align} &&-\frac{1}{2}\left[\frac{\partial^2 \phi_m}{\partial r_\perp ^2} +\frac{1}{r_\perp} \frac{\partial \phi_m}{\partial r_\perp} - \frac{m^2}{r_\perp^2}\phi_m+\partial_z^2\phi_m\right] +\nonumber\\ &&~~~~~~~~~~~~~~ V_T(r_\perp,z)\phi_m=\mu\phi_m.\nonumber \end{align}\]

Near the symmetry axis, the dominant contribution is the centrifugal term \(-m^2\phi_m/r_\perp^2\), so to leading order the equation reduces to \(\partial_{r_\perp}^2\phi_m+\frac{1}{r_\perp}\partial_{r_\perp}\phi_m-\frac{m^2}{r_\perp^2}\phi_m\simeq0\).

Using the ansatz \(\phi_m(r_\perp,z)\sim A(z)r_\perp^p\), and substituting gives \(p(p-1)r_\perp^{p-2}+pr_\perp^{p-2}-m^2r_\perp^{p-2}=0\), which implies \(p^2-m^2=0\), and finally \(p=\pm m\).

The branch \(p=-m\) is singular at the axis and it must be discarded, which results in regular behavior \(\phi_m(r_\perp,z)=A(z)r_\perp^m+\mathcal{O}(r_\perp^{m+2})\). This implies that \(\lim_{r_\perp\to0}\phi_m(r_\perp,z)=0\), so that the vortex boundary condition reduces to the condition \(\phi_m(0,z)=0\).

Finally, the density behaves as \(\rho_m=|\phi_m|^2\sim r_\perp^{2m}\), which vanishes on the axis, \(\lim_{r_\perp\to0}\rho_m=0\).

References↩︎