Primordial Axion Stars and Galaxy Halos Formation

1. Cosmological Production of Axions

We assume the universe is described by $\Lambda$CDM model [35,36]. In this model, the universe comprises photons, neutrinos, ordinary (baryon) matter, cold dark matter, and dark energy. The latter is defined by the cosmological constant, $\Lambda$, and is responsible for the observable acceleration in the expansion of our universe. The axion DM can be described by a model with the action given by [11,37,38],

$$\begin{array}{c}\hfill S=\int \sqrt{-g}{d}^{4}x\left({\displaystyle \frac{1}{2}}{g}^{\mu \nu }{\varphi }_{,\mu }{\varphi }_{,\nu }-V\left(\varphi \right)\right).\end{array}$$

(3)

The corresponding equation of motion is the Klein-Gordon equation:

$$\square \varphi =-{\displaystyle \frac{\partial V\left(\varphi \right)}{\partial \varphi }}.$$

(4)

Qualitatively the potential $V\left(\varphi \right)$, generated by QCD instantons, may be written as [5,11]

$$\begin{array}{c}\hfill V\left(\varphi \right)={\left(m_a{f}_a\right)}^2\left(1-cos(\varphi /f_a)\right),\end{array}$$

(5)

where $m_a$ is the axion mass. The potential $V\left(\varphi \right)$ is the periodic function with a period defined by a shift-invariance, $\varphi \to \varphi +2\pi f_a$. There are N degenerate vacua, and the axion field has a range $0<\varphi <2\pi f_a$ [11].

1.1. Non-Relativistic Limit of the Klein-Gordon Equation

In the non-relativistic limit, the Klein-Gordon equation has the form of the Schrödinger equation with the effective potential describing the interacting forces in the system in the non-relativistic limit [20,21,24,39,40]. Writing,

$$\varphi ={\displaystyle \frac{1}{2\sqrt{m_a}}}\left(\psi e^{-i{m}_{a}t}+{\psi }^{*}{e}^{i{m}_{a}t}\right),$$

(6)

one can show that the non-relativistic theory is described by the effective action for the complex field $\psi$ (see Appendix for detail):

$$\begin{array}{c}\hfill S_{\mathrm{eff}}=\int d^{4}x\left({\displaystyle \frac{i}{2}}\left({\psi }^{*}\dot{\psi }-{\dot{\psi }}^{*}\psi \right)-{\displaystyle \frac{1}{2m_a}}\nabla {\psi }^{*}·\nabla \psi -V_{\mathrm{eff}}\left({\psi }^{*}\psi \right)-{\displaystyle \frac{1}{8\pi G}}\nabla \Phi ·\nabla \Phi -m_a{\psi }^{*}\psi \Phi \right),\end{array}$$

(7)

where $\Phi \left(\mathbf{r}\right)$ is the Newtonian gravitational potential, and $V_{eff}\left(\right|\psi \left|\right)$ is the effective potential obtained as the non-relativistic limit of the instanton potential Eq. (5). The equations of motion are given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & i{\displaystyle \frac{\partial \psi }{\partial t}}=-{\displaystyle \frac{1}{2m}}{\nabla }^2\psi +m_a\Phi \psi +{\displaystyle \frac{\partial V_{\mathrm{eff}}}{\partial {\psi }^{*}}}\psi ,\hfill \\ \hfill & {\nabla }^2\Phi =4\pi G{m}_a{\left|\psi \right|}^2.\hfill \end{array}\end{array}$$

(8)

For the axion field, the effective instanton non-relativistic potential potential was derived in Refs. [39,41] and takes the form

$$\begin{array}{c}\hfill V_{\mathrm{eff}}\left({\psi }^{*}\psi \right)={\displaystyle \frac{1}{2}}{m}_a{\psi }^{*}\psi +m_a^2{f}_a^2\left[1-J_0\left(n^{1/2}\right)\right],\end{array}$$

(9)

where $J_0\left(z\right)$ is a Bessel function, $n={\psi }^{*}\psi /n_0$ and $n_0=m_a{f}_a^2/2$.

Figure 1 represents the normalized effective potential, $V_0\left(x\right)=1-J_0\left(x\right)$, with the first term subtracted:

$$\begin{array}{c}\hfill V_0={\displaystyle \frac{1}{2m{n}_0}}\left({\mathcal{V}}_{\mathrm{eff}}-{\displaystyle \frac{1}{2}}{m}_a{\psi }^{*}\psi \right)\end{array}$$

(10)

The stable solutions for the axion stars with $n\gg 1$ (dense stars) correspond to the minima of the effective potential [20].

