Diverse dark matter haloes in Two-field Fuzzy Dark Matter

Cosmology in the presence of multiple axions, namely “the axiverse”, is generically expected in string theory [1, 2]. Predictions for the axion mass spectrum depend on the details of the precise dynamics of the dimensional compactification that describes our universe, which so far remains unknown. Knowing that axion masses are generated only by non-perturbative (instanton) effects, typical axion masses are exponentially smaller than the standard model scale, and the axion spectrum covers a wide range in masses, which may even reach as low as the Hubble scale of $10^{-33}~{\rm eV}$. The fuzzy dark matter (FDM) proposal composed of ultralight axions of mass $\sim 10^{-22}$ eV [3, 4, 5, 6] is well motivated from this perspective.

Cosmological simulations of FDM reveal rich wave-like structures of constructive and destructive interference on the de Broglie scale of $\sim 1~{\rm kpc}$, including the NFW-like haloes of galaxies and a stable “soliton” ground state at the core of every collapsed structure [7, 8, 9, 10]. This soliton formation and wave-like behaviour are unique properties that distinguish FDM from cold dark matter (CDM) [11, 12] and provide an alternative solution to the small-scale issues of the standard $\Lambda$CDM model [13, 14, 15, 16], without invoking the sub-grid physics of “baryonic feedback” [17, 18, 19] that is yet to be understood.

However, the canonical model of FDM faces several challenges from observational Lyman-$\alpha$ data, which excludes either the particle mass up to $10^{-20}~{\rm eV}$ [20, 21, 22] or the cosmological abundance to less than $30\%$ [23]. Stellar dynamics inside galaxies and satellites also places stringent constraints on the mass spectrum [24, 25, 26, 27, 28, 29]. What has become increasingly clear is that FDM haloes are not able to simultaneously fit two distinguishable classes of dwarf galaxies, namely the classical dwarf spheroidals (dSphs) [30] and the physically much smaller and less luminous of ultra-faint dwarfs (UFDs) [31]. The primary concern is that star clusters in UFDs tend to be overheated by density fluctuations from the soliton or interference granules [25, 29], unless they are significantly stripped away via tidal evolution [32, 33]. This paradoxical situation of FDM is reminiscent of the diversity problem in galactic rotation curves encountered in the standard CDM paradigm [34].

Meanwhile, it has been recently suggested that a cosmological model with multiple ultralight axions, or equivalently multiple FDM species [35, 36, 37], may accommodate diverse profiles of dark matter haloes [38, 39, 40, 41, 42], hence circumvents the aforementioned issue with dwarf galaxies. The two-field fuzzy dark matter (2FDM) model with particle masses separated by 1-2 orders of magnitude provides the most simplified construction of the axiverse where “halo diversity” is anticipated [35, 43]. Specifically, density fluctuations in separate patches of the universe with varying fractions of the 2FDM species would form haloes of different inner core profiles at late times. Despite the feasibility of this idea, cosmological simulations for the 2FDM model have not been able to capture the relatively wide range of de Broglie wavelengths spanned by both FDM species due to the enhanced dynamical range required by a large particle mass hierarchy.

In this paper, we make the first attempt to simulate the structure formation for two axion species differing in particle mass by a factor of 5, the highest ratio by far for cosmological simulations of this type. We develop an optimized spectral solver for 2FDM equations of motion with self-consistent initial conditions solved from 2FDM linear perturbations. The improved simulation technique combined with the large mass factor allows$-$for the first time$-$investigation of dark matter haloes and their diverse structures with unprecedented resolution.

In the non-relativistic limit, the dynamics of the 2FDM model is governed by the coupled Schrödinger-Poisson (SP) equations, which can be written explicitly in the comoving coordinates as

Here $\psi_{1}$ and $\psi_{2}$ represent wavefunctions of the 2FDM fields with a mass $m_{1}$ and $m_{2}$, respectively. $\Phi$ denotes a gravitational potential sourced by the total density, $\rho=\rho_{1}+\rho_{2}\equiv|\psi_{1}|^{2}+|\psi_{2}|^{2}$. $\bar{\rho}$ is the average density over the comoving volume of interest and $a$ is the scale factor. For clarity reasons, let us refer to the light field and the heavy field as $\psi_{1}$ and $\psi_{2}$ from now on, assuming $m_{2}>m_{1}$. We will also refer to the total field as the sum of $\psi_{1}$ and $\psi_{2}$ in terms of dark matter density.

