Nested solitons in two-field fuzzy dark matter

Dark matter as scalar particles consisting of multiple species is well motivated in string theory where axion fields are ubiquitous. A two-field fuzzy dark matter (FDM) model features two species of ultralight axion particles with different masses, $m_1 \neq m_2$, which is extended from the standard one-field model with $m_a \sim 10^{-22}\,{\rm eV}$. Here we perform numerical simulations to explore the properties of two-field FDM haloes. We find that the central soliton has a nested structure when $m_2 \gg m_1$, which is distinguishable from the generic flat-core soliton in one-field haloes. However, the formation of this nested soliton is subject to many factors, including the density fraction and mass ratio of the two fields. Finally, we study non-linear structure formation in two-field cosmological simulations with self-consistent initial conditions and find that the small-scale structure in two-field cosmology is also distinct from the one-field model in terms of DM halo counts and soliton formation time.


INTRODUCTION
The nature of dark matter (DM) remains elusive since its first evidence was found almost a century ago by Fritz Zwicky (Zwicky 1937).The cold dark matter (CDM) hypothesis, where DM is made of non-relativistic, collision-less particles, has successfully explained the anisotropies of the Cosmic Microwave Background and the largescale structure of the universe; see Aghanim et al. (2020) and Alam et al. (2021) for the most recent observational data from Planck and eBOSS, see Vogelsberger et al. (2020) for a review on stateof-the-art cosmological simulations of galaxy formation.From a microscopic perspective, the models of weakly interacting massive particles (WIMPs), where DM particles are presumably heavy and have interactions of electroweak scales with Standard Model particles, have been the leading candidates for CDM in the past decades.However, CDM is also facing several challenges due to inconsistencies between -body simulations and observations on small scales, including the "core-cusp problem" (de Blok 2010), the "too-big-tofail problem" (Garrison-Kimmel et al. 2014;Boylan-Kolchin et al. 2011), the diversity problem of rotation curves (Oman et al. 2015) -together with the null results of WIMP signals in many laboratory searches (Roszkowski et al. 2018).
Even though baryonic feedback, e. g. from supernova explosions or stellar winds (Governato et al. 2012;Di Cintio et al. 2014), may partially explain these CDM small-scale crises, other dark matter ★ E-mail: hnluu@connect.ust.hkmodels, e. g. warm dark matter (Lovell et al. 2012), fuzzy dark matter (Hu et al. 2000), or self-interacting dark matter (Spergel & Steinhardt 2000), have emerged as promising alternatives to CDM.Fuzzy dark matter (FDM), also known as wave dark matter (Schive et al. 2014a;Hui 2021) or ultralight axions (Marsh 2016;Hui et al. 2017;Ferreira 2021), is a particularly interesting DM candidate to resolve the small-scale problems but still maintain CDM predictions on large scales.FDM is typically comprised of extremely light bosons with mass   ∼ 10 −22 eV, which have de Broglie wavelengths of astronomical scales,  dB ∝ (  ) −1 ∼ kpc.
Therefore, in FDM cosmology, bound structures below kiloparsec scales are washed out due to an effective "quantum" pressure acting against gravity.For the same reason, FDM haloes form a central flat core called soliton (Schive et al. 2014a;Schwabe et al. 2016;Mocz et al. 2017) rather than the cuspy Navarro-Frenk-White (NFW) profile of CDM (Navarro et al. 1996).A typical density profile of FDM haloes also features sophisticated interference granules of density fluctuations (Schive et al. 2014a;Mocz et al. 2017) around the central core.However, the simplest model of FDM, where   is the only parameter, has been tightly constrained by observational data from the Lyman- forest (Armengaud et al. 2017;Iršič et al. 2017;Kobayashi et al. 2017), ultrafaint dwarf galaxies (Marsh & Niemeyer 2019;Hayashi et al. 2021;Dalal & Kravtsov 2022), and strong lensing (Powell et al. 2023).The current consensus is that a single FDM species of mass 10 −22 eV is unlikely to account for the total dark matter budget while still addressing the small-scale issues of CDM at the same time.
Several extensions of the most simplified FDM model have been widely explored in literature.A few notable studies are selfinteracting FDM (Mocz et al. 2023) where the quartic self-coupling term is included, vector fuzzy dark matter (Amin et al. 2022) where FDM is a higher-spin field, and mixed cold and fuzzy dark matter (Schwabe et al. 2020).From a theoretical perspective, the ultralight mass of the FDM field can be generated in the context of string compactifications.It is expected that there may exist more than one species of ultralight axions with masses spanning many orders of magnitudes as in the string axiverse scenario (Svrcek & Witten 2006;Arvanitaki et al. 2010).Thus, the multiple-field FDM model is naturally motivated by string theory.Along this direction, the multiple-field FDM model with a specific focus on the two-field case has been proposed in Luu et al. (2020).
We previously found that the ground state of a time-independent two-field FDM system is a nested configuration of two solitons, each from one component field, which is distinct from the usual soliton in the one-field FDM model.Ever since our work, there has been increasing interest in the two-field FDM model, starting with Eby et al. (2020) and Guo et al. (2021) where the authors studied the semianalytical profile of two-field solitons.More recently, Huang et al. (2023) performed the first two-field cosmological simulations with CDM initial conditions; Gosenca et al. (2023) investigated stellar heating in two-and higher-field FDM haloes; Glennon et al. (2023) simulated soliton collision of two-field FDM with self-interactions and inter-field interactions; Jain et al. (2023) studied Bose-Einstein condensation in the kinetic regime.These studies showed that the phenomenology of the two-field FDM model can be complex, novel, and compelling.Importantly, it has the potential to retrieve the initial promises of the one-field FDM model.
In this work, we carry out non-cosmological and cosmological simulations of the two-field FDM model to study the final state of virialized haloes.Our main goal is to verify whether the nested soliton, previously predicted in Luu et al. (2020) can emerge and become stable in numerical simulations.The paper is organized as follows.In Sec. 2, we review the theoretical background of the twofield Schrödinger-Poisson equations.Sec. 3 describes the numerical method for our simulations and ground-state solver in detail.Sec. 4 and Sec. 5 set up several simulations and then report the main results and new findings of our paper.Finally, conclusions and potential directions for future work are discussed in Sec. 6.

