Cosmological structure formation and soliton phase transition in fuzzy dark matter with axion self-interactions

ABSTRACT

1 INTRODUCTION

Ultra-light bosons continue to be a popular candidate for the dark matter in our Universe (Hu, Barkana & Gruzinov 2000; Guzmán & Ureña-López 2003; Hui et al. 2017; Mocz et al. 2019; Burkert 2020; Niemeyer 2020; Hui 2021). The so-called fuzzy dark matter (FDM) model postulates a particle mass of m ∼ 10⁻²² eV, which introduces wave dynamics in the dark matter on the de Broglie wavelength λ_dB ∼ 1 kpc – the scale of galaxies. The arising quantum ‘pressure’ (really a pressure tensor) suppresses small-scale power in the initial dark matter power spectrum (Hu et al. 2000), modifies the halo mass function (Schutz 2020), and creates soliton cores at the centres of dark matter haloes (Schive, Chiueh & Broadhurst 2014a; Marsh & Pop 2015). Solitons are quasi-stable cored objects with total mass scaling inversely with radius, unique to the FDM model. FDM has seen a rise in direct numerical simulations that investigate non-linear and small-scale features of the model (Schive et al. 2014a; Mocz et al. 2017; Du et al. 2018; Mocz et al. 2019; Laguë et al. 2020; Mocz et al. 2020; Schwabe et al. 2020; Veltmaat, Schwabe & Niemeyer 2020; Li, Hui & Yavetz 2021; May & Springel 2021; Nori & Baldi 2021). A challenge for the FDM model continues to understand whether an ultra-light particle mass can simultaneously predict the Lyman α forest power spectrum extracted from high-redshift quasars (Iršič et al. 2017; Nori et al. 2019), as well as explain the core sizes of satellite galaxies (Burkert 2020; Safarzadeh & Spergel 2020). This is a Catch-22 problem (Davies & Mocz 2020) of sorts, in which smaller boson masses lead to larger, less dense cores, but also less structure in the Lyman α forest. Dalal & Kravtsov (2022) use sizes and stellar kinematics of ultra-faint dwarf galaxies to place a strict lower limit of m > 3 × 10⁻¹⁹ eV in the simple FDM model, arguing the one-parameter family of soliton solutions at the centre of haloes cannot fit the observational data at lower particle masses (neglecting dynamical heating).

This Catch-22 may be resolved by the introduction of a second relevant scale, determined by the way that FDM particles interact with one another. This could arise naturally in one of the main candidate models for the FDM boson: the hypothetical axion arising from the symmetry breaking, needed to solve the strong CP problem (Peccei & Quinn 1977; Weinberg 1978). In this model, the axion would have a decay constant (or symmetry-breaking scale) f associated with it, which would give rise to an attractive self-interaction (SI; Desjacques, Kehagias & Riotto 2018; Arvanitaki et al. 2020). Such an ultra-light axion may constitute a considerable fraction of the present-day critical density of the Universe (e.g. Marsh 2016; Hui et al. 2017; Desjacques et al. 2018):

For fiducial values m ∼ 10⁻²² eV and f ∼ 10¹⁷ GeV, the attractive SI is tiny: the dimensionless strength of the quartic coupling m²/f² ∼ 10⁻⁹⁶, and hence the attractive SI has so far been ignored in most numerical simulations. Despite this tiny value, the analytical findings of Desjacques et al. (2018) indicate that the cosmic web is influenced by a small, non-vanishing self-coupling among ultra-light axions. Desjacques et al. (2018) show that attractive SI can have a significant impact on the stability of cosmic structures at low redshift, including filaments and soliton cores. A noticeable effect on cosmological scales is likely to be seen at f ≲ 10¹³ GeV. Other analytical studies have also indicated that attractive SI would only allow solitons to remain stable below some critical maximum mass (Vakhitov & Kolokolov 1973; Chavanis 2011, 2016). Below that mass, solitons are in a dilute phase (Chavanis & Delfini 2011). Above that mass, the solitons collapse and form dense solitons (Braaten, Mohapatra & Zhang 2016), which are stabilized by higher order repulsive terms in the expansion of the SI potential (Eby et al. 2016; Chavanis 2018).¹ The collapse of the solitons may be accompanied by a sort of ‘explosion’ (a burst of relativistic axions), leading to a bosenova (Levkov, Panin & Tkachev 2017). The bosenova phenomenon occurs in the case of a relatively strong SI |$f\lt M_P\sim 10^{19}\, {\rm GeV}$|⁠. In certain regimes, not relevant here, it is necessary to take general relativity into account and the collapse rather leads to a black hole (Helfer et al. 2017).

