# Constraining Ultralight Axions with Galaxy Surveys

A. Laguë,11footnotetext: Corresponding author. J. R. Bond  R. Hložek  K. K.
Rogers  D. J. E. Marsh  and D. Grin

###### Abstract

Ultralight axions and other bosons are dark matter candidates present in many
high energy physics theories beyond the Standard Model. In particular, the
string axiverse postulates the existence of up to $\mathcal{O}(100)$ light
scalar bosons constituting the dark sector. Considering a mixture of axions
and cold dark matter, we obtain upper bounds for the axion relic density
$\Omega_{a}h^{2}<0.004$ for axions of mass $10^{-31}\;\mathrm{eV}\leq
m_{a}\leq 10^{-26}\;\mathrm{eV}$ at 95% confidence. We also improve existing
constraints by a factor of over 4.5 and 2.1 for axion masses of $10^{-25}$ eV
and $10^{-32}$ eV, respectively. We use the Fourier-space galaxy clustering
statistics from the Baryon Oscillation Spectroscopic Survey (BOSS) and
demonstrate how galaxy surveys break important degeneracies in the axion
parameter space compared to the cosmic microwave background (CMB). We test the
validity of the effective field theory of large-scale structure approach to
mixed ultralight axion dark matter by making our own mock galaxy catalogs and
find an anisotropic ultralight axion signature in the galaxy quadrupole. We
also observe an enhancement of the linear galaxy bias from 1.8 to 2.4 when
allowing for 5% of the dark matter to be composed of a $10^{-28}$ eV axion in
our simulations. Finally, we develop an augmented interpolation scheme
allowing a fast computation of the axion contribution to the linear matter
power spectrum leading to a 70% reduction of the computational cost for the
full Monte Carlo Markov chains analysis.

##  1 Introduction

Multiple hypotheses have been put forward to explain the nature of dark
matter. Among the most promising candidates from particle physics are axions
which were initially theorized as a solution to the charge parity problem in
quantum chromodynamics (QCD). A cosmological population of such axions
produced by vacuum realignment [1, 2, 3] would have a relic density given by
[4]

$\displaystyle\Omega_{a}^{\mathrm{QCD}}h^{2}\approx
0.12\left(\frac{6\;\mu\mathrm{eV}}{m_{a}}\right)^{1.165}\theta_{I}^{2},$ (1.1)

where $\theta_{I}$ is the initial misalignment angle and $h$ is the Hubble
constant today in units of 100 km/s/Mpc, and $m_{a}$ is the axion mass. We
note that the QCD axion has a lower mass bound of around $m_{a}\sim 6\;\mu$eV
for $\theta_{I}=\mathcal{O}(1)$ to avoid the relic density from exceeding the
cosmological dark matter density $\Omega_{d}h^{2}\approx 0.12$. Escapes from
this can be found in the so-called “anthropic” axion window by fine-tuning the
initial axion field. Indeed, assuming the PQ symmetry is broken after
inflation, it is possible to prevent the lighter QCD axion from exceeding the
total dark matter relic density by appropriately reducing the initial
misalignment angle [5].

On the other hand, axions also appear naturally in string theory where they
can be many orders of magnitude lighter than their QCD counterparts leading to
the concept of an axiverse [6]. The relic density in this regime follows
(assuming again a vacuum realignment formation mechanism) [7]

$\displaystyle\Omega_{a}h^{2}=0.12\left(\frac{m_{a}}{4.7\times
10^{-19}\;\mathrm{eV}}\right)^{1/2}\left(\frac{f_{a}}{10^{16}\;\mathrm{GeV}}\right)^{2}\left(\frac{\Omega_{m}h^{2}}{0.15}\right)^{3/4}\left(\frac{3.4\times
10^{3}}{1+z_{\mathrm{eq}}}\right)^{3/4}\theta_{I}^{2},$ (1.2)

where $f_{a}$ is the symmetry breaking scale of the axion, $\Omega_{m}$ is the
total matter density, and $z_{\mathrm{eq}}$ is the redshift of matter-
radiation equality. In the latter case, they are theorized to form in a
plenitude of $\mathcal{O}(100)$ logarithmically distributed masses [6, 8].
Ref. [9] studied the axion mass and decay constant distributions in Type IIB
compactifications on Calabi-Yau manifolds and found and almost log-flat
distributions for $m_{a}$, and log-normal distributions for $f_{a}$. The
implications of a single ultralight scalar boson of mass $m_{a}\sim 10^{-22}$
eV have been explored in the context of fuzzy dark matter (FDM) (see e.g. Ref.
[10] for a review). This scenario assumes however that a single light boson
composes exactly all of the dark matter. Given the form of Eq. (1.2), we can
ask what happens if we relax this assumption and allow multiple axions each
composing a fraction of the dark matter. Many studies have been devoted to
probing the existence of such particles as a subdominant component [11, 12]
and any evidence for such a multi-component dark sector would be considered a
smoking gun for the string landscape [6].

There is also motivation for the existence of light bosonic dark matter from a
cosmological perspective as tensions arise between the predictions of cold
dark matter (CDM)-only simulations and observations. This is known as the
small-scale crisis of CDM. Notably, there is the missing satellite problem
where the number of satellite galaxies predicted by CDM is much larger than
what has been historically observed [13, 14]. There is also the cusp-core
problem where halos in CDM simulations have a diverging density profile at
their center when observations seem to suggest the existence of central cores
with an approximately constant density [15]. Finally, there is the too-big-to-
fail problem which describes the fact that some CDM subhalos are theoretically
too massive not to have formed stars and should be observationally detectable
[16]. Many modifications to CDM have been proposed to alleviate these tensions
such as self-interacting dark matter [17] and warm dark matter [18]. Also,
detailed hydrodynamical simulations have shown that baryonic effects such as
feedback and active galactic nuclei solve many of these issues [19, 20, 21,
22]. For ultralight axions, it has been shown that their lack of small-scale
clustering and formation of solitonic cores could offer a solution to all of
the above tensions [23, 24]. It has been argued in previous work that the
existence of ultralight bosons of mass $m_{a}\sim 10^{-22}$ eV could
potentially exacerbate the cusp-core problem by creating overly dense cores
[25] and predict dwarf galaxies that are too massive [26]. However, these
conclusions are based on the assumption that the dark matter is composed of a
single axion species whereas high energy physics naturally allows a number of
axion species $N_{\mathrm{ax}}\gg 1$.

Existing bounds on ultralight axion dark matter ($m_{a}\lesssim
10^{-20}\,\mathrm{eV}$) exploit their distinctive effect on the cosmic large-
scale structure. This is most prominently a suppression in the growth of
structure, relative to the CDM limit, below a characteristic scale owing to
so-called “quantum pressure.” This manifests as a cut-off in the linear matter
power spectrum, where the cut-off wavenumber is a monotonically increasing
function of the axion mass [27]. On the largest observable scales, the cosmic
microwave background (CMB) rules out axions being all the dark matter for
$10^{-33}\,\mathrm{eV}\leq m_{a}\leq 10^{-24}\,\mathrm{eV}$ and places
percent-level bounds on their fractional contribution to the dark matter
[28].222The CMB is also sensitive to axion-induced isocurvature modes and the
mass regime where ultralight axions behave like dark energy. On the smallest
scales currently accessible in the quasi-linear matter power spectrum [$k\sim
20\,h\,\mathrm{Mpc}^{-1}$; 29, 30], the Lyman-alpha forest [31, 32, 33] sets
the strongest bound on the mass scale ($m_{a}\sim 10^{-22}\,\mathrm{eV}$)
motivated above, allowing axions to be all the dark matter only if
$m_{a}>2\times 10^{-20}\,\mathrm{eV}$ [34, 35]. In setting this bound, it is
important to marginalise carefully the partly degenerate suppression scales
arising from the temperature and pressure of the intergalactic medium [36, 37,
35]. Competitive bounds also arise from not observing an axion-induced
suppression in the Milky Way sub-halo mass function ($m_{a}>2.9\times
10^{-21}\,\mathrm{eV}$), with a different set of nuisance parameters and
systematic effects [38]. Also note that the Lyman-$\alpha$ forest bounds the
axion fraction in the mixed dark matter scenarios, but only at the $\approx
30$% level [32]. X-ray observations have also been used to make constraints on
ultralight axions as dark matter reaching a constraint of $m_{a}>7\times
10^{-23}$ eV [39]. In this work, we exploit a different regime in the matter
power spectrum, as well as, for the first time, the anisotropic effect on
matter clustering.

Figure 1: Constrained region of the $m_{a}-\Omega_{a}/\Omega_{d}$ parameter
space when combining galaxy clustering and CMB data. The light and dark red
regions represent the $2\sigma$ and $1\sigma$ allowed regions respectively.
The grey shaded area denotes the boundary of the $2\sigma$ allowed region
using the Planck 2015 data release from the analysis of Ref. [40] and the grey
error bars denote the 2$\sigma$ bound found using only galaxy surveys as part
of the present work. Previous work combining the CMB and galaxy clustering
data found $\Omega_{a}/\Omega_{d}\leq 0.048$ for a similar axion mass range
[40] while we find that a $10^{-28}$ eV axion has $\Omega_{a}/\Omega_{d}\leq
0.014$ at 95% confidence. The wider range of allowed $\Omega_{a}/\Omega_{d}$
when adding CMB prior at $10^{-25}$ eV stems from a degeneracy with the power
spectrum amplitude (see Fig. 16).

Galaxies are tracers of the matter power spectrum and deviations from the CDM
power spectrum due to the presence of ultralight axions has a characteristic
signature which we can probe using galaxy clustering data. In order to
constrain small concentrations of ultralight axions, we measure the maximal
amount of suppression of the galaxy power spectrum from the BOSS DR12 survey
for various axion scenarios. Without loss of generality, we consider only one
ultralight axion species while the remainder of the dark matter is treated as
pure CDM. This model is equivalent to a multi-axion scenario where only one of
the axion masses falls in the ultralight regime. Deviations from CDM would be
even larger in the presence of multiple ultralight fields (assuming the same
relic density) which implies that the single axion + CDM model in this study
will lead to conservative constraints. In Section 2, we describe the physics
of ultralight axions and their behaviour in more detail. Next, in Section 3,
we present the galaxy clustering model used in this work while the full set of
parameters and priors are described in Section 4. In Section 5, we validate
our model with large-scale structure simulations both with and without an
axionic component. Finally, we consider the results of the full likelihood
analysis in Section 6, and discuss the implications of our findings in Section
7. In Appendix A we develop our model in greater detail and examine the impact
of the wave effects on higher order corrections to the power spectrum, in
Appendix B we study the impact of the power spectrum suppression on the galaxy
quadrupole, and Appendix C we describe our optimized axion transfer function
interpolation scheme.

##  2 Ultralight Axion Physics

Ultralight axions are scalar bosons described by a non-relativistic field
$\phi$ with a cosine potential which we approximate as a harmonic quadratic
potential of the form

$\displaystyle V(\phi)$ $\displaystyle\propto 1-\cos(\phi/f_{a})$ (2.1)
$\displaystyle\approx\frac{1}{2}m_{a}^{2}\phi^{2}.$ (2.2)

This assumes $\phi\ll f_{a}$ which is consistent with the treatment of Ref.
[40, 28]333See Ref. [41] for certain cases where $\phi~{}\sim f_{a}$ could be
advantageous.. Across the mass range of interest, this approximation is valid
for $f_{a}$ slightly above the Grand Unified Theory (GUT) scale (see e.g. Ref.
[42]). The corresponding density and pressure of the background axion field
are then given respectively by

$\displaystyle\rho_{a}$
$\displaystyle=\frac{1}{2a^{2}}\dot{\phi_{0}}^{2}+\frac{m_{a}^{2}}{2}\phi_{0}^{2},$
(2.3) $\displaystyle P_{a}$
$\displaystyle=\frac{1}{2a^{2}}\dot{\phi_{0}}^{2}-\frac{m_{a}^{2}}{2}\phi_{0}^{2},$
(2.4)

where the dot denotes the derivative with respect to conformal time and
$\phi_{0}$ is the homogeneous background field. The equation of motion of the
field for the field reads

$\displaystyle\ddot{\phi_{0}}+2\mathcal{H}\dot{\phi_{0}}+\frac{1}{2}a^{2}m_{a}^{2}\phi_{0}=0,$
(2.5)

where $\mathcal{H}=aH$ is the conformal Hubble factor. At early times ($H\gg
m_{a}$), the field is slowly rolling with $\dot{\phi}\approx 0$, the axion
equation of state is $w_{a}\equiv P_{a}/\rho_{a}\approx-1$, and the field
behaves as a dark energy component. At late times ($H\ll m_{a}$), the field’s
equation of state oscillates around zero and the energy density of the field
obeys the scaling $\rho_{a}\propto a^{-3}$. It is therefore a dark matter
component after the field started oscillating around the potential minimum.
The value of the scale factor when this transition happens is denoted
$a_{\mathrm{osc}}$ and defined such that

$\displaystyle 3H(a_{\mathrm{osc}})\approx m_{a}.$ (2.6)

This implies that the time of transition of the field is given as a function
of its mass. For $m_{a}\lesssim 10^{-28}$ eV, and
$a_{\mathrm{osc}}>a_{\mathrm{eq}}$ the field does not behave as dark matter at
matter-radiation equality and cannot constitute the entirety of the dark
matter. We refer to those axions as dark-energy-like axions. For any mass, the
mean axion density depends on $a_{\mathrm{osc}}$ and is given by

$\displaystyle\rho_{a}=a_{\mathrm{osc}}^{3}\bigg{[}\frac{1}{2a^{2}}\dot{\phi_{0}}^{2}+\frac{m_{a}^{2}}{2}\phi_{0}^{2}\bigg{]}_{m_{a}=3H}.$
(2.7)

One of the main signatures of axions is their lack of clustering on small
scales. Axions have a characteristic Jeans scale $k_{J}$ below which the
growth of structure is suppressed444See Ref. [43] for a study of the Jeans
instability in the presence of self-interactions.. For the linear matter power
spectrum with dark-matter-like axions, the suppression is frozen in at matter-
radiation equality and we have [27]