THEORETICAL BACKGROUND
In this section, we introduce the time-dependent and timeindependent Schrödinger-Poisson equations of the two-field FDM model in Sec.2.1 and Sec.2.2, respectively.

Two-field Schrödinger-Poisson equations
In the non-relativistic limit, the minimal (non-self-interacting) system of two-field FDM with different masses is governed by gravitationally coupled Schrödinger-Poisson equations.In physical coordinates in the Newtonian gauge, these equations are given by where Φ is the gravitational potential,  1 and  2 are the wavefunctions of the FDM field with masses  1 and  2 , respectively.The first field  1 is always chosen as the lighter FDM field in our convention throughout this study.The density of each field is given by   = |  | 2 , where  = 1, 2. The total density is obtained by summing individual densities,  =  1 +  2 , while ρ is the corresponding mean of the total density.For convenience, we define a parameter to quantify the ratio of the heavy-field mean density to the total density as and   denotes the total mass of each field enclosed within a volume with length  box along each spatial dimension.As the two-field Schrödinger-Poisson equations (1) are highly non-linear due to gravitational couplings, numerical simulations are mandatory to study the system's dynamics.Simulation techniques for the Schrödinger-Poisson equations will be reviewed later in Sec.3.1.

Time-independent Schrödinger-Poisson equations
Some prominent features of FDM haloes, such as the stable groundstate solution, can be found in an approximate limit.Once the field dynamics are stabilized, the time-dependent part of the FDM fields can be factorized out,   (x, ) = exp(−  /ℏ)   (x), where   is the "eigenenergy" of each field.Let us also assume that the system is spherically symmetric,   (x) =   (), Φ(x) = Φ(), and omit the Hubble flow as well as the mean density.The time-independent Schrödinger-Poisson equations are then derived as follows (3) The ground-state solution of this system is called a soliton, which is a cored density profile that satisfies the following boundary conditions and Φ(0) can be arbitrarily chosen.Integrating the density of each field yields the corresponding total mass of the soliton In practice, the set of equations (3) can be further simplified thanks to the scaling symmetry of the Schrödinger-Poisson equations.The Schrödinger-Poisson equations are unchanged under a transformation of variables with parameter  such that (Ji & Sin 1994) This scaling relation suggests that Eqs.
(3) can be re-scaled into where a new set of dimensionless variables has been introduced with a tilde sign.They are related to the original variables by Here,  is a free scaling parameter, corresponding to  in (6). is the normalization factor of the soliton mass, ∫ | ψ | 2 r2  r =  , /, and  represents an effective FDM mass (which can, however, be freely chosen), which is usually set to  = √  1  2 in our computations.

NUMERICAL METHOD
In this section, we review the pseudo-spectral method used for our simulations in Sec.3.1, then provide detailed guidelines to solve for the ground state of two-field equations in Sec.3.2.