The SP equations in (1) can be solved numerically by the pseudo-spectral method which evolves the system unitarily via a series of “kick” and “drift” steps as described in Ref. [41]. To achieve desired performance and scalability, we develop a fully parallelized solver optimized for the 2FDM model in this work.

We perform a high-resolution cosmological simulation in the 2FDM model where $m_{1}=10^{-22}~{\rm eV}$ and $m_{2}=5\times 10^{-22}~{\rm eV}$, with a density ratio of $\beta_{2}\equiv\Omega_{2}/\Omega_{m}=0.7$. Although the chosen values of $m_{1}$ and $m_{2}$ are in tension with observations, the relevant physics still applies for higher particle masses of the same ratio $m_{2}/m_{1}$ thanks to the scale invariance of SP equations. The simulation volume has a side length of $1.7~{\rm Mpc}/h$ in each dimension. The spectral resolution is $2048^{3}$ for the 2FDM fields. Periodic boundary condition is automatically applied in the spectral solver. The background cosmology is set up with cosmological parameters from the most updated Planck data [44] except for $A_{s}=10^{-8}$. This enhanced value of $A_{s}$ is to compensate for the density fluctuations from large-scale modes in a small simulation volume.

To obtain self-consistent initial conditions, at first we employ N-GenIC [45] to generate a realization of initial random phases in the momentum space. The initial power spectra are then computed with the Boltzmann solver CAMB [46] specifically modified for the 2FDM model. Finally, the initial wavefunctions of the 2FDM fields are solved from Madelung (fluid) equations. More details about how the initial conditions are set up can be found in the Supplemental Material.

We evolve the 2FDM system from the starting redshift of $z=127$ to the final redshift of $z=3.4$. In terms of convergence, we find that some smallest features are slightly under-resolved at the final redshift, but they do not affect the main results discussed below.

Fig. 1 shows how structures evolve in the 2FDM cosmology. We observe in total 13 haloes formed by the end of our simulation. These haloes are numbered based on their (approximate) formation history where smaller numbers indicate haloes that form earlier.

To understand the halo evolution, we compare the radial profiles of a few haloes from their formation redshift ($z=z_{f}$) to the final redshift ($z=3.4$) in Fig. 2. Only four haloes are displayed here as they represent different viable evolution patterns found among the other ones. Generally, every halo starts off with a density fraction of $\psi_{2}$ higher than that of $\psi_{1}$ as expected from their initial cosmological abundance. Each halo, however, evolves to a distinct final state depending on its formation time.

In a few haloes that form first in the simulation, such as Halo 1, $\psi_{2}$ reaches a stable configuration shortly while the density of $\psi_{1}$ keeps growing until $z=3.4$. As a result, $\psi_{1}$ becomes comparable to $\psi_{2}$ in terms of mass and density content at the last redshift. Interestingly, Halo 1 is just a typical example to show that the higher cosmological density of $\psi_{2}$, i.e., $\beta_{2}>0.5$, does not always translate to its dominance inside individual haloes.

In other haloes that form at later redshifts, such as Halo 8, the 2FDM fields experience a similar growth but $\psi_{1}$ ends up a sub-dominant component compared to $\psi_{2}$. It seems that the distribution of $\psi_{2}$ in Halo 8 is not massive enough to support the clustering of $\psi_{1}$. In addition, the mass accumulation rate of $\psi_{1}$ here is relatively slow, which means $\psi_{1}$ might have already settled in its virialized state without further growth.

There are also extreme cases such as Halo 12 which is severely deficient of $\psi_{1}$. The lack of $\psi_{1}$ in this halo can be explained by its late formation time ($z_{f}=3.6$). Since $\psi_{1}$ only accounts for $30\%$ of the total DM budget, most of $\psi_{1}$ abundance has already clustered in haloes that form earlier. Thus, we anticipate that any haloes forming later than Halo 12 are also completely dominated by $\psi_{2}$.