1.2. Formation of Axion Stars

Axion stars are characterized by their high density and stability, maintained by a balance between gravitational attraction and quantum pressure due to the axions’ wave-like nature. They can be categorized into two types:

• Dilute Axion Stars ($n_{core}\ll 1$). These stars are less dense and more extended, and insert "their size is determined by "their size is determined by the axion mass and self-interaction strength.
• Dense Axion Stars ($n_{core}\gtrsim 1$). These are more compact and can undergo complex dynamical processes, including potential collapse into black holes and phase transition to the dilute stars under certain conditions.

The stability and structure of axion stars are determined by the balance between gravitational attraction and the quantum pressure arising from the Heisenberg uncertainty principle.

2. Model

The non-relativistic axions can form a Bose-Einstein condensate (BEC). The equations that describe a static BEC of axions can be obtained by considering an axion field $\psi (\mathbf{r},t)$ whose time dependence is $\psi \left(\mathbf{r}\right)exp(-i\mu t)$, where $\mu$ is the chemical potential [20]. They can be expressed conveniently as coupled equations for $\psi \left(\mathbf{r}\right)$ and the gravitational potential $\Phi \left(\mathbf{r}\right)$:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\nabla }^2\widehat{\psi }=& 2m_a^2\left(\Phi -{\displaystyle \frac{\mu }{m_a}}-{\displaystyle \frac{1}{2}}+{\displaystyle \frac{J_1\left(\right|\widehat{\psi }\left|\right)}{|\widehat{\psi }|}}\right)\widehat{\psi },\hfill \\ \hfill {\nabla }^2\Phi =& 4\pi G{m}_a{n}_0{\left|\widehat{\psi }\right|}^2,\hfill \end{array}\end{array}$$

(11)

where $\widehat{\psi }=\psi /\sqrt{n_0}$, $G=1/m_p^2$, and $m_p$ denotes the Planck mass ($m_p=1.2·{10}^{19}\mathrm{GeV}$). To further simplify matters, we rescale coordinates and functions: $r\to r^{*}=r/r_0$, $\psi \to \psi \left(r^{*}\right)$ and $\Phi \to \gamma \Phi \left(r^{*}\right)$, where $r_0=\gamma /m_a$ and $\gamma =m_p/f_a$.

Hereafter, we drop *’s in rescaled functions and parameters. Then, the dimensionless equations can be recast as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\nabla }^2\widehat{\psi }=& 2\left(\Phi -\gamma \epsilon \right)\widehat{\psi },\hfill \\ \hfill {\nabla }^2\Phi =& 4\pi |\widehat{\psi }{|}^2,\hfill \end{array}\end{array}$$

(12)

where

$$\begin{array}{c}\hfill \epsilon =\overline{\mu }+{\displaystyle \frac{1}{2}}-{\displaystyle \frac{J_1\left(\right|\widehat{\psi }\left|\right)}{|\widehat{\psi }|}},\end{array}$$

(13)

and we set $\overline{\mu }=\mu /m_a$.

In Table 1, we present various length and mass scales to facilitate recalculation in physical units. In particular, the mass of the axion star can be obtained using the formula:

$$\begin{array}{c}\hfill M=3M_0\int {\left|\widehat{\psi }\right|}^2{d}^3{V}^{*},\end{array}$$

(14)

where the integration is performed over dimensionless variables ${\mathbf{r}}^{*}$, and $M_0=4\pi n_0{r}_0^3/3$.

Our model is based on the previous studies of the rotating axion stars made in [20,22,23,40,42]. In general $\psi$ may consist of a superposition of various spherical harmonics $Y_l^m$,

$$\begin{array}{c}\hfill \widehat{\psi }(r,\theta ,\phi )=\sum _{\ell =0}^{\infty }\sum _{m=-l}^{\ell }{f}_{\ell }^m{r}^{\ell }{Y}_{\ell }^m(\theta ,\phi ).\end{array}$$

(15)

For simplicity, we assume the following ansatz:

$$\widehat{\psi }(r,\theta ,\varphi )=R\left(r\right)Y_l^m(\theta ,\phi ).$$

(16)

After substitution of $\widehat{\psi }(r,\theta ,\varphi )$, one can rewrite Eq. (12) as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & {\Delta }_{r}R-{\displaystyle \frac{l(l+1)}{r^2}}R=2\left(\Phi -\gamma \overline{\mu }-{\displaystyle \frac{\gamma }{2}}+{\displaystyle \frac{\gamma J_1\left(R\right|Y_l^m\left|\right)}{R|Y_l^m|}}\right)R,\hfill \\ \hfill & {\nabla }^2\Phi =4\pi R^2{\left|Y_l^m\right|}^2,\hfill \end{array}\end{array}$$

(17)

where ${\Delta }_r$ denotes the radial part of the Laplacian,

$$\begin{array}{c}\hfill {\Delta }_r={\displaystyle \frac{1}{r^2}}{\displaystyle \frac{\partial }{\partial r}}\left(r^2{\displaystyle \frac{\partial }{\partial r}}\right).\end{array}$$

(18)

Note that (17) still depends on angular variables $\theta$ and $\phi$. Therefore, following [40], we take the angular average of these equations. rgb]1,0,0We define the angular averaged gravitational potential and use the normalization of the spherical harmonics:

$$\begin{array}{cc}\hfill & \tilde{\Phi }={\displaystyle \frac{1}{4\pi }}\int \Phi (r,\theta ,\phi )d\Omega ,\int Y_{\ell }^m{Y}_{{\ell }^{\prime }}^{m^{\prime }*}d\Omega ={\delta }_{\ell {\ell }^{\prime }}{\delta }_{m{m}^{\prime }}.,\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill & {\mathcal{J}}_l^m={\displaystyle \frac{1}{4\pi }}\int {\displaystyle \frac{J_1\left(R\right|Y_l^m\left|\right)}{R|Y_l^m|}}d\Omega .\hfill \end{array}$$

(20)

Then, we obtain for angular averaged potentials the following equations of motion:

$$\begin{array}{cc}\hfill {\Delta }_{r}R=& 2\left(\tilde{\Phi }+{\displaystyle \frac{l(l+1)}{2r^2}}-\gamma \overline{\mu }-{\displaystyle \frac{\gamma }{2}}+\gamma {\mathcal{J}}_l^m)\right)R,\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill {\Delta }_r\tilde{\Phi }=& R^2.\hfill \end{array}$$

(22)

For $l=0,1$, we have obtained the analytical expressions for ${\mathcal{J}}_l^m$. The computation yields:

$$\begin{array}{cc}\hfill & {\mathcal{J}}_1^0=J_0\left(a_{0}R\right)-{\displaystyle \frac{J_1\left(a_{0}R\right)}{a_{0}R}}+{\displaystyle \frac{\pi }{2}}\left(J_1\left(a_{0}R\right)H_0\left(a_{0}R\right)-J_0\left(a_{0}R\right)H_1\left(a_{0}R\right)\right),\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill & {\mathcal{J}}_1^{\pm 1}={\displaystyle \frac{1-cos\left(a_{1}R\right)}{a_1^2{R}^2}},{\mathcal{J}}_0^0={\displaystyle \frac{J_1\left(R\right)}{R}},\hfill \end{array}$$

(24)

where $a_0=\sqrt{3/4\pi }$, $a_1=\sqrt{3/8\pi }$ and $H_n\left(z\right)$ denotes the Struve function.

For a given function $R\left(r\right)$, the gravitational potential can be obtained due to the integration of Eq. (). The solution can be written as follows:

$$\begin{array}{c}\hfill \tilde{\Phi }\left(r\right)=-{\displaystyle \frac{1}{r}}{\int }_0^r{x}^2{R}^2\left(x\right)dx-{\int }_r^{r_h}x{R}^2\left(x\right)dx,\end{array}$$

(25)

where $r_h$ is the radius of the star. From here it follows that $\tilde{\Phi }\left(0\right)={\tilde{\Phi }}_0$ and ${\Phi }^{\prime }\left(0\right)=0$, where

$$\begin{array}{c}\hfill {\tilde{\Phi }}_0=-{\int }_0^{r_h}x{R}^2\left(x\right)dx.\end{array}$$

(26)

We impose the following boundary conditions to solve numerically the coupling system of the differential equations (21) - ():

$$\begin{array}{c}\hfill R\left(0\right)=R_0,R{\left(0\right)}^{\prime }=R_1,\tilde{\Phi }\left(0\right)={\Phi }_0,{\tilde{\Phi }}^{\prime }\left(0\right)=0.\end{array}$$

(27)

The radius of the axion star is obtained as the solution of the equation:

$$\begin{array}{c}\hfill \tilde{\Phi }\left(r_h\right)=-{\displaystyle \frac{M^{*}\left(r_h\right)}{r_h}},\end{array}$$

(28)

where $M^{*}\left(r_h\right)={\int }_0^{r_h}{x}^2{R}^2\left(x\right)dx$ is the scaled mass of the star, $M^{*}=M/M_0$, and

$$\begin{array}{c}\hfill \tilde{\Phi }\left(r_h\right)=-{\displaystyle \frac{1}{r_h}}{\int }_0^{r_h}{x}^2{R}^2\left(x\right)dx.\end{array}$$

(29)

The chemical potential is subject to the boundary conditions:

$$\begin{array}{c}\hfill \overline{\mu }=\left[{\displaystyle \frac{1}{\gamma }}\left(\tilde{\Phi }+{\displaystyle \frac{l(l+1)}{2r^2}}\right)+{\mathcal{J}}_l^m-{\displaystyle \frac{1}{2}}\right]{|}_{r=r_h}=0\end{array}$$

(30)

3. Results and Discussion

Following the procedure outlined in Sect. Section 2 we have solved numerically the equations of motion (19) - () for the axion mass, $m_a=5.3·{10}^{-}19$ eV, and the boundary conditions presented below.

• Dense stars:

  $$\begin{array}{c}\hfill R\left(0\right)=6,R{\left(0\right)}^{\prime }=0,\Phi \left(0\right)=-1.5,{\Phi }^{\prime }\left(0\right)=0,(l=0,m=0).\end{array}$$
• Dilute stars:

  $$\begin{array}{cc}\hfill & R\left(0\right)=0.1,R{\left(0\right)}^{\prime }=0,\Phi \left(0\right)=-0.065,{\Phi }^{\prime }\left(0\right)=0,(l=0,m=0),\hfill \\ \hfill & R\left(0\right)=0,R{\left(0\right)}^{\prime }=10,\Phi \left(0\right)=-0.001,{\Phi }^{\prime }\left(0\right)=0,(l=1,m=0),\hfill \\ \hfill & R\left(0\right)=0,R{\left(0\right)}^{\prime }=0.12,\Phi \left(0\right)=-0.5,{\Phi }^{\prime }\left(0\right)=0,(l=1,m=\pm 1).\hfill \end{array}$$

The results are summarized in Figure 2 and Table 2. To obtain the redshift, we use the relation $z=n_h/n_{DM}$, where $n_{DM}={\overline{\rho }}_{\mathrm{DM}}/m_a$ is the average density of DM obtained from the data presented in Refs. [39,43]: ${\overline{\rho }}_{\mathrm{DM}}=1.3\times {10}^{-6}\mathrm{GeV}/{\mathrm{cm}}^3$. As follows from Table 2, the formation of the axion stars begins in the radiation-dominated era, which, according to the $\Lambda$CDM model, corresponds to redshift $z\approx (5.244·{10}^3-1.47·{10}^{11}$).

The mass of galaxies can vary significantly depending on their type and size:

• Dwarf Galaxies: $M\approx ({10}^7-{10}^9)M_{\odot }$.
• Spiral Galaxies: $M\approx ({10}^{10}-{10}^{12})M_{\odot }$.
• Elliptical Galaxies: $M\approx ({10}^{12}-{10}^{13})M_{\odot }$.

Considering the distribution of different types of galaxies, the average mass of galaxies in our universe can be roughly estimated as follows: $M\approx {10}^{11}{M}_{\odot }$. As one can see, the mass of axion stars presented in Table 2 is inside the range of the observable data.

The current size of the axion stars can obtained using the relation: $d=2z{r}_h$. We obtain

• $d=52.2$ kpc ($M=1.57·{10}^{10}{M}_{\odot }$)
• $d=76.4$ kpc ($M=3.43·{10}^9{M}_{\odot }$)
• $d=2.52$ kpc ($M=4·{10}^7{M}_{\odot }$)
• $d=8.3$ kpc ($M=3.92·{10}^{10}{M}_{\odot }$)

These results are in agreement with an estimation of the halo size for typical galaxies:

• Spiral Galaxies: 61.3 to 92 kpc in diameter
• Elliptical Galaxies: 92 to 306.6 kpc in diameter
• Dwarf Galaxies: up to 9.2 kpc in diameter

Primordial ULA stars offer potential insights into the nature of dark matter and its role in galaxy formation. As axion stars accrete additional axions, they contribute to the overall density of dark matter halos. This increased density can enhance the gravitational potential of the halo, making it more effective at attracting and retaining baryonic matter. By acting as seeds for structure formation, primordial axion stars could play a pivotal role in the formation of galaxies and shape the universe as we observe it today. Further research is needed to understand their properties, but their potential impact on our understanding of the universe is profound.