## Diverse dark matter haloes in Two-field Fuzzy Dark Matter

Hoang Nhan Luu , 1, 2, ∗ Philip Mocz , 3, † Mark Vogelsberger , 2, ‡ Alvaro Pozo , 1, 4 Tom Broadhurst , 1, 4, 5 S.-H. Henry Tye , 6, 7 Tao Liu , 6 Leo W.H. Fung , 6 George F. Smoot , 1, 6, 8, 9 Razieh Emami , 10 and Lars Hernquist 10

1 Donostia International Physics Center, Basque Country UPV/EHU, San Sebastian, E-48080, Spain

2 Department of Physics, Kavli Institute for Astrophysics and Space Research,

Massachusetts Institute of Technology, Cambridge, MA 02139, USA

3 Flatiron Institute, 162 5th Ave, New York, NY, 10010, USA

4 University of the Basque Country UPV/EHU, Department of Theoretical Physics, Bilbao, E-48080, Spain

5 Ikerbasque, Basque Foundation for Science, Bilbao, E-48011, Spain

6 Department of Physics and Institute for Advanced Study,

The Hong Kong University of Science and Technology, Hong Kong

7 Department of Physics, Cornell University, Ithaca, NY 14853, USA

8 Paris Centre for Cosmological Physics, APC, AstroParticule et Cosmologie,

Universit´ e de Paris, CNRS/IN2P3, CEA/lrfu, 10, rue Alice Domon et Leonie Duquet, 75205 Paris CEDEX 13, France emeritus

9 Physics Department, University of California at Berkeley, CA 94720, Emeritus

10 Center for Astrophysics, Harvard &amp; Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

(Dated: August 5, 2024)

Fuzzy dark matter (FDM) is a compelling candidate for dark matter, offering a natural explanation for the structure of diffuse low-mass haloes. However, the canonical FDM model with a mass of 10 -22 eV encounters challenges in reproducing the observed diversity of dwarf galaxies, except for possibly scenarios where strong galactic feedback is invoked. The introduction of multiple-field FDM can provide a potential resolution to this diversity issue. The theoretical plausibility of this dark matter model is also enhanced by the fact that multiple axion species with logarithmicallydistributed mass spectrum exist as a generic prediction of string theory. In this paper we consider the axiverse hypothesis and investigate non-linear structure formation in the two-field fuzzy dark matter (2FDM) model. Our cosmological simulation with an unprecedented resolution and selfconsistent initial conditions reveals the diverse structures of dark matter haloes in the 2FDM model for the first time. Depending on the formation time and local tidal activities, late-time haloes can host solitons of nested cores or solitons of one dominant species.