$\displaystyle k_{J,\mathrm{eq}}\approx
9\bigg{(}\frac{m_{a}}{10^{-22}\;\mathrm{eV}}\bigg{)}^{1/2}\mathrm{Mpc}^{-1}.$
(2.8)

Figure 2: Linear power spectrum for different axion masses and densities. We
note that the level of suppression increases with the axion density while the
scale at which the axion power spectrum deviates from the CDM ($\Omega_{a}=0$)
case is dependent on the axion mass. The grey shaded area highlights the
Fourier modes probed by our analysis between $k_{\mathrm{min}}=0.025\;h/$Mpc
and $k_{\mathrm{max}}=0.23\;h/$Mpc. The black arrow indicates the Jeans scale
for the only axion mass for which $a_{\mathrm{osc}}<a_{\mathrm{eq}}$ among the
three showed. For reference, the massive neutrino characteristic scale of
$k_{\mathrm{nr}}$ for a neutrino with a mass of $150$ eV would correspond to
the Jeans scale of axions of mass $10^{-26}$ eV.

This holds when assuming that the axions behave as dark matter at matter-
radiation equality or that $m_{a}\gtrsim 3H(a_{\mathrm{osc}})$. In the
alternate case, the field behaves as a dark energy component until it starts
oscillating. The Jeans scale at matter-radiation equality is not the optimal
reference scale for these axions. We need instead to use the scale at horizon
crossing when the field started oscillating which we can calculate as [44, 6,
42]

$\displaystyle k_{m}=a_{\mathrm{osc}}H\left(a_{\mathrm{osc}}\right).$ (2.9)

We note that in both cases, the loss in power on small scales is dependent on
the axion fraction $\Omega_{a}/\Omega_{d}$ where $\Omega_{d}$ is the total
dark matter density (combining CDM and axions). The level of suppression may
then only be partial when fixing the axion mass and varying the axion
concentration with the $\Omega_{a}\to 0$ limit recovering the pure CDM result.
This deviation from CDM in the linear power spectrum illustrated in Fig. 2 is
the main signature we use in this study to constrain axions. The matter power
spectra including axionic effects are computed using the adapted Boltzmann
code axionCAMB [40]. It is worth noting that the code solves for the evolution
of the axion perturbations using an effective fluid approach. It has been
shown to match closely the exact solution, although biases could occur if the
sensitivity of the experiments used approaches the cosmic variance limit [45].

##  3 Model

Our datasets consist of the multipole moments of the CMASS
($z_{\mathrm{eff}}=0.57$) and LOWZ ($z_{\mathrm{eff}}=0.32$) catalogs from the
BOSS survey [46]. We use the North and South Galactic Caps denoted NGC and SGC
for our analysis. Given galaxy positions, the local galaxy overdensity is
defined as

$\displaystyle\delta_{g}(\mathbf{x})\equiv\frac{n_{g}(\mathbf{x})}{\bar{n}_{g}}-1,$
(3.1)

where $n_{g}(\mathbf{x})$ is the number density of galaxies at the comoving
position $\mathbf{x}$ and $\bar{n}_{g}$ is the mean number density for the
survey volume $V$. After taking the Fourier transform $\delta_{g}(\mathbf{k})$
of the real-space overdensity, the multipole moments $P_{\ell}$ of the power
spectrum are calculated by integrating over the angle between the line-of-
sight and the wavevector $\mathbf{k}$. It is customary to decompose the latter
into $k$ which is the scalar amplitude of the vector and $\mu$ which is the
cosine of angle between the wavevector and the line-of-sight. The integral
expression for the galaxy multipoles can then be written as (using the
Yamamoto estimator [47])

$\displaystyle
P_{\ell}(k)=\left\langle\frac{2\ell+1}{V}\int\frac{d\Omega_{k}}{4\pi}\delta_{g}(\mathbf{k})\delta_{g}(-\mathbf{k})\mathcal{P}_{\ell}\left(\hat{\mathbf{k}}\cdot\hat{\mathbf{z}}\right)\right\rangle,$
(3.2)

where the hat denotes unit vectors and $\mathcal{P}_{\ell}$ are the Legendre
polynomials of degree $\ell$. In the data, the shot noise
$P_{\mathrm{shot}}=1/\bar{n}_{g}$ is subtracted from the monopole. For the
purpose of this study, we will focus on the monopole ($\ell=0$) and quadrupole
($\ell=2$) moments of the power spectrum.

Figure 3: Monopoles and quadrupoles of the galaxy clustering obtained from the
linear power spectra of Fig. 2. The model is evaluated with fixed cosmological
and bias parameters at redshift $z=0.57$. We highlight a slight increase in
the quadrupole on small scales indicating that the structure suppression from
axions also has anisotropic effects.

Our model for the multipoles begins with the linear matter power spectrum
generated from a set of standard flat $\Lambda$CDM cosmological parameters and
a set of two axion parameters. We first obtain the power spectrum without the
axion effects using the Boltzmann code Class [48] and the axion transfer
function with the adapted code axionCAMB. We refer to this power spectrum as
the CDM power spectrum. The reason for this use of code combination is to
separate the calculations for the axion effects from the rest of the
cosmological calculations since they take significantly longer to complete. We
opt to interpolate over the axion transfer function rather than the full
matter power spectrum since the axion transfer function is independent of
$A_{s}$. This implies that we can reduce the dimensionality of our
interpolation tables for the axion transfer function by interpolating only
over the axion fraction, the Hubble constant, the baryon density, and the dark
matter density. This saves considerable computational resources. We obtain
percent-level errors compared to the use of the full axionCAMB when the axion
fraction is less than 10%. As for the choice of Boltzmann code for the
$\Lambda$CDM part of the calculation, we can use any code of preference since
the CDM power spectrum has been factored out of the axion transfer function
and since the latter is only weakly dependent on the choice of cosmology at
low axion fractions (we explore this more rigorously in Appendix C).
Ultimately, the choice of the Class code was made due to the fact that it is
the one implemented with the galaxy multipole calculator which we will use for
our analysis.

We then load the interpolation tables and use them to compute the axion
transfer function so that the linear power spectrum with axions is given by

$\displaystyle
P_{\mathrm{lin}}(k,z)=T_{\mathrm{ax}}^{2}\left(k,z;m_{a},\Omega_{a}/\Omega_{d}\right)P_{\mathrm{lin,\;CDM}}(k,z).$
(3.3)

Examples of $P_{\mathrm{lin}}$ and $T_{\mathrm{ax}}^{2}$ are given in the top
and bottom rows of Fig. 2 respectively. The computation of the axion transfer
function interpolation tables is given in more detail in Appendix C.

The redshift-space galaxy power spectrum can be summarized as a sum of
corrections to the (biased) linear matter power spectrum (including axion
effects)

$\displaystyle
P_{g}(\mathbf{k})=b_{1}^{2}P_{\mathrm{lin}}(k)+P_{\mathrm{1-loop}}(\mathbf{k})+P_{\mathrm{counter}}(\mathbf{k})+P_{\mathrm{stoch}}(k),$
(3.4)

where the contributions are respectively from the 1-loop corrections, the
counterterms, and the stochastic terms. This model is based on the Effective
Field Theory of Large Scale Structure (EFTofLSS) [49] which allows high
accuracy prediction of the power spectrum up to semi-linear scales. This
approach consists in perturbatively expanding the equations for the power
spectrum corrections in powers of $k/k_{\mathrm{NL}}$ where $k_{\mathrm{NL}}$
is the scale at which non-linear effects begin to dominate. It accounts for
physics on smaller scales through the counterterms while stochastic terms
model the variations between the expectation values of the stress tensor and
its value in a given realization. For more details on the EFTofLSS, we direct
the readers to Ref. [49] and to subsequent publications.

The EFTofLSS was designed with the assumption of collisionless dark matter. In
the case of ultralight axions of mass $m_{a}\leq 10^{-25}$ eV, the dark matter
particles do not cluster on small scales due to their scale-dependent sound
speed [50, 51, 52]. For this reason, we may treat them like fast moving
neutrinos even if the axions are strictly speaking non-relativistic. The
inclusion of fast moving neutrinos in the EFTofLSS was performed in Ref. [53].
It was found that at first order, the supplementary counterterms that need to
be added to account for neutrinos have the same functional form as the
counterterms already present. Given that we vary the prefactors of these
counterterms with a loose prior in our analysis, we assume that the functional
form of our counterterms accounts for the non-clustering axions. Therefore the
EFTofLSS can be extended to ultralight axions when marginalizing over the
amplitude of the counterterms (which are defined in more detail in Appendix
A). Furthermore, we explore the impact of the wave effects on the clustering
dynamics on higher order corrections in Appendix A while we leave the
computation of the counterterms including axion wave effects for future work.
We conclude that the deviations from pure CDM dynamics are a small
contribution to already subdominant higher order corrections, especially at
low axion fractions.

We then use the public PyBird code to apply the following transformations to
the linear matter power spectrum (with the corresponding reference in the
literature):

  1. 1.

compute and add the non-linear one-loop correction to the power spectrum [54],

  2. 2.

perform the IR-resummation of the power spectrum [55],

  3. 3.

apply window function effects [56],

  4. 4.

apply the Alcock-Paczynski effect [57] ,

  5. 5.

account for power spectrum binning [57],

  6. 6.

correct for fiber collisions [57].

In the following section, we will briefly discuss each of the steps performed
on the linear matter power spectrum to obtain a realistic model for the galaxy
multipole data. The calculations for each of these terms from the parameters
of Table 1 are explained in greater detail in Appendix A. Examples of galaxy
multipoles obtained with our model including axions with fixed nuisance
parameters are shown in Fig. 3.

The PyBird code has the convenient feature of computing the effect of window
functions from surveys. We use this when running the tests on the simulated
data. However, the dataset we use from BOSS combines the North and South
Galactic Caps for each redshift slice. We therefore use a custom window
function computation. We obtain the monopole and quadrupole with the publicly
available window matrices $W_{ij}$ using Ref. [58] Eq. (15)

$\displaystyle P_{\ell}(k_{i})=\sum_{j}W_{ij}^{\ell 0}P_{0}^{\mathrm{no-
window}}(k_{j})+\sum_{j}W_{ij}^{\ell 2}P_{2}^{\mathrm{no-window}}(k_{j}),$
(3.5)

where $k_{i}$ are the $k$-bins of the BOSS dataset and $k_{j}$ are 1000
selected $k$-bins up to 0.5 $h$/Mpc.

##  4 Methodology

Our model based on the EFTofLSS was described in Section 3 and has 10 free
parameters. The monopole and quadrupole are both calculated at redshifts
$z=0.32$ and $z=0.57$. For each likelihood calculation, the cosmological
parameters are assumed to be the same for the two redshift slices, but the
bias, counterterms, and stochastic parameters are allowed to vary. The axion
fraction is varied like the cosmological parameters while the axion mass is
fixed for each individual calculation. We thus make sixteen independent
likelihood calculations: for each of the eight mass bins between $10^{-32}$ eV
and $10^{-25}$ eV, we compute the posterior distribution for two choices of
priors. For the first half of the runs, we include only information about
galaxy clustering and use a Big Bang Nucleosynthesis (BBN) prior on the baryon
density $\omega_{b}$. For the second set of runs which combine our setup with
existing constraints, we instead impose a CMB prior on all cosmological
parameters along with a prior on the axion fraction also from the same CMB
analysis. We thus use the pre-computed chains from the Planck likelihood
including an axion component directly. This method has been used as an
alternative to joint likelihood analysis when combining large-scale structure
and CMB data with the EFTofLSS (see Ref. [59, 57] where a CMB prior is put on
the baryon density and sound horizon).

Parameter | Type | Prior
---|---|---
$100\omega_{b}$ | Cosmology | $\mathcal{N}(2.214,\;0.038)$
$\omega_{d}$ | Cosmology | $\mathcal{U}(0.05,\;0.2)$
$\ln(10^{10}A_{s})$ | Cosmology | $\mathcal{U}(\log(5),\;\log(50))$
$h$ | Cosmology | $\mathcal{U}(0.55,\;.91)$
$m_{a}$ | Axion | Discrete Bins
$\Omega_{a}/\Omega_{d}$ | Axion | $\mathcal{U}(0,\;1)$
$b_{1}$ | Bias | $\mathcal{U}(0,\;4)$
$b_{2}$ | Bias | $\mathcal{U}(-2\sqrt{2},\;2\sqrt{2})$
$b_{3}$ | Bias | $\mathcal{N}(0,\;2)$
$b_{4}$ | Bias | Fixed (equal to $b_{2}$)
$c_{\mathrm{ct}}$ | Counterterm | $\mathcal{N}(0,\;2)$
$\tilde{c}_{r,1}$ | Counterterm | $\mathcal{N}(0,\;8)$
$\tilde{c}_{r,2}$ | Counterterm | Fixed to 0
$\tilde{c}_{\epsilon,1}/1000$ | Stochastic | $\mathcal{U}(-10,\;10)$
$\tilde{c}_{\epsilon,2}/1000$ | Stochastic | $\mathcal{N}(0,\;3)$
Table 1:  Set of cosmological, axion, bias, and counterterm parameters with
their corresponding priors. $\mathcal{U}(a,b)$ indicates that the prior is a
uniform distribution with lower and upper bounds $a$ and $b$ and
$\mathcal{N}(\mu,\sigma)$ indicates that the prior is a Gaussian distribution
with mean $\mu$ and variance $\sigma^{2}$. The parameter $h$ is the Hubble
constant today in units of $100$ km/s/Mpc and the stochastic terms are quoted
in units of $[\mathrm{Mpc}/h]^{3}$.