Pseudo-spectral solver
The dynamical two-field FDM system in Eqs.
(1) can be simulated using the pseudo-spectral method extended for two fields with different masses.The simulation volume is a cubic box discretized on a fixed Cartesian grid with uniform resolution.Periodic boundary conditions are implicitly imposed in the spectral solver.The wavefunction of each field is evolved unitarily for each timestep Δ as follows: • Step 1: Evolve  1 and  2 with the potential "kick" operator by half a timestep • Step 2: Evolve  1 and  2 with the kinetic "drift" operator by a full timestep in Fourier space • Step 3: Calculate the total density and update the potential (also in Fourier space) • Step 4: Repeat Step 1 with the updated potential.
Here "fft" and "ifft" represent the Fast Fourier Transform and its inverse.We employ the FDM module of the magneto-hydrodynamics code AREPO (Springel 2010;Pakmor et al. 2016;Weinberger et al. 2020) originally implemented by the authors of Mocz et al. (2017) and modify it accordingly for two-field FDM simulations.
Since the pseudo-spectral algorithm converges at second order in time and exponentially in space, the heavy-field (smaller) soliton is fully resolved even when it overlaps only a few cells along each dimension.On the one hand, the spacing requirement of the grid, Δ, needs to cover the finest features of the heaviest field (Mocz et al. 2018;May & Springel 2021) where  ,max denotes the maximum velocity of the corresponding field.On the other hand, the timestep criterion needs to be satisfied for the lightest field (Schwabe et al. 2016;May & Springel 2021), i. e., These two criteria often require simulations with very high temporal and spatial resolution to resolve the inner core of both fields.Therefore, multiple-field simulations are much more computationally challenging than one-field simulations.

Two-field ground state solver
The dimensionless equations ( 7) are particularly simple to solve and useful to scale up to different physical systems.For instance, in the one-field FDM model, the scaling parameter  can be chosen as the FDM mass   and we obtain a completely scale-free set of equations where the field indices have been omitted.As the above equations are independent of FDM parameters, the ground state corresponds to a unique solution of ψ.This solution can be freely scaled up (via the parameter ), yielding the soliton profile approximated by (Schive et al. 2014a) where  c is the core radius defined as the radius where the core density drops by half, ( c ) = (0)/2.The total mass of each soliton is given approximately by   ≃  ,0 ( c /kpc) (  /10 −22 eV) 2 where  ,0 ≡ 2.2 × 10 8 M ⊙ .We note that this profile only describes the ground-state solution of a one-field FDM system.In contrast, the ground-state solution of a two-field FDM system is not well approximated by this profile.As such, we will refer to the ground state of Eqs. ( 14), and equivalently the profile (15), as one-field soliton, or simply soliton.The ground state of Eqs.
For the two-field equations in (7), we cannot eliminate the FDM masses with scaling relations.Thus, the ground-state solution is dependent on at least two parameters: the FDM mass ratio  2 / 1 and the central wavefunction ratio defined as since ψ1 (0) can always be fixed to a constant, e. g., ψ1 (0) = 1.We apply the shooting method to solve these equations: (i) pick a not-too-large endpoint rend and arbitrarily choose some initial values for Ẽ ; (ii) solve ( 7) as an initial condition problem where ψ1 (0) = 1, ψ2 (0) =  2 ; (iii) adjust Ẽ1 and Ẽ2 simultaneously such that ψ ( rend ) → 0 with multi-dimensional Newton-Raphson algorithm; (iv) repeat step (i) to (iii) until rend is large enough.There may exist several excited states of different energies satisfying regular boundary conditions.Thus, we must impose the condition ψ () > 0 ∀ to guarantee the ground-state solution.We list in Tab. 1 some eigenenergies corresponding to the ground state of the two-field systems for different values of  2 and  2 as an example.Although these steps are straightforward, we cannot apply them to solve the two-field system with a high mass ratio,  2 / 1 ≫ 1.In these limits, we need to tune Ẽ2 extremely precisely compared to Ẽ1 so that the solution does not blow up.As Ẽ2 reaches floating-point precision, the shooting method stops converging without arbitraryprecision arithmetic.The high-mass limit is, however, an ideal condition to apply the flat density approximation, where the density distribution of the light field is treated as constant with respect to the heavy field assuming that their cores are very different in sizes, r,2 ≪ r,1 .This allows the use of what are effectively the onefield equations ( 14), with an additional constant external source.The approximation procedure also follows four steps: • Step 1: Solve for the heavy wavefunction with a constant density distribution of the light field as an external source, setting  =  1 and  = M (arbitrary), i. e.
where we have adopted an approximation ψ1 ( r) ≃  −1 2 .• Step 2: Rescale the variables by setting  =  2 and  = • Step 3: Solve for the light wavefunction with an additional source of the heavy field just found: • Step 4: To return to our "usual choice" of  = √  1  2 , rescale the variables again by additionally setting We then obtain an approximate solution of Eqs. ( 7).
The eigenenergies are not crucial in this method, but they can be straightforwardly retrieved from the scaling relations: These values provide a reasonable estimate of the exact energies (see Tab. 1) and they can serve as an initial guess for the precise shooting method.Figure 1 compares the ground-state wavefunctions solved by both methods for Scenario B with  2 = 0.5.Even though the mass ratio is not high in this case, the flat-density approximation (dashed curves) yields virtually indistinguishable profiles from the precise ones (solid curves)., are randomly placed in the simulation box.These solitons are not the same as the two-field soliton, which may (or may not) form as the final state in each scenario.The last column shows whether the two-field soliton is found at the end of each simulation.Regardless of the method used, we can restore the two-field soliton density by computing the scaling factor  based on the relation between dimensionless and physical variables.More specifically, we require