Introduction. 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 eV. The fuzzy dark matter (FDM) proposal composed of ultralight axions of mass ∼ 10 -22 eV [3-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 ∼ 1 kpc, including the NFW-like haloes of galaxies and a stable 'soliton' ground state at the core of every collapsed structure [7-10]. This soliton formation and wave-like behaviour are unique properties that distinguish the FDM from the cold dark matter (CDM) [11, 12] and provide an alternative solution to the small-scale issues of the standard ΛCDM model [1316], without invoking the sub-grid physics of 'baryonic feedback' [17-19] that is yet to be understood.

challenges from observational Lymanα data, which excludes either the particle mass up to 10 -20 eV [20-22] or the cosmological abundance to less than 30% [23]. Stellar dynamics inside galaxies and satellites also placed stringent constraints on the mass spectrum [24-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 disruptions [32, 33]. Nevertheless, this paradoxical situation of FDM is reminiscent of the diversity problem in galactic rotation curves encountered in the standard CDM paradigm [34].

However, the canonical model of FDM faces several

Meanwhile, it has been recently suggested that a cosmological model with multiple ultralight axions, or equivalently multiple FDM species [35-37], may accommodate diverse profiles of dark matter haloes [38-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. Although this idea is theoretically viable, 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 improved spectral solver for 2FDM equations of motion with self-consistent initial conditions solved from 2FDM linear perturbations. The improved simulation combined with the large mass factor allows -for the first time -investigation of dark matter haloes and their diverse structures with unprecedented resolution.

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

$$& \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \\ & i \hbar { \frac { \partial \psi _ { 1 } } { \partial t } } = - \frac { \hbar { 2 } } { 2 m _ { 1 } a ^ { 2 } } \nabla ^ { 2 } \psi _ { 1 } + \frac { m _ { 1 } } { a } \Phi \psi _ { 1 }, \quad & \text{fig. } \text{tion} \\ & i \hbar { \frac { \partial \psi _ { 2 } } { \partial t } } = - \frac { \hbar { 2 } } { 2 m _ { 2 } a ^ { 2 } } \nabla ^ { 2 } \psi _ { 2 } + \frac { m _ { 2 } } { a } \Phi \psi _ { 2 }, \quad & \text{(1)} \\ & \nabla ^ { 2 } \Phi = 4 \pi G \left ( | \psi _ { 1 } | ^ { 2 } + | \psi _ { 2 } | ^ { 2 } - \bar { \rho } \right ). \quad & \text{that}$$

Here, ψ 1 and ψ 2 represent wavefunctions of the 2FDM fields with a mass m 1 and m 2 , respectively. Φ denotes the gravitational potential sourced by the total density, ρ = ρ 1 + ρ 2 ≡ | ψ 1 | 2 + | ψ 2 | 2 . ¯ 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 ψ 1 and ψ 2 from now on, assuming m &gt; m 2 1 . We will also refer to the total field as the sum of ψ 1 and ψ 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 [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 eV and m 2 = 5 × 10 -22 eV, with a density ratio of β 2 ≡ Ω 2 / Ω m = 0 7. . The chosen values of m 1 and m 2 are in tension with observations, but the relevant physics still applies for higher particle masses of the same ratio m /m 2 1 thanks to the scale invariance of SP equations. The simulation volume has a side length of 1 7 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.

We employ N-GenIC [44] to generate initial random

FIG. 1. Projected densities of the total field in the simulation volume at several redshifts. The white circles mark all haloes formed with a solitonic structure at the corresponding redshift denoted in each panel. The number adjacent to each halo indicates (approximately) the order of formation where smaller numbers correspond to earlier-forming haloes. Note that Halo 6 and 7 have merged in the last redshift.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00827v1-figure-ee5e33e3-6bec-4b08-b1d9-aa333cdefb93.jpg)

phases in the momentum space. The initial power spectra are computed with a modified version of the Boltzmann solver CAMB [45]. The initial wavefunctions of the 2FDM fields are then solved from Madelung (fluid) equations. The background cosmology is set up with cosmological parameters from the most updated Planck data [46] 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. More details about how the initial conditions are set up can be found in the supplemental materials.

We then 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.

Structure formation in the 2FDM cosmology. 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 earlier-forming haloes.

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

FIG. 2. Radial profiles of some representative haloes at several redshifts. In each panel, z f denotes the first redshift at which the corresponding halo forms with a solitonic structure in the simulation. Blue, green and red curves correspond to profiles of ψ , ψ 1 2 and the total field respectively. H x -ψ 1 and H x -ψ 2 denote ψ 1 and ψ 2 profiles of Halo x . Since the density of ψ 1 in Halo 12 is too low, its profiles (in faded blue) are multiplied by 10 for better illustration.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00827v1-figure-c53d77c1-14d1-4fd5-8722-223a68866b95.jpg)

four haloes are displayed here as they represent different viable evolution patterns found among the other ones. Generally, every halo starts off with the density fraction of ψ 2 higher than that of ψ 1 as expected from their initial cosmological abundance. Each halo, however, yields a distinct final state depending on its formation time.

In a few haloes that form first in the simulation, such as Halo 1, ψ 2 reaches a stable configuration shortly while the density of ψ 1 keeps growing until the last redshift. As a result, ψ 1 becomes comparable to ψ 2 in terms of mass and density content at z = 3 4. . Interestingly, Halo 1 is just a typical example to show that the higher cosmological density of ψ 2 , i.e. β 2 , 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 ψ 1 ends up a sub-dominant component compared to ψ 2 . It seems that the distribution of ψ 2 in Halo 8 is not massive enough to support the clustering of ψ 1 . In addition, the mass accumulation rate of ψ 1 here is relatively slow, which means ψ 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 ψ 1 . The lack of ψ 1 in this halo can be explained by its late formation time ( z f = 3 6). . Since ψ 1 only accounts for 30% of the total DM budget, most of ψ 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 ψ 2 .

Halo 3 is somewhat special as it undergoes a separate evolution from the remaining haloes. Since this halo locates in the neighborhood of another massive object, i.e. , Halo 1, it is constantly exerted by a strong gravitational potential. Thus, we observe tidal disruptions in the density profile of ψ 1 and ψ 2 when Halo 3 approaches Halo 1. Most notably, only the mass of ψ 2 is significantly stripped away while most of ψ 1 mass remains intact. The reason why there exists such asymmetry in mass stripping is unknown at the moment, but this phenomenon indicates a possible mechanism to create ψ 1 -dominated haloes.

Another curious case 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 ψ 1 and ψ 2 similarly follow those of Halo 1, i.e., they also yield an equivalent amount of ψ 1 and ψ 2 at the end.

Diverse structures of 2FDM haloes. We have shown that the virialized configuration of haloes are not universal the 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 ψ 1 and ψ 2 on large scales and the lower density distribution of ψ 1 compared to the one of ψ 2 . The sliced densities (bottom-left panels) give a closeup view of individual haloes with their central solitons surrounded by density granules from wave interference. Most importantly, the radial profiles (right panels) reveal the 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 ψ 1 compared to ψ 2 with the presence of two central cores. In Halo 1, 2, 5 where the cores of ψ 1 and ψ 2 are approximately concentric (∆ r 12 &lt; 1 5 kpc), the . simulation profiles perfectly match the so-called nested soliton , namely the ground-state solution of the timeindependent SP equations [35, 41]. On contrary, in Halo 4,6(7) where the two cores are not aligned (∆ r 12 &gt; 2 kpc) due to soliton random walk [32], the nested soliton does not yield a good fit as expected. The 2FDM haloes in this group would be observed as haloes with centrally nested structures (from their total-field profiles).

The second population consists of Halo 8, 10, 11, 12

FIG. 3. (Left top) Projected densities of ψ 1 and ψ 2 at z = 3 9. . The white circles mark the haloes at this redshift as in Fig. 1. (Left bottom) Sliced densities of the haloes marked in the above panels. Each slice displays a cross section through the center of the associated halo with a side length of 200 kpc /h . (Right) Radial profiles of all haloes at z = 3 4. . The solid curves show the simulation profiles of the total field (red), ψ 1 (blue), and ψ 2 (green). The dashed curves show the best-fit soliton profiles of ψ 1 (blue) and ψ 2 (green) in each halo. The halo center is chosen as the barycenter of the total field. As the barycenters of ψ 1 and ψ 2 are not always aligned, ∆ r 12 denotes the distance between them in the units of kpc. The soliton profile of only one field is shown in Halo 3, 9 and 12 because the other field does not host a soliton. ∆ 12 is also undefined for these haloes.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00827v1-figure-ff228591-438d-47bc-9783-e427dfe61141.jpg)

and 13. These haloes include an extremely low to moderate amount of ψ 1 . In this case, even when the two cores of ψ 1 and ψ 2 form in a concentrically nested configuration, only the soliton of ψ 2 is visible as it overwhelmingly dominates the one of ψ 1 . As the outer regions of these haloes are also dominated by ψ 2 , they would be most likely observed as ψ 2 -only haloes.

dSphs and UFDs can be explained by two axion species with m 1 ∼ 10 -22 eV and m 2 ∼ 10 -20 eV, respectively.

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 ψ 2 soliton, so that Halo 3 is eventually dominated by ψ 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 ψ 2 . Although the density of ψ 2 is still comparable to the one of ψ 1 here, there is no ψ 2 soliton formed in this halo (see the sliced and projected densities). Eventually, Halo 3 and 9 would be observed as ψ 1 -only haloes.

Halo diversity in observation. Can we seek the aforementioned halo populations in observational data? Previous studies [35, 43] suggested that the DM profiles of

If we assume that the current simulation results can be extrapolated to m 2 within this mass range, the ψ 1 -only and ψ 2 -only populations would satisfy the observational constraints of dSphs and UFDs, respectively. However, there are two problems with this explanation. Firstly, the ψ 2 -only haloes are more common than the ψ 1 -only ones, but the number of presently detected UFDs are of the same order of magnitude with dSphs. Secondly, ψ 2 -only haloes can form at later times than ψ 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ψ 1 haloes like Halo 12, which would appear much earlier in a higherm 2 2FDM simulation. However, at the moment we have feasible but yet definitive connections between our simulation haloes and 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 exists, 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 ψ 1 soliton of 10 9 M . Additionally, in an analy-⊙ sis of FDM in the nuclear star cluster of MW, the authors of [48] found a positive hint for a soliton of an FDM field with a mass 10 -20 5 . eV, which is approximately the boson mass that are expected for ψ 2 . As such, we may hope to continue searching for nested solitons in other galaxies in the near future.

Conclusion and outlook. 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 ψ 2 and they only see ψ 2 solitons as the central major components. Old and more massive haloes, on the other hand, incorporate comparable amounts of ψ 1 and ψ 2 . Hence, these haloes host solitons with observable nested structures. Lastly, haloes dominated by ψ 1 may emerge via tidal disruption in rare 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 ψ 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 /m 2 1 = 5 and β 2 = 0 7. Even though we also examine simulations with . other input parameters of m /m 2 1 and β 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 materials.

With the first compelling evidence for halo diversity discovered in this study, future work could aim at simulating larger cosmological volumes with higher resolution at lower redshifts, if possible. It turns out that 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.

Acknowledgments. HN acknowledges several discussion from Jens St¨cker, Prof. Ra´l Angulo and Prof. Siu u mon White who have provided valuable insights to the final results of this paper. HN also thanks Prof. Ra´l Anu gulo 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).

