Dynamical galactic effects induced by solitonic vortex structure in bosonic dark matter

The nature of dark matter (DM) remains one of the unsolved mysteries of modern physics. An intriguing possibility is to assume that DM consists of ultralight bosonic particles in the Bose–Einstein condensate (BEC) state. We study stationary DM structures by using the system of the Gross–Pitaevskii and Poisson equations, including the effective temperature effect with parameters chosen to describe the Milky Way galaxy. We have investigated DM structure with BEC core and isothermal envelope. We compare the spherically symmetric and vortex core states, which allows us to analyze the impact of the core vorticity on the halo density, velocity distribution, and, therefore, its gravitational field. Gravitational field calculation is done in the gravitoelectromagnetism approach to include the impact of the core rotation, which induces a gravimagnetic field. As a result, the halo with a vortex core is characterized by smaller orbital velocity in the galactic disk region in comparison with the non-rotating halo. It is found that the core vorticity produces gravimagnetic perturbation of celestial body dynamics, which can modify the circular trajectories.


I. INTRODUCTION
The nature of DM particles remains one of the most fascinating puzzles of modern physics.The DM large-scale properties consistent with astrophysical observations are successfully explained by the cold dark matter model (CDM), which describes DM as a collisionless sufficiently cold perfect fluid.However, at smaller scales, the CDM encounters the cusp-core, missing satellites, and too-big-to-fail problems.One possibility to solve them is to assume that DM particles are ultra-light bosons as it is assumed in ultralight dark matter (ULDM) models [1].Generically, these models are characterized by the suppression of the small-scale structures, the presence of cores, and dynamic effects which arise from the BEC formed in the central regions of galaxies.Such DM halo proposals were investigated in [2][3][4][5][6].
The ULDM model is supported indirectly by observations.For example, in cosmological sim-ulations [7] it was found that the bosonic DM can indeed reproduce the observed distribution of matter at very large scales [8,9], though the mass of such bosons should be extremely small.There have been also studies on some tensions of the ULDM with observational data from the rotation curves of galaxies including the Milky Way, which could probe the particle mass in the range m = 10 −22 −10 −21 eV [10,11].Furthermore, the viability of the ULDM model was studied with the stellar kinematics measurements in dwarf galaxies [12].Another important piece of evidence is the DM nongravitational self-interaction, which has been recently reported for collisions of the clusters [13,14].In addition, the DM halo model must ensure the stability of a predicted halo.The stability of compact astrophysical objects which may be formed due to the Bose-Einstein condensation of ULDM was shown numerically [15].
In the present paper, we discuss DM, which consists of ultra-light bosons with repulsive self-interaction.
Such models make use of two macroscopic quantum phenomena: Bose-Einstein condensation and superfluidity.Bose-Einstein condensate in the mean-field approximation is described by the Gross-Pitaevskii equation.By adding dissipation in the Gross-Pitaevskii equation one obtains a more general model, which includes the effective temperature effect and predicts that the ULDM halo consists of a BEC core and an isothermal envelope [16].Such core-envelope structure in the ULDM model was also discussed in [17][18][19][20].Another important property, superfluidity, allows the quantization of the circulation and thus the possibility of the formation of vortices in the core of the halo.The central object of our study, the vortex, has a vanishing wavefunction at the vortex line, with a quantized circular flow around the vortex line [1].According to the recent numerical studies [21,22] only the non-rotating soliton and single-charged vortex are stable, even being strongly perturbed.In the present work, we consider a DM halo, which consists of two regions -core and isothermal envelope, while the core could be either a soliton or a single-charged vortex.
Most of our knowledge about DM is based on its gravitational interaction with baryonic matter objects.Thus, testing the validity of the UDM theory requires a detailed investigation of the DM gravitational field.The DM density distribution, predicted by ULDM models, has been extensively studied in numerical simulations and applied in studies aimed at reconstructing the gravitational potential of DM halos for the Milky Way [23] and dwarf galaxies [24].In general, one can determine the gravitational field of the ULDM by solving the Einstein equations with the DM density and rotation flow as sources of the gravitational field, where rotation flow is induced by the BEC superfluidity.Thus, in the ULDM model, we should be able to deduce the impact of the superfluid DM rotation on the observations.The dominant effect of the vortex existence is due to the different core density distributions.Moreover, rotation flows produce v/c and higher order effects, which can be taken into account in the gravitoelectromagnetism approach discussed in [25][26][27][28][29] and used in our calculations below.The gravitoelectro-magnetic formulation of a slowly rotating, selfgravitating, and dilute BEC intended for astrophysical applications in the context of DM halos was discussed in [30].As a rule, the gravimagnetic force is quite weak and does not affect significantly the dynamics of astrophysical systems.However, in the central region of the BEC core, the DM density vanishes while the vortex flow velocity dramatically increases, which can affect the dynamics of luminous matter in the central region of the galaxies.
In the present work, we calculate the DM gravitational field, which is needed for analysis of the observable predictions of the DM model, namely, to study how DM affects the movement of luminous matter.In our study, DM is the only source of a gravitational field, while luminous matter moves along geodesics, induced by DM.A more precise description of galactic kinematics is given by modeling the baryonic contribution to the gravitational potential which can distort the BEC soliton structures [31,32].Such a contribution was found to be significant for the Milky Way (MW) but not essential for the SPARC LSB galaxies [33].In this paper, we will limit ourselves to some simple consequences of the ULDM model on the galactic kinematics, namely, rotation curves and deviation of circular trajectory, induced by the gravimagnetic force.The more detailed study in this direction is beyond the scope of the current paper, though it is an interesting perspective on further work.
The paper is organized as follows.In Sec.II, we develop the key parameters of our model, define the equations for halo structure, and formulate the gravitoelectromagnetism ansatz.In Sec.III, we discuss the halo density profile for two stable core configurations and define the corresponding hydrodynamical velocity.In Sec.IV, the gravielectric (Newtonian) field of the halo is calculated and the rotational curves are obtained.Sec.V provides gravimagnetic field calculations and our estimates of the gravimagnetic effect on circular trajectory.The results are summarized in Sec.VI.