Scenario
to fit a two-field FDM soliton with central density Once  is determined, the physical radius and the physical wavefunction of both fields can be derived according to (8).

SOLITON COLLISION WITH TWO-FIELD FDM
We perform non-cosmological simulations to investigate the features of two-field haloes.Initially, the soliton cores with density profiles given by ( 15) are randomly placed in a box of length The bottom panels are the corresponding radial profiles of each field.The red dashed curve in the left panel is shown as the sum of the light (blue dashed) and heavy (green dashed) one-field solitons.The FDM mass of the heavy field is five times that of the light field,  2 = 5 1 = 4 × 10 −22 eV, while its abundance is four times less than the light field,  2 = 0.25, as indicated in the figures.We distinguish the core radius of the best-fit one-field soliton (dashed), denoted as FDM, and two-field soliton (solid), denoted as 2FDM, as the former is overshooting the simulation profile of the light field.In contrast, the core radius of the heavy field is almost identical in both profiles (with dashed and solid green curves almost overlapping).The orange dashed curve denotes the best fit of a NFW profile with the scale density  ,0 ∼ 9 × 10 7 M ⊙ kpc −3 and the scale radius   ∼ 1.85 kpc to the outer halo of the total field.
box = 50 kpc with zero velocities.We then let them gravitationally interact, merge, and form a virialized halo.The final state is captured at time  = 2.4−2.6 Gyr, when the two-field halo has become relatively stationary.Table 2 summarizes the setup of the soliton collision simulations in the following sections.

Nested soliton profile
Firstly, we collide 15 solitons of the light field,  1 = 8 × 10 −23 eV, with 45 solitons of the heavy field with  2 = 4 × 10 −22 eV in a simulation box with a resolution of 1024 3 cells.The core radius of the light-field solitons is manually set to  ,1 = 2 kpc.The core radius of the heavy-field solitons is inferred from  ,1 to match the mass of the light-field ones,  ,2 =  ,1 ( 1 / 2 ) 2 , so that the mean density of the heavy field accounts for 25% of the total density in the system, i. e.,  2 = 0.25.We refer to this simulation as Scenario A.
Figure 2 shows slices of the density distribution and corresponding radial profiles that are averaged over angular directions in Scenario A. The origin is chosen at the location with maximum total density.
As expected, the halo of the total field and the individual components show a central core surrounded by interference granules in the simulation.Most notably, the core of the total field exhibits a nested structure of two constituent solitons with a transition around 0.16 kpc.This nested profile is fitted well with the ground-state soliton (red solid curve) of the time-independent two-field Schrödinger-Poisson equations (3).Similarly, the predicted two-field soliton profiles of the light field (blue solid curve) and heavy field (green solid curve) are in excellent agreement with the inner cores of the associated simulation fields, although the heavy-field soliton is only resolved by three cells.The outer halo follows an NFW-like profile, similar to the one-field FDM model, as expected.
The one-field profile of the light field (blue dashed curve) with a larger core radius does not fit the simulation data, unlike its two-field counterpart.As this profile is determined by solving the Schrödinger-Poisson equations for the light field only, the corresponding core radius is noticeably wider due to the absence of gravitational attraction from the heavy-field soliton.On the other hand, the one-field and two-field soliton profiles of the heavy field are virtually identical since it only experiences gravitational effects from an almost flat and low-density distribution of the light field, i. e.,  2 (0) ≫  1 (0).The one-field profile of the total field (red dashed curve), defined as the sum of two component one-field solitons, is not compatible with the simulation data in a similar manner to the light field.In summary, the two-field halo profile contains a nested soliton with a visible transition from one core to another and a second transition from the soliton to the remaining NFW-like halo.Yet, the soliton-soliton transition is not always visible due to time-dependent fluctuations of the heavy-field granules which induce a slight density excess over the light-field soliton.The soliton-halo transition is found at  ∼ 3 − 3.5 , for the light and heavy field, which agrees with previous studies (Schive et al. 2014b;Mocz et al. 2017).As for the total field, although the soliton-halo transition can be seen around 0.5 kpc, it is generally ill-defined when the light-field soliton blends with the interference granules of the heavy field as mentioned above.