## Supplemental Materials

The supplemental materials include additional information that may be useful for readers interested in technical details. In Sec. A we provide technical details about how to derive the initial conditions for the 2FDM cosmology. In Sec. B we compare and discuss 2FDM simulations having different set of input parameters.

## A. Initial conditions

In single-field FDM simulations, the matter power spectrum is typically derived using numerical Boltzmann codes such as axionCAMB [49, 50] or using the approximate axion transfer function (at z = 0) given by [3]

$$T ( k ) = \sqrt { \frac { P _ { \text{FDM} } } { P _ { \text{CDM} } } } = \frac { \cos x _ { J } ^ { 3 } } { 1 + x _ { J } ^ { 8 } }, \text{ \quad \ \ } ( 2 )$$

where x J = 1 61 . m 1 / 18 22 k/k J, eq ; k J, eq = 9 m 1 / 2 22 Mpc -1 and m 22 = m/ 10 -22 eV. Note that the transfer function (2) has a scale-dependent growth and assumes the total dark matter budget composed of axions.

So far, FDM initial conditions are usually extrapolated from the N-body particle distribution via density assignment such as the cloud-in-cell algorithm because it is particularly convenient to generate N-body initial conditions via publicly available codes [44, 51] in the standard ΛCDM cosmology, even for an arbitrary power spectrum. However, this approach is rather unreliable for ICs generation of field-based simulations for two reasons. Firstly, the velocities of N-body particles are typically computed with the Zel'dovich approximation or the second-order Lagrangian Perturbation Theory that is tailored for CDM perturbations. Secondly, particle discreteness might introduce abnormal mesh noise on the smallest scales of FDM simulations when being converted to fields. These effects might cause an artificial formation of compact objects at late times, especially for FDM with a high particle mass. As such, the most robust and selfconsistent approach is to derive the initial wavefunction of FDM directly from its linear perturbations in the early universe, which is implemented as follows.