The attractive SI arises as a leading term from the scalar field dark matter potential, which may have various forms. A number of these have been studied by Boltzmann-type cosmological simulations, which have been found to have excess power in the dark matter power spectrum compared against cold dark matter (CDM) (Cedeño, González-Morales & Ureña-López 2017; Ureña-López 2019; Linares Cedeño & Ureña-López 2021; Linares Cedeño, González-Morales & Ureña-López 2021; Medellín-González, Ureña-López & González-Morales 2021).

The goal of this paper is to offer the first cosmological simulation of FDM with attractive SI, and detailed and accurate treatment of the quantum pressure via a spectral method, in order to study the impact of instabilities on structure formation in the post-recombination universe. Local numerical simulations with attractive SI have been performed recently (Chen et al. 2021; Glennon & Prescod-Weinstein 2021) at the scale of one cluster in a static background. Cosmological simulations including gravity and attractive SIs in an expanding universe were carried out in Amin & Mocz (2019) using the Schrödinger–Poisson (SP) system. In that work, however, the focus was on soliton formation and their gravitational clustering rather than late-time structure formation.²

This paper is organized as follows. In Section 2, we lay out the non-relativistic limit for the axion dark matter model, relevant for our cosmological simulations. In Section 3, we describe the simulations. In Section 4, we discuss the impact of SI on the dark matter power spectrum. In Section 5, we explore the phase transition that affects dark matter solitons due to the SI. We offer our concluding remarks in Section 6.

2 FDM WITH ATTRACTIVE SELF-INTERACTION

We assume a real scalar field ϕ in the weak-field limit in an expanding universe, with an instantonic axion potential (Peccei & Quinn 1977; Di Vecchia & Veneziano 1980; Witten 1980) |$\mathcal {V}(\phi)$|⁠:

Such a system is governed by the Klein–Gordon–Einstein (KGE) equations.

In the non-relativistic limit (c → ∞), making the Klein transformation

to separate the fast oscillations (with a pulsation ω = mc²/ℏ ≫ H) from the slow evolution of the complex wavefunction ψ, the KGE equations reduce to the Gross–Pitaevski–Poisson (GPP) equations in an expanding universe:

where H is the Hubble constant, V is the gravitational potential seeded by the density ρ ≡ |ψ|², and a_s < 0 is an effective s-scattering length of the SI related to the axion decay constant f via:

For a detailed derivation of equations (4) and (5), see Chavanis (2018). In the above, the Hubble constant |$H\equiv \dot{a}/a$| encodes cosmological expansion, where a ≡ 1/(1 + z) is the cosmological expansion factor and z is the redshift.

Equations (4) and (5) with a_s = 0 are the SP equations (e.g. Schive et al. 2014a; Mocz et al. 2018) commonly used to simulate FDM neglecting SI. The |ψ|² and |ψ|⁴ terms in the equation come from a Taylor expansion of the non-relativistic limit of the instantonic axion potential equation (2). The |ψ|² is an attractive SI term. The next-order |ψ|⁴ term, only relevant at very high densities, is repulsive.

2.1 Soliton instability

The SP equations admit a well-known stable ground state soliton solution, approximated analytically by (Schive et al. 2014a):³

where r is the spherical coordinate, r_c is the core radius, and ρ₀ is the central density:

The soliton core has total mass M_c:

With attractive SI added, the soliton becomes unstable above a maximum critical mass (Chavanis 2011, 2018):