Soliton formation and instabilities
Even though we have seen the existence of a stable nested soliton as the ground-state solution of two-field Schrödinger-Poisson equations, it may not always be able to form with any value of the mass fraction.Recently, one particular cosmological simulation of twofield FDM (Huang et al. 2023) with  2 = 3 1 = 3 × 10 −22 eV and  = 0.25 has pointed out that the heavy-field soliton cannot form under deep gravitational potential fluctuations induced by the light field.We are, therefore, motivated to study the formation of two-field solitons in Scenario B as follows.
We collide 10 solitons of each field with  1 = 8 × 10 −23 eV and  2 = 1.6 × 10 −22 eV in a simulation box with a resolution of 512 3 cells.The heavy-field solitons are placed concentrically with the light-field ones to enhance collision efficiency.The total mass is fixed,  total =  1 +  2 = 3 × 10 9 M ⊙ .The heavy-field density ratio varies from 30% to 70%,  2 = {0.3,0.5, 0.7}, in three simulations.Thus, the core radii of the initial solitons need to be adjusted in each simulation to satisfy the density ratio required where  , denotes the number of initial solitons from each FDM species and  ,1 =  ,2 = 10 in this case.Figure 3 (upper row) shows radial profiles of the light and heavy field in different simulations with  2 = 2 1 = 1.6 × 10 −22 eV and different mass fractions  2 .We observe distinctive results in each simulation, which are discussed in order as follows.When the heavy field is dominant,  2 = 0.7, the two-field solitons fit the numerical simulation well for both component fields.The one-field profile of the heavy field is still aligned with its two-field profile, but a similar alignment is not seen in the case of the light field.These findings completely agree with the previous results in Scenario A. When the two fields are comparable in abundance,  2 = 0.5, the halo profile of each field still respects the corresponding two-field ground state.However, both one-field profiles deviate from the two-field ones.
When the light field is dominant,  2 = 0.3, the inner core of this field is consistent with both one-field and two-field solitons.In contrast, the inner profile of the heavy-field halo fails to fit either of these solitons, which implies that the second soliton cannot form in this case.Our results independently confirm that there exists a limit of a mass fraction where the formation of the heavy-field soliton is inhibited, as pointed out by Huang et al. (2023).
We estimate this threshold by running additional simulations which are set up identically as above with  2 varying continuously from 0.36 to 0.54.We then calculate the goodness-of-fit of the heavyfield soliton which describes how well it fits the corresponding simulation data, defined as where  sol,2 and  sim,2 are the two-field soliton and the simulation profile of the heavy field, respectively;  fit is the number of simulation data points used counting from the halo center since only these points are expected to fit the soliton.The inset of the upper right panel in Fig. 5 shows − 2 with respect to  2 for  fit = 3 (orange) and  fit = 4 (blue).A sudden drop from  2 ≃ 0 implies no soliton found in the simulation data.It is, hence, straightforward to infer the mass fraction threshold to be approximately  ≃ 0.44 for Scenario B.
On the other hand, we must highlight two important properties of two-field FDM haloes that are found in this set of simulations.Firstly, the one-field profiles fail to describe two-field solitons, especially for the light field.These discrepancies are particularly prominent when the heavy field becomes more dominant ( 2 = 0.5 compared to  2 = 0.7 in Scenario A) or when the heavy-field mass is closer to the light-field mass ( 2 = 2 1 in Scenario B compared to  2 = 5 1 in Scenario A).As the heavy-field soliton is more extended and massive in these two limits, its gravitational pull on the light field is also larger and vice versa.The case  2 = 0.7 provides an example when the light soliton is so compressed that its core is no longer flat while the case  2 = 0.5 has the solitons of both fields being squeezed under the gravitational potential of each other.Thus, fitting solely one-field solitons to the numerical simulation data could lead to the non-detection of solitons, which is a misidentification of the real twofield ground state.Secondly, the threshold of a mass fraction where the heavy-field core becomes disrupted is determined by several characteristics of the two-field FDM system.The nested soliton can form with a mass fraction of  2 = 0.25 as seen in Scenario A, while the same phenomenon does not happen in the case  2 = 0.3 of Scenario B with a higher mass fraction.
However, comparing soliton formation in Scenario A and Scenario B is not completely fair as there are other varying factors, in addition to  2 and  2 , in these two scenarios such as the total mass  (or equivalently the mean density ρ) and the initial conditions (core radius, number of colliding solitons).Therefore, we set up Scenario C similar to Scenario B but with the heavy-field mass varying,  2 = {4.0,4.8, 5.6} × 10 −22 eV, and the density fraction fixed to  2 = 0.3 in three simulations.The core radius of initial solitons also changes according to Eq. ( 23).
The final profiles are shown in the lower row of Fig. 3, where the heavy-field soliton can only form with  2 = 7 1 = 5.6 × 10 −22 eV, but not with  2 = 5 1 = 4 × 10 −22 eV and  2 = 6 1 = 4.8 × 10 −22 eV, for a fixed  2 = 0.3.The only difference in the last case with  2 = 7 1 is that its initial state at  = 0 consists of much denser solitons (see Eq. 15 and Tab. 2) which enables the system to reach gravitational equilibrium faster.Once the stable soliton of the heavy field has formed, it is no longer susceptible to tidal disruption from density fluctuations of the light field.Thus, it is likely that the formation of a nested soliton in the two-field FDM model is viable if the heavy-field soliton formation has completed relatively earlier than the light-field one.This condition can be realized even if the heavy-field density is sub-dominant as long as the two FDM masses differ by orders of magnitude,  2 ≫  1 .