For a set of random phases { θ k } in the momentum space, the initial density fluctuations of the FDM field can be solved via an inverse Fourier transform

$$\delta ( { \mathbf x } ) \leftarrow \text{ifft} \left [ \Delta _ { k } e ^ { i \theta _ { \mathbf k } } \right ], \ \Delta _ { k } ^ { 2 } = k ^ { 3 } P _ { k } / ( 2 \pi ) ^ { 3 }. \quad ( 3 )$$

Here P k is computed at the initial redshift, e.g., z ini = 127 in our simulation, which is then related to the present spectrum of CDM by the linear growth factor and the

FIG. A. (Left) Linear power spectra at the initial redshift z ini = 127. The spectra of ψ 1 and ψ 2 in the 2FDM model are shown as the blue-solid and green-solid curves. The spectra of the light ( m 22 = 1) and heavy ( m 22 = 5) FDM species are shown as the blue-dashed and green-dashed curves. The spectrum of CDM is also shown as the black-solid curve for comparison. The dotted vertical line indicates the minimum wavenumber (the largest scale) resolved by the simulation box. (Right) Projected densities of initial fluctuations with matching phases in various dark matter models. The 2FDM/FDM initial conditions are generated as wavefunctions on a uniform grid, while the CDM ones are generated as N-body particles with N-GenIC [44].

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00827v1-figure-eca685fb-ec8f-4b74-b216-feb52d2aa96a.jpg)

transfer function, P k = D 2 ( z ini ) /D 2 (0) T 2 ( k P ) CDM . This procedure can be repeated to find the (conformal) time derivatives of the density field δ ′ ( x ) = ∂δ/∂η , which are analogous to the particle velocities in N-body simulations. Both δ ( x ) and δ ′ ( x ) are necessary degrees of freedom to derive the complex-valued wavefunction of the FDM field in the Schr¨dinger equation o

$$i \hbar { \frac { \partial \psi } { \partial t } } = - \frac { \hbar { ^ { 2 } } } { 2 m ^ { 2 } } \nabla ^ { 2 } \psi + \frac { m } { a } \Phi \psi \,. \quad \quad ( 4 ) \quad \text{ice} \quad \text{in c}$$

If we decompose ψ = ψ e r iα and define the Madelung velocity as v M ≡ ℏ /m ∇ α , the continuity equation yields a relation of α and δ ′ at the first perturbative order as

$$\frac { \partial \rho } { \partial t } + \nabla \cdot ( \rho v _ { M } ) = 0 \ \rightarrow \ \nabla ^ { 2 } \alpha = - \frac { m } { \hbar } a \delta ^ { \prime } \,. \ \ ( 5 ) \quad \text{gravi} \\ \quad \text{that} \\ \quad \text{--web}.$$

We then solve Eq. (5) to obtain the phase α while the modulus ψ r can be easily inferred from the density contrast by ρ = Ω DM cr ρ (1 + δ ) = | ψ r | 2 , which makes the wavefunction ψ fully specified.

These computations can be generalized for two axion species in the 2FDM model, providing δ 1 , δ 2 and δ ′ 1 , δ ′ 2 . However, the transfer function of the two axion species is no longer given by (2). Instead, we need to solve the linear perturbation equations of both fields

$$\dot { \delta } _ { i } & = - k u _ { i } - \dot { h } / 2 - 3 \mathcal { H } c _ { s, i } ^ { 2 } \delta _ { i } - 9 \mathcal { H } ^ { 2 } c _ { s, i } ^ { 2 } u _ { i } / k \,, \quad \text{(6)} \quad \text{Fij} \\ \dot { u } _ { i } & = - \mathcal { H } u _ { i } + k c _ { s, i } ^ { 2 } \delta _ { i } + 3 \mathcal { H } c _ { s, i } ^ { 2 } u _ { i } \,, \quad \text{(6)} \quad \text{(k)} \quad \text{that}$$

where u i is the heat flux and h is the trace of the metric perturbations in the synchronous gauge. The effective sound speed and the conformal Hubble function are defined as