Halo 3 is somewhat special as it undergoes a separate track from the remaining haloes. Since this halo locates in the neighborhood of another massive object, i.e., Halo 1, it evolves under a gravitational potential during the entire lifetime. Thus, we observe strong tidal disruption in the density profile of $\psi_{1}$ and $\psi_{2}$ when Halo 3 approaches Halo 1. Most notably, only the mass of $\psi_{2}$ is significantly stripped away while most of $\psi_{1}$ mass remains intact. The reason why there exists such an asymmetry is unknown at the moment, but this phenomenon indicates a possible mechanism to create $\psi_{1}$-dominated haloes.

Another curious scenario is the coalescence of Halo 6 and 7 into a more massive halo, as seen in Fig. 1. Halo mergers such as this pair are expected to be common in a larger volume. Here the evolution of $\psi_{1}$ and $\psi_{2}$ similarly follow those of Halo 1, i.e., they also yield an equivalent amount of $\psi_{1}$ and $\psi_{2}$ at the end.

As the virialized configuration of haloes are not universal in a 2FDM cosmology, it is important to examine the demographics of each halo population in more details.

Fig. 3 provides a complete landscape of the comoving volume and 12 haloes found at $z=3.4$. The projected densities (top-left panels) show an almost identical filamentary structure of $\psi_{1}$ and $\psi_{2}$ on large scales and the lower density distribution of $\psi_{1}$ compared to the one of $\psi_{2}$. The sliced densities (bottom-left panels) give a close-up view of individual haloes with their central solitons surrounded by density granules from wave interference. Most important of all, the radial profiles (right panels) clearly illustrate a diversity of haloes with distinct solitonic core structures. We find that 2FDM haloes can be separated into three populations as follows.

The first population consists of Halo 1, 2, 4, 5 and 6(7). The radial profiles of these haloes show a high fraction of $\psi_{1}$ compared to $\psi_{2}$ with the presence of two central cores. In Halo 1, 2, 5 where the cores of $\psi_{1}$ and $\psi_{2}$ are approximately concentric ($\Delta r_{12}<1.5~{\rm kpc}$), the simulation profiles perfectly match the so-called nested soliton, namely the ground-state solution of the time-independent SP equations [35, 41]. On the contrary, in Halo 4 and 6(7) where the two cores are not aligned ($\Delta r_{12}>2~{\rm kpc}$) due to soliton random walk [32], the nested soliton does not yield a good fit as expected. Overall, the haloes belonging to this population would be observed as ones with centrally nested structures (from the total field perspective).

The second population consists of Halo 8, 10, 11, 12 and 13. These haloes include an extremely low to moderate amount of $\psi_{1}$. In this case, even when the two cores of $\psi_{1}$ and $\psi_{2}$ form in a concentrically nested configuration, only the soliton of $\psi_{2}$ is visible as it overwhelmingly dominates the one of $\psi_{1}$. As the outer regions of these haloes are also dominated by $\psi_{2}$, they would be most likely observed as $\psi_{2}$-only haloes.

The third population consists of Halo 3 and 9. As previously mentioned, Halo 3 has been subject to tidal interactions since its formation. This effect causes a complete disruption of the $\psi_{2}$ soliton, so that Halo 3 is eventually dominated by $\psi_{1}$ with its soliton. A similar phenemenon also occurs with Halo 9 during its evolution under the gravitational potential of Halo 4 (see Fig. 1), resulting in a considerable mass loss of $\psi_{2}$. Although the density of $\psi_{2}$ is still comparable to the one of $\psi_{1}$ here, there is no $\psi_{2}$ soliton formed in this halo (see the sliced and projected densities). Eventually, Halo 3 and 9 would be observed as $\psi_{1}$-only haloes.

Can we seek the aforementioned halo populations in observational data? Previous studies [35, 43] suggested that the DM profiles of dSphs and UFDs can be explained by two axion species with $m_{1}\sim 10^{-22}~{\rm eV}$ and $m_{2}\sim 10^{-20}~{\rm eV}$, respectively.