triggering a phase transition between dilute (equation 7) and dense solitons. The precise outcome of a dense soliton requires reverting back to the relativistic version of the governing physical equations; however, in the non-relativistic version of the equations, the repulsive |ψ|⁴ term may regularize and balance the attractive |ψ|² term and form a compact object of approximately constant density (Braaten et al. 2016; Chavanis 2018)

2.2 Linear instability scales

For non-relativistic self-gravitating Bose–Einstein Condensates with an attractive SI in an expanding universe, the equation for the density contrast in the linear regime is (Chavanis 2012)

where δ is the overdensity parameter. Equation (12) can be obtained from the hydrodynamic representation of the GPP equations. Structure formation results from the competition between the quantum pressure, the attractive SI, and the self-gravity. The competition between the quantum pressure and the self-gravity defines a (comoving) quantum Jeans wavenumber (Khlopov, Malomed & Zeldovich 1985):

The competition between the quantum pressure and the SI defines a (comoving) SI wavenumber (Chavanis 2011):

When all effects (gravity, quantum pressure, and SI) are taken into account, the critical wavenumber is obtained by putting the term in parenthesis in equation (12) equal to zero. This condition can be written as

Therefore, the critical wavenumber can be expressed in terms of k_I and k_J as

Jeans-type instability occurs for k < k_c.

In a cosmological context, with density ρ ∝ a⁻³ and a Hubble parameter h = 0.7, Desjacques et al. (2018) rewrites the comoving instability scales as

which allows us see their fiducial values and scaling. These scales imply that in the simple FDM model structure formation happens at physical scales which is larger than the Jeans scale k < k_J, and with SI, there is a secondary instability mode at k < k_I.

3 NUMERICAL SIMULATIONS

In this work, we consider five dark matter-only simulations: a reference CDM set-up, a 0-SI FDM run, and FDM runs with SI characterized by a scattering length a_s = −{1, 4, 8} × 10⁻⁷⁵ cm. Table 1 summarizes the simulation parameters and set-up, as well as calculates some relevant instability scales and masses described throughout the text. The axion mass is fixed to m = 10⁻²² eV. The FDM runs have a resolution of 1024³ and a box size of L_box = 1.5 h⁻¹Mpc. The box size is limited because we use a spectral method (Mocz et al. 2020) to evolve the wavefunction on a uniform grid and accurately resolve small-scale features and interference patterns that arise in solving Schrödinger-type systems. The numerical method is implemented as a module in the arepo code (Springel 2010), which is a state-of-the-art high-performance cosmological code for dark matter and baryonic simulations. The CDM simulation is performed using the N-body technique with a resolution of 512³. Our simulations use cosmological parameters of Ω_m = 0.3089, Ω_Λ = 0.6911, Ω_b = 0.0486, andh = 0.6774 consistent with the Planck observations of temperature and polarization anisotropies of the comic microwave background (Planck Collaboration XIII 2016).

Initial conditions are created as a random realization of a Gaussian field, with initial radial 1D power spectrum at redshift z = 127 calculated by axioncamb (Lewis, Challinor & Lasenby 2000; Hlozek et al. 2015). All simulations are generated with the same initial random seed for phases and amplitudes, allowing for direct comparison of structures across the simulations. In contrast to CDM, which is a scale-free theory where dark matter structure exists on all physical scales, in FDM there is a cut-off in the dark matter power above a wavenumber (Hu et al. 2000; Hui et al. 2017):

The simulations are run down to a redshift of z = 2, after which the uniform resolution is insufficient to resolve small-scale structures.

The SI strengths in our simulations are set to be stronger than the fiducial value that would predict the natural abundance of dark matter via equation (1), given our choice of axion particle mass m = 10⁻²² eV, i.e.f ≲ 10¹⁷ GeV. This choice was made for a few reasons: (1) Desjacques et al. (2018), estimate that large-scale structure is impacted at lower decay constants: f ∼ 10¹³ GeV. (2) Numerical limitations make stronger SI easier to resolve on cosmological scales, and we wish to numerically verify relevant instability scales. (3) Results may be interpolated between the SI and no SI cases. (4) The dark matter abundance (equation 1), does not necessarily have to hold for all axion-like particle models.