$$c _ { s, i } ^ { 2 } \equiv \frac { k ^ { 2 } / ( 4 m _ { i } ^ { 2 } a ^ { 2 } ) } { 1 + k ^ { 2 } / ( 4 m _ { i } ^ { 2 } a ^ { 2 } ) }, \quad \mathcal { H } \equiv \frac { a ^ { \prime } } { a } = a H \,. \quad ( 7 )$$

We note that Eqs. (6) only show the effective description of axion perturbations when the axion oscillations become much faster than the Hubble timescales. In practice, the exact and effective fluid treatment are combined in our calculations for optimal speed and accuracy, following the procedure in the previous works from [49, 52].

At the first glance, the linear equations of each field in (6) seem independent of each other. They are, however, gravitationally coupled to each other via the 'metric' h that are governed by Einstein equations. As a result, perturbations of the sub-dominant field is modulated by the dominant one on all scales. Fig. A (left panels) clearly illustrates this feature of the power spectra of the 2FDM fields. In the 2FDM model, as ψ 2 perturbations accumulate earlier to form potential wells that attract ψ 1 , the spectrum of ψ 1 closely traces that of ψ 2 on all scales. On the other hand, these spectra clearly separate in the single-field FDM model, i.e, the curve of the lighter field (10 -22 eV) has a higher cut-off scale (lower k ) than to the one of the heavier field (5 × 10 -22 eV).

Fig. A (right panels) also shows that the initial 'clumpiness' of the 2FDM field is distinguishable from that of the FDM field. If ψ 1 and ψ 2 were initially set up with the FDM (instead of 2FDM) power spectra, we would have observed a significantly lower concentration

of ψ 1 in every halo, which might weaken the argument for halo diversity. Thus, it is important to emphasize that initial conditions in the 2FDM model must be properly inferred from 2FDM linear perturbations.

## B. Other scenarios

In the main text, our discussions revolve around a specific scenario of the 2FDM model with m /m 2 1 = 5 and β 2 = 0 7. . If this scenario is treated as the 'fiducial' model, a natural question is whether halo diversity can be realized in a more generic model. As such, we have performed other simulations (with a resolution of 1024 3 ) to examine the impacts of input parameters on the 2FDM cosmology. Fig. B (left panels) compares structure formation in these scenarios at z = 5 9 (due to a lower . resolution). By convention we always keep the mass of ψ 1 constant at m 1 = 10 -22 eV and only change m 2 when varying m /m 2 1 .

When the abundance of ψ 2 decreases, as for the case of β 2 = 0 5 or . β 2 = 0 3 with a fixed . m /m 2 1 = 5, the halo number reduces considerably, with three in the former and only one in the latter case, because the initial matter spectrum is more suppressed with more ψ 1 in the dark matter budget. Among the haloes that already form, the central region is completely dominated by a soliton core of ψ 1 , e.g. see the radial profiles of Halo 1 in Fig. B (right panels). Compared to the single-field FDM model, these 2FDM models still have an enhanced small-scale density distribution, but the intrinsic halo structures are almost identical, i.e., only ψ 1 -dominated solitons would be observed.

On the other hand, in the model with a lower mass ratio m /m 2 1 = 3 and a fixed β 2 = 0 7, we only find the . formation of haloes with nested solitons, which can be seen in Halo 1 but also in two other haloes not shown here. This result seems consistent with what was found by [36] using the adaptive mesh refinement method in an equivalent set up. Although the nested soliton is the most unique signature of the 2FDM model, its presence in every halo does not explain the diversity of dwarf galaxies.

There are certainly other aspects that we are unable to study comprehensively due to the limited resolution. Firstly, it is straightforward to notice that the model with m /m 2 1 = 5 and β 2 = 0 5 looks . (almost) statistically identical to the one with m /m 2 1 = 3 and β 2 = 0 7 on . large scales, which may imply a mass-abundance degeneracy. If that is the case, a 2FDM model defined by β 2 and α 2 ≡ m /m 2 1 is analogous to another model with β ′ 2 &lt; β 2 and α ′ 2 &gt; α 2 , and vice versa.

Secondly, it seems halo diversity can only be archived by a combination of a large mass hierarchy and an appropriate density ratio. In case β 2 is too low or too high, one population of FDM would dominate all dark matter haloes at late times. From our simulations, the thresh- old of β 2 seems to fall within 0 6-0 7 for . . m /m 2 1 = 5. Assume that we can make a naive extrapolation from the mass-abundance degeneracy above, a higher-mass ψ 2 would require less β 2 for halo diversity. In other words, ψ 1 would not dominate and suppress the soliton formation of ψ 2 at β 2 &lt; 0 7 thanks to an earlier accumulation . of ψ 2 in every halo, which also agrees with what we found from idealized simulations in the previous study [41].

