On the cosmic web elongation in fuzzy dark matter cosmologies: Effects on density profiles, shapes, and alignments of haloes

ABSTRACT

The fuzzy dark matter (FDM) scenario has received increased attention in recent years due to the small-scale challenges of the vanilla Lambda cold dark matter (ΛCDM) cosmological model and the lack of any experimental evidence for any candidate particle. In this study, we use cosmological N-body simulations to investigate high-redshift dark matter haloes and their responsiveness to an FDM-like power spectrum cutoff on small scales in the primordial density perturbations. We study halo density profiles, shapes, and alignments in FDM-like cosmologies (the latter two for the first time) by providing fits and quantifying departures from ΛCDM as a function of the particle mass m. Compared to ΛCDM, the concentrations of FDM-like haloes are lower, peaking at an m-dependent halo mass and thus breaking the approximate universality of density profiles in ΛCDM. The intermediate-to-major and minor-to-major shape parameter profiles are monotonically increasing with ellipsoidal radius in N-body simulations of ΛCDM. In FDM-like cosmologies, the monotonicity is broken, haloes are more elongated around the virial radius than their ΛCDM counterparts and less elongated closer to the centre. Finally, intrinsic alignment correlations, stemming from the deformation of initially spherically collapsing haloes in an ambient gravitational tidal field, become stronger with decreasing m. At z ∼ 4, we find a 6.4σ-significance in the fractional differences between the isotropized linear alignment magnitudes D_iso in the m = 10⁻²² eV model and ΛCDM. Such FDM-like imprints on the internal properties of virialized haloes are expected to be strikingly visible in the high-z Universe.

dark matter, large-scale structure of Universe, cosmology: theory

1 INTRODUCTION

1.1 Wave dark matter

The six-parameter standard Lambda cold dark matter (ΛCDM) cosmological scenario that emerged in the late 1990s can boast about a wide range of observational successes, from the Cosmic Microwave Background (CMB; Planck Collaboration XIII 2016), the Lyman-α forest (Iršič et al. 2017b), and galaxy clustering (Nuza et al. 2013) to weak gravitational lensing (Murata et al. 2018).

While this ‘minimal’ approach might be justified by the good fit to the CMB data and by model selection criteria like Occam’s razor, some of the theoretical arguments underpinning the ΛCDM approach lack sufficient justification (Di Valentino, Melchiorri & Silk 2015): The total neutrino mass ∑m_ν = 0.06 eV c⁻² is set to the minimal value allowed by solar and terrestrial neutrino oscillation experiments while cosmological data sets are sensitive to ∼100 meV variations and the dark energy equation of state is fixed to w = −1 so as to represent a cosmological constant rather than allow for time- and position-dependent components. Further, the effective number of light neutrino species is set to the Standard Model value of N_eff ∼ 3.046 which e.g. inflationary reheating or non-standard decoupling could alter.

CMB data assuming a ΛCDM cosmology gives rise to an anomalous weak lensing amplitude A_lens (Calabrese et al. 2008) and is in tension with low-redshift observations of Type Ia supernovae, standardizable quasars and cosmic chronometers (Ó Colgáin et al. 2022). Even worse, ΛCDM faces several small-scale crises such as the ‘missing satellite problem’ (Moore et al. 1999) and the ‘core-cusp problem’ (Wyse & Gilmore 2007): Galaxies predicted by CDM extend to much lower masses, well below the observed dwarf galaxies, with steeper, singular mass profiles. Dwarf spheroidal (dSph) galaxies in particular exhibit a surprising uniformity of their central masses, |$M(\lt 300 \ \text{pc}) \sim 10^7 \, \mathrm{M}_{\odot }$| (Strigari et al. 2008), and shallow density profiles (Bullock & Boylan-Kolchin 2017). Cf. Perivolaropoulos & Skara (2021) for a review on ΛCDM challenges.

In this work, we focus on fuzzy dark matter (FDM; Hui et al. 2017). The constituent particles in this dark matter (DM) model are drawn from an entire family (Arvanitaki et al. 2010; Visinelli & Vagnozzi 2019) of ultra-light bosons in the mass range 10⁻³³−10⁻¹⁹ eV. The lower bound, m_H ∼ 10⁻³³ eV, is the Hubble scale and is equivalent to saying that the Compton wavelength of FDM is smaller than the visible Universe, while the upper bound depends on what masses are to be considered ultra-light. The uncertainty principle that these quantum particles obey counters gravity below a particle mass-dependent Jeans scale or de Broglie wavelength λ_dB(m) (Khlopov, Malomed & Zeldovich 1985). The largely redshift-insensitive comoving de Broglie wavelength |$\lambda _{\text{dB},c} \sim (1+z)^{\frac{1}{4}}m^{-\frac{1}{2}}$| simultaneously suppresses small-scale structure and limits the central density of collapsed haloes (Schive, Chiueh & Broadhurst 2014a). The masses of cosmologically interesting FDM are traditionally expected to lie around m ∼ 10⁻²²–10⁻²⁰ eV: A dwarf galaxy in a halo of mass |$M_h \sim 4\times 10^9 \, \mathrm{M}_{\odot }\, h^{-1}$| and a virial radius |$R_{\text{vir}} \sim 40 \ \text{kpc}\, h^{-1}$| will have a virial velocity |$v _{\text{vir}} \sim 20$| km s⁻¹, which corresponds to a de Broglie wavelength of |$\lambda _{\text{dB}} \sim 7 \ \text{kpc}\, h^{-1}$| assuming m ∼ 5 × 10⁻²² eV.

Since the interparticle separation becomes smaller than the de Broglie wavelength for a DM particle mass of m < 30 eV, the ensuing wave character of FDM warrants the alternative term wave dark matter (Hui et al. 2017). The wave nature is reflected in the formation of interference patterns and solitonic cores in centres of haloes. Such cores match the observed flat cores of dwarf galaxies very well, and together with the suppression of small-scale structure solve some small-scale problems of ΛCDM.

However, experimental observations provide increasingly stringent bounds on the admissible mass ranges of FDM. Here, we disregard laboratory searches for axion-like signatures such as those involving haloescopes and helioscopes, some of which might have non-negligible sensitivity to ultra-light axions, i.e. FDM (Marsh 2016; Irastorza & Redondo 2018). Instead, we focus on particle mass constraints from astrophysical observations. In Fig. 1, we collate a non-exhaustive list of recent constraints ordered by publication date. While the density profiles in the central regions of dwarf spheroidals can be explained with an FDM soliton core potential provided that m ≲ 1.1 × 10⁻²² eV (Marsh & Pop 2015), Ly α observations find a lower bound m ≳ 3.8 × 10⁻²¹ eV to have enough power on Mpc-scales in the Ly α forest (Iršič et al. 2017a), constituting a Catch-22 problem. Recently, Safarzadeh & Spergel (2020) showed that a pure FDM cosmology is inconsistent with dwarf spheroidals across the entire parameter space by at least 3σ. More sophisticated models of reionization are needed though to make some of these observational constraints more reliable: For instance, most forest analyses assume that the frequency-averaged amplitude of the ionizing background J in the neutral hydrogen density n_H i ∝ (1 + δ_b)²/(T^0.7J) has negligible spatial fluctuations (Hui et al. 2017). Some FDM constraints might also be biased due to poor star formation modelling.

Figure 1.

Recent particle mass constraints for warm dark matter (WDM) / FDM from astrophysical observations in the FDM window of 10⁻²⁶ eV < m < 10⁻¹⁶ eV. The red shading indicates a disfavoured parameter space, though not necessarily a 2σ-constraint. The orange shading indicates forecasts for future observatories. Note that most constraints assume a pure WDM / FDM cosmology with a fixed DM particle mass, rather than DM cocktails or a distribution of DM particle masses. The resolution / fidelity of the simulations underlying some of the constraints differ widely. UFDs = ultrafaint dwarf galaxies; SMBHs = supermassive black holes.

High-redshift observations are most promising to provide evidence for/against FDM. The main reason can be traced back to non-linear structure formation, which redistributes power between scales so as to make the low-redshift FDM matter power spectrum resemble the ΛCDM one. For instance, in N-body simulations of the m ∼ 1.6 × 10⁻²² eV model, by z ∼ 1 discrepancies between the respective power spectra larger than 10 per cent exist only for scales k > 20 h c Mpc⁻¹ ∼2k_1/2 (Viel et al. 2012), corresponding to about twice the half-mode scale k_1/2. As we will see in Section 4, FDM imprints on the internal properties of virialized haloes persist for somewhat longer than suggested by the power spectrum alone.

1.2 Shapes of haloes and galaxies

Since the accretion of matter can be clumpy and directional, halo growth is anisotropic, resulting in the formation of non-spherical triaxial haloes, see e.g. Jeeson-Daniel et al. (2011). In ΛCDM theory, a blue-nugget transition for low-mass galaxies at z ∼ 2–4 was identified, a term introduced in Dekel & Burkert (2014). Before the transition, a DM-dominated central body causes the majority of galaxies to be aligned with their DM haloes. Both the galaxy and its host DM halo are on average fairly elongated, though as we will see in Section 4.3, their shapes strongly depend on halocentric distance. The transition itself is mediated by compaction events such as mergers, counter-rotating streams, or violent disc instabilities. Once complete, a self-gravitating baryonic core emerges, which scatters small-pericentre orbits of DM and stars, leading to a more spherical configuration. The outer DM halo remains less affected by baryonic self-gravitation in the centre.

Is this transition motivated by semi-analytical arguments confirmed by simulations? Chua et al. (2019) and Tomassetti et al. (2016) noted that the baryonic core emerges by virtue of radiative processes such as radiative heating and cooling, star formation, chemical evolution as well as strong supernova and AGN feedback. These allow a proper condensation of baryons into the centre of haloes, scattering DM particles that approach the halo centre and modifying box orbits into rounder passages. Eventually, a high star-formation rate (SFR) around cosmic noon (z ∼ 2) and a lower rate of gas inflow eventually leads to an inside-out quenching process into a compact, passive red nugget, the likely progenitor of the centre of an early-type galaxy today.

Given certain assumptions concerning the distribution of the 2D axis ratio Q for a sample of observed galaxies, one can use the inferred values of Q to put constraints on the distribution of the 3D axis ratios. For instance, Chang et al. (2013) and van der Wel et al. (2014) assumed that the triaxiality and the edge-on ellipticity are Gaussian distributed and concluded that low-mass galaxies (i.e. M_⋆ ≤ 7 × 10⁹ M|$_{\odot }\, h^{-1}$|⁠) are more elongated at z ≳ 1, while they are consistent with a population of flattened spheroids at the current epoch, hinting again at the blue-nugget transition. In a nutshell, the observational basis for this picture is becoming solid, both for the red-nugget phenomenon (Whitaker et al. 2012) and for their potential blue-nugget progenitors (Tacchella et al. 2015; Takeuchi et al. 2015). Since we are focusing on high redshifts of z ≳ 3.4, baryonic effects on the shapes of DM haloes are insignificant, further justifying our use of N-body simulations.

1.3 Correlated alignments of shapes / spin

While the high-z statistical properties of halo shapes can already hint at the nature of DM as we will see in Section 4.3, we gain additional insight when considering alignment correlations between haloes. The simplest approach is to consider geometrical alignment measures such as shape–position and shape–shape alignments, defined in Section 3.2. While these are easy to define, their measurement can be challenging, as exemplified in Pandya et al. (2019; CANDELS Collaboration). They assigned each of the observed galaxies probabilities of being intrinsically elongated, flattened, or spheroidal by essentially interpolating the results of Zhang et al. (2019). Yet they did not detect significant geometrical alignment signals in CANDELS observations, possibly owing to spectroscopic incompleteness.

More apt alignment measures can be motivated via the mechanisms underlying them: Correlated alignments are typically driven by stretching or compression of initially spherically collapsing mass distributions in some gravitational gradient. Alternatively, alignments can result from the mutual acquisition of angular momentum through tidal torquing of aspherical protogalactic mass distributions during galaxy formation, for instance due to the proximity of a large-scale cosmic filament. On large scales, the former effect is well captured by the Linear Alignment Model (LAM; Hirata & Seljak 2004) while the latter is described by tidal torquing theory (TTT; Doroshkevich 1970; White, McAdam & Jones 1984), a second-order effect. Starting in the 2000s, intrinsic alignment as found in the SDSS and WiggleZ DES data sets proved to be consistent with LAM and TTT predictions (Mandelbaum et al. 2006, 2011; Hirata et al. 2007; Lee 2011). However, the assumptions of TTT appear to be violated at lower redshift as hinted at by both simulations (Zjupa, Schäfer & Hahn 2021) and observations (KiDS+GAMA; Johnston et al. 2019). Given that we are post-processing N-body simulations, we will calculate LAM best fits in Section 4.4 while ignoring TTT.