II. MODEL
A. Ultra-light dark matter model and halo structure In this section, we briefly discuss the model, suggested in [16].The structure of the DM halo is described by the Gross-Pitaevskii-Poisson (GPP) equations, which define the dynamical evolution of self-gravitating BEC field ψ where Xdr is the spatial average over halo, m is the bosonic particle mass, denotes reduced Planck constant, k B is Boltzmann constant.The first equation can be obtained by incorporating dissipative effects into the Schrödinger equation by means of the theory of scale relativity.This generalization of the Schrödinger equation means basically taking into account the interaction of the system with the external environment.The model Eqs.
(1),(2) were derived in [34], and we follow this approach in our current work.
We consider the BEC model with parameters γ = 2 and K = 2πas 2 m 3 , where a s denotes the s-wave scattering length of the selfinteraction.The parameter η 0 determines the equation of state of DM [35].The first term on the right-hand side of Eq.( 1) is the kinetic term, and the second describes the interaction with the condensate gravitational potential Φ g .The third term takes into account the bosonic selfinteraction (we will consider only the case γ = 2 which corresponds to binary collisions).The fourth term accounts for the core, and the fifth term describes an isothermal envelope with effective temperature T which surrounds the core.These terms can be derived from the Lynden-Bell theory of violent relaxation [16].The last term with ξ < 0 is a damping term and ensures that the system relaxes towards the equilibrium state.
An important feature of the Gross-Pitaevskii (GP) equation is that it satisfies the H-theorem, i.e., the free energy F of the system decreases Ḟ = −ξ ρu 2 dr ≤ 0.
where ρ = |ψ| 2 denotes BEC density and u = ∇S(r, t)/m is the velocity field.These quantities are obtained by application of Madelung transformation, according to the expression ψ(r, t) = ρ(r, t)e iS(r,t)/ , where S(r, t) is the action.The negative sign of ξ implies that the system relaxes towards the state with zero hydrodynamical velocity u = 0. Therefore, a stationary vortex solution with nonzero u can be found only if we set ξ = 0.
The free energy F = E − T S B is expressed through the total energy E, the effective temperature T , and the Boltzmann entropy S B = −k B (ρ/m)(ln ρ − 1)dr.The total energy consists of the classical kinetic energy Θ c = 1/2 ρu 2 dr, the quantum kinetic energy Θ Q = 1/m ρQdr, the gravitational potential energy W = 1/2 ρΦ g dr, and the internal energy of the self-interaction U = K ρ 2 dr, is the quantum potential.A stable equilibrium state corresponds to the minimum of the free energy F at fixed total mass M of BEC.This gives the following condition of quantum hydrostatic equilibrium [16]: where P = Kρ 2 + ρ kB T m is pressure due to the self-interaction and effective temperature.Taking into account the Poisson equation ( 2) and neglecting the quantum pressure term Q, we obtain the following equation of state where G is the gravitational constant.The solution of this equation is discussed in Sec.III.