- ∗ hoang.luu@dipc.org
- † pmocz@flatironinstitute.org
- ‡ mvogelsb@mit.edu
- [1] P. Svrcek and E. Witten, Axions in string theory, Journal of High Energy Physics 2006 , 051 (2006), arXiv:hepth/0605206 [hep-th].
- [2] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, String Axiverse, Phys. Rev. D 81 , 123530 (2010), arXiv:0905.4720 [hep-th].
- [3] W. Hu, R. Barkana, and A. Gruzinov, Cold and fuzzy dark matter, Phys. Rev. Lett. 85 , 1158 (2000), arXiv:astro-ph/0003365.
- [4] D. J. E. Marsh, Axion Cosmology, Phys. Rept. 643 , 1 (2016), arXiv:1510.07633 [astro-ph.CO].
- [5] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Ultralight scalars as cosmological dark matter, Phys. Rev. D 95 , 043541 (2017), arXiv:1610.08297 [astro-ph.CO].
- [6] E. G. M. Ferreira, Ultra-light dark matter, Astron. Astrophys. Rev. 29 , 7 (2021), arXiv:2005.03254 [astro-ph.CO].
- [7] H.-Y. Schive, T. Chiueh, and T. Broadhurst, Cosmic Structure as the Quantum Interference of a Coherent Dark Wave, Nature Phys. 10 , 496 (2014), arXiv:1406.6586 [astro-ph.GA].
- [8] B. Schwabe, J. C. Niemeyer, and J. F. Engels, Simulations of solitonic core mergers in ultralight axion dark matter cosmologies, Phys. Rev. D 94 , 043513 (2016), arXiv:1606.05151 [astro-ph.CO].
- [9] P. Mocz, M. Vogelsberger, V. H. Robles, J. Zavala, M. Boylan-Kolchin, A. Fialkov, and L. Hernquist, Galaxy formation with BECDM - I. Turbulence and relaxation of idealized haloes, Mon. Not. Roy. Astron. Soc. 471 , 4559 (2017), arXiv:1705.05845 [astro-ph.CO].
- [10] S. May and V. Springel, Structure formation in largevolume cosmological simulations of fuzzy dark matter: impact of the non-linear dynamics, Mon. Not. Roy. Astron. Soc. 506 , 2603 (2021), arXiv:2101.01828 [astroph.CO].
- [11] P. Mocz et al. , First star-forming structures in fuzzy cosmic filaments, Phys. Rev. Lett. 123 , 141301 (2019), arXiv:1910.01653 [astro-ph.GA].
- [12] P. Mocz et al. , Galaxy formation with BECDM - II. Cosmic filaments and first galaxies, Mon. Not. Roy. Astron. Soc. 494 , 2027 (2020), arXiv:1911.05746 [astro-ph.CO].
- [13] A. A. Klypin, A. V. Kravtsov, O. Valenzuela, and F. Prada, Where are the missing Galactic satellites?, Astrophys. J. 522 , 82 (1999), arXiv:astro-ph/9901240.
- [14] W. J. G. de Blok, The Core-Cusp Problem, Advances in Astronomy 2010 , 789293 (2010), arXiv:0910.3538 [astroph.CO].
- [15] M. Boylan-Kolchin, J. S. Bullock, and M. Kapling-

## Radial profiles of Halo 1

FIG. B. (Left) Projected densities of the (total) field in several dark matter models at redshift z = 5 9. . The input parameters are indicated in each panel. Initial conditions are calculated with the same phases (same random seed) but with different power spectra. The colored density scale is the same as in Fig. 3. 2FDM and FDM simulations shown here have a spectral resolution of 1024 3 , except for the CDM simulation with an N-body resolution of 512 3 particles and a softening length of 1 66 kpc. . Halo 1 is marked by the white circle in each panel. (Right) The radial profiles of Halo 1 in each dark matter model. The total field, ψ 1 and ψ 2 correspond to the red, blue and green curves. We note that the FDM and CDM model only have one dark matter component, hence only one density profile plotted in red. The vertical line (dashed grey) in the last panel denotes a radius equal to three times of the softening length in the (N-body) CDM simulation.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00827v1-figure-dbd75f50-c960-4168-9583-3f83909cb590.jpg)

hat, Too big to fail? The puzzling darkness of massive Milky Way subhaloes, MNRAS 415 , L40 (2011), arXiv:1103.0007 [astro-ph.CO].