SI additionally affects the growth of perturbations in the early Universe, but can be safely neglected in the linear regime if the axion SI is f ≲ 3 × 10¹⁵ GeV (Desjacques et al. 2018; Chavanis 2021). We have neglected the effect of SI on structure formation in the linear regime: we have used identical initial conditions for all our FDM simulations, given by axioncamb which does not include SI. We point out that we have chosen strong SI strengths that would actually have some moderate effect on build-up of small-scale dark matter power in the linear regime prior to the epoch our simulations are started, which we have ignored. This approach makes it more straightforward to interpolate the simulations to weaker SI, whose effects are more difficult to resolve in our cosmological volume.

Since there is a cut-off of power in the initial power spectrum (equation 19), linear theory predicts that dark matter haloes form only down to a particular mass (Hui et al. 2017)

Non-linear structure formation may support less massive quantum ‘pressure’-supported haloes (solitons) of mass (Hui et al. 2017)

which has indeed been verified by numerical simulations (Mocz et al. 2019).

In our study, we investigate the effect of the SI instability scale (equation 14) on non-linear structure formation, as well as its impact on a soliton phase transition above M_max (equation 10).

4 DARK MATTER POWER SPECTRUM

Fig. 1 shows the evolving dark matter power spectrum in our five simulations at redshifts z = 127, 31, 15, 7, 3, 2. FDM (with and without SI), shows a reduction of power compared to CDM across all redshifts, due to the initial cut-off scale k_1/2. However, as seen in the power spectra, the inclusion of SI leads to the growth of additional small-scale power on the instability scale k_c (equation 16). Fig. 1 marks the location of the combined instability scale k_c, as well as the individual components: the Jeans instability scale k_J and the SI instability scale k_I. It can be seen that for our parameters, at high redshifts z ≳ 20, the instability occurs on the SI scale k_c ≃ k_I, while at lower redshifts z ≲ 20, the instability is set by the Jeans scale k_c ≃ k_J (Chavanis 2021). SI becomes less important at lower redshifts, as also seen in its scaling with a in equation (14).

Figure 1.

Evolving dark matter power spectra of our numerical simulations, with instability scales indicated. Shown also on the bottom row, are projected dark matter densities at z = 2, with blue arrows denoting formed solitons.

Open in new tabDownload slide

Interestingly, the dark matter power spectra for the various SI scenarios approximately converge by redshift z = 2 as the field evolves non-linearly. In terms of observational consequences and constraints, one could place on the model, SI has an effect of creating excess power at higher redshifts, z > 3 while resulting in similar power at lower redshifts. Most significant effects may be seen at z ∼ 7.

Fig. 1 also shows the projected dark matter density field in the bottom row. CDM is strikingly different since it forms dark matter sub-haloes on all spatial scales down to the numerical resolution limit. The FDM simulations (with and without SI) have reduced structure below k_1/2 and resemble each other more closely. However, the inclusion of SI has slightly accelerated structure formation, which has made filaments thicker and voids less dense at the <10 per cent level. For future work, it would be of interest to study with baryonic simulations how the change in filament potentials affects star formation. Insight from our previous study comparing baryons in CDM and FDM with no SI suggests that in these high-redshift first objects, baryons trace dark matter well, rather than having baryonic effects dwarf the signatures (Mocz et al. 2019). It would also be of interest to study stacked void profiles in larger-scale simulations to see how they differ between CDM/FDM/WDM.

5 SOLITONS