TWO-FIELD FDM IN COSMOLOGY
We perform DM-only cosmological simulations of the two-field FDM model with self-consistent initial conditions.

Simulation setup
Four simulations are conducted in a comoving volume that has a length of  box = 1.4 ℎ −1 Mpc with a spatial resolution of 1024 3 cells.The four simulations are set up as follows: (i) CDM via body dynamics described by Vlasov-Poisson equations, (ii) one-field FDM with   =  1 = 10 −22 eV for the light field only (i.e.  2 = 0), (iii) one-field FDM with   =  2 = 2 × 10 −22 eV for the heavy field only (i.e.  2 = 1), (iv) two-field FDM with  1 and  2 as given above and  2 = 0.5.We note that the mass fraction of two fields in this case is defined as the ratio of cosmological density parameters, i. e., where Ω a, denotes the present-day density fraction of each FDM field.Initial conditions are generated from the same random seed for comparison at  = 127 using N-GenIC from Springel (2015) with a two-field matter power spectrum calculated by the publicly available Boltzmann code aHCAMB from Luu (2023) 1 .The heavy-field mass is assumed to be  2 = 2 × 10 −22 eV, two times heavier than the light-field mass of  1 = 10 −22 eV.Cosmological parameters are chosen as constrained by Planck (Aghanim et al. 2020), i. e., Ω m = 0.314, Ω Λ = 0.686, ℎ = 0.6732,   = 0.966.One exception is the primordial spectrum amplitude, which is boosted to   = 7.5 × 10 9 to enhance small-scale fluctuations.Thus, the simulation box is more representative of an overdense region than an average cosmological volume in the two-field FDM universe.The physical system in these simulations evolved until redshift  ≃ 3.5, after which halo cores become too small in the comoving grid to resolve.