As our dataset, we use the CMASS combined NGC and SGC as well as the LOWZ
combined NGC and SGC (made public by the BOSS collaboration [60]). The
covariance matrices are computed from the MultiDark-Patchy mock catalogs [61,
62]. These simulations were run using an augmented Lagrangian perturbation
theory algorithm combined with a bias calculation scheme to populate the
simulated halos with galaxies [63]. We fit the monopole ($\ell=0$) and
quadrupole ($\ell=2$) moments of the anisotropic galaxy power spectrum as
defined in Eq. (3.2). To calculate the posterior distribution for the
parameters of Table 1, we use the MCMC sampler emcee [64] with which we create
our own likelihood code adapted for axion parameters. The priors on the
cosmological parameters are chosen to match closely those of similar analyses
and are generic for $\Lambda$CDM parameters estimation using BOSS data [59].
The priors on the nuisance parameters are theoretical priors motivated by the
fact that EFT parameters are expected to be $\mathcal{O}(1)$ [65, 66].

We employ the same fitting method as Ref. [40] to obtain valid $\Lambda$CDM +
ultralight axions priors instead of the regular $\Lambda$CDM priors from the
Planck analysis. We check fo MCMC convergence using the spectral method of Ref
[67]. Finally, we compile all of the posterior distributions for the axion
fraction as a function of the axion mass for the galaxies-only and
galaxies+CMB scenarios and obtain a constraint plot shown in Fig. 1 for the
$m_{a}-\Omega_{a}/\Omega_{d}$.

##  5 Tests on Simulations

We test our model on two sets of simulated data. In the first, we use the
MultiDark-Patchy Mocks [61, 62] set of simulations provided with the BOSS
datasets. These are the standard simulations used to validate galaxy
clustering models in redshift space. We run our model on the monopoles and
quadrupoles extracted from the simulated CMASS NGC data in order to compare
the posterior distribution of our cosmological parameters with the true values
used in the simulations. For the second set, we use the large scale Lagrangian
based simulation code Peak-Patch [68] in order to generate our own set of
simulations. We introduce axions in one set of simulations and assume
$\Lambda$CDM conditions for the other. We extract the multipoles from the
simulations before running our model on this set of simulations as well. This
is to determine if axionic components in the dark matter can be revealed by
our model if they exist.

###  5.1 MultiDark-Patchy Mocks

In order to validate the model presented in Sec 3, we make use of the
MultiDark-Patchy mock simulations. These simulations are based on a
$\Lambda$CDM cosmology and do not include ultralight axions. We use the priors
in Table 1 and we test using two priors on $\omega_{b}$, with one from the CMB
and one from BBN. Using these two choices of priors, we make two analyses with
the axion fraction fixed to zero (assuming pure CDM) and two analyses allowing
the axion fraction to vary (assuming an axion mass of $m_{a}=10^{-27}$ eV),
for a total of four runs. We also set the sum of the neutrino masses to zero
since the simulations were run without massive neutrinos.

(a)

(b)

Figure 4: Model validation on high-$z$ NGC mock catalogs with a Planck (a) or BBN (b) prior on $\omega_{b}$. The dashed lines denote the input parameters of the simulations. CDM denotes the use of the $\Lambda$CDM model while the CDM + Axions denotes the extended CDM model where we also vary the axion fraction. The NGC mocks were run with a CDM cosmology. $\omega_{b}$ Prior | Planck Prior | BBN Prior
---|---|---
Parameter | Best-Fit | Mean & 68% C.L. | Best-Fit | Mean & 68% C.L.
$100\omega_{b}$ | $2.218$ | $2.215\pm 0.015$ | $2.220$ | $2.215\pm 0.038$
$\omega_{d}$ | $0.110$ | $0.1171^{+0.0082}_{-0.0094}$ | $0.118$ | $0.1163^{+0.0089}_{-0.010}$
$\ln(10^{10}A_{s})$ | $3.171$ | $3.03\pm 0.15$ | $3.055$ | $3.04^{+0.18}_{-0.16}$
$h$ | $0.665$ | $0.677^{+0.017}_{-0.019}$ | $0.676$ | $0.674\pm 0.019$
$b_{1}$ | $1.828$ | $1.97^{+0.15}_{-0.18}$ | $1.928$ | $1.96^{+0.16}_{-0.21}$
$b_{2}$ | $0.345$ | $0.59^{+0.23}_{-0.43}$ | $0.442$ | $0.53^{+0.27}_{-0.44}$
$b_{3}$ | $-1.667$ | $-0.7\pm 1.8$ | $-0.659$ | $-0.8^{+1.8}_{-1.6}$
$c_{\mathrm{ct}}$ | $-0.725$ | $-0.4^{+1.9}_{-2.1}$ | $0.615$ | $-0.4\pm 2.0$
$\tilde{c}_{r,1}$ | $-13.445$ | $-16.3^{+6.9}_{-4.4}$ | $-14.903$ | $-15.4^{+6.1}_{-5.2}$
$\tilde{c}_{\epsilon,1}/1000$ | $3.942$ | $2.7^{+1.8}_{-2.7}$ | $2.910$ | $3.0^{+2.0}_{-2.9}$
$\tilde{c}_{\epsilon,2}/1000$ | $-0.640$ | $-0.8\pm 2.4$ | $-0.659$ | $-0.5\pm 2.5$
$\Omega_{m}$ | $0.300$ | $0.304^{+0.012}_{-0.013}$ | $0.306$ | $0.304^{+0.012}_{-0.015}$
$\sigma_{8}$ | $0.830$ | $0.806\pm 0.046$ | $0.821$ | $0.804^{+0.052}_{-0.042}$
Table 2: Marginalized maximal likelihood, mean and $1\sigma$ bounds on the
direct and derived parameters in our analysis from the mock high-$z$ NGC
catalogs. The results are for the pure CDM runs with a prior on $\omega_{b}$
taken from Planck (left) or BBN (right). The units are the same as for Table
1.

The marginalized posterior for $h$, $\ln(10^{10}A_{s})$, and $\omega_{d}$ are
given in Fig. 4(a) for two values of the prior on $\omega_{b}$ and with the
inclusion or exclusion of axions. The full results for all parameters
including the derived $\Omega_{m}$ and $\sigma_{8}$ are given in Table 2 for
the axion-free runs. The results for the latter tests using the BBN prior are
compatible with analyses run for the PyBird code [69]. The main difference we
have observed is that our estimated value of $A_{s}$ is unbiased contrary to
the $2.3\sigma$ tension found in the original PyBird analysis. We attribute
this to the only difference which is in our treatment of the stochastic
parameters (see Appendix A). We specifically include stochastic contributions
to the monopole and find that the inclusion of a stochastic term for the
quadrupole does not affect our results. We note finally that we do not observe
the aforementioned bias in our analysis of the BOSS DR12 data using
$\Lambda$CDM.

###  5.2 Peak-Patch Mocks

(a)

(b)

(c)

(d)

Figure 5: Analysis of the Peak-Patch mock galaxy catalogs where the CDM
scenario denotes an axion-free cosmology and where CDM+Axions denotes a 5%
concentration of $10^{-28}$ eV axions. Such high axion fraction is disfavoured
by observations and is meant to serve as an illustrative case. (a) Simulated
galaxy monopoles with the BOSS simulated catalogs results used for HOD
calibration, (b) galaxy monopoles when removing the contribution of the linear
bias showing a closer similarity with the linear results shown in Fig. 2 with
the error bars removed for clarity, (c) same as (a) for the quadrupoles, (d)
halo mass functions of the halos in the simulation boxes with the Tinker
result used for the CDM calibration. Figure 6: Posterior distributions of the
cosmological parameters, including the axion density for the Peak-Patch mock
galaxy catalogs. The covariance matrix is taken from the MultiDark-Patchy mock
simulations and scaled by 16 with a Planck prior on
$\\{\omega_{d},\omega_{b},\ln(10^{10}A_{s}),h\\}$. The CDM+Axions scenario
denotes simulations containing 5% of $10^{-28}$ eV axions.

We generate a set of 88 (1024 Mpc)3 boxes at a grid resolution of 1 Mpc. The
boxes have a mean redshift of $z\approx 0.57$, and half of them are based on a
$\Lambda$CDM cosmology with parameters $\\{h=0.67,\ln
10^{10}A_{s}=3.091,\omega_{b}=0.02222,\omega_{d}=0.119\\}$ and the other half
with the same cosmological parameters but with a non-zero axion concentration
of $\Omega_{a}/\Omega_{d}=0.05$ for an axion mass of $m_{a}=10^{-28}$ eV. We
use the large-scale structure simulation code Peak-Patch [68] which finds
collapsed halos from a density field on a uniform grid and displaces them to
their final positions using second order Lagrangian perturbation theory. The
algorithm has been shown to reproduce very closely the results of $N$-body
simulations. To generate realistic halo catalogs, we only allow the CDM
component to cluster into halos given that the Jeans scale of axions at
$10^{-28}$ eV forbids them to cluster even at later times. This lack of
clustering combined with the suppression of power on most scales in the
simulation lead to a strong deficit in very massive halos. This is again
mirrored by results involving massive neutrinos where $N$-body simulations
show the same type of halo mass function (HMF) difference as that obtained in
our simulations [70, 71]. The HMF obtained in the presence of ultralight
axions is shown in Fig. 5(d) also agrees with the analytical result of Ref.
[42] for the same axion mass.

We then populate the halos with galaxies using a halo occupation distribution
(HOD) technique. We use a similar prescription to that employed in Refs. [72,
73] where the expectation values of the number of central and satellite
galaxies for a halo of mass $M$ are given by

$\displaystyle\left\langle N_{\text{cen }}\right\rangle_{M}$
$\displaystyle=\frac{1}{2}\left[1+\operatorname{erf}\left(\frac{\log M-\log
M_{\min}}{\sigma_{\log M}}\right)\right],$ (5.1) $\displaystyle\left\langle
N_{\text{sat }}\right\rangle_{M}$ $\displaystyle=\left\langle N_{\text{cen
}}\right\rangle_{M}\left(\frac{M}{M_{\text{sat
}}}\right)^{\alpha}\exp\left(\frac{-M_{\text{cut }}}{M}\right).$ (5.2)

We set $M_{\mathrm{min}}=10^{13.09}\;M_{\odot}/h$, $\sigma_{\log M}=0.19$,
$M_{\mathrm{cut}}=10^{14.2}\;M_{\odot}/h$,
$M_{\mathrm{sat}}=10^{12.9}\;M_{\odot}/h$, $\alpha=0.25$. These values are
obtained by computing the galaxy multipoles for the HOD parameters over a
6-dimensional grid near the values found in the literature. The final
parameters are chosen to minimize the $\chi^{2}$ between the BOSS mock galaxy
catalog multipoles and the average of the CDM-only simulations. For the BOSS
mocks, we use the CMASS sample given the redshift of our simulated boxes.

As our goal is to test our model in a simulated setting, we do not reproduce
all of the characteristics of the BOSS survey. For instance, we do not account
for the limited survey window. To remain self-consistent, we do not apply a
window function when running our MCMC analysis with this dataset. The galaxy
multipoles for the simulated boxes are shown in the two top panels of Fig 5.
We can see that the approximate fitting procedure results in a very close
monopole to the BOSS simulations. The quadrupole is not as well reproduced but
still gives a reasonable approximation on large scales. The small-scale
discrepancy may be due to attribution of velocities to satellite galaxies. The
amplitude of the monopoles initially seems to contradict the linear
predictions, but Fig. 5(b) shows that the expected suppression of power due to
axions is recovered when taking into account the biasing effects of galaxies.
The value of the bias used is the best-fit value from our statistical analysis
of the galaxy multipoles with our EFTofLSS model.

Given that our number of simulations is much smaller than the $\sim 1000$ of
the BOSS simulated dataset, we use the MultiDark-Patchy mocks covariance
matrix scaled down by a factor of 16. We also limit the maximal scale to
$k=0.2\;h/$Mpc so that we have a reliable modelling of the quadrupole. The
main impediment in designing these simulations is that mock galaxy catalogs
are calibrated on $N$-body simulations with galaxy formation. In our case,
such a calibration was not possible since a mixed axion dark matter simulation
suite with baryons is not yet available. We assume that the HOD parameters are
unchanged from the pure CDM case when including a 5% component of ultralight
axions. The study of the impact of subdominant ultralight axionic components
on galaxy formation will be left for future studies once the computational
tools have been developed. For this reason, we prefer to keep the error margin
on the monopole and quadrupole moments for these simulations high, to use
these results conservatively, and to impose Planck-like priors to all
cosmological parameters to isolate the effects of axions.

A supplementary application of our simulations is to investigate what kind of
biases the existence of axions would bring in the measurements of cosmological
parameters. Since we impose strong priors on the varied parameters and since
the covariance is large, we do not observe any bias. However, derived and
nuisance parameters may be biased even with CMB priors.

Figure 7: Marginalised posterior distributions of the $S_{8}$ parameter
measured for a simulated universe with 5% of $10^{-28}$ eV axions. We observe
that assuming that dark matter is purely CDM in a universe containing axions
leads to a biased measurement by many standard deviations. The $1\sigma$
intervals for each of the posteriors are highlighted by shaded areas with the
true value of $S_{8}$ denoted by the dashed line. $S_{8}$ is very sensitive to
power spectrum suppression making it a an ideal probe of ultralight axions.

To illustrate this, in Fig. 7, we show the posterior distribution of the
$S_{8}\equiv\sigma_{8}(\Omega_{m}/0.3)^{0.5}$ parameter which is a measure of
large-scale structure clustering. We observe that the distribution for the
total matter density ($\Omega_{m}$) is unchanged due to strong priors, but
obtain a biased estimate on the derived $S_{8}$. We find that using a model
assuming no axions in a universe where they exist as a subdominant component
leads to biased measurements. The model accounts for the lack of clustering
(lower $S_{8}$) by boosting the value of the linear galaxy bias. Also, a model
with axions biases $S_{8}$ in galaxies high compared to the true value if you
analyse assuming CDM. However, the effect is driven by the galaxy bias, which
causes axions to increase the galaxy power spectrum relative to CDM. This
enhancement is consistent with the same effect seen in the HI halo model [42].
Given our large covariance matrix, the $\chi^{2}$ difference between the two
models is only about a factor of 2, but this is sufficient to obtain a lower
Bayes information criterion (BIC) [74] despite the addition of a free
parameter. However the difference in BIC is too small to definitively favour
one model over the other in this context. Better calibrated simulations with a
more realistic covariance may shed more light on model selection in the
presence on ultralight axions.