- Phys. Rev. Lett. 119 , 031302 (2017), arXiv:1703.04683 [astro-ph.CO].
- [16] J. S. Bullock and M. Boylan-Kolchin, Small-Scale Challenges to the ΛCDM Paradigm, Ann. Rev. Astron. Astrophys. 55 , 343 (2017), arXiv:1707.04256 [astro-ph.CO].
- [17] F. Governato, A. Zolotov, A. Pontzen, C. Christensen, S. H. Oh, A. M. Brooks, T. Quinn, S. Shen, and J. Wadsley, Cuspy no more: how outflows affect the central dark matter and baryon distribution in Λ cold dark matter galaxies, MNRAS 422 , 1231 (2012), arXiv:1202.0554 [astro-ph.CO].
- [18] A. Di Cintio, C. B. Brook, A. V. Macci` o, G. S. Stinson, A. Knebe, A. A. Dutton, and J. Wadsley, The dependence of dark matter profiles on the stellar-to-halo mass ratio: a prediction for cusps versus cores, Mon. Not. Roy. Astron. Soc. 437 , 415 (2014), arXiv:1306.0898 [astro-ph.CO].
- [19] M. Vogelsberger, F. Marinacci, P. Torrey, and E. Puchwein, Cosmological Simulations of Galaxy Formation, Nature Rev. Phys. 2 , 42 (2020), arXiv:1909.07976 [astroph.GA].
- [20] E. Armengaud, N. Palanque-Delabrouille, C. Y` eche, D. J. E. Marsh, and J. Baur, Constraining the mass of light bosonic dark matter using SDSS Lymanα forest, Mon. Not. Roy. Astron. Soc. 471 , 4606 (2017), arXiv:1703.09126 [astro-ph.CO].
- [21] V. Irˇiˇ, s c M. Viel, M. G. Haehnelt, J. S. Bolton, and G. D. Becker, First constraints on fuzzy dark matter from Lymanα forest data and hydrodynamical simulations,
- [22] K. K. Rogers and H. V. Peiris, Strong Bound on Canonical Ultralight Axion Dark Matter from the LymanAlpha Forest, Phys. Rev. Lett. 126 , 071302 (2021), arXiv:2007.12705 [astro-ph.CO].
- [23] T. Kobayashi, R. Murgia, A. De Simone, V. Irˇ siˇ, and c M. Viel, Lymanα constraints on ultralight scalar dark matter: Implications for the early and late universe, Phys. Rev. D 96 , 123514 (2017), arXiv:1708.00015 [astroph.CO].
- [24] S.-R. Chen, H.-Y. Schive, and T. Chiueh, Jeans Analysis for Dwarf Spheroidal Galaxies in Wave Dark Matter, Mon. Not. Roy. Astron. Soc. 468 , 1338 (2017), arXiv:1606.09030 [astro-ph.GA].
- [25] D. J. E. Marsh and J. C. Niemeyer, Strong Constraints on Fuzzy Dark Matter from Ultrafaint Dwarf Galaxy Eridanus II, Phys. Rev. Lett. 123 , 051103 (2019), arXiv:1810.08543 [astro-ph.CO].
- [26] T. Broadhurst, I. de Martino, H. N. Luu, G. F. Smoot, and S. H. H. Tye, Ghostly Galaxies as Solitons of BoseEinstein Dark Matter, Phys. Rev. D 101 , 083012 (2020), arXiv:1902.10488 [astro-ph.CO].
- [27] M. Safarzadeh and D. N. Spergel, Ultra-light Dark Matter Is Incompatible with the Milky Way's Dwarf Satellites, ApJ 893 , 21 (2020), arXiv:1906.11848 [astroph.CO].
- [28] K. Hayashi, E. G. M. Ferreira, and H. Y. J. Chan, Narrowing the Mass Range of Fuzzy Dark Matter with Ultrafaint Dwarfs, Astrophys. J. Lett. 912 , L3 (2021),

arXiv:2102.05300 [astro-ph.CO].