We have previously demonstrated in Mocz et al. (2019) via direct numerical simulation that the first structures that form in FDM are filamentary and undergo an instability to form solitons with mass as low as M_min (equation 21), which is below the cut-off scale predicted from linear theory, M_1/2 (equation 20). We observe a similar situation in our FDM simulation, where we form an M = 1.6 × 10⁸ M_⊙ soliton at redshift z = 2.2, which can be fit by the analytic soliton profile given by equation (7). Fig. 2 shows the measured radial profile and the analytic model, which provides a reasonable fit to the core size and central density. The soliton mass M is greater than M_min = 1.4 × 10⁷ M_⊙ but well below M_1/2 = 5 × 10¹⁰ M_⊙. The radial profiles are plotted in terms of physical units, rather than comoving units, as solitons are physical objects detached from cosmological expansion.

Figure 2.

Radial profile (physical units) of soliton at time of formation for the FDM simulations of various SI strengths. Note the presence of a dense soliton at the two highest SI strengths. For reference, a dilute soliton profile of radius 1.4 kpc is shown (thick grey line). The inset shows the redshift of formation of the solitons as a function of the SI strength.

Open in new tabDownload slide

It is interesting to observe the behaviour of solitons when SI is activated in the simulations. In the weakest SI case, a_s = −1 × 10⁷⁵ cm, solitons above M_max = 1.6 × 10⁸ M_⊙ are analytically expected to go unstable. The soliton is below this threshold and thus maintains its cored shape (Fig. 2) and is just slightly more compact and centrally concentrated due to the impact of the attractive SI. For simulations with stronger SI, the soliton mass M is now above the critical stable mass: M > M_max. Collapse is seen here, and the radial profiles are cuspy (Fig. 2). That is, the soliton has phase transitioned from a dilute to a dense state (Chavanis 2018). The simulation lacks the spatial resolution to fully resolve the final compact object with central density given by equation (11), which would be parsec-sized. The critical transition from dilute to dense solitons is also further corroborated with idealized simulations of a single quasi-stationary halo in Appendix  A.

Fig. 2 also demonstrates that SI also leads to earlier formation of solitons. The redshift z = 2.2 soliton in the no SI case forms before z = 5 in the strongest SI simulation. The formation of the soliton is defined as either the point in time that the filament forms an overdense ∼kpc core that can be approximated by the analytic soliton model, or forms a compact cusp (<kpc).

Finally, SI leads to the formation of additional solitons. The no SI FDM simulation forms just a single soliton in the filament by z = 2 in our 1.5 h⁻¹ Mpc box. However, as indicated in Fig. 1 by arrows, SI can cause constructive interference overdensities to collapse into dense solitons. This handful of additional dense solitons are difficult to resolve due to our limited spatial resolution.

6 CONCLUDING REMARKS

We have investigated ultra-light m = 10⁻²² eV FDM simulations with an attractive SI added, in order to explore the effect of the axion decay constant f on cosmic structure formation. We found that an axion decay constant of f ≲ 10¹⁴ GeV leads to a noticeable increase in small-scale power. This finding is consistent with analytic expectations for the instability scale due to attractive SI (Desjacques et al. 2018). SI also leads to the formation of dense rather than dilute solitons above a critical mass threshold; and thus, the prediction of the FDM model with SI is that the Universe would be populated with ‘bosenova’ that results from cosmological initial conditions (Levkov et al. 2017). Our simulations also show that increased SI leads to the formation of additional solitons in cosmic filaments, where interference patterns can cause over-densities that may be unstable under the SI. Given the above, our work highlights the important changes to the model predictions of FDM, if the boson is associated with an axion and the self-coupling is taken into account.

Our work has investigated a relatively low axion decay constant (f = 5.0 × 10¹³–1.4 × 10¹⁴ GeV), where the effects of SI are more noticeable. Such a low value would need a physical motivation beyond the simplest models. For a value of f ≃ 10¹⁷ GeV, which is the fiducial value that predicts the total dark matter abundance in the simplest models, the attractive SI would not have a significant impact on the structure of cosmic filaments. The critical mass for soliton collapse would also be significantly larger: M_max ≃ 10¹¹ M_⊙, which would not be cosmologically relevant to alter soliton core shapes, given the soliton core–halo mass relation (Schive et al. 2014b; Chavanis 2021). Hence, cosmological structure can place useful constraints on m and f simultaneously, which we leave for upcoming future work. Qualitatively, the inclusion of attractive SI goes in the right direction of solving the Catch-22 problem (Davies & Mocz 2020) that FDM currently faces: namely that a low particle mass m is needed to predict large, low-density cores, but that erases too much structure in the high-redshift Lyman α forest – which may be recovered to an extent with SI, without the need to invoke baryonic feedback physics.