Results and Discussion
Figure 4 provides a glimpse of the evolution of DM clustering in cosmological simulations of four different DM models, as set up in Sec.5.1.Each panel shows the projected density distribution of DM structures which consists of various DM haloes connected by cosmic filaments.As expected from previous FDM simulations (Mocz et al. 2019;May & Springel 2023), small-scale fluctuations are smoothed out in all FDM models, compared to the scale-independent clustering in the CDM model.The number of FDM haloes is, thus, strongly suppressed due to a sharp cut-off below the Jeans scale in the initial matter power spectra and the "quantum pressure" from FDM dynamics that counteracts gravity.In terms of one-field simulations, the heavy-field structures are more abundant and form earlier than the light-field ones throughout cosmic history.The first structure of the heavy field emerges at  ≃ 6.4, while the first light-field halo is not seen until  ≃ 3.5.We also identify four heavy-field haloes in total compared to only one light-field halo at the final redshift.
On the other hand, the two-field FDM cosmology (with  2 = 0.5) looks like an intermediate state of the one-field models.Two haloes have formed by the end of our simulation with the first one becoming visible at  ≃ 4.8.Zooming into this two-field halo in Fig. 5 (upper panel) reveals the soliton core and wave interference patterns on scales of the de Broglie wavelength of both fields.Since the total field is a superposition of two-component fields, its core and granular structures are slightly denser and fuzzier.We found that the radial profile of the central soliton exactly follows the two-field ground state (lower panel), as in the non-cosmological simulations, even though the solitons are just barely resolved.The nested structure is not apparent in this case, owing to the low mass ratio of two fields, i. e.,  2 = 2 1 .More importantly, the haloes in our two-field simulation consist of solitons of both light and heavy fields.These findings show that the light-field soliton formation in the two-field simulation is greatly enhanced, as no soliton appeared until the last redshift in the onlylight-field simulation.The light-field density likely started clustering faster thanks to deep gravitational potentials induced by overdense regions of heavy-field haloes which had formed beforehand in the two-field cosmology.Therefore, in the limit of high mass ratios,  2 ≫  1 , we speculate that heavy-field soliton formation may not be susceptible to the oscillating potential of the light-field soliton in a realistic cosmological context, even with low-density fractions  2 between two fields, as its clustering commences significantly earlier in the cosmic history.On the other hand, late-time DM haloes may predominantly incorporate solitons of both fields in this context, which poses a new challenge to explain the diverse core size of dwarf galaxies by multiple FDM species (Luu et al. 2020;Pozo et al. 2023), where soliton formation of each FDM species should be mostly independent.

CONCLUDING REMARKS
As evidence has begun to build up against the one-field FDM model with mass   ∼ 10 −22 eV, the two-field FDM model with different masses  1 and  2 offers an interesting and plausible alternative for the FDM paradigm.In this study, we have conducted state-of-theart simulations to study structure formation in the two-field FDM model.With simulations of soliton collision, we discovered that the virialized haloes feature a nested soliton and an NFW-like outer profile.The nested soliton is a stable configuration of two concentric cores of different sizes, which can be solved as the ground state of two-field time-independent Schrödinger-Poisson equations.
Unlike the one-field soliton profile, component solitons in the two-field model are more compact due to mutual gravitational attraction.Depending on the density ratio of the component fields, soliton formation may not occur in case the heavy field is excessively dominated by the light field.Deep gravitational potentials induced by soliton oscillations may be responsible for this phenomenon (Huang et al. 2023).However, this threshold of soliton instabilities is also sensitive to initial conditions and mass differences between the two fields.A general pattern is that as the second field becomes heavier, the nested soliton can form at a lower density fraction.
We also had a preliminary look into the large-scale structure of the two-field FDM model.Despite the limited resolution, two-field soliton profiles are found in our cosmological simulation.Since there are three viable profiles in two-field FDM cosmology: two-field nested soliton, one-field light soliton, and one-field heavy soliton, one might expect to find these different types of haloes in a larger volume, which would explain the separate hierarchy of core sizes in spheroidal dwarf galaxies and ultrafaint dwarf galaxies (Luu et al. 2020).Unfortunately, every halo in our small-box simulation only contains a two-field soliton, so we leave this interesting problem to a future study.
The two-field model can potentially alleviate several tensions between FDM predictions and observational data.Hayashi et al. (2021) argued that ultrafaint dwarf galaxy cores as FDM solitons imply   ≳ 10 −21 eV, incompatible with the 10 −22 eV mass for dwarf spheroidal galaxies, which can be precisely explained by the heavy field.Dalal & Kravtsov (2022) suggested that orbiting stars are perturbed as they interact with the gravitational potential of density granules in FDM haloes, which is particularly prominent for lighter FDM fields, so the mass range   ≲ 3 × 10 −19 eV is strongly excluded.As indicated by Gosenca et al. (2023), the presence of one (or more) heavier field(s) helps to smooth out the density fluctuations in FDM haloes, hence relaxing the bound.Lastly, FDM suppresses the linear matter power spectrum at small scales, which alters the structure of the intergalactic medium (IGM) on the scales probed by the Lyman- forest of early-time quasars.By relating the one-dimensional flux spectrum with the FDM power spectrum, Armengaud et al. (2017) and Iršič et al. (2017) constrained   ≲ ×10 −21 eV.Kobayashi et al. (2017) also derived similar constraints for the mixed DM model for any amount of FDM more than 30% of the total DM.Although there still exist several uncertainties of IGM gas physics assumed in these studies, which should be better understood (Hui et al. 2017), we expect a sufficiently large density fraction of the heavy field will resolve the Lyman- constraint.Nevertheless, more detailed studies on the two-field FDM model in theory, phenomenology, and simulation are needed to settle these issues.
Although our work provides a fairly comprehensive study of twofield FDM regime, the results are restricted to the low mass ratio  2 / 1 = 2 − 7 and the minimal model with uncoupled fields.It is, therefore, interesting to explore different variations of the two-field FDM in a broader context.For instances, Schwabe et al. (2020) have extensively studied the mixed DM model and pointed out that the soliton cannot form in cosmological simulations with FDM constituting less than 10% of the total DM.Provided CDM is effectively treated as the heavy FDM field with much higher mass than the light field,  2 ≫  1 , there seems to exist a density threshold that inhibits the soliton formation of the light field, which deserves more investigation.Another example is Cyncynates et al. (2022) who proposed a coupled "friendly" two-axion model with a cosine potential.As long as the heavy-field mass is sufficiently close to the light-field one,  2 / 1 ∼ 0.75, the coupling triggers autoresonance mechanism which synchronizes the background oscillations of both fields.Furthermore, small-scale density fluctuations are enhanced in this model, which can alleviate the Lyman- constraints mentioned above, a feature shared by the extreme axion framework with an initial large displacement angle & Chiueh 2017; Arvanitaki et al. 2020).
At the end of the day, since the nested soliton is exclusive to two (or multiple) species of FDM, potential observations of this profile in future galaxy surveys will be a smoking gun signature of the multiple-axion scenario.