- [29] N. Dalal and A. Kravtsov, Excluding fuzzy dark matter with sizes and stellar kinematics of ultrafaint dwarf galaxies, Phys. Rev. D 106 , 063517 (2022), arXiv:2203.05750 [astro-ph.CO].
- [30] M. G. Walker, M. Mateo, E. W. Olszewski, J. Pe˜arrubia, n N. W. Evans, and G. Gilmore, A Universal Mass Profile for Dwarf Spheroidal Galaxies?, ApJ 704 , 1274 (2009), arXiv:0906.0341 [astro-ph.CO].
- [31] J. D. Simon, The Faintest Dwarf Galaxies, ARA&amp;A 57 , 375 (2019), arXiv:1901.05465 [astro-ph.GA].
- [32] H.-Y. Schive, T. Chiueh, and T. Broadhurst, Soliton Random Walk and the Cluster-Stripping Problem in Ultralight Dark Matter, Phys. Rev. Lett. 124 , 201301 (2020), arXiv:1912.09483 [astro-ph.GA].
- [33] B. T. Chiang, H.-Y. Schive, and T. Chiueh, Soliton Oscillations and Revised Constraints from Eridanus II of Fuzzy Dark Matter, Phys. Rev. D 103 , 103019 (2021), arXiv:2104.13359 [astro-ph.CO].
- [34] K. A. Oman et al. , The unexpected diversity of dwarf galaxy rotation curves, Mon. Not. Roy. Astron. Soc. 452 , 3650 (2015), arXiv:1504.01437 [astro-ph.GA].
- [35] H. N. Luu, S. H. H. Tye, and T. Broadhurst, Multiple Ultralight Axionic Wave Dark Matter and Astronomical Structures, Phys. Dark Univ. 30 , 100636 (2020), arXiv:1811.03771 [astro-ph.GA].
- [36] H. Huang, H.-Y. Schive, and T. Chiueh, Cosmological simulations of two-component wave dark matter, Mon. Not. Roy. Astron. Soc. 522 , 515 (2023), arXiv:2212.14288 [astro-ph.CO].
- [37] M. Gosenca, A. Eberhardt, Y. Wang, B. Eggemeier, E. Kendall, J. L. Zagorac, and R. Easther, Multifield ultralight dark matter, Phys. Rev. D 107 , 083014 (2023), arXiv:2301.07114 [astro-ph.CO].
- [38] H.-K. Guo, K. Sinha, C. Sun, J. Swaim, and D. Vagie, Two-scalar Bose-Einstein condensates: from stars to galaxies, Journal of Cosmology and Astroparticle Physics 2021 , 028 (2021), arXiv:2010.15977 [astro-ph.CO].
- [39] J. Eby, M. Leembruggen, L. Street, P. Suranyi, and L. C. R. Wijewardhana, Galactic condensates composed of multiple axion species, Journal of Cosmology and Astroparticle Physics 2020 , 020 (2020), arXiv:2002.03022 [hep-ph].
- [40] A. Maleknejad and E. McDonough, Ultralight pion and superheavy baryon dark matter, Phys. Rev. D 106 ,
- 095011 (2022), arXiv:2205.12983 [hep-ph].
- [41] H. N. Luu, P. Mocz, M. Vogelsberger, S. May, J. Borrow, S. H. H. Tye, and T. Broadhurst, Nested solitons in two-field fuzzy dark matter, Mon. Not. Roy. Astron. Soc. 527 , 4172 (2024), [Erratum: Mon.Not.Roy.Astron.Soc. 528, 2882 (2024)], arXiv:2309.05694 [astro-ph.CO].
- [42] F. van Dissel, M. P. Hertzberg, and J. Shapiro, Core and halo properties in multi-field wave dark matter, Journal of Cosmology and Astroparticle Physics 2024 , 077 (2024), arXiv:2310.19762 [astro-ph.CO].
- [43] A. Pozo, T. Broadhurst, G. F. Smoot, T. Chiueh, H. N. Luu, M. Vogelsberger, and P. Mocz, Dwarf galaxies united by dark bosons, Phys. Rev. D 109 , 083532 (2024), arXiv:2302.00181 [astro-ph.CO].
- [44] V. Springel, N-GenIC: Cosmological structure initial conditions, Astrophysics Source Code Library, record ascl:1502.003 (2015).
- [45] A. Lewis, A. Challinor, and A. Lasenby, Efficient computation of CMB anisotropies in closed FRW models, Astrophys. J. 538 , 473 (2000), arXiv:astro-ph/9911177.
- [46] N. Aghanim et al. (Planck), Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641 , A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO].
- [47] I. De Martino, T. Broadhurst, S. H. H. Tye, T. Chiueh, and H.-Y. Schive, Dynamical Evidence of a Solitonic Core of 10 9 M ⊙ in the Milky Way, Phys. Dark Univ. 28 , 100503 (2020), arXiv:1807.08153 [astro-ph.GA].
- [48] F. Toguz, D. Kawata, G. Seabroke, and J. I. Read, Constraining ultra light dark matter with the Galactic nuclear star cluster, MNRAS 511 , 1757 (2022), arXiv:2106.02526 [astro-ph.GA].
- [49] R. Hlozek, D. Grin, D. J. E. Marsh, and P. G. Ferreira, A search for ultralight axions using precision cosmological data, Phys. Rev. D 91 , 103512 (2015), arXiv:1410.2896 [astro-ph.CO].
- [50] D. Grin, D. J. E. Marsh, and R. Hlozek, axionCAMB: Modification of the CAMB Boltzmann code, Astrophysics Source Code Library, record ascl:2203.026 (2022).
- [51] O. Hahn and T. Abel, MUSIC: MUlti-Scale Initial Conditions, Astrophysics Source Code Library, record ascl:1311.011 (2013).
- [52] H. N. Luu, Axion-Higgs cosmology: Cosmic microwave background and cosmological tensions, Phys. Rev. D 107 , 023513 (2023), arXiv:2111.01347 [astro-ph.CO].
Hoang Nhan Luu Philip Mocz Mark Vogelsberger Alvaro Pozo Tom Broadhurst S. -H. Henry Tye Tao Liu Leo W. H. Fung George F. Smoot Razieh Emami Lars Hernquist
Diverse dark matter haloes in Two-field Fuzzy Dark Matter
Fuzzy dark matter (FDM) is a compelling candidate for dark matter, offering a
natural explanation for the structure of diffuse low-mass haloes. However, the
canonical FDM model with a mass of $10^{-22}~{\rm eV}$ encounters challenges in
reproducing the observed diversity of dwarf galaxies, except for possibly
scenarios where strong galactic feedback is invoked. The introduction of
multiple-field FDM can provide a potential resolution to this diversity issue.
The theoretical plausibility of this dark matter model is also enhanced by the
fact that multiple axion species with logarithmically-distributed mass spectrum
exist as a generic prediction of string theory. In this paper we consider the
axiverse hypothesis and investigate non-linear structure formation in the
two-field fuzzy dark matter (2FDM) model. Our cosmological simulation with an
unprecedented resolution and self-consistent initial conditions reveals the
diverse structures of dark matter haloes in the 2FDM model for the first time.
Depending on the formation time and local tidal activities, late-time haloes
can host solitons of nested cores or solitons of one dominant species.