B. Gravitoelectromagnetic approach
To determine the gravitational field of DM halo we employ the well-known gravitoelectromagnetism (GEM) approach [29] which was previously applied to galactic structures in [25,27,28].According to the GEM formalism, in the case of a test particle (which is luminous matter in our case) moving much slower than the speed of light c, it is convenient to represent the spacetime metric in the form where Φ g and A g are the GEM scalar (gravielectric) and vector (gravimagnetic) potentials.
For the gravitoelectromagnetic fields E g and B g the Einstein equations imply the following relations: Here sources of the gravitational field are mass density ρ and matter current j = ρu (u is the matter velocity).
Since these equations are clearly analogous to those in the electromagnetic theory, their solutions have a form similar to Maxwell's theory where integration proceeds over r ′ occupied by DM particles, ρ(r ′ ) is the condensate density and u(r ′ ) is the BEC velocity at r ′ .r are coordinates associated with the test particle, moving along geodesics in the BEC gravitational field.Finally, the geodesic movement for a test particle, which corresponds to the spacetime metric in the GEM form, ∂A g ∂x i dx dt can be equivalently described as the classical motion mẍ = F g in the gravitoelectromagnetic analog of the Lorentz force where v is the particle velocity and m is its mass.
Here we introduced gravielectric a E = −E g and gravimagnetic a B = − 2 c v × B g components of acceleration.

III. HALO DENSITY PROFILE
The model, based on the generalized GPP equations (see Eqs. ( 1), ( 2)) describes the coreenvelope structure of DM halo with a dense core and diffuse isothermal envelope.The model yields the following equation of state for the ULDM P = Kρ 2 + ρ kB T m (see Sec. II).Thus, one can conclude, that in the core region equation of state is approximately P = Kρ 2 , because the weak self-interaction dominates over effective temperature impact due to large density.That is why the latter will be neglected in the discussion of the core states.On the contrary, in the isothermal envelope region, we have the equation of state P = ρ kB T m , which means that the effective temperature term plays a crucial role there.
Based on these considerations, we calculate the halo density in two steps.Firstly, we reproduce the numerical result for the total density of the non-rotating halo (see the original result in [16]), which defines density distribution in the isothermal envelope region.This step is needed as a starting point to define isothermal enve-lope density distribution and to compare the discussed in [16] s = 0 solitonic core with the new case of vortex core s = 1.Secondly, under the assumption that core and envelope do not interact, we discuss the core density profile separately by means of variational ansatz [36].This way we will study the spherically symmetric (s = 0) and the single-charged vortex (s = 1) solutions for the core density distribution.