Recently, the focus of alignment studies has shifted to isolating the contaminating contribution of large-scale intrinsic alignments to the weak gravitational lensing signal by the large-scale structure. In fact, intrinsic alignments are estimated to be one of the most serious physical systematic effects to the cosmic shear signal in the new era of systematics-dominated cosmology. The magnitude of the impact of intrinsic alignment on the observed lensing spectrum is comparable to the changes in cosmology even for a deep, stage IV weak lensing survey (Troxel & Ishak 2015). This next generation of weak lensing surveys comprising the Euclid mission, the Vera Rubin Observatory , and the Nancy Roman Space Telescope will all produce unprecedented weak lensing measurements in the coming decades and are all relying on proper alignment modelling. Understanding how FDM modifies intrinsic alignment correlations is an important ingredient in such models.

1.4 Outline

We focus on intermediate scales of k ∼ 2–18 h Mpc⁻¹ to study halo density profiles, shapes, and alignments at high redshifts of z > 3.4 using a suite of N-body simulations of CDM and FDM-like scenarios. The paper is laid out as follows: In Section 2, we describe our large-scale simulations. We present relevant definitions and the post-processing tools in Section 3. Results for the high-z statistics of FDM-like haloes are given in Section 4. A summary follows in Section 5. The Supplementary Materials are dedicated to initial conditions and pre-initial conditions in simulations with a power spectrum cutoff, to convergence tests and to the impact of quantum pressure on profiles of density and shape.

2 CDM AND cFDM SIMULATIONS

We search for imprints of a small-scale cutoff in the initial DM power spectrum using cosmological N-body simulations performed with the state-of-the-art moving-mesh code arepo (Springel 2010). Gravitational forces are computed using a Tree-PM method where long-range forces are calculated on a particle mesh and the short-range forces are calculated using a tree-like hierarchical multipole expansion scheme.

Bona fide FDM simulations are much more challenging than CDM simulations as the oscillations with highest frequency |$\omega \propto m^{-1}\lambda _{\text{dB}}^{-2}$| occur in the densest regions, requiring very fine temporal resolution even for moderate spatial resolution. For instance, a de-Broglie wavelength of λ_dB = 0.6 kpc h⁻¹ and velocity of |$v =200$| km s⁻¹ translate to an oscillation time-scale of τ_osc = 2π/ω ∼ 3 × 10⁶ yr h⁻¹. In this work, we use a widely adopted approximation for FDM which we shall call classical FDM (cFDM). As opposed to a bona fide superfluid, cFDM approximates FDM as a classical collision-less fluid, governed by the Vlasov-Poisson equation, but with FDM initial conditions. The initial small-scale suppression in the power spectrum is modelled using axioncamb (Hložek et al. 2015). We will refer to this exponential-like suppression as a cutoff. By construction, cFDM ignores dynamical effects of the quantum pressure that are apparent on small scales (Mocz et al. 2019, 2020; May & Springel 2021). Note that cFDM particles do not possess thermal velocities, i.e. they are not a valid WDM implementation.

Why is cFDM a valid approximation? The scales of interest in this work are the intermediate ones of k ∼ 2–18 h Mpc⁻¹, corresponding to halo virial masses M_h = 4π(π/k)³ρ_m/3 of |$M_h \sim 2\times 10^9\!-\!10^{12} \ \mathrm{M}_{\odot }\, h^{-1}$|⁠. On these scales, quantum pressure is not very significant. The reason is that the absolute fractional difference between growth rates in FDM versus cFDM does not exceed 10 per cent for particle masses¹ around m ∼ 10⁻²² eV and mass scales around |$M \sim 2\times 10^9 \ \mathrm{M}_{\odot }\, h^{-1}$| (Schive et al. 2016; Corasaniti et al. 2017). In fact, this fractional difference is less than 5 per cent for mass scales beyond |$M \sim 4\times 10^9 \ \mathrm{M}_{\odot }\, h^{-1}$|⁠, representing the vast majority of the haloes in our inventory, described in Section 3.3. Note that these estimates are based on the solutions to the linearized governing equations (Schrödinger-Poisson versus Vlasov-Poisson) extrapolated down to z_end ∼ 4 in an EdS-like universe without baryons, as in this work.

The N-body suite offers a much cleaner platform to assess the imprints of a cutoff in the initial power spectrum, as in full hydrodynamical runs subgrid baryonic physics uncertainties would make it challenging to disentangle the resolution effects that are due to baryonic physics from e.g. the resolution convergence of the iterative shape procedure outlined in Section 3.1. For instance, Chua et al. (2019) finds poor shape convergence in illustris simulations that can be traced back to the underprediction of galaxy stellar mass and galaxy formation efficiency at lower resolutions (Vogelsberger et al. 2013).

We use cosmological volumes with two different box side lengths, L_box = 10 and 40 cMpc h⁻¹, for each of three DM resolutions, N = 256³, 512³, and 1024³. This set of simulation parameters best balances the competing demands of high resolution (for convergence in e.g. halo shapes) and large volume (to obtain accurate statistical distributions). It also allows us to test convergence both in statistical distributions and with respect to mass resolution. Our N-body simulations have initial conditions generated at z = 127 using ngenic(Springel et al. 2005), an initial conditions code which relies on the Zeldovich approximation. The cosmological boxes are evolved until z= 3.4. We run simulations over a range of bosonic particle masses, m = 10⁻²² eV, 7 × 10⁻²² eV, 2 × 10⁻²¹ eV. The fixed comoving Plummer-equivalent gravitational softening length is set to ϵ = 0.19 ckpc h⁻¹. We adopt a Planck-like cosmology with Ω_m = 0.3089, Ω_Λ = 0.6911, h = 0.6774, σ₈ = 0.81 (and n_s = 0.9665 for the N-body simulations of ΛCDM), cf. Planck Collaboration XIII (2016). See Table 1 for an overview of the N-body simulations.

Table 1.

Overview of the N-body simulations and the parameters used: (1) simulation type, including the particle mass m in case of cFDM; (2) side length of simulation box; (3) number of DM resolution elements; (4) mass per DM resolution element. While the 10 cMpc h⁻¹ boxes are used for convergence tests, the bulk of this work analyses the 40 cMpc h⁻¹ boxes.

Simulation type	Box side length (cMpc h−1)	DM particles	mDM (⁠|$10^6 \, \mathrm{M}_{\odot }$|⁠)
ΛCDM	10	2563 / 5123 / 10243	5.11 / 0.64 / 0.08
ΛCDM	40	2563 / 5123 / 10243	327 / 40.9 / 5.11
cFDM, m = 2 × 10−21 eV	10	2563 / 5123 / 10243	5.11 / 0.64 / 0.08
cFDM, m = 2 × 10−21 eV	40	2563 / 5123 / 10243	327 / 40.9 / 5.11
cFDM, m = 7 × 10−22 eV	10	2563 / 5123 / 10243	5.11 / 0.64 / 0.08
cFDM, m = 7 × 10−22 eV	40	2563 / 5123 / 10243	327 / 40.9 / 5.11
cFDM, m = 10−22 eV	10	2563 / 5123 / 10243	5.11 / 0.64 / 0.08
cFDM, m = 10−22 eV	40	2563 / 5123 / 10243	327 / 40.9 / 5.11

Substructure is identified using the subfind algorithm explained in Springel et al. (2001). After identifying DM haloes (called groups) using the friends-of-friends (fof) algorithm with a standard linking length of b = 0.2 × mean interparticle separation, the algorithm searches for overdense regions inside each fof group and prunes them according to a gravitational boundedness criterion. In the case of cFDM, subfind identifies many particles beyond the virial radius R_vir as part of the halo if they are gravitionally bound to it, providing an opportunity to study the ambient cosmic environment of haloes in Section 4.3.

3 POST-PROCESSING TOOLS

3.1 Shape determination algorithm

Ideally, one would like to find the shapes of isodensity surfaces for a given 3D particle distribution. With only |$\lesssim \mathcal {O}(10^4)$| particles, however, the next-best choice is to assume that the distribution can be reasonably fit by an ellipsoid, see Zemp et al. (2011). The axis ratios can be computed from the unweighted shape tensor q_ij, which is the mode-centred second moment of the mass distribution:

$$\begin{eqnarray*} q_{ij} = \frac{1}{\sum _k m_k} \sum _k m_k r^{\text{mode}}_{k,i}r^{\text{mode}}_{k,j}. \end{eqnarray*}$$

(1)

Here, m_k is the mass of the k-th particle, and |$r^{\text{mode}}_{k,i}$| is the i-th component of its position vector with respect to the densest location within the cloud, i.e. mode. We have compared the potential centre definition with the result of the ‘shrinking sphere’ algorithm that we use to find the densest location following Power et al. (2003), and have found very good agreement. The overall shape could be described by the axis ratios

$$\begin{eqnarray*} q := \frac{b}{a}, \ \ s := \frac{c}{a}, \end{eqnarray*}$$

(2)

where a, b, and c are the eigenvalues of q_ij corresponding to the major, intermediate, and minor axes, respectively. The ratio of the minor-to-major axis s has traditionally been used as a canonical measure of the sphericity. A perfectly spherical cloud has q = s = 1. The triaxiality parameter, first defined by Franx et al. (1991) as

$$\begin{eqnarray*} T := \frac{1-q^2}{1-s^2}, \end{eqnarray*}$$

(3)

measures the prolateness / oblateness of a halo. T = 1 describes a completely prolate (colloquially elongated) halo, while T = 0 describes a completely oblate (colloquially flattened) halo. In practice, haloes with T > 0.67 are considered prolate while those with T < 0.33 are dubbed oblate. Haloes with 0.33 < T < 0.67 are said to be triaxial.

Note that q and s are in fact a function of distance from the centre of the particle cloud. A proper distance measure to choose is the ellipsoidal radius

$$\begin{eqnarray*} r_{\text{ell}} = \sqrt{x_{\text{pf}}^2+\frac{y_{\text{pf}}^2}{(b/a)^2}+\frac{z_{\text{pf}}^2}{(c/a)^2}}, \end{eqnarray*}$$

(4)

measured from the mode. Here, (x_pf, y_pf, z_pf) are the coordinates of the particle in the eigenvector coordinate system of the ellipsoid (also called principal frame), i.e. r_ell corresponds to the semimajor axis a of the ellipsoidal surface through that particle.

We use an iterative process (Dubinski & Carlberg 1991; Katz 1991) to compute the axis ratios q(r_ell) and s(r_ell). First, we discretize the major axis into logarithmically spaced radial bins. At a given semimajor axis r_ell ≡ a, we start by considering particles inside a sphere of radius a. Next, the shape tensor q_ij(r_ell), equation (1), is diagonalized which yields the eigenvectors and eigenvalues at distance r_ell. The eigenvectors give the directions of the semiprincipal axes. As opposed to the differential version of this algorithm, the geometrical meaning of the eigenvalues remains somewhat uncertain and dependent on the mass density profile (Zemp et al. 2011). It is generally assumed though that the eigenvalues of q_ij are still proportional to the semimajor axes squared. With this approximation, we calculate the axis ratios

$$\begin{eqnarray*} q = \sqrt{\frac{q_{yy}}{q_{xx}}}, \ \ s = \sqrt{\frac{q_{zz}}{q_{xx}}}, \end{eqnarray*}$$

(5)

where q_zz < q_yy < q_xx are the principal components of the tensor. Now, we keep the length of the semimajor axis fixed – though the orientation can change – and calculate q_ij again by summing over all particles within the new, deformed integration volume (= ellipsoid) with semimajor axis a = r_ell and axis ratios q and s. In other words, particle (x_pf, y_pf, z_pf) enters the summation in equation (1) if

$$\begin{eqnarray*} \sqrt{x_{\text{pf}}^2+\frac{y_{\text{pf}}^2}{q^2}+\frac{z_{\text{pf}}^2}{s^2}} \le a. \end{eqnarray*}$$

(6)

This iteration is repeated until convergence is reached. As a convergence criterion, we typically require that the fractional difference between two iteration steps in both axis ratios is smaller than 10⁻².

Even if the particle resolution exceeds |$\mathcal {O}(10^4)$| particles, we refrain from running the differential version of this algorithm which replaces ellipsoids with ellipsoidal shells of width Δa (Zemp et al. 2011). While some of our haloes of the intermediate- and high-mass end would satisfy this criterion for the 1024³ resolution runs, many shells at large ellipsoidal radii struggle to converge with the differential version.