If we assume that the current simulation results can be extrapolated to a larger value of $m_{2}$, the $\psi_{1}$-only and $\psi_{2}$-only populations mentioned above would satisfy the observational constraints of dSphs and UFDs, respectively. However, there are two problems with this identification. Firstly, the $\psi_{2}$-only haloes are more common than the $\psi_{1}$-only ones, but the number of presently detected UFDs are of the same order of magnitude with dSphs. Secondly, $\psi_{2}$-only haloes generally form later than $\psi_{1}$-only haloes whereas UFDs are believed to be much older than dSphs [31]. It is possible that UFDs should be identified with extremely-low-$\psi_{1}$ haloes like Halo 12, which would appear much earlier in a higher-$m_{2}$ 2FDM simulation. However, we need cosmological simulations with a proper implementation of star formation and baryonic physics to address these problems. At the moment there exist feasible but yet definitive connections between our simulation haloes and the observed dwarf galaxies.

On the other hand, the first population with nested soltions can be identified with normal galaxies such as Milky Way (MW). There are, in fact, some observational evidence of such nested profile in the MW center. For instance, Ref. [47] found that the excess velocity dispersion of the central MW bulge stars could be accounted for by a $\psi_{1}$ soliton of a mass $10^{9}~{\rm M}_{\odot}$. Meanwhile, a similar analysis of FDM in the nuclear star cluster of MW [48] found a positive hint for a soliton of an FDM field with a mass $10^{-20.5}~{\rm eV}$, which is approximately the boson mass expected for $\psi_{2}$. As such, we may hope to continue searching for nested solitons in other galaxies in the near future.

The existence of $\psi_{1}$-deficient galaxies such as Halo 12 for a 2FDM cosmology with non-negligible $\Omega_{1}$ also affects the existing limits on canonical FDM derived from data of isolated systems [29, 49]. Since the considered system may be hosted by a 2FDM halo like this one, their constraints should apply to the particle mass of the heavier species of FDM, ${\it i.e.}$, $m_{2}$, instead of $m_{1}$ as the most conservative estimation. In addition, it is necessary to revisit these constraints with tidal stripping taken into account because the stellar heating effect caused by transient density fluctuations generic to FDM are highly suppressed in that context [32, 33].

We have demonstrated that the 2FDM cosmology can provide rich and complex structure formation, with the diversity of dark matter haloes being one of its most distinguishing features. Most young and small haloes are dominated by $\psi_{2}$ and they only see $\psi_{2}$ solitons as the central major components. Old and more massive haloes, on the other hand, incorporate comparable amounts of $\psi_{1}$ and $\psi_{2}$. Hence, these haloes host solitons with observable nested structures. Lastly, haloes dominated by $\psi_{1}$ may emerge via tidal disruption in natural encounters with nearby haloes.

It is important to note that these results can be subject to some limitations of our simulation. For instance, the major restriction of the spectral method (uniform resolution) is that we can only evolve the system to a certain point before features become smaller than what can be represented by the maximum wavenumber in the simulation. As a consequence, the peaked (central) density of $\psi_{2}$ in each halo can be underestimated when its soliton is not completely resolved. This drawback, however, does not affect any qualitative conclusions about halo diversity. Another caveat is that the above results only apply to the 2FDM system where $m_{2}/m_{1}=5$ and $\beta_{2}=0.7$. Even though we also examine simulations with other input parameters of $m_{2}/m_{1}$ and $\beta_{2}$, dark matter haloes here are typically dominated by a single population of haloes, i.e., no halo diversity is realized. Extended discussions on these scenarios can be found in the Supplemental Material.

With the first compelling evidence for halo diversity clearly shown in this study, future work could aim at simulating larger cosmological volumes with higher resolution at lower redshifts, if possible. The string axiverse, as represented by the 2FDM model, still holds surprisingly many yet-to-be-discovered potentials, which will hopefully open a new window to the long-standing dark matter mystery of our time.

HN appreciates Dr. Jens Stücker, Prof. Raúl Angulo and Prof. Simon White for valuable discussions and insights to the final results of this paper. HN also thanks Prof. Raúl Angulo and Donostia International Physics Center (DIPC) for providing a simulation storage on DIPC Supercomputing Center. HN, TB, HT, TL and GS acknowledge support from the Research Grants Council of Hong Kong through the Collaborative Research Fund C6017-20G. LH acknowledges support by the Simons Collaboration on “Learning the Universe”. The simulations in this work were run on the Engaging cluster sponsored by Massachusetts Institute of Technology (MIT). TB is supported by the Spanish grant PID2023-149016NBI00 (funded by MCIN/AEI/10.13039/501100011033 and by ”ERDF A way of making Europe”.