##  6 Results

###  6.1 Galaxy Surveys Only

(a)

(b)

Figure 8: Maximum likelihood curve of the monopole and quadrupole for the (a)
LOWZ and (b) CMASS datasets. The model displayed has an axion mass of
$m_{a}=10^{-27}$ eV. The bottom panels show the residuals. Figure 9: Maximum
likelihood (M.L.) monopole for $m_{a}=10^{-27}$ eV. The full line denotes the
overall maximum likelihood curve identical to the CMASS result of Fig. 8. The
dashed line denotes the maximum likelihood imposing that axions compose at
least 5% of the dark matter. We note that this value is about double the 95%
C.L. upper bound for ultralight axions from CMB data and is used for
illustration purposes. We also highlight the axion Jeans scale as it is the
point at which the model begins to deviate from the data. Even if the linear
spectra do match on small scales, the galaxy power spectra do not converge for
$k<k_{J,\mathrm{eq}}$ since the values of the nuisance parameters are
different for both fits.

Using the model described in Section 3, we fit the monopoles and quadrupoles
of the CMASS and LOWZ datasets simultaneously. For each run, we allow the
axion density to vary but fix the axion mass. The best-fit result for the case
with $m_{a}=10^{-27}$ eV is plotted along with the full dataset and the
residuals in Fig. 8. From the fits, we obtain an upper bound on the axion
density at the 68% and 95% confidence levels. We give an example of a value of
the axion fraction in tension with the data in Fig. 9. In this case, the axion
fraction is set to 5% and the model deviates from observations at scales below
the Jeans scale of the axions. This illustrates how the linear power spectrum
suppression (see Fig. 2) leads to detectable effects in the galaxy clustering
multipoles allowing us to extract information about the axion density from
these data.

The combined constraints as a function of mass are displayed in Fig. 1 where
the allowed regions are highlighted in shades of red. Over-plotted is the CMB
analysis for ultralight axions for the Planck 2015 data release taken from
Ref. [28]. We first note that the constraints obtained from the galaxy surveys
are entirely independent on the CMB data with the exception that the power
spectrum tilt has been fixed to the Planck [75] best-fit value of
$n_{s}=0.9611$ throughout this work with the exception of the analysis of
simulations which have been run with $n_{s}=0.9655$. This matches the approach
taken by other BOSS data analyses using similar models [57]. Also, the sum of
the neutrino masses has been fixed to $\sum m_{\nu}=0.06$ eV as in part of the
analysis conducted in Ref. [28]. The prior for the baryon density was taken
from BBN and the analysis in this section can therefore serve as a
verification of CMB studies.

(a)

(b)

Figure 10: Posterior distributions on the axion fraction for the range of
axion masses studied. We note having much better constraint of the heavier
axions $m_{a}\gtrsim 10^{-28}$ eV with the exception of $m_{a}=10^{-25}$ eV
where $k_{J,\mathrm{eq}}\approx k_{\mathrm{max}}$. The numerical values of the
95% C.L. bounds from these marginalized posterior distributions are listed in
Table 3. We note that only the axions with mass $m_{a}\gtrsim 10^{-27}$ eV on
the right hand side are considered to be dark-matter-like since $m_{a}\lesssim
H(a_{\mathrm{osc}})$.

For masses below $10^{-28}$ eV, this study provides an important independent
verification of the bounds on $\Omega_{a}$ found in previous work. Our
$1\sigma$ exclusion bounds match those obtained from CMB very closely. In the
heavier mass range, the constraints we obtain are improved by up to a factor
of 3.7 relative to the CMB-only constraints. We show the full posterior
distributions for the axion fraction in Fig. 10. There we see more clearly the
difference between the highest and lowest masses considered and we observe the
increase in variance for the masses at the lower end. This is due to the
choice of range in scales of galaxy surveys and the axion mass range. Given
any choice of lower $k$ for the dataset, all axion masses for which
$k_{m}(m_{a})\ll k_{\mathrm{min}}$ will be indistinguishable from a simple
shift in amplitude in the power spectrum. The scales below the Jeans scale
will be suppressed by an amount proportional to the axion fraction. However,
if all scales probed by the galaxy survey fall far below this scale, then all
modes are suppressed by the same value. In this setting the axion fraction is
degenerate with the primordial power spectrum amplitude ($A_{s}$).

Figure 11: Joint posterior distribution for the power spectrum amplitude, the
total matter density, and the axion density for four axion masses highlighting
the dependence of the degeneracies on the particle mass.

This can be clearly seen in Fig. 11 where we show the joint posterior
distribution for the axion fraction and the power spectrum amplitude. We note
that for a mass of $10^{-27}$ eV with
$k_{J,\mathrm{eq}}(m_{a})/k_{\mathrm{min}}\approx 1.7$, the parameters appear
independent while in the case with a mass of $10^{-30}$ eV and
$k_{m}(m_{a})/k_{\mathrm{min}}<0.2$, we see a very strong degeneracy between
the two parameters. Also in Fig. 10, we observe a slight preference for a non-
zero axion fraction in the marginalized posterior. This is caused by opening
the $A_{s}$ degeneracy and marginalising over $\Omega_{d}$ and it was also
observed in Ref. [28]. We note that the axion fraction for a $10^{-32}$ eV
axion is not degenerate with either $\omega_{m}$ or $A_{s}$ which explains why
a peak in the posterior distribution is not observed for that mass in Fig. 10.

Figure 12: Area of joint posterior ellipses for three pairs of parameters. The
areas are rescaled by the average taken over axion masses to allow for a
better comparison. We note a simultaneous weakening of the constraints on
$\Omega_{a}/\Omega_{d},\;\ln(10^{10}A_{s}),\;\omega_{m}$ at a mass of
$10^{-30}$ eV.

Heavier axion masses ($10^{-27}-10^{-28}$ eV) also suffer from another
degeneracy with the total matter density as shown in Fig. 11. This is similar
to the degeneracy observed for massive neutrinos as the power spectrum below a
characteristic scale is partially suppressed. In the case of axions, this
suppression as a function of axion concentration and axion mass is shown in
Fig. 2 while for a single massive neutrino species, the suppression can be
approximated by [76]

$\displaystyle\frac{\Delta
P(k)}{P(k)}\approx-0.08\omega_{m}^{-1}\left(\frac{m_{\nu}}{1\;\mathrm{eV}}\right)\;\;\;\mathrm{for}\;k>k_{\mathrm{nr}},$
(6.1)

where $m_{\nu}$ is the neutrino mass and where $k_{\mathrm{nr}}$ is the
horizon scale at the time the neutrinos become non-relativistic which can be
expressed as

$\displaystyle k_{\mathrm{nr}}\approx
0.026\omega_{m}^{1/2}\left(\frac{m_{\nu}}{1\;\mathrm{eV}}\right)^{1/2}\;\mathrm{Mpc}^{-1}.$
(6.2)

Although there exists no analytical approximation of the suppression level of
the power spectrum due to axions there exists a few semi-analytic descriptions
of Eq. (6.1) which can be found in Ref. [44, 77, 78]. The similarities between
massive neutrinos and axions have been the subject of extensive studies using
a collection of other cosmological probes [78]. We note that the axions’ Jeans
scale and the neutrinos $k_{\mathrm{nr}}$ are identical in their particle mass
dependence. Also, we can observe from the BOSS data that axions seem to suffer
from a similar relationship between the total matter density and the axion
concentration where an increase in the matter density can (partially)
compensate for the loss of clustering on small scales. This is evident in the
case of neutrinos when taking the limit $\omega_{m}\gg m_{\nu}/\mathrm{eV}$
and for axions this can be read from Fig. 11. Finally, it is worth noting that
neutrinos are also degenerate with the power spectrum amplitude when studying
the Lyman-$\alpha$ forest [79] as are the lightest axions in our study. It is
worth noting that allowing the sum of the neutrino masses to vary may also
weaken some of the constraints in this work for some axion masses for which
$k_{m}\sim k_{\mathrm{nr}}$.

Fig. 12 summarizes the degeneracies between the axion fraction and the
cosmological parameters which have been varied in this study. In order to
quantify the impact of varying a given parameter on our constraints on the
axion fraction, we use the ellipse area given by [80]

$\displaystyle A(X,Y)=\pi\sigma_{X}\sigma_{Y}\sqrt{1-\rho(X,Y)^{2}}$ (6.3)

where
$X\in\left\\{\Omega_{a}/\Omega_{d},\ln(10^{10}A_{s}),\omega_{m}\right\\}$,
$\sigma_{i}$ is the 1$\sigma$ error on the parameter $i$, and where $\rho$ is
the usual correlation coefficient

$\displaystyle\rho(X,Y)=\frac{\sigma_{XY}}{\sigma_{X}\sigma_{Y}}.$ (6.4)

This derived quantity captures degeneracies with axions but it also accounts
for the capacity of the data to constrain the parameters correlated with
axions. This is justified since a degeneracy with a very well constrained
parameter will not lead to a wider posterior distribution for the axion
fraction. We note a strong increase in area for axions with a mass of
$m_{a}\sim 10^{-30}$ eV for the ellipse taken from the joint posterior
distribution of the matter density and the power spectrum amplitude. This
indicates a degeneracy which loosens the constraints on the axion fraction at
that mass. This can be seen in Fig. 1 as a wider permitted region at
$10^{-30}$ eV. This explains why the one and two sigma contours are not
monotonically decreasing from $10^{-32}$ eV as one would naively expect.

Finally, we look for degeneracies between the axion fraction and the nuisance
parameters of the model. We show the marginalised posterior distribution of
the nuisance parameters along with the axion fraction in Fig. 13. Given the
shape of the galaxy monopoles in Fig. 3, we expected to find strong
degeneracies between the axion fraction and the linear galaxy bias $b_{1}$.
However, we do not observe such degeneracy for either of the redshift slices
in the data. Another study based on a similar model finds no strong degeneracy
between the galaxy biases and the power spectrum amplitude [59] which is
itself degenerate with the axion fraction for low axion masses ($\leq
10^{-29}$ eV). We conclude that our choice of model based on the EFTofLSS was
well suited for probing ultralight axions using galaxy clustering despite the
high number of nuisance parameters associated with working with such datasets.

Figure 13: Marginalized 2D posterior distributions for the model’s nuisance
parameters (see Table 1) and the axion fraction (leftmost column) for an axion
mass of $10^{-27}$ eV. We note no severe degeneracies between the axion
fraction and the nuisance parameters for either low or high redshifts. Similar
results are obtained for all other axion masses in our analysis.

###  6.2 Galaxy Surveys with a CMB Prior

We now recalculate the constraints obtained in the previous section by
replacing the BBN prior on $\omega_{b}$ with a CMB prior on all cosmological
parameters and the axion fraction. The CMB prior is obtained from the chains
from the analysis in Ref. [28]. This work obtained constraints on axions using
a combination of temperature, polarization and lensing data from the CMB. It
also considered adiabatic initial conditions and repeated the analysis
including axion-sourced isocurvature perturbations. In this work, we use the
former as our prior since we assume adiabatic initial conditions for the
density perturbations. The chains are once again run for each mass bin between
$10^{-25}\;\mathrm{eV}$ and $10^{-32}\;\mathrm{eV}$. The resulting updated
constraint plot is shown in Fig. 1.

Figure 14: Joint posterior distributions for an axion with a mass of $10^{-32}$ eV for three experimental setups. We note an improvement on the constraint on the axion fraction when breaking the degeneracy with $H_{0}$ present with the CMB data. The gray shaded area represent the confidence interval for $h$ from the SH0ES measurement [81]. Figure 15: 68% (dark-colored) and 95% (light-colored) confidence level bounds on the axion density from the CMB data, galaxy clustering and the combined measurements. Mass [$\log(m_{a}/\mathrm{eV})$] | CMB (Planck 2015) | Galaxy Clustering | Combined Measurements
---|---|---|---
$-25$ | $<0.101$ | $<0.022$ | $<0.041$
$-26$ | $<0.006$ | $<0.005$ | $<0.004$
$-27$ | $<0.003$ | $<0.004$ | $<0.002$
$-28$ | $<0.003$ | $<0.005$ | $<0.002$
$-29$ | $<0.003$ | $<0.010$ | $<0.002$
$-30$ | $<0.004$ | $<0.011$ | $<0.002$
$-31$ | $<0.006$ | $<0.011$ | $<0.003$
$-32$ | $<0.015$ | $<0.011$ | $<0.007$
Table 3:  Numerical values of the 95% C.L. upper bounds on $\Omega_{a}h^{2}$
from three different experimental configurations as shown in Fig. 15.

The joint analysis allows for much tighter constraints on the axion fraction
by breaking degeneracies between $\Omega_{a}$ and other cosmological
parameters. See Fig. 14 for an example of broken degeneracy in the case of
$H_{0}$ with $10^{-32}$ eV axions. When considering only galaxy clustering
data, the presence of axions increases the error bars on $H_{0}$ thus
initially alleviating the tension with late time measurements. This is very
similar to the impact of the inclusion of a higher effective number of
relativistic species of Ref. [82]. When including the CMB prior however, we
observe that a large fraction of axions at a mass of $10^{-32}$ eV leads to an
even higher discrepancy with SN1a measurements of $H_{0}$. This contrasts with
other ultralight scalars with modified field potentials which have been
proposed as solutions to this tension [83]. There is one case however where
the addition of a CMB prior does not lead to stronger constraints on the axion
fraction.

Figure 16: Degeneracy between the power spectrum amplitude and the axion
fraction for an axion mass of $10^{-25}$ eV. We observe that the CMB prior
favours a region of parameter space for which the galaxy clustering data has a
wider uncertainty on the axion fraction.

Indeed, the constraint on the axion fraction with the mass of $10^{25}$ eV
worsens when adding a CMB prior. This is due to the fact that the CMB prior
favours a higher value of $A_{s}$ which is slightly degenerate with the axion
fraction at that mass as shown in Fig. 16. Another contributing factor is that
the CMB prior does not constrain the axion fraction as well as for the axion
masses below $10^{-25}$ eV. Performing a joint likelihood analysis rather than
imposing a prior on the cosmological parameters may allow for stronger
constraints for this mass bin and is left for future work. We note however
that galaxy clustering measurements alone improve existing constraints on the
axion fraction at that mass by over a factor of 4.5 (see Table 3).