It is instructive to estimate the convergence radius r_conv beyond which we deem our shape calculations reliable. Power et al. (2003) demonstrated that deviations from convergence depend (for appropriate choices of other numerical parameters) solely on the number of particles, and scale roughly with the collisional ‘relaxation’ time, t_relax. The latter quantifies after how much time the actual velocity of test particles significantly differs from the velocity that they would have had if the mass of the other particles were smoothly distributed, rather than finitely sampled. Since t_relax(r) scales roughly like the enclosed number of particles times the local orbital time-scale, one often expresses t_relax in units of the circular orbit time-scale t_circ at R₂₀₀. According to |$t_{\text{circ}}(R_{200})= 2\pi R_{200}/v _{R_{200}}=\sqrt{3\pi /(G\bar{\rho })} = \sqrt{3\pi /(200G\rho _{\text{crit}}\Omega _m)}$|⁠, the latter is of the order of the age of the Universe. The rescaled relaxation time may then be written as

$$\begin{eqnarray*} \kappa (r) = \frac{t_{\text{relax}}(r)}{t_{\text{circ}}(R_{200})} = \frac{\sqrt{200}}{8} \frac{N(r)}{\ln N(r)}\left[\frac{\bar{\rho }(r)}{\rho _{\text{crit}}}\right]^{-1/2}, \end{eqnarray*}$$

(7)

where N(r) is the enclosed number of particles and |$\bar{\rho }(r)$| is the mean enclosed density within r. While Power et al. (2003) found that to achieve deviations of no more than 10 per cent in the local density profile requires κ ∼ 1, Vera-Ciro et al. (2011) found that κ = 7 marked the convergence radius where the axis ratios of individual relaxed haloes agreed in different resolution runs of a quiescent halo in the Aquarius simulations of Milky Way-mass dark matter haloes. In other words, as our inner convergence radius r_{ρ, conv} for density profiles we thus adopt the halocentric radius r_{ρ, conv} = 〈r〉, typically averaged over a mass bin, that satisfies κ(r_{ρ, conv}) = 1, and analogously for shape profiles using the ellipsoidal radius.

Predicting halo shapes around the virial radius R_vir is of importance due to the success of the spherical and ellipsoidal collapse models. However, there is also considerable interest in studying halo shapes at radii that correspond to the outskirts of the spiral galaxies which they may host, including our own Milky Way. We will thus put some focus on the specific radius of R₁₅ = 0.15 R_vir. The choice of R₁₅ is motivated by observational measurements of the properties of the Milky Way halo, which using tidally disrupted dwarf satellites have been performed out to 40 kpc h⁻¹ from the Galactic Centre (Law & Majewski 2010).

3.2 Alignment correlations

Following Pandya et al. 2019, we consider two distinct types of alignment: shape-position alignments and shape-shape alignments. As for the former, we are interested in how the major axis of a particle distribution relates to the vector pointing at the mode of another particle distribution in the sample. As for the latter measure, we quantify how much the major axes of two particle distributions deviate from one another.

We compute the alignment angles trivially using

$$\begin{eqnarray*} |\cos \theta | = \frac{|\mathbf {v}_1\cdot \mathbf {v}_2|}{\Vert \mathbf {v}_1 \Vert \ \Vert \mathbf {v}_2 \Vert } , \end{eqnarray*}$$

(8)

where θ is the alignment angle between the two vectors |$\mathbf {v}_1$| and |$\mathbf {v}_2$|⁠, and the absolute value is taken since shapes in the axisymmetric ellipsoidal approximation are spin-2 objects. For shape-position alignments, |$\mathbf {v}_1$| is the major axis of one particle group while |$\mathbf {v}_2$| is the position difference vector connecting the two modes. For shape–shape alignments, both |$\mathbf {v}_1$| and |$\mathbf {v}_2$| are major axes.

Central subhaloes can be identified using subfind: The subhalo in each fof group with the lowest potential resolution element is classified as central. They are expected to exert a continuous gravitational torque on non-centrals, colloquially called satellites, aligning them partially with the halo. This is exemplified by the halo alignment model by Schneider & Bridle (2010). Since we are interested in alignments mediated by the large-scale structure, we discard satellites in our alignment studies. In fact, we discard satellites throughout this entire work, including for density and shape profile estimations, all of which they would bias otherwise.

3.3 Halo inventory

As a minimum halo resolution we reject all haloes consisting of fewer than 400 DM resolution elements, i.e. |$M_{\text{min}} = 400 \times m_{\mathrm{DM}} = 2.0 \times 10^9 \, \mathrm{M}_{\odot }\, h^{-1}$| for the runs with 1024³ resolution, L_box = 40 cMpc h⁻¹. This resolution floor is sufficient to discard the vast majority of small-mass clumps that stem from the artificial fragmentation of filaments (Lovell et al. 2014), even for the m = 10⁻²² eV halo sample. We conclude this by comparing M_min to the location of the artificial small-mass upturn in the inferred halo mass function (not shown). At redshift z = 4.38, this leaves us with ∼46 000 haloes for CDM, ∼39 000 haloes for cFDM with m = 2 × 10⁻²¹ eV, ∼32 000 haloes for cFDM with m = 7 × 10⁻²² eV and ∼5000 haloes for cFDM with m = 10⁻²² eV. The ratios of these numbers are in good agreement with the numbers obtained when extrapolating the Schneider et al. (2012a) WDM-to-CDM recalibration fit

$$\begin{eqnarray*} \frac{\mathrm{d}n^{\text{WDM}}/\mathrm{d}M}{\mathrm{d}n^{\text{CDM}}/\mathrm{d}M}(M) = \left(1+\frac{M_{1/2}}{M}\right)^{-\beta } \end{eqnarray*}$$

(9)

with β = 1.16 and half-mode mass M_1/2 (see equation 11) to higher redshifts. Here, we use the half-mode matching formula

$$\begin{eqnarray*} m_X = 0.84 \left(\frac{m}{10^{-22}\text{eV}}\right)^{0.39} \ \text{keV} \end{eqnarray*}$$

(10)

to translate from a cFDM particle mass m to a thermal relic WDM mass m_X (Marsh 2016). When considering mass density profiles, we impose additional relaxedness conditions outlined in Section 4.2 which reduce the number of haloes by a factor of ∼0.2. Finally, for the construction of median shape profiles we additionally discard haloes whose shape determination fails at the virial radius R_vir, leading to a reduction in the number of haloes by a factor of ∼0.25.

4 DENSITY, SHAPE, AND ALIGNMENT STATISTICS

To assess the high-z statistics of dark matter haloes in the cFDM cosmologies, we investigate four statistical measures of the large-scale structure in order of complexity.

4.1 Inferred power spectra

The simplest such measure is the two-point statistics of the DM density field, which is traced by DM resolution elements in the arepo simulations. The theoretical axioncamb power spectra at z = 127 drop by many orders of magnitudes at high wavenumbers in an approximately exponential fashion, accompanied by various axion wiggles. The resolvability of the cutoff with a finite number of DM resolution elements is an indicator of the quality of the initial conditions. We find that the inferred power spectra from the simulations replicate the cutoff down to very low power spectra magnitudes of P(k) ∼ 10⁻¹³ h⁻³ Mpc³. For a cFDM cosmology with m = 10⁻²² eV, the latter value is attained at a spatial scale of k ∼ 30 h Mpc⁻¹.

The ratio of cFDM to CDM power spectra for different axion masses are shown in Fig. 2 for a wide range of redshifts.² To obtain such high-k power spectrum estimation fidelity, we specified a non-canonical Piecewise Cubic Spline (PCS) mass assignment scheme for the kernel of the window function |$W(\mathbf {x})$| to obtain the density field on a mesh of resolution 1024³. PCS is a third-order window function, i.e. there are p = 3 + 1 grid points per dimension to which each particle is assigned (Sefusatti et al. 2016). As a comparison, the widely adopted cloud-in-cell (CIC, cf. Hockney & Eastwood 1981) mass assignment scheme is a first-order window function. We compensate for the window function |$W(\mathbf {x})$| by trivially deconvolving the kernel in Fourier space.

Figure 2.

Evolved power spectrum ratios P(k)_cFDM/P(k)_CDM of the comoving DM density field, colour-coded by redshift. Results are shown for N-body runs with 1024³ resolution and L_box = 40 cMpc h⁻¹. The vertical lines mark the FDM half-mode scale k_1/2 (orange) and the Nyquist frequency k_Ny (magenta). The DM model is indicated on top of each panel and the redshift is specified by the colourbar labelling. Ratios of linear power spectra as obtained with axioncamb at the corresponding redshifts are traced by the dashed curves.

The colour-coding in Fig. 2 reflects the logarithmically spaced scale factors a that were chosen for the snapshots. Constant spacing in the logarithm of the scale factor has desirable properties when considering galaxy populations and cosmology, see Baldry (2018). Since the differences between two scale factors thus translate to approximately (exactly in Einstein-de-Sitter space) constant comoving distances, the power spectra amplitudes undergo approximately constant multiplicative shifts as time evolves, at least in the regime of linear structure formation. We show curves for 34.3 ≥ z ≥ 3.4.

The half-mode scale (vertical orange line in Fig. 2) is illustrative to consider since the cutoff occurs over a range of scales, and following Marsh (2016) the definition we adhere to is

$$\begin{eqnarray*} T(k_{1/2}) = \left(\frac{P_{\text{lin}}^{\text{FDM}}(k_{1/2},z)}{P_{\text{lin}}^{\text{CDM}}(k_{1/2},z)}\right)^{\frac{1}{2}} = \frac{1}{2}, \ \forall z \ge 0, \end{eqnarray*}$$

(11)

where T(k) is the redshift-independent FDM-to-CDM linear transfer function ratio.³

The linear theory predictions in Fig. 2 are obtained with axioncamb and thus bona fide linear FDM theory. The fact that the linear theory power spectrum ratios of FDM to CDM do not all overlap is indicative of scale-dependent growth, a feature of FDM even when its governing equations (Schrödinger-Poisson) are linearized (Marsh 2016). For CDM and cFDM, linear growth always corresponds to scale-independent growth due to the properties of the Vlasov-Poisson system of equations. The non-linear predictions as obtained with the N-body code show that at high redshifts of z ≳ 25 and large scales both CDM and cFDM follow the standard linear evolution of the power spectrum P(k, a) ∝ a²P(k) with growing scale factor a. This is indeed predicted for the matter-dominated era, with distinct modes evolving independently as a distribution of Gaussian random fields. It is the smaller scales for which this simple approximation breaks down first. Instead, non-linear structure formation successfully redistributes power to smaller scales.

For the cFDM cosmologies, we observe a wiggly pattern of power at high k which goes beyond the mere axion wiggles and is thus indicative of numerical lattice features. As the grid pre-initial condition (cf. Supplementary Materials) that we employ simply places N³ particles onto the grid points of a 3D Cartesian lattice, the regular grid spacing introduces periodic signals into the pre-initial condition. The periodic signals get reinforced by gravity at high redshift while actual structure formation eventually overtakes the N-body discreteness effects in amplitude. The lattice features can persist in low-density regions of very low-resolution simulations even at redshift z = 0, cf. Wang & White (2007) and Liao (2018). It is evident from Fig. 2 that the redshift z_* at which lattice features become insignificant on all resolvable scales is dependent on the particle mass m, or equivalently the half-mode scale k_1/2. More precisely, z_* grows with increasing m. Amongst the 1024³ resolution, L_box = 40 cMpc h⁻¹ runs, lattice features are most profound for the m = 10⁻²² eV scenario for which they become unrecognizable on all resolvable scales around z_* ∼ 12.

We admit that in our convergence tests in which we also employ simulations of resolution 256³ and 512³, some unavoidable lattice-induced small-scale power is still present, especially for the 256³ runs with m = 10⁻²² eV, L_box = 10 cMpc h⁻¹, and n.b. L_box = 40 cMpc h⁻¹. However, halo density and shape distributions and by extension shape alignments are robust against very low-magnitude small-scale power. The question of ideal initial loads for simulations with power spectra cutoffs is an intricate one. We try to answer this question in the Supplementary Materials, in favour of grid pre-initial conditions.

4.2 Density profiles

Before we continue in Section 4.3 with shape statistics results, we first remind the reader of some well-known results for FDM-like halo density profiles and then generalize them to the new parameter space.

It is well established that quasi-equilibrium CDM haloes have an approximately universal density profile that can be reproduced by rescaling the NFW profile (Navarro, Frenk & White 1996)

$$\begin{eqnarray*} \frac{\rho (r)}{\rho _{\text{crit}}} = \frac{\delta _c}{(r/r_s)(1+r/r_s)^2}. \end{eqnarray*}$$

(12)

