# Improved Mixed Dark Matter Halo Model for Ultralight Axions

Sophie M. L. Vogt1,2<EMAIL_ADDRESS>David J. E. Marsh3
<EMAIL_ADDRESS>Alex Laguë4,5,6,7<EMAIL_ADDRESS>1Institut für
Astrophysik und Geophysik, Universität Göttingen, Germany
2Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians-Universität
München, Scheinerstr. 1, 81679 München, Germany
3Theoretical Particle Physics and Cosmology, King’s College London, Strand,
London WC2R 2LS, United Kingdom
4Department of Physics and Astronomy, University of Pennsylvania,
Philadelphia, PA 19104, USA
5David A. Dunlap Department of Astronomy & Astrophysics, University of
Toronto, Canada
6Canadian Institute for Theoretical Astrophysics, University of Toronto,
Canada
7Dunlap Institute for Astronomy and Astrophysics, University of Toronto,
Canada

###### Abstract

We present a complete halo model for mixed dark matter composed of cold dark
matter (CDM) and ultralight axion-like particles (ULAs). Our model treats ULAs
as a biased tracer of CDM, in analogy to treatments of massive neutrinos and
neutral hydrogen. The model accounts for clustering of ULAs around CDM host
halos, and fully models the cross correlations of both components. The model
inputs include the ULA Jeans scale, and soliton density profile. Our model can
be used to predict the matter power spectrum, $P(k)$, on non-linear scales for
sub-populations of ULAs across the mass range $10^{-33}\text{ eV}\leq m\leq
10^{-21}\text{ eV}$, and can be calibrated against future mixed DM simulations
to improve its accuracy. The mixed DM halo model also allows us to assess the
importance of various approximations. The code is available at
https://github.com/SophieMLV/axionHMcode.

cosmology, large-scale structure

##  I Introduction

The standard $\Lambda$ cold dark matter ($\Lambda$CDM) cosmological model is
highly successful at describing the Universe Aghanim  _et al._ (2020a, b), yet
the microphysical nature of DM remains a mystery. DM candidates’ masses span
orders of magnitude from the heaviest primordial black holes to the lightest
sub-eV particles (see Ref. Alves Batista  _et al._ (2021)). Ultralight axion-
like particles (ULAs) in the mass range $10^{-33}\text{ eV}\leq m\leq
10^{-21}\text{ eV}$ are the lightest DM candidates, yet they are highly
constrained, with current data allowing only a few percent contribution to the
energy density Hlozek  _et al._ (2015a, 2018); Kobayashi  _et al._ (2017);
Laguë  _et al._ (2022). On the other hand, such particles are abundant in the
string theory landscape Witten (1984); Svrcek and Witten (2006); Conlon
(2006); Arvanitaki  _et al._ (2010), as has been shown recently in
increasingly explicit compactifications Mehta  _et al._ (2021); Demirtas  _et
al._ (2021); Cicoli  _et al._ (2022). Thus, one expects a sub-population of
the cosmic DM to be composed of ULAs. Next generation cosmological surveys
will increase the precision of ULA searches by orders of magnitude, and could
detect a sub-population of ULAs as small as $\mathcal{O}(0.1\%)$ Hložek  _et
al._ (2017); Bauer  _et al._ (2020); Farren  _et al._ (2022); Abazajian  _et
al._ (2016); Dvorkin  _et al._ (2022); Flitter and Kovetz (2022).

Exploiting to the full the next generation of cosmological data requires
considering cosmological statistics beyond the linear regime, and parameter
estimation from such data requires the non-linear physics to be computable
rapidly. Non-linear physics can be modelled extremely accurately using
$N$-body and hydrodynamical simulations (e.g. Refs. Vogelsberger  _et al._
(2020)), but such methods are not appropriate for parameter estimation. Two
methods that allow for fast estimation of non-linear observables are
_emulators_ and the _halo model_ (HM). Emulators are machine learning inspired
methods to interpolate accurately on large grids of simulations, and have been
employed to great effect in a variety of cases in $\Lambda$CDM and beyond
Heitmann  _et al._ (2013, 2016); Rogers  _et al._ (2019); Pedersen  _et al._
(2021).

The HM Cooray and Sheth (2002), on the other hand, has the advantage of
maintaining the speed of an emulator while being physics-inspired, and
elements of the HM can be calibrated on a smaller number of simulations.
HMCode is one example of the HM that competes with emulator accuracy on the
power spectrum, $P(k)$, with a small number of parameters calibrated to
simulation Mead  _et al._ (2015a). As an example, HMCode has been successfully
calibrated on fixed cosmology models of active galactic nucleus (AGN)
feedback, and then used to predict the cosmological parameter dependence Mead
_et al._ (2021a). The method of accurate HM calibration is also used by the
_Euclid_ consortium to model the non-linear power spectrum for cosmological
parameter estimation from clusters Castro  _et al._ (2022). Being physics-
inspired, the HM is thus highly suitable to apply to beyond $\Lambda$CDM
models, where suitable simulations might be limited in dynamic range, sparse
in parameter space, or non-existent.

In the following, we develop the mixed DM halo model for ULAs, greatly
improving on the early work in Ref. Marsh and Silk (2014). Our mixed HM draws
inspiration from the treatment of neutrino clustering around halos LoVerde and
Zaldarriaga (2014); Massara  _et al._ (2014), and the neutral hydrogen HM
Padmanabhan and Refregier (2017); Bauer  _et al._ (2020). Our model contains
physically motivated elements that are suitable to calibrate against mixed DM
simulations, when they reach the appropriate scales.

Throughout this work we refer to the ultralight sub-component of DM as
“axions”, since the existence of such particles is motivated in the string
theory landscape Arvanitaki  _et al._ (2010); Mehta  _et al._ (2021). However,
our treatment makes no assumption about the $CP$ properties of the ultralight
component, so it is equally valid for scalars (e.g. dilaton-like fields
Hamaide  _et al._ (2022)). Ultralight bosonic components with more degrees of
freedom, such as a complex scalar Li  _et al._ (2014) or vector Gorghetto  _et
al._ (2022), are amenable to the same modelling (possessing a Jeans scales,
and forming solitons) but differ in details that require treating separately
(e.g. different linear transfer function and early Universe behaviour).

This paper is organised as follows. In Section II, we briefly outline the
relevant aspects of ULA cosmology. In Section III we develop the theory of the
HM in the case of mixed DM. In Section IV we give the results of our model in
terms of the non-linear power spectrum. We conclude in Section V. Appendix A
discusses the modifications to the HM of HMCode, which we adopt as baseline.
Appendix B shows the mixed DM halo model for massive neutrinos, following Ref.
Massara  _et al._ (2014), and assesses the effect of neutrino clustering in
halos compared to the approximate treatment in HMCode (we find differences of
order a few percent at $k\approx 1\,h\text{ Mpc}^{-1}$). Appendix C shows some
convergence checks on our numerical implementation. We take baseline
cosmological parameters as shown in Table 1.

Parameter | Fiducial value
---|---
$h$ | 0.674
$\omega_{\mathrm{d}}$ | 0.12
$\omega_{\mathrm{b}}$ | 0.02237
$f_{\mathrm{ax}}$ | 0.1
$N_{\textnormal{eff}}$ | 3.046
$A_{s}$ | $2.1\times 10^{-9}$
$n_{s}$ | 0.97
$k_{\textnormal{piv}}$ [$\textnormal{Mpc}^{-1}$] | 0.05
Table 1: Fiducial cosmological parameters and their values for our flat mixed
dark mater (MDM) cosmology to compute the linear and non-linear power
spectrum. $h$ is the Hubble parameter; $\omega_{\mathrm{d}}$ and
$\omega_{\mathrm{b}}$ are the dark matter and baryon reduced density
parameters, respectively,
$f_{\mathrm{ax}}=\omega_{\mathrm{ax}}/\omega_{\mathrm{d}}$ the axion fraction
, $N_{\textnormal{eff}}$ is the effective number of neutrinos, $A_{s}$ is the
scalar amplitude, $n_{s}$ is the scalar spectral index and
$k_{\textnormal{piv}}$ is the pivotal scale.

##  II Axion Physics

The halo model aims to predict the non-linear power spectrum given as input
the linear power spectrum. We begin this section describing linear
perturbation theory for ULAs, and key results (many more details can be found
in e.g. Refs. Marsh (2016a); Hlozek  _et al._ (2015a) and references therein).
We then briefly describe the non-linear theory that underlies the mixed CDM-
ULA simulations of Refs. Schwabe  _et al._ (2020); Lague  _et al._ (2022),
which gives the form of the halo density profiles we adopt.

###  II.1 Perturbation Theory

Relic ULAs can be produced by the misalignment mechanism of a classical field
Preskill  _et al._ (1983); Abbott and Sikivie (1983); Dine and Fischler
(1983). The scalar field $\phi$ obeys the classical Klein-Gordon-Equation:

$\Box\phi-m_{\textnormal{ax}}^{2}\phi=0\,,$ (1)

where $\Box$ is the D’Alembertian operator for spacetime metric $g_{\mu\nu}$,
and we have assumed small displacements from the vacuum. Consider cosmological
perturbation theory in the Newtonian gauge Ma and Bertschinger (1995). At
zeroth order in perturbations, the axion field obeys the ordinary differential
equation:

$\ddot{\phi}_{0}+3H\dot{\phi}_{0}+m_{\textnormal{ax}}^{2}\phi_{0}=0\,,$ (2)

where $\phi_{0}$ is a function only of time. Here $H$ denotes the Hubble
parameter defined by $H=\dot{a}/a$, dots denote derivatives with respect to
cosmic time, $t$, and $a$ is the cosmic scale factor.

Eq. (2) is the equation for a damped harmonic oscillator. As long as $m\ll H$
(early times) the field is “frozen”, i. e. the harmonic oscillator is
overdamped. As the universe evolves $H\sim m_{a}$ and the field starts to
oscillate. Since the axion field is overdamped to begin with, the initial
conditions at time $t_{i}$ can be set to $\phi_{0}(t_{i})=\phi_{i}$ and
$\dot{\phi}(t_{i})=0$.

From the full Klein-Gordon-Equation, Eq. (1), we can compute the equation of
motion for the axion overdensity $\delta_{\textnormal{ax}}$ and thus find the
(linear) matter power spectrum with ULAs. In the Newtonian gauge the perturbed
field equation reads Marsh (2016a)

$\delta\phi^{\prime\prime}+2\mathcal{H}\delta\phi^{\prime}+(k^{2}+m_{\textnormal{ax}}^{2}a^{2})\delta\phi=(\Psi^{\prime}+3\Phi)\phi^{\prime}-2m_{\textnormal{ax}}^{2}a^{2}\phi\Psi\,,$
(3)

with the two Newtonian scalar potentials $\Phi$ and $\Psi$ and primes denotes
derivative with respect to conformal time $\tau$
($\textnormal{d}t=a(t)\textnormal{d}\tau$). We have assumed
$\delta_{\textnormal{ax}}\ll 1$ and $\delta\phi\ll\phi_{0}$.

When the axion field $\phi$ starts to oscillate (when the damping term in Eq.
2 becomes less than the mass, i.e. $H\lesssim m_{\rm ax}$) one can find with
the WKB-Ansatz an expression for the effective axion sound speed chan Hwang
and Noh (2009):

$c_{\textnormal{s}}^{2}=\frac{\frac{k^{2}}{4m_{\textnormal{ax}}^{2}a^{2}}}{1+\frac{k^{2}}{4m_{\textnormal{ax}}^{2}a^{2}}}\,.$
(4)

With the sound speed we obtain for the equation of motion for the axion
overdensity in the Newtonian gauge Marsh (2016a)

$\delta^{\prime\prime}+\mathcal{H}\delta^{\prime}+c_{\textnormal{s}}^{2}k^{2}\delta-3\mathcal{H}\Phi^{\prime}+k^{2}\Psi-3\Phi^{\prime\prime}=0\,.$
(5)

Compared with the equation of motion for the CDM overdensity, Eq. (5) has an
extra term proportional to the sound speed. This term goes to zero as $k\to 0$
and thus ULAs behave like CDM on large scales. For small scales the sound
speed is no longer negligible and the ULA overdensity oscillates instead of
growing. This behaviour is different from CDM, since the CDM overdensity has a
growing solution on all scales. The Jeans scale is the approximate scale where
the transition between the two regimes takes place Khlopov  _et al._ (1985).

If $k\ll m_{\textnormal{ax}}a$ the sound speed reads:

$c_{\textnormal{s}}^{2}=\frac{k^{2}}{4m_{\textnormal{ax}}^{2}a^{2}}\,.$ (6)

Then the Jeans scale reads Marsh (2016a)

$k_{\textnormal{J}}=66.5a^{1/4}\left(\frac{\Omega_{\textnormal{ax}}h^{2}}{0.12}\right)^{1/4}\left(\frac{m_{\textnormal{ax}}}{10^{-22}\,\textnormal{eV}}\right)^{1/2}\,\textnormal{Mpc}^{-1}\,,$
(7)

with $\bar{\rho}_{\textnormal{ax}}=\rho_{{\textnormal{ax}},0}a^{-3}$.
Suppression of the ULA power spectrum relative to CDM begins at the Jeans
scale at matter-radiation equality, $k_{\rm J,eq}$ Hu  _et al._ (2000); Bauer
_et al._ (2020).

With the above derivation we can summarise the effect on the matter power
spectrum if the DM is a mixture of CDM and ULAs:

  1. 1.

The axion overdensity behaves exactly as the CDM overdensity for large scales,
i. e. $k<k_{\textnormal{J}}$ and thus there is no change in the matter power
spectrum for these scales.

  2. 2.