Combining CMB data with galaxy clustering allows us to overcome the
degeneracies present in the CMB analysis. The improvements on the bounds on
the axion density are displayed in Fig. 15. There we see a factor of two
improvement on CMB alone for most of the axion masses considered. A joint
analysis of CMB and galaxy survey data was conducted in Ref [40] with the
Planck and the WiggleZ [84] surveys and gave $\Omega_{a}/\Omega_{d}\leq 0.048$
for a similar mass range. This work was based on the 2013 data release and
considered adiabatic initial conditions. Two factors contribute to improved
constraints in the present work. The first is that the WiggleZ data were
limited to the galaxy power spectrum while the BOSS data used includes both
the monopoles and quadrupoles. This allows us to detect the axions’
anisotropic effects on galaxies and amplifies the size of our galaxy
clustering dataset by about a factor of four. The second factor is the
development of modeling techniques allowing us to model the power spectrum
multipoles up to a scale of $k=0.23\;h/$Mpc while the initial work of Ref.
[40] was limited to $k_{\mathrm{max}}=0.2\;h/$Mpc. This difference is most
notable on the heaviest mass of $10^{-25}$ eV where our 95% C.L. upper bound
is improved by a factor of $\sim 1.6$ since the power spectrum cut-off
manifests on the smallest scales.

##  7 Discussion

We introduced ultralight axion physics in a galaxy clustering model by
modifying the linear power spectrum and then calculating higher order
corrections using the EFTofLSS approach through the use of the public code
PyBird. To increase calculation efficiency, we also developed a way to obtain
the linear matter power spectrum in the presence of ultralight axions using
interpolations thus significantly decreasing the total computation time. We
first ran the canonical tests for galaxy modeling using the MultiDark-Patchy
mocks simulation suite and verified that our model was able to recover
$\Lambda$CDM parameters in the absence of axions. We also generated our own
mock catalogs using an HOD approach and the Lagrangian based simulation code
Peak-Patch. This allowed us to assert that our model is sensitive to even
small axion fractions of order $\Omega_{a}/\Omega_{d}\sim 0.05$. Results from
simulations including $10^{-28}$ eV axions also showed that they could bias
our estimates of $S_{8}$ if they are present but unaccounted for in our model.
In the present work, we have not considered interactions between axions and
the standard model. A typical interaction takes the form $\phi{\bf E}\cdot{\bf
B}$. Such an interaction affects the polarisation of the CMB, and for
$10^{-33}\text{ eV}\leq m_{a}\leq 10^{-28}\text{ eV}$ could explain the
recently observed _birefringence_ in CMB anisotropies [85]. If this
measurement is confirmed, this motivates a drive to improve searches for
axions in this mass window using structure formation, which could allow an
independent measurement of $m_{a}$ and $f_{a}$ impossible with birefringence
alone. We note that birefringence may also be induced by a network cosmic
strings if the Peccei-Quinn symmetry of the axion field is broken after
inflation [86, 87].

Having tested our model in various axion and axion-free cosmologies, we ran
multiple studies on the BOSS galaxy clustering dataset in redshift space. We
performed 16 different runs using eight logarithmically spaced mass bins and
two choices of priors. Our axion mass range extended over seven orders of
magnitude from $10^{-32}$ eV to $10^{-25}$ eV. Axions lighter than this would
behave as a dark energy component in the relevant regime. Axions heavier than
the upper bound would have a Jeans scale much smaller than the semi-linear
scales to which we are sensitive making them identical to CDM for $k\leq
k_{\mathrm{max}}$. Constraining axions heavier than this mass scale would
require a fully non-linear treatment which has not yet been fully explored in
the mixed axion dark matter case. See Ref. [88] for a semi-linear treatment at
high redshift for initial conditions and Ref. [89] for the spherical collapse
in mixed dark matter case. Using a BBN prior on the baryon density only, we
arrived at constraints matching that of the CMB analysis for the lighter
masses, but improved them significantly for masses above $10^{-29}$ eV. Using
a CMB prior on all cosmological parameters including the axion parameters, we
arrived at the strongest constraint on the axion density with
$\Omega_{a}h^{2}<0.022$ for the heaviest axions in this study. We also arrived
at the bounds $\Omega_{a}h^{2}<0.002$ for the mass window
$10^{-30}\;\mathrm{eV}\leq m_{a}\leq 10^{-27}\;\mathrm{eV}$ and
$\Omega_{a}h^{2}<0.007$ for axions with $m_{a}=10^{-32}\;\mathrm{eV}$. These
results are expressed in terms of the axion fraction in Fig. 1 and in terms of
the axion relic density in Table 3 and Fig. 15. Combining CMB and large-scale
structure data, we were able to break important degeneracies between axions
and $\Lambda$CDM parameters and we arrived at the most stringent bounds on the
axion density for this mass range to date. The 95% C.L. bounds on the axion
fraction obtained in this work approach those projected for the CMB-S4
observatory [90] alone but were attained using existing data thus showcasing
the power of multi-probe approaches in cosmology.

## Acknowledgments

We would like to thank Mona Dentler, Colin Hill, Evan McDonough, and Mikhail
Ivanov for useful comments and suggestions. RH is a CIFAR Azrieli Global
Scholar, Gravity & the Extreme Universe Program, 2019, and a 2020 Alfred. P.
Sloan Research Fellow. RH is supported by Natural Sciences and Engineering
Research Council of Canada. RB is a CIFAR Fellow. AL, RB and RH are supported
by Natural Sciences and Engineering Research Council of Canada. The Dunlap
Institute is funded through an endowment established by the David Dunlap
family and the University of Toronto. We acknowledge that the land on which
the University of Toronto is built is the traditional territory of the
Haudenosaunee, and most recently, the territory of the Mississaugas of the New
Credit First Nation. We are grateful to have the opportunity to work in the
community, on this territory. Computations were performed on the Niagara
supercomputer [91, 92] at the SciNet HPC Consortium. SciNet is funded by: the
Canada Foundation for Innovation; the Government of Ontario; Ontario Research
Fund - Research Excellence; and the University of Toronto. The data
visualization in this work used the software package getdist [93]. The MCMC
chains were run in parallel with MPI through the use of the schwimmbad library
[94].

##  Appendix A More on Perturbation Theory

###  A.1 Higher Order Terms with CDM Dynamics

The first loop correction has two contributions and takes the form of [95]

$\displaystyle P_{1-\mathrm{loop}}(k)=P_{22}(k)+P_{13}(k),$ (A.1)

where

$\displaystyle P_{22}(k)$ $\displaystyle\equiv
2\int\frac{d^{3}q}{(2\pi)^{3}}\left[F_{2}(\mathbf{q},\mathbf{k}-\mathbf{q})\right]^{2}P_{\mathrm{lin}}(|\mathbf{k}-\mathbf{q}|)P_{\mathrm{lin}}(q),$
(A.2) $\displaystyle P_{13}(k)$ $\displaystyle\equiv
6\int\frac{d^{3}q}{(2\pi)^{3}}F_{3}(\mathbf{q},-\mathbf{q},\mathbf{k})P_{\mathrm{lin}}(k)P_{\mathrm{lin}}(q),$
(A.3)

where $P_{\mathrm{lin}}(k)$ is the linear matter power spectrum555The linear
matter power spectrum is often denoted as $P_{11}(k)$. and where the functions
$F_{2,3}$ are standard Eulerian perturbation theory kernels which can be found
by recursion relations [95]. To model the data, we must also consider the
effects redshift space distortions and the fact that galaxies are biased
tracers. We can model redshift space distortions by mapping the redshift space
positions $\mathbf{x}_{r}$ from the real-space positions following

$\displaystyle\mathbf{x}_{r}=\mathbf{x}+\frac{\hat{\mathbf{z}}\cdot\hat{\mathbf{v}}}{\mathcal{H}}\hat{\mathbf{z}},$
(A.4)

where $\mathbf{v}$ is the tracer’s velocity. The impact of redshift space
distortions on biased tracers has been studied in detail in Ref. [96]. The
relationship between the halo overdensities in redshift space are then
obtained from the real space results by

$\displaystyle\delta_{h,r}(\mathbf{k})=\delta_{h}(\mathbf{k})+\int
d^{3}xe^{-i\mathbf{k}\cdot\mathbf{x}}\left[\exp\left(-i\frac{k_{z}}{\mathcal{H}}v_{h,z}\right)-1\right](1+\delta_{h}(\mathbf{x})).$
(A.5)

Expanding the halo overdensity up to third order, we get

$\displaystyle\delta_{h,r}=\delta_{h,r}^{(1)}+\delta_{h,r}^{(2)}+\delta_{h,r}^{(3)}+\delta^{(3)}_{h,r,\mathrm{counter}}+\delta^{(3)}_{h,r,\mathrm{stoch}}+....$
(A.6)

where the two last terms account for the counter and stochastic terms arising
from the renormalization of the power spectrum (need more details). The halo
power spectrum in redshift space is computed with loop-order corrections from
the linear matter power spectrum by the inclusion of higher order terms and
with the addition of counter and stochastic corrections. The full expression
from Ref. [54] reads

$\displaystyle\left\langle\delta_{h,r}(\mathbf{k})\delta_{h,r}(\mathbf{k})\right\rangle$
$\displaystyle=\left\langle\delta_{h,r}^{(1)}\delta_{h,r}^{(1)}\right\rangle+\left\langle\delta_{h,r}^{(2)}\delta_{h,r}^{(2)}\right\rangle+2\left\langle\delta_{h,r}^{(1)}\delta_{h,r}^{(3)}\right\rangle+\left\langle\delta_{h,r}\delta_{h,r}\right\rangle_{\mathrm{ct}}+\left\langle\delta_{h,r}\delta_{h,r}\right\rangle_{\epsilon}$
(A.7) $\displaystyle=b_{1}^{2}P_{\mathrm{lin}}(k)+2\int
d^{3}q\left[K_{h,r}^{(2)}(b_{1},b_{2},b_{4};\mathbf{q},\mathbf{k}-\mathbf{q})_{\mathrm{sym}}\right]^{2}P_{\mathrm{lin}}(|\mathbf{k}-\mathbf{q}|)P_{\mathrm{lin}}(q)$
$\displaystyle+6b_{1}\int
d^{3}qK_{h,r}^{(3)}(b_{1},b_{3};\mathbf{q},-\mathbf{q},\mathbf{k})_{\mathrm{sym}}P_{\mathrm{lin}}(q)P_{\mathrm{lin}}(k)+\left\langle\delta_{h,r}\delta_{h,r}\right\rangle_{\mathrm{ct}}+\left\langle\delta_{h,r}\delta_{h,r}\right\rangle_{\epsilon},$
(A.8)

where $b_{1,2,3,4}$ are the bias parameters which will be varied in our
analysis and $K_{h,r}^{(2,3)}$ are the symmetrized versions of the redshift-
space halo kernels whose complete derivation can be found in Ref. [54]. The
redshift-space kernels are

$\displaystyle\begin{split}K_{h,r}^{(2)}\left(\mathbf{q}_{1},\mathbf{q}_{2}\right)=K_{\delta_{h}}^{(2)}&\left(\mathbf{q}_{1},\mathbf{q}_{2}\right)+f\mu^{2}K_{\theta_{h}}^{(2)}\left(\mathbf{q}_{1},\mathbf{q}_{2}\right)\\\
&+\frac{1}{2}b_{1}\mu
f\left(\frac{kq_{2z}}{q_{2}^{2}}+\frac{kq_{1z}}{q_{1}^{2}}\right)K_{\theta_{h}}^{(1)}\left(\mathbf{q}_{1}\right)+\frac{1}{2}\mu^{2}f^{2}\frac{k^{2}q_{1z}q_{2z}}{q_{1}^{2}q_{2}^{2}}K_{\theta_{h}}^{(1)}\left(\mathbf{q}_{1}\right)K_{\theta_{h}}^{(1)}\left(\mathbf{q}_{2}\right)\end{split}$
(A.9)
$\displaystyle\begin{split}K_{h,r}^{(3)}\left(\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{q}_{3}\right)=&K_{\delta_{h}}^{(3)}\left(\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{q}_{3}\right)+f\mu^{2}K_{\theta_{h}}^{(3)}\left(\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{q}_{3}\right)\\\
&+b_{1}\mu
f\left(\frac{kq_{3z}}{q_{3}^{2}}\right)K_{\delta_{h}}^{(2)}\left(\mathbf{q}_{1},\mathbf{q}_{2}\right)K_{\theta_{h}}^{(1)}\left(\mathbf{q}_{3}\right)+\mu
f\left(\frac{k\left(q_{1z}+q_{2z}\right)}{\left(\mathbf{q}_{1}+\mathbf{q}_{2}\right)^{2}}\right)K_{\theta_{h}}^{(2)}\left(\mathbf{q}_{1},\mathbf{q}_{2}\right)\\\
&+\frac{1}{2}\mu^{2}f^{2}\left(\frac{kq_{1z}}{q_{1}^{2}}\frac{kq_{2z}}{q_{2}^{2}}\right)K_{\theta_{h}}^{(1)}\left(\mathbf{q}_{1}\right)K_{\theta_{h}}^{(1)}\left(\mathbf{q}_{2}\right)K_{\delta_{h}}^{(1)}\left(\mathbf{q}_{3}\right)\\\
&+\mu^{2}f^{2}\left(\frac{k\left(q_{1z}+q_{2z}\right)}{\left(\mathbf{q}_{1}+\mathbf{q}_{2}\right)^{2}}\frac{kq_{3z}}{q_{3}^{2}}\right)K_{\theta_{h}}^{(2)}\left(\mathbf{q}_{1},\mathbf{q}_{2}\right)K_{\theta_{h}}^{(1)}\left(\mathbf{q}_{3}\right)\\\
&+\frac{1}{6}\mu^{3}f^{3}\left(\frac{kq_{1z}}{q_{1}^{2}}\frac{kq_{2z}}{q_{2}^{2}}\frac{kq_{3z}}{q_{3}^{2}}\right)K_{\theta_{h}}^{(1)}\left(\mathbf{q}_{1}\right)K_{\theta_{h}}^{(1)}\left(\mathbf{q}_{2}\right)K_{\theta_{h}}^{(1)}\left(\mathbf{q}_{3}\right),\end{split}$
(A.10)