ACKNOWLEDGEMENTS

Support (PM) for this work was provided by NASA through Einstein Postdoctoral Fellowship grant number PF7-180164 awarded by the Chandra X-ray Centre, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. AF was supported by a Royal Society University Research Fellowship. MV acknowledges support through NASA ATP 19-ATP19-0019, 19-ATP19-0020, 19-ATP19-0167, and NSF grants AST-1814053, AST-1814259, AST-1909831, AST-2007355, and AST-2107724. MBK acknowledges support from NSF CAREER award AST-1752913, NSF grants AST-1910346 and AST-2108962, NASA grant 80NSSC22K0827, and HST-AR-15809,HST-GO-15658, HST-GO-15901,HST-GO-15902, HST-AR-16159, and HST-GO-16226 from the Space Telescope Science Institute, which was operated by AURA, Inc., under NASA contract NAS5-26555. MA was supported by a NASA grant 80NSSC20K0518. SB was supported by the UK Research and Innovation (UKRI) Future Leaders Fellowship (grant number MR/V023381/1). TD acknowledges support from the Isaac Newton Studentship and the Science and Technology Facilities Council under grant number ST/V50659X/1. This work was performed using the Cambridge Service for Data Driven Discovery (CSD3), part of which was operated by the University of Cambridge Research Computing on behalf of the STFC DIRAC HPC Facility (www.dirac.ac.uk). The DIRAC component of CSD3 was funded by BEIS capital funding via STFC capital grants ST/P002307/1 and ST/R002452/1 and STFC operations grant ST/R00689X/1. DIRAC is part of the National e-Infrastructure. This material was based upon work supported by the National Science Foundation under grant DGE-2108962. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

References

APPENDIX: IDEALIZED SIMULATIONS OF SOLITON PHASE TRANSITION

We perform additional simulations of an idealized FDM halo with SI, to confirm the transition from dilute to dense solitons above the critical mass M_max (Chavanis 2018) in an idealized setting with higher effective resolution of the core. The set-up follows Mocz et al. (2017), where random solitons are merged to form a single quasi-stationary halo. The simulation has a box size of L = 20 kpc, resolution 400³, axion mass |$m=10^{-22} \, {\rm eV}$|⁠, and is run for 4 Gyr. In the reference case with SI switched off (a_s = 0), the result is a quasi-stationary halo with a soliton core of mass |$M=1.2\times 10^9 \, \mathrm{M}_\odot$| (radius 0.2 kpc). We consider additional simulation cases with SI strengths: |$a_s = -\lbrace 0.5,0.9,1.4,1.5,1.6,2 \rbrace \times 10^{-77} \, {\rm cm}$|⁠, corresponding to |$f=\lbrace 2, 1.5, 1.2, 1.15, 1.125, 1 \rbrace \times 10^{15} \, {\rm GeV}$|⁠.

Fig. A1 shows the resulting radial profiles of the halo for each SI strength. As the attractive SI strength increases, the soliton becomes more dense and compact. The phase transition is observed when the soliton mass is M > M_max, which is the case for the two strongest SI strengths simulated. The outer radial profile of the dark matter halo is close to an r⁻² isothermal profile, as analytically predicted in Chavanis (2019), and is largely unaffected by the collapse of the central soliton. A more detailed study of idealized collapse will be presented by Painter et al. (in preparation).

Figure A1.

Radial profiles and projected densities for idealized FDM halo with SI. A phase transition is observed to occur in the central soliton at large attractive SI strengths. For reference, a dilute soliton profile of radius 2 kpc is shown (thick grey line).

Open in new tabDownload slide