Axionic Dark Matter Halos in the Gravitational Field of Baryonic Matter
11 Axionic Dark Matter Halos in the Gravitational Field of Baryonic Matter
Gennady P. Berman , Vyacheslav N. Gorshkov , and Vladimir I. Tsifrinovich , Theoretical Division, T-4, Los Alamos National Laboratory, Los Alamos, NM 87545, USA National Technical University of Ukraine “Igor Sikorsky Kyiv Polytechnic Institute”, Kyiv, 03056, Ukraine Department of Applied Physics, NYU Tandon School of Engineering, Brooklyn, NY 11201, USA
LA-UR-20-24164
Abstract
We consider a dark matter halo (DMH) of a spherical galaxy as a Bose-Einstein condensate (BEC) of the ultra-light axions (ULA) interacting with the baryonic matter. In the mean-field (MF) limit, we have derived the integro-differential equation of the Hartree-Fock type for the spherically symmetrical wave function of the DMH component. This equation includes two independent dimensionless parameters: (i) - the ratio of baryon and axion total mases and (ii) - the ratio of characteristic baryon and axion spatial parameters. We extended our “dissipation algorithm” for studying numerically the ground state of the axion halo in the gravitational field produced by the baryonic component. We calculated the characteristic size, c x , of DMH as a function of and . and obtained an analytical approximation for c x . Introduction
Ultralight axions (ULA) have attracted much attention as a possible candidate for the dark matter halo (DMH) in galaxies (see, for example, [1-5]). Typical models suggest that axionic DMH contains a Bose-Einstein condensate (BEC) core and a quasi-classical axionic envelope (see, for example, [3,4]). In our previous work [5], we suggested that in some galaxies the whole DMH can be described by a simple model of the spherically symmetrical self-gravitating BEC completely ignoring the baryonic component. Certainly, it is possible only when the mass of the baryonic component is negligible compared to the mass of DMH. It is known that in some galaxies the mass of the baryonic matter is essential. Moreover, recent observations show that the mass of the baryonic matter may exceed the mass of the DMH [6]. In this work we suggest that our simple model of spherically symmetrical axionic BEC [5] can be applied to the description of the DMH of some galaxies with the spherically symmetrical distribution of the baryonic mass for an arbitrary relation between the masses of the baryonic and axionic components. We consider a spherically symmetrical BEC of self-gravitating axions which experiences the gravitational field of a given spherically symmetrical distribution of the baryonic matter. We note that a closely related problem of axionic BEC, which experiences the gravitational field of a black hole, was considered in ref. [4]. In Section II, we derive, in the mean-field (MF) limit, the non-linear integro-differential equation for the axionic BEC in the presence of the baryonic component. In Section III, we describe our dissipation method for numerical solution of this equation, and in Section IV, we describe the results of our numerical simulations. The main objective of our paper is to calculate the characteristic size of DMH as a function of two independent parameters, the ratio of masses of the baryonic and axionic components and the size of the baryonic components. Our results may be relevant for description of recently discovered galaxies with the dominant baryonic component.
II. Non-Linear Integro-Differential Equation for the axionic DMH in the presence of baryonic mass
We start with the non-relativistic Schrödinger equation, ,..., ; ,..., ;
N N N i r r t t H r r t , (2.1) for N axions which interact gravitationally between themselves and with the given distribution of the baryonic matter. The system is described by the Hamiltonian, i N NN a b i ari i i j Na i j
H m U r Gmm r r (2.2) In (2.2), b U r is a given baryonic potential, ' ',' bb b r d rU r GM r r (2.3) where a m and , b M are the mass of an axion and the total baryonic mass, correspondingly, and b r is the probability density for the baryonic matter. We use the normalization condition for , b r b r d r (2.4) For definiteness, we take the density of the baryonic matter for the spherical galaxies from ref. [7], and obtain the following normalized probability density, ,4 bb b rr r r r (2.5) where b r is a characteristic dimensional parameter. Using, the approach, described in our previous work [5], we derive a single-particle equation of the Hartree-Fock type, in the mean-field (MF) limit, ', '' ', , ,2 ' ' br a b a aa r t d rr d ri r t Gm M Gm M r tt m r r r r (2.6) Here a a M m N is the total mass of the axionic halo. We will use the following notations,
222 3 23 - characteristic size of the axionic halo with no baryons;- dimensionless time;- dimensionless coordinate;, , - dimensionless wave function;. a a aa aa ab a b r Gm Mtm rrx rx r r tx r r (2.7) In the dimensionless notations (2.7), Eq. (6) becomes, ', '' '1, , ,2 ' ' ,4 , 1, 1. bxb b x d xx d xi x xx x x xx x xx d x x d x (2.8) In (2.8), two independent dimensionless parameters are introduced, , ba MM , ba rr (2.9) By using the explicit expression for the baryonic density, the baryonic potential in (2.8) becomes, ' ' 1 ' ' ln 1 .' ' ' ' xbb x x d x dx dxU x x x x xx x x (2.10) Eigenvalue problem in the MF limit.
Below, we are interested in the spherically symmetric ground-state eigenfunction, , x of Eq. (2.8). The corresponding eigenvalue problem for Eq. (2.8) is described by the equation, ' '1 ln 1 ,2 ', xa a x d xx xx x xm r E (2.11) where E and are the dimensional and dimensionless, eigenenergies of an axion, correspondingly. Eq. (2.11) can be rewritten as, x x x x x x x s s ds x s s dsx x x x
0, 4 1. x x x dx (2.12) For numerical simulations, it is convenient to use the equation for, x x x (2.13) We obtain from (2.12) the equations which include the effective potentials, xx U x xx (2.14) x x sU x s ds dsx x s (2.15) x dx U x x dxx (2.16) x dx (2.17) III. Numerical Protocol
The numerical solution of the system of equations (2.14) - (2.17) is based on an iterative method for finding a stationary solution to the equation, , , ( ) , ,2 x t x ti i t U x t x tt x , (3.1) in the region, x X . The value of the upper boundary of the computational domain, X , guarantees a high accuracy of the solution to the problem, if in the integrals (2.14) - (2.17), it is used instead of infinity. One iteration cycle includes the following computational procedures. If the wave function, ( ) , k x t , is known at the end of the performed iteration cycle with the number, k , then, based on this function, the potential, ( ) ( ) k U x , and the energy, ( ) k , are redefined (see Eqs. (2.15) and (2.16)). The calculated ( ) ( ) k U x and ( ) k , are used in Eq. (3.1) for calculation of ( 1) ( ) , k k x t T at the end of the ( 1) k -th iteration cycle ( T is the iteration period). This procedure continues up to establishing of a stationary state. Now, we discuss the important factors that provide fast convergence of iterations. The calculations of ( 1) , k x t , with the initial condition, ( ) ( ) , , k k x t t x are performed using the implicit difference method of the first order of accuracy in time [8], j j j j j k k ji i i i i i i i U xt x , (3.2) where t and x are the discretization steps in time and space, correspondingly. We used the values, t and / , 4 10 x X N N . One iteration cycle assumes calculations of J steps in time ( J T J t ) in order to get the next approximation of the wave, ( 1) ( ) , k k x t T . In the general case, this method of solving the Schrödinger equation is imperfect, since it leads to an unphysical attenuation of the wave function in time. However, in this iterative method, it is this attenuation that ensures the convergence of the proposed iteration method. Let us explain this statement in detail. Represent the wave function, ( ) ( ) , , k k x t in the form of a series, ( ) ( ) ˆ, ( ) k k n n x t a x , (3.3) where ˆ ( ) n x are exact solutions of Eq. (2.14) for energy levels, n . The calculation of evolution of ( ) , k x t , according to the scheme (3.3), leads to the expression [9], ( ) [ ]( ) ˆ, exp( ) ( ) kn t tk n n n x t a i E t e x , / 2 n n E t , ( ) kn n E . (3.4) From (3.4), it is seen that the decomposition spectrum narrows over time and concentrates in the zone E , that covers the value, ( ) k . Therefore, the calculation of the next approximate value of the energy, ( 1) k , by the function, ( 1) ( ) , , k k x t x t T , on the ( 1) k -th iteration-cycle, will give a value closer to one of the levels n E than on the k -cycle. As a result of a sequence of iterations, the solution of equation (3.1) converges to one of the exact solutions of Eq. (2.14): , / 0 x t t in (3.1). To which particular level - it depends on the initial approximation, chosen on the basis of qualitative physical considerations. Obviously, to search for the ground state, an iterative procedure should begin with a function of the form, (0) 2 2 exp( / ) x x x . From this function, one should determine the potential, (0) ( ) U x . The energy value, (0) , should be set deliberately negative, since the relation (2.16) does not guarantee negativity of (0) when using the function, (0) x , in it. In subsequent iterations, when calculating ( ) k and ( ) ( ) k U x , it is preferable to use not the function ( ) k x , but its module, ( ) k x . If the solution is close to the exact one, then the real and imaginary parts of the wave function are similar to each other, so a replacing, ( ) ( ) k k x x , does not change the result. However, this technique accelerates the convergence of iterations to the exact solution. Since the difference scheme (3.2) introduces an attenuation, the wave function should be renormalized to unity at the beginning of each iteration cycle according to (2.17). When approaching the stationary state, the corresponding normalization factor approaches unity, and the sequence of energy values, ( ) k , also saturates. The control of these quantities determines the number of necessary iterations, iter N . In our problem, it is convenient to calculate the characteristics of the system on the plane ( , ) along the trajectories, const , by changing the value, . At the beginning of such trajectories ( , baryons are absent), the wave function of axions is easily calculated in [5], and then it is used as an initial approximation for subsequent growth of . Changing the parameter by 0.3 requires 7 iterations ( iter N ) to obtain a new solution (the error in the energy value, , if we use the solution for the previous value of as the initial approximation). In conclusion, we note that the use of the Crank and Nicolson difference scheme [9], for solving equation (3.1), does not lead to convergence of iterations, since such a scheme does not introduce attenuation into the wave function and does not narrow the expansion spectrum (3.4) from iteration to iteration. IV. Results
Our first objective was to calculate the characteristic dimensionless size, c x , of the axionic DMH as a function of two dimensionless parameters: the ratio, , of the total baryonic and axionic masses, and the size, , of the baryonic component (in terms of the parameter, a r ). In [5], where the baryonic component was ignored , the size of the DMH was estimated under the condition that the probability distribution of the axionic component, , P x x approaches zero. Namely, ( 8) 0.0091 ( 0). c P x P x
With no baryonic component, according to [5], c x . The same criterion, for calculation of , c x was used in this paper, when the baryonic component is present In Fig. 1, we show the dependence of c x on two parameters, and . As expected, the characteristic size of DMH drops with increase of the relative mass of the baryonic component, , and decrease of its size, . Fig. 1, also shows the eigenenergy of an axionic ground state, , as the function of the parameters, and . In Fig. 2, we show the MF potential energy, , U x and the probability density, , P x for three concrete values of and . We have obtained an approximate analytical estimate for the characteristic size of the DMH, , c x in terms of the corresponding size, ,0 c c x x with no baryonic component,
121 (1 2.5 )(1 ) c c x x . (4.1) Fig. 3, demonstrates the validity of our analytical estimate (4.1). The solid lines represent the numerical solutions. The circles represent the approximation (4.1). V. Conclusion
In this paper, we follow our previous suggestion in ref. [5] that the DMH in some galaxies can be described as a BEC of ULA. We consider the spherical galaxies with an arbitrary relation between the masses of the baryonic and axionic components which may be especially relevant in connection to recently discovered galaxies with dominant baryonic component [6]. To take into consideration the baryonic component, we extended both our mathematical model based on MF limit and our numerical code based on the iterative approach. Using these techniques, we calculated the characteristic size of the DMH as a function of the ratio of masses of the baryonic and axionic components and characteristic size of the baryonic component. We have found an accurate analytical estimate for the DMH size. The results obtained in this paper may be relevant for the description of the DMH in galaxies with significant or dominant baryonic components.
Acknowledgements
The work by G.P.B. was done at Los Alamos National Laboratory managed by Triad National Security, LLC, for the National Nuclear Security Administration of the U.S. Department of Energy under Contract No. 89233218CNA000001
References
1. W. Hu, R. Barkana, and A. Gruzinov. PRL, , 1158 (2000), arXiv:astro-ph/0003365. 2. C.G. Boehmer and T. Hartko. and J. Cosmol. Astropart. Phys., , 025 (2007), arXiv: 0705.4158 [astro-ph].
3. L. Hui, J.P. Ostriker, S. Tremaine, and E. Witten, Phys. Rev. D, , 043541 (2017). 4. P.H. Chavanis, Eur. Phys. J. Plus., , 352 (2019). 5. G. P. Berman, V. N. Gorshkov, V.I. Tsifrinovich, M. Merkli and X. Wang. Modern Physics Letters A, , 1950361 (2019). 6. Q. Guo, H. Hu, Z. Zheng, S. Liao, W. Du, S. Mao, L. Jiang, J. Wang, Y. Peng, L. Gao, J. Wang and H. Wu, Nature Astronomy, , 246 (2020). 7. M. Nagashima and Y. Yoshii, Mon. Not. R. Astron. Soc., , 509 (2003). 8. D. Potter, Computational Physics, J. Wiley & Sons, 1973, ISBN-10: 0471695556 9. J. Crank and P. Nicolson, Proc. Camb. Philos. Soc., , 50 (1947). Figures
Fig. 1:
The characteristic size, , c x (A) end energy, , (B) of the DMH as a function of two dimensionless parameters, and : , When baryonic component is absent c x . Fig. 2:
Left – MF potential energy of DM axions (blue curves) and baryons (red curves). Right – probability function of axions, ( ) . P x x
Upper insert – P x - probability function of axions for , when baryons are absent ( ; the depth of the MF potential well - U ). ( A ) - ( B ) - ( C ) - . Fig. 3:
The size of the DMH as a function of the ratio, , of masses of baryonic and axionic components, at three values of the characteristic dimensionless size, ,