where the dependence on the bias parameters is mae explicit in the higher
order real-space kernels. For the halo overdensity $\delta_{h}$ they are given
by

$\displaystyle\begin{split}K_{\delta_{h}}^{(2)}(k,q,x)_{\mathrm{sym}}=&\frac{b_{1}-2q^{3}+k^{3}x+4kq^{2}x-k^{2}q\left(1+2x^{2}\right)}{k^{2}+q^{2}-2kqx}\\\
&+\frac{b_{2}}{7}\frac{7q^{2}-14kqx+k^{2}\left(5+2x^{2}\right)}{k^{2}+q^{2}-2kqx}+b_{4}\end{split}$
(A.11)
$\displaystyle\begin{split}K_{\delta_{h}}^{(3)}(k,q)_{\mathrm{UV}-\mathrm{sub},\mathrm{sym}}=\frac{b_{1}}{504k^{3}q^{3}}\left(-38k^{5}q+48k^{3}q^{3}-18kq^{5}+9\left(k^{2}-q^{2}\right)^{3}\log\left[\frac{k-q}{k+q}\right]\right)\\\
+\frac{b_{3}}{756k^{3}q^{5}}\left(2kq\left(k^{2}+q^{2}\right)\left(3k^{4}-14k^{2}q^{2}+3q^{4}\right)+3\left(k^{2}-q^{2}\right)^{4}\log\left[\frac{k-q}{k+q}\right]\right).\end{split}$
(A.12)

And for the halos’ velocity divergence $\theta_{h}$, we have

$\displaystyle K_{\theta_{h}}^{(1)}(k,q,x)_{\mathrm{sym}}$ $\displaystyle=1,$
(A.13) $\displaystyle K_{\theta_{h}}^{(2)}(k,q,x)_{\mathrm{sym}}$
$\displaystyle=\frac{k^{2}\left(7kx-q\left(1+6x^{2}\right)\right)}{14q\left(k^{2}+q^{2}-2kqx\right)},$
(A.14) $\displaystyle
K_{\theta_{h}}^{(3)}(k,q,x)_{\mathrm{UV}-\mathrm{sub},\mathrm{sym}}$
$\displaystyle=\frac{12k^{7}q-82k^{5}q^{3}+4k^{3}q^{5}-6kq^{7}+3\left(k^{2}-q^{2}\right)^{3}\left(2k^{2}+q^{2}\right)\log\left[\frac{k-q}{k+q}\right]}{504k^{3}q^{5}},$
(A.15)

where UV-sub denotes the term with the high $k$ modes subtracted. We once
again refer readers to Ref. [54] for more detail. The counterterm contribution
is given by

$\displaystyle\left\langle\delta_{h,r}(\mathbf{k})\delta_{h,r}(\mathbf{k})\right\rangle_{\mathrm{ct}}$
$\displaystyle=2P_{\mathrm{lin}}(k)\left(b_{1}+f\mu^{2}\right)$
$\displaystyle\;\;\;\times\left(\mu^{2}\left(\frac{k}{k_{\mathrm{M}}}\right)^{2}\tilde{c}_{r,1}+\mu^{4}\left(\frac{k}{k_{\mathrm{M}}}\right)^{2}\tilde{c}_{r,2}+c_{\mathrm{ct}}\left(\frac{k}{k_{\mathrm{NL}}}\right)^{2}\right),$
(A.16)

where $f\equiv d\ln D/d\ln a$ is the logarithmic derivative of the linear
growth factor $P(k,t)=\frac{D^{2}(t)}{D^{2}(t_{0})}P(k,t_{0})$, where $k_{m}$
is the scale of halos (more details here), and where $k_{\mathrm{NL}}$ is the
scale at which non-linearities begin to dominate. The cosine of the angle
between the line-of-sight and and the wavevector is denoted $\mu$ and
$\\{\tilde{c}_{r,1},\tilde{c}_{r,2},c_{\mathrm{ct}}\\}$ is the set of
counterterms which will be varied in our analysis.

$\displaystyle\left\langle\delta_{h,r}\delta_{h,r}\right\rangle_{\epsilon}=\frac{1}{\bar{n}_{W}}\left[c_{\epsilon,1}+(c_{\epsilon,2}+c_{\epsilon,2}f\mu^{2})\left(\frac{k}{k_{\mathrm{M}}}\right)^{2}\right],$
(A.17)

where $\bar{n}_{W}$ is the typical halo density which we approximate using
[97]

$\displaystyle\bar{n}_{W}=\left(\int
dM\frac{M^{2}}{\bar{\rho}^{2}}\frac{dn}{dM}\right)^{-1}\approx
330\;h^{3}/\mathrm{Mpc}^{3}.$ (A.18)

Given that $k_{\mathrm{M}}\approx 0.7\;h/\mathrm{Mpc}$ and
$k_{\mathrm{max}}\leq 0.23\;h/\mathrm{Mpc}$, we find little scale dependence
on the stochastic terms and replace them with a single constant shot noise
contribution to the monopole term along with a quadratic high-momentum term
and absorb the contribution of the inverse halo density giving

$\displaystyle\left\langle\delta_{h,r}\delta_{h,r}\right\rangle_{\epsilon}\approx\tilde{c}_{\epsilon,1}+\tilde{c}_{\epsilon,2}\left(\frac{k}{k_{\mathrm{M}}}\right)^{2}.$
(A.19)

And this completes our model. We thus have initially nine free parameters
aside from the axion and cosmological parameters which are
$\\{b_{1,2,3,4},c_{\mathrm{ct}},\tilde{c}_{r,1},\tilde{c}_{r,2},\tilde{c}_{\epsilon,1},\tilde{c}_{\epsilon,2}\\}$.
We identify which of the parameters we fix and the priors used for each of the
varied ones in Section 4.

###  A.2 Wave Effects on Loop Corrections

In this section, we tackle the importance of the dynamical effects of
ultralight axions. Due to their quantum diffusion, axions do not cluster in
the same way as collisionless dark matter and care must be taken if we wish to
use the EFTofLSS which was designed for CDM. The wave effects on the
clustering of axions is encoded at first order in the linear growth factor
which is found by solving [52]

$\displaystyle\ddot{D}(k,a)+2H\dot{D}(k,a)+\left(\frac{\hbar^{2}k^{4}}{4m_{a}^{2}a^{2}}-4\pi
G\bar{\rho}\right)D(k,a)=0,$ (A.20)

where the overdot denotes a derivative with respect to time. For simplicity,
we will assume an Einstein deSitter cosmology from now on. We note that we
recover the scale-independent CDM growth factor $D_{\mathrm{CDM}}(a)\propto a$
as we go to large scales. On small scales, the growth of perturbations is
suppressed and rapidly oscillates around zero. Therefore, the assessment the
importance of wave effects will have to be conducted for each scale.
Integrating the differential equation for the growth factor from an initial
$a_{\mathrm{in}}$, we have the solution [50, 51, 98]

$\displaystyle
D(k,a)=\left(\frac{a_{\mathrm{in}}}{a}\right)^{1/4}\frac{J_{-5/2}\left(\hbar
k^{2}/m_{a}H_{0}\sqrt{a}\right)}{J_{-5/2}\left(\hbar
k^{2}/m_{a}H_{0}\sqrt{a_{\mathrm{in}}}\right)},$ (A.21)

where the Bessel function of the first kind is given by

$\displaystyle J_{-5/2}(x)=\sqrt{\frac{2}{\pi x}}\bigg{(}\frac{3\cos
x}{x^{2}}+\frac{3\sin x}{x}-\cos x\bigg{)}.$ (A.22)

Taking the $k\to 0$ limit, we get a constant growth equal to
$a/a_{\mathrm{in}}$ as expected. Deviations from the CDM growth on small
scales is what alters the kernels of in the integrals for the 1-loop
corrections for the matter power spectrum. However, these deviations only
affect the part of the matter which is composed of axions and it was shown
earlier that the axion perturbations at the linear level are suppressed below
their Jeans scale at matter-radiation equality. We argue here that _the loss
of power in the perturbations on small scales erases the wave effects before
they become dominant._ Let us define the quantum corrections to the growth
factor as

$\displaystyle D_{Q}\equiv D(k,a)-D_{\mathrm{CDM}}(a).$ (A.23)

The final amplitude of total quantum corrections depend on the size of the
axion perturbations which they affect. Therefore, the amplitude of the quantum
corrections at linear order are proportional to the product of the growth
factor and the perturbations at each scale. If we define the relative
corrections normalized with the CDM perturbations we have

$\displaystyle\mathcal{Q}^{(1)}(k,a)=\left(\frac{D_{Q}(k,a)}{D_{\mathrm{CDM}}(a)}\right)^{2}\frac{\left\langle\delta^{(1)}_{a}\delta^{(1)}_{a}\right\rangle}{\left\langle\delta^{(1)}_{c}\delta^{(1)}_{c}\right\rangle},$
(A.24)

where $\delta_{a,c}^{(i)}$ is the $i^{\mathrm{th}}$ order perturbations and
where $a$ and $c$ denote the axion and CDM components respectively.

Figure 17: Amplitude of the wave effects at linear order. We see that the wave
effects begin to dominate in a regime where the axion perturbations are
severely suppressed and that the relative size of the quantum corrections are
bounded from above. This was generated for an axion mass of $10^{-27}$ eV at
redshift $z=1$.

We plot $\mathcal{Q}^{(1)}$ in Fig. 17 where we can see that the sharp cutoff
in the variance of the axion perturbations keeps the amplitude of the total
quantum corrections bounded. Note that we used the fitting formula of Ref [88]
Eq. (23) to remove numerical artifacts when computing the growth factor
directly from Bessel functions.

It remains to show that this reasoning applies to higher order perturbations
e.g. $\delta^{(2)},\;\delta^{(3)}$, etc. This is due to the fact that higher
order perturbations are computed using lower order ones through recursion
relations. Therefore if $\delta^{(1)}_{a}(k)=0$ for $k\gg k_{J,\mathrm{eq}}$,
then so will $\delta_{a}^{(2)}(k)$ and $\delta_{a}^{(3)}(k)$. With regards to
the growth, the second order growth factor can be found by solving the
inhomogeneous equation [95]

$\displaystyle\ddot{D}_{2}(k,a)+2H\dot{D}_{2}(k,a)+\left(\frac{\hbar^{2}k^{4}}{4m_{a}^{2}a^{2}}-4\pi
G\bar{\rho}\right)D_{2}(k,a)=-4\pi G\bar{\rho}D(k,a)^{2}.$ (A.25)

We anticipate $D_{2}$ to behave in a similar way to the linear growth on small
scales ($D_{2}\to 0$) and on large scales ($D_{2}\to D_{2,\mathrm{CDM}}$). We
thus expect that the quantum corrections to the 1-loop power spectrum (which
are the higher order equivalent to $\mathcal{Q}^{(1)}$) to be small as well.
This is corroborated by the analysis of Ref. [98] where it was shown that the
wave effects on higher order corrections become important at scales at which
the 1-loop power spectrum is heavily suppressed. It is worth pointing out that
this study was conducted for a single component dark matter of
$m_{a}=10^{-23}$ eV. In the case of a two-component dark matter, the total
dark matter perturbations ($\delta_{d}$) obey

$\displaystyle\left\langle\delta_{d}(k)\delta_{d}(k)\right\rangle$
$\displaystyle=\left\langle\left[\frac{\Omega_{a}}{\Omega_{d}}\delta_{a}(k)+\left(1-\frac{\Omega_{a}}{\Omega_{d}}\right)\delta_{c}(k)\right]\left[\frac{\Omega_{a}}{\Omega_{d}}\delta_{a}(k)+\left(1-\frac{\Omega_{a}}{\Omega_{d}}\right)\delta_{c}(k)\right]\right\rangle$
(A.26)
$\displaystyle=\left(\frac{\Omega_{a}}{\Omega_{d}}\right)^{2}\left\langle\delta_{a}(k)\delta_{a}(k)\right\rangle+2\left(1-\frac{\Omega_{a}}{\Omega_{d}}\right)\left(\frac{\Omega_{a}}{\Omega_{d}}\right)\left\langle\delta_{a}(k)\delta_{c}(k)\right\rangle$
$\displaystyle\;\;\;\;+\left(1-\frac{\Omega_{a}}{\Omega_{d}}\right)^{2}\left\langle\delta_{c}(k)\delta_{c}(k)\right\rangle.$
(A.27)

Therefore the total one-loop contribution to the dark matter gives

$\displaystyle
P_{\mathrm{1-loop}}^{dd}(k)=\left(\frac{\Omega_{a}}{\Omega_{d}}\right)^{2}P_{\mathrm{1-loop}}^{aa}(k)+2\left(\frac{\Omega_{a}}{\Omega_{d}}\right)\left(1-\frac{\Omega_{a}}{\Omega_{d}}\right)P_{\mathrm{1-loop}}^{ac}(k)+\left(1-\frac{\Omega_{a}}{\Omega_{d}}\right)^{2}P_{\mathrm{1-loop}}^{cc}(k).$
(A.28)

Since only $P^{aa}$ and $P^{ac}$ have contributions from the wave effects, we
see that the amplitude of these effects in the mixed case are multiplied by at
least one power of $\Omega_{a}/\Omega_{d}$ which is less than 0.05 in the vast
majority of the cases in our analysis. We conclude that the corrections from
the wave effects on axion dynamics are even smaller in the mixed case and
vanish in the limit $\Omega_{a}\ll 1$, as expected.

##  Appendix B Axions’ Anisotropic Effects

In the rightmost panel of Fig. 3, we observe an enhancement of the quadrupole
moment despite a suppression of structure for the $\ell=0$ multipole. This is
a unique feature as it indicates possible anisotropic effects of the structure
suppression of axions and constitutes a completely new signature beyond the
well-known structure suppression. To investigate this, we make use of a very
simple redshift space model for the galaxy power spectrum where we approximate
the galaxy power spectrum as