Figure 2 .
Figure2.The top panels show slices through a snapshot of the total, light, and heavy fields (from left to right) of a virialized two-field FDM halo at  ∼ 2.6 Gyr.The bottom panels are the corresponding radial profiles of each field.The red dashed curve in the left panel is shown as the sum of the light (blue dashed) and heavy (green dashed) one-field solitons.The FDM mass of the heavy field is five times that of the light field,  2 = 5 1 = 4 × 10 −22 eV, while its abundance is four times less than the light field,  2 = 0.25, as indicated in the figures.We distinguish the core radius of the best-fit one-field soliton (dashed), denoted as FDM, and two-field soliton (solid), denoted as 2FDM, as the former is overshooting the simulation profile of the light field.In contrast, the core radius of the heavy field is almost identical in both profiles (with dashed and solid green curves almost overlapping).The orange dashed curve denotes the best fit of a NFW profile with the scale density  ,0 ∼ 9 × 10 7 M ⊙ kpc −3 and the scale radius   ∼ 1.85 kpc to the outer halo of the total field.

Figure 3 .
Figure 3. Radial density profiles from several non-cosmological simulations calculated at  ∼ 2.4 Gyr.(Upper row) Profiles of the light,  1 = 8 × 10 −23 eV, and heavy field,  2 = 1.6 × 10 −22 eV, in three simulations with varying mass fractions  2 ≡  2 / total (Scenario B).The inset of the last panel shows the goodness-of-fit of the heavy-field profile (see definition in (24)) in other simulations with different  2 and at the same moment.(Lower row) Profiles of the light and heavy fields in similar simulations with varying heavy-field mass  2 and a fixed mass fraction  2 = 0.3 (Scenario C).The heavy field cannot reach a stable ground state in some simulations where neither the two-field soliton nor the one-field soliton fits the simulation profile.

Figure 4 .
Figure 4. Snapshots of projected (comoving) density in four simulations at  = 15.1, 6.4, 4.8, 3.5 to illustrate structure formation of different DM models throughout cosmological history.These values of redshift are chosen to fully capture the onset of soliton formation in each FDM simulation.The simulation volume has a comoving side length of  box = 1.4 ℎ −1 Mpc.The white box at  = 4.8 in the bottom row marks the first structure in a two-field FDM cosmology with  2 = 2 1 = 2 × 10 −22 eV and  2 = 0.5.The close-up view of this region is shown in Fig. 5.

Figure 5 .
Figure5.The upper panels show close-up slices through the center of the first two-field halo formed at  = 4.8 as seen in Fig.4(all quantities in comoving coordinates).The lower panel shows the radial profile of this halo in comoving coordinates, including a two-field soliton at the center (red) and an outer NFW-like profile.The profiles of the light field (blue) and the heavy field (green) are also displayed for comparison.Their solitonic cores are just marginally resolved due to the limited resolution of our simulation.