The latter is fully determined by two parameters, the characteristic or scale radius r_s ≡ r₋₂ at which the logarithmic slope has the isothermal value of −2, i.e. |$\mathrm{d} \ln {\rho } / \mathrm{d} \ln r|_{r_{-2}} = -2$|⁠, and a characteristic overdensity δ_c. This universality of spherically averaged density profiles holds over many orders of magnitudes in mass, cosmological parameters, magnitude and shape of the initial density fluctuation power spectrum and even in several modified gravity models (Wang & White 2009; Hellwing et al. 2013). However, halo density profiles are better approximated with the Einasto profile (Einasto 1965)

$$\begin{eqnarray*} \ln \left(\frac{\rho _E(r)}{\rho _{-2}}\right) = -\frac{2}{\alpha }\left[\left(\frac{r}{r_{-2}}\right)^{\alpha }-1\right]. \end{eqnarray*}$$

(13)

While it is often reported that the third parameter is needed to more accurately capture the halo-to-halo variation in profile shape as well as its mass-dependence, this conclusion is not correct. For even if the third parameter α is fixed as a function of halo mass and redshift to around α ∼ 0.17, the Einasto profile provides a better fit to simulation results (Gao et al. 2008).

It has been found (González-Samaniego, Avila-Reese & Colín 2016) that low-redshift WDM halo profile shapes can be well accounted for by the NFW and Einasto profiles. We thus use the Einasto profile as our fitting model for cFDM cosmologies and assess the fidelity thereof a posteriori. We obtain ρ(r) curves by brute-force binning of DM resolution elements into 20 mode-centred radial bins, equally spaced in log r spanning −2.0 ≤ log (r/R₂₀₀) ≤ 0.0. We also compute the total enclosed mass M(r) and mean inner density profiles |$\bar{\rho }(r) = M(r)/(4/3)\pi r^3$| in order to compute κ-profiles according to equation (7).

We construct the c(M, z) relation by largely following Ludlow et al. (2016): First, we remove substructure contamination by focusing on particles identified by subfind as part of central subhaloes. We then fit the median mode-centred density profiles after averaging over logarithmic mass bins⁴ of width Δlog M = 0.2. The advantage of this approach is that it smooths out any features unique to individual systems and dampens the influence of outliers. Best-fitting Einasto profiles are determined by adjusting the three parameters of equation (13) in order to minimize a figure-of-merit defined as

$$\begin{eqnarray*} \psi ^2 = \frac{1}{N_{\text{bin}}} \sum _{i = 1}^{N_{\text{bin}}} \left[\ln \rho _i - \ln \rho _E(r_i; \rho _{-2}, r_{-2}, \alpha)\right]^2. \end{eqnarray*}$$

(14)

Here, the number of radial bins is N_bin < 20 since we discard those that fall below the inner convergence radius estimate r_{ρ, conv} given by κ(r_{ρ, conv}) = 1. This choice for r_{ρ, conv} differs from the choice in Ludlow et al. (2016): They choose r_{ρ, conv} = 3ϵ, where ϵ is the gravitational softening length. While both estimates for r_{ρ, conv} were first introduced in (Power et al. 2003), the latter is less robust since it constitutes a halo mass-independent and thus halo resolution-independent estimate. In addition, radial bins that exceed the outer limit of 0.8 R₂₀₀ are also discarded in the fit since they may correspond to radii that are not fully relaxed (Ludlow et al. 2016).

To account for smooth accretion, minor mergers, and occasional major mergers with systems of comparable mass that drive individual haloes out of quasi-equilibrium, we choose to discard unrelaxed haloes by combining two conditions. The first is to request haloes to have a normalized offset between the mode |$\mathbf {r}^{\text{mode}}$| of the halo and its centre of mass |$\mathbf {r}^{\text{CoM}}$| of no more than |$d_{\text{off}} = |\mathbf {r}^{\text{mode}}-\mathbf {r}^{\text{CoM}}|/R_{200}\lt 0.07$|⁠. As mentioned in Section 2, some cFDM haloes are deeply embedded into nearby cosmic filaments, especially at high redshift. Consequently, we choose to evaluate d_off by excluding DM resolution elements outside of R₂₀₀ to avoid filament-induced biases. The second condition we impose is that the virial parameter η = (2K − S_p)/|W| < 1.35. Here, K is the kinetic energy of the halo, W its potential energy, and S_p the surface pressure term

$$\begin{eqnarray*} S_p = \oint \rho (\mathbf {x}\cdot \mathbf {v})(\mathbf {v}\cdot \mathrm{d}\mathbf {S}) = 4\pi R_{\text{vir}}^3\rho _{\text{vir}}v _r^2 \end{eqnarray*}$$

(15)

that follows from relaxing the assumption that the density at the halo boundary is zero (Davis, D’Aloisio & Natarajan 2011; Klypin et al. 2016). The surface integral is taken over the surface of a sphere of virial radius R_vir, and the density ρ_vir and radial velocity dispersion |$v _r$| at the virial radius can be estimated from the halo density profile. The surface pressure correction is of considerable importance in the high redshift, rapid merging phase of the Universe, and upon its exclusion one would categorize many relaxed haloes as unrelaxed (Klypin et al. 2016).

We avoid using alternative relaxedness conditions such as the one based on the substructure mass fraction f_sub = M_sub(< R₂₀₀)/M₂₀₀ (Neto et al. 2007). It is often constrained to f_sub < 0.1, yet the inferred f_sub is resolution-dependent and thus heavily underestimated for low-mass systems. Ludlow et al. (2016) resorts to a dynamical age requirement to flag haloes undergoing rapid accretion and equal-mass merging, a largely redundant condition in view of our imposed η < 1.35. In other words, our set of two conditions proves sufficient to remove the ‘upturn’ in the concentration of high-mass haloes reported by e.g. Klypin et al. (2011).

The results of the Einasto fitting procedure for four mass bins in the N-body, 1024³ resolution, L_box = 40 cMpc h⁻¹ runs at redshift z = 4.38 are shown in Fig. 3. The curves show the median density profile in the respective mass bin while the shaded area delineates the 25–75th percentile, measuring the halo-to-halo variation. We find that both CDM and cFDM haloes for the m = 10⁻²² eV scenario have density profiles that can be well fit with an Einasto profile. The red and blue dashed verticals demarcate the convergence radius r_{ρ, conv} defined by κ(r_{ρ, conv}) = 1, as explained around equation (7). As is evident from said equation, larger mass haloes have smaller convergence radii. In the plotted range of radii, their density profiles also do not exhibit spuriously large densities close to the halo centres, as opposed to small-mass haloes (first row).

Figure 3.

Density profiles for N-body, 1024³ resolution, L_box = 40 cMpc h⁻¹ runs at redshift z = 4.38, for cFDM with particle mass m = 10⁻²² eV (red), versus CDM haloes (blue). The mass bins span |$10^9-10^{10}\, \mathrm{M}_{\odot }\, h^{-1}$| (top left-hand panel) all the way to |$10^{12}-10^{13}\, \mathrm{M}_{\odot }\, h^{-1}$| (bottom right-hand panel). See legends for details. The shaded areas delineate the 25–75th percentiles. The red and blue dashed verticals denote κ = 1 lower convergence radii, naturally moving towards smaller halocentric radii for higher mass haloes. The κ = 1 lines for CDM and cFDM overlap in most panels. The r⁻¹ and r⁻² power laws trace two characteristic regimes of an NFW-like profile. The cyan dashed verticals show the (Schive et al. 2014b) estimates for the soliton core radius that one would find in bona fide FDM simulations. For the two higher mass bins, the core radius is too small to be shown, r_c/R₂₀₀ = 4.5 × 10⁻³ for |$10^{11}-10^{12}\, \mathrm{M}_{\odot }\, h^{-1}$| and r_c/R₂₀₀ = 1.1 × 10⁻³ for |$10^{12}-10^{13}\, \mathrm{M}_{\odot }\, h^{-1}$|⁠.

The core radii of solitons that one would encounter in bona fide FDM simulations are given by cyan dashed verticals following (Schive et al. 2014b). The well-known trend of core radii growing with decreasing halo mass is recovered. In the mass bin |$10^{11}-10^{12}\, \mathrm{M}_{\odot }\, h^{-1}$|⁠, the soliton core radius r_c/R₂₀₀ ∼ 4 × 10⁻³ normalized by the median R₂₀₀ value is not shown and similarly for the mass bin |$10^{12}-10^{13}\, \mathrm{M}_{\odot }\, h^{-1}$| for which r_c/R₂₀₀ ∼ 10⁻³. Even in the lowest mass bin of |$10^{9}-10^{10}\, \mathrm{M}_{\odot }\, h^{-1}$|⁠, the core radius is smaller than the convergence radius, demonstrating that we would not expect solitons to contaminate our density profile inferences at this level of resolution. In the Supplementary Materials, we show in detail how quantum pressure in FDM induces alterations of the halo density profile.

Finally, note that in CDM the median scale radius r₋₂ migrates towards larger radii as the halo mass grows. This is of great importance. In fact, the major advantage of NFW and Einasto profiles is that their scale radius r_s, often called concentration radius, can be used to define the halo concentration

$$\begin{eqnarray*} c := \frac{R_{200}}{r_{-2}} \end{eqnarray*}$$

(16)

as the ratio of R₂₀₀ to that of the scale radius. This definition can be extended straightforwardly to alternative DM scenarios such as WDM and FDM. As discussed by NFW for the case of CDM, c and M₂₀₀ do not take on arbitrary values, but correlate in a way that reflects the mass-dependence of halo formation times (Bullock et al. 2001). Earlier assembly corresponds to higher characteristic densities, reflecting the larger background density at that epoch.

Ludlow et al. (2014) later showed that the concentration of a CDM halo can be inferred from the critical density of the Universe at a characteristic time along its mass accretion history. The insight was used to construct an analytic model for the mass-concentration-redshift relation c(M, z) for CDM haloes, reproducing the inferred relation from many CDM simulations. In Ludlow et al. (2016), the generalization of the model has been introduced and applied to WDM, relying on the full ‘collapsed mass history’ of a halo instead of its main progenitor. As such, the model is based on extended Press-Schechter theory, allowing us to assess the dependence of halo concentrations on cosmological parameters and on the shape of the linear matter power spectrum alike. It is this model whose validity for cFDM and high redshifts we would like to comment on.

Fig. 4 shows parts of the c(M, z) relation as inferred from the N-body, 1024³ resolution, L_box = 40 cMpc h⁻¹ runs at redshift z = 4.38. Square markers indicate concentrations as obtained from Einasto best fits in the respective mass bin, while the colour of each hexbin is proportional to the number of haloes whose concentration falls therein. The solid blue curve labelled ‘LB 16, CDM’ is a broken power-law fitting formula for c(M, z) in the Planck cosmology from Ludlow et al. (2016) (equation C1). At fixed redshift z, the concentrations of CDM haloes decrease monotonically with increasing halo mass, reflecting the lower background density at the epoch of their formation. At fixed mass, concentrations c decrease monotonically with increasing redshift (not shown), indicative of the redshift evolution of the reference density, also called the pseudo-evolution of mass (Diemer, More & Kravtsov 2013).

Figure 4.

Concentration–mass relation c(M, z) in different cosmologies, for N-body, 1024³ resolution, L_box = 40 cMpc h⁻¹ runs at redshift z = 4.38. We compare the CDM haloes (blue) with cFDM haloes for particle mass m = 2 × 10⁻²¹ eV (top), m = 7 × 10⁻²² eV (middle), and m = 10⁻²² eV (bottom). The square markers indicate median concentrations obtained from Einasto best fits while the colour of each hexbin is proportional to the number of haloes whose Einasto concentration falls therein. ‘LB 16, CDM’ and ‘LB 16, cFDM’ show an analytical model prediction from Ludlow et al. (2016) for model parameters f = 0.02 and C = 650 (see the text).

For cFDM, the concentration-mass relation is non-monotonic. In solid red labelled ‘LB 16, cFDM’, we add the model prediction from Ludlow et al. (2016). To maintain consistency, as input to the model we use linear FDM matter power spectra obtained with axioncamb rather than the Viel et al. (2005) parametrization suitable for WDM. The model parameter f that is used to define the halo collapse redshift is set to f = 0.02 while the proportionality factor C relating mean inner halo densities to the critical density of the Universe at the collapse redshift is set to C = 650. We find that their model slightly overpredicts the cFDM concentrations at the small-mass end, while performing better for higher masses. At given z, concentrations peak at around two orders of magnitude above the half-mode mass scale M_1/2, with concentrations declining above and below this peak mass scale M_peak ∼ 100 × M_1/2. This is a known result which has also been described by e.g. Bose et al. (2015). However, both the latter work and that of Ludlow et al. (2016) focus on WDM simulations at low redshifts of z ≤ 3 and equivalent cFDM particle masses of 3.3 × 10⁻²¹ and 4.5 × 10⁻²² eV. We extend the result to cFDM of lower particle mass of m = 10⁻²² eV and to redshifts in the range 3.4 ≤ z ≤ 6.3, beyond which our halo catalogues become too sparse and noisy to make statistically robust conclusions. We concede that for the m = 1 × 10⁻²² eV run (lower panel), our lack of high-mass haloes makes the extrapolation towards M_peak only tentative.

As in WDM, the non-monotonicity of c(M, z) is due to the suppression of gravitational collapse below the half-mode scale, breaking the scale-invariance of the assembly process. It imprints a preferred scale on the mass accretion histories. While Ludlow et al. (2016) apply their model of the c(M, z) relation to WDM, here we show that it also fares well for cFDM simulations with initial conditions based on axioncamb, reproducing the non-monotonicity and even the constancy of M_peak with cosmic time.

Another viewpoint on halo concentrations can be obtained by looking at the inferred probability distribution function (PDF) for c(M, z). To this end, Fig. 5 juxtaposes normalized histograms of CDM concentrations versus cFDM ones, again at redshift z = 4.38. Following Jing & Suto (2002), we fit the PDFs with a phenomenological lognormal distribution

$$\begin{eqnarray*} p(c)\mathrm{d}c = \frac{1}{\sqrt{2\pi }\sigma _{c}}\exp \left[-\frac{(\ln c - \ln \bar{c})^2}{2\sigma _c^2}\right] \mathrm{d}\ln c, \end{eqnarray*}$$

(17)

which we find to perform well across all mass bins, cFDM particle masses, and again for redshifts in the range 3.4 ≤ z ≤ 6.3. The scatter in c(M, z) around the median relations has been investigated by other authors. While Diemer & Kravtsov (2015) reports a CDM scatter of 0.16 dex at all redshifts and masses using a halo inventory that includes unrelaxed haloes, Neto et al. (2007) finds |$\sigma _{\log _{10}c}\sim 0.1$| in a selection of mass bins for relaxed haloes at z = 0. At z = 4.38, we identify a CDM scatter of σ_c ∼ 1.32 in the highest mass bin, which translates to |$\sigma _{\log _{10}c}\sim 0.12$| and is thus consistent. The monotonicity of concentrations in CDM is evident from the PDFs, as in the analogous, bottom panel of Fig. 4.

Figure 5.

Concentration probability distribution functions in different cosmologies, for N-body, 1024³ resolution, L_box = 40 cMpc h⁻¹ runs at redshift z = 4.38. We compare CDM haloes (blue) for mass bin |$10^9-10^{10}\, \mathrm{M}_{\odot }\, h^{-1}$| (top), |$10^{10}-10^{11}\, \mathrm{M}_{\odot }\, h^{-1}$| (middle), and |$10^{11}-10^{12}\, \mathrm{M}_{\odot }\, h^{-1}$| (bottom) with cFDM haloes for particle mass m = 10⁻²² eV. The solid curves are lognormal best fits.

We take a moment to note that the agreement we find with Ludlow et al. (2016) is conditional on adopting (most of) their selection functions, their choices for calculating halo density profiles and their Einasto fitting procedure. Most importantly, by discarding unrelaxed haloes the median c(M, z) relation in CDM becomes strictly monotonically decreasing with halo mass. If we were to include unrelaxed haloes, we would encounter an ‘upturn’ in the concentration of high-mass haloes. Rather than ‘LB 16, CDM’, a much better fitting formula would then be provided by Diemer & Joyce (2019). While their analytic framework naturally incorporates unrelaxed haloes in CDM, its extension to cFDM would necessitate a new parameter in the model. Since the upturn at the high-mass end is largely insensitive to a primordial power spectrum truncation, we remove the upturn and accept a bias towards dynamically older systems.

Another controversy concerns the very definition of halo concentration for a three-parameter model such as the Einasto profile. For the latter, the parameter α co-determines not only the shape of the profile but also the density of the halo at its centre. An intuitive definition of halo concentration c mapping more dense centres to higher concentrations must make c a function of both R₂₀₀/r₋₂ and α. For instance, the ratio of maximum to virial circular velocities can be cast into a concentration by assuming a profile parametrization as is done in Klypin et al. (2016), but it can also be used as an alternative method for characterizing the halo concentration that does not assume any density profile (Gao et al. 2004). Other refinements could be added by considering ellipsoidal profiles rather than spherical profiles (cf. Section 7), in which case we would take into account the shape of the dark matter particle distribution, constraining the bins to follow the main directions of this distribution. When using ellipsoidal profiles, fitted masses, and concentrations thus tend to be higher (Gonzalez et al. 2022).

Since our goal is to generalize and extend density profile results found for WDM to cFDM, we choose the simple c = R₂₀₀/r₋₂ definition, contend with spherically averaged profiles, and exclude unrelaxed haloes, all in agreement with Ludlow et al. (2016).

4.3 Shape statistics

The lack of a fully self-consistent theory of DM halo structure formation let alone galaxy formation motivates the use of cosmological simulations. This is even more true for alternative DM scenarios, for which many seminumerical methods such as a generalized halo model have only recently been introduced (Schneider et al. 2012a).

It is known that halo growth is anisotropic by virtue of accretion being clumpy and directional, such as along filaments and sheets. While the anisotropy is known to result in non-spherical haloes (Macciò, Dutton & van den Bosch 2008), statistical properties of cFDM halo shapes have not yet been investigated, to the best of our knowledge. Here, we fill this gap by focusing on how halo shapes are affected by gradually decreasing the cFDM particle mass m.

Shape profiles out to 5R_vir are calculated by applying the mode-centred Katz-Dubinski algorithm of Section 3.1 on the central subhaloes. Fig. 6 compares CDM and cFDM median shape profiles for N-body, 1024³ resolution, L_box = 40 cMpc h⁻¹ runs at redshift z = 4.38. We highlight beyond-virial radius regions in the plot as they ought to be interpreted as part of the cosmic environment rather than the haloes per se. Let us begin by outlining some shape characteristics of CDM haloes. We find that haloes are least spherical near the halo centre, with axis ratios 〈q〉 ∼ 0.65 and 〈s〉 ∼ 0.45 for intermediate-mass haloes at R₁₅ = 0.15 R_vir. As a comparison, the same group of haloes is significantly more spherical near the virial radius R_vir with 〈q〉 ∼ 0.70 and 〈s〉 ∼ 0.50.

Figure 6.

Shape profiles in different cosmologies, for N-body, 1024³ resolution, L_box = 40 cMpc h⁻¹ runs at redshift z = 4.38. We show the intermediate-to-major axis ratio q (first row), the minor-to-major axis ratio s (second row), and the triaxiality T (third row) profiles for CDM (blue) as compared to cFDM (red) with particle mass m = 10⁻²² eV for halo mass bins |$10^{9}-10^{10} \, \mathrm{M}_{\odot }\, h^{-1}$| (first column), |$10^{10}-10^{11} \, \mathrm{M}_{\odot }\, h^{-1}$| (second column), and |$10^{11}-10^{12} \, \mathrm{M}_{\odot }\, h^{-1}$| (third column). Median difference curves are drawn in solid purple. Halo-to-halo variations are represented by shaded regions which enclose the 25th–75th percentiles. The dashed blue and red lines delineate κ = 7 lower convergence radii. The shaded light-grey region denotes the ambient cosmic environment of the haloes.

We speculate that the monotonicity of q and s as a function of ellipsoidal radius is a property of hierarchical structure formation, i.e. occurs in cosmological models in which the dimensionless matter power spectrum Δ²(k) = k³P(k)/(2π²) increases with wavenumber k. In theory, it should thus be predictable using (semi-)analytical arguments. Similar to halo density profiles, hierarchical clustering clearly does not produce objects with a universal shape, but the distributions may be close to universal.⁵ Spherically averaged halo density profiles ρ(r) are considered universal if they are purely a function of total halo mass M. Current efforts to understand the c-M relation are precisely tackling the question of universality since the latter is broken unless concentration can itself be described as a function of mass. Klypin et al. (2016) find that in the so-called plateau regime, concentration does not depend on halo mass. A potential universality for shapes is underinvestigated as of yet, most likely due to its higher level of complexity arising from 3D considerations as opposed to 1D. We advocate for more emphasis thereon as it is beyond the scope of this work.

As with density profiles, special care is needed to understand how non-relaxed haloes modify any such monotonicity. For our shape analysis, including Fig. 6, we do not impose any relaxedness condition, which might be the reason why we find the monotonicity in q and s to break for the highest mass haloes, which are expected to undergo many mergers and are still in the fast-accretion regime in the nomenclature of Diemer & Joyce (2019). On a similar note, how does monotonicity in q and s change in the low-redshift Universe of z ≲ 4? As N-body simulations including our own ones indicate, the monotonicity becomes a steeper one as time evolves.⁶

Another consequence of the higher rate of mergers of high-mass CDM haloes that we just hinted at is the enhanced prolateness of such haloes, already identified by Allgood et al. (2006). At z = 4.38, low-mass haloes with |$M = 10^{9}-10^{10} \, \mathrm{M}_{\odot }\, h^{-1}$| have a triaxiality of 〈T〉 ∼ 0.70 at the virial radius while high-mass haloes are more elongated with 〈T〉 ∼ 0.80.

4.3.1 Impact of cosmic filaments on cFDM haloes

How does a power spectrum cutoff affect halo shapes? As we see in Fig. 6 for the rather extreme m = 10⁻²² eV model, cFDM haloes are less oblate than their CDM counterparts beyond and around the virial radius R_vir. The lower values of q and s are indicative of the increased role cosmic filaments play in cFDM, into which most haloes at z = 4.38 are embedded. As we will show in future work, cosmic filaments in cFDM feature higher mean overdensities than in CDM and in relative terms contribute more to the build-up of large-scale tidal forces. This leads to an elongation (higher T-values) of DM haloes at large ellipsoidal radii. Their profound gravitational influence renders small-mass haloes less spherical (reduced q and s) even down to deeper regions around R₁₅. High-mass cFDM haloes of mass |$M = 10^{11}-10^{12} \, \mathrm{M}_{\odot }\, h^{-1}$| are less disrupted by cosmic filaments as their deeper gravitational potential wells provide greater shielding.

To appreciate the redshift- and cFDM particle mass-dependence of shape profiles, Fig. 7 shows triaxiality curves over a wider range in parameter space. We find that the influence of cosmic filaments on shape profiles fades away as time evolves. Beyond-virial radius profiles of cFDM haloes for m = 2 × 10⁻²¹ eV closely match their CDM counterparts around z = 5.56, while the profiles do not show much resemblance even at z = 4.38 for the more extreme m = 10⁻²² eV run. This suggests that with power redistributing among the various scales due to non-linear structure formation, cFDM / FDM filaments eventually fragment over time into less prolate structures (Smith & Markovic 2011; Mocz et al. 2019). This complicated process takes place over a redshift range that is hard to pinpoint, as the fragmentation time of a filament depends on its core radius r₀, amongst others. The latter is a good proxy for where the simple model of isothermal cylinders in hydrostatic equilibrium breaks down (Ramsøy et al. 2021). Some filaments prove stable enough so as to not break up at all. Overall though, the effect of the power spectrum cutoff becomes less noticeable at lower redshifts. Since structure formation is predominantly linear and thus more pristine at higher redshift, we conclude that with increasing redshift it is easier to distinguish cFDM from CDM. This is insofar as the theoretically relevant but not directly observable halo shapes are concerned. An accurate high-z determination of galaxy shapes and a thorough understanding of their dependence on parent halo shapes thus might shed light on the nature of DM.

Figure 7.

Triaxiality profiles in different cosmologies, for N-body, 1024³ resolution, L_box = 40 cMpc h⁻¹ runs. We compare CDM (blue) to cFDM with particle mass m = 2 × 10⁻²¹ eV (magenta), m = 7 × 10⁻²² eV (dark-red), and m = 10⁻²² eV (orange-red) across redshifts indicated above the columns. The shaded regions enclose the 25th–75th percentiles. The results are for intermediate-mass haloes in the range |$10^{10}-10^{11} \, \mathrm{M}_{\odot }\, h^{-1}$|⁠. The vertical dashed lines delineate κ = 7 lower convergence radii. The cyan-coloured dashed horizontal line at T = 2/3 denotes the triaxial-prolate transition.

Enhanced prolateness of haloes is expected to be observed in WDM simulations as well, since thermal velocities are incapable of modifying halo shapes (in addition to the linear effects stemming from the matter power spectrum suppression) beyond a very small core with size of a few parsecs for m_X ∼ 1 keV (Macciò et al. 2012). Hence, we find that a combination of anisotropic gravity and a primordial power spectrum cutoff results in an elongation of the cosmic web at high-z, which is most noticeable for haloes and filaments. Crudely speaking, the cutoff-mediated elongation of filaments is passed on to the embedded haloes via tidal shearing, until the tenuous filaments eventually break up around z ∼ 6 for m ∼ 7 × 10⁻²² eV and elongation fades away with non-linear structure formation, cf. Mocz et al. (2019).

4.3.2 Anticorrelation between sphericity and halo mass

We confirm the positive correlation between triaxiality T and halo mass M at the virial radius R_vir. It is more customary though to quantify the correlation of the median sphericity 〈s〉 with halo mass M. Similar to the results of Chua et al. (2019), Butsky et al. (2016), we find a negative correlation of 〈s〉 with M for CDM haloes. Following Allgood et al. (2006), we express the sphericity–mass relation as

$$\begin{eqnarray*} \langle s \rangle = a\left(\frac{M_{\text{vir}}}{M_{\ast }(z)}\right)^b, \end{eqnarray*}$$

(18)

where M_vir represents the virial mass and M_*(z) is the cosmology-dependent characteristic non-linear mass for redshift z such that the RMS real-space top-hat smoothed overdensity at scale R(M_*) = (3M_*/(4πρ_critΩ_m))^1/3 is δ_c = 1.68. In other words,

$$\begin{eqnarray*} \sigma ^2(M_{\ast }(z), z) = \int _0^{\infty }\frac{{\rm d}k}{k}\left[\frac{3j_1(kR(M_{\ast }))}{kR(M_{\ast })}\right]^2\Delta ^2(k,z) \stackrel{!}{=} \delta _c^2, \end{eqnarray*}$$

(19)

where j₁(x) is a spherical Bessel function of the first kind and Δ²(k, z) = k³P^lin(k, z)/(2π²). Allgood et al. (2006) found that the dependence of the distributions of halo shapes on the amplitude of density perturbations, σ₈, was well described by the cosmology dependence of M_* alone. We remind the reader that the power spectrum P^lin(k, z) that appears in the definition of the variance σ² refers to the linear power spectrum in the corresponding cosmology, extrapolated to possibly low redshift z. The same remark holds true for the definition of |$\sigma _{8} := \left[\sigma ^2(8\ \text{Mpc}\, h^{-1}, z=0)\right]^{\frac{1}{2}}$|⁠. In ΛCDM, |$M_{\ast }(z=0) = 4.66\times 10^{12}\, \mathrm{M}_{\odot }\, h^{-1}$|⁠, |$M_{\ast }(z=1) = 1.54 \times 10^{11}\, \mathrm{M}_{\odot }\, h^{-1}$|⁠, and |$M_{\ast }(z=3.41) = 8.84 \times 10^{7}\, \mathrm{M}_{\odot }\, h^{-1}$|⁠.

For cFDM cosmologies, one could resort to the Viel et al. (2005) parametrization of the linear WDM-to-CDM transfer function ratio, since scales beyond k = 10 h Mpc⁻¹ are largely suppressed in the integrand of equation (19) by the Bessel factor, while the linear cFDM power spectrum below k = 10 h Mpc⁻¹ can be well approximated by the Viel et al. (2005) parametrization. However, for pure cFDM and WDM cosmologies |$\mathcal {C}$|⁠, |$M_{\ast }(z, \mathcal {C})$| ceases to exist for larger redshifts, e.g. z ≳ 2 for m = 10⁻²² eV or equivalently m_WDM = 0.84 keV. The reason is that σ²(M, z) flattens off for small M, instead of diverging as lim_{M → 0}σ²(M, z) → ∞ as happens for cosmologies with a high-k power index n ≥ −3 such as ΛCDM.

In addition, |$M_{\ast }(z, \mathcal {C} = \Lambda \text{CDM})$| drops quickly below the resolution limit |$m_{\text{DM}} = 5.11 \times 10^6 \, \mathrm{M}_{\odot }\, h^{-1}$| of the N = 1024³, L_box = 40 cMpc h⁻¹ runs at around z ≳ 4. Table 2 thus shows best-fitting parameters for z = 3.41 only, with |$M_{\ast }(z, \mathcal {C} = \Lambda \text{CDM})$| values adopted for the cFDM cosmologies as well.

Table 2.

Fitting parameters to the mass-dependence equation 〈p〉 = a(M_vir/M_*(z = 3.41))^b for two different radii R_vir and R₁₅, with p = q or s. Results are shown for N-body, 1024³ resolution, L_box = 40 cMpc h⁻¹ runs.

	CDM		m = 2 × 10−21 eV		m = 7 × 10−22 eV		m = 10−22 eV
	a	b	a	b	a	b	a	b
qvir	0.759	−0.012	0.768	−0.009	0.740	−0.004	0.437	0.069
q15	0.648	0.006	0.652	0.005	0.631	0.012	0.494	0.044
svir	0.541	−0.011	0.521	−0.006	0.511	−0.005	0.201	0.124
s15	0.289	0.081	0.276	0.080	0.267	0.089	0.233	0.098

We confirm the anticorrelation (b < 0) between sphericity and halo mass that has been found previously by other authors (Allgood et al. 2006; Chua et al. 2019) and extend the result to moderate values of the cFDM particle mass m. What is the physical origin of this phenomenon? Again, we expect haloes with masses above M_* to undergo a higher rate of mergers than haloes with masses below M_*. Since it has been shown that this merging happens along preferred directions (Zentner et al. 2005), the enhanced prolateness and reduced sphericity is primarily due to merging. Halo merging thus plays a significant role in the distribution of shapes.

For particle mass m = 10⁻²² eV, however, we observe that the anticorrelation between sphericity and halo mass breaks down and turns into a positive correlation. This is true for both b = 0.069 at q_vir and b = 0.124 at s_vir. We find no clear anticorrelation for CDM in the central regions around R₁₅ as opposed to some low-redshift results, e.g. Allgood et al. (2006), Chua et al. (2019). However, the magnitude of the power index b is still preferentially increased by the power spectrum cutoff, up to b = 0.098 at s₁₅ for cFDM with m = 10⁻²² eV. Some authors have reported no significant change in the magnitude of the power index b across redshift, with some even claiming a steepening with redshift (Jing & Suto 2002). However, Fig. 2 in Allgood et al. (2006) hints at the opposite trend, which we can confirm. In fact, the magnitude of b is so low at z ≳ 4 that the fitting exponent is consistent with zero in CDM, another reason against pursuing such a fitting procedure at too high redshifts.

4.3.3 PDF of halo shapes

For the PDFs of the axis ratios, we recover the CDM trends already identified in Schneider, Frenk & Cole (2012b) at lower redshifts: As Fig. 8 indicates for z = 4.38, sphericity s follows an approximately symmetrical distribution, with cFDM haloes featuring lower values especially beyond and around R_vir. Lowest mass haloes in cFDM, however, peak at sphericities as low as s ∼ 0.15, indicating that the majority of these haloes is highly non-spherical. We note that T is negatively skewed, suggesting that there is a heavy tail of oblate and triaxial haloes in both CDM and cFDM scenarios. As in Fig. 7, we find a reduced prolateness (smaller median of T) of intermediate- and especially high-mass cFDM haloes compared to the CDM ones around R₁₅. We see that for intermediate-mass cFDM haloes, this can be traced to an increased median in the q-distribution and a reduced median in the s-distribution. For high-mass cFDM haloes, it primarily stems from an increased median in the q-distribution.

Figure 8.

Probability density functions of the shape parameters in different cosmologies, for N-body, 1024³ resolution, L_box = 40 cMpc h⁻¹ runs, at redshift z = 4.38. We compare the intermediate-to-major axis ratio q (top row), the minor-to-major axis ratio s (middle row) and the triaxiality T (bottom row) PDFs in cFDM with particle mass m = 10⁻²² eV (red) to CDM (blue) at R₂₀₀ (solid lines) and R₁₅ (dashed lines). Mass bins are |$10^{9}-10^{10} \, \mathrm{M}_{\odot }\, h^{-1}$| (left column), |$10^{10}-10^{11} \, \mathrm{M}_{\odot }\, h^{-1}$| (middle column) and |$10^{11}-10^{12} \, \mathrm{M}_{\odot }\, h^{-1}$| (right column). The cyan-coloured dashed vertical line at T = 2/3 characterizes the triaxial-prolate transition. While the symbols denote normalized histogram values, the curves plot Gaussian kernel density estimates thereof.

The nature of DM thus has strong imprints on halo shapes, most pronounced for the smallest explored haloes at R₁₅. Specifically, such haloes in cFDM are less spherical (smaller s and q) while having an overall higher triaxiality (larger T). It is possible that given a primordial power spectrum and expansion history of the Universe, one could estimate the shape distribution of virialized objects at a given redshift z (semi-)analytically and thus theorize, for instance, the smaller median of T for intermediate- and high-mass haloes in cFDM cosmologies. To the best of our knowledge, such an analysis has yet to be performed.

4.4 Intrinsic alignment

The reason why intrinsic alignments have gained considerable attention in recent years is twofold. First, the next generation of galaxy weak lensing surveys such as the Euclid mission will suffer from highly biased cosmological inferences without a proper alignment modelling informed by simulations and seminumerical methods. At the same time, intrinsic alignments of galaxies provide complementary cosmological information on their own, assuming they can be disentangled from the lensing signal by means of either galaxy colour (Yao et al. 2020) or polarization data (Brown & Battye 2011).

To remain consistent throughout this work, we focus on the linear intrinsic alignment strengths of DM haloes and the impact of a power spectrum cutoff thereon. One should be careful when carrying over our conclusions for the linear alignment strengths of haloes to the population of elliptical galaxies they may host. Ellipticals have more spherical gas distributions than the underlying DM halo, which is a corollary of the X-ray shape theorem (Buote & Canizares 1994), since the gravitational potential represents an overall average of the local density profile. Gas pressure is also more isotropic than the anisotropic velocity ellipsoids of a DM halo, further contributing to more spherical gas distributions. These theoretical arguments are corroborated by various observations (Buote & Canizares 1998; Morandi, Pedersen & Limousin 2010). Moreover, the misalignment between ellipticals and their DM haloes can pose additional problems in case such a misalignment is correlated with the local tidal fields rather than being purely random. To make matters worse, the linear intrinsic alignment strength of spiral galaxies is expected to be about an order of magnitude smaller than those for ellipticals (Zjupa et al. 2021).

4.4.1 Geometric alignments

The simplest way to quantify the intrinsic alignment of haloes is via geometry. One measure is shape-position alignment, which we have introduced in Section 3.2. It can be considered a special case of the more general notion of ellipticity-direction cross-correlation (Lee et al. 2008). Likewise, the shape–shape alignment measure can be mapped onto the ellipticity–ellipticity correlation function.

In Fig. 9 we present shape–position alignment statistics for CDM as compared to cFDM haloes at redshift z = 4.38, restricted to mass bin |$10^{10}-10^{11} \, \mathrm{M}_{\odot }\, h^{-1}$|⁠. We find that for large pair separations, the median value for |cos θ| asymptotically approaches 0.5, a value consistent with a uniform random shape orientation and halo clustering. For small pair separations, intermediate-mass CDM haloes conspire to a median |cos θ| value of around |$0.653^{+0.006}_{-0.002}$|⁠, which is in very close agreement with the value Pandya et al. (2019) found for real-space 3D mock light cones that they generated for the CANDELS survey at lower redshifts of z = 1.0–2.5. What we can also confirm (not shown) is that the shape–shape alignment strength barely reaches values of |cos θ| ∼ 0.6 for smallest separation halo pairs, which is systematically lower than the shape–position alignment magnitude.

Figure 9.

High-z shape–position alignment statistics in different cosmologies, for N-body, 1024³ resolution, L_box = 40 cMpc h⁻¹ runs, at redshift z = 4.38. We compare CDM haloes (blue) to cFDM haloes with particle mass m = 2 × 10⁻²¹ eV (magenta), m = 7 × 10⁻²² eV (dark-red), and m = 10⁻²² eV (orange-red). The results are for intermediate-mass haloes in the range |$10^{10}-10^{11} \, \mathrm{M}_{\odot }\, h^{-1}$|⁠. Top: Shape–position alignments versus 3D pair separation, out to L_box/2 = 20 cMpc h⁻¹ to avoid geometric bias induced by the simulation box. The shaded regions enclose the 25th–75th percentiles. Bottom: PDF of shape–position alignments, with results of the Kolmogorov–Smirnov test of independence between the CDM and cFDM samples.

For cFDM cosmologies, Fig. 9 clearly demonstrates how the geometric alignment is gradually enhanced as the cFDM particle mass m is decreased. For our most extreme cFDM model with m = 10⁻²² eV, the median |cos θ| curve climbs up to |$0.753^{+0.021}_{-0.010}$|⁠, the wider-spread percentiles effected by the smaller number of haloes compared to CDM. The enhancement is due to the lack of small-scale structure such as smaller-mass haloes that would regularize the tidal fields generated by quasi-linear cosmic filaments and to a lesser extent 2D cosmic sheets. Shape–position alignments are significantly different from the uniform random result even out to 3D pair separations of 10 cMpc h⁻¹, where |cos θ| ∼ 0.54. Upon inspection of the |cos θ| distribution after marginalizing over 3D pair separations, we find that the Kolmogorov–Smirnov p-value, p_KS, for rejecting the hypothesis that the cFDM and CDM signals are drawn from the same distribution is very low. In other words, the hypothesis can be ruled out with very high significance. As expected, p_KS increases with m.

4.4.2 Intrinsic alignment modelling

While the linear alignment model (LAM) introduced in Hirata & Seljak (2004) and Catelan et al. (2001) is traditionally applied to elliptical galaxies, the condition of a virialized, velocity dispersion-stabilized system is also satisfied by haloes. To keep our focus on the properties of the cosmic web and to harness greater statistical power, we investigate intrinsic alignment correlations of haloes within the LAM. In this model, the halo is assumed to be spherically symmetric in isolation, perturbed by the presence of a cosmic tidal field. Let σ²(r) denote the isotropic velocity dispersion of DM resolution elements in the unperturbed halo with |$r=|\mathbf {r}| = |\mathbf {x}-\mathbf {x}_0|$| being the distance of the particle to the centre of mass of the halo. In a steady-state Jeans equilibrium one has

$$\begin{eqnarray*} \frac{1}{\rho } \frac{\mathrm{d}(\sigma ^2\rho)}{\mathrm{d}r} = - \frac{\mathrm{d}\Phi _{\text{eq}}}{\mathrm{d}r}, \end{eqnarray*}$$

(20)

with density ρ(r) and gravitational potential Φ_eq(r). The Jeans equilibrium density profile is thus

$$\begin{eqnarray*} \rho (r) = \rho _0 \cdot \exp \left(-\frac{\Phi _{\text{eq}}(r)}{\sigma ^2}\right), \end{eqnarray*}$$

(21)

with a constant ρ₀ if σ² is constant with radius.

In the presence of an anisotropic gravitational tidal field surrounding the halo, the halo will find a new dynamical equilibrium within the free fall time-scale |$\varpropto1/\sqrt{G\rho }$|⁠. As the halo is virialized, setting ρ = 200Ω_mρ_crit with |$\rho _{\mathrm{crit}} = 3H_0^2/(8\pi G)$| yields the fraction |$\sim \sqrt{4\pi /300\Omega _m}$| of the Hubble time 1/H₀ as the upper limit for the free fall time-scale. Taylor-expanding the anisotropic large-scale gravitational potential |$\Phi (\mathbf {r})$| to second order around the centre of mass of the halo gives

$$\begin{eqnarray*} \Phi (\mathbf {r}) = \Phi (\mathbf {x}_0) + \Phi _{,a}(\mathbf {x}_0)r_a + \frac{1}{2}\Phi _{,ab}(\mathbf {x}_0)r_ar_b + ..., \end{eqnarray*}$$

(22)

where we use the Einstein summation and _{, a} denotes derivation with respect to the position coordinate r_a. While the first derivative in equation (22) accelerates the halo, providing it with a non-zero peculiar velocity, shape modifications are effected by higher order derivatives (Piras et al. 2018; Ghosh, Durrer & Schäfer 2021; Giesel, Ghosh & Schäfer 2022). Provided σ² remains unscathed by the perturbation and the halo is small compared to the curvature scale |$\mathcal {S}$|⁠,

$$\begin{eqnarray*} \frac{1}{\mathcal {S}^2} := \bigg \vert \frac{1}{\sigma ^2}\frac{\mathrm{d}^2}{\mathrm{d}r^2} \Phi (r)\big \vert _{\mathbf {x}_0} \bigg \vert , \end{eqnarray*}$$

(23)

one can perturb the steady-state density profile ρ(r) of equation (21) to lowest order to obtain

$$\begin{eqnarray*} \tilde{\rho }(\mathbf {r}) \sim \rho (r)\left(1-\frac{\Phi _{,ab}(\mathbf {x}_0)r_ar_b}{2\sigma ^2}\right). \end{eqnarray*}$$

(24)

The first three even moments of the density distribution are given by

$$\begin{eqnarray*} S_0 := \int \mathrm{d}^3 r \rho (\mathbf {r}), \end{eqnarray*}$$

(25)

$$\begin{eqnarray*} q_{ij} := \frac{1}{S_0} \int \mathrm{d}^3 \rho (\mathbf {r})r_ir_j, \end{eqnarray*}$$

(26)

$$\begin{eqnarray*} s_{ijkl} := \frac{1}{S_0} \int \mathrm{d}^3 \rho (\mathbf {r})r_ir_jr_kr_l. \end{eqnarray*}$$

(27)

Equation (26) is the generalization of equation (1) for the continuous case. The perturbation of the second moment q_ij around |$\mathbf {r}_0$| due to cosmic tides is thus

$$\begin{eqnarray*} \Delta q_{ij}(\mathbf {x}_0) := \tilde{q}_{ij}(\mathbf {x}_0)-q_{ij} = -\frac{1}{2\sigma ^2}\Phi _{,ab}(\mathbf {x}_0)s_{ijkl}. \end{eqnarray*}$$

(28)

The (in the case of haloes) unobservable complex ellipticity in projection along the Cartesian z-axis follows from |$\tilde{q}_{ij} \equiv \tilde{q}_{ij}(\mathbf {x}_0)$| as

$$\begin{eqnarray*} \tilde{\epsilon } = \frac{\tilde{q}_{xx}-\tilde{q}_{yy}}{\tilde{q}}+i\frac{2\tilde{q}_{xy}}{\tilde{q}}, \end{eqnarray*}$$

(29)

where the halo size |$\tilde{q} = \tilde{q}_{xx}+\tilde{q}_{yy} = q$| is assumed to be constant. By inserting the tidal shape perturbation (28) we obtain

$$\begin{eqnarray*} \tilde{\epsilon } = \epsilon _0 - \frac{1}{2\sigma ^2q}\sum _{i,j \in \lbrace x,y\rbrace }\Phi _{,ij}(s_{ijxx} - s_{ijyy} + 2is_{ijxy}). \end{eqnarray*}$$

(30)

We note here that by symmetry (‘matching the indices’), equation (28) must hold to lowest order even if we do not assume an unperturbed halo that is spherically symmetric, |$q_{ij} = \frac{q}{2}\delta _{ij}$|⁠. For now we generalize the unperturbed halo to being perturbatively isotropic, i.e. |$q_{ij} = \frac{q}{2}\delta _{ij} + \eta _{ij}$|⁠, in which case the leading order contribution becomes

$$\begin{eqnarray*} \tilde{\epsilon } = \epsilon _0 - D(\Phi _{,xx} - \Phi _{,yy} + 2i\Phi _{,xy}) =: \epsilon _0 - D(T_{+} + iT_{\times }). \end{eqnarray*}$$

(31)

The constant

$$\begin{eqnarray*} D := \frac{c/q+q/4}{\sigma ^2} \end{eqnarray*}$$

(32)

absorbs the halo properties, i.e. the concentration or peakiness |$\tilde{c} = c$| which is the fourth cumulant of the density |$\rho (\mathbf {r})$|⁠, the velocity dispersion σ² and its size q.

In a statistical sample such as haloes from a simulation box, the intrinsic ellipticity ϵ can be modelled by a random variable with zero mean and some amount of dispersion, such that after dropping the tilde one obtains

$$\begin{eqnarray*} \epsilon = \epsilon _{+} + i\epsilon _{\times } = -D(T_{+} + iT_{\times }), \end{eqnarray*}$$

(33)

which constitutes the essence of LAM. Equation (33) is reminiscent of gravitational lensing, and indeed the analogy is justified as intrinsic alignment is governed by the second derivatives of Φ/σ² whereas lensing is governed by the second derivatives of Φ/c². D is strictly positive and is traditionally quoted in units of c², since it is beneficial to work with the dimensionless gravitational potential Φ/c².

To obtain D for our simulation samples, we first compute the discrete version of the ellipticity (29) for each halo. Next, we calculate the local gravitational tidal field

$$\begin{eqnarray*} \frac{\Phi _{,ij}}{c^2}(\mathbf {x}) = \frac{\partial ^2\Phi (\mathbf {x})}{c^2\partial x_i\partial x_j} \end{eqnarray*}$$

(34)

at the centre of mass of each halo. This can be achieved by first determining the 3D overdensity field |$\delta (\mathbf {x})$| for the full simulation volume using a CIC interpolation. We can then use a discrete Fourier transform to solve Poisson’s equation and obtain the Hessian of the potential algebraically in Fourier space via

$$\begin{eqnarray*} \frac{\Phi _{,ij}}{c^2}(\mathbf {x}) = \frac{3\Omega _m}{2\chi _H^2a}\mathcal {F}^{-1}\left[ \frac{k_ik_j}{|\mathbf {k}|^2} \exp \left(-\frac{1}{2}|\mathbf {k}|^2\lambda ^2\right)\mathcal {F}[\delta (\mathbf {x})]\right], \end{eqnarray*}$$

(35)

where |$\mathbf {k}$| is the comoving wave vector, |$\chi _H := \frac{c}{H_0}$| is the Hubble-distance and a the cosmic scale factor. There is a non-removable uncertainty on the scale the effective tidal field relevant for LAM in equation (33) should be evaluated on. Here, we consider a Gaussian smoothing scale of λ = 1 cMpc h⁻¹. This scale corresponds to haloes of mass |$M = 4\pi /3 \Omega _m \rho _{\mathrm{crit}} \lambda ^3 \sim 10^{11} \, \mathrm{M}_{\odot }\, h^{-1}$|⁠. Finally, to obtain the tidal shear at the halo positions, we perform an inverse CIC interpolation, effectively increasing the smoothing scale λ by about one grid cell.

Note that the derivation of equation (33) assumes that the objects under consideration are in a steady-state Jeans equilibrium with small perturbations. When fitting for D, we thus discard DM resolution elements that lie outside of the overdensity radius R₂₀₀ of the parent fof halo in an attempt to better satisfy the assumption. To reduce the bias in determining the complex ellipticity ϵ, we further restrict ourselves to central subhaloes, as we did when investigating mass density and shape profiles.

4.4.3 Isotropization of the ellipticity frame

As ϵ₊ and ϵ_× transform as components of a (spin-2) tensor, they are sensitive to the absolute orientation of structure. Since our simulation boxes of side lengths 10 and 40 cMpc h⁻¹ are not large enough to contain isotropic large-scale structure, cosmic variance will induce statistical differences in the parameters derived from the two components ϵ₊ and ϵ_×. To circumvent this, we follow a randomization technique outlined in Zjupa et al. (2021): It consists of randomizing the local xy-frame in which the alignment parameter D is measured by rotating the frame by a random angle between 0 and π, and averaging over the obtained results for D. The more randomizations one performs the more D as obtained from ϵ₊ will match the one obtained from ϵ_×. In practice it is sufficient to perform 10³ randomizations, the results of which we will denote by D_iso.

4.4.4 Linear alignment model fitting

In Fig. 10 we present some selected best-fitting results for the alignment parameter D. Note that choosing a smoothing scale of λ = 1 cMpc h⁻¹ while not imposing any fof halo mass cuts is slightly inconsistent, in that a smoothing scale of λ = 1 cMpc h⁻¹ would suggest imposing a halo mass floor of |$M_{\text{min}} \sim 10^{11} \, \mathrm{M}_{\odot }\, h^{-1}$|⁠. In turn, M_min determines the smoothing scale below which the matter power spectrum P(k) used for computing tidal fields should be cut off, since tidal field fluctuations on scales smaller than the scale corresponding to M_min cannot be relevant to the alignment process. The choice of the lower mass limit M_min matters when considering averaged values of virial quantities and, consequently, of the alignment parameter, the primary reason why the employment of a mass scale of |$10^{12}\, \mathrm{M}_{\odot }\, h^{-1}$| led Tugendhat & Schäfer (2018) to measure a large D-value of 10⁻⁴ (Mpc h⁻¹)². Since our focus lies on analysing trends with the cFDM particle mass m and redshift z, rather than exploring the whole parameter space that includes the smoothing scale λ and halo mass M, our lower mass limit M_min is merely determined by the requirement to reasonably resolve a halo. In other words, by discarding central subhaloes that are comprised of fewer than 400 particles, we obtain |$M_{\text{min}} = 400 \times m_{\mathrm{DM}} = 2.0 \times 10^9 \, \mathrm{M}_{\odot }\, h^{-1}$|⁠. For a comprehensive analysis on the subtle λ-dependence cf. Zjupa et al. (2021).

Figure 10.

Intrinsic alignment strengths in different cosmologies. We show the correlation of the two components of the halo ellipticity ϵ₊, ϵ_× with the respective tidal field components T₊, T_×, for N-body, 1024³ resolution, L_box = 40 cMpc h⁻¹ runs at redshift z = 4.38. Each dot represents one halo colour-coded by its triaxiality T. The blue dashed lines correspond to a linear fit to the binned data points, the shaded bands displaying the standard error on the mean (SEM) in each bin. The blue (CDM) and red solid lines (cFDM) depict the fits to the anisotropy-corrected data. All values are in units of c²(cMpc h⁻¹)².

We find that the D-values in Fig. 10 that are obtained from the real part of the halo ellipticity differ from the imaginary one by at most 2.6σ. In contrast, the isotropized values D_iso differ by at most 0.24σ. The latter ones are thus clearly more reliable, hence we present results for the isotropized values D_iso in Table 3, varying the redshift and the cosmology. We observe two important intrinsic alignment trends which we describe in Sections 4.4.5 and 4.4.6.

Table 3.

Results for the isotropized alignment parameter D_iso for various redshifts z (first column), in CDM (second column), and cFDM of various particle masses m (columns three, four, and five), obtained by averaging the isotropized linear fits from the two components ϵ₊, ϵ_×. Results are for N-body, 1024³ resolution, L_box = 40 cMpc h⁻¹ runs. The second column quotes D_iso values for CDM in units of c²(Mpc h⁻¹)², while the last three columns show the fractional differences |$\left(D^{\mathrm{cFDM}}_{\text{iso}}-D^{\mathrm{CDM}}_{\text{iso}}\right)/D^{\mathrm{CDM}}_{\text{iso}}$|⁠. The redshift spacings correspond to equal logarithmic spacings in the scale factor a, i.e. a_{z = 7.01}/a_{z = 10.90} = a_{z = 4.38}/a_{z = 7.01} ∼ 1.49.

Redshift	CDM	m = 2 × 10−21 eV	m = 7 × 10−22 eV	m = 10−22 eV
z = 10.90	1.01 × 10−6 ± 4 × 10−8	0.07 ± 0.05	0.08 ± 0.06	0.45 ± 0.19
z = 7.01	9.59 × 10−7 ± 2.4 × 10−8	0.06 ± 0.04	0.07 ± 0.03	0.58 ± 0.11
z = 4.38	8.67 × 10−7 ± 1.9 × 10−8	0.09 ± 0.03	0.12 ± 0.04	0.51 ± 0.08

4.4.5 Increase of alignment strength with redshift

As Table 3 suggests, we discern a strong dependence of D_iso on cosmic time, decreasing from D_iso ∼ 1.01 × 10⁻⁶c² (cMpc h⁻¹)² at redshift z = 10.90 to about D_iso ∼ 8.67 × 10⁻⁷c² (cMpc h⁻¹)² at redshift z = 4.38 in case of CDM. The same trend is exhibited by the cFDM cosmologies as well, with reductions in D_iso by about 15 per cent from z = 10.90 to z = 4.38. The origin of the z-dependence lies in the intricate responsivity of DM dynamics to changes in the surrounding tidal field. Such responsivity is believed to be impacted by the concentration c of the halo and other halo properties, but we defer such analysis to future work. At lower redshifts, the z-dependence is more accentuated, with D_iso decreasing by a factor of ∼3.7 from z = 1 to z = 0, see Zjupa et al. (2021) for an analogous analysis of the intrinsic alignment of the luminosity distribution of ellipticals in the V-band. The trend of an increasing alignment strength with redshift is also observed at low-z by Samuroff, Mandelbaum & Blazek (2021) in both MassiveBlack-II and the TNG300 simulation, as well as by Tenneti et al. (2015) for the non-linear alignment model in the MassiveBlack-II simulation.

4.4.6 Increase of alignment strength with cFDM particle mass

As we decrease the cFDM particle mass m, one noticeable change in Fig. 10 comes from the higher average in the triaxiality T of haloes. This is in agreement with our results in Section 4.3, where we have seen that cFDM haloes are more prolate at the virial radius R₂₀₀ than their CDM counterparts. By the same token, the distributions of both the real and the imaginary part of the ellipticity ϵ_× exhibit more outliers / heavier tails at larger absolute values than for CDM. This is expected, since prolate haloes tend to look – on average – prolate in projection as well. This circumstance together with the distribution of local tidal field values T₊, T_×, conspire to a strong dependence of the alignment strength on the cFDM particle mass m. The smaller m and thus the smaller the half-mode scale k_1/2, the more D_iso grows. The fractional difference in D_iso between the most extreme cFDM model with m = 10⁻²² eV and CDM is as high as 0.45 ± 0.19 at z = 10.9. At z = 4.38, the fractional difference reads 0.51 ± 0.08, which is different from zero at a level of about 6.4σ.

5 CONCLUSIONS

In this work, we show that the impact of a small-scale cutoff in the primordial power spectrum on high-redshift virialized structures such as haloes is profound, altering their density profiles, shapes, and intrinsic alignment. We analyse dark matter haloes in a suite of cosmological N-body simulations with the following specifications: two box sizes of L_box = 10 and 40 cMpc h⁻¹, DM resolutions of 256³, 512³, and 1024³, ΛCDM and classical FDM (cFDM) with particle masses m = 10⁻²², 7 × 10⁻²², 2 × 10⁻²¹ eV. cFDM ignores quantum pressure which is justified on the intermediate scales of k ∼ 2–18 Mpc h⁻¹ we are interested in assuming that gravitational interactions do not re-thermalize axions. We focus on three distinct properties of central subhaloes in this work:

• The spherically averaged internal mass distribution of dark matter haloes reflects the broken hierarchy of structure formation in cFDM in that the monotonicity in the concentration–mass relation c–M is also broken. Instead, the concentration peaks at around two orders of magnitudes above the half-mode mass M_1/2, a result that is fairly constant with redshift and in agreement with Bose et al. (2015), Ludlow et al. (2016). Universality of density profiles in cFDM is at best an approximation, since concentration does not depend on halo mass close to the concentration peak in cFDM cosmologies. This is similar to the plateau regime at the high-mass end (not captured with relaxed haloes) that is apparent in both CDM and cFDM cosmologies (Klypin et al. 2016). Concentration PDFs can be well fit by lognormal distributions in CDM as well as cFDM cosmologies in the investigated redshift range of 3.4 ≤ z ≤ 6.3.

• The asymmetry-sensitive shape parameters q (intermediate-to-major axis ratio) and s (minor-to-major axis ratio) exhibit a monotonicity relation in N-body simulations of ΛCDM: Both q and s are monotonously increasing as a function of ellipsoidal radius for pure N-body simulations, a monotonicity that is thus well maintained at high redshifts in ΛCDM. It is known that baryonic physics in the centres of haloes breaks the monotonicity below z ∼ 3, with radiative feedback processes sphericalizing the central regions (Chua et al. 2019, 2021). We now find that monotonicity is also broken by a small-scale cutoff in the primordial density fluctuations, in models such as FDM and cFDM. In the latter, using a large sample of DM haloes we show that they are significantly less spherical (lower q and s) and more prolate (higher triaxiality T) than their CDM counterparts close to and beyond the virial radius R_vir. Small-mass haloes experience the strongest shape distortions, mediated by the tidal fields of cosmic filaments into which many of them are embedded. Closer to the halo centre, we observe higher q-values for high-mass haloes and in addition lower s-values for intermediate mass haloes, translating into lower triaxialities T in both cases. Likewise, the anticorrelation between sphericity and halo mass that is observed in simulations of CDM is broken in cFDM cosmologies.

• The intrinsic alignment of haloes/ galaxies is both a blessing and a curse, the former due to the cosmological information that it contains and the latter by constituting a significant contaminant for past and future weak lensing surveys such as Euclid. We find that geometric measures for intrinsic alignment such as the shape–position and the shape–shape alignment statistics are sensitive to the cosmological model. In particular, for small 3D halo pair separations we show that median shape–position correlations |cos θ| increase with decreasing particle mass m and in the extreme cFDM model with m = 10⁻²² eV climb up to |$0.753^{+0.021}_{-0.010}$| as opposed to |$0.653^{+0.006}_{-0.002}$| in ΛCDM, at redshift z = 4.38 for intermediate halo masses of |$10^{10}-10^{11}\, \mathrm{M}_{\odot }\, h^{-1}$|⁠. We also carry out a linear alignment model analysis on the haloes across the full mass range in the various cosmologies. For the inferred isotropized linear alignment magnitudes D_iso, we confirm the ΛCDM-trend of larger D_iso values with increasing redshift and extend the result to cFDM. More importantly, we find very significant differences in D_iso between the cFDM models and ΛCDM, with the m = 10⁻²² eV model differing from ΛCDM by as much as 6.4σ at redshift z = 4.38.

We have seen that imprints of a primordial power spectrum cutoff on the internal structure and alignment of DM haloes are strikingly visible upon closer look at the properties of haloes in the high-z Universe. In this spirit, this work attempts to improve our understanding of the high-z Universe in the alternative FDM scenario. Except for the central regions of haloes in which quantum pressure gives rise to de Broglie wavelength-sized solitonic cores as demonstrated in the Supplementary Materials, the conclusions in this work are expected to hold for string theory-motivated FDM cosmologies. High-redshift observations are thus a promising avenue to provide evidence for / against FDM cosmologies.

In particular, on scales slightly larger than those considered in this work, the next generation of line intensity mapping (LIM) surveys is expected to be instrumental in putting constraints on FDM and mixed dark matter models, especially in view of upcoming experiments such as SPHEREx. While this spacecraft will observe near-infrared hydrogen recombination lines, 21 cm cosmology promises new insights at even higher redshifts well into Cosmic Dawn with the upcoming SKA telescope. Its low-frequency component SKA1-Low will be able to rule out a m ∼ 2.6 × 10⁻²¹ eV FDM-like cosmology at more than 2σ with 1000 h of observations at z ∼ 5 (Carucci et al. 2015). Combining an SKA1-Mid-like LIM survey with the CMB measurements at the future Simons Observatory, an m = 10⁻²² eV FDM model could be constrained at a few per cent (Bauer et al. 2021). Upcoming 21 cm power spectrum measurements by HERA will help determine the FDM particle mass to within 20 per cent at 2σ level (Jones et al. 2021), with even better sensitivity in the mass range 10⁻²⁵ eV ≤m ≤ 10⁻²³ eV (Flitter & Kovetz 2022). However, it remains to show how some of these forecasts can be extended to higher redshifts and how they are modified by non-linear FDM dynamics.

SUPPORTING INFORMATION

SupplementaryMaterials.pdf

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

It is a pleasure to thank Shihong Liao, who has provided powerful insights into low-magnitude power spectrum estimation techniques and the setting up of initial loads for simulations with a primordial power spectrum cutoff. We further thank Debora Sijacki and Nina Sartorio for fruitful discussions. We are also grateful to Sophie Koudmani and Martin Bourne for providing help on the ins and outs of arepo. TD acknowledges support from the Isaac Newton Studentship and the Science and Technology Facilities Council (STFC) under grant number ST/V50659X/1. MBK acknowledges support from NSF CAREER award AST-1752913, NSF grants AST-1910346 and AST-2108962, 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 is operated by AURA, Inc., under NASA contract NAS5-26555. AF is supported by the Royal Society University Research Fellowship. The simulations were performed under DiRAC project number ACSP253 using the Cambridge Service for Data Driven Discovery (CSD3), part of which is operated by the University of Cambridge Research Computing on behalf of the STFC DiRAC HPC Facility. 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.

DATA AVAILABILITY

We have published cosmicprofiles, a fast, pip- and conda-installable open-source implementation of various shape and density profile estimation and fitting algorithms in Cython. Its documentation can be found here. High-level data products are available upon reasonable request.

Footnotes

1

We will use the FDM particle mass m_FDM and the cFDM particle mass m_cFDM interchangeably and use a shorthand notation, m = m_FDM = m_cFDM.

2

The power spectra were obtained using nbodykit, a versatile open-source Python package (Hand et al. 2018).

3

Note that a redshift-independent T(k) is only a good approximation for m ≳ 10⁻²⁴ eV (Marsh 2016).

4

All of the mass bins in this work refer to fof masses.

5

Strictly speaking, the monotonicity can only be considered a universality if the dependence on at least the redshift and cosmological parameters of the q- and s-profiles can be subsumed via an appropriate rescaling. Since we do not investigate this claim thoroughly, we will henceforth refrain from calling the monotonicity a universality.

6

If we were to include baryons and radiative feedback processes, the steepening would be strongly impeded and the monotonicity eventually broken at lower redshifts (Chua et al. 2019, 2021), see also Section 1.2.