A. Isothermal envelope
In the first case of a non-rotating core, we can set u = 0 , and then the Eq. ( 3) simplifies where we took into account that K = 4πas 2 m 3 .It is convenient to introduce the density function and the radial coordinate ρ = ρ c e −f , y = r/r 0 , where and ρ c defines the density at the center.The equation of state can be rewritten in the following form: where χ = 4πa s 2 ρ c /(m 2 k B T ).The boundary conditions are f (0) = 0 and df dy (0) = 0, which define the boundary conditions for the density function ρ(0) = ρ c and dρ(0) dr = 0. We solve Eq. ( 11) numerically for different values of χ and present solutions in Fig. 1 (a).The isothermal envelope density distribution is defined as ρ = ρ 0 e −f = ρ 0 f N (r), where f is a numerical solution of Eq. (11).
The profile has a solitonic core and an isothermal envelope whose density decreases as ρ(r ∞ /(4πGr 2 ) [16] in agreement with observations (here v ∞ is the constant rotational velocity in the large distance limit).The existence of a BEC core in the ULDM model was also discussed in [17][18][19][20].
The possible physical origin of the coreenvelope structure could be the merger of twostate configurations when the total system tends to a virialized state, and the obtained averaged profile has a core and a tail structure [37].The process of halo formation usually undergoes gravitational cooling [38], which is discussed in [38,39].Gravitational cooling process for inreitially quite arbitrary density profiles leads to relaxation and virialization through the emission of scalar field particles [40].The resulting profile has the same dense core and diffuse envelope structure.
In the case s = 1, the hydrodynamical velocity u does not vanish in the inner region due to the existence of the vortex.The definition of the velocity profile in the isothermal halo region is a complicated task.One would expect that there is an intermediate region between the core and isothermal envelope, where the hydrodynamical velocity is small but nonzero, and at large enough distances, we should have u = 0.This is due to the divergent mass of the isothermal envelope, which therefore cannot rotate in order for kinetic energy to be finite.For an estimate, we simply put u = 0 in the whole isothermal envelope region.This approximation can be justified by the negligibly small density of the isothermal envelope in comparison with the core density, so its rotation would have no sufficient impact on the system.Hence the density profile in the envelope region remains unchanged.Thus, to define isothermal envelope density distribution we use the numerical solution for ρ = ρ 0 f N (r), obtained earlier in the case of non-rotating core.The density profile in the core region will be discussed in the next section in detail.
To reproduce the Milky Way halo mass M = 1.3 × 10 12 M ⊙ and radius R halo = 287kpc [41], taking into account the model described in [16], we choose the following values of the particle mass m = 2.92 × 10 −22 eV /c 2 = 0.52 × 10 −57 kg, scattering length a s = 8.17 and distance scaling parameter r 0 = 0.071kpc.Then χ = 20 and the temperature of BEC of such ultralight bosons is much larger than the effective temperature.For the spherically symmetric case, this yields the core with mass M c = 6.39 × 10 10 M ⊙ and radius R c = 1kpc.

B. Core stationary states
The dynamics of self-gravitating BEC of N weakly interacting bosons with mass m is described by the GPP system of equations with the term, corresponding to the effective temperature impact: where g = 4π 2 as m is the coupling strength that corresponds to the two-particle interaction, a s is the s-wave scattering length, Φ is the gravitational potential and G is gravitational constant.
The GPP system of Eqs.( 12) and ( 13) includes three crucial physical parameters: particle mass m, the total number of particles N (or, equivalently, total mass M ) and coupling strength g (or, equivalently, self-interaction constant λ 8π = as λc , where λ c = mc is the Compton wavelength of bosons) [36].
The GPP system of equations is invariant under the transformation t , where λ * > 0, which allows us to scale-out the coupling constant to g = 1.
In order to simplify calculations, it is convenient to introduce dimensionless variables and wave function where the dimensional variables are related to the dimensionless ones as follows: and will be neglected in the following discussion because the corresponding term T eff ln |ψ| is negligibly small in the core region.Therefore, we neglect the temperature effects in the analysis of the BEC core density distribution .
For the BEC core mass M c = 6.39 × 10 10 M ⊙ and radius R c = 1kpc, we solve the GPP equations ( 14) and ( 15) by using the variational ansatz in cylindrical coordinates r, z Here R and η are variational parameters, which will be fixed later.Constant A is fixed by the normalization condition the cases s = 0, 1 are considered, and N 0 is defined by the core mass where m Pl = c G is the Planck mass and λ/(8π) = 1.21×10 −91 is the self-interaction coupling constant.
The dimensionless quantities and the physically observed ones are related as follows: where R 99 is the dimensionless radius which contains 99 percent of the mass of the core (the vari-ational analysis gives R 99 ≈ 2.38R in the case of solitonic core and R 99 ≈ 2.58R in the case of vortex core), ρ is the condensate density, and ρ 0 = M A 2 /(L 3 N 0 ) is the density scaling parameter.R c denotes the total radius of the core in physical units.
Using the variational ansatz for the BEC wave function (16), we obtain the energy [36]: where Γ(x) denotes the Gamma function, Erfc(x) is the complementary error function and L s (x) is the Laguerre polynomials.Here ǫ = ( 2 /4πm Pl λ 2 c )(8π/λ) 3/2 is characteristic energy, which does not depend on variational parameters.
In what follows, we will use r 0 = 2.18 × 10 18 m = 0.071kpc = 0.22L as the distance scaling parameter.
In the subsection below, we investigate the case s = 0.

Non-rotating spherically-symmetric core
In this case, the BEC wave function in Eq. ( 16) depends only on radial distance r in spherical coordinates (21) and the density function (see Eq. ( 19)) equals In what follows, r will denote spherical distance, when the s = 0 case is discussed.
We should relate R and the BEC core radius R c which is defined through , the numerical result for halo density (see Fig. 1 a) gives Rc LR = 1.64 or R = 8.66 in the r 0 scale.It is interesting to compare the obtained R with its value in the variational analysis method used in [36].Substituting η = 1 and s = 0 in the energy functional in Eq. ( 20), we get Its extremum is defined by the equation that gives R = 1.73 or R = 7.86 in the r 0 scale.Thus, R c = 0.9 kpc (see Eq. ( 18)) and, therefore, the variational analysis method and numerical calculation (see Fig. 1 (a)) are in a good agreement.

Rotating axially-symmetric core
In the case s = 1 (see Eq.( 16)), we have a wave function, which depends on cylindrical coordinates r, z, φ and the density function equals where A is given by Eq. ( 17).The dimensionless total energy in Eq. ( 20) for s = 1 reads Equations of an extremum of the total energy with respect to η and R yield the solution η = 1.464, and R = 1.226 in the L scale.In the r 0 scale, we have R = 5.57.To determine the core density distribution, we use the variational analysis result.We assume that the core interacts negligibly weakly with the isothermal envelope.Therefore, for the isothermal envelope region, we use the numerical distribution f N (r sph ) = f N ( √ r 2 + z 2 ) (see Fig. 1 (a)), derived under u = 0 condition.Thus, we obtain (see Fig. 1 (b)) where r = x 2 + y 2 and z are cylindrical coordinates and r sph = x 2 + y 2 + z 2 .Here ρ 0 is the spherical halo central density.The spherically symmetric isothermal envelope density ρ(r, z) = ρ 0 f N (r sph /r 0 ) is found numerically by solving Eq. (11).The total core radius R c is defined by Eq. (18).
By using u = j ph ρ and the particle current we find the velocity distribution u(r) of DM particles where α = /(mr 0 c) = 0.31 • 10 −3 .Obviously, the velocity of condensate particles increases while approaching the center of the vortex.Note that there is an inner region where the velocity becomes of the order of c and, therefore, this region cannot be described by making use of the gravitoelectromagnetism ansatz (see Appendix A for explanation).This region is limited by the radial distance r = αr 0 = 2.2 × 10 −5 kpc.
In the two following sections, by using the for-malism of GEM, we describe particle movement in the gravitational field of DM in the s = 0, 1 states aiming to understand how baryonic matter particles interact with the proposed DM .

IV. GRAVIELECTRIC FIELD AND ROTATION CURVES
In this section, we obtain numerical results for the gravielectric (Newtonian) component of the DM halo gravitational field.Having calculated the field, we analyze the rotation curves, predicted by the model in the cases of soliton and vortex core.
To determine the gravielectric potential in the case of a non-rotating halo we use the numerically obtained density distribution (see Fig. 1  (a)).In the case of a rotating axially symmetric halo, the mass density distribution is shown in Fig. 1 (b).In the spherically symmetric case of nonrotating halo (s = 0), only the radial component of the gravielectric field is not zero (see Eq.( 7)) and the corresponding gravielectric acceleration a E = −E g = a E e r (see Eq.( 9)) is presented in Fig. 2. The acceleration at large distances behaves like a E /a 0 = 82.66r0 /r, i.e., a E = 9.3 × 10 −29 kpc 2 s 2 × 1/r.Here a 0 = Gρ 0 r 0 = 5.38 × 10 −13 km/s 2 .In the core region, where the density distribution is described by the variational ansatz (22), the gravielectric potential and the corresponding acceleration can be found analytically The general solution is given by where Erf(x) denotes the error function and c 1 , c 2 are constants.We can set c 2 = 0.At a large distance, the gravielectric potential of the halo must be equal to the potential of a body with the same mass M = π 3/2 ρ 0 R 3 .This implies that c 1 = 0. Thus, Φ(r) is completely determined and we have the radial acceleration  Clearly, a E has a maximum at r = R = 8.66 in the r 0 scale in agreement with the radial gravielectric acceleration shown in Fig. 2.
Gravielectric field in the case of vortex core has radial and z components in cylindrical coordinates, namely a E = a Er e r + a Ez e z .They are illustrated in Figs. 4 and 5, respectively.The radial dependence of the gravielectric radial acceleration in the z = 0 plane is shown in Fig. 3. Notice that at r ≈ 0.81r 0 = 0.058 kpc the acceleration projection changes sign, hence test particles are repelled in the interior region and attracted in the exterior region.This result stems from the geometry of the considered doughnut-shaped halo with a hole.Now we aim to determine the impact of the gravitational field of the DM halo on the movement of celestial bodies in the Milky Way galaxy.According to our model (see Sec. III) density distribution depends on the state of the core, which must lead to a difference between the rotation curves, which they induce.To demonstrate how the gravielectric acceleration induces rotation in the s = 0 and s = 1 cases, we present the rotation velocity v in the z = 0 plane as a function of the radial distance r in Fig. 6.The new result here is the curve in the case s = 1, while s = 0 case was discussed earlier in [16].The two halos with s = 0 and s = 1 core have equal mass, which is the observed mass of DM halo in the Milky Way, according to the model discussed in Sec.III.The numerical results indeed show that at large distances the corresponding rotational curves have the same asymptotic.Note that, the gravielectric force in s = 1 case changes its sign at r = 0.81r 0 = 0.058 kpc.Hence, at distances less than 0.058 kpc there are no stable rotation orbits in the rotating halo model.However, the stable orbits are possible if one includes not only DM but also the other sources of the gravitational field, namely, the baryonic galactic bulge and the supermassive black hole in the central region of the galaxy.

V. GRAVIMAGNETIC FIELD IN THE BEC CORE
In this section, we obtain numerical results for the gravimagnetic (first post-Newtonian) component of the DM halo gravitational field (see Subsec.II B of Sec.II).The component is induced by a moving source, hence, it is nonzero only in the second case of the DM halo with a vortex core.
To determine the gravimagnetic potential in the case of a rotating axially symmetric halo we use the mass density and velocity distributions given by Eqs. ( 26) and ( 27).The calculation is based on Eqs. ( 8) and ( 6).The results of numerical integration for radial and z-components of the gravimagnetic field, B = B r e r + B z e z , are shown in Figs.7 and 8, respectively.Fig. 9 displays the z-component of the gravimagnetic field B z in the z = 0 plane (the radial component of the gravimagnetic field equals zero in this plane).Having determined the gravimagnetic field, we can calculate the corresponding acceleration of the test particle.Using Eq.( 9), we have = a Br (a, b, k)e r + a Bz (a, b, k)e z .This allows us to estimate the impact of the gravimagnetic field on stars' motion.In the case of the Milky Way galaxy, v = v 0 + γ 0 r 0 a if a < a break and v = v 1 + γ 1 r 0 a for a ≥ a break [42].Constants γ 0 , γ 1 , a break , and v 1 are different for the thick and thin galactic disks' velocity profiles.Setting R break = r 0 a break = 5 kpc and v 0 =0 in both cases gives the values of parameters presented in Table I.This approximation is valid up to 13 kpc = 180r 0 [42].We should emphasize here that v includes only the component of velocity directed along e φ and does not include the component along e r .It is important to distinguish the φ-component and the absolute value of the whole velocity when dealing with sufficiently non-circular elliptic orbits.
According to Eq.( 9), the gravimagnetic acceleration in galactic plane c = 0 can be estimated as where i = 0 for a < a break and i = 1 for a ≥ a break .The corresponding plot is shown in Fig. 10.The spike on the red curve, which shows the modulus of the ratio of the gravimagnetic acceleration to the gravielectric one |a Br /a Er | appears because the gravielectric acceleration changes sign at r = 0.81r 0 = 0.058 kpc.It is interesting that a Br (a) tends to a constant in the a ≪ 1 limit (see Fig. 11).This directly follows from the analytical expression.In the interior region r < 5 kpc, we have a Br a 0 = −1.38αγ 0 r 0 a c B r (a, b, 0).
Table I.Parameters of the Milky Way's rotational velocity profiles [42].In the a ≪ 1 limit, we find The last integral can be calculated numerically which yields We see that in the case under consideration the gravimagnetic acceleration indeed tends to be a constant in the a ≪ 1 limit.
The gravimagnetic field calculations performed in this section allow us to obtain some testable predictions of the model.According to numerical results for B g and E g , the gravielectric force changes its sign at r = 0.81r 0 = 0.058 kpc, and the gravimagnetic force component is attractive or repulsive, depending on the direction of the motion.The acceleration in the polar coordinates (r, φ) is given by a = (r − r φ2 )e r + (r φ + 2 ṙ φ)e φ .Then the equations of motion for a star take the form Since the gravielectric acceleration dominates over the gravimagnetic one, it suffices to take the latter into account as a perturbation.Therefore, we treat B g as the first-order perturbation and expand φ(t) and r(t) around the solution r c and φ c determined by the gravielectric acceleration.For r = r c + δr and φ = φ c + δφ, in the zeroth-order, we have the Kepler problem equations with E r (r c ) calculated numerically in Sec.IV.The corresponding solutions are elliptic orbits.For simplicity, we will consider only the case of circular orbits r c (φ) = r c = const.By substituting E r (r c + δr) ≈ E(r c ) + dE dr (r c )δr in Eqs. ( 28) and ( 29), we obtain where w 0 = dφc dt is angular frequency, induced by gravielectric field.Thus, it can be explicitly written as w where we set the integration constant to zero.Substituting this relation in the first equation, we find From the numerical result, we see that f (E) is positive and tends to zero at a large distance.

VI. CONCLUSIONS
We investigated the model of DM halo with BEC core composed of ultra-light bosonic particles.Solving the generalized GPP equations for self-gravitating BEC we obtained the density profile of the DM halo and analyzed its core and envelope structure.The density and velocity profiles were found for two types of stable structures with topological charges (s = 0 and s = 1) of the BEC core.
Using this DM halo description, we investi-gated its gravitational field and the impact of this field on the baryonic matter.The key result of our paper is that the observable effects, predicted by the ULDM halo model, depend on the state of the core.In particular, solitonic and vortex cores yield different density and velocity distributions and thus different gravitational fields.The doughnut-like density distribution (vanishing at the vortex core) and vortex flows (rapidly increasing at the vortex axis) of the BEC core can significantly modify both gravielectric and gravimagnetic components of the gravitational field.We described the gravitational fields of these two core configurations by using the gravimagnetism approach.A dominant component of the gravitational field is the gravielectric (Newtonian) one, which generates the rotation of celestial bodies in the galaxy.The rotational velocity induced by the halo with vortex is smaller close to the core region but has the same asymptotics at large distances in comparison with the non-rotating halo.
The first post-Newtonian component of the gravitational field, which is called gravimagnetic, is induced by the rotation of the BEC vortex core and appears only in the model of a rotating halo.Although, as expected, the gravimagnetic acceleration is much weaker than the gravielectric one, it can affect the dynamics of baryonic matter in the halo, especially in its inner region.In our simplified perturbation approach for circular orbit gravimagnetic field yields radius and frequency shift, and can also induce trajectory oscillations, depending on initial conditions.There are several possible directions in which the present study could be extended.An analysis of gravitational fields beyond the gravimagnetic approach is required in the central region of the galaxy, due to the high rotational velocity of BEC there.Furthermore, according to astrophysical observations, there is a supermassive black hole in the center of our galaxy whose presence should be taken into account.Finally, the gravitational effects of baryonic matter should be included in further studies.

VII. ACKNOWLEDGMENTS
The authors are grateful to Yelyzaveta Nikolaieva, Sebastian Ulbricht, Stanislav Vilchinskii, and Luca Salasnich for useful discussions and comments.A.Y. acknowledge support from BIRD Project "Ultracold atoms in curved geometries" of the University of Padova.

APPENDIX A
Let us discuss the self-consistency of our model, which makes use of the GEM approach to describe the first post-Newtonian contribution to the gravitational field potential.We assumed that a test particle (celestial body acted upon by the gravitational field) propagates with a nonrelativistic speed v so that all terms of higher than linear order in O(v/c) can be neglected in the equations of motion.As to DM, we describe it by using the nonlinear Schrödinger equation with gravitational potential Φ g .
Since the hydrodynamical velocity in the vortex (the state with s = 1) is u(r) = αcr 0 /r, it increases at small r and attains at r ∼ αr 0 values of the order of c. Obviously, the Newtonian treatment is not applicable in this region.Therefore, we use the Klein-Gordon equation in order to describe the relativistic equation of motion of bosons, as follows: where U = 2m 2 gN |φ| 2 and φ is the scalar field.We neglect the effective temperature because only the core region is investigated (the hydrodynamical velocity u(r) is nonzero only in the core region) and ∇ α denotes covariant derivative in curved space-time.
The metric in the GEM approach reads (here all notations are the same as in Subsec.II B) and the Laplace operator is given by where g = det(g µν ) ≈ −1.Then we have where fields Φ g and A g are time-independent.
Taking into account the gauge condition ∂ i A i g = 0, we find To obtain a nonrelativistic approximation of the Klein-Gordon equation we represent the scalar field in the form φ = e imc 2 t/ ψ.Substituting this expression in the Klein-Gordon equation and multiplying by e −imc 2 t/ we get Neglecting terms of order of (u/c) 2 and higher (A g ∼ u/c), we obtain Finally, after some straightforward simplifications, we derive the Schrödinger equation in the Thus, we conclude that the model is self-consistent if we take into account only terms up to u/c, or, equivalently, in the region, where the hydrodynamical velocity of vortex is not relativistic (u ≪ c).

Figure 1 .
Figure 1.Halo density profile ρ/ρ0 as a function of dimensionless r/r0 coordinate in the plane z = 0, both x and y axes have log scale.The cyan insets in both plots show 3D density isosurfaces of the corresponding BEC cores.Left panel (a) shows the halo with the BEC core in a soliton state (s = 0).Three curves correspond to different values of parameter χ = 4πas 2 ρc/(m 2 kBT ), so that while increasing χ one decreases effective temperature T and vice versa.Right panel (b) shows the halo with the core in a vortex state (s = 1, χ = 20).Note, we investigate in detail the isothermal envelope for χ = 20, which is consistent with observations for the Milky Way.The black dashed line divides the distribution into two parts: the inner region with a rotating core and the outer region composed of an isothermal envelope.

Figure 2 .
Figure 2. The radial component of gravielectric field aE/a0 (blue dashed line) and density (red solid line) of the non-rotating halo (s = 0 core) as functions of the dimensionless r/r0 coordinate, both x and y axes have log scale.Here a0 = 5.38 × 10 −13 km s 2 , r0 = 71pc.

Figure 3 .
Figure 3.The radial component of gravielectric acceleration aEr/a0 (blue dashed line) and density (red line) of the rotating halo (s = 1 core) as functions of dimensionless r/r0 coordinate in the z = 0 plane, both x and y axes have log scale.Here a0 = 5.38 × 10 −13 km s 2 , r0 = 71pc.

Figure 4 .
Figure 4.The radial component of gravielectric acceleration aEr/a0 induced by the rotating halo (s = 1 core) as a function of dimensionless r/r0 and z/r0 coordinates.Here a0 = 5.38 × 10 −13 km s 2 , r0 = 71pc.The left panel shows the isothermal envelope region with the three axes in the log scale and the right panel is a zoom-in of the core region.

Figure 5 .
Figure 5.The z-component of gravielectric acceleration aEz/a0 induced by the rotating halo (s = 1 core) as a function of dimensionless r/r0 and z/r0 coordinates.Here a0 = 5.38 × 10 −13 km s 2 , r0 = 71pc.The left panel shows the isothermal envelope region with the three axes in the log scale and the right panel is a zoom-in of the core region.

Figure 6 .
Figure 6.The rotation (Kepler) velocity v in the z = 0 plane as a function of the radial distance r.The pink dashed line corresponds to the nonrotating spherical halo (s = 0 core) and the cyan solid line to the rotating halo (s = 1 core).The background represents a gradient plot of the density distribution in s = 1 case

Figure 10 .
Figure 10.The radial component of gravimagnetic acceleration (solid and dashed blue lines) aBr/a0 for thin and thick disks, respectively, and the absolute value of the ratio of gravimagnetic acceleration to gravielectric |aBr/aEr| (both for thin and for thick disks) as a function of r/r0.Here a0 = 5.38 × 10 −13 km s 2 , r0 = 71pc.