$\displaystyle P_{g}(k,\mu)\approx e^{-(k\mu
f\sigma_{v})^{2}}\left(1+f\mu^{2}\right)^{2}b_{g}^{2}P_{\mathrm{lin}}(k),$
(B.1)

where $b_{g}$ is the galaxy bias and where $\sigma_{v}$ is the galaxy velocity
dispersion. This model is based on the Kaiser approximation [99] with a
Gaussian kernel for the finger-of-God effects. The velocity dispersion can be
roughly approximated at linear order with (see Ref. [100] and references
therein)

$\displaystyle\sigma_{v,\mathrm{lin}}^{2}=\frac{1}{6\pi^{2}}\int
dqP_{\mathrm{lin}}(q).$ (B.2)

Using Eq. (3.2), we have that the multipoles of the power spectrum are

$\displaystyle
P_{\ell}(k)=\frac{2\ell+1}{2}b_{g}^{2}P_{\mathrm{lin}}(k)\underbrace{\int_{-1}^{1}d\mu\;e^{-(k\mu
f\sigma_{v})^{2}}\left(1+f\mu^{2}\right)^{2}\mathcal{P}_{\ell}(\mu)}_{\equiv\mathcal{B}_{\ell}(k;\sigma_{v})}.$
(B.3)

Figure 18: Function $\mathcal{B}_{\ell}$ as defined in Eq. (B.3) as a function
of $k$ for the monopole and quadrupole. We note an enhancement on small scales
when reducing the velocity dispersion (e.g. by suppressing the matter power
spectrum).

From this simple model, we find that the increase in the quadrupole moment is
attributable to a decrease in the velocity divergence which arises when the
linear matter power spectrum is suppressed and the value of the integral in
Eq. (B.2) decreases. This decrease in $\sigma_{v}$ leads to a higher value of
the $B_{\ell}$, especially for the $\ell=2$ as shown in Fig. 18.

##  Appendix C Axion Transfer Function Interpolation

The axion transfer function defined in Eq. 3.3 captures the deviation from
$\Lambda$CDM due to axions in the matter power spectrum. It is most often
obtained through semi-analytic approximations [27] or numerically with adapted
Boltzmann codes. In the present study however, none of those options are
optimal. The approach given in Ref. [27] was only for the case of
$\Omega_{a}/\Omega_{d}=1$ and the axionCAMB has a long runtime of up to five
seconds. Although the code was used extensively for CMB studies and Fisher
forecasts analyses, the higher dimensionality of the current study makes the
high runtime problematic as larger number of samples is required for the MCMC
chains to converge.

To overcome the computation time issue, we use an interpolating approach where
the result of the axionCAMB code is calculated for a predefined grid in
parameter space. The the RegularGridInterpolator of the Scipy
package666https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.RegularGridInterpolator.html
is used to linearly interpolate the axion transfer function for the choice of
parameters given by the MCMC sampler. Two approaches are considered when
building the interpolator. The first is to assume the axion transfer function
depends only on the axion parameters. In other words, at a fixed mass $m_{a}$,
one only need to interpolate over different values of $\Omega_{a}/\Omega_{d}$.
This approximation is valid in the limit where the axion fraction is very
small since

$\displaystyle
T_{\mathrm{ax}}^{2}(\Omega_{a}/\Omega_{d},\mathbf{b}_{\mathrm{cosmo}};k)=\frac{P(\Omega_{a}/\Omega_{d},\mathbf{b}_{\mathrm{cosmo}};k)}{P(\Omega_{a}/\Omega_{d}=0,\mathbf{b}_{\mathrm{cosmo}};k)}\to
1\;\;\;\mathrm{when}\;\;\;\Omega_{a}/\Omega_{d}\to
0\;\;\;\forall\;\mathbf{b}_{\mathrm{cosmo}},\;k,$ (C.1)

where $\mathbf{b}_{\mathrm{cosmo}}$ denotes the vector of cosmological
parameters without including $A_{s}$777Although the power spectrum amplitude
is varied in this study, it does not affect the transfer function which is a
ratio of power spectra.. So we expect the axion transfer function to be
independent of cosmological parameters at sufficiently small concentrations.
We note that the prior given for the axion fraction includes a wide range of
values for which this does not necessarily hold. This observation is
illustrated by the difference in transfer function variation between Fig.
19(a) and Fig. 19(b).

(a)

(b)

Figure 19: Variation of the axion transfer function when fixing axion mass and
density while perturbing cosmological parameters about their fiducial values.
Top: Variation with an axion density of 5% with $m_{a}=10^{-27}$ eV. Bottom:
Variation with an axion density of 30% at the same axion mass.

The interpolator can be corrected by including prefactors accounting for the
correlation between the axion density and the cosmological parameters.

$\displaystyle
T_{\mathrm{ax,\;corrected}}(\mathbf{b}_{\mathrm{axion}},\mathbf{b}_{\mathrm{cosmo}};k)=\prod_{b_{i}\in\mathbf{b}_{\mathrm{cosmo}}}\mathcal{C}_{b_{i}}(b;k)\times
T_{\mathrm{ax,\;naive}}(\mathbf{b}_{\mathrm{axion}};k),$ (C.2)

where

$\displaystyle\mathcal{C}_{b_{i}}(b,k)\equiv\frac{T_{\mathrm{ax}}(\mathbf{b}_{\mathrm{axion}},b_{i}=b;k)}{T_{\mathrm{ax,\;fid}}(\mathbf{b}_{\mathrm{axion}};k)}.$
(C.3)

Figure 20: Relative error from the two interpolation method described in this
section with respect to the axionCAMB result. The initial interpolation is
done assuming the transfer function depends only on axion parameters and is
independent of cosmological ones.

In other words, we create four distinct interpolators. The first three vary a
single cosmological parameter along with the axion fraction and evaluates the
axion transfer function as a function of scale up to $k=0.5\;h/$Mpc. We then
divide the result by the axion transfer function with the same axion fraction
and at the fiducial values of all the cosmological parameters. This
calculation captures the effect of varying both a cosmological parameter and
the axion fraction simultaneously. The normalization ensures that
$\mathcal{C}_{b_{i}}(b,k)\to 1$ when $b\to b_{\mathrm{fid}}$. After
calculating the correction factors for $\\{\omega_{b},\omega_{d},H_{0}\\}$, we
calculate the naive axion transfer function which accounts for the impact of
axions while keeping the cosmological parameters constant. The product of the
four captures the effects of axions and the correlation between axions and
cosmological parameters. Note that this process is repeated for each axion
mass bin as the axion mass is fixed for each run. Accounting for these small
correction improves the accuracy of the interpolation scheme as shown in Fig.
20 where the transfer function was calculated for three randomly sampled
choices of cosmological parameters and compared to the axionCAMB result. This
improved scheme allows the interpolation to maintain percent level error with
the axionCAMB output at an axion concentration of
$\Omega_{a}/\Omega_{d}\lesssim 0.1$ which, as shown in Fig. 1, is where more
than 95% of the samples are taken.

The final step in using the interpolator is to multiply the resulting transfer
function with the Class matter power spectrum which interfaces well with the
PyBird functionalities. By not requiring the high-precision parameter in the
input file, we get an average runtime per execution of 2.7 seconds. The
combination of the Class and the interpolation scheme is .7 s per execution
which represents an improvement by a factor of 3.85 per iteration. We run the
axionCAMB for the interpolator for 40 different choices of each cosmological
parameters and for 80 axion fractions. Combined with the naive transfer
function which requires 80 executions, we arrive at a total of $\sim 10^{4}$
executions of the Boltzmann code for the creation of the interpolation tables.
Assuming an average $\sim 2\times 10^{5}$ samples for the chains to converge,
we get that the total runtime on a single core for the Boltzmann code approach
is of $5.4\times 10^{5}$ s per chain while the interpolation scheme is at
$1.7\times 10^{5}$ s which represents a reduction of about 70% in total
computational cost.

## References

  * [1] L. F. Abbott and P. Sikivie.  A cosmological bound on the invisible axion.  Physics Letters B, 120(1-3):133–136, January 1983.
  * [2] J. Preskill, M. B. Wise, and F. Wilczek.  Cosmology of the invisible axion.  Physics Letters B, 120(1-3):127–132, January 1983.
  * [3] M. Dine and W. Fischler.  The Not So Harmless Axion.  Physics Letters B, 120:137–141, 1983.
  * [4] G. Ballesteros, J. Redondo, A. Ringwald, and C. Tamarit.  Standard Model—axion—seesaw—Higgs portal inflation. Five problems of particle physics and cosmology solved in one stroke.  JCAP, 2017(8):001, August 2017.
  * [5] F. Wilczek.  A Model of Anthropic Reasoning, Addressing the Dark to Ordinary Matter Coincidence.  arXiv e-prints, pages hep–ph/0408167, August 2004.
  * [6] A. Arvanitaki et al.  String axiverse.  Physical Review D, 81(12):123530, 2010.
  * [7] M. Tanabashi et al.  Review of particle physics.  Phys. Rev. D, 98:030001, Aug 2018.
  * [8] T. C. Bachlechner, K. Eckerle, O. Janssen, and M. Kleban.  Multiple-axion framework.  Phys. Rev. D, 98(6):061301, 2018.
  * [9] V. M. Mehta, M. Demirtas, C. Long, D. J. E. Marsh, L. McAllister, and M. J. Stott.  Superradiance in String Theory.  arXiv e-prints, page arXiv:2103.06812, March 2021.
  * [10] L. Hui.  Wave Dark Matter.  arXiv e-prints, page arXiv:2101.11735, January 2021.
  * [11] D. Grin et al.  Gravitational probes of ultra-light axions.  arXiv preprint arXiv:1904.09003, 2019.
  * [12] J. Diehl and J. Weller.  Constraining Ultra-light Axions with Galaxy Cluster Number Counts.  arXiv e-prints, page arXiv:2103.08674, March 2021.
  * [13] A. Klypin, A. V. Kravtsov, O. Valenzuela, and F. Prada.  Where Are the Missing Galactic Satellites?  APJ, 522(1):82–92, September 1999.
  * [14] B. Moore et al.  Dark Matter Substructure within Galactic Halos.  APJ Letters, 524(1):L19–L22, October 1999.
  * [15] W. J. G. de Blok.  The Core-Cusp Problem.  Advances in Astronomy, 2010:789293, January 2010.
  * [16] M. Boylan-Kolchin, J. S. Bullock, and M. Kaplinghat.  Too big to fail? The puzzling darkness of massive Milky Way subhaloes.  Monthly Notices of the Royal Astronomical Society, 415(1):L40–L44, July 2011.
  * [17] D. N. Spergel and P. J. Steinhardt.  Observational Evidence for Self-Interacting Cold Dark Matter.  Physical Review Letters, 84(17):3760–3763, April 2000.
  * [18] J. Sommer-Larsen and A. Dolgov.  Formation of Disk Galaxies: Warm Dark Matter and the Angular Momentum Problem.  The Astrophysical Journal, 551(2):608–623, April 2001.
  * [19] A. Pontzen and F. Governato.  How supernova feedback turns dark matter cusps into cores.  Monthly Notices of the Royal Astronomical Society, 421(4):3464–3471, April 2012.
  * [20] A. Pontzen and F. Governato.  Cold dark matter heats up.  Nature, 506(7487):171–178, February 2014.
  * [21] E. Tollet et al.  NIHAO - IV: core creation and destruction in dark matter density profiles across cosmic time.  Monthly Notices of the Royal Astronomical Society, 456(4):3542–3552, March 2016.
  * [22] A. Fitts et al.  fire in the field: simulating the threshold of galaxy formation.  Monthly Notices of the Royal Astronomical Society, 471(3):3547–3562, November 2017.
  * [23] D. J. E. Marsh and J. Silk.  A model for halo formation with axion mixed dark matter.  Monthly Notcies of the Royal Astronomical Society, 437(3):2652–2663, January 2014.
  * [24] D. J. E. Marsh and A.-R. Pop.  Axion dark matter, solitons and the cusp-core problem.  Monthly Notices of the Royal Astronomical Society, 451(3):2479–2492, August 2015.
  * [25] R. H. Victor, J. S. Bullock, and M. Boylan-Kolchin.  Scalar field dark matter: helping or hurting small-scale problems in cosmology?  Monthly Notices of the Royal Astronomical Society, 483(1):289–298, 2019.
  * [26] M. Safarzadeh and D. N. Spergel.  Ultra-light dark matter is incompatible with the milky way’s dwarf satellites.  The Astrophysical Journal, 893(1):21, 2020.
  * [27] W. Hu, R. Barkana, and A. Gruzinov.  Fuzzy Cold Dark Matter: The Wave Properties of Ultralight Particles.  Physics Review Letters, 85(6):1158–1161, August 2000.
  * [28] R. Hložek, D. J. E. Marsh, and D. Grin.  Using the full power of the cosmic microwave background to probe axion dark matter.  Monthly Notices of the Royal Astronomical Society, 476(3):3063–3085, May 2018.
  * [29] E. Boera, G. D. Becker, J. S. Bolton, and F. Nasir.  Revealing Reionization with the Thermal History of the Intergalactic Medium: New Constraints from the Ly$\alpha$ Flux Power Spectrum.  The Astrophysical Journal, 872(1):101, February 2019.
  * [30] S. Chabanier, M. Millea, and N. Palanque-Delabrouille.  Matter power spectrum: from Ly $\alpha$ forest to CMB scales.  Monthly Notices of the Royal Astronomical Society, 489(2):2247–2253, October 2019.
  * [31] V. Iršič, M. Viel, M. G. Haehnelt, J. S. Bolton, and G. D. Becker.  First Constraints on Fuzzy Dark Matter from Lyman-$\alpha$ Forest Data and Hydrodynamical Simulations.  Physical Review Letters, 119(3):031302, July 2017.
  * [32] T. Kobayashi et al.  Lyman-$\alpha$ constraints on ultralight scalar dark matter: Implications for the early and late universe.  Physical Review D, 96(12):123514, 2017.
  * [33] E. Armengaud, N. Palanque-Delabrouille, C. Yèche, D. J. E. Marsh, and J. Baur.  Constraining the mass of light bosonic dark matter using SDSS Lyman-$\alpha$ forest.  Monthly Notices of the Royal Astronomical Society, 471(4):4606–4614, November 2017.
  * [34] K. K. Rogers and H. V. Peiris.  Strong Bound on Canonical Ultralight Axion Dark Matter from the Lyman-Alpha Forest.  Physical Review Letters, 126(7):071302, February 2021.
  * [35] K. K. Rogers and H. V. Peiris.  General framework for cosmological dark matter bounds using N -body simulations.  Physical Review D, 103(4):043526, February 2021.
  * [36] Z. Lukić, C. W. Stark, P. Nugent, M. White, A. A. Meiksin, and A. Almgren.  The Lyman $\alpha$ forest in optically thin hydrodynamical simulations.  Monthly Notices of the Royal Astronomical Society, 446(4):3697–3724, February 2015.
  * [37] J. Oñorbe, J. F. Hennawi, and Z. Lukić.  Self-consistent Modeling of Reionization in Cosmological Hydrodynamical Simulations.  The Astrophysical Journal, 837(2):106, March 2017.
  * [38] DES Collaboration, E. O. Nadler, et al.  Constraints on Dark Matter Properties from Observations of Milky Way Satellite Galaxies.  Physical Review Letters, 126(9):091101, March 2021.
  * [39] A. Maleki, S. Baghram, and S. Rahvar.  Investigation of two colliding solitonic cores in fuzzy dark matter models.  Physical Review D, 101(2):023508, January 2020.
  * [40] R. Hlozek, D. Grin, D. J. E. Marsh, and P. G. Ferreira.  A search for ultralight axions using precision cosmological data.  Physics Review D, 91(10):103512, May 2015.
  * [41] K.-H. Leong, H.-Y. Schive, U.-H. Zhang, and T. Chiueh.  Testing extreme-axion wave-like dark matter using the BOSS Lyman-alpha forest data.  Monthly Notices of the Royal Astronomical Society, 484(3):4273–4286, April 2019.
  * [42] J. B. Bauer et al.  Intensity mapping as a probe of axion dark matter.  Monthly Notices of the Royal Astronomical Society, 500(3):3162–3177, 2021.
  * [43] M. I. Khlopov, B. A. Malomed, and I. B. Zeldovich.  Gravitational instability of scalar fields and formation of primordial black holes.  Monthly Notices of the Royal Astronomical Society, 215:575–589, August 1985.
  * [44] L. Amendola and R. Barbieri.  Dark matter from an ultra-light pseudo-Goldsone-boson.  Physics Letters B, 642(3):192–196, November 2006.
  * [45] J. Cookmeyer, D. Grin, and T. L. Smith.  How sound are our ultralight axion approximations?  Physical Review D, 101(2):023501, January 2020.
  * [46] F. Beutler et al.  The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: anisotropic galaxy clustering in Fourier space.  Monthly Notices of the Royal Astronomical Society, 466(2):2242–2260, April 2017.
  * [47] K. Yamamoto, M. Nakamichi, A. Kamino, B. A. Bassett, and H. Nishioka.  A Measurement of the Quadrupole Power Spectrum in the Clustering of the 2dF QSO Survey.  Publications of the Astronomical Society of Japan, 58:93–102, February 2006.
  * [48] D. Blas, J. Lesgourgues, and T. Tram.  The Cosmic Linear Anisotropy Solving System (CLASS). Part II: Approximation schemes.  JCAP, 2011(7):034, July 2011.
  * [49] J. J. M. Carrasco, M. P. Hertzberg, and L. Senatore.  The effective field theory of cosmological large scale structures.  Journal of High Energy Physics, 2012(9):82, 2012.
  * [50] P. H. Chavanis.  Growth of perturbations in an expanding universe with Bose-Einstein condensate dark matter.  Astronomy and Astrophysics, 537, 2012.
  * [51] A. Suárez and P.-H. Chavanis.  Hydrodynamic representation of the Klein-Gordon-Einstein equations in the weak field limit: I. General formalism and perturbations analysis.  Phys. Rev. D, 92(2), 3 2015.
  * [52] D. J. E. Marsh.  Axion cosmology.  Physics Reports, 643:1–79, July 2016.
  * [53] L. Senatore and M. Zaldarriaga.  The Effective Field Theory of Large-Scale Structure in the presence of Massive Neutrinos.  arXiv e-prints, page arXiv:1707.04698, July 2017.
  * [54] A. Perko, L. Senatore, E. Jennings, and R. H. Wechsler.  Biased Tracers in Redshift Space in the EFT of Large-Scale Structure.  arXiv e-prints, page arXiv:1610.09321, October 2016.
  * [55] M. M. Ivanov and S. Sibiryakov.  Infrared resummation for biased tracers in redshift space.  JCAP, 2018(7):053, July 2018.
  * [56] J. Zheng et al.  The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: a tomographic measurement of structure growth and expansion rate from anisotropic galaxy clustering in Fourier space.  Monthly Notices of the Royal Astronomical Society, 484(1):442–450, March 2019.
  * [57] G. d’Amico et al.  The cosmological analysis of the SDSS/BOSS data from the Effective Field Theory of Large-Scale Structure.  JCAP, 2020(5):005, May 2020.
  * [58] H. Gil-Marín et al.  The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: RSD measurement from the LOS-dependent power spectrum of DR12 BOSS galaxies.  Monthly Notices of the Royal Astronomical Society, 460(4):4188–4209, August 2016.
  * [59] M. M. Ivanov, M. Simonović, and M. Zaldarriaga.  Cosmological parameters from the BOSS galaxy power spectrum.  Journal of Cosmology and Astroparticle Physics, 2020(5):042, May 2020.
  * [60] H. Gil-Marín et al.  The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: BAO measurement from the LOS-dependent power spectrum of DR12 BOSS galaxies.  Monthly Notices of the Royal Astronomical Society, 460(4):4210–4219, August 2016.
  * [61] F.-S. Kitaura et al.  The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: mock galaxy catalogues for the BOSS Final Data Release.  Monthly Notices of the Royal Astronomical Society, 456(4):4156–4173, 01 2016.
  * [62] S. A. Rodríguez-Torres et al.  The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: modelling the clustering and halo occupation distribution of BOSS CMASS galaxies in the Final Data Release.  Monthly Notices of the Royal Astronomical Society, 460(2):1173–1187, August 2016.
  * [63] F. S. Kitaura, G. Yepes, and F. Prada.  Modelling baryon acoustic oscillations with perturbation theory and stochastic halo biasing.  Monthly Notices of the Royal Astronomical Society, 439:L21–L25, March 2014.
  * [64] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman.  emcee: The mcmc hammer.  PASP, 125:306–312, 2013.
  * [65] A. Chudaykin, M. M. Ivanov, O. H. E. Philcox, and M. Simonović.  Nonlinear perturbation theory extension of the Boltzmann code CLASS.  Physical Review D, 102(6):063533, September 2020.
  * [66] T. Nishimichi, G. D’Amico, M. M. Ivanov, L. Senatore, M. Simonović, M. Takada, M. Zaldarriaga, and P. Zhang.  Blinded challenge for precision cosmology with large-scale structure: Results from effective field theory for the redshift-space galaxy power spectrum.  Physical Review D, 102(12):123541, December 2020.
  * [67] J. Dunkley et al.  Fast and reliable Markov chain Monte Carlo technique for cosmological parameter estimation.  Monthly Notices of the Royal Astronomical Society, 356(3):925–936, January 2005.
  * [68] G. Stein, M. A. Alvarez, and J. R. Bond.  The mass-Peak Patch algorithm for fast generation of deep all-sky dark matter halo catalogues and its N-body validation.  Monthly Notices of the Royal Astronomical Society, 483(2):2236–2250, February 2019.
  * [69] G. D’Amico and P. Zhang L. Senatore.  Limits on wCDM from the EFTofLSS with the PyBird code.  JCAP, 2021(1):006, January 2021.
  * [70] E. Castorina, E. Sefusatti, R. Sheth, F. Villaescusa-Navarro, and M. Viel.  Cosmology with massive neutrinos ii: on the universality of the halo mass function and bias.  Journal of Cosmology and Astroparticle Physics, 2014:049–049, 2013\.
  * [71] E. Castorina, C. Carbone, J. Bel, E. Sefusatti, and K. Dolag.  DEMNUni: the clustering of large-scale structures in the presence of massive neutrinos.  Journal of Cosmology and Astroparticle Physics, 2015(7):043, July 2015.
  * [72] J. L. Tinker et al.  Cosmological constraints from galaxy clustering and the mass-to-number ratio of galaxy clusters.  The Astrophysical Journal, 745(1):16, 2011.
  * [73] F. Beutler et al.  The clustering of galaxies in the sdss-iii baryon oscillation spectroscopic survey: testing gravity with redshift space distortions using the power spectrum multipoles.  Monthly Notices of the Royal Astronomical Society, 443(2):1065–1089, 2014.
  * [74] G. Schwarz et al.  Estimating the dimension of a model.  Annals of statistics, 6(2):461–464, 1978.
  * [75] Planck Collaboration, N. Aghanim, et al.  Planck 2018 results. VI. Cosmological parameters.  AAP, 641:A6, September 2020.
  * [76] W. Hu, D. J. Eisenstein, and M. Tegmark.  Weighing neutrinos with galaxy surveys.  Physical Review Letters, 80(24):5255, 1998.
  * [77] D. J. E. Marsh and P. G. Ferreira.  Ultralight scalar fields and the growth of structure in the Universe.  Physical Review D, 82(10):103528, November 2010.
  * [78] D. J. E. Marsh, E. Macaulay, M. Trebitsch, and P. G. Ferreira.  Ultralight axions: Degeneracies with massive neutrinos and forecasts for future cosmological observations.  Physical Review D, 85(10):103514, May 2012.
  * [79] P. Christian et al.  Massive neutrinos and degeneracies in lyman-alpha forest simulations.  Journal of Cosmology and Astroparticle Physics, 2020(04):025, 2020\.
  * [80] D. Coe.  Fisher Matrices and Confidence Ellipses: A Quick-Start Guide and Software.  arXiv e-prints, page arXiv:0906.4123, June 2009.
  * [81] M. J. Reid, D. W. Pesce, and A. G. Riess.  An Improved Distance to NGC 4258 and Its Implications for the Hubble Constant.  The Astrophysical Journal, 886(2):L27, December 2019.
  * [82] J. L. Bernal, L. Verde, and A. G. Riess.  The trouble with H0.  Journal of Cosmology and Astroparticle Physics, 2016(10):019, October 2016.
  * [83] V. Poulin, T. L. Smith, D. Grin, T. Karwal, and M. Kamionkowski.  Cosmological implications of ultralight axionlike fields.  Physical Review D, 98(8):083525, October 2018.
  * [84] C. Blake et al.  The WiggleZ Dark Energy Survey: the selection function and z = 0.6 galaxy power spectrum.  Monthly Notices of the Royal Astronomical Society, 406(2):803–821, August 2010.
  * [85] Y. Minami and E. Komatsu.  New Extraction of the Cosmic Birefringence from the Planck 2018 Polarization Data.  Physical Review Letters, 125(22):221301, November 2020.
  * [86] P. Agrawal, A. Hook, and J. Huang.  A CMB Millikan experiment with cosmic axiverse strings.  Journal of High Energy Physics, 2020(7):138, July 2020.
  * [87] M. Jain, A. J. Long, and M. A. Amin.  CMB birefringence from ultra-light axion string networks.  arXiv e-prints, page arXiv:2103.10962, March 2021.
  * [88] A. Laguë, J. R. Bond, R. Hložek, D. J. E. Marsh, and L. Söding.  Evolving ultralight scalars into non-linearity with lagrangian perturbation theory.  Monthly Notices of the Royal Astronomical Society, 03 2021.  stab601.
  * [89] B. Schwabe et al.  Simulating mixed fuzzy and cold dark matter.  Physics Review D, 102(8):083518, October 2020.
  * [90] R. Hložek, D. J. E. Marsh, D. Grin, R. Allison, J. Dunkley, and E. Calabrese.  Future CMB tests of dark matter: Ultralight axions and massive neutrinos.  Physical Review D, 95(12):123511, June 2017.
  * [91] L. Chris et al.  SciNet: Lessons learned from building a power-efficient top-20 system and data centre.  Journal of Physics: Conference Series, 256:012026, nov 2010.
  * [92] M. Ponce et al.  Deploying a Top-100 Supercomputer for Large Parallel Workloads: the Niagara Supercomputer.  arXiv e-prints, page arXiv:1907.13600, July 2019.
  * [93] A. Lewis.  GetDist: a Python package for analysing Monte Carlo samples.  2019\.
  * [94] A. M. Price-Whelan and D. Foreman-Mackey.  schwimmbad: A uniform interface to parallel processing pools in python.  The Journal of Open Source Software, 2(17), sep 2017.
  * [95] F. Bernardeau, S. Colombi, E.Gaztañaga, and R. Scoccimarro.  Large-scale structure of the universe and cosmological perturbation theory.  Physics Reports, 367(1-3):1–248, Sep 2002.
  * [96] L. Senatore and M. Zaldarriaga.  Redshift Space Distortions in the Effective Field Theory of Large Scale Structures.  arXiv e-prints, page arXiv:1409.1225, September 2014.
  * [97] T. Baldauf, E. Schaan, and M. Zaldarriaga.  On the reach of perturbative methods for dark matter density fields.  Journal of Cosmology and Astroparticle Physics, 2016(03):007, 2016\.
  * [98] X. Li et al.  Numerical and perturbative computations of the fuzzy dark matter model.  Physical Review D, 99(6), 3 2019.
  * [99] N. Kaiser.  Clustering in real space and in redshift space.  Monthly Notices of the Royal Astronomical Society, 227(1):1–21, 1987.
  * [100] A. J. S. Hamilton.  Linear redshift distortions: A review.  The evolving universe, pages 185–275, 1998.