A suppression by a factor of
$(1-\Omega_{\textnormal{ax}}/\Omega_{\textnormal{m}})^{2}$ for
$k>k_{\textnormal{J}}$ because the axion oscillates on these scales and the
overdensity is negligible compared to $\delta_{\textnormal{c}}$ and
$\delta_{\textnormal{b}}$.

  3. 3.

An extra suppression again for $k>k_{J}$ due to the fact that the axion field
suppresses the growing solution of the cold dark matter field,
$\delta_{\textnormal{c}}\propto
a^{\frac{-1+\sqrt{25-24\Omega_{\textnormal{ax}}/\Omega_{\textnormal{m}}}}{4}}$.

These effects lead to the presence of a step-like feature in the power
spectrum, which begins near the Jeans scale at equality, and has an amplitude
determined by the ULA fraction relative to CDM (see Refs. Amendola and
Barbieri (2006); Marsh and Ferreira (2010); Arvanitaki  _et al._ (2010)).

To illustrate the described behaviour, Fig. 1 shows the matter power spectrum
with an axion mass $m_{\textnormal{ax}}=10^{-28}$ eV and an axion fraction
$0.05<\Omega_{\textnormal{ax}}/\Omega_{\textnormal{d}}\leq 0.25$ (top) and the
ratio to the power spectrum in a $\Lambda$CDM universe (bottom) computed with
the Boltzmann code axionCAMB Hlozek  _et al._ (2015b), which solves first
order perturbation theory for ULAs in synchronous gauge coupled to all other
$\Lambda$CDM components, taken here with adiabatic initial conditions in a
radiation dominated Universe.

Figure 1: The linear matter power spectrum computed with axionCAMB with axions
of a mass $m_{\textnormal{ax}}=10^{-28}$ eV and an axion fraction
$0.05<\Omega_{\textnormal{ax}}/\Omega_{\textnormal{d}}\leq 0.25$ on the top.
On the bottom the ratio between the different matter power spectra with axions
to a $\Lambda$CDM power spectrum are shown. The dashed, vertical grey lines
indicate $k_{\textnormal{J}}$ for the lowest axion fraction (right line) and
highest fraction (left line) at $a_{\textnormal{eq}}$.

###  II.2 Non-Linear Regime

Taking Eq. (1), we insert the WKB ansatz:

$\phi=\frac{1}{\sqrt{2m_{\rm ax}^{2}}}\left(\psi e^{im_{\rm
ax}t}+\psi^{*}e^{-im_{\rm ax}t}\right)\,,$ (8)

working to linear order in perturbations of the metric, and taking the non-
relativistic limit, we find the Schrödinger-Poisson equations describing the
complex field $\psi$:

$\displaystyle i\frac{\partial}{\partial{t}}\psi$
$\displaystyle=-\frac{1}{2m}\nabla^{2}\psi+m_{\rm ax}\Phi\psi,$ (9)
$\displaystyle\nabla^{2}\Phi$ $\displaystyle=4\pi G\left[m_{\rm
ax}\left(|\psi|^{2}-\langle|\psi|\rangle^{2}\right)+\delta\rho_{\rm
f}\right],$ (10)

where $G$ is Newton’s constant, $\nabla$ is the flat space Laplacian, angle
brackets denote spatial average, and $\delta\rho_{\rm f}$ denotes fluid
density perturbations (e.g. CDM and baryons). This system of equations makes
no assumptions about the smallness of density perturbations, and can be used
to evolve ULAs into the non-linear regime while fully capturing wavelike
dynamics Widrow and Kaiser (1993); Schive  _et al._ (2014a); Mocz  _et al._
(2017a).

The key features of ULAs in the non-linear regime are: the persistence of the
Jeans scale, i.e. effective pressure opposing gravitational collapse, leading
to the existence of stable solitons Schive  _et al._ (2014a), condensation due
to wave scattering Levkov  _et al._ (2018), relaxation Mocz  _et al._ (2017b);
Dalal and Kravtsov (2022), and interference effects in the multi streaming
regime Schive  _et al._ (2014a); Mocz  _et al._ (2019); Gough and Uhlemann
(2022). Many different numerical approaches have been adopted to solve the
Schrödinger-Poisson equations, with different methods being useful in
different regimes of physical interest (e.g. Refs. Nori and Baldi (2018);
Veltmaat and Niemeyer (2016); Chen  _et al._ (2021)).

The code axioNyx (based on the Nyx code Almgren  _et al._ (2013)) was
presented and released in Ref. Schwabe  _et al._ (2020), and has been used to
simulate mixed ULA DM. These simualtions inform our MDM halo model. The MDM
simulations of Ref. Schwabe  _et al._ (2020) studied spherical collapse, and
noted that the ULAs follow CDM on large scales, while forming solitons on
smaller scales. Ref. Lague  _et al._ (2022) has extended the study of mixed DM
to cosmological initial conditions, although box sizes limit the ability to
measure the power spectrum and halo mass function. In order to properly evolve
the ULA wavefunction, the simulation grid spacing must be smaller than half
the de Broglie wavelength $\lambda_{\rm dB}\propto 1/m_{\rm ax}v$ of the ULAs
May and Springel (2021). In the presence of a dominant CDM component
($\Omega_{\textnormal{ax}}/\Omega_{\textnormal{d}}\lesssim 0.5$), the
potential wells in which the wavefunction evolves become steeper, which
increases the ULA velocity dispersion and decreases its wavelength. The mixed
DM simulations thus require a higher resolution than their pure ULA
counterparts. This makes running large cosmological volumes involving solving
the full Schrödinger-Poisson with the MDM system difficult. Alternative MDM
simulation algorithms based on a particle treatment of the ULAs (such as the
one implemented in ax-gadget Nori and Baldi (2018)) could capture the large-
scale dynamics of MDM and be a complement to the spectral and finite
difference methods at a reduced computational cost.

##  III Halo Model Theory

In this section we explore the halo model (HM) as a theoretical approach to
calculate the non-linear matter power spectrum. We first introduce the halo
model in a standard $\Lambda$CDM cosmology, as reviewed by Ref. Cooray and
Sheth (2002), in Section III.1. We then extend the HM to a mixed DM cosmology
in Section III.2, following a biased tracer treatment. The specific HM
ingredients in the case of ULAs are shown in Section III.3.

###  III.1 The $\Lambda$CDM Halo Model

In a universe where all matter is assumed to be cold we can assume that all
matter is contained in halos and thus the matter power spectrum is the sum of
two terms, the one halo term $P^{\mathrm{1h}}$ (correlation in the same halo)
and the two halo term $P^{\mathrm{2h}}$ (correlation between two different
halos) Cooray and Sheth (2002)

$P(k)=P^{\mathrm{1h}}(k)+P^{\mathrm{2h}}(k)\,,$ (11)

where the one halo and two halo terms have the following forms

$P^{\mathrm{1h}}(k)=\frac{1}{\bar{\rho}^{2}}\int\mathrm{d}MM^{2}n(M)|\tilde{u}(k,M)|^{2}\,,$
(12)

and

$P^{\mathrm{2h}}(k)=P^{\mathrm{L}}(k)\left[\frac{1}{\bar{\rho}}\int\mathrm{d}MMn(M)b(M)|\tilde{u}(k,M)|\right]^{2}\,.$
(13)

The halo mass function (HMF) is $n(M)$, $b(M)$ is the halo bias,
$\tilde{u}(k,M)$ is the Fourier transform of the halo density profile and
$P^{\mathrm{L}}(k)$ is the linear matter power spectrum. These are defined in
the following.

The improper integral in the two halo term, Eq. (13), should go to unity if
$k\to 0$ because matter is unbiased with respect to itself Cooray and Sheth
(2002). So, to ensure the correct behaviour for low $k$’s the mass interval
should be chosen large enough. This numerical problem was studied in Ref. Mead
_et al._ (2020) and was solved by adding some correction factors. The
implications on the halo model are discussed in Appendix C.

The Fourier transform of a (normalised) radial density profile
$u(r,M,z)=\rho(r,M,Z)/M$ is given by

$\tilde{u}(k,M,z)=4\pi\int_{0}^{r_{\textnormal{v}}}u(r,M,z)\frac{\sin(kr)}{kr}r^{2}\textnormal{d}r\,.$
(14)

Here the profile is truncated at the virial radius $r_{\textnormal{v}}$, i. e.
we assume that the density profile is zero for $r>r_{\textnormal{v}}$ and that
the mass of the halo is given by

$M=\frac{4\pi}{3}\bar{\rho}\Delta_{\textnormal{v}}(z)r_{\textnormal{v}}^{3}\,,$
(15)

with $\Delta_{\textnormal{v}}$ the virial overdensity (see below). Still
assuming a $\Lambda$CDM universe the density profile of a dark matter halo can
be described by the Navarro-Frenk-White (NFW) profile Navarro  _et al._ (1997)

$\rho_{\textnormal{NFW}}(r,M)=\frac{\rho_{\mathrm{char}}}{r/r_{\textnormal{s}}\left(1+r/r_{\textnormal{s}}\right)^{2}}\,,$
(16)

with $r_{\textnormal{s}}$ the scale radius and $\rho_{\textnormal{char}}$ the
characteristic density of the profile which ensures that the integral over the
NFW density profile gives the enclosed mass in Eq. (15) and can be computed to
be:

$\rho_{\mathrm{char}}=\rho_{\textnormal{crit}}\frac{\Omega_{\textnormal{m}}(z)\Delta_{\textnormal{v}}c^{3}}{3f(c)}\,,$
(17)

with $f(x)=-\frac{x}{1+x}+\ln(x+1)$ and
$c=r_{\textnormal{v}}/r_{\textnormal{s}}$ is the halo concentration parameter.

Evaluating the Fourier transformation, Eq. (14), with the NFW-profile gives
Scoccimarro  _et al._ (2001):

$\displaystyle\tilde{u}(k,M,z)$
$\displaystyle=\frac{1}{f(c)}\bigg{(}\cos(b)(\mathrm{Ci}(b+kr_{\textnormal{v}})-\mathrm{Ci}(b))$
$\displaystyle\hskip
2.84544pt+\sin(b)(\mathrm{Si}(b+kr_{\textnormal{v}})\left.-\mathrm{Si}(b))-\frac{\sin(kr_{\textnormal{v}})}{b+kr_{\textnormal{v}}}\right)\,.$
(18)

Here $b=kr_{\textnormal{v}}/c$ and Si$(x)$ and Ci$(x)$ are the sine and cosine
integrals. To calculate the halo mass function (HMF) $n(M)$ and the halo bias
$b(M)$ we need the variance of the linear power spectrum

$\displaystyle\sigma(R)^{2}=\frac{1}{2\pi^{2}}\int_{0}^{\infty}P^{\mathrm{L}}(k,z)\tilde{W}(Rk)^{2}k^{2}\textnormal{d}k$
(19) $\displaystyle\hskip
28.45274pt\mathrm{with}\quad\tilde{W}(x)=\frac{3}{x^{3}}(\sin\,x-x\cos\,x).$

Here we assumed a spherical top hat window function, $W$, in real space. The
variance above can be transformed to a function of the halo mass by
$M=4/3\pi\bar{\rho}R^{3}$ . The halo mass function, $n(M)$, is given by Press
and Schechter (1974):

$n(M,z)=\frac{1}{M}\frac{\textnormal{d}\tilde{n}}{\textnormal{dln}M}=\frac{1}{2}\frac{\bar{\rho}(z)}{M^{2}}f(\nu)\left|\frac{\mathrm{dln}\sigma^{2}}{\mathrm{dln}M}\right|\,.$
(20)

where $\tilde{n}$ is the halo number density, $\nu=\delta_{\rm
crit}(z)/\sigma(M,z)$, with $\delta_{\rm crit}$ the critical linear density
threshold for halo collapse, and $f(\nu)$ is the multiplicity function, which
for ellipsoidal collapse is given by the Sheth-Tormen (ST) multiplicity
function Sheth and Tormen (2002)

$f_{\textnormal{ST}}(\nu)=A\sqrt{\frac{2}{\pi}}\sqrt{q}\nu(1+(\sqrt{q}\nu)^{-2p})e^{-\frac{q\nu^{2}}{2}}\,,$
(21)

with $A=0.3222,\ p=0.3,\ q=0.707$. The last term in Eq. (20) can be calculated
by inserting the definition of the variance $\sigma(M,z)$ and using
$x\mathrm{d}x=\mathrm{dln}x$:

$\displaystyle\frac{\mathrm{d\,ln}\sigma^{2}}{\mathrm{d\,ln}M}$
$\displaystyle=\frac{3}{\sigma^{2}R^{4}\pi^{2}}\int_{0}^{\infty}\mathrm{d}k\frac{P^{\mathrm{L}}(k)}{k^{2}}\tilde{I}(k,R)$
(22) $\displaystyle\tilde{I}(k,R)$ $\displaystyle=(\sin(kR)-kR\cos(kR))$
$\displaystyle\hskip
2.84544pt\left[\sin(kR)\left(1-\frac{3}{(kr)^{2}}\right)+\frac{3}{kR}\cos(kR)\right].$
(23)

The halo bias can computed with the theory of Sheth and Tormen (2002) to be

$b(m,z)=1+\frac{1}{\delta_{\rm
crit}}\left(q\nu^{2}-1+\frac{2p}{1+(\sqrt{q}\nu)^{2p}}\right)\,.$ (24)

It remains to specify two quantities: the virial overdensity
$\Delta_{\mathrm{v}}$, and the concentration parameter, $c(M)$, of the NFW
profile. The first one is found for a $\Lambda$CDM cosmology by simulations
and a fitting formula was constructed which reads Bryan and Norman (1998)

$\Delta_{\textnormal{v}}(z)=\frac{18\pi^{2}+82x-39x^{2}}{\Omega_{\textnormal{m}}(z)}\,,$
(25)

with $x=\Omega_{\textnormal{m}}(z)-1$.

To find a functional form of the concentration parameter we follow Ref.
Bullock  _et al._ (2001) and assign for each halo of mass $M$ at redshift $z$
a formation redshift $z_{\textnormal{f}}$ by the equation

$M_{*}(z_{\textnormal{f}})=0.01M\,,$ (26)

where $M_{*}(z)$ is the collapsing mass defined by
$\sigma(M_{*}(z))=\delta_{\rm crit}(z)$. Here we assumed that $\delta_{\rm
crit}$ is constant and is given by $\delta_{\rm crit}=1.686$. An equivalent
definition for the formation redshift is given by:

$\frac{D(z_{\textnormal{f}})}{D(z)}\sigma(0.01M,z)=\delta_{\rm crit}(z)\,,$
(27)

here $D(z)$ is the linear growth factor
$D(z)/D(0)=\delta_{L}(\mathbf{x},z)/\delta_{L}(\mathbf{x},0)$. With the
equations above and fits to simulations the concentration parameter is Bullock
_et al._ (2001)

$c(M,z)=4\left(\frac{1+z_{\textnormal{f}}(M,z)}{1+z}\right)\,.$ (28)

Note that the defining equation for the formation redshift Eq. (27) can give
$z_{\textnormal{f}}<z$. In this case the authors of Bullock  _et al._ (2001)
forced $z_{\textnormal{f}}=z$ and the concentration parameter is $c(M,z)=4$.
Thus, the prefactor in Eq. (28) is called the minimum halo concentration.

###  III.2 Mixed Dark Matter Halo Model

In any MDM model where the clustering properties of the components differ, the
HM is more complex. As we saw in Section II.1 ULAs cannot cluster on small
scales, i. e. smaller than the Jeans scale and thus the assumption that all
matter is contained in halos is no longer valid. We also expect that the
Schrodinger-Poisson equation will cause the internal structure of ULA density
profiles to depart from the CDM profile on small scales. We consider MDM
models where the ULA component is a sub-dominant component of the DM. We
follow closely the model for massive neutrinos plus CDM (mixed hot and cold
DM) of Ref. Massara  _et al._ (2014), which we reproduce in more detail in
Appendix B (the authors of this paper developed a general approach of a MDM
halo model with a sub-domiant component as a biased tracer of the CDM and thus
their method is applicable in our analysis).

One of the main assumption in Ref. Massara  _et al._ (2014) is that the non-
cold component of the DM, i. e. massive neutrinos, clusters in the potential
wells of the cold matter halos, see e. g. Refs. Singh and Ma (2003); Ringwald
and Wong (2004); LoVerde and Zaldarriaga (2014). Thus the non-cold component
is treated as a biased tracer of the cold matter (for simplicity we neglect
the non-trivial clustering of baryons on the scales of interest). A similar
biased tracer technique is used in the halo model of neutral hydrogen, e. g.
see Refs. Villaescusa-Navarro  _et al._ (2018); Bauer  _et al._ (2020), which
can match well to simulations. ULAs are expected to behave in the same way:
tracing the dominant cold matter in halos above the Jeans scale with some
characteristic internal density profile.

Since the matter power spectrum is proportional to the matter overdensity
squared, $P(k)\propto\delta_{\textnormal{m}}^{2}$, the total matter
overdensity in a MDM cosmology is a sum of the cold matter
$\delta_{\textnormal{c}}$ and ULAs $\delta_{\textnormal{a}}$:

$\delta_{\textnormal{m}}=\frac{\Omega_{\textnormal{c}}}{\Omega_{\textnormal{m}}}\delta_{\textnormal{c}}+\frac{\Omega_{\textnormal{a}}}{\Omega_{\textnormal{m}}}\delta_{\textnormal{a}}\,,$
(29)

with $\Omega_{\textnormal{c}}$ the weighted sum of the density parameters of
CDM and baryons and $\Omega_{\textnormal{a}}$ the density parameter of axions.
Then the power spectrum reads: 111Note that we consider only adiabatic
perturbations, where the cross correlation $P_{\rm c,a}\propto\delta_{\rm
c}\delta_{\rm a}$.

$P_{\textnormal{m}}(k)=\left(\frac{\Omega_{\textnormal{c}}}{\Omega_{\textnormal{m}}}\right)^{2}P_{\textnormal{c}}(k)+\frac{2\Omega_{\textnormal{c}}\Omega_{\textnormal{a}}}{\Omega^{2}_{\textnormal{m}}}P_{\textnormal{c,a}}(k)+\left(\frac{\Omega_{\textnormal{a}}}{\Omega_{\textnormal{m}}}\right)^{2}P_{\textnormal{a}}(k)\,,$
(30)

where $P_{\textnormal{c}}(k)$, $P_{\textnormal{c,a}}(k)$ and
$P_{\textnormal{a}}(k)$ are the cold matter, cross and ULA power spectrum
respectively. We can already see from the prefactors of the above equation
that the main contribution to the non-linear total matter power spectrum comes
from the cold part, because we assume
$\Omega_{\textnormal{a}}\ll\Omega_{\textnormal{c}}$.

For the cold matter power spectrum we can use the standard halo model
described in Section III.1 and thus $P_{\textnormal{c}}(k)$ is given by
Equation (11), since for cold matter we still assume that all (cold) matter is
bound into halos.

Next we have to find expressions for the cross and non-cold non-linear power
spectrum. As explained above ULAs have a component that cannot cluster and
thus evolve approximately linearly. Thus, the axion overdensity has a
component in halos and a linear component, and can be written as:

$\delta_{\textnormal{a}}=F_{\textnormal{h}}\delta_{\textnormal{h}}+(1-F_{\textnormal{h}})\delta_{\textnormal{L}}\,.$
(31)

Here $\delta_{\textnormal{h}}$ and $\delta_{\textnormal{L}}$ are the halo and
linear overdensities, respectively, and $F_{\textnormal{h}}$ is the ULA
fraction in halos, i. e. $F_{\textnormal{h}}\in[0,1]$. With this overdensity
the cross and ULA power spectra read Massara  _et al._ (2014)

$\displaystyle P_{\textnormal{c, a}}(k)=F_{\textnormal{h}}P_{\textnormal{c,
a}}^{\textnormal{h}}(k)+(1-F_{\textnormal{h}})\sqrt{P_{\textnormal{c}}(k)P_{\textnormal{a}}^{\textnormal{L}}(k)}\,,$
(32) $\displaystyle
P_{\textnormal{a}}(k)=F_{\textnormal{h}}^{2}P_{\textnormal{a}}^{\textnormal{h}}(k)+2F_{\textnormal{h}}(1-F_{\textnormal{h}})\sqrt{P_{\textnormal{a}}^{\textnormal{h}}(k)P_{\textnormal{a}}^{\textnormal{L}}(k)}$
$\displaystyle\hskip
28.45274pt+(1-F_{\textnormal{h}})^{2}P_{\textnormal{a}}^{\textnormal{L}}(k)\,,$
(33)

where $P_{\textnormal{c, a}}^{\textnormal{h}}$ and
$P_{\textnormal{a}}^{\textnormal{h}}$ are the cross and axion non-linear power
spectra, respectively, and $P_{\textnormal{c,a}}^{\textnormal{L}}$ and
$P_{\textnormal{a}}^{\textnormal{L}}$ the corresponding linear power spectra.
The linear power spectra can be computed directly with axionCAMB Hlozek  _et
al._ (2015b) by using the transfer functions and the primordial power
spectrum, whereas for the non-linear part we use the HM.

To find the expression for the cross and axion non-linear power spectrum the
biased tracer technique plays an important role. As mentioned above this
method assumes that axion halos only form around cold matter halos (c-halos)
and thus the halo mass function for ULAs is the same as for the cold field,
$n(M_{\rm ax})\textnormal{d}M_{\rm
ax}=n(M_{\textnormal{c}})\textnormal{d}M_{\textnormal{c}}$, and the linear
axion halo bias corresponds to the c-halo bias, $b(M_{\rm
ax})=b(M_{\textnormal{c}})$ (i. e. the axion halo mass is itself a function of
the c-halo mass).

Another quantity which also will change is the exact form of the halo mass
function in Equation 20 which can be rewritten in the following form

$n_{\textnormal{c}}\textnormal{d}M_{\textnormal{c}}=\frac{\bar{\rho}_{\textnormal{tot}}}{M}f(\nu)\frac{1}{\nu}\textnormal{d}\nu\,,$
(34)

with $M=M_{\textnormal{c}}+M_{\textnormal{a}}$ and
$\nu=\delta_{\textnormal{c}}/\sigma(M)$. In the case of a MDM with massive
neutrinos N-body simulations showed that the HMF is fitted much better if only
the cold matter field is used in the peak height variable $\nu$ Villaescusa-
Navarro  _et al._ (2014), Castorina  _et al._ (2014), Costanzi  _et al._
(2013). The authors of this paper series also made predictions for the halo
bias and showed that also for the bias only the cold matter field should be
used. We assume that this continues to be true for other biased tracers such
as ULAs.

Computing the cross and axion non-linear power spectra with the halo model
with the biased tracer technique as well as the cold matter description gives
for the one and two halo term of the cross power

$\displaystyle P_{\textnormal{c,a}}^{\textnormal{1h}}(k)$
$\displaystyle=\frac{1}{F_{\textnormal{h}}\bar{\rho}_{\textnormal{c}}\bar{\rho}_{\textnormal{a}}}\int_{M_{\textnormal{cut}}}^{\infty}\textnormal{d}M_{\textnormal{c}}n(M_{\textnormal{c}})M_{\textnormal{c}}M_{\textnormal{a}}(M_{\textnormal{c}})$
$\displaystyle\hskip
71.13188pt\tilde{u}_{\textnormal{c}}(k,M_{\textnormal{c}})\tilde{u}_{\textnormal{a}}(k,M_{\textnormal{c}})\,,$
(35) $\displaystyle P_{\textnormal{c, a}}^{\textnormal{2h}}(k)$
$\displaystyle=\frac{1}{F_{\textnormal{h}}\bar{\rho}_{\textnormal{c}}\bar{\rho}_{\textnormal{a}}}\int_{0}^{\infty}\textnormal{d}M_{\textnormal{c}}n(M_{\textnormal{c}})b(M_{\textnormal{c}})M_{\textnormal{c}}\tilde{u}_{\textnormal{c}}(k,M_{\textnormal{c}})$
$\displaystyle\times\int_{M_{\textnormal{cut}}}^{\infty}\textnormal{d}M_{\textnormal{c}}n(M_{\textnormal{c}})b(M_{\textnormal{c}})M_{\textnormal{a}}(M_{\textnormal{c}})\tilde{u}_{\textnormal{a}}(k,M_{\textnormal{c}})P_{\textnormal{c}}^{\textnormal{L}}(k)\,$
(36)

and for the one and two halo of the ULA power:

$\displaystyle P_{\textnormal{a}}^{\textnormal{1h}}(k)$
$\displaystyle=\frac{1}{(F_{\textnormal{h}}\bar{\rho}_{\textnormal{a}})^{2}}\int_{M_{\textnormal{cut}}}^{\infty}\textnormal{d}M_{\textnormal{c}}n(M_{\textnormal{c}})M_{\textnormal{a}}(M_{\textnormal{c}})^{2}\tilde{u}_{\textnormal{a}}(k,M_{\textnormal{c}})^{2}\,,$
(37) $\displaystyle P_{\textnormal{a}}^{\textnormal{2h}}(k)$
$\displaystyle=\left[\frac{1}{F_{\textnormal{h}}\bar{\rho}_{\textnormal{a}}}\int_{M_{\textnormal{cut}}}^{\infty}\textnormal{d}M_{\textnormal{c}}n(M_{\textnormal{c}})b(M_{\textnormal{c}})M_{\textnormal{a}}(M_{\textnormal{c}})\times\right.$
$\displaystyle\hskip
14.22636pt\left.\tilde{u}_{\textnormal{a}}(k,M_{\textnormal{c}})\right]^{2}P_{\textnormal{c}}^{\textnormal{L}}(k)\,.$
(38)

Here we have introduced $M_{\textnormal{cut}}$: the cut-off mass below which
the axions can no longer cluster. The cut-off mass is also involved in the
clustered fraction $F_{\textnormal{h}}$ which is defined via Massara  _et al._
(2014)

$F_{\textnormal{h}}=\frac{1}{\bar{\rho}_{\textnormal{a}}}\int_{M_{\textnormal{cut}}}^{\infty}\textnormal{d}M_{\textnormal{c}}n(M_{\textnormal{c}})b(M_{\textnormal{c}})M_{\textnormal{a}}(M_{\textnormal{c}})\,.$
(39)

This means the three new quantities we have to specify to complete the MDM HM
are: the cut-off mass, $M_{\textnormal{cut}}$, the ULA halo mass relation, i.
e. the ULA halo mass as a function of the c-halo mass, and the ULA halo
density profile.

###  III.3 Cut-off Mass, Axion Halo Mass Relations and Axion Halo Density
Profile

Cut-off Mass To find an expression for the cut-off mass we will follow Ref.
Marsh and Silk (2014), where we proposed that in a pure ULA cosmology no halo
will form if the halo Jeans scale is larger than the virial radius. The halo
Jeans scale is similar to the linear Jeans scale, but depends on the halo
profile of the c-halo. The halo Jeans scale is given by Hu  _et al._ (2000):

$\displaystyle k_{\textnormal{hJ}}=$ $\displaystyle
66.5(1+z)^{-1/4}\left(\frac{\Omega_{\textnormal{m}}h^{2}}{0.12}\right)^{1/4}\left(\frac{m_{a}}{10^{-22}\,\textnormal{eV}}\right)^{1/2}$
$\displaystyle\left(\frac{\rho_{\textnormal{NFW}}(r_{\textnormal{hJ}})}{\bar{\rho}_{m}}\right)^{1/4}\,\textnormal{Mpc}^{-1}\,.$
(40)

Here $r_{\textnormal{hJ}}$ is the halo Jeans length calculated by converting
$r=\pi/k$. Note that here the radius is set to half the wavelength as in Ref.
Marsh and Silk (2014) instead of $r=2\pi/k$ as in Ref. Hu  _et al._ (2000).

Since we assume that $r_{\textnormal{hJ}}\leq r_{\textnormal{v}}$ the NFW
profile can be approximated as

$\rho_{\textnormal{NFW}}(r_{\textnormal{hJ}})\approx\frac{\bar{\rho}_{\mathrm{m}}\Delta_{\textnormal{v}}c^{2}r_{\textnormal{v}}}{3f(c)r_{\textnormal{hJ}}}\,$
(41)

with $f(x)$ the function as in Eq. (17). Inserting this approximation back
into Equation (III.3) and using $k_{\textnormal{hJ}}=\pi/r_{\textnormal{hJ}}$,
the halo Jeans length becomes

$\displaystyle r_{\textnormal{hJ}}\approx$ $\displaystyle
2.2a^{-1/3}\left(\frac{m_{a}}{10^{-22}\,\textnormal{eV}}\right)^{-2/3}\left(\frac{100f(c)}{f(10)c^{2}}\right)^{1/3}$
$\displaystyle\left(\frac{M}{10^{10}M_{\odot}}\right)^{-1/9}\left(\frac{\Omega_{\textnormal{m}}h^{2}}{0.12}\right)^{-2/9}\,\textnormal{kpc}\,.$
(42)

Fig. 2 shows the ratio of the halo Jeans length to the virial radius as a
function of the dark matter halo mass for different axion masses. The
horizontal dashed line in the figure indicates the point where the ratio of
the two length is equal to one and thus gives the cut-off mass
$M_{\textnormal{cut}}$ for the ULA halo.

Figure 2: The ratio of the halo Jeans scale $r_{\textnormal{hJ}}$ to the
virial radius $r_{\textnormal{v}}$ as a function of the c-halo mass for
different axion masses colour coded as indicated in the colourbar. The
horizontal dashed line shows the point for which the two quantities are equal
and determines the cut-off mass.

ULA Halo Mass Relation In this paper we assume, that for a ULA halo which is
located around a cold matter halo with mass $M_{\textnormal{c}}$ the mass is
given by the cosmic abundance, i. e. the ULA halo mass relation is
$M_{\textnormal{ax}}=\frac{\Omega_{\textnormal{ax}}}{\Omega_{\textnormal{c}}}M_{\textnormal{c}}$.
Since ULAs cannot cluster into halos on small scales, this relation is only
given for cold matter halos above the cut-off mass $M_{\textnormal{cut}}$ and
below $M_{\textnormal{cut}}$ the ULA halo mass is assumed to be zero. In Fig.
3 the halo mass relation for axion halos is plotted for axions in the mass
range of $10^{-28}\,\textnormal{eV}\leq m_{\textnormal{ax}}\leq 10^{-21}\,$eV
and an axion fraction of 0.1. In the next section we will see how the axion
halo mass relation helps to find the axion halo density profile.

Figure 3: The mass of the axion halo as a function of the corresponding c-halo
mass for axions in the mass range $10^{-28}\,\textnormal{eV}\leq
m_{\textnormal{ax}}\leq 10^{-21}\,$eV as colour coded in the colourbar and an
axion fraction of $f_{\textnormal{ax}}=0.1$. Below the cut-off mass for each
axion mass the axion halo mass is assumed to be zero.

Axion Halo Density Profile In a cosmology with ULAs we have no fitting
function of the axion density profile. There are simulations with mixture of
ULAs and cold dark matter, though, that tell us something about the shape of
the profile Schwabe  _et al._ (2020). However, in the case of a pure axion DM
cosmology a density profile is found by simulations and a fitting formula was
determined in Schive  _et al._ (2014a). The high resolution simulations showed
that the core of the axion density profile is given by a soliton whereas the
outer regions follow a CDM NFW-profile as in Eq. (16). Ref. Schive  _et al._
(2014a) found that the soliton in a pure ULA cosmology is well fitted by:

$\displaystyle\rho_{\textnormal{c}}(r)=$
$\displaystyle\frac{1.9(1+z)}{(1+9.1\times
10^{-2}(r/r_{\textnormal{c}})^{2})^{8}}\left(\frac{r_{\textnormal{c}}}{\textnormal{kpc}}\right)^{-4}$
$\displaystyle\left(\frac{m_{\textnormal{ax}}}{10^{-23}\,\textnormal{eV}}\right)^{-2}\,M_{\odot}\textnormal{pc}^{-3}\,,$
(43)

with $r_{\textnormal{c}}$ the core radius where the density drops to one half
of the central density. Further, Ref. Schive  _et al._ (2014b) determined this
core radius to be:

$\displaystyle r_{\textnormal{c}}=$ $\displaystyle
1.6(1+z)^{-1/2}\left(\frac{m_{\textnormal{ax}}}{10^{-22}\,\textnormal{eV}}\right)^{-1}\left(\frac{\Delta_{\textnormal{v}}(z)}{\Delta_{\textnormal{v}}(0)}\right)^{-1/6}$
$\displaystyle\left(\frac{M_{\textnormal{h}}}{10^{9}\,M_{\odot}}\right)^{-1/3}\,\textnormal{kpc}\,.$
(44)

Here $M_{\textnormal{h}}$ is the mass of the ULA halo.

Ref. Schwabe  _et al._ (2020) have simulated spherical collapse of halos in
mixed CDM-ULA models, which showed that also in this case the ULA halo density
profile is given by a soliton core and an NFW profile in the outer regions.
The soliton core forms only as long as $f_{\textnormal{ax}}\geq 0.1$. For
lower axion fractions the simulations showed that strong fluctuations in the
central density profile do not allow a fit to the soliton profile Schwabe  _et
al._ (2020). Therefore, we restricted the axion fraction to the range
$f_{\textnormal{ax}}\in[0.1,0.5]$ (the upper bound is given by the biased
tracer approach), and consider only halos that host soliton cores. We are
conservative with the lower bound, but new high resolution simulations of
mixed CDM-ULA model could show that we can remove the lower bound and that we
only have to set the upper bound of the axion fraction, i. e.
$f_{\textnormal{ax}}\leq 0.5$. The NFW profile for the axion halo in the outer
regions was found to be the same as the surrounding c-halo scaled by the
cosmic abundance $\Omega_{\textnormal{ax}}/\Omega_{\textnormal{c}}$. Since the
axion halos form in the potential wells of the c-halos the cold matter
influenced the ground state solution of the ULA and the solution given by the
soliton above has to be modified. We keep the same shape, Eq. (III.3), but
determine the core radius from the CDM virial velocity, and rescale the
central soliton density.
To construct a new soliton radius we used the radius defined by the
characteristic velocity, $v_{\textnormal{c}}$, which scales like the de
Broglie radius. Fits to simulations Mocz  _et al._ (2017a) give

$r_{\textnormal{c}}=\frac{2\pi}{7.5}\frac{\hbar}{m_{\textnormal{ax}}v_{\textnormal{c}}}\,,$
(45)

if the core is in equilibrium with its host halo. If we take the
characteristic velocity equal to the virial velocity Niemeyer (2020);
Eggemeier  _et al._ (2022),
$v_{\textnormal{v}}=(GM_{\textnormal{h}}/r_{\textnormal{v}})^{1/2}$, then the
soliton radius becomes:

$\displaystyle r_{\textnormal{c}}=$ $\displaystyle
1.2(1+z)^{-1/2}\left(\frac{m_{\textnormal{ax}}}{10^{-22}\,\textnormal{eV}}\right)^{-1}\left(\frac{\Delta_{\textnormal{v}}(z)}{\Delta_{\textnormal{v}}(0)}\right)^{-1/6}$
$\displaystyle\left(\frac{M_{\textnormal{h}}}{10^{9}\,M_{\odot}}\right)^{-1/6}\left(\frac{\Omega_{\textnormal{m}}h^{2}}{0.12}\right)^{-1/6}\,\textnormal{kpc}\,.$
(46)

In addition to the modified core radius we have to change the central density
of the soliton such that the ULA halo has the correct mass given by the ULA
halo mass relation described above. We thus rescale the soliton density with a
factor $A$, which is set by fixing:

$M_{\textnormal{ax}}=4\pi\int_{0}^{r_{\textnormal{v}}}\textnormal{d}r(\Theta(r_{\textnormal{i}}-r)A\rho_{\textnormal{c}}(r)+\Theta(r-r_{\textnormal{i}})\rho_{\textnormal{NFW}}(r))r^{2}\,,$
(47)

with $\Theta$ the Heaviside step function and $r_{\textnormal{i}}$ the radius
where the two profiles cross.

Figure 4: The axion halo density profiles for three different c-halos for
$m_{\textnormal{a}}=10^{-22}\,$eV with $f_{\textnormal{ax}}=0.1$ as solid
lines. The profile has a soliton core, see Eq. (III.3), with a soliton radius
as in Eq. (III.3) and an NFW profile on the outer region (see the text for the
exact construction). For comparison also the corresponding NFW profiles of the
c-halos are plotted in dashed lines.

With this at hand we can now determine the ULA halo density profile. By
computing the central density scaling, $A$, we found that Eq. (47) has no
solution if the axion halo mass is very close to the cut-off mass. This is
because the virial radius of a halo with a mass equal to the cut-off mass is
very similar to the soliton radius of Eq. (III.3) where the soliton density
falls rapidly with increasing radius. Therefore, we decide to set the new cut-
off mass to a little higher value where no solution is found for the central
density of the soliton profile. The axion halo profile is shown in Figure 4
for three different halo masses and an axion mass of
$m_{\textnormal{a}}=10^{-22}\,\textnormal{eV}$ with $f_{\textnormal{ax}}=0.1$.
The profiles show the soliton core and an NFW profile for large $r$. For
comparison also the NFW profiles of the corresponding c-halo are shown in
dashed lines. We can see that the axion profile is less massive in the core
and shows a flat core rather than a cusp like the NFW profile. Our constructed
profiles resemble closely the simulated profiles of Ref. Schwabe  _et al._
(2020).

##  IV Results

Figure 5: The non-linear power spectrum from the halo model in a MDM cosmology
with axions of $m_{\textnormal{ax}}=10^{-25}\,$eV and
$f_{\textnormal{ax}}=0.1$. The three terms of the power spectrum as in Eq.
(30) are shown as well as the resulting total non-linear matter power spectrum
in dashed black. For comparison also the total matter power spectrum in a
$\Lambda$CDM cosmology is represented as a black dotted line.

###  IV.1 Power Spectrum from MDM Halo Model

The non-linear power spectrum with the extended halo model described in the
previous Section is shown in Fig. 5 for a MDM cosmology with
$m_{\textnormal{ax}}=10^{-25}\,$eV and $f_{\textnormal{ax}}=0.1$. To
understand the influence of ULAs on the power spectrum in more detail, we
compare the MDM halo model with the $\Lambda$CDM halo model in Fig. 6 for a
ULA fraction of 0.1 and mass range $10^{-28}\,\textnormal{eV}\leq
m_{\textnormal{a}}\leq 10^{-21}\,$eV. As expected the non-linear power
spectrum shows a suppression on large wavenumbers compared to pure CDM,
asymptoting to a constant step-size. The size of the step is fixed by the
relative abundance, $\Omega_{\rm ax}/\Omega_{\rm d}$, and transitions with
$m_{\rm ax}$ as $k_{\rm J,eq}$ crosses through $1\,h\text{ Mpc}^{-1}$.

The suppression scale in the non-linear power can start at very different
wavenumbers compared to the linear theory, with the difference depending on
the ULA mass (a similar effect has also been seen before in Refs. Hložek  _et
al._ (2017); Marsh (2016b)). In linear theory, suppression relative to CDM
comes from the mass dependence of the axion Jeans scale. But why does the
suppression wavenumber in the non-linear power spectrum change? This can be
understood when we look at the formula for the non-linear power spectrum which
is given by a one and two halo term, see Eq. (11), and the transition between
these two terms is around $k_{\textnormal{t}}\sim 1\,h\textnormal{Mpc}^{-1}$.
Furthermore, the two halo term is proportional to the linear power spectrum
and thus, as long as the suppression of the linear power spectrum starts below
the transition wavenumber, i. e. $k_{\textnormal{J,eq}}<k_{\textnormal{t}}$,
the suppression of the non-linear power spectrum starts at the same scale as
the linear one. This occurs for masses $m_{\rm ax}\lesssim 10^{-23}\text{
eV}$.

For wavenumbers higher than the transition wavenumber the one halo term starts
to dominate and the non-linear power spectrum departs strongly from the linear
one. Hence, the location of the suppression no longer depends simply on the
linear Jeans scale. For higher mass ULAs with linear Jeans scale
$k_{\textnormal{J,eq}}>k_{\textnormal{t}}$, the difference to the pure CDM
case is driven by the one-halo term, which in turn is dominated by the cold
halos themselves. The power spectrum drops when halos contributing to it
become less massive than the cut-off mass. Additional suppression is driven
“passively” by the effect of ULAs in reducing the clustering of CDM on small
scales, which reduces the variance of CDM density fluctuations.

Figure 6: The ratio of the non-linear (solid lines) and linear (dashed lines)
power spectrum in a MDM cosmology with axions in the mass range
$10^{-28}\,\textnormal{eV}\leq m_{\textnormal{ax}}\leq 10^{-21}\,$eV and 10 %
axions to the $\Lambda$CDM case. The non-linear power spectra are calculated
with the halo model described in this work and the linear ones are calculated
with axionCAMB.

Another feature which can be seen in Fig. 6 is that the lowest masses
considered, $m_{\textnormal{ax}}\leq 10^{-27}\,$eV, the ratio has a spoon-like
shape. A similar shape was found if one compares the non-linear power spectrum
with massive neutrinos to the power spectrum of a $\Lambda$CDM cosmology
computed with the halo model Refs. Massara  _et al._ (2014); Hannestad  _et
al._ (2020), or also from simulation, see e. g. Refs. Brandbyge  _et al._
(2008); Bird  _et al._ (2012). The appearance of a similar feature for ULAs is
thus not unexpected.

The last feature we can see in Fig. 6 is an _enhancement_ in power for the MDM
model compared to pure CDM around $k\sim 1\,h\textnormal{Mpc}^{-1}$ for the
higher ULA masses. We can understand this when we look at the ratio of the
three different parts of the power spectrum in Eq. (11) to the $\Lambda$CDM
halo model, as shown in Fig. 7. In the ratios of the cross and axion parts
(bottom left and right panels) we see a strong enhancement at the scales
mentioned above and this comes from the shape of the ULA halo density profile
which is different from the cold one, i.e. this is caused by the coherence of
the soliton, which increases the correlation function of the ULA field on
small scales. The enhancement is not present for all ULA masses, since it
requires a conincidence between the one-to-two halo transition in $P(k)$, and
the size of the soliton in the halo mass dominating the power at this scale,
which occurs for $m_{\rm ax}\approx 10^{-22}\text{ eV}$. This prediction of
our model is in complete agreement with the simulations of Ref. Nori  _et al._
(2019), who observed a small increase in the power for $m_{\rm ax}\approx
10^{-22}\text{ eV}$ in pure ULA cosmologies only after accounting for the
effect of the “quantum pressure” terms in the effective fluid description of
the Schrödinger-Poisson equation.

In Fig. 6 we have shown for illustration the relative power over a wide $k$
range, and in particular for very large wavenumbers. This means, however, that
the halo model is evaluated up to very small lengths where internal properties
of individual halos have to be taken into account, e. g. baryonic feedback
from star formation or active galactic nuclei Freundlich  _et al._ (2016).
Thus our model is idealised on these scales, and should not be considered
realistic. See Refs. Mead  _et al._ (2015b, 2021a); Massara  _et al._ (2014)
for more discussion. We expect our model to be relatively accurate up to
approximately $k=10\,h\textnormal{Mpc}^{-1}$.

Figure 7: The ratio of the halo model in a MDM cosmology with axions (top
left) and all three terms as in Eq. (11) (top right: cold part, bottom left:
cross part and bottom right: axion part) to the $\Lambda$CDM halo model. Here
the axions have a mass of $10^{-28}\,\textnormal{eV}\leq
m_{\textnormal{ax}}\leq 10^{-21}\,$eV and 10 % of the DM are axions.

###  IV.2 DE-like Axions

Our halo model should work extremely well for dark energy (DE) like ULAs, as
defined by Ref. Hlozek  _et al._ (2015a) with $m_{\rm ax}\leq 10^{-28}\text{
eV}$, where no simulations are available at the moment. As mentioned above in
the discussion of Fig. 6 a spoon like shape is seen for very light ULAs, i. e.
DE-like axions. We want investigate this feature further comparing to the
$\Lambda$CDM power spectrum for different ULA fractions, i. e. $0.05\leq
f_{\mathrm{ax}}\leq 0.25$, and an axion mass of
$m_{\textnormal{ax}}=10^{-28}\,$eV, shown in Fig. 8. Here the we use our mixed
Halo model for axion fractions below the discussed lower bound of
$f_{\mathrm{ax}}$ in Sec. III.3. But at this low fractions and small ULA mass
the effect of the axions on the non-linear power spectrum is very small and
thus we decided to extend the mixed halo model for lower ULA mass also to
smaller axion fractions. We observe that the spoon like shape is more dominant
for smaller axion fraction and faded away when the fraction is raised.

For lower masses still, $m_{\rm ax}<10^{-28}\text{ eV}$, we see from Fig. 3
that we do not expect such ULAs to reside in any cosmologically known halos
$M_{h}\lesssim 10^{15}M_{\odot}$. This justifies the approximation taken in
Ref. Bauer  _et al._ (2020) to remove ULAs entirely from the HM at low masses,
and has the same effect as the removal of neutrinos from halos in the case of
HMCode (the effect of this approximation compared to the full mixed halo model
of neutrinos is discussed in Appendix B). This suggests that for $m_{\rm
ax}<10^{-28}\text{ eV}$ one can leverage the accuracy of HMCode for ULAs at
any density fraction allowed by current constraints from linear scales Hlozek
_et al._ (2018), although Lagrangian perturbation theory Laguë  _et al._
(2022) will also be accurate in this regime. This is due to the fact that the
axion perturbations $\delta_{\rm ax}$ have a scale-dependent growth which is
suppressed on small scales. In the case of DE-like axions, the perturbations
on scales $k\lesssim 0.1\,h\text{ Mpc}^{-1}$ will not grow until the present
day and remain in the linear regime.

Figure 8: The ratio of the non-linear (solid lines) and linear (dashed lines)
power spectrum in a MDM cosmology with axions of a mass $10^{-28}\,$eV and for
an axion fraction in the range $0.05\leq f_{\mathrm{ax}}\leq 0.25$ to the
$\Lambda$CDM case.

##  V Discussion and Outlook

We presented in this paper an improved halo model for a cosmology composed of
CDM and ULAs. The standard pure CDM halo model assumes that all matter is
bound in halos and that the two-point correlation between the matter is given
by the correlations inside one halo and between two separate halos. However,
in a mixed dark matter cosmology with a sub-dominant ULA component, because of
the ULA Jeans scale we can no longer assume that all matter is bound into
halos. More generally, due to their different clustering properties CDM and
ULAs must be treated differently. In the power spectrum we have a cold-cold
part, a cross correlation between the cold and the non-cold matter, and an ULA
self-clustering part. The ULA power spectrum has a “smooth” component from the
matter that is not contained in halos and a component that can cluster inside
halos on large enough scales

Our model accounts for all of these effects. Our model assumes (in a manner
that is consistent with observations) that ULAs make up only a sub-dominant
component of the total matter, and thus we treat ULAs as a biased tracer of
the cold matter. In the spirit of the biased tracer model (which has been
applied successfully to neutrinos and neutral hydrogen), we have proposed a
density profile for ULAs inside halos, as well as a relationship between the
ULA and CDM halo masses, including a cut-off mass below which ULAs do not
cluster inside halos. For the ULA density profile in a mixed halo, no fitting
formulae are available in the literature, and we proposed a model based on
Ref. Schive  _et al._ (2014a) fitting formula in a pure ULA cosmology, and
observations from simulations in mixed ULA-CDM spherical collapse of Ref.
Schwabe  _et al._ (2020).

For ULAs in the mass range $10^{-28}\text{ eV}\lesssim m_{\rm ax}\lesssim
10^{-23}\text{ eV}$ our model predicts that the suppression scale between the
MDM power spectrum and pure CDM occurs near the linear theory Jeans scale,
with an additional spoon-like suppression similar to the case of neutrinos,
before asymptoting to a constant fixed by the relative DM density fractions.
At larger ULA masses, the suppression scale is controlled by the one-halo
term, and moves out to larger wavenumbers. For ULAs with $m_{\rm ax}\approx
10^{-22}\text{ eV}$ we observe a small enhancement in the power relative to
CDM on intermediate scales, which we attribute to the role of solitons in the
power spectrum. This prediction of our model is in qualitative agreement with
the simulations of Ref. Nori  _et al._ (2019). Finally, we also considered DE-
like ULAs with $m_{\rm ax}\lesssim 10^{-28}\text{ eV}$ and also found a spoon-
like feature in the power spectrum.

Our model is inspired by the mixed ULA-CDM simulations of Ref. Schwabe  _et
al._ (2020), who propose the density profile we use based on the Schroödinger-
Poisson equation, and also many simulations that observe a minimum halo mass
in pure ULA cosmologies in N-body (e.g. Ref. Schive  _et al._ (2016)).
However, we have not been able to calibrate our model on cosmological mixed
ULA-CDM simulations. Some such simulations are in preparation Lague  _et al._
(2022). However, the box size of these simulations is too small to have a
large number of halos. Thus we cannot test the cut-off mass in our model. The
relative power spectrum is in principle well predicted in relatively small
boxes Corasaniti  _et al._ (2017), however this is only true on scales where
there are many halos contributing to the power such that the one-halo term is
well sampled. We have also found that mixed DM simulations are limited in
resolution of the density profile Lague  _et al._ (2022). Nonetheless, this
shows us the way forward to calibrating our halo model in future.

An interesting extension of our work would be to include multiple ULA sub-
components, as one might expect in an “axiverse” Arvanitaki  _et al._ (2010).
The principles of the biased tracer model could be adopted for each component,
and one might expect e.g. “nested solitons” in the density profile Luu  _et
al._ (2020). Indeed, mixed ULA simulations were recently reported Gosenca  _et
al._ (2023). The linear transfer function in such a model is expected to
demonstrate multiple step features, and one would expect this to be reflected
in our model also in the halo mass function and non-linear power. An accurate
linear transfer function has, however, not been computed in such a case.
Extending e.g. axionCAMB to multiple fields would be evidently worthwhile.

We expect our halo model to be extremely useful in future analyses of
cosmological data. We demonstrated recently that the halo model for pure ULAs
can be used in analysis of Dark Energy Survey data Dentler  _et al._ (2022).
Ref. Dentler  _et al._ (2022) was only able to constrain the ULA mass, and
could not vary the fraction at low masses due to the lack of an appropriate
halo model. The model presented here fulfils that purpose and will allow for a
combined CMB+DES analysis covering all masses and fractions across the ULA
parameter space. Such an analysis will plug an important gap in current ULA
constraints between $10^{-25}\text{ eV}$ and $10^{-23}\text{ eV}$. This has
the potential to probe new parameter space of string theory models Arvanitaki
_et al._ (2010); Demirtas  _et al._ (2021); Mehta  _et al._ (2021), some of
which now make specific predictions in this region Cicoli  _et al._ (2022),
and has wider implications for the understanding of DM at the low-mass
frontier Grin  _et al._ (2019); Alves Batista  _et al._ (2021). Our model will
continue to be useful for constraining ULAs with next generation cosmological
data such as Simons Observatory Ade  _et al._ (2019), CMB-S4 Dvorkin  _et al._
(2022), and _Euclid_ Amendola  _et al._ (2018).

###### Acknowledgements.

We acknowledge useful discussions with Mona Dentler, Jens Niemeyer, and Bodo
Schwabe. We acknowledge Jurek Bauer for collaboration in an early stage of
this project. DJEM is supported by an Ernest Rutherford Fellowship from the
Science and Technologies Facilities Council (UK). This work made use of the
open-source libraries matplotlib Hunter (2007), numpy Harris  _et al._ (2020),
scipy Virtanen  _et al._ (2020), astropy Astropy Collaboration  _et al._
(2022).

##  Appendix A HMcode Parameters and Comparison

The halo model is a very good and quick model to find the non-linear power
spectrum and is used in a lot of different codes. One of the most frequently
used codes. the HMcode, is provided by Ref. Mead  _et al._ (2021b), with
previous versions Refs. Mead  _et al._ (2015a, 2016). HMCode introduces a
number of parameters to improve the model in its fit to simulations over the
standard HM. The parameters were fit to simulations of Ref. Heitmann  _et al._
(2013) such that the model accurately matches these simulations as well as the
simulations from Ref. Heitmann  _et al._ (2016). The newest code HMCode shows
excellent agreement to simulations for $k\lesssim 10\,h\textnormal{Mpc}^{-1}$
and $z<2$ with a root mean square of at most 2.5% Mead  _et al._ (2021b).

Since HMCode is a frequently used halo model code and shows excellent results,
we decide to implement the parameter in our halo model code, called
axionHMcode. When cosmological mixed DM ULA simulations become available we
can compare them with our MDM halo model with the parameters from HMCode to
see if these parameters also improve the HM with ULAs. In total there are six
new parameter and we will discuss them in this Appendix.

Halo bloating term: The NFW profile from Eq. (16) is modified in HMcode by the
parameter $\eta$ such that

$\rho_{\textnormal{NFW}}(r,M,\eta)=\frac{\rho_{\mathrm{char}}}{r/(\nu^{\eta}r_{s})\left(1+r/(\nu^{\eta}r_{s})\right)^{2}},\qquad
r\leq r_{\textnormal{v}}\nu^{\eta}\,,$ (48)

with

$\eta(z)=0.1281\sigma_{8,\,\textnormal{cc}}^{-0.3644}(z)>0\,,$ (49)

where $\sigma_{8,\,\textnormal{cc}}$ refers to the variance of the cold matter
linear power spectrum for $R=8\,h^{-1}\,\mathrm{Mpc}$ in Eq. (19). It can be
shown that the halo bloating term influences the Fourier transformation of the
NFW profile, Eq. (18), by scaling the wavenumber $k$ in the following way
Copeland  _et al._ (2020):

$\tilde{u}(k,M,z)\to\tilde{u}(\nu^{\eta}k,M,z)\,.$ (50)

Since $\eta>0$ the $k$-space halo density profile is pushed to higher $k$’s
for small halos and the profile for massive halos is transformed to smaller
wavenumbers. The exact shape of $\tilde{u}(k,M,z)$ is very important for the
one halo term and smaller halos play an increasingly large role for high
$k$’s. This means that at a scale around $k\sim 1\,h\textnormal{Mpc}^{-1}$,
where the one halo term becomes the dominant component, the halo bloating term
shows an effect on the non-linear power spectrum. Since on larger scales the
contribution from halos with high masses dominates we expect a suppression in
the power spectrum, because $\nu^{\eta}<1$. Then for higher and higher $k$’s
the shapes of smaller halos are more important and the $\eta$ parameter will
give an enhancement in the non-linear power spectrum ($\nu^{\eta}>1$). This
described effect can be seen in Fig. 9, where the ratio between the halo model
to the halo bloating term and the standard model is shown in blue.

One halo term damping: In the standard halo model approach the one halo term
is constant on large scales. However, this is not the correct behaviour due to
mass and momentum conservation. It was shown in Smith  _et al._ (2003) that
the one halo term should grow like $P^{1\textnormal{h}}(k)\propto k^{4}$ at
small $k$ (i.e. damp away compared to constant at small $k$). In HMCode this
is implemented trough the modification

$P^{1\textnormal{h}}(k)\to
P^{1\textnormal{h}}(k)\frac{(k/k_{*})^{4}}{1+(k/k_{*})^{4}}\,.$ (51)

The modified one halo term grows as expected and is suppressed on large
scales. This also ensures that on large scales the non-linear power spectrum
is given by the two halo term and hence is equal to the linear power spectrum.
The suppression depends on the free parameter $k_{*}$ which is fitted from
simulations to be Mead  _et al._ (2021b)

$k_{*}=0.05618\sigma_{8,\,\textnormal{cc}}^{-1.013}(z)\,h\textnormal{Mpc}^{-1}\,.$
(52)

Since the one halo term is suppressed on large scales, the total non-linear
power is also suppressed on large scales if the one halo term is modified as
above. This can be seen in Fig. 9 in orange.

Two halo term damping: Like the one halo term also the two halo term is damped
on some scales. These lengths are scales larger than a particular wavenumber
$k_{d}$. The damping takes the form

$P^{2\textnormal{h}}(k)\to
P^{2\textnormal{h}}(k)\left(1-f\frac{(k/k_{\textnormal{d}})^{n_{\textnormal{d}}}}{1+(k/k_{\textnormal{d}})^{n_{\textnormal{d}}}}\right)\,,$
(53)

where $n_{\textnormal{d}},\ k_{\textnormal{d}}$ and $f$ are fitting
parameters. The three parameters in Eq. (53) are fitted and given by Mead  _et
al._ (2021b)

$\displaystyle k_{\textnormal{d}}$
$\displaystyle=0.05699\sigma_{8,\,\textnormal{cc}}^{-1.089}(z)\,h\textnormal{Mpc}^{-1}\,,$
(54) $\displaystyle f$
$\displaystyle=0.2696\sigma_{8,\,\textnormal{cc}}^{0.9403}(z)\,,$ (55)
$\displaystyle n_{\textnormal{d}}$ $\displaystyle=2.853\,.$ (56)

We see that the damping power $n_{\textnormal{d}}>0$ and hence there is a
damping for $k\gg k_{\textnormal{d}}$ as long as the two halo term dominates
in the non-linear power spectrum. The explained behaviour can be seen in the
Fig. 9 in green.

Smoothing parameter: In the standard halo model the power spectrum is the sum
of the one halo and two halo term, see Eq. (11). However, if the two halo term
is of comparable size to the one halo term (transition region), the assumption
of a purely additive behaviour is too simple. The scale of the transition
region is also known as the quasi-linear regime and Ref. Mead  _et al._
(2015a) found that the halo model performed quite poorly at the quasi-linear
regime. Hence the HMcode introduces a transition parameter by modelling:

$P(k)=(P^{\mathrm{1h}}(k)^{\alpha}+P^{\mathrm{2h}}(k)^{\alpha})^{1/\alpha}\,,$
(57)

where $\alpha$ is the parameter that shapes the transition. If $\alpha<1$ the
transition is smoothed whereas for $\alpha>1$ the transition is sharper. The
general form of the smoothing parameter is assumed to be

$\alpha=ab^{n_{\textnormal{eff}}(z)}\,.$ (58)

Here $a$ and $b$ are the fitting parameters and $n_{\textnormal{eff}}(z)$ is
the effective spectral index at the non-linear length $R_{\textnormal{nl}}$,
where
$\sigma_{\textnormal{cc}}(R_{\textnormal{nl}},z)=\delta_{\textnormal{c}}$:

$n_{\textnormal{eff}}(z)=-\frac{\textnormal{d}\,\textnormal{ln}\sigma_{\textnormal{cc}}^{2}(R,z)}{\textnormal{d}\,\textnormal{ln}R}\bigg{|}_{\sigma=\delta_{\rm
crit}}-3\,.$ (59)

In the newest version HMCode the smoothing parameter is fitted to be Mead  _et
al._ (2021b)

$\alpha=1.876\times 1.603^{n_{\textnormal{eff}}(z)}\,.$ (60)

The effect of the smoothing parameter can be seen in Fig. 9 in red and a clear
enhancement around the transition region is visible since $\alpha<1$ in the
analysed cosmology. Note that the halo model with the smoothing parameter is
not equal to the standard halo model on large scales because only the
smoothing parameter is used and the one halo term is not damped in this plot
and thus the constant one halo term influences the power spectrum on large
scales.

In Figure 9 also the ratio with all parameters of the HMCode is shown by a
solid black line and the ratio to the linear power spectrum is shown in dotted
grey. On large scales the non-linear power spectrum with all parameters
described above and the linear power spectrum coincide and the same is true if
only the one halo damping is used. Therefore, we use only the one halo damping
in axionHMcode as default to ensure the correct behaviour on large scales. The
other parameters are optional in our halo model, and free parameters in the
code.

Figure 9: Ratio between the halo model with one of the parameters in HMCode
and the standard halo model. In the ratio with the two halo damping three
parameters are involved (see the text for more details). Each parameter shows
a different effect on the non-linear power spectrum. The black line shows the
ratio when all parameters and a non-$\Lambda$CDM approach for the
overdensities are used.

##  Appendix B Massive Neutrinos

The full HM for massive neutrinos was developed similarly to our model in Ref.
Massara  _et al._ (2014), and calibrated to a small range of simulations.
Crucially, the simulations in this case include the distinct dynamics of
massive neutrinos compared to CDM, and allow them to non-linearly cluster
around halos.

In contrast, the MiraTitan Heitmann  _et al._ (2016) emulator simulates only
the _linear_ evolution of neutrinos and does not allow them to cluster inside
halos. Thus, the massive neutrino model in HMCode, which is calibrated to
MiraTitan, takes an approximate treatment of massive neutrinos that removes
them from halos. The accuracy of this treatment, as applied to MiraTitan, is
discussed in depth in Ref. Mead  _et al._ (2021b). We developed our own
implementation of the massive neutrino halo model of Ref. Massara  _et al._
(2014) 222Available at https://github.com/SophieMLV/nuHMcode. Comparing this
treatment to the HMCode approximation can be used as an estimate of the
theoretical error in HMCode, and consequently in MiraTitan, of neglecting full
neutrino dynamics. We briefly present these results here. A more complete
comparison, using e.g. larger neutrino simulations such as Ref. Emberson  _et
al._ (2017), would be interesting, but is beyond the scope of the present
work.

Neutrino Halo model ingredients: The halo model with massive neutrinos uses
the same formulae as the model with ULAs described in Section III.2. The
difference is that the massive neutrino halos have a different shape, halo
mass relation and the cut-off mass. The advantage for the model with massive
neutrinos is that we have fitting functions for the neutrino halo ($\nu$-halo)
profiles for $\sum m_{\nu}=0.3\,$eV and $\sum m_{\nu}=0.6\,$eV from
simulations in Ref. Villaescusa-Navarro  _et al._ (2013). The authors found
the following fitting function

$\rho_{\nu}(r,M_{\textnormal{c}})=\frac{\delta_{\textnormal{core}}}{1+\left(r/r_{\textnormal{core}}\right)^{\alpha}}\bar{\rho}_{\nu}\,.$
(61)

Here $\delta_{\textnormal{core}}$, $r_{\textnormal{core}}$ and $\alpha$ are
fitting parameters and depend on the mass of the corresponding c-halo
$M_{\textnormal{c}}$.

However, the resolution of the N-body simulations in Ref. Massara  _et al._
(2014) was not high enough to resolve the core of the neutrino halos for
c-halos of mass below $M_{\textnormal{c}}\sim 10^{13.5}\,h^{-1}\,M_{\odot}$.
Hence, the profile in this case is chosen to behave like

$\rho_{\nu}(r)=\frac{\kappa}{r^{\beta}}\bar{\rho}_{\nu}\,$ (62)

and to reproduce the outer region of the neutrino density profile as in Eq.
(61). Here $\kappa$ and $\beta$ are again determine by fitting to the
simulation results. For a total neutrino masses of $\sum m_{\nu}=0.3\,$eV and
$\sum m_{\nu}=0.6\,$eV the parameters are given in Table 2.

parameter |  $\sum m_{\nu}$ = 0.3 eV |  $\sum m_{\nu}$ = 0.6 eV
---|---|---
$\delta_{\textnormal{core}}$ | $6.056\times 10^{-8}M_{\textnormal{c}}^{0.58}$ | $3.7478\times 10^{-8}M_{\textnormal{c}}^{0.64}$
$r_{\textnormal{core}}$ [$h^{-1}$kpc] | $4.029\times 10^{-8}M_{\textnormal{c}}^{0.68}$ | $2.046\times 10^{-4}M_{\textnormal{c}}^{0.43}$
$\alpha$ | $-6.69+0.24\log(M_{\textnormal{c}})$ | $-4.62+0.19\log(M_{\textnormal{c}})$
$\beta$ | $-2.06+0.09\log(M_{\textnormal{c}})$ | $-3.64+0.15\log(M_{\textnormal{c}})$
$\kappa$ [$(h^{-1}\textnormal{kpc})^{-\beta}$] | $0.19+3.242\times 10^{-19}M_{\textnormal{c}}^{1.5}$ | $0.24+1.144\times 10^{-20}M_{\textnormal{c}}^{1.7}$
Table 2: Fitting functions for the parameters $\delta_{\textnormal{core}}$,
$r_{\textnormal{core}}$, $\alpha$, $\beta$ and $\kappa$ in Eq. (61) and Eq.
(62) as in figure 10 of Villaescusa-Navarro  _et al._ (2013).

Different approaches can be made to find a cut-off mass. We continue to follow
Ref. Massara  _et al._ (2014) who defined $M_{\textnormal{cut}}$ as the c-halo
for which the corresponding $\nu$-halo has a mass of at least 10 % of the
background neutrino density enclosed in the same radius:

$M_{\nu}(M_{\textnormal{cut}})=0.1\frac{4\pi\bar{\rho}_{\nu}}{3}r_{\textnormal{v}}^{3}(M_{\textnormal{cut}})\,.$
(63)

with the cut-off mass, $M_{\textnormal{cut}}$, the $\nu$-halo density profile,
$\rho_{\nu}$, and the halo mass relation for the $\nu$-halo,
$M_{\nu}(M_{\textnormal{c}})$. Fig. 10 shows the mass of the $\nu$-halo as a
function of the corresponding cold matter halo. In the figure we also see the
transition between the two neutrino halo profiles around
$M_{\textnormal{c}}=10^{13.5}\,h^{-1}\,M_{\odot}$ as a small discontinuity in
$M_{\nu}(M_{\textnormal{c}})$.

Figure 10: The mass of the neutrino halo as a function of the corresponding
c-halo for a total neutrino mass of $\sum m_{\nu}=0.3\,$eV in blue and for
$\sum m_{\nu}=0.6\,$eV in orange. At the cut-off mass the halo mass relation
drops to zero, i. e. there is no $\nu$-halo for a c-halo with mass below the
cut-off mass. The jump in the halo mass relation around
$M_{\textnormal{c}}=10^{13.5}\,h^{-1}\,M_{\odot}$ comes from the change in the
density profile, Eq. (61) and (62). Note that these lines do not have a slope
equal to the cosmic abundance $\Omega_{\nu}/\Omega_{\textnormal{c}}$.

Comparison to HMCode: We want to compare the results from our halo model code
with massive neutrinos, $\nu$HMCode, with the HMCode predictions. The massive
neutrinos are implemented such that they evolve linearly and do not cluster in
halos at all. This means Ref. Mead  _et al._ (2021b) used the normal
$\Lambda$CDM halo model as in Section III, but removed the massive neutrino
from the halo density profile, Eq. (18), by transforming the profile as

$\tilde{u}(k,M)\to(1-f_{\nu})\tilde{u}(k,M)\,.$ (64)

The modifications of the halo model for HMCode in Appendix A can also be
adopted to our $\nu$HMCode. Note that for the smoothing parameter $\alpha$ the
smoothing is only applied to the cold matter part of the fully extended halo
model. We decided to smooth only the cold part because the expression for the
cold part is given by the standard halo model, Eq. (11), by using the cold
matter quantities rather than the total matter terms and thus the expression
is very similar to the one where the smoothing parameter is used in HMcode.

So, we compare here the following models first our full $\nu$HMcode without
the parameters with the HMcode massive neutrino approximation also with no
parameters and second the $\nu$HMcode with all parameters with the HMCode
approximation also with all parameters. The difference between these models
can be understood as the effect of the clustered treatment of the massive
neutrinos on the non-linear power spectrum. The ratios of the models for both
massive neutrino masses can be seen in Figure 11 and we see that the full
treatment in $\nu$HMcode has an enhancement for large wavenumbers in both
configurations and for both masses. The extra power comes from the additional
clustering of massive neutrinos inside halos, which were taken into account in
our halo model. The difference for a total neutrino mass $\sum m_{\nu}=0.3$ eV
is not larger than 1 % and for $\sum m_{\nu}=0.3$ eV the discrepancy never
reaches 3 % for wavenumbers below $k\sim 10\,h\textnormal{Mpc}^{-1}$. Thus the
advantage of using the full treatment of massive neutrinos is very small and
the question is whether it is worthwhile to use the more computationally
intensive but more accurate model or to work with the simplified model which
achieves comparable results in much less time. Moreover, the difference of the
two models as shown in Fig. 11 can be used for an approximation of the error
of simulations which treat massive neutrinos only linearly.

Figure 11: The ratio of the $\nu$HMcode without the HMCode parameters (dashed
lines) and with the HMCode parameters (dashed-dotted lines) to the
corresponding approximation of the power spectrum with massive neutrinos as in
HMCode for $\sum m_{\nu}=0.3$ eV in green and $\sum m_{\nu}=0.6$ eV in red.

##  Appendix C Convergence Checks

In Section III we mentioned that we have to check if the integral involved in
the two halo term, see Eq. (13), converges. In theory this is an improper
integral, but in the numerical computation the integral has to be evaluated on
a finite interval $[M_{\textnormal{min}},M_{\textnormal{max}}]$. With the
finite interval we can no longer ensure that to $k\to 0$ the integral
converges to

$\frac{1}{\bar{\rho}}\int_{0}^{\infty}\mathrm{d}MMn(M)b(M)=1\,,$ (65)

which implies that on large scales the non-linear power spectrum is equal to
the linear one (this makes sense because we live in an isotropic and
homogeneous universe). So we have to be careful when choosing the boundaries
for the integral of the two halo term. This problem was discussed in the
Appendix A of Ref. Mead  _et al._ (2020) and the authors showed that a
relatively simple modification of the two halo term integral solves this
problem.

Note that the problem of the finite interval is not the upper bound, since the
HMF $n(M)$ has an exponential cut-off approximately around
$10^{16}M_{\odot}\,h^{-1}$ and thus if we take $M_{\textnormal{max}}$
sufficiently large it can be understood as infinite. Instead, the lower bound
is a problem due to the constantly rising HMF which states that a large amount
of the matter is contained in low mass halos Mead  _et al._ (2020) and thus
the correction of the two halo term integral depends on the lower bound
$M_{\textnormal{min}}$.

So the solution of Ref. Mead  _et al._ (2020) is defining a new function

$A(M_{\textnormal{min}})=1-\frac{1}{\bar{\rho}}\int_{M_{\textnormal{min}}}^{\infty}\mathrm{d}MMn(M)b(M)\,,$
(66)

which gives the missing part of the integral below $M_{\textnormal{min}}$.
With this at hand we can define a new halo mass function which takes care of
the missing part from halos with masses below $M_{\textnormal{min}}$ Schmidt
(2016). The new HMF is given by the transformation

$n(M)\to
n(M)+\frac{A(M_{\textnormal{min}})\delta(M-M_{\textnormal{min}})}{b(M_{\textnormal{min}})M_{\textnormal{min}}/\bar{\rho}}\,,$
(67)

where $\delta$ is the Dirac delta distribution. This new HMF gives a
correction to the two halo term integral:

$\frac{1}{\bar{\rho}}\int_{M_{\textnormal{min}}}^{\infty}\mathrm{d}MMn(M)b(M)|\tilde{u}(k,M)|+\frac{A(M_{\textnormal{min}})\tilde{u}(k,M_{\textnormal{min}})}{M_{\textnormal{min}}/\bar{\rho}}\,.$
(68)

So as $k\to 0$ the first term goes to
$1/\bar{\rho}\int_{M_{\textnormal{min}}}^{\infty}\mathrm{d}MMn(M)$ and the
second term $\tilde{u}(k,M_{\textnormal{min}})\to
M_{\textnormal{min}}/\bar{\rho}$. Using Eq. (66) we see that the new two halo
term integral goes to one. So, we always use this new two halo term integral
when we compute the two halo term in our code.

For the one halo term we do not need this correction since in the one halo
term low mass halos contribute very little Mead  _et al._ (2020).

## References

  * Aghanim  _et al._ (2020a) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A1 (2020a), arXiv:1807.06205 [astro-ph.CO] .
  * Aghanim  _et al._ (2020b) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020b), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Alves Batista  _et al._ (2021) R. Alves Batista _et al._ , (2021), arXiv:2110.10074 [astro-ph.HE] .
  * Hlozek  _et al._ (2015a) R. Hlozek, D. Grin, D. J. E. Marsh, and P. G. Ferreira, Phys. Rev. D 91, 103512 (2015a), arXiv:1410.2896 [astro-ph.CO] .
  * Hlozek  _et al._ (2018) R. Hlozek, D. J. E. Marsh, and D. Grin, Mon. Not. Roy. Astron. Soc. 476, 3063 (2018), arXiv:1708.05681 [astro-ph.CO] .
  * Kobayashi  _et al._ (2017) T. Kobayashi, R. Murgia, A. De Simone, V. Iršič, and M. Viel, Phys. Rev. D 96, 123514 (2017), arXiv:1708.00015 [astro-ph.CO] .
  * Laguë  _et al._ (2022) A. Laguë, J. R. Bond, R. Hložek, K. K. Rogers, D. J. E. Marsh, and D. Grin, JCAP 01, 049 (2022), arXiv:2104.07802 [astro-ph.CO] .
  * Witten (1984) E. Witten, Phys. Lett. B 149, 351 (1984).
  * Svrcek and Witten (2006) P. Svrcek and E. Witten, JHEP 06, 051 (2006), arXiv:hep-th/0605206 .
  * Conlon (2006) J. P. Conlon, JHEP 05, 078 (2006), arXiv:hep-th/0602233 .
  * Arvanitaki  _et al._ (2010) A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, Phys. Rev. D 81, 123530 (2010), arXiv:0905.4720 [hep-th] .
  * Mehta  _et al._ (2021) V. M. Mehta, M. Demirtas, C. Long, D. J. E. Marsh, L. McAllister, and M. J. Stott, JCAP 07, 033 (2021), arXiv:2103.06812 [hep-th] .
  * Demirtas  _et al._ (2021) M. Demirtas, N. Gendler, C. Long, L. McAllister, and J. Moritz, (2021), arXiv:2112.04503 [hep-th] .
  * Cicoli  _et al._ (2022) M. Cicoli, V. Guidetti, N. Righi, and A. Westphal, JHEP 05, 107 (2022), arXiv:2110.02964 [hep-th] .
  * Hložek  _et al._ (2017) R. Hložek, D. J. E. Marsh, D. Grin, R. Allison, J. Dunkley, and E. Calabrese, Phys. Rev. D 95, 123511 (2017), arXiv:1607.08208 [astro-ph.CO] .
  * Bauer  _et al._ (2020) J. B. Bauer, D. J. E. Marsh, R. Hložek, H. Padmanabhan, and A. Laguë, Monthly Notices of the Royal Astronomical Society 500, 3162 (2020), arXiv:2003.09655.
  * Farren  _et al._ (2022) G. S. Farren, D. Grin, A. H. Jaffe, R. Hložek, and D. J. E. Marsh, Phys. Rev. D 105, 063513 (2022), arXiv:2109.13268 [astro-ph.CO] .
  * Abazajian  _et al._ (2016) K. N. Abazajian _et al._ (CMB-S4), (2016), arXiv:1610.02743 [astro-ph.CO] .
  * Dvorkin  _et al._ (2022) C. Dvorkin _et al._ , in  _2022 Snowmass Summer Study_ (2022) arXiv:2203.07064 [hep-ph] .
  * Flitter and Kovetz (2022) J. Flitter and E. D. Kovetz, (2022), arXiv:2207.05083 [astro-ph.CO] .
  * Vogelsberger  _et al._ (2020) M. Vogelsberger, F. Marinacci, P. Torrey, and E. Puchwein, Nature Rev. Phys. 2, 42 (2020), arXiv:1909.07976 [astro-ph.GA] .
  * Heitmann  _et al._ (2013) K. Heitmann, E. Lawrence, J. Kwan, S. Habib, and D. Higdon, The Astrophysical Journal 780, 111 (2013).
  * Heitmann  _et al._ (2016) K. Heitmann, D. Bingham, E. Lawrence, S. Bergner, S. Habib, D. Higdon, A. Pope, R. Biswas, H. Finkel, N. Frontiere, and S. Bhattacharya, The Astrophysical Journal 820, 108 (2016).
  * Rogers  _et al._ (2019) K. K. Rogers, H. V. Peiris, A. Pontzen, S. Bird, L. Verde, and A. Font-Ribera, JCAP 02, 031 (2019), arXiv:1812.04631 [astro-ph.CO] .
  * Pedersen  _et al._ (2021) C. Pedersen, A. Font-Ribera, K. K. Rogers, P. McDonald, H. V. Peiris, A. Pontzen, and A. Slosar, JCAP 05, 033 (2021), arXiv:2011.15127 [astro-ph.CO] .
  * Cooray and Sheth (2002) A. Cooray and R. Sheth, Physics Reports 372, 1 (2002), arXiv: astro-ph/0206508.
  * Mead  _et al._ (2015a) A. J. Mead, J. A. Peacock, C. Heymans, S. Joudaki, and A. F. Heavens, Monthly Notices of the Royal Astronomical Society 454, 1958 (2015a).
  * Mead  _et al._ (2021a) A. Mead, S. Brieden, T. Tröster, and C. Heymans, Monthly Notices of the Royal Astronomical Society 502, 1401 (2021a), arXiv: 2009.01858.
  * Castro  _et al._ (2022) T. Castro _et al._ (Euclid), (2022), arXiv:2208.02174 [astro-ph.CO] .
  * Marsh and Silk (2014) D. J. E. Marsh and J. Silk, Monthly Notices of the Royal Astronomical Society 437, 2652 (2014), arXiv:1307.1705.
  * LoVerde and Zaldarriaga (2014) M. LoVerde and M. Zaldarriaga, Phys. Rev. D 89, 063502 (2014), arXiv:1310.6459 [astro-ph.CO] .
  * Massara  _et al._ (2014) E. Massara, F. Villaescusa-Navarro, and M. Viel, Journal of Cosmology and Astroparticle Physics 2014, 053 (2014), arXiv:1410.6813.
  * Padmanabhan and Refregier (2017) H. Padmanabhan and A. Refregier, Mon. Not. Roy. Astron. Soc. 464, 4008 (2017), arXiv:1607.01021 [astro-ph.CO] .
  * Hamaide  _et al._ (2022) L. Hamaide, H. Müller, and D. J. E. Marsh, Phys. Rev. D 106, 123509 (2022), arXiv:2210.03705 [astro-ph.CO] .
  * Li  _et al._ (2014) B. Li, T. Rindler-Daller, and P. R. Shapiro, Phys. Rev. D 89, 083536 (2014), arXiv:1310.6061 [astro-ph.CO] .
  * Gorghetto  _et al._ (2022) M. Gorghetto, E. Hardy, J. March-Russell, N. Song, and S. M. West, JCAP 08, 018 (2022), arXiv:2203.10100 [hep-ph] .
  * Marsh (2016a) D. J. E. Marsh, Physics Reports 643, 1 (2016a), arXiv:1510.07633.
  * Schwabe  _et al._ (2020) B. Schwabe, M. Gosenca, C. Behrens, J. C. Niemeyer, and R. Easther, Physical Review D 102, 083518 (2020), arXiv:2007.08256.
  * Lague  _et al._ (2022) A. Lague _et al._ , In preparation (2022).
  * Preskill  _et al._ (1983) J. Preskill, M. B. Wise, and F. Wilczek, Phys. Lett. B 120, 127 (1983).
  * Abbott and Sikivie (1983) L. Abbott and P. Sikivie, Phys. Lett. B 120, 133 (1983).
  * Dine and Fischler (1983) M. Dine and W. Fischler, Phys. Lett. B 120, 137 (1983).
  * Ma and Bertschinger (1995) C.-P. Ma and E. Bertschinger, Astrophys. J. 455, 7 (1995), arXiv:astro-ph/9506072 .
  * chan Hwang and Noh (2009) J. chan Hwang and H. Noh, Physics Letters B 680, 1 (2009).
  * Khlopov  _et al._ (1985) M. Khlopov, B. A. Malomed, and I. B. Zeldovich, Mon. Not. Roy. Astron. Soc. 215, 575 (1985).
  * Hu  _et al._ (2000) W. Hu, R. Barkana, and A. Gruzinov, Physical Review Letters 85, 1158 (2000), arXiv:astro-ph/0003365.
  * Amendola and Barbieri (2006) L. Amendola and R. Barbieri, Phys. Lett. B 642, 192 (2006), arXiv:hep-ph/0509257 .
  * Marsh and Ferreira (2010) D. J. E. Marsh and P. G. Ferreira, Phys. Rev. D 82, 103528 (2010), arXiv:1009.3501 [hep-ph] .
  * Hlozek  _et al._ (2015b) R. Hlozek, D. Grin, D. J. E. Marsh, and P. G. Ferreira, Physical Review D 91, 103512 (2015b).
  * Widrow and Kaiser (1993) L. M. Widrow and N. Kaiser, Astrophys. J. Lett. 416, L71 (1993).
  * Schive  _et al._ (2014a) H.-Y. Schive, T. Chiueh, and T. Broadhurst, Nature Physics 10, 496 (2014a), arXiv:1406.6586.
  * Mocz  _et al._ (2017a) P. Mocz, M. Vogelsberger, V. H. Robles, J. Zavala, M. Boylan-Kolchin, A. Fialkov, and L. Hernquist, Monthly Notices of the Royal Astronomical Society 471, 4559 (2017a).
  * Levkov  _et al._ (2018) D. Levkov, A. Panin, and I. Tkachev, Phys. Rev. Lett. 121, 151301 (2018), arXiv:1804.05857 [astro-ph.CO] .
  * Mocz  _et al._ (2017b) P. Mocz, M. Vogelsberger, V. H. Robles, J. Zavala, M. Boylan-Kolchin, A. Fialkov, and L. Hernquist, Monthly Notices of the Royal Astronomical Society 471, 4559 (2017b), https://academic.oup.com/mnras/article-pdf/471/4/4559/19609125/stx1887.pdf .
  * Dalal and Kravtsov (2022) N. Dalal and A. Kravtsov, Phys. Rev. D 106, 063517 (2022), arXiv:2203.05750 [astro-ph.CO] .
  * Mocz  _et al._ (2019) P. Mocz _et al._ , Phys. Rev. Lett. 123, 141301 (2019), arXiv:1910.01653 [astro-ph.GA] .
  * Gough and Uhlemann (2022) A. Gough and C. Uhlemann,  (2022), 10.21105/astro.2206.11918, arXiv:2206.11918 [astro-ph.CO] .
  * Nori and Baldi (2018) M. Nori and M. Baldi, Mon. Not. Roy. Astron. Soc. 478, 3935 (2018), arXiv:1801.08144 [astro-ph.CO] .
  * Veltmaat and Niemeyer (2016) J. Veltmaat and J. C. Niemeyer, Phys. Rev. D 94, 123523 (2016), arXiv:1608.00802 [astro-ph.CO] .
  * Chen  _et al._ (2021) J. Chen, X. Du, E. W. Lentz, D. J. E. Marsh, and J. C. Niemeyer, Phys. Rev. D 104, 083022 (2021), arXiv:2011.01333 [astro-ph.CO] .
  * Almgren  _et al._ (2013) A. Almgren, J. Bell, M. Lijewski, Z. Lukic, and E. Van Andel, Astrophys. J. 765, 39 (2013), arXiv:1301.4498 [astro-ph.IM] .
  * May and Springel (2021) S. May and V. Springel, Monthly Notices of the Royal Astronomical Society 506, 2603 (2021), arXiv:2101.01828 [astro-ph.CO] .
  * Mead  _et al._ (2020) A. J. Mead, T. Tröster, C. Heymans, L. Van Waerbeke, and I. G. McCarthy, Astronomy & Astrophysics 641, A130 (2020).
  * Navarro  _et al._ (1997) J. F. Navarro, C. S. Frenk, and S. D. M. White, The Astrophysical Journal 490, 493 (1997), arXiv: astro-ph/9611107.
  * Scoccimarro  _et al._ (2001) R. Scoccimarro, R. K. Sheth, L. Hui, and B. Jain, The Astrophysical Journal 546, 20 (2001), arXiv: astro-ph/0006319v2.
  * Press and Schechter (1974) W. H. Press and P. Schechter, Astrophys. J. 187, 425 (1974).
  * Sheth and Tormen (2002) R. K. Sheth and G. Tormen, Monthly Notices of the Royal Astronomical Society 329, 61 (2002), arXiv:astro-ph/0105113.
  * Bryan and Norman (1998) G. L. Bryan and M. L. Norman, The Astrophysical Journal 495, 80 (1998).
  * Bullock  _et al._ (2001) J. S. Bullock, T. S. Kolatt, Y. Sigad, R. S. Somerville, A. V. Kravtsov, A. A. Klypin, J. R. Primack, and A. Dekel, Monthly Notices of the Royal Astronomical Society 321, 559 (2001), arXiv:astro-ph/9908159.
  * Singh and Ma (2003) S. Singh and C.-P. Ma, Phys. Rev. D 67, 023506 (2003).
  * Ringwald and Wong (2004) A. Ringwald and Y. Y. Y. Wong, Journal of Cosmology and Astroparticle Physics 2004, 005 (2004).
  * Villaescusa-Navarro  _et al._ (2018) F. Villaescusa-Navarro, S. Genel, E. Castorina, A. Obuljen, D. N. Spergel, L. Hernquist, D. Nelson, I. P. Carucci, A. Pillepich, F. Marinacci, B. Diemer, M. Vogelsberger, R. Weinberger, and R. Pakmor, The Astrophysical Journal 866, 135 (2018), arXiv:1804.09180.
  * Villaescusa-Navarro  _et al._ (2014) F. Villaescusa-Navarro, F. Marulli, M. Viel, E. Branchini, E. Castorina, E. Sefusatti, and S. Saito, Journal of Cosmology and Astroparticle Physics 2014, 011 (2014), arXiv:1311.0866.
  * Castorina  _et al._ (2014) E. Castorina, E. Sefusatti, R. K. Sheth, F. Villaescusa-Navarro, and M. Viel, Journal of Cosmology and Astroparticle Physics 2014, 049 (2014), arXiv:1311.1212.
  * Costanzi  _et al._ (2013) M. Costanzi, F. Villaescusa-Navarro, M. Viel, J.-Q. Xia, S. Borgani, E. Castorina, and E. Sefusatti, Journal of Cosmology and Astroparticle Physics 2013, 012 (2013), arXiv:1311.1514.
  * Schive  _et al._ (2014b) H.-Y. Schive, M.-H. Liao, T.-P. Woo, S.-K. Wong, T. Chiueh, T. Broadhurst, and W.-Y. P. Hwang, Physical Review Letters 113, 261302 (2014b), arXiv:1407.7762.
  * Niemeyer (2020) J. C. Niemeyer, Progress in Particle and Nuclear Physics 113, 103787 (2020).
  * Eggemeier  _et al._ (2022) B. Eggemeier, B. Schwabe, J. C. Niemeyer, and R. Easther, Physical Review D 105, 023516 (2022).
  * Marsh (2016b) D. J. E. Marsh, (2016b), arXiv:1605.05973 [astro-ph.CO] .
  * Hannestad  _et al._ (2020) S. Hannestad, A. Upadhye, and Y. Y. Wong, Journal of Cosmology and Astroparticle Physics 2020, 062 (2020).
  * Brandbyge  _et al._ (2008) J. Brandbyge, S. Hannestad, T. Haugbølle, and B. Thomsen, Journal of Cosmology and Astroparticle Physics 2008, 020 (2008).
  * Bird  _et al._ (2012) S. Bird, M. Viel, and M. G. Haehnelt, Monthly Notices of the Royal Astronomical Society 420, 2551 (2012).
  * Nori  _et al._ (2019) M. Nori, R. Murgia, V. Iršič, M. Baldi, and M. Viel, Mon. Not. Roy. Astron. Soc. 482, 3227 (2019), arXiv:1809.09619 [astro-ph.CO] .
  * Freundlich  _et al._ (2016) J. Freundlich, A. El-Zant, and F. Combes, SF2A-2016: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics , 153 (2016).
  * Mead  _et al._ (2015b) A. J. Mead, J. A. Peacock, C. Heymans, S. Joudaki, and A. F. Heavens, Monthly Notices of the Royal Astronomical Society 454, 1958 (2015b).
  * Schive  _et al._ (2016) H.-Y. Schive, T. Chiueh, T. Broadhurst, and K.-W. Huang, Astrophys. J. 818, 89 (2016), arXiv:1508.04621 [astro-ph.GA] .
  * Corasaniti  _et al._ (2017) P. S. Corasaniti, S. Agarwal, D. J. E. Marsh, and S. Das, Phys. Rev. D 95, 083512 (2017), arXiv:1611.05892 [astro-ph.CO] .
  * Luu  _et al._ (2020) H. N. Luu, S. H. H. Tye, and T. Broadhurst, Phys. Dark Univ. 30, 100636 (2020), arXiv:1811.03771 [astro-ph.GA] .
  * Gosenca  _et al._ (2023) M. Gosenca, A. Eberhardt, Y. Wang, B. Eggemeier, E. Kendall, J. L. Zagorac, and R. Easther, (2023), arXiv:2301.07114 [astro-ph.CO] .
  * Dentler  _et al._ (2022) M. Dentler, D. J. E. Marsh, R. Hložek, A. Laguë, K. K. Rogers, and D. Grin, Mon. Not. Roy. Astron. Soc. 515, 5646 (2022), arXiv:2111.01199 [astro-ph.CO] .
  * Grin  _et al._ (2019) D. Grin, M. A. Amin, V. Gluscevic, R. Hlǒzek, D. J. E. Marsh, V. Poulin, C. Prescod-Weinstein, and T. L. Smith, (2019), arXiv:1904.09003 [astro-ph.CO] .
  * Ade  _et al._ (2019) P. Ade _et al._ (Simons Observatory), JCAP 02, 056 (2019), arXiv:1808.07445 [astro-ph.CO] .
  * Amendola  _et al._ (2018) L. Amendola _et al._ , Living Rev. Rel. 21, 2 (2018), arXiv:1606.00180 [astro-ph.CO] .
  * Hunter (2007) J. D. Hunter, Computing in Science Engineering 9, 90 (2007).
  * Harris  _et al._ (2020) C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del R’ıo, M. Wiebe, P. Peterson, P. G’erard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Nature 585, 357 (2020).
  * Virtanen  _et al._ (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. Jarrod Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. Carey, İ. Polat, Y. Feng, E. W. Moore, J. Vand erPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and S. . . Contributors, Nature Methods 17, 261 (2020).
  * Astropy Collaboration  _et al._ (2022) Astropy Collaboration, A. M. Price-Whelan, P. L. Lim, N. Earl, N. Starkman, L. Bradley, D. L. Shupe, A. A. Patil, L. Corrales, C. E. Brasseur, M. N”othe, A. Donath, E. Tollerud, B. M. Morris, A. Ginsburg, E. Vaher, B. A. Weaver, J. Tocknell, W. Jamieson, M. H. van Kerkwijk, T. P. Robitaille, B. Merry, M. Bachetti, H. M. G”unther, T. L. Aldcroft, J. A. Alvarado-Montes, A. M. Archibald, A. B’odi, S. Bapat, G. Barentsen, J. Baz’an, M. Biswas, M. Boquien, D. J. Burke, D. Cara, M. Cara, K. E. Conroy, S. Conseil, M. W. Craig, R. M. Cross, K. L. Cruz, F. D’Eugenio, N. Dencheva, H. A. R. Devillepoix, J. P. Dietrich, A. D. Eigenbrot, T. Erben, L. Ferreira, D. Foreman-Mackey, R. Fox, N. Freij, S. Garg, R. Geda, L. Glattly, Y. Gondhalekar, K. D. Gordon, D. Grant, P. Greenfield, A. M. Groener, S. Guest, S. Gurovich, R. Handberg, A. Hart, Z. Hatfield-Dodds, D. Homeier, G. Hosseinzadeh, T. Jenness, C. K. Jones, P. Joseph, J. B. Kalmbach, E. Karamehmetoglu, M. Kaluszy’nski, M. S. P. Kelley, N. Kern, W. E. Kerzendorf, E. W. Koch, S. Kulumani, A. Lee, C. Ly, Z. Ma, C. MacBride, J. M. Maljaars, D. Muna, N. A. Murphy, H. Norman, R. O’Steen, K. A. Oman, C. Pacifici, S. Pascual, J. Pascual-Granado, R. R. Patil, G. I. Perren, T. E. Pickering, T. Rastogi, B. R. Roulston, D. F. Ryan, E. S. Rykoff, J. Sabater, P. Sakurikar, J. Salgado, A. Sanghi, N. Saunders, V. Savchenko, L. Schwardt, M. Seifert-Eckert, A. Y. Shih, A. S. Jain, G. Shukla, J. Sick, C. Simpson, S. Singanamalla, L. P. Singer, J. Singhal, M. Sinha, B. M. SipHocz, L. R. Spitler, D. Stansby, O. Streicher, J. Sumak, J. D. Swinbank, D. S. Taranu, N. Tewary, G. R. Tremblay, M. d. Val-Borro, S. J. Van Kooten, Z. Vasovi’c, S. Verma, J. V. de Miranda Cardoso, P. K. G. Williams, T. J. Wilson, B. Winkel, W. M. Wood-Vasey, R. Xue, P. Yoachim, C. Zhang, A. Zonca, and Astropy Project Contributors, apj 935, 167 (2022), arXiv:2206.14220 [astro-ph.IM] .
  * Mead  _et al._ (2021b) A. Mead, S. Brieden, T. Tröster, and C. Heymans, Monthly Notices of the Royal Astronomical Society 502, 1401 (2021b), arXiv: 2009.01858.
  * Mead  _et al._ (2016) A. J. Mead, C. Heymans, L. Lombriser, J. A. Peacock, O. I. Steele, and H. A. Winther, Monthly Notices of the Royal Astronomical Society 459, 1468 (2016).
  * Copeland  _et al._ (2020) D. Copeland, A. Taylor, and A. Hall, Monthly Notices of the Royal Astronomical Society 493, 1640 (2020).
  * Smith  _et al._ (2003) R. E. Smith, J. A. Peacock, A. Jenkins, S. D. M. White, C. S. Frenk, F. R. Pearce, P. A. Thomas, G. Efstathiou, and H. M. P. Couchman, Monthly Notices of the Royal Astronomical Society 341, 1311 (2003).
  * Emberson  _et al._ (2017) J. D. Emberson _et al._ , Res. Astron. Astrophys. 17, 085 (2017), arXiv:1611.01545 [astro-ph.CO] .
  * Villaescusa-Navarro  _et al._ (2013) F. Villaescusa-Navarro, S. Bird, C. Peña-Garay, and M. Viel, Journal of Cosmology and Astroparticle Physics 2013, 019 (2013).
  * Schmidt (2016) F. Schmidt, Physical Review D 93, 063512 (2016).