Diversity of Fuzzy Dark Matter Solitons

The existence of dark matter (DM) can be inferred from the rotation curves of galaxies Rubin:1982kyu , the evolution of large scale structure Davis:1985rj and the gravitational lensing observations Clowe:2006eq . According to the free streaming length of DM, DM can be divided into cold, warm and hot categories. Although the standard Lambda cold DM ($\Lambda$CDM) cosmological model is very successful and the latest cosmic microwave background (CMB) observations Planck:2018vyg suggests that CDM accounts for about $26\%$ of today’s energy density in the Universe, we don’t know what particle or object CDM is. Especially, as one of the most promising candidates for CDM particle, the weak interacting massive particles (WIMPs) grounded on supersymmetric theories of particle physics still have not been detected PandaX-II:2016vec ; LUX:2015abn ; ATLASCMS ; AMS:2014bun ; Fermi-LAT:2011baq . Moreover, primordial black holes (BHs) can also serve as CDM Carr:2016drx but they still have not been identified. These null results coupled with the failure of CDM particles on sub-galactic scales Primack:2009jr imply that DM may not be cold.

There is a promising alternative to CDM which is an ultralight scalar field with spin-$0$, extraordinarily light mass ($\sim 10^{-22}\rm{eV}/c^{2}$) and de Broglie wavelength comparable to few kpc, namely fuzzy DM (FDM) Hu:2000ke . Due to the large occupation numbers in galactic halos, FDM behaves as a classical field obeying the coupled Schrödinger–Poisson (SP) system of equations

where $m$ is the mass of FDM which is described by the wavefunction $\Psi$, the gravitational potential $\Phi$ is sourced by the FDM density $\rho=|\Psi|^{2}$. The SP system also can govern the evolution of FDM in an expanding Universe Schive:2014dra ; Schive:2014hza ; DeMartino:2017qsa ; Mocz:2019pyf . It is worth noting that the SP system is just the weak field limit of the general relativistic Einstein–Klein–Gordon (EKG) system Kaup:1968zz ; Ruffini:1969qy ; Ma:2023vfa and the SP system should be replaced by the Gross-Pitaevskii-Poisson (GPP) system when the self-interactions between FDM particles exist Chavanis:2011zi ; Chavanis:2011zm .

For simplicity, in this parer, we confine ourself to the SP system in a non-expanding Universe. The evolution of this system is usually simulated numerically, including the formation, the perturbation, the interference/collision and the tidal disruption/deformation of FDM solitons Guzman:2004wj ; Paredes:2015wga ; Edwards:2018ccc ; Munive-Villa:2022nsr . If one just care about stationary solutions, one can turn to the shooting method to find the eigenvalues of equilibrium configurations Guzman:2004wj ; Davies:2019wgi .

Since FDM solitons are only dependent on the mass of FDM, the SP system has the scaling symmetry. As a result, the Universe should be very monotonous and full of similar structures. But in fact there is a great diversity of structures in the Universe. Therefore, the real astrophysics systems must follow the different variants of the exact SP system, even though the modified gravity and the different interactions between FDM are not considered. In this paper, we will consider several variants due to different extra terms, such as the gravitoelectric field, the gravitomagnetic field, an extra FDM soliton and the baryon profile. Our motivations are as follows. First, as pointed out in Guzman:2004wj , the SP system is the weak field limit of the EKG system. Therefore, there must be a variant system between these two ends where the gravitoelectromagnetism is a good approximation. While the gravitoelectric field or the gravitational potential of a supermassive BH has been discussed Edwards:2018ccc ; Munive-Villa:2022nsr ; Davies:2019wgi , the gravitomagnetic field due to the system’s angular momentum has not been considered before. Secondly, although the interference/collision with an extra FDM soliton also has been studied Paredes:2015wga ; Edwards:2018ccc ; Munive-Villa:2022nsr , they just dealt with the cases where the density-ratio of two solitons is $\sim 1$. For the extreme cases where the density-ratio of the denser solitons to another one is $\gtrsim 10^{4}$, numerical simulations need a larger simulation box to contain the larger soliton with lower density or need a higher resolution to depict the smaller soliton with higher density. However, a simulation with a larger simulation box and higher resolution at the same time is prohibitively expensive and almost impossible. Finally, a non-spherical baryon profile is located in galaxies which results in non-spherical FDM solitons, for example the impact of a baryon profile with the cylindrical symmetry and the parity symmetry has been quantified Bar:2019bqz . In fact, the baryon profile in the Milky Way may be ellipsoidal Iocco:2015xga ; Lin:2019yux .

This paper is organized as follows. In section II, we briefly review the effects of the gravitoelectric field due to a supermassive BH on the FDM solitons. In section III, we prove the effects of the gravitomagnetic field due to the rotation velocity of FDM, the supermassive black hole spin and the orbital motion of supermassive BH binary on the FDM solitons are negligible. In section IV, we study the interaction between soliton binary with the extreme density-ratio. In section V, we calculate the FDM solitons in a given ellipsoidal baryon profile. Finally, a brief summary and discussions are included in section VI.

The effects of the gravitoelectric field includes the tidal disruption/deformation of FDM solitons by a nearby supermassive BH Edwards:2018ccc ; Munive-Villa:2022nsr and the modified formation of FDM solitons by a central supermassive BH Davies:2019wgi . While the former non-spherical impacts are more complicated and have to be studied by numerical simulations, the latter spherical ones are simple and can be studied by the shooting method. In this section, we briefly review the latter cases.

Firstly, the gravitoelectric field $\mathbf{E}_{\rm g}$ or the gravitational potential $\Phi_{\rm e}$ sourced by a central supermassive BH contributes to a variant of the exact SP system as

where the supermassive BH with mass $M_{\rm bh}$ has a point mass potential

When the system in question features the spherical symmetry, the ansatz of $\Psi(r,t)=e^{-i\gamma t/\hbar}\psi(r)$ means the FDM soliton density is $\rho(r)=|\Psi|^{2}=\psi^{2}(r)$ and the FDM soliton mass is $M=\int_{0}^{\infty}4\pi r^{2}\rho(r)dr$. After defining a number of dimensionless variables as

the dimensionless version of Eq. (2) is

Fulfilling the arbitrary normalization $\tilde{\psi}(\tilde{r}=0)=1$ and $\tilde{\Phi}(\tilde{r}=0)=0$¹¹1Although the normalization $\tilde{\Phi}(\tilde{r}=0)=0$ for $\Phi$ is different from the normalization $\Phi_{\rm e}(r=\infty)=0$ for $\Phi_{\rm e}$ (Eq. (3)), this difference does not affect the solutions with the normalization $\tilde{\psi}(\tilde{r}=0)=1$., the boundary conditions $\tilde{\psi}(\tilde{r}=\infty)=0$, $\frac{\partial\tilde{\psi}}{\partial\tilde{r}}|_{\tilde{r}=0}=0$ and $\frac{\partial\tilde{\Phi}}{\partial\tilde{r}}|_{\tilde{r}=0}=0$, choosing the supermassive BH mass $\tilde{M}_{\rm bh}=\{0.0,0.5,1.0,1.5,2.0,2,5\}$ as the input parameter and adjusting the quantized eigenvalue $\tilde{\gamma}$, we can calculate the equilibrium configurations from Eq. (10) by the shooting method. Given the input parameter $\tilde{M}_{\rm bh}$, only the solution from the smallest $\tilde{\gamma}$ is stable and the ground state. In Tab. 1, we list the input values of $\tilde{M}_{\rm bh}$ and the corresponding eigenvalue $\tilde{\gamma}$ and soliton mass $\tilde{M}$ of the ground state solutions. In Fig. 1, the corresponding soliton profiles are plotted. We find that the larger $\tilde{M}_{\rm bh}$ leads to a more denser and compact soliton, hence a smaller dimensionless soliton mass $\tilde{M}$ due to the fixed normalization $\tilde{\psi}(\tilde{r}=0)=1$. As for its physical mass $M$, in fact it increases with $\tilde{M}_{\rm bh}$.

The SP system is the weak field limit of the general relativistic EKG system Guzman:2004wj . Therefore, there must be a variant system between these two extreme ends by taking some general relativistic corrections into consideration, for example the gravitomagnetic field $\mathbf{B}_{\rm g}$. Then there is a variant of the exact SP system as

where the gravitational potential $\Phi_{\rm m}$ is related to the gravitomagnetic field $\mathbf{B}_{\rm g}$ and the angular momentum of FDM particle $\hat{\mathbf{L}}$ as

Generally speaking, the linear term $\Phi_{\rm m1}$ is much larger than the quadratic one $\Phi_{\rm m2}$. In this section, we will talk about some possible sources of $\mathbf{B}_{\rm g}$ and their effects on the FDM solitons.

For the cases where the density-ratio of two solitons is $\sim 1$, this system obviously does not feature the spherical symmetry and can be studied by numerical simulations Paredes:2015wga ; Edwards:2018ccc ; Munive-Villa:2022nsr . On the contrary, for the extreme case where the density-ratio of two solitons is $\gtrsim 10^{4}$, there are two stages when this system do feature an approximate and local spherical symmetry: stage 1) the smaller soliton with higher density is just immerse in the boundary of the larger soliton with lower density which is relatively flat and can be considered as a background; stage 2) the larger soliton with lower density shares the same center with the smaller soliton with higher density which are more stable and can be considered as a background. The collisional dynamics include these two stages and we can still use the shooting method to deal with them. Finally, we can use the exact SP system (Eq. (1)) to describe the background soliton and the following variant of the exact SP system to describe the changed soliton during the above two stages

whose dimensionless version is

where $\rho_{\rm bg}$ is the density of the background soliton and its dimensionless counterpart is

Since this variant SP system (Eq. (31)) has the scaling symmetry

we can set up the background soliton with $\lambda=10^{-2}$ for stage 1 and with $\lambda=10^{2}$ for stage 2. Given the background density $\rho_{\rm bg}$, the changed soliton can be found by solving Eq. (31) with the shooting method. By assuming that its density profile $\tilde{\psi}^{2}(\tilde{r})$ is changing but its dimensionless mass $\tilde{M}$ is fixed as the initial value listed in the first row of Tab. 1, we get the eigenvalue $\tilde{\gamma}=0.6495$ for stage 1 and the eigenvalue $\tilde{\gamma}=69.2491$ for stage 2. For stage 1, the density profile $\tilde{\psi}^{2}(\tilde{r})$ of each soliton is shown in the left subplot of Fig. 5; for stage 2, each soliton’s $\tilde{\psi}^{2}(\tilde{r})$ is shown in the right subplot of Fig. 5.

We can find that the smaller soliton with higher density is almost not affected by the background soliton for stage 1 and the larger soliton with lower density shrink dramatically for stage 2. After stage 2, there are two comparable solitons and the system becomes complicated whose further evolution should be studied by numerical simulation.

Instead of a background soliton shown in section IV, we consider a background baryon profile in this section. More precisely, we consider an ellipsoidal baryon profile similar to the Milk Way’s Iocco:2015xga ; Lin:2019yux

where $\rho_{\rm bulge}$ is a triaxial and bar-shaped bulge of stars at the inner few kpc of the Milky Way Lopez-Corredoira:2005anz ; Ryu:2008dx

the disk density $\rho_{\rm disk}\approx 0.05M_{\odot}{\rm{pc}}^{-3}$ Lopez-Corredoira:2005anz can be neglected but $\rho_{\rm gas}$ is a non-negligible part of the baryons in the Milky Way in the form of molecular, atomic and ionised hydrogen and heavier elements. Here we assume that $\rho_{\rm gas}\approx m_{H_{2}}\left\langle n_{H_{2}}\right\rangle$, where $m_{H_{2}}$ is the mass of molecular hydrogen and $\left\langle n_{H_{2}}\right\rangle$ is the number density of molecular hydrogen in central molecular zone Ferriere:2007yq ,

For this non-spherical system, the shoot method is failed and numerical simulations should be performed. As an approximation, however, we consider three background baryon profiles respectively: $\rho_{\rm bg}^{x}(r)=\rho_{\rm bg}(x=r,y=0,z=0)$, $\rho_{\rm bg}^{y}(r)=\rho_{\rm bg}(x=0,y=r,z=0)$ and $\rho_{\rm bg}^{z}(r)=\rho_{\rm bg}(x=0,y=0,z=r)$. Then the shoot method still works for $\rho_{\rm bg}^{i}(r)$. If the results by different $\rho_{\rm bg}^{i}(r)$ are very different, it means the original $\rho_{\rm bg}(x,y,z)$ is too non-spherical to be described by the approximation $\rho_{\rm bg}^{i}(r)$.

To find the equilibrium configurations, as shown in section II, we usually use the arbitrary normalization $\tilde{\psi}(\tilde{r}=0)=1$. It means that the normalized soliton mass $\tilde{M}=\int_{0}^{\infty}\tilde{r}^{2}\tilde{\psi}^{2}d\tilde{r}$ is deviate from the true dimensionless mass by a scaling factor $\lambda$. For the Milk Way, the soiliton physical mass $M\approx 1.25\times 10^{9}M_{\odot}$ is obtained with the soliton-halo mass relation Schive:2014hza . On the other hand, taking this scaling factor into consideration, $\tilde{\rho}_{\rm bg}^{i}\rightarrow\lambda^{2}\tilde{\rho}_{\rm bg}^{i}$ means the background baryon profile is modified. To fix both of the soltion mass and the background baryon profile as the Milk Way’s, we introduce a scaling factor in advance and obtain a new dimensionless variant as

where $\lambda^{-2}_{0}$ is a parameter and not changing with the scaling factor $\lambda$. Therefore, this system still has the scaling symmetry with $\lambda^{-2}_{0}\tilde{\rho}_{\rm bg}^{i}\rightarrow\lambda^{2}\lambda^{-2}_{0}\tilde{\rho}_{\rm bg}^{i}$. Given an initial $\lambda_{0}$, we can obtain the normalized soliton profile $\tilde{\psi}^{2}$. And then we can obtain the scaling factor $\lambda$ by requiring the soliton mass to be the Milk Way’s. If $\lambda/\lambda_{0}\neq 1$, we replace the initial $\lambda_{0}$ with the newest $\lambda$ and repeat the former two steps. When $\lambda/\lambda_{0}-1<0.0001$, we stop the iteration.

In Tab. 2, Tab. 3 and Tab. 4, we list the eigenvalues $\tilde{\gamma}$, the soliton mass $\tilde{M}$ of the ground state solutions and the final values of $\lambda_{0}$ for the different background baryon profile $\tilde{\rho}_{\rm bg}^{i}$ and the mass of FDM $m$ respectively. Meanwhile, we plot the normalized and the physical density profiles of the ground state solutions with different background baryon profile $\tilde{\rho}_{\rm bg}^{i}$ and the mass of FDM $m$ in Fig. 6. Comparing the soliton profiles with different $\tilde{\rho}_{\rm bg}^{i}$ plotted in every subplot, we find the fine-structure of the background baryon profile does not matter and our treatment for the ellipsoidal baryon profile is reasonable. And again, comparing the soliton profiles with different $m$, we find the heavier FDM particles form a more compact denser soltion which is more insensitive to the background baryon profile.

In this paper, we enumerate some reasonable variants of the exact SP system, such the variants due to a central supermassive BH (section II), the system’s own angular momentum (section III), an extra denser and compact soliton (section IV) and an ellipsoidal baryon background (section V). All of them can be considered as an almost spherical system and solved with the shooting method. We reach the same conclusion Davies:2019wgi that the larger central supermassive BH makes the soliton more denser and compact. For the first time, we prove the gravitomagnetic field derived from the angular momentum can be neglected compared to the contributions from the gravitoelectric field which is the only survival in the Newtonian regime. Unlike the interference or collison of two similar solitons, the beginning of the interaction between two solitions with an extreme density-ration is relatively simple: the smaller soliton with higher density is almost not affected but the larger soliton with lower density shrink dramatically. Finally, we revisit the topic discussed in Bar:2019bqz that the baryon background is important for the solition profile. But we propose a reasonable treatment for the ellipsoidal baryon background and then highlight the importance of the baryon background for lighter FDM particle.

Of course, we can enlarge the parameter space of the underlying theoretical model and produce a greater diversity in the solution of soliton. For example, the self-interactions between FDM particles Chavanis:2019bnu or the modified gravitational potential due to the modified gravity Chen:2024pyr would result in a different coupled equations, hence a different soliton profile. Although the Universe has picked out the only underlying rule, these variants of the exact SP system presented in this paper are common to the theoretical models.