# Fuzzy Dark Matter and the Dark Energy Survey Year 1 Data

Mona Dentler1, David J. E. Marsh2, Renée Hložek3,4, Alex Laguë3,4,5, Keir K.
Rogers3, Daniel Grin6.
1Institut für Astrophysik, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany
2King’s College London, Strand, London, WC2R 2LS, United Kingdom
3Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50
St. George Street, Toronto, Ontario, M5S 3H4, Canada
4Department of Astronomy and Astrophysics, University of Toronto, 50 St.
George Street, Toronto, Ontario, M5S 3H4, Canada
5Canadian Institute for Theoretical Astrophysics, University of Toronto, 60
St. George St., Toronto, ON, M5S 3H8, Canada
6Department of Physics and Astronomy, Haverford College, 370 Lancaster Avenue,
Haverford, PA 19041, United States Corresponding Author:
<EMAIL_ADDRESS>
###### Abstract

Gravitational weak lensing by dark matter halos leads to a measurable imprint
in the shear correlation function of galaxies. Fuzzy dark matter (FDM),
composed of ultralight axion-like particles of mass $m\sim 10^{-22}\text{
eV}$, suppresses the matter power spectrum and shear correlation with respect
to standard cold dark matter. We model the effect of FDM on cosmic shear using
the optimised halo model HMCode, accounting for additional suppression of the
mass function and halo concentration in FDM as observed in $N$-body
simulations. We combine Dark Energy Survey year 1 (DES-Y1) data with the
_Planck_ cosmic microwave background anisotropies to search for shear
correlation suppression caused by FDM. We find no evidence of suppression
compared to the preferred cold DM model, and thus set a new lower limit to the
FDM particle mass. Using a log-flat prior and marginalising over uncertainties
related to the non-linear model of FDM, we find a new, independent 95% C.L.
lower limit $\log_{10}m>-23$ combining _Planck_ and DES-Y1 shear, an
improvement of almost two orders of magnitude on the mass bound relative to
CMB-only constraints. Our analysis is largely independent of baryonic
modelling, and of previous limits to FDM covering this mass range. Our
analysis highlights the most important aspects of the FDM non-linear model for
future investigation. The limit to FDM from weak lensing could be improved by
up to three orders of magnitude with $\mathcal{O}(0.1)$ arcmin cosmic shear
angular resolution, if FDM and baryonic feedback can be simultaneously
modelled to high precision in the halo model.

###### keywords:

cosmology: theory, dark matter, elementary particles

††pagerange: Fuzzy Dark Matter and the Dark Energy Survey Year 1
Data–LABEL:lastpage

##  1 Introduction

Dark matter (DM) is one of the most pressing issues in modern particle physics
and cosmology, with its existence confirmed by a range of observations
covering scales from the galactic neighbourhood to the cosmos (Aghanim et al.,
2020b; de Salas & Widmark, 2021). At one extreme end of DM candidates reside
primordial black holes (PBHs), with macroscopic masses measured in solar
masses, $M_{\odot}$. Observations constrain the parameter space of PBHs
tightly, leaving a single window around $M\approx 10^{-12}M_{\odot}$ where
they can compose all of the DM with a monochromatic mass function, and with
tight constraints on the PBH fraction across the rest of the parameter space
(Green & Kavanagh, 2021). At the other extreme end of DM parameter space are
ultralight bosons, including pseudoscalars such as the axion (Peccei & Quinn,
1977; Weinberg, 1978; Wilczek, 1978; Abbott & Sikivie, 1983; Preskill et al.,
1983; Dine & Fischler, 1983), scalars (e.g. Turner, 1983; Li et al., 2014),
vectors (e.g. Graham et al., 2016), and tensors (e.g Babichev et al., 2016).

In this work we are concerned with ultralight scalars and pseudoscalars that
can be treated as having vanishing non-gravitational self interactions at
early times (Marsh & Hoof, 2021). We consider DM composed of a scalar field,
$\phi$, with potential $V(\phi)=m^{2}\phi^{2}/2$, where $m$ is the particle
mass, and neglect all other DM interactions (hence the model applies to
scalars and pseudoscalars equally). We consider masses in the range of
$10^{-27}\text{ eV}\lesssim m\lesssim 10^{-19}\text{eV}$. We refer to such
particles as Fuzzy Dark Matter (FDM), or “ultralight bosonic dark matter”
(UBDM). Ultralight axions (ULAs) provide one particle physics model for FDM,
valid in the mass range of interest for “decay constants” in the range
$f_{a}\approx 10^{17}\text{ GeV}$. Such ULAs are a feature of the string
theory landscape (Arvanitaki et al., 2010; Marsh, 2016b; Hui et al., 2017;
Mehta et al., 2021; Cicoli et al., 2022), and can also be realised in field
theory models (Kim & Marsh, 2016; Davoudiasl & Murphy, 2017).

In the context of cosmology, the key phenomenological feature of FDM/ULAs is
the existence of a scale depending on the particle mass $m$, called the Jeans
scale $k_{J}(m)$ (Khlopov et al., 1985). At smaller wavenumbers $k<k_{J}$, FDM
is phenomenologically virtually indistinguishable from cold DM (CDM). However,
at higher wavenumbers $k>k_{J}$, the scalar field gradient energy leads to an
effect called “quantum pressure”, 111Erroneously, since it is neither quantum,
nor a pressure. which is absent in the case of CDM. The presence of quantum
pressure below the Jeans scale counteracts gravity and results in a
suppression of the matter density power spectrum, $P(k)$ (Hu et al., 2000;
Marsh & Ferreira, 2010; Marsh, 2016b; Hui et al., 2017; Marsh & Hoof, 2021).

Our evidence for DM comes entirely from its gravitational interactions. So too
come the most widely applicable and generic constraints on the properties of
DM, in particular the lower bound to the particle mass, $m$. Cosmic microwave
background (CMB) temperature, polarisation, and gravitational lensing
anisotropies measured by the _Planck_ satellite (Aghanim et al., 2020a)
establish the lower limit $m\gtrsim\mathcal{O}(\text{few})\times
10^{-25}\text{ eV}$ for the dominant component of DM (Hlozek et al., 2015,
2018; Poulin et al., 2018; Linares Cedeño et al., 2021). The CMB bound is
extremely robust, since it relies only on linear physics of cosmological
perturbation theory and decoupling of the photon-baryon plasma, and
marginalises over uncertainties in the cosmological parameters (baryon
density, DM density, primordial power, Hubble parameter, and CMB optical
depth).

Extending the lower limit to the FDM mass requires probes of the power
spectrum on scales smaller than those accessible from the CMB, moving towards
the quasi-linear and fully non-linear regimes of structure formation. Many
independent probes have been considered, including from the epoch of
reionisation (e.g. Bozek et al., 2015; Sarkar et al., 2016; Schive et al.,
2016; Corasaniti et al., 2017). In these studies, cosmological and systematic
uncertainties are left largely unaccounted for. A more rigorous approach to
small and non-linear scales is required.

One such rigorous approach has been developed for the Lyman-$\alpha$ forest
flux power spectrum (see e.g. Narayanan et al., 2000; McDonald et al., 2006;
Viel et al., 2013; Iršič et al., 2017a). The highest-resolution Lyman-$\alpha$
forest probes the matter power spectrum as traced by the intergalactic medium
(IGM) for line-of-sight velocity wavenumbers
$k_{\mathrm{f}}<0.2\,\text{s}\,\text{km}^{-1}$ and redshifts $4.2<z<5$ (Boera
et al., 2019). The most accurate model on such small scales is a cosmological
hydrodynamical simulation (e.g. Lukić et al., 2015), ideally run in pairs in
order to suppress sample variance (Anderson et al., 2019).

The computational cost of such simulations necessitates, for parameter
inference, the use of _emulators_ : interpolation methods in high-dimensional
parameter spaces222This is despite the computational efficiency in
Lyman-$\alpha$ forest simulations where hot and dense gas particles are
converted into collisionless particles. This has negligible impact on
Lyman-$\alpha$ forest statistics, which are sourced in the lower-density IGM
(Viel et al., 2004). (e.g. Lawrence et al., 2010a). In particular, Bird et al.
(2019); Rogers et al. (2019); Rogers & Peiris (2021b, a) built Lyman-$\alpha$
forest emulators capable of marginalising the astrophysical uncertainties of
IGM temperature, energy deposition in the IGM, and the strength of the
ionising background, in addition to cosmological parameters, and the FDM mass,
leading to the rigorous bound $m>2\times 10^{-20}\text{ eV}$ at 95%
credibility. The methodology to marginalise uncertainties improves on previous
Lyman-$\alpha$ forest bounds on ULAs (Iršič et al., 2017b; Armengaud et al.,
2017; Kobayashi et al., 2017).

A shortcoming of the constraints to FDM from the Lyman-$\alpha$ forest is that
the full power of a joint likelihood analysis with the CMB has not yet been
exploited (although see e.g. Seljak et al., 2006, for the case of other
cosmological parameters). Studies using the Lyman-$\alpha$ forest neglect the
interplay of the quantum pressure and baryons in ways which may weaken bounds
(Zhang et al., 2018). However, large-scale structure simulations including FDM
dynamics indicate that, for current data, the effect is only on the order of a
few percent (Nori et al., 2019). Further, existing Lyman-$\alpha$ forest
bounds neglect the large-scale fluctuations ($\sim 40\,\mathrm{Mpc}$) in IGM
temperature and ionisation arising from a spatially-inhomogeneous reionisation
(Hui et al., 2017). Although, initial studies (Wu et al., 2019; Molaro et al.,
2021) suggest that the effect on the small scales (sub-Mpc) that drive dark
matter bounds are currently statistically negligible. Nonetheless, this
motivates complementary and rigorous probes of the dark matter nature, in
particular exploiting scales intermediate to those traced by the CMB primary
anisotropies ($\sim\mathrm{Gpc}$) and the high-resolution Lyman-$\alpha$
forest ($\sim\mathrm{Mpc}$).

Figure 1: One dimensional posterior on the FDM particle mass, $m$, assuming a
log-flat prior. The DES cosmic shear correlation function, $\xi_{\pm}$ probes
scales one order of magnitude smaller than those probed by the _Planck_ CMB
anisotropies. Thus, given the absence of evidence for scale dependence of
$\xi_{\pm}$ induced by the FDM Jeans scale, $k_{J,eq}\propto m^{1/2}$, Eq.
(1), we are able to place a new lower limit to $m$ that is almost two orders
of magnitude tighter than the limit from the CMB alone. This new limit is
unaffected by a variety of uncertainties associated to the non-linear model.
The small peak in the posterior using DES alone arises due to degeracies that
are broken by the combination with CMB data, and is not statistically
significant.

An alternative and complementary observable that probes small and non-linear
scales is provided by the galaxy shear weak lensing power spectrum, and is the
subject of the present work. In particular, we make use of the Dark Energy
Survey Year 1 (DES-Y1, Abbott et al., 2018a) weak lensing likelihood
implemented in the publicly available cosmosis (Zuntz et al., 2015) analysis
framework. 333While this work was in preparation the DES-Y3 results were
released as preprints (Pandey et al., 2021; Porredon et al., 2021). We will
study these in a future publication, once the DES-Y3 likelihoods have been
publicly released. A virtue of weak lensing is that the theoretical input is
the total matter power spectrum, $P(k)$, offering a direct probe of the
clustering not just of a tracer such as hydrogen gas in the IGM, but of the DM
itself.

The standard approach to the non-linear scales in $P(k)$ adopted in weak
lensing analyses by DES and the Kilo Degree Survey (KiDS, Joudaki et al.,
2020) is the _Halo Model_ (reviewed in e.g. Cooray & Sheth, 2002). The primary
DES-Y1 analysis uses the halo model inspired fitting function halofit (Smith
et al., 2003; Bird et al., 2012; Takahashi et al., 2012). halofit is, however,
not calibrated for use with ULA/FDM cosmologies (Hložek et al., 2017), and so
a more physical approach is required. The halo model implementation HMCode
(Mead et al., 2015) provides such a model, which can be adapted to both warm
DM (WDM) and FDM (Marsh, 2016a). 444HMCode modifications in Marsh (2016a)
accurately reproduce the WDM non-linear $P(k)$ ratio to CDM in Corasaniti et
al. (2017) to $\mathcal{O}(\rm few)\%$ accuracy.  HMCode is a complementary
and alternative method to power spectrum emulators. It is useful for weak
lensing, which is dominated by high density regions, where an analytical model
saves considerable computational cost. Using the halo model, weak lensing bias
parameters and other uncertainties in the model are easily accounted for in
Bayesian parameter estimation, i.e. marginalised over.

The halo model in HMCode is modified with physically motivated fitting
parameters in order to match emulators and simulations of the non-linear power
spectrum to within around 2% accuracy. HMCode also includes a model for the
important effect of baryonic feedback from Active Galactic Nuclei (AGNs),
which affects $P(k)$ at the low redshifts $0\lesssim z\lesssim 2$ of the
DES-Y1 galaxy sample. In the present analysis, following Abbott et al.
(2018a), we mask small scales in the DES-Y1 data that are affected by AGN
feedback. We comment in our conclusions on the effect of including these
scales in future analyses.

As described below, the halo model uses physical inputs of the linear matter
power spectrum, the halo mass function (HMF), $n(M)$, and halo density profile
parameterised by the halo concentration-mass relation, $c(M)$. This makes the
halo model adaptable to different DM theories, and can account for physics
observed in single simulations, but not performed in significant numbers to
build an emulator around. HMCode further introduces a phenomenological model
for the transition between the one halo and two halo terms in the quasi-linear
regime, calibrated to $N$-body simulations with different initial transfer
functions that cover a range of scales of interest for FDM.

In the following we introduce physically motivated modifications to HMCode
based on observations from $N$-body simulations to model the non-linear scales
of FDM. We then propagate (approximated) systematic uncertainties in the
models adopted for $n(M)$, $c(M)$, and the quasi-linear smoothing, along with
all relevant DES nuisance parameters. This allows us to perform a search for
evidence of the FDM Jeans scale in the galaxy shear correlation function.

We furthermore perform a _global fit_ incorporating both DES and _Planck_
data, which are highly complementary in the cosmological parameter space.
_Planck_ effectively anchors all the primary cosmological parameters, leaving
the DES correlation functions as a lever arm into small scales, increasing the
sensitivity to the FDM Jeans scale (and thus $m$) without any degeneracy with
other parameters of the model.

We can anticipate our results by considering how weak lensing data can be used
as a probe of the linear power spectrum, $P(k)$ (Tegmark & Zaldarriaga, 2002).
As shown in Chabanier et al. (2019) the _Planck_ CMB constrains $P(k)$ for
wavevectors $k<0.3\,h\,\text{Mpc}^{-1}$, while the final DES bin covers
$1\,h\,\text{Mpc}^{-1}<k<5\,h\,\text{Mpc}^{-1}$. The ULA power spectrum is
suppressed below the Jeans scale at matter-radiation equality, given by (Hu et
al., 2000):

$k_{\rm J,eq}=12.9\,h\text{
Mpc}^{-1}\left(\frac{0.7}{h}\right)\left(\frac{m}{10^{-22}\text{
eV}}\right)^{1/2}$ (1)

Thus, extending constraints to $P(k)$ by an order of magnitude in $k$ using
DES compared to _Planck_ alone should leverage two orders of magnitude in
sensitivity to the FDM mass. This improvement in the limit to $m$ driven by
DES is evident in the one dimensional posterior derived from our analysis,
shown in Fig. 1 (these results are discussed further in Section 5).

This paper is organised as follows. We start by introducing cosmic shear as a
probe of dark matter and outline the halo model formalism including the
relevant changes to the model in the case of FDM in Section 2. We outline our
statistical methodology and the DES-Y1 and Planck data in Section 3. We
present our results in Section 4 and conclude in Section 5. The Supplementary
Material considers the introduction of the galaxy correlation function, which
we omit from our main analysis due to uncertainty in the use of a scale
independent galaxy bias, discusses aspects of the ULA halo model not included
in our analysis, and discusses massive neutrinos in more detail.

##  2 Observables and Models

###  2.1 DES-Y1 Observables

We now introduce the observables analyzed in this study. We show how these
observables depend on the total matter density power spectrum $P(k)$, which
ultimately allows us to set constraints on the viable range of FDM particle
mass $m$.

The DES-Y1 survey Abbott et al. (2018a) measures the distribution of the
number and ellipticites of galaxies, covering a total redshift range of
$0\lesssim z\lesssim 2$. The data are divided in tomographic redshift bins
$i$, where the number of bins is four or five, depending on the respective
observable. Each redshift bin $i$ contains a distribution of galaxies
$n^{i}(z)$, and the total number of galaxies in each bin is given
$\bar{n}^{i}=\int{\rm d}z\,n^{i}(z)$. The measured distributions of galaxies
and their ellipticities yield a set of three two-point observables, which we
refer to as ‘3x2pt’ below. It consists of the individual observables galaxy
clustering $w(\theta)$, galaxy-galaxy lensing $\gamma_{t}(\theta)$, and the
two components of the cosmic shear $\xi^{\pm}(\theta)$. In detail, $w(\theta)$
measures the distribution of pairs of angular separations of galaxies as
compared to what is expected assuming a random distribution. The observables
$\gamma_{t}(\theta)$ and $\xi_{\pm}(\theta)$ quantify shape distortions due to
lensing by foreground mass distributions. $\gamma_{t}(\theta)$ measures, for
pairs of source and lens galaxies, the ellipticities of the the source
galaxies tangential to the line connecting to the lens galaxy.
$\xi_{+}(\theta)$ and $\xi_{-}(\theta)$, on the other hand, express for pairs
of source galaxies the sum (+) and the difference (-) of the product of the
tangential ellipticities and the cross ellipticities.

The effects measured by DES-Y1 are all sourced by the distribution of matter
along the light of sight, and hence the corresponding 2pt observables are
predicted via projections of the power spectrum $P(\ell,\chi)$ along this
direction (Abbott et al., 2018a; Krause et al., 2017):

$\displaystyle P^{ij}_{\kappa\kappa,\kappa
g,gg}(\ell)=\int\rm{d}\chi\,\frac{q^{i}_{g/\kappa}q^{j}_{g/\kappa}}{\chi^{2}}\times
P\left(\frac{\ell+1/2}{\chi},z(\chi)\right),$ (2)

where the weighting functions $q^{i}_{g/\kappa}$ are defined as

$\displaystyle
q^{i}_{g}(k,\chi)=b^{i}(k,z(\chi))\frac{n^{i}_{g}(z(\chi))}{\bar{n}^{i}_{g}}\frac{\rm
dz}{\rm d\chi}$ (3)

and

$\displaystyle
q^{i}_{\kappa}(\chi)=\frac{3H_{0}^{2}\Omega_{m}}{2c^{2}}\frac{\chi}{a(\chi)}\int_{\chi}^{\chi_{h}}\,\frac{n^{i}_{\kappa}(z(\chi^{\prime}))\rm
dz/\rm
d\chi^{\prime}}{\bar{n}^{i}_{\kappa}}\frac{\chi^{\prime}-\chi}{\chi^{\prime}}\rm
d\chi^{\prime}.$ (4)

The comoving distance is $\chi$, and $b$ is the linear bias. The final
observables are calculated by transforming to real-space as follows:

$\displaystyle w^{i}(\theta)$ $\displaystyle=\int\frac{\rm
d\ell\,\ell}{2\pi}J_{0}(\theta\ell)P^{ij}_{gg}(\ell)$ (5)
$\displaystyle\gamma^{ij}_{t}(\theta)$ $\displaystyle=(1+m^{j})\int\frac{\rm
d\ell\,\ell}{2\pi}J_{2}(\theta\ell)P^{ij}_{g\kappa}(\ell)$ (6)
$\displaystyle\xi^{ij}_{+/-}(\theta)$
$\displaystyle=(1+m^{i})(1+m^{j})\int\frac{\rm
d\ell\,\ell}{2\pi}J_{0/4}(\theta\ell)P^{ij}_{\kappa\kappa}(\ell),$ (7)

where $J_{\nu}$ denotes the Bessel function of the first kind, and the
parameter $m^{i}$ is the multiplicative shear bias.

When galaxies are used as tracers for the matter, we need the linear bias
$b^{i}(k,z(\chi))$ to relate the distribution of galaxies to the underlying
total matter distribution, resulting in the radial weight function Eq. (3). By
contrast, the lensing efficency Eq. (4) only depends on the matter density and
the galaxy distribution along the line of sight. Consequently, the projection
$P^{ij}_{\kappa g}(\ell)$ and hence the observable $\gamma_{t}$ depends
linearly on the galaxy bias, while the projection $P^{ij}_{gg}(\ell)$ and
hence the observable $w$ depends quadratically on the galaxy bias. Thus,
measuring both $w$ and $\gamma_{t}$ allows the linear bias parameters to be
measured (Abbott et al., 2018a).

Figure 2: Cosmic shear data and model predictions for CDM and FDM with four
different values for the mass parameter $m$. We show our model predictions for
the linear (muted colours) and non-linear (saturated colours) power spectrum
at redshift $z=1$ in the top left panel, the respective predictions for CDM
are shown in dashed, black. The redshift is chosen near the peak of the galaxy
distribution (inset top right pannel) in the fourth DES-Y1 redshift bin. We
show the line-of-sight projected power spectrum for the $4,4$ bin in the top
right panel, with the galaxy distribution used in the integral kernel for the
projection inset. Finally we compare our predictions for the shear observables
$\xi_{+}$ (left) and $\xi_{-}$ (right) to the data in the $4-4$ bin in the
bottom row. We mark DES-Y1 data points in black, indicating the square-roots
of the diagonal entries of the covariance matrix as vertical bars. The gray
shaded regions mask angular scales excluded from the analysis due to modeling
uncertainties. The case $m=10^{-22}\text{ eV}$ is indistinguishable from CDM
over the scales shown. Our model predictions compared to the full DES-Y1 data
are shown in the Supplementary Material.

In FDM models, the bias parameter is expected to have additional scale
dependence due to the Jeans scale (Hlozek et al., 2015). Recently Laguë et al.
(2022) simulated Baryon Oscillation Spectroscopic Survey (BOSS) data in
cosmologies with a sub-dominant fraction of ULAs with $m\leq 10^{-24}\text{
eV}$, modelling the galaxy bias in order to derive rigorous constraints. Such
an analysis for DES, with heavy ULAs composing all the DM, would require fully
non-linear simulations and is beyond the scope of the present work. However,
we expect that for the relatively heavy ULAs considered in the present work
($m\gtrsim 10^{-24}\text{ eV}$ consistent with CMB and BOSS lower bounds) the
galaxy bias for ULAs will be linear and (almost) scale-independent, just like
the CDM galaxy bias. Nonetheless, given our relative ignorance of how to model
the galaxy bias for ULAs in our mass range of interest, we consider the shear
only analysis (using only $\xi_{\pm}$), to be the most conservative, and
present it as our main analysis.

We show the cosmic shear data together with the predictions for CDM as well as
for four different FDM masses $m$ for our fiducial cosmology for the $4,4$ bin
in Fig. 2. The full set of bin combinations is illustrated in Fig. 12. The
full 3x2pt analysis is considered in the Supplementary Material, where we show
how the respective observables $\gamma_{t}$ and $w$ are affected by FDM (Fig.
13). As anticipated in the introduction, for masses $\lesssim 0.5\times
10^{-23}\rm{eV}$, the predictions from FDM model differ visibly from the CDM
case.

For the cosmic shear observables $\xi_{\pm}$, the angular scales $\theta$
measured by DES-Y1 and shown in Fig. 12 are related to the radial distances
$k$ in the power spectrum $P(k)$ through the projection Eq. 2 and the
filtering with the Bessel functions $J_{0/4}$ in Eq. (7). As noted above, we
can use the results by Chabanier et al. (2019) to deduce that the DES-Y1 data
are sensitive to values of $k\lesssim 5\,h\mathrm{Mpc}^{-1}$ and hence to the
quasi-linear and mildly non-linear scales of the power spectrum $P(k)$. As we
discuss below, this sensitivity makes the DES-Y1 data set an interesting probe
of DM. On the other hand, it requires modelling of the power spectrum beyond
linear theory.

###  2.2 The FDM Power Spectrum

Because the DES-Y1 observables defined by Eqs. (5), (6) and (7) are related to
the total matter density power spectrum, they are all sensitive to the
distribution of DM. As motivated above, in our analysis, we focus on cosmic
shear, which depends on the total power spectrum of the gravitational
potential. The gravitational potential is related to the total matter density
by the Poisson equation. It is an unbiased tracer of the total matter content
of the Universe, including DM. Since the power spectrum measures the
clustering of DM under gravity it is sensitive to any departures from CDM.

In the case of FDM, these departures come at different scales. At linear
scales, we predict two different effects. Firstly, the background energy
density evolves differently. For the scales of interest in the present study,
$m\gg 10^{-27}\text{ eV}$, this effect does not play a role (Hlozek et al.,
2015). Secondly, gradients in the FDM field play the role of pressure, and
lead to a Jeans scale suppressing structure formation (Khlopov et al., 1985).
Both of these effects are captured by solving the perturbed Klein-Gordon-
Einstein equations, coupled to the rest of the contents of the Universe via
the metric and the Einstein equations. In the present work, the linear power
spectrum is evaluated with camb (Lewis et al., 2000), modified in the case of
ULAs/FDM to axionCAMB (Hlozek et al., 2015, 2014). 555axionCAMB is based on
the 2013 release of camb. We checked via a visual comparison of the posterior
distributions that our results in the case of CDM are unchanged between the
2013 version of camb and the more recent versions available within cosmosis .
Work on an updated version of axionCAMB is ongoing. axionCAMB solves the
Klein-Gordon equation in the Madelung fluid variables, using a WKB
approximation for the sound speed at late times.

The suppression of the FDM linear power with respect to the CDM case
congruously leads to a different distribution of halos, in particular a
reduction in the number and inner density of halos near to and below the Jeans
mass. We implement these effects on the non-linear power spectrum within a
modified version of the Halo Model as described in Section 2.4. Finally, the
specific properties of FDM leads to wave-like effects on very small scales, in
particular the formation of solitonic cores in the centers of halos. Since
these effects are restricted to scales far below the scales probed by the
DES-Y1 data sets, we do not discuss them in the context of this analysis. The
justification to neglect solitons and certain other effects in our analysis
are discussed in detail in the Supplementary Material.

###  2.3 Halo Model and HMCode

The starting point of all cosmological analyses is linear perturbation theory.
Beyond the linear regime, different variants of the Halo Model provide a semi-
analytical construction of the non-linear power spectrum from the
corresponding power spectrum obtained from linear theory. We summarise the
general idea behind this approach below. We furthermore give a short
description of HMCode (Mead et al., 2015), the specific implementation of the
Halo Model which we use as basis for our adaption of the Halo Model. In a
final step, we specify the modelling choices we made for our ULA/FDM version
of HMCode in S ec. 2.4.

####  2.3.1 Halo Model – Basic Concept

Calculations in cosmological models are often more conveniently done in $k$
space and hence we usually define observables in terms of the power spectrum
$P(k)$. However, the Halo Model is based on physical concepts in real space
like the halo density profile and the corresponding halo mass contained within
a certain radius. We therefore find it more favourable to introduce the Halo
Model in terms of the matter density correlation function
$\xi({\mathbf{x}}_{1},{\mathbf{x}}_{2})$, the Fourier transform of the power
spectrum $P(k)$

Assuming that all matter is contained with halos of density
$\rho({\mathbf{x}})$, we can factorise
$\xi({\mathbf{x}}_{1},{\mathbf{x}}_{2})$ as (Cooray & Sheth, 2002; Hayashi &
White, 2008a)

$\displaystyle\xi({\mathbf{x}}_{1},{\mathbf{x}}_{2})$
$\displaystyle=\langle\delta({\mathbf{x}}_{1}),\delta({\mathbf{x}}_{2})\rangle$
$\displaystyle=\langle\frac{\rho_{h}({\mathbf{x}}_{1})}{\bar{\rho}}\frac{\rho_{h^{\prime}}({\mathbf{x}}_{2})}{\bar{\rho}}\rangle+\langle\left(\frac{\rho_{h}({\mathbf{x}}_{1})-\bar{\rho}}{\bar{\rho}}\right)\left(\frac{\rho_{h}({\mathbf{x}}_{2})-\bar{\rho}}{\bar{\rho}}\right)\rangle$
$\displaystyle\equiv\xi^{2h}({\mathbf{x}}_{1},{\mathbf{x}}_{2})+\xi^{1h}({\mathbf{x}}_{1},{\mathbf{x}}_{2}),$
(8)

where we define the two-halo term $\xi^{2h}$ as the correlation function
between densities attributed to two different halos $h$, $h^{\prime}$, while
the one-halo term $\xi^{1h}$ is defined as the correlation between densities
attributed to the same halo. It is commonly assumed that the density
distribution within a halo, $\rho_{h}$, is a function of space with the total
mass $M$ of the halo as a single parameter, i.e.
$\rho_{h}=\rho_{h}({\mathbf{x}}|M)$. The ensemble average in Eq. (8) can be
replaced by the average over space and the average over the halo mass $M$, as
follows (Cooray & Sheth, 2002)

$\displaystyle\xi^{1h}({\mathbf{x}}_{1},{\mathbf{x}}_{2})$
$\displaystyle=\frac{1}{\bar{\rho}^{2}}\int\mathrm{d}M\,n(M)\times$
$\displaystyle\int\mathrm{d}{\mathbf{r}}\,$
$\displaystyle\left(\rho_{h}({\mathbf{x}}_{1}-{\mathbf{r}}|M)-\bar{\rho}\right)\left(\rho_{h}({\mathbf{x}}_{2}-{\mathbf{r}}|M)-\bar{\rho}\right)$
(9) $\displaystyle\xi^{2h}({\mathbf{x}}_{1},{\mathbf{x}}_{2})$
$\displaystyle=\frac{1}{\bar{\rho}^{2}}\int\mathrm{d}M_{1}\,n(M_{1})\int\mathrm{d}M_{2}\,n(M_{2})\,\times$
$\displaystyle\int\\!\int\mathrm{d}{\mathbf{r}}_{1}\mathrm{d}{\mathbf{r}}_{2}\,$
$\displaystyle\rho_{h_{1}}({\mathbf{x}}_{1}-{\mathbf{r}}_{1}|M_{1})\rho_{h_{2}}({\mathbf{x}}_{2}-{\mathbf{r}}_{2}|M_{2})\,\xi_{hh}({\mathbf{r}}_{1},{\mathbf{r}}_{2}|M_{1},M_{2}).$
(10)

In the above equation, $n(M)$ denotes the halo mass function, which describes
the number of halos per unit volume, as a function of their total mass $M$.
The term $\xi_{hh}({\mathbf{r}}_{1},{\mathbf{r}}_{2}|\allowbreak M_{1},M_{2})$
encodes the correlation between two halos of mass $M_{1}$ at point
${\mathbf{r}}_{1}$ and mass $M_{2}$ at point ${\mathbf{r}}_{2}$, respectively.
If the correlation between halos varies slowly on scales of the order of the
halo sizes, the halo density distributions $\rho_{h_{1}}$, $\rho_{h_{2}}$ in
Eq. (10) can be approximated by Dirac-$\delta$-functions times the respective
masses, i.e. $\rho_{h}({\mathbf{x}}-{\mathbf{r}}|M)\sim
M\,\delta^{D}({\mathbf{x}}-{\mathbf{r}})$. The correlation function $\xi_{hh}$
itself can as first approximation be replaced by the correlation function
$\xi_{lin}$ obtained from linear theory, times the linear bias $b(M)$

$\displaystyle\xi^{2h}({\mathbf{x}}_{1},{\mathbf{x}}_{2})\sim
b(M_{1})b(M_{2})\xi^{lin}({\mathbf{x}}_{1},{\mathbf{x}}_{2})\sim\xi^{lin}({\mathbf{x}}_{1},{\mathbf{x}}_{2}).$
(11)

which is justified at large scales $|{\mathbf{x}}_{1}-{\mathbf{x}}_{2}|\gg 1$
because at large scales we can again assume
$\rho_{h}({\mathbf{x}}-{\mathbf{r}}|M)\sim
M\,\delta^{D}({\mathbf{x}}-{\mathbf{r}})$, and
$\int\mathrm{d}M\,n(M)\,b(M)\,M=\bar{\rho}$ by definition of the linear bias.
At smaller scales, in principal we need to account for the linear bias $b(M)$.
However, at these scales, the total correlation function is dominated by the
one-halo term (this approximation as implemented in HMCode is discussed
further below).

Hence, within the framework of the halo model introduced above, two components
need to be predicted: the halo mass function $n(M)$ and the halo density
profile $\rho_{h}({\mathbf{r}}|M)$. A conventional choice for the halo mass
function is obtained from the ellipsoidal collapse model (Bond & Myers, 1996;
Sheth et al., 2001; Sheth & Tormen, 1999; Cooray & Sheth, 2002) within the
Press & Schechter (1978) approach. The result is

$\displaystyle n^{\rm ST}(M)=$
$\displaystyle\frac{\bar{\rho}}{M}f(\nu)\frac{\mathrm{d}\nu}{\mathrm{d}M}\,,$
(12) $\displaystyle f(\nu)=$ $\displaystyle
A\left[1+\frac{1}{(a\,\nu^{2})^{p}}\right]e^{-a\,\nu^{2}/2}$ (13)

where the parameters $a=0.707$, $p=0.3$ and $A=0.2162$ are obtained from
empirical fits to $\Lambda$CDM simulations (Mead et al., 2015). The variable
$\nu$ is defined by the ratio

$\displaystyle\nu$ $\displaystyle\equiv\frac{\delta_{c}}{\sigma}.$ (14)

In the above equation, the variable $\delta_{c}$ denotes the critical density
defined as the minimum density for which a region collapses under its own
gravity. The variable $\sigma^{2}$, defined by

$\displaystyle\sigma^{2}(M)=\langle\delta_{s}(M)^{2}\rangle,$ (15)

denotes the autocorrelation function measuring the correlation between
smoothed overdensities $\delta_{s}({\mathbf{r}},M)$, which are averaged over a
sphere with radius $R$ as follows:

$\displaystyle\delta_{s}({\mathbf{r}},M)=\frac{3}{4\pi
R(M)^{3}}\int_{0}^{{\mathbf{R}}}\mathrm{d}^{3}x\,\delta({\mathbf{x}}-{\mathbf{r}}),$
(16)

where the smoothing radius $R=R(M)$ is chosen such that the corresponding
sphere contains on average the mass $M=4\pi/3\bar{\rho}R^{3}$ (Mo et al.,
2010).

The conventional model for the second component of the halo model, the halo
density profile, $\rho_{h}({\mathbf{r}}|M)$, is the Navarro-Frenk-White (NFW)
profile, defined as (Navarro et al., 1997; Hayashi & White, 2008b)

$\displaystyle\rho_{h}(r|\rho_{s},r_{s})-\bar{\rho}=\frac{\rho_{s}}{(r/r_{s})(1+r/r_{s})^{2}},$
(17)

which depends on the two parameters scale density $\rho_{s}$ and scale radius
$r_{s}$. In order to obtain a finite halo mass $M=\int\rho_{h}(r)d^{3}r$, the
NFW profile, Eq. (17) needs to be truncated. Conventionally, the truncation
radius is set to the virial radius $R_{vir}$. Defining the concentration
parameter $c$ as

$\displaystyle c\equiv\frac{R_{vir}}{r_{s}},$ (18)

and noting

$\displaystyle\rho_{s}=\frac{c^{3}M}{4\pi
R_{vir}^{3}\left(\log(1+c)-c/(1+c)\right)},$ (19)

we can express the NFW profile in terms of the variables $M$, $R_{vir}$ and
$c$ as follows (Cooray & Sheth, 2002)

$\displaystyle\rho_{h}(r|\rho_{s},r_{s})\bigg{|}_{r<R_{vir}}\longrightarrow\;\rho_{h}(r|M,R_{vir},c)\bigg{|}_{r<R_{vir}}.$
(20)

With the halo density profile and the halo mass function defined, we have
constructed a general halo model.

####  2.3.2 Implementation in HMCode

We now describe HMCode , an implementation of the Halo Model by Mead et al.
(2015). The authors match the Halo Model parameters to a set of simulations
with varying cosmological parameters Lawrence et al. (2010b). Furthermore, to
improve the predictive power of the model, they introduce a range of
additional parameters which are fitted to these simulations. A summary of all
parameters as well as the best fit values reported by Mead et al. (2015) is
given in Table 1. Below, we explain the physical motivation as well as the the
effect of each parameter on the power spectrum. 666While this work was in
preparation, HMCode was updated Mead et al. (2021), however we kept our
analysis pipeline with the 2015 version.

For the Halo Model parameters halo mass $M$, concentration parameter $c$, and
critical density $\delta_{c}$, the authors use the following prescriptions:
The halo mass parameter $M$ is defined such that the average density of the
halo $\bar{\rho}_{h}(M)$ is greater than the background density $\bar{\rho}$
by a factor $\Delta_{v}=\bar{\rho}_{h}/\bar{\rho}$. $\Delta_{v}$ is fitted to
simulations according to the prescription in Table 1.

The concentration parameter $c$ is calculated according to the prescription of
Bullock et al. (2001)

$c^{\rm B}(M,z,P^{lin}(k))=A\frac{1+z_{\rm coll}(M,P^{lin}(k))}{1+z},$ (21)

where the minimum halo concentration parameter $A$ in the above equation is
again fitted to simulations (Mead et al., 2015). The redshift $z_{\rm
coll}(M)$ of collapse for halos of mass $M$ is defined as the redshift where
the mass of the halo was lighter by a factor of $f_{\rm coll}=0.01$ as
compared to the current mass at redshift $z$, estimated according to the
relation

$\displaystyle\frac{g(z_{\rm coll})}{g(z)}\sigma(f_{\rm
coll}M,P^{lin}(k))=\delta_{c},$ (22)

where $g$ is the linear growth factor. Thus, the concentration-mass relation
$c(M,z)$ is determined by the input linear theory power spectrum that is used
to compute $g$ and $\sigma$, which we denote explicitly in Eqs. (21)) and
(22)). The critical density for collapse, $\delta_{c}$, has the functional
form given in Table 1, with co-efficients determined by best fit to
simulations.

In addition to tuning the generic Halo Model parameters $M$, $c$ and
$\delta_{c}$ to simulations, Mead et al. (2015) introduce a range of
modifications, which are intended to parametrize physical effects which are
neglected in the basic Halo Model introduced in Section 2.3.1. These
modifications affect the one-halo term, the two-halo term and the combination
of both terms as described below.

HMCode approximates the two-halo term as given in Eq. (11), neglecting the
bias. This argument is confirmed by Mead et al. (2021), who find that
neglecting the linear bias, $b(M)$, in Eq. (11) only changes the total
correlation function at a significant level (compared to the desired accuracy)
once the one halo term is already dominant. Thus, $b(M)$, can be neglected
within the desired accuracy of HMCode . In fact, this approximation over-
predicts the correlation function on smaller scales of the order of the halo
radii. Hence, Mead et al. (2015) adopt an empirical damping model as follows

$\displaystyle\Delta^{2}_{2h}(k)$
$\displaystyle=\left[1-f\tanh^{2}(k\sigma_{v}/\sqrt{f})\right]\Delta_{lin}^{2}(k),$
(23) $\displaystyle\sigma^{2}_{v}$
$\displaystyle=\frac{1}{3}\int_{0}^{\infty}\frac{\Delta_{lin}^{2}(k)}{k^{3}}\mathrm{d}k,$
(24)

where $f$ fixed by an empirical fit, and $\Delta^{2}_{lin}(k)=4\pi
V(k/2\pi)^{3}P(k)$, where the power spectrum $P(k)$ is the Fourier transform
of the correlation function $\xi(x)$.

For the one-halo term, Mead et al. (2015) introduce two modifications. The
first modification is supposed to correct for the problem that the simple
expression Eq. (9) does respect that halos are mutually exclusive, which
results in too large contributions of the one-halo term to the correlation
function at very large scales. The correction is given by

$\displaystyle\Delta^{2}_{1h}=\left[1-e^{-(k/k_{*})^{2}}\right]\Delta^{2}_{1h,gen},$
(25)

where $\Delta^{2}_{1h,gen}=4\pi V(k/2\pi)^{3}P_{1h,gen}(k)$, with
$P_{1h,gen}(k)$ being the Fourier transform of the generic one-halo term
$\xi^{1h}(r)$ defined in Eq. (9). The parameter $k_{*}$ is again fixed by an
empirical fit. The authors furthermore modify the halo density profile in
Fourier space:

$\displaystyle\rho_{h}(k|M,R_{vir})\longrightarrow\;\rho(\nu^{\eta}k,M,R_{vir})$
(26)

where $\eta$ is again a heuristic fit parameter. With this modification of the
density profile, the total halo mass $M$ is unaltered. Furthermore, the
profile of halos with $\nu=1$ does not change, while for halos with $\nu>1$
($\nu<1$) the profile flattens (sharpens), as $\rho_{s}$ decreases (increases)
with $\nu^{-3\eta}$ while the scale radius $r_{s}$ and the cut-off increase
(decrease) as $r_{s}\,\nu^{\eta}$ and $R_{vir}\,\nu^{\eta}$, respectively.

Another shortcoming of the generic halo model is the crude transition between
the one-halo term appropriate for small scales and the two-halo term
appropriate for large scales. Physically, one expects in the so-called quasi-
linear regime a smooth transition. Mead et al. (2015) fit this smooth
transition regime with the following functional form

$\displaystyle\Delta^{2}(k)=\left[(\Delta^{2}_{2h})^{\alpha}+(\Delta^{2}_{1h})\right]^{1/\alpha},$
(27)

with $\alpha$ an empirical fit parameter given in Table 1, and
$\Delta^{2}_{2h}$ and $\Delta^{2}_{1h}$ defined according to Eqs. (23) and
(25), respectively. The fit parameter $\alpha$ is determined by the effective
spectral index, $n_{\rm eff}$, given by:

$3+n_{\rm eff}(z)=-\left.\frac{d\ln\sigma^{2}(R,z)}{d\ln R}\right|_{R=R_{\rm
nl}}\,,$ (28)

where $\sigma^{2}$ is the variance of the linear power spectrum, and the non-
linear scale $R_{\rm nl}$ is defined by $\sigma^{2}(R_{\rm nl})=1$. We discuss
the smoothing Eq. (27) and its role in FDM models in depth in Section 2.4.2.

Finally, Mead et al. (2015) fit the critical density $\delta_{c}$ the
virialised overdensity $\Delta_{v}=3m/4\pi R_{vir}^{3}\bar{\rho}$ used to
define $R_{vir}$, as well as the parameter $A$ used to define the
concentration parameter $c(A)$. See Table 1.

Parameter | Description |  Mead et al. (2015) fit value
---|---|---
$\Delta_{v}$ | Virialised halo overdensity | 418$\,\Omega_{m}(z)^{-0.352}$
$\delta_{c}$ | Linear collapse threshold | $1.59+0.0314\ln{\sigma_{8}(z)}$
$\eta$ | Halo bloating | $0.603-0.3\sigma_{8}(z)$
$f$ | Linear spectrum transition damping | $0.199\sigma_{8}(z)^{4.29}$
$k_{*}$ | One-halo damping wavenumber | $0.584\sigma_{v}(z)^{-1}$
$A$ | Min. halo concentration | $3.13$
$\alpha$ | Quasi-linear one- to two-halo term softening | $2.93\times 1.77^{n_{\mathrm{eff}}}$
Table 1: The halo model parameters and fitted values from Mead et al. (2015).
The additional non-linear nuisance parameters varied in this work are
specified in Eq. (36).

###  2.4 Halo Model and Lensing Observables for FDM

The lensing observables for FDM are shown in Figs. 12 and 13. In the following
we describe the modifications to the halo model that lead to the observed
effects, which cause ULAs to be observationally distinct from CDM. Effects we
did not include in our halo model, and which are deemed to be unimportant for
DES-Y1 observables, are discussed in the Supplementary Material.

####  2.4.1 Halo Mass Function and Concentration

It is known that for models with a truncated linear power spectrum,
$P^{lin}(k)$, such as models with FDM, but also with warm dark matter (WDM),
the standard predictions for $n(M,P^{lin}(k))$ and $c(M,P^{lin}(k))$ described
above disagree with the results of $N$-body simulations when using a real-
space top-hat window function. Specifically, in simulations $n(M)$ and $c(M)$
display additional suppression and a cut-off compared to the predictions from
the real space spherically averaged variance $\sigma(M,P^{lin}(k))$ computed
from the linear power spectrum (e.g Schneider et al., 2012; Schneider et al.,
2013; Schive et al., 2016; Corasaniti et al., 2017). The additional
suppression in $n(M)$ in an $N$-body simulation is visible after removing
numerical artefacts (“spurious haloes” Wang & White, 2007). However, such a
suppression is expected from basic physical principles applied to models with
a truncated linear power spectrum: we expect there to be no haloes below some
cut-off scale, and a component of unbound DM. Similarly, a turnover in $c(M)$
at low $M$ can be interpreted as a turnover in the redshift of collapse, with
low mass halos formed due to fragmentation of larger ones (see Eq. 21). HMCode
imposes a minimum value of $c=A>1$ at large masses, and we keep the same
minimum value also at low masses, implying $z_{\rm coll}>0$ for all halos.
This is consistent with the assumption that all halos are described by NFW
profiles, even if they are rare with low $n(M)$.

The exact shape of $n(M)$ and $c(M)$ below the cut-off is not well known,
since the spurious haloes are sometimes removed in a somewhat _ad hoc_ way.
This motivates exploring a range of models for the cut-offs, and marginalising
over this uncertainty when deriving cosmological constraints. We retain the
real space window function since it allows physical mass assignment, and is
already implemented in HMCode . We model the cut-offs in $n(M)$ and $c(M)$
within HMCode building on the results reported from $N$-body simulations, and
Marsh (2016a).

In Fig. 3 we show different model predictions for the HMF. In solid, black, we
show the Sheth-Tormen model $n^{\rm ST}(P_{\rm CDM})$, introduced in Eq. (13),
for a CDM linear power spectrum. In dashed, blue we show the Sheth-Tormen
model $n^{\rm ST}(P_{\rm FDM})$ for a FDM linear power spectrum. Below some
halo mass, the Sheth-Tormen model predicts a suppression of the $n^{\rm
ST}(P_{\rm FDM})$ with respect to $n^{\rm ST}(P_{\rm CDM})$, which is rather
uniform for varying $m$, reflecting the suppression of the linear power
spectrum. However, as discussed above, the prediction from the Sheth-Tormen
model, $n^{\rm ST}(P_{\rm FDM})$, deviates from simulation results which
feature an additional suppression of the HMF. To assess the impact on
prediction of the observables, we estimate the maximum visible halo masses
$M$. Assuming spherical collapse with a critical mass $\delta_{c}$ as
implemented in HMCode , we relate the the halo mass $M$ to a virial radius
$R_{\rm vir}$, which can be converted to an angular scale via the angular
diameter distance $D_{A}$. To make a conservative estimate we take the
redshift of $z=0.09$, below which are only 1% of the DES-Y1 galaxy sample in
the nearest redshift bin. We observe that differences between the CDM and FDM
HMF need to occur for relatively large halo masses $M$ to impact the DES-Y1
observables, implying that this effect will be only relevant for very low
particle masses. For consistency, we aim at matching the cut-off in the HMF
and in the concentration mass relation to be discussed below to simulation
results. Schive et al. (2016) showed that the suppression with respect to the
CDM case can be described quite accurately by the following two-parameter
model

$\displaystyle n(P_{\rm CDM})=n^{\rm ST}(P_{\rm CDM})\,\Delta_{n}^{\rm
CDM}(M_{0},\alpha_{1},\alpha_{2}),$ (29)

where $n(P_{\rm CDM})$ is the HMF taking into account the additional
suppression observed in simulations and $n^{\rm ST}(P_{\rm CDM})$ is the
original HMF for CDM in the Sheth-Tormen model. The correction function with
respect to the CDM case $n^{\rm ST}(P_{\rm CDM})$ is defined as

$\displaystyle\Delta_{n}^{\rm
CDM}(M^{(n)}_{0},\alpha_{1},\alpha_{2})=\left[1+\left(\frac{M}{M^{(n)}_{0}}\right)^{-\alpha_{1}}\right]^{-(\alpha_{2}/\alpha_{1})}$
(30)

with a scaling mass $M^{(n)}_{0}=1.6\times
10^{10}(m_{Ax}/10^{-22}\rm{eV})^{-4/3}M_{\odot}$. For halo masses $M\gtrsim
M^{(n)}_{0}$ the HMF agrees for FDM and CDM. For halo masses $M\lesssim
M^{(n)}_{0}$, the FDM HMF is suppressed with respect to the CDM case. The
steepness of the suppression is controlled by the parameter $\alpha_{2}$ and
the sharpness of the transition at $M\sim M_{0}$ is controlled by
$\alpha_{1}$. Fitting this model to $N$-body simulations results in
$\alpha_{1}=1.1$, $\alpha_{2}/\alpha_{1}=2.2$ (Schive et al., 2016). We show
this model in thin, black dashed in Fig. 3.

To implement the additional suppression into HMCode , we want to correct the
Sheth-Tormen model for FDM, $n^{\rm ST}(P_{FDM})$. We use the ansatz

$\displaystyle n(M,P_{\rm FDM})=n^{\rm ST}(P_{\rm FDM})\Delta_{n}^{\rm
FDM}(M^{(n)}_{0},\alpha_{1},\alpha_{2}),$ (31)

where we choose the correction function to be of the same functional form as
Eq. (30), and $n^{\rm ST}(P_{\rm FDM})$ is the Sheth-Tormen mass function
computed with the FDM linear power spectrum. This is justified because for
small halo masses $M\lesssim M^{(n)}_{0}$, the Sheth-Tormen prediction for
both, the CDM and the ULA case are well described by a power-law. Re-fitting
the slope-parameter $\alpha_{2}$ within our fiducial model, we find
$\alpha_{2}=1.86$ for our correction function in Eq. (31). We show our model
in solid, green in Fig. 3.

Figure 3: Suppression of the Halo Mass Function in our model. We show the
Sheth-Tormen prediction Eq. 12 (Sheth & Tormen, 1999) based on a CDM (black,
solid) and a FDM (orange, dashed) linear power spectrum, for a particle mass
$m=10^{-24}$ eV. In simulations, an additional suppression of the HMF is
observed, parametrized by Schive et al. (2016) (thin black, dashed). Our model
of the additional suppression is shown in green. To account for modeling
uncertainties we introduce two nuissance parameters $\alpha_{1}$, $\alpha_{2}$
(cf. Eq. 31). We show the effect of varying $\alpha_{1}$ ($\alpha_{2}$),
keeping $\alpha_{2}$ ($\alpha_{1}$) fixed in light green (dark green). The
grey-shaded regions correspond to Halo Masses whose virial radius cannot be
resolved on the scales present in the DES-Y1 data. Note that we pick different
values for the particle mass parameter $m$ in Fig. 3 and Fig. 4, because
differences between FDM and CDM affect the HMF and the concentration mass
relation at different scales. Figure 4: Suppression of the concentration mass
relation in our model. We show the case for CDM and FDM directly from the
linear power spectrum using the Bullock et al. (2001) relation in HMCode , Eq.
(21). Schneider et al. (2013) found additional suppression below the half-mode
mass in related WDM models. We include this extra suppression using two
nuisance parameters, $\gamma_{1}$ and $\gamma_{2}$. We show the effects within
each prior, holding the other one fixed, by the light green ($\gamma_{1}$) and
dark green ($\gamma_{2}$) shaded regions. We denote the mass scale $M$ above
which formally halos have not yet formed, i.e. $z_{\rm coll}<0$. For these
halos, which are very rare according to the respective HMF, we assume the
minimum halo concentration $A$. The vertical shaded region indicates
approximately halos that are below the DES resolution. Note that we pick
different values for the particle mass parameter $m$ in Fig. 3 and Fig. 4 to
highlight the differences between FDM and CDM affect the HMF and the
concentration mass relation, which occur on different scales.

Our approach for the concentration parameter is very similar to our treatment
of the HMF. In Fig. 4, we illustrate the different models for the
concentration parameter. In solid, black, we show the Bullock et al. (2001)
prediction, as introduced in Eq. (21), for a CDM linear power spectrum,
$c^{\rm B}(P_{\rm CDM})$. In thick blue, dashed we show the Bullock et al.
(2001) prediction for the FDM case, $c^{\rm B}(P_{\rm CDM})$, with the
modified collapse redshift, $z_{\rm coll}(P_{\rm FDM})$ computed from the FDM
linear theory power spectrum from axionCAMB . As in the case of the HMF,
$c^{\rm B}(P_{\rm FDM})$ is suppressed with respect to $c^{\rm B}(P_{\rm
CDM})$, turning flat for small halo masses $M$. However, the suppression
predicted by the Bullock et al. (2001) model is less than measured in
simulations. Note that we show our models for the concentration parameter for
a larger particle mass $m=10^{-23}$ eV as compared to Fig. 4. This is because
the differences between FDM and CDM affect the concentration parameter at
higher mass scales $M$, because the halo density profile, and hence the
concentration parameter, trace the matter density at the time of collapse.
Therefore, the concentration parameter is sensitive to a _higher_ particle
mass scale $m$.

We follow the approach of Schneider et al. (2012), who parameterise the
additional suppression of the concentration parameter within their WDM
simulations (which includes a free streaming scale very similar to FDM) as

$\displaystyle c(M,P_{\rm CDM}))=c^{\rm B}(P_{\rm CDM})\,\Delta_{c}^{\rm
CDM}(M^{(c)}_{0},\gamma_{1},\gamma_{2}),$ (32)

where $c(M,P_{\rm CDM})$ is the concentration parameter taking into account
the additional suppression observed in simulations and $c^{\rm B}(M)$ is the
original concentration parameter for CDM (Eq. 21) specified in Bullock et al.
(2001). Schneider et al. (2012) find that the two-parameter model

$\displaystyle\Delta_{c}^{\rm
CDM}(M^{(c)}_{0},\gamma_{1},\gamma_{2})=\left[1+\gamma_{1}\left(\frac{f\,M^{(c)}_{0}}{M}\right)\right]^{-\gamma_{2}}$
(33)

provides a good fit to their simulations.

In Eq. (33) the parameter $\gamma_{2}$ controls the slope of the cut-off. The
second parameter, $\gamma_{1}$, shifts the position of the cut-off scale. Note
the slightly different form of the suppression function $\Delta_{c}^{\rm CDM}$
as compared to the HMF suppression function $\Delta_{n}$. As above,
$M^{(c)}_{0}$ denotes the scale parameter. Note that we need to adapt a larger
scale mass $M^{(c)}_{0}=1/f\,M^{(n)}_{0}$ compared to the HMF, where the
parameter $f_{\rm coll}=0.01$ is used in the Bullock et al. (2001)
prescription implemented in HMCode . Therefore, to be consistent with
Schneider et al. (2012), we introduce an additional factor of $f$ in Eq. (33).
In thin, dashed, we show the prediction using the additional suppression
formula Eq. (33), with the best-fit parameters $\gamma_{1}=15$ and
$\gamma_{2}=0.3$. We see that the cut-off scale is moved by a factor
$\sim\gamma_{1}=15$ to the left compared to the scale $M^{(c)}_{0}$. To
implement this additional suppression within HMCode , we use the following
prescription

$\displaystyle c(M,P_{\rm FDM})=c^{\rm B}(P_{\rm FDM})\,\Delta^{\rm
FDM}_{c}(M^{(c)}_{0},\gamma_{0},\gamma_{1},\gamma_{2})$ (34)

with

$\displaystyle\Delta^{\rm
FDM}_{c}={\left[1+\left(\frac{M^{(c)}_{0}}{M}\right)\right]^{-\gamma_{0}}}\Delta_{c}^{\rm
CDM}(M^{(c)}_{0},\gamma_{1},\gamma_{2}),$ (35)

where $\gamma_{0}=\mathrm{d}\log c^{\rm B}(P_{\rm FDM})/\mathrm{d}\log M$,
where we evaluate the derivative at $4\cdot M^{(c)}_{0}$, slightly before
$c^{\rm B}(P_{\rm CDM})$ and $c^{\rm B}(P_{\rm FDM})$ start to deviate. The
new term in brackets in Eq. (35) counteracts the flattening of the
concentration curve in the Bullock et al. (2001) model for FDM, such that
without the suppression term $\Delta_{c}^{\rm CDM}$, we would approximately
recover the CDM curve, as demonstrated by the thin, green, dashed curve in
Fig. 4. The suppression must first be counteracted, then suppressed again,
since during a Monte Carlo parameter “scan” HMCode does not have access to the
CDM power spectrum and the CDM concentration parameter, but only the axion
concentration $c(M)$ (note that we drop the implicit dependence on
$P_{\mathrm{FDM}}$ for simplicity).

The suppression of $n(M)$ and $c(M)$ in Eqs. (31) and (34) is adopted to
account for a number of physical factors. As discussed above, $N$-body
simulations with a free streaming scale are better fit using a sharp-$k$
window function than a real space one (Schneider et al., 2013). The shape
including the sharp-$k$ function can be minimised by a real space space window
with an additional suppression, while retaining the advantage of a well
defined halo mass (Schive et al., 2016). HMCode has access to the linear
theory power spectrum from axionCAMB for FDM, and so the suppression function
must be modified to be with respect to $n^{\rm ST}(P_{\rm FDM})$. Such a
function can also account for the additional suppression of halo formation
caused by the Jeans scale and modified collapse barrier (Marsh & Silk, 2014;
Du, 2018). To account for the Jeans scale, the collapse barrier, $\delta_{c}$
should in principle also be modified, where pressure increases the overdensity
required for spherical collapse (as in the case of WDM velocities as shown in
Benson et al., 2013). A modified barrier for FDM was proposed, and implemented
in an approximate manner, in Marsh & Silk (2014); Marsh (2016a), while the
excursion set for the modified barrier was solved by Du et al. (2017); Du
(2018). A modified collapse barrier leads to additional suppression of the
mass function with respect to CDM, and also with respect to $N$-body
simulations with a free streaming scale. The halo mass function has not been
measured in simulations including the “quantum pressure” on all scales, and is
likely beyond the realm of present computational ability, although see e.g.
Schive et al. (2014a); Nori et al. (2019); Li et al. (2019); Mina et al.
(2022); Veltmaat & Niemeyer (2016); Mocz et al. (2019). Thus we adopt the
suppression function in Eq. (31) with the additional uncertainty shown in Fig.
3. On the other hand, Marsh (2016a) found that the modified barrier had little
effect on the non-linear power spectrum on observationally relevant scales.

In addition to these physical reasons to allow for variation in the fitted
cut-offs for $n(M)$ and $c(M)$, there is also unaccounted for possible
dependence on cosmological parameters, and on the method for removing spurious
haloes. The required simulations to account for variability in $n(M)$ and
$c(M)$ below the cut-off for ULAs are not available in the literature, and
performing such simulations is beyond the scope of the present work (and is
likely beyond present computational ability to resolve the required large and
small scales simultaneously). We thus adopt the following wide, flat, priors
on the nuisance fitting parameters:

$\displaystyle\alpha_{1}:\,$ $\displaystyle[0.6,2.0]$ (36)
$\displaystyle\alpha_{2}:\,$ $\displaystyle[1.43,2.54]$
$\displaystyle\gamma_{1}:\,$ $\displaystyle[5.0,45.0]$
$\displaystyle\gamma_{2}:\,$ $\displaystyle[-0.37,-0.23]$

We illustrate the variation induced in the FDM model sampling the nuisance
parameters from our priors in Figs. 3 and 4 by the green, shaded areas around
our fiducial suppression models described by Eqs. (31) and (34). Using this
method allows us in the following to present rigorous, Bayesian limits on the
axion mass from weak lensing, marginalising over theoretical uncertainty in
the underlying halo model.

In modelling the suppression of the HMF and the concentration mass relation,
we assume that the correction functions $\Delta^{\rm FDM}_{n}$ and
$\Delta^{\rm FDM}_{c}$, do not change with redshift. The fitting functions on
which we base our work from Schive et al. (2016) and Schneider et al. (2012)
are reported at $z=0$ and argued to be redshift independent.

For the HMF, the cut-off at large masses, when the asymptotic value of
$\sigma>1$, can be seen to be redshift independent by appeal to the known
redshift independence of the $P(k)$ cut-off (Hu et al., 2000), and origin of
the cut-off in the sharp-$k$ filtering model. We also compared our HMF cut-off
to the fit reported by Corasaniti et al. (2017): this shows very mild redshift
evolution within the uncertainties of the methods employed and accounted for
by different approaches to remove spurious halos. There is also redshift
dependence in the high mass end of the HMF: this is accounted for entirely in
our model by using the exact $P(k,z)$ for FDM from axionCAMB.

For $c(M,z)$ we investigated possible redshift dependence using the
alternative concentration calculation presented by Ludlow et al. (2016),
calibrated to simulations of warm DM. The method of Ludlow et al. (2016)
correctly reproduces the cut-off in $c(M,z)$ seen in simulations. We applied
this calculation to our FDM power spectra and observed no redshift dependence
in the location of the cut-off. In future, it would be highly desirable to
implement the method of Ludlow et al. (2016) into HMCode, but this is beyond
the scope of the present work, and we find the fits of Schneider et al.
(2012), with the adopted uncertainty, to be sufficient, and to have little
effect on our constraints.

####  2.4.2 One-to-Two Halo Term Smoothing

Figure 5: Evolution of $n_{\rm eff}$ as function of redshift $z$ as compared
to the CDM case, which is calibrated to simulations as described in Mead et
al. (2015, 2021). On the far left we show a histogram of the values of $n_{\rm
eff}$ determined on the nodes of the calibration emulator. In the top right
panel, we show in black the evolution of $n_{\rm eff}$ in the case of CDM for
our fiducial cosmology. Colour-coded by the mass $m_{\rm Ax}$ we show the
evolution of $n_{\rm eff}$ for the ULA case. In dashed we show models which
are strongly excluded by our analysis. In solid we show models which are
allowed or marginally excluded by our analysis. In the bottom panel, we show
$n^{i}_{\rm src}$, the normalised redshift distribution of the numbers of
galaxies in four redshift bins $i=1\dots 4$. We expect the bins to be affected
by capping over the redshifts that contain the most source galaxies. This
implies that while the parameter $n_{\rm eff}$ leaves the calibrated range
(shown by the green band in the top panels) at some $z$ for all our ULA
models, this does not necessarily impact the model prediction for all the
bins. As a guide to the eye, we mark the redshift at which a scale of $1\,\rm
Mpc$ near the quasi-linear regime is smaller than the smallest measured
angular bin at $2.8^{\prime}$ for our fiducial cosmology.

HMCode adopts a smoothing of the transition between the one halo term and the
two halo term (Eq. 27), which captures aspects of clustering in the quasi-
linear regime. The presence of this smoothing term leads to FDM displaying
_enhanced_ power, and thus shear correlation, relative to CDM over a narrow
range of masses near $m\approx 10^{-24}\text{ eV}$ (see Fig. 12, and also in
the the Supplementary Material, Fig. 13). This behaviour is counter-intuitive
from the point of view of the linear theory, where FDM only suppresses power.
An enhancement of power in the quasi linear regime caused by models with
suppressed linear power on small scales can, however, be explained by appeal
to the bias. Models with suppressed linear power generically have enhanced
bias on scales above the suppression scale (see e.g. Carucci et al., 2015;
Bauer et al., 2021; Laguë et al., 2022) due to conservation of mass. On
intermediate scales, this enhancement of the bias can compensate for the
linear theory suppression, leading to an overall enhancement of the power. It
is possible that HMCode indirectly accounts for this effect via the one-to-two
halo smoothing dependence on $n_{\rm eff}$ (see Eq. 28 and Table 1).

However, since the HMCode treatment is entirely phenomenological, we must
check to what extent the effect is properly calibrated for the FDM model, and
to what extent an erroneous or extrapolated treatment might influence our
reults. We investigate this in the following.

The underlying parameter, $n_{\rm eff}$, that controls the one-to-two halo
smoothing is calibrated to cosmological emulator nodes, where the power
spectrum is computed directly from $N$-body simulations. Varying the
cosmological parameters in these $N$-body simulations, can lead to a range of
values for the smoothing. This is illustrated by the yellow histogram in Fig.
5. The lowest value of $n_{\rm eff}$ simulated in the $N$-body simulations is
$n_{\rm eff}\approx-2.6$. Fig. 5 also shows the value of $n_{\rm eff}(z)$ for
a reference CDM cosmology (solid black curve), and for the same cosmology with
various ULA masses. When the ULA linear suppression scale is larger than the
non-linear scale (i.e. $k_{J,eq}<k_{nl}$), this leads to the asymptotic value
$n_{\rm eff}=-3$, which is outside the range of values calibrated by
simulations. We expect such models to be strongly excluded by DES, which
requires the formation of non-linear cosmological structure. 777Indeed, we
could even exclude such models by prior, since non-linear structure formation
is required for galaxy formation and the presence of life (Weinberg, 1987;
Tegmark et al., 2006).

Fig. 5 shows in the bottom panel the normalised distribution of source
galaxies for DES-Y1 as a function of redshift. The grey line indicates the
redshift $z\approx 0.5$ at which a scale of $1\,\rm Mpc$ is less than
$2.8^{\prime},$ the smallest angular bin in the DES-Y1 data. The value of
$n_{\rm eff}(z)$ computed from the FDM linear power spectrum remains fully
within the calibrated range for DES source galaxies with redshifts $z>0.5$ for
$m\gtrsim 10^{-23}\text{ eV}$, which, as we will see, covers all of the 95%
confidence region of our constraints. Furthermore, the parameter $n_{\rm eff}$
impacts the quasi-linear scale. For higher redshifts, this scale cannot be
resolved within the range of measured angular sizes. For higher redshifts, the
effect of $n_{\rm eff}$ on the model of the observables is expected to be
smaller, because more of the linear part of the power spectrum and less of the
quasi-linear contribute to the observables.

To further investigate the effect on cosmic shear observables of $n_{\rm
eff}(z)$ leaving the calibrated range, we computed the effect on the shear
correlation function $\xi^{+}$ in the $i$-$i$-bins for two possibilities for
$n_{\rm eff}$ (and thus smoothing parameter $\alpha$), shown in Fig. 6. The
first takes $\alpha$ at its value specified by $n_{\rm eff}$ according to the
fits in HMCode . The second possibility we investigate applies a hard cap to
$\alpha$, forcing it to lie in the range covered by the emulator nodes. We
show the relative difference between the two treatments compared to the data
error for an axion mass of $\sim 3\times 10^{-24}\,$eV (this mass is strongly
excluded by our analysis, while the enhanced shear correlation caused by the
smoothing parameter is maximal). As illustrated by Fig. 6, the difference
between the two approaches (capping or not capping the smoothing parameter) is
small compared to the uncertainty in the data.

For higher values of $m$, the difference between capping or not capping
becomes smaller, and for $m>10^{-24}\text{ eV}$ is less than about 10% the
size of the DES error. Since our lower limit to $m$ will turn out to be an
order of magnitude larger than this, we conclude that the effect of $n_{\rm
eff}$ leaving the calibrated range of HMCode on our results is likely small,
and will not affect our conclusions. Nonetheless in the analysis in Section 4,
we consider three possibilities for $n_{\rm eff}$: 1) leaving $n_{\rm eff}$
_free_ , as implemented in HMCode , 2) _capping_ $n_{\rm eff}$ to lie in the
calibrated region, and 3) interpolating between those two possibilities by
imposing an additional parameter which we _marginalise_ over.

We can thus conclude that, for the purposes of limiting the FDM mass while it
is fixed to be all of the dark matter, we can safely take the prescription for
one-to-two halo smoothing as specified by HMCode , since the power spectrum
slope falls within the range of the cosmologies covered by the emulator nodes
for FDM masses within the bulk of the posterior probability, and the effect of
capping the smoothing parameter has an effect on the correlation function
smaller than the observational error bars. Nonetheless, the calibration of
HMCode for beyond CDM models in the quasi-linear regime is an important open
problem and should be investigated further using simulations. It will be
particularly important to search for evidence of sub-dominant ULAs with
$m<10^{-23}\text{ eV}$ using high precision lensing observables from e.g.
_Euclid_ (Marsh et al., 2012; Amendola et al., 2018).

Finally, we note that the smoothing term adopted in HMCode is not present in
the comparably accurate halo model/perturbation theory hybrid of Sullivan et
al. (2021). The hybrid model does, however, contain elements that might
account for the same features in $P(k)$ captured by the HMCode smoothing term.
In the hybrid model the two halo term is modified to follow the Zel’dovich
approximation (Lagrangian perturbation theory), while the one halo term is
modified to include a “broadband expansion” of the halo profile, which has
free coefficients to be fitted to simulations, and can be used to model
baryonic feedback. The hyrbid model can account for enhancement of power in
the quasi-linear regime (as caused by HMCode smoothing when $n_{\rm eff}$ is
small) via a polynomial transfer function with further fitting parameters
augmenting the Zel’dovich approximation. The coefficients of all these fits to
simulations are found to vary with cosmological parameters, but Sullivan et
al. (2021) do not vary $\sum m_{\nu}$ or $n_{s}$, both of which would affect
$n_{\rm eff}$ dependence of the fitting parameters in the transition region.
Thus it is not possible to assess whether the hybrid model would similarly
predict the observed effects of FDM in HMCode based on the fits presented thus
far by Sullivan et al. (2021). It would be useful in further application of
the hybrid model to investigate parameter dependences of the fitting based on
model independent phenomenological parameters such as $n_{\rm eff}$, which can
be applied to any input linear power spectrum.

Figure 6:  Discrepancy in the model predictions for $\xi_{+}$ (solid lines)
and $\xi_{-}$ (dashed lines) for different treatments of the $n_{\rm eff}$
parameter, relative to the error (square-root of the diagonal entries of the
covariance matrix). We show the diagonal bins for three different axion masses
below our limit. We compare the case where we use $n_{\rm eff}$ as originally
implemented in HMCode , which will leave the CDM uncalibrated region for some
$z$ (see. Fig. 5 and text for details) and the case when $n_{\rm eff}$ is
capped to remain within the CDM calibrated region. The grey band (grey, dashed
line) shows $\xi_{+}$ ($\xi_{-}$) data points that are excluded from the
analysis by the small scale cut imposed by the DES collaboration. For axion
masses at and reasonably below our bounds, represented by the green curve, the
effect of capping versus not capping is small relative to the error in the
data. For $m$ far smaller than our lower limit, the effect becomes stronger,
with the biggest impact on the 3-3 bin. However, these masses are in any case
excluded.

##  3 Data and Methodology

###  3.1 Statistical Methods

We perform Bayesian analysis using the publicly available code cosmosis to
implement the DES-Y1 and _Planck_ likelihoods. We use multinest (Feroz et al.,
2009) as implemented within cosmosis to derive parameter constraints. We used
a stopping criterion of $\log{Z}=10^{-5}$ and 2100 (3500 for the full “3x2pt"
analysis, see the Supplementary Material) livepoints to ensure smooth
confidence intervals from the multinest chains.

On a subset of runs, we also compared the multinest constraints with Markov
Chain Monte Carlo (MCMC) analysis using the affine invariant emcee (Foreman-
Mackey et al., 2013) sampler as a cross-check. We used these tests to also
perform convergence tests of our result using the spectral method (Dunkley et
al., 2005), and to compute the autocorrelation time of our chains.

The autocorrelation method is useful for chains with multiple walkers, because
the error on the MCMC estimate for a given parameter is proportional to
$\tau/N$ where $N$ is the total number of samples (i.e. the number of samples
per walker times the number of walkers) and where $\tau$ is the integrated
autocorrelation time (Goodman & Weare, 2010). We found that the standard
autocorrelation time convergence criterion described in the emcee
documentation to be overly stringent for our cases with a large number of
parameters and walkers, and time consuming likelihood and theory calculations.
Specifically, finding a reliable estimate of the autocorrelation time would
have required a prohibitively large number of steps and computation time. We
thus used to an autoregressive model to get an estimate of the autocorrelation
time with fewer samples. After verifying that the autoregressive model did not
underestimate the autocorrelation time, we ensured that all the chains
satisfied $N\gg\tau$ using that estimate. An example comparison of the
autocorrelation time estimates from the three methods described above can be
found in Fig. 7.

Figure 7: Estimation methods for the chain autocorrelation times. We show that
the autoregressive model approach gives a more stable estimate of the
autocorrelation time with short chains. The fact that the autoregressive model
slightly overestimates the autocorrelation time means that this method to
evaluate convergence is conservative since more samples are needed to reach
$N\gg\tau$.

###  3.2 DES Data Sets

Our analysis is based on the DES-Y1 3x2pt data set consisting of the
observables $\xi_{\pm}$, $\gamma_{t}$ and $w$ described in Section (2.1).
Consistent with the discussion there, we focus on the cosmic shear data in our
main analysis, because it does not depend on the galaxy bias and is therefore
less sensitive to modelling uncertainties. The DES data on the cosmic shear
observables in the $4,4$ bin was shown already in Fig. 2 illustrating the
effect of the FDM mass holding all other parameters fixed. For the respective
figures the other bin combinations as well as the figures of $\gamma_{t}$ and
$w$ we refer to the Supplementary Material.

We validated our pipeline for the DES-Y1 data with a series of cross checks,
described briefly here. First, we investigate the effects of our non-linear
model on CDM by comparing our constraints on the cosmological parameters using
HMCode and the alternative non-linear model Halofit (Smith et al., 2003),
using both 3x2pt and shear only, finding no significant differences. As noted
in Hložek et al. (2017), Halofit cannot be applied to ULAs/FDM. Next, we
compared the results of our CDM analysis to those of the publicly available
chains from DES 888http://desdr-
server.ncsa.illinois.edu/despublic/y1a1_files/chains/, finding visual
agreement in the two-dimensional posteriors.

As discussed by Joudaki et al. (2020), the marginalised limits on the main
parameters constrained by weak lensing, ($\sigma_{8}$,$\Omega_{m}$), are prior
dependent: the upper limit on $\sigma_{8}$ depends on how $A_{s}$ is sampled,
and on the assumptions made about neutrino mass. Furthermore, the limit on
$\Omega_{m}$ relies on the 3x2pt data combination, and thus depends on the
galaxy bias model. Using weak lensing as an anchor to the non-linear scales to
constrain DM thus requires further inputs for the cosmological parameters.

###  3.3 Planck Data and Data Combination

The full _Planck_ data set of two point statistics includes temperature $T$,
$E$ and $B$ mode polarisation, and lensing deflection $\phi$ power spectra and
cross-spectra. The _Planck_ lensing deflection is correlated with DES lensing
observables. We do not model this correlation, and so we omit $\phi$ auto and
cross spectra from our _Planck_ data. Small angle polarisation $B$ modes are
generated primarily from lensing, and so we omit these also from our analysis.
We neglect the covariance between the _Planck_ $TT$ angular power spectrum and
DES weak lensing arising from the lensing-induced smoothing of small-scale
$TT$ peaks and troughs. We consider only the adiabatic mode of initial
conditions, as described in Hlozek et al. (2015), omitting isocurvature modes.
A complete analysis of _Planck_ data in the ULA model including lensing and
isocurvature was performed by Hlozek et al. (2018). For the particle masses of
interest, omitting isocurvature modes is equivalent to a prior on the Hubble
scale during inflation $H_{I}\lesssim 10^{12}\text{ GeV}$ (Marsh, 2016b;
Hlozek et al., 2018).

In our analysis we combine the _Planck_ CMB data with DES-Y1 for both shear
only and 3x2pt. In this combination, the _Planck_ CMB data constrains the six
standard cosmological parameters extremely well, and for CDM the effect of
DES-Y1 is only to slightly tighten the error bars by relatively small amounts.
Applied to constraining the nature of DM, however, the utility of the
combination is greater: _Planck_ anchors the large scales and the cosmological
parameters, while DES anchors smaller scales and thus constrains the nature of
DM further.

Parameter | Prior range
---|---
DM particle mass $\log{m}/\mathrm{eV}$ | $\mathcal{U}[-25.5,-22]$
DM density $\Omega_{d}$ | $\mathcal{U}[0.1,0.9]$
Hubble parameter $H_{0}$ [km/s/Mpc] | $\mathcal{U}[55,90]$
Baryon density $\Omega_{b}$ | $\mathcal{U}[0.03,0.07]$
Scalar spectral index $n_{s}$ | $\mathcal{U}[0.87,1.07]$
Scalar amplitude $A_{s}$ | $\mathcal{U}[5e^{-10},5e^{-9}]$
Neutrino density $\Omega_{\nu}h^{2}$ | $\mathcal{U}[0,0.01]$
Table 2: Cosmological parameters varied in the analysis.

###  3.4 Parameters, Priors, and Models

We adopt primary cosmological parameters:

$\\{A_{s},n_{s},H_{0},\Omega_{b},\Omega_{m},\sum m_{\nu}\\}\,,$ (37)

where $A_{s},n_{s}$ are the primordial power spectrum amplitude and spectral
index, $H_{0}$ is the present day Hubble parameter, $\Omega_{b}$ is the baryon
density, and $\Omega_{m}$ is the total matter density,
$\Omega_{m}=\Omega_{b}+\Omega_{d}$ and the DM ($\Omega_{d}$) is composed of
either CDM (for test cases) or FDM, but _not_ both. Massive neutrinos are
included with the parameter $\sum m_{\nu}$ for which we assume a single
massive neutrino and $N_{\rm eff}=2.04$ massless neutrinos. The neutrino mass
is related to the neutrino density parameter by $\Omega_{\nu}h^{2}=\sum
m_{\nu}/94.\text{ eV}$. For detailed discussion of how massive neutrinos are
treated in the halo model, we refer to Mead et al. (2015, 2021).

The parameters are varied with flat priors as specified in Table 2. These
priors are consistent with the choices made by the DES collaboration in their
own analysis. These are different from the priors adopted in standard CMB
analyses, which use a log prior on $A_{s}$, vary the physical densities
$\Omega_{c}h^{2}$ and $\Omega_{b}h^{2}$, and $H_{0}$ is a derived parameter
from the angular scale of the first acoustic peak. The parameter choices of
Eq. (37) are not optimal for a CMB analysis. When analysing galaxy data only,
Joudaki et al. (2020) showed how these different priors affect the resulting
bounds on $S_{8}$ and $\Omega_{m}$. When we combine DES-Y1 with _Planck_ , the
prior dependence largely vanishes due to the strong constraining power of
_Planck_ for the power spectrum amplitude and matter density. When using the
_Planck_ data, we vary the optical depth $\tau$ with a uniform prior, and the
_Planck_ absolute calibration parameter, $A_{\rm Planck}$, with a Gaussian
prior.

The DES-Y1 analysis contains a large number of nuisance parameters, which are
listed in Table 3. These are varied, and marginalised over in any data
combination including DES.

Parameter | Sampling range | Prior
---|---|---
Bin Bias parameters
$b_{1}$ | $\mathcal{U}[0.8,3.0]$ | $\mathcal{U}[0.8,3.0]$
$b_{2}$ | $\mathcal{U}[0.8,3.0]$ | $\mathcal{U}[0.8,3.0]$
$b_{3}$ | $\mathcal{U}[0.8,3.0]$ | $\mathcal{U}[0.8,3.0]$
$b_{4}$ | $\mathcal{U}[0.8,3.0]$ | $\mathcal{U}[0.8,3.0]$
$b_{5}$ | $\mathcal{U}[0.8,3.0]$ | $\mathcal{N}[0.8,3.0]$
Shear calibration parameters
$m_{1}$ | $\mathcal{U}[-0.1,1]$ | $\mathcal{N}[0.012,0.023]$
$m_{2}$ | $\mathcal{U}[-0.1,1]$ | $\mathcal{N}[0.012,0.023]$
$m_{3}$ | $\mathcal{U}[-0.1,1]$ | $\mathcal{N}[0.012,0.023]$
$m_{4}$ | $\mathcal{U}[-0.1,1]$ | $\mathcal{N}[0.012,0.023]$
Intrinsic alignment parameters parameters
$a$ | $\mathcal{U}[-5,5]$ | $\mathcal{U}[-5,5]$
$\alpha$ | $\mathcal{U}[-10,10]$ | $\mathcal{U}[-10,10]$
Weak lensing photo-$z$ error parameters
$wl-pz_{1}$ | $\mathcal{U}[-0.1,1]$ | $\mathcal{N}[-0.001,0.016]$
$wl-pz_{2}$ | $\mathcal{U}[-0.1,1]$ | $\mathcal{N}[-0.019,0.013]$
$wl-pz_{3}$ | $\mathcal{U}[-0.1,1]$ | $\mathcal{N}[0.009,0.011]$
$wl-pz_{4}$ | $\mathcal{U}[-0.1,1]$ | $\mathcal{N}[-0.018,0.022]$
Lensing photo-$z$ error parameters
$l-pz_{1}$ | $\mathcal{U}[-0.05,0.05]$ | $\mathcal{N}[0.008,0.007]$
$l-pz_{2}$ | $\mathcal{U}[-0.05,0.05]$ | $\mathcal{N}[-0.005,0.007]$
$l-pz_{3}$ | $\mathcal{U}[-0.05,0.05]$ | $\mathcal{N}[0.006,0.006]$
$l-pz_{4}$ | $\mathcal{U}[-0.05,0.05]$ | $\mathcal{N}[0.00,0.01]$
$l-pz_{5}$ | $\mathcal{U}[-0.05,0.05]$ | $\mathcal{N}[0.0,0.01]$
CMB Parameters
$\tau$ | $\mathcal{U}[0.02,0.2]$ | $\mathcal{U}[0.02,0.2]$
$A_{\rm Planck}$ | $\mathcal{U}[0.9,1.1]$ | $\mathcal{N}[1,0.0025]$
Table 3: The nuisance parameters for the DES-Y1 and Planck18 dataset, and
their sampling priors.

The axion mass is varied as:

$\mathcal{P}(\log_{10}[m/\text{eV}])=\mathcal{U}[-25.5,-22]\,.$ (38)

Such a log-flat prior, while arguably less informative than alternatives, must
have upper and lower limits to keep the prior volume finite. The lower bound
in Eq. (38) is motivated by the CMB constraints varying the axion fraction
(Hlozek et al., 2018), who find that at our lower limit the CMB forbids axions
from being all the DM at high significance. Lower masses have vanishing
posterior weight, and our exclusion limit is not affected by the lower limit
of the prior extending to smaller values. Since we expect a one sided lower
limit to $\log_{10}(m/\text{eV})$, however, the upper bound of the prior will
affect the exclusion. In such a case, the upper bound of the prior should be
chosen to be close to the expected experimental sensitivity, which can in
principle be computed _a priori_ knowing only the experimental specifications
(e.g. from a Fisher matrix forecast).

We set the upper limit of our prior following the discussion in the
introduction, as follows. We wish to set our prior upper limit to be
equivalent to CDM, but to establish in a data-driven manner. Chabanier et al.
(2019) project the DES-Y1 data into wavenumber bins on the linear power
spectrum, the locations of which can be computed without specifying the
measurement value in the bin. The highest bin for DES-Y1 covers
$1\,h\,\text{Mpc}^{-1}<k<5\,h\,\text{Mpc}^{-1}$. Setting $k_{\rm
J,eq}(m)=5\,h\,\text{Mpc}^{-1}$ in Eq. (1) gives $m=1.5\times 10^{-23}\text{
eV}$ as the approximate DES sensitivity to $m$. Due to the rough nature of
this estimate, we choose our maximum axion mass to be an order of magnitude
larger than the expected sensitivity. A log-flat prior on $m$ was also used in
Rogers & Peiris (2021b), where the lower limit was also set according to
existing bounds, and the upper limit was set according to the smallest scale
(i.e. heaviest axion mass) to which the Lyman-$\alpha$ forest data were
sensitive.

Figure 8: Posterior density for the DM particle mass $m$, the mass fluctuation
parameter $S_{8}$, and the matter density $\Omega_{m}$ for different
combinations of data sets.

The disadvantage of a log-flat prior is that there is infinite prior volume
between any chosen upper limit and the CDM limit $m\rightarrow\infty$. The
alternative prior:

$\mathcal{P}(10^{-23}\text{ eV}/m)=\mathcal{U}[0,200]\,,$ (39)

has finite prior volume up to the CDM limit, $\mathcal{P}(0)$. Such a prior
has been considered constraining the dark matter nature from the Lyman-alpha
forest (e.g. Iršič et al., 2017a). However, this inverse mass prior requires
the choice of a reference mass scale set by the lower limit and, unlike the
log-flat prior, it does not apply equal prior probability to different
logarithmic particle mass scales within the limits set by existing bounds and
projected sensitivity. In the following, our main analysis uses the log-flat
prior Eq. (38), but we also present a comparison to the case Eq. (39) for a
subset of our analyses.

We adopt various nuisance parameters to account for uncertainty in the non-
linear model of the HMF and halo concentration, as given in Eq. (36). In a
subset of our analysis we also investigated a nuisance parameter for the halo
model smoothing parameter, $\alpha$, as described in Section 2.4.2.

##  4 Results

Our baseline analysis uses the combination of _Planck_ $TT$, $TE$, and $EE$
spectra and DES-Y1 shear. Fig. 8 shows marginalised two dimensional posteriors
on the FDM mass and cosmological parameters for our baseline analysis. We show
the total matter density $\Omega_{m}=\Omega_{d}+\Omega_{b}$, and the “matter
fluctuation parameter” $S_{8}=\sigma_{8}\,(\Omega_{m}/0.3)^{0.5}$. These are
the parameters of most interest for weak lensing surveys (Abbott et al.,
2018a), and the other parameters are well constrained in combination with
_Planck_.

For _Planck_ alone, there is a slight degeneracy between $S_{8}$ and $m$ at
very low masses, due to the suppression of the linear power caused by FDM.
Fig. 8 demonstrates that, due to the tight constraints on the primary
parameters provided by the _Planck_ data, there is no degeneracy between the
ULA mass and other cosmological parameters at higher mass (the same is true of
the cosmological parameters that are not shown in this figure). For DES alone
there are stronger degeneracies, and wider limits to $S_{8}$ and $\Omega_{m}$.

Compared to _Planck_ alone, adding in the DES-Y1 data provides an anchor to
measure the power spectrum into the non-linear regime, leading to an improved
lower bound to the FDM mass:

$\displaystyle\log_{10}(m/\text{eV})$ $\displaystyle\geq-23.0\quad\text{(95\%
C.L.)}$ $\displaystyle\Rightarrow m$ $\displaystyle\geq 1\times 10^{-23}\text{
eV}\,.$ (40)

The marginalized posterior distributions for the axion mass for Planck and DES
are shown in Fig. 1. Both probes show an agreement with $m>10^{-22}$ eV, and
are therefore consistent with DM being cold in a model allowing FDM. DES alone
has a peak probability at $m\approx 10^{-22.5}$ eV which we attribute to the
lower values of $S_{8}$ and $\Omega_{m}$ measured by this survey in isolation.
In the posterior for DES alone, we observe a peak probability at $m\approx
10^{-22.5}$ eV which we attribute to a multinest sampling artefact, which has
no effect on the value of the exclusion limit (the limits agree between
multinest and emceee and no corresponding peak is present with emcee).

Figure 9: Cosmological parameters and axion mass posteriors with three
different treatments of the two halo smoothing parameter, $n_{\rm eff}$, for
_Planck_ +DES-Y1 shear (see text for description). Consistency between the
treatements demonstrates that the additional effects caused by $n_{\rm eff}$
leaving the range calibrated to simulation do not affect our conclusions
related to the cosmological parameters and FDM mass. The vertical dashed lines
show the best-fit parameter values from the published Planck TT,EE,TE+lowE
results Aghanim et al. (2020b), and the 95% upper bound for the neutrino mass.
Shifts in the central values compared to _Planck_ , though consistent at
$1\sigma$ (shaded region), are caused by our adopting priors consistent with
the DES-Y1 analysis, but different from the _Planck_ analysis (see discussion
in Joudaki et al., 2020).

As discussed in Section 2.4.2, the two-halo smoothing term in HMCode is
controlled by a parameter, $n_{\rm eff}$, which for FDM leaves the range
calibrated to simulation. To estimate the dependence of our cosmological
results on the smoothing term, we ran three separate analyses treating $n_{\rm
eff}$ differently: one in which it was _free_ to take the value outside of the
range calibrated through simulations and predicted by the FDM linear power
spectrum, one in which $n_{\rm eff}$ was _capped_ to remain within the
calibrated range, and one in which it was _marginalised_ with an additional
nuisance parameter if it left the calibrated range.

We compared the maximum likelihood in a range of mass bins for the different
models for $n_{\rm eff}$. As anticipated in our discussion in sec. 2.4.2, we
find that the different models for $n_{\rm eff}$ impact the likelihood at
small masses to a small extent. At higher masses, near to where we set our
exclusion limits, the effect is negligible for the combined DES and Planck18
likelihood and very small for the DES likelihood alone. We then proceeded to
check the effect of the three different treatments of $n_{\rm eff}$ on the
posterior. For DES alone, we found a weak dependence on the $n_{\rm eff}$
parameter in the 95% region of the two dimensional marginal posteriors, which
is not present in the combined analysis. The one-dimensional posteriors for
_Planck_ +DES-Y1, for all three treatments, are shown in Fig. 9. These results
are consistent between all three treatments, demonstrating that our
conclusions are not sensitive to the treatment of the smoothing term outside
the range calibrated to simulation.

Figure 10: Cosmological parameters and axion mass posteriors with two
different treatments of the halo mass function, $n(M)$, and concentration mass
relation, $c(M)$, cut-offs (see Figs. 4 and 3). In the upper panels, the cut-
offs are fixed to their fiducial values following simulations in Schneider et
al. (2012); Schive et al. (2016), while the lower panels additional
uncertainty in the low mass shape is marginalised. Consistency between the
treatments demonstrates that the additional effects caused by uncertainty in
$n(M)$ and $c(M)$ below the cut-off scale do not affect our conclusions
related to the cosmological parameters and FDM mass. The vertical dashed lines
show best-fit values from the Planck (TT,EE,TE+lowE) data combination Aghanim
et al. (2020b), and the 95% confidence limit for the neutrino mass. Shifts in
the central values compared to _Planck_ , though consistent at $1\sigma$
(shaded region), are caused by our adopting priors consistent with the DES-Y1
analysis, but different from the _Planck_ analysis (see discussion in Joudaki
et al., 2020).

In addition, Fig. 9 also demonstrates the consistency between our results with
FDM and the _Planck_ analysis, which assumes CDM. We also verified this
independently with our own CDM analysis, finding that the marginalised
posteriors on the cosmological parameters align for CDM and ULAs. This
implies, qualitatively, that there is little to no effect on $\sigma_{8}$ and
$H_{0}$ tensions including FDM as all the DM in the mass range we have
considered.

Similarly to the uncertainty in the halo model smoothing parameter, we found
that our analysis was also insensitive to the other nuisance parameters listed
in Eq. 36, which are introduced to account for uncertainty in the HMF and
$c(M)$. The nuisance parameters introduced to account for this uncertainty
show no degeneracy with any cosmological parameters, in particular they are
not degenerate with the FDM mass. The marginalised posterior on the axion mass
does not change when these parameters are varied compared to the case when
they are held fixed, Fig. 10.

From these runs varying the parameters controlling the uncertainty on the FDM
non-linear power, we conclude that these uncertainties lead to no significant
effect on the FDM mass constraint derived from DES-Y1 weak lensing
measurements at the current level of precision in the data.

Figure 11: Prior dependence. Purple histograms show the case of a log-flat
prior on the particle mass, while green histograms show the case of a uniform
prior on the inverse mass. The left panel shows the resulting posterior in
log-space, while the right panel shows the posterior in $1/m$ space. In each
case we find a lower/upper limit to the variable in the space over which the
prior was uniform, but we find a different shaped posterior after a change of
variables. For each prior, we mark the 95% limit by a vertical line to
illustrate the noticeable difference in the derived mass bound. We also show
as un-filled histograms (solid lines) the case of binning our chains in the
original parameter and transforming the histogram bounds to the complementary
parameter. After this rescaling (effectively a change of prior), these
transformed histograms agree very well with the posteriors obtained from the
complementary prior. This illustrates that the different results for the mass
parameter $m$ are indeed caused by the different priors and are not a sampling
artefact or an effect in the likelihood. This illustrates that care must be
taken to include knowledge of the prior volume when producing mass limits
after a changeof variable.

Fig. 11 demonstrates the effect on the posterior distribution from the prior
on the FDM mass (see § 3.4 for more discussion). Our multinest chains
correspond to samples drawn from the posterior distribution of the parameters,
in particular the mass parameter $m$, sampled in $\log_{10}(m/\text{eV})$ and
$10^{-23}\,\text{eV}/m$, respectively. When we calculate the 95% percentile of
these distributions, we find a mass bound of $m=9.2\times 10^{-24}\,$eV for
the log-flat prior and $m=6.7\times 10^{-24}\,$eV for the prior flat in $1/m$
when transformed to limits in $m$ itself. We can see how this difference comes
about by binning our chains evenly in both mass parameters, $\log(m)$ and
$1/m$, where we convert from one mass parameter to the other _before_ binning
the samples from our chains. The resulting histograms, illustrated by filled
histograms in Fig. 11, differ notably when a change of variables is performed.
In particular, we notice a peak in the posterior around $10^{-23}\text{ eV}$
when transforming the $1/m$ prior to log space. This is caused entirely by the
choice of scale inherent in the $1/m$ prior, and is not a feature of the data.

To demonstrate that this is caused by the choice of prior, rather than a
feature in the data/likelihood or a sampling artefact, we furthermore evenly
bin our chains in the respective original mass parameter and _subsequently_
convert the resulting histogram bounds to the complementary mass parameter.
This process is equivalent to importance re-sampling, where we would re-weight
each element of the Multinest chain by the ratio of the respective priors. For
our special case of two flat priors, the respective posteriors agree up to
normalization within the range of the two priors, therefore we can directly
transform the histogram, which represents our estimate of the posterior.
Indeed, we find that after suitable normalisation, the resulting histogram
agrees well with the posterior distribution obtained from the complementary
parameter. This illustrates that there are no features in the data leading to
the different shapes of distribution. In a Bayesian analysis one must always
present limits to the parameter which is varied and be careful comparing
results using different priors. Our log-uniform prior selects a scale only in
the data-driven upper limit, as discussed above, and therefore has physical
motivation.

##  5 Discussion and Conclusions

###  5.1 Summary

We have developed a rigorous pipeline to analyse the effects of FDM on the
DES-Y1 data, and thus derived a new lower bound on the FDM particle mass
combining DES-Y1 and _Planck_ data. Our pipeline uses likelihoods and samplers
from cosmosis . The linear power spectrum was computed with axionCAMB , while
the non-linear power was modelled using the halo model as implemented in
HMCode . HMCode is calibrated to simulations of CDM, including effects that
lead to moderate suppression of the linear $P(k)$ on small scales, thus
capturing some of the effects of FDM. We modify the halo model for FDM by
introducing additional cut-offs in the HMF and halo concentration, seen in
$N$-body simulations of FDM. Uncertainties on the non-linear model are
propagated through our analysis pipeline by marginalising over nuisance
parameters. We found that the additional uncertainty on the non-linear model
of FDM did not affect our results, since the additional parameters of the
model are uncorrelated with the cosmological parameters, in particular to the
FDM mass, and are unconstrained by the data at the current level of precision.

We also presented a bound on the FDM mass from the CMB alone, demanding that
FDM is all of the DM, we find $\log_{10}(m/\text{eV})\geq-24.6\quad\text{(95\%
C.L.)}$ The combination of DES-Y1 and _Planck_ data leads to an improvement of
this bound by approximately two orders of magnitude with respect to the CMB
alone: $\log_{10}(m/\text{eV})\geq-23.0\quad\text{(95\% C.L.)}$ The
improvement is driven by the fact that DES data on small angular scales probes
the cut-off in the linear power spectrum induced by FDM on scales inaccessible
to the CMB. The CMB plays the role of fixing the primary cosmological
parameters on large scales, leaving no room for degeneracy between the FDM
induced cut-off and e.g. the power spectrum amplitude or spectral index in
their effects on the small scales probed by DES. These results are shown in
Fig. 1.

A striking result of HMCode applied to FDM is that the one-to-two halo
smoothing parameter leads to FDM _enhancing_ the power relative to CDM over a
small range of scales in the quasi-linear regime and particle masses around
$10^{-24}\text{ eV}$. We checked that this effect is largely within the
calibrated range of the linear transfer function parameters of HMCode , within
the redshift range probed by DES galaxies. When the FDM effect leaves the
calibrated range, we tested that this does not drive constraints by comparing
a variety of treatments for the smoothing parameter, and finding that our
bound on the particle mass was unaffected. The increase in power can be
explained qualitatively by noting that halo bias in models with reduced linear
power increases on scales near the suppression scale, leading in some cases to
slightly increased non-linear power on larger scales (“moving power around”).
HMCode does not directly compute the halo bias, however effects such as this
can be captured in the calibration to simulations by the phenomenological
modifications to the basic halo model parameterised by an effective slope at
the non-linear scale.

In our main analysis, we chose to use only the shear correlation functions in
the DES-Y1 data, $\xi_{\pm}$. This is because there are possible unaccounted
for aspects of the galaxy bias in FDM in our model, which might call an
analysis of the galaxy power spectrum, and galaxy-lensing cross correlation
into question. However, we also presented an analysis including these data
with a simple linear bias model for FDM, leading to a slight improvement in
the lower bound on the particle mass to
$\log_{10}(m/\text{eV})\geq-22.8\quad\text{(95\% C.L.)}$ (see the
Supplementary Material).

We also considered varying the statistical aspects of our analysis. We first
compared different sampling techniques of MCMC and nested sampling, finding
agreement between both. We next considered the effect of different priors on
the FDM mass. Our main analysis adopts a log-flat prior between rough limits
imposed by previous analyses, and set by the expected sensitivity of the data.
We compared these results to the case with a uniform prior on the inverse
mass, which allows finite prior volume at the CDM limit. Accounting for the
different prior volumes, the results are consistent (as expected).

Having developed such a rigorous analysis, we can be confident that the
presented bound on the FDM mass is statistically and theoretically robust.

###  5.2 Comparison to Other Results

Our model approximates FDM as having vanishing self-interactions, i.e. scalar
potential $V(\phi)\approx m^{2}\phi^{2}/2$. Axions are expected to have
periodic potentials, and so if the FDM is an ultra-light axion, its self
interactions are described by the potential,
$V=m^{2}f_{a}^{2}[1-\cos{(\phi/f_{a})}]$, where $f_{a}$ is the axion decay
constant. Linares Cedeño et al. (2021) presented constraints on ULAs including
such interactions parameterised by a coupling constant
$\lambda=3M_{pl}^{2}/f_{a}^{2}$, and showed that the limit on $m$ is not
affected significantly by the self-interaction strength for parameter values
where the implementation of linear perturbation theory is numerically
reliable.

Using the CMB, Linares Cedeño et al. (2021) find the bound
$\log_{10}(m/\text{eV})>-23.99$ at 95% C.L., which is slightly stronger than
our bound form the CMB alone. However, Linares Cedeño et al. (2021) took a
maximum value of $\log_{10}(m/\text{eV})=-16$ in their log-flat prior. This is
far beyond the sensitivity of the CMB, and so their limit is dominated by the
choice of prior.

Many complementary cosmological probes can also be used to measure the FDM
mass and density fraction relative to CDM including galaxy formation (Schive
et al., 2016; Bozek et al., 2015; Corasaniti et al., 2017), X-ray observations
(Maleki et al., 2020), the Lyman-$\alpha$ effective opacity (Sarkar et al.,
2021), and the galaxy power spectrum multipoles (which can be used in
conjunction with CMB data) (Laguë et al., 2022). Upcoming cosmological surveys
will further improve sensitivity to FDM. Forecasts from line intensity mapping
surveys project that the ULA fraction could be measured at the percent level
(Bauer et al., 2021), and is sensitive to masses on the order of
$10^{-22}\text{ eV}$ and below. There is similar sensitivity for future CMB
experiments, especially so-called “high definition CMB” (Hložek et al., 2017;
Sehgal et al., 2019). Measurements of the cluster pairwise velocity dispersion
through the kinetic Sunyaev-Zel’dovich (kSZ) effect could probe ULA mass
fractions of $\sim 5\%$ in a window near $m\sim 10^{-27}~{}{\rm eV}$, while
the high-$\ell$ Ostriker-Vishniac anisotropies could reach mass fractions of
$0.1\%$ near $m\sim 10^{-27}~{}{\rm eV}$, and is sensitive to pure FDM up to
masses as high as $m\sim 10^{-22}~{}{\rm eV}$ (Farren et al., 2022). Future
pulsar timing array measurements could probe axion masses around $10^{-22}$ eV
(PPTA Collaboration, 2018), and percent level fractions for $m\approx
10^{-23}\text{ eV}$ (Khmelnitsky & Rubakov, 2014).

Additional bounds on the FDM mass can be obtained using astrophysical data on
the local distribution of DM. Rotation curves and velocity dispersions of many
types of galaxies are affected by the “solitonic core” formed in FDM density
profiles (Schive et al., 2014a; Marsh & Pop, 2015; Chen et al., 2017;
González-Morales et al., 2017; Wasserman et al., 2019; Hayashi et al., 2021;
Bar et al., 2018). We note, however, that these probes are affected by
important systematic errors. Indeed it has been shown that stars, gas and
black holes have significant impacts on the small scale ULA density profile
(Veltmaat et al., 2020; Davies & Mocz, 2020; Chan et al., 2018). Stringent
constraints have been obtained using black hole superradiance (Arvanitaki et
al., 2010; Arvanitaki & Dubovsky, 2011; Stott & Marsh, 2018; Stott, 2020)
which excludes narrow bands in $m$ above the region where FDM has significant
effects on cosmology and galaxy formation (the black hole in M87 reaches
furthest down into the FDM range, Davoudiasl & Denton, 2019) . Heating of
stars by small scale fluctuations of FDM also provides a probe of the particle
mass (Hui et al., 2017). Heating of the Milky Way disk leads to the limit
$m\gtrsim 10^{-22}\text{ eV}$ (Church et al., 2019), while the survival of the
Eridanus-II star cluster excludes a range of masses near $10^{-20}\text{ eV}$
(Marsh & Niemeyer, 2019; Chiang et al., 2021). The Milky Way satellite
population, including Eridanus-II (Marsh & Niemeyer, 2019), and the satellites
surveyed by DES (Nadler et al., 2019), appears to disfavour $m\lesssim
10^{-21}\text{ eV}$, although it has been argued that this bound can be
avoided if the true masses are underestimated, and ultra faint galaxies are in
fact tidally stripped and isolated solitons (Chiang et al., 2021; Broadhurst
et al., 2020).

Finally, we note the recent results of Blum & Teodori (2021) and Allali et al.
(2021), who find that the introduction of a partial ULA dark matter component
may help alleviate the various “tensions” in recent cosmological data
combinations. Blum & Teodori (2021) finds that $10^{-25}$ eV ULAs composing
10% of the dark matter could alleviate the $H_{0}$ tension from gravitational
lensing. Allali et al. (2021) proposes a mixed DM model with ULAs paired with
a decaying dark energy component that can simultaneously resolve the $H_{0}$
and $\sigma_{8}$ tensions.

###  5.3 Looking Forward

On small angular scales, $\theta\lesssim\mathcal{O}(1)$ arcmin, baryonic
effects in the form of Active Galactic Nucleus (AGN) feedback can affect weak
lensing correlation functions. HMCode has been calibrated to simulations of
AGN feedback for a wide range of cosmologies, and in the weak lensing analysis
analysis of KiDS in Joudaki et al. (2020) the uncertainty in the AGN feedback
model is marginalised over. AGN feedback is not expected to be precisely the
same for CDM and FDM. The DES-Y1 analysis masks small angular scales in order
to avoid uncertainty due to AGN feedback (indicated by the grey regions in
Figs. 2, 12 and 13), and so our present analysis is not affected by the model
of AGN feedback. Taking full advantage of weak lensing on small angular scales
to constrain the nature of DM requires simulations including AGN feedback
beyond CDM in order to calibrate e.g. HMCode and emulators.

If such calibrations were performed, how could weak lensing tests of FDM
improve in the future? We can estimate this using the bounds to the effective
linear power spectrum (Tegmark & Zaldarriaga, 2002; Chabanier et al., 2019).
Tegmark & Zaldarriaga (2002) showed that the effective linear theory
wavenumber $k\propto 1/\theta$ for weak lensing correlation functions measured
on an angular scale $\theta$. In our analysis, the smallest value of $\theta$
was $\theta\approx 4$ arcmin in $\xi_{+}$. The smallest scales used in the
DES-Y3 analysis (Abbott et al., 2022) are around 1 arcmin. Using again that
$k_{\rm J,eq}\propto m^{1/2}$ (Eq. 1), substantial improvement in ULA
constraints is possible with existing data, to $m\sim 10^{-22}\text{ eV}$ (at
the time of writing the DES-Y3 likelihood is not yet available). The _Euclid_
angular resolution is around 0.1 arcmin (Laureijs et al., 2011), which by the
same logic might be sensitive to $m\sim 10^{-20}\text{ eV}$, comparable to the
Lyman-alpha forest constraints of Rogers & Peiris (2021b).

When smaller scales are included in weak lensing analyses of FDM, the effects
on the power spectrum from the models of $c(M)$ and $n(M)$ will become
increasingly important. A rigorous analysis of such data will require the
tools we have developed to account for systematic uncertainty in the low mass
behaviour of $c(M)$ and $n(M)$. We also advocate more in depth studies of the
one-to-two halo transition region in simulations of FDM and mixed models of
CDM and FDM.

## Acknowledgments

MD and DJEM were supported at the University of Göttingen by the Alexander von
Humboldt Foundationa and the German Federal Ministry of Education and
Research. DJEM is supported at King’s College London by an Ernest Rutherford
Fellowship of the Science and Technologies Facilities Council (UK). RH is a
CIFAR Azrieli Global Scholar in the Gravity and the Extreme Universe Program,
and acknowledges funding from the NSERC Discovery Grants program, the Alfred
P. Sloan Foundation and the Connaught Fund. The Dunlap Institute is funded
through an endowment established by the David Dunlap family and the University
of Toronto. The authors at the University of Toronto 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. They are grateful to have the opportunity to work in
the community, on this territory. DG acknowledges support in part by NASA ATP
grant 17-ATP17-0162 and the Provost’s office at Haverford College. We would
like to thank the following colleagues for useful discussions: Sebastian Hoof,
Shahab Joudaki, Alex van Engelen. We thank Aaron Ludlow for supplying code to
reproduce the work of Ludlow et al. (2016). MD thanks members of the Munich
cosmology group for helpful suggestions. We analysed some of our results using
ChainConsumer (Hinton, 2016), and made use of the following open source
libraries: numpy Harris et al. (2020), matplotlib Hunter (2007). Our cosmosis
analysis relies, in addition to what is cited in the text, on the following
works: Kirk et al. (2012); Bridle & King (2007); Kilbinger et al. (2009). The
_Planck_ data used in this article was accessed through the likelihoods in
cosmosis, the original data are available here: Planck Legacy Archive.

This project used public archival data from the Dark Energy Survey (DES),
based on Abbott et al. (2018b); Morganson et al. (2018); Flaugher et al.
(2015), Data Products. Funding for the DES Projects has been provided by the
U.S. Department of Energy, the U.S. National Science Foundation, the Ministry
of Science and Education of Spain, the Science and Technology
FacilitiesCouncil of the United Kingdom, the Higher Education Funding Council
for England, the National Center for Supercomputing Applications at the
University of Illinois at Urbana-Champaign, the Kavli Institute of
Cosmological Physics at the University of Chicago, the Center for Cosmology
and Astro-Particle Physics at the Ohio State University, the Mitchell
Institute for Fundamental Physics and Astronomy at Texas A&M University,
Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à
Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento
Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação,
the Deutsche Forschungsgemeinschaft, and the Collaborating Institutions in the
Dark Energy Survey. The Collaborating Institutions are Argonne National
Laboratory, the University of California at Santa Cruz, the University of
Cambridge, Centro de Investigaciones Energéticas, Medioambientales y
Tecnológicas-Madrid, the University of Chicago, University College London, the
DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische
Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the
University of Illinois at Urbana-Champaign, the Institut de Ciències de
l’Espai (IEEC/CSIC), the Institut de Física d’Altes Energies, Lawrence
Berkeley National Laboratory, the Ludwig-Maximilians Universität München and
the associated Excellence Cluster Universe, the University of Michigan, the
National Optical Astronomy Observatory, the University of Nottingham, The Ohio
State University, the OzDES Membership Consortium, the University of
Pennsylvania, the University of Portsmouth, SLAC National Accelerator
Laboratory, Stanford University, the University of Sussex, and Texas A&M
University. Based in part on observations at Cerro Tololo Inter-American
Observatory, National Optical Astronomy Observatory, which is operated by the
Association of Universities for Research in Astronomy (AURA) under a
cooperative agreement with the National Science Foundation.

## Data Availability Statement

The data underlying this article, in the form of chains, will be shared on
reasonable request to the corresponding author.

## References

  * Abbott & Sikivie (1983) Abbott L. F., Sikivie P., 1983, Phys. Lett. B, 120, 133
  * Abbott et al. (2018a) Abbott T. M. C., et al., 2018a, Phys. Rev. D, 98, 043526
  * Abbott et al. (2018b) Abbott T. M. C., et al., 2018b, ApJS, 239, 18
  * Abbott et al. (2022) Abbott T. M. C., et al., 2022, Phys. Rev. D, 105, 023520
  * Aghanim et al. (2020a) Aghanim N., et al., 2020a, Astron. Astrophys., 641, A5
  * Aghanim et al. (2020b) Aghanim N., et al., 2020b, Astron. Astrophys., 641, A6
  * Allali et al. (2021) Allali I. J., Hertzberg M. P., Rompineve F., 2021, arXiv e-prints, p. arXiv:2104.12798
  * Amendola & Barbieri (2006) Amendola L., Barbieri R., 2006, Phys. Lett. B, 642, 192
  * Amendola et al. (2018) Amendola L., et al., 2018, Living Rev. Rel., 21, 2
  * Anderson et al. (2019) Anderson L., Pontzen A., Font-Ribera A., Villaescusa-Navarro F., Rogers K. K., Genel S., 2019, Astrophys. J., 871, 144
  * Armengaud et al. (2017) Armengaud E., Palanque-Delabrouille N., Yèche C., Marsh D. J. E., Baur J., 2017, Mon. Not. Roy. Astron. Soc., 471, 4606
  * Arvanitaki & Dubovsky (2011) Arvanitaki A., Dubovsky S., 2011, Phys. Rev. D, 83, 044026
  * Arvanitaki et al. (2010) Arvanitaki A., Dimopoulos S., Dubovsky S., Kaloper N., March-Russell J., 2010, Phys. Rev. D, 81, 123530
  * Babichev et al. (2016) Babichev E., Marzola L., Raidal M., Schmidt-May A., Urban F., Veermäe H., von Strauss M., 2016, JCAP, 09, 016
  * Bar et al. (2018) Bar N., Blas D., Blum K., Sibiryakov S., 2018, Phys. Rev. D, 98, 083027
  * Bauer et al. (2021) Bauer J. B., Marsh D. J. E., Hložek R., Padmanabhan H., Laguë A., 2021, MNRAS, 500, 3162
  * Benson et al. (2013) Benson A. J., et al., 2013, MNRAS, 428, 1774
  * Bird et al. (2012) Bird S., Viel M., Haehnelt M. G., 2012, MNRAS, 420, 2551
  * Bird et al. (2019) Bird S., Rogers K. K., Peiris H. V., Verde L., Font-Ribera A., Pontzen A., 2019, JCAP, 02, 050
  * Blum & Teodori (2021) Blum K., Teodori L., 2021, arXiv e-prints, p. arXiv:2105.10873
  * Boera et al. (2019) Boera E., Becker G. D., Bolton J. S., Nasir F., 2019, ApJ, 872, 101
  * Bond & Myers (1996) Bond J. R., Myers S. T., 1996, ApJS, 103, 1
  * Bozek et al. (2015) Bozek B., Marsh D. J. E., Silk J., Wyse R. F. G., 2015, Mon. Not. Roy. Astron. Soc., 450, 209
  * Bridle & King (2007) Bridle S., King L., 2007, New J. Phys., 9, 444
  * Broadhurst et al. (2020) Broadhurst T., de Martino I., Luu H. N., Smoot G. F., Tye S. H. H., 2020, Phys. Rev. D, 101, 083012
  * Bullock et al. (2001) Bullock J. S., Kolatt T. S., Sigad Y., Somerville R. S., Kravtsov A. V., Klypin A. A., Primack J. R., Dekel A., 2001, MNRAS, 321, 559
  * Carucci et al. (2015) Carucci I. P., Villaescusa-Navarro F., Viel M., Lapi A., 2015, JCAP, 07, 047
  * Chabanier et al. (2019) Chabanier S., Millea M., Palanque-Delabrouille N., 2019, Mon. Not. Roy. Astron. Soc., 489, 2247
  * Chan et al. (2018) Chan J. H. H., Schive H.-Y., Woo T.-P., Chiueh T., 2018, MNRAS, 478, 2686
  * Chen et al. (2017) Chen S.-R., Schive H.-Y., Chiueh T., 2017, MNRAS, 468, 1338
  * Chen et al. (2021) Chen J., Du X., Lentz E. W., Marsh D. J. E., Niemeyer J. C., 2021, Phys. Rev. D, 104, 083022
  * Chiang et al. (2021) Chiang B. T., Schive H.-Y., Chiueh T., 2021, Phys. Rev. D, 103, 103019
  * Church et al. (2019) Church B. V., Mocz P., Ostriker J. P., 2019, MNRAS, 485, 2861
  * Cicoli et al. (2022) Cicoli M., Guidetti V., Righi N., Westphal A., 2022, JHEP, 05, 107
  * Cooray & Sheth (2002) Cooray A., Sheth R., 2002, Physics reports, 372, 1
  * Corasaniti et al. (2017) Corasaniti P. S., Agarwal S., Marsh D. J. E., Das S., 2017, Phys. Rev. D, 95, 083512
  * Davies & Mocz (2020) Davies E. Y., Mocz P., 2020, MNRAS, 492, 5721
  * Davoudiasl & Denton (2019) Davoudiasl H., Denton P. B., 2019, Phys. Rev. Lett., 123, 021102
  * Davoudiasl & Murphy (2017) Davoudiasl H., Murphy C. W., 2017, Phys. Rev. Lett., 118, 141801
  * Dine & Fischler (1983) Dine M., Fischler W., 1983, Phys. Lett. B, 120, 137
  * Du (2018) Du X., 2018, PhD Thesis. Georg-August Universität, Göttingen
  * Du et al. (2017) Du X., Behrens C., Niemeyer J. C., Schwabe B., 2017, Phys. Rev. D, 95, 043519
  * Dunkley et al. (2005) Dunkley J., Bucher M., Ferreira P. G., Moodley K., Skordis C., 2005, Mon. Not. Roy. Astron. Soc., 356, 925
  * Farren et al. (2022) Farren G. S., Grin D., Jaffe A. H., Hložek R., Marsh D. J. E., 2022, Phys. Rev. D, 105, 063513
  * Feroz et al. (2009) Feroz F., Hobson M. P., Bridges M., 2009, Mon. Not. Roy. Astron. Soc., 398, 1601
  * Flaugher et al. (2015) Flaugher B., et al., 2015, AJ, 150, 150
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, Publ. Astron. Soc. Pac., 125, 306
  * González-Morales et al. (2017) González-Morales A. X., Marsh D. J. E., Peñarrubia J., Ureña López L. A., 2017, Mon. Not. Roy. Astron. Soc., 472, 1346
  * Goodman & Weare (2010) Goodman J., Weare J., 2010, Communications in Applied Mathematics and Computational Science, 5, 65
  * Graham et al. (2016) Graham P. W., Mardon J., Rajendran S., 2016, Phys. Rev. D, 93, 103520
  * Green & Kavanagh (2021) Green A. M., Kavanagh B. J., 2021, J. Phys. G, 48, 4
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hayashi & White (2008a) Hayashi E., White S. D., 2008a, Monthly Notices of the Royal Astronomical Society, 388, 2
  * Hayashi & White (2008b) Hayashi E., White S. D. M., 2008b, Monthly Notices of the Royal Astronomical Society, 388, 2
  * Hayashi et al. (2021) Hayashi K., Ferreira E. G. M., Chan H. Y. J., 2021, ApJ, 912, L3
  * Helfer et al. (2017) Helfer T., Marsh D. J. E., Clough K., Fairbairn M., Lim E. A., Becerril R., 2017, JCAP, 03, 055
  * Hinton (2016) Hinton S. R., 2016, The Journal of Open Source Software, 1, 00045
  * Hložek et al. (2017) Hložek R., Marsh D. J. E., Grin D., Allison R., Dunkley J., Calabrese E., 2017, Phys. Rev. D, 95, 123511
  * Hlozek et al. (2014) Hlozek R., Grin D., Marsh D. J. E., Ferreira P. G., 2014, https://github.com/dgrin1/axionCAMB
  * Hlozek et al. (2015) Hlozek R., Grin D., Marsh D. J. E., Ferreira P. G., 2015, Phys. Rev. D, 91, 103512
  * Hlozek et al. (2018) Hlozek R., Marsh D. J. E., Grin D., 2018, Mon. Not. Roy. Astron. Soc., 476, 3063
  * Hu et al. (2000) Hu W., Barkana R., Gruzinov A., 2000, Phys. Rev. Lett., 85, 1158
  * Hui et al. (2017) Hui L., Ostriker J. P., Tremaine S., Witten E., 2017, Phys. Rev. D, 95, 043541
  * Hunter (2007) Hunter J. D., 2007, Computing in Science Engineering, 9, 90
  * Iršič et al. (2017a) Iršič V., et al., 2017a, Phys. Rev. D, 96, 023522
  * Iršič et al. (2017b) Iršič V., Viel M., Haehnelt M. G., Bolton J. S., Becker G. D., 2017b, Phys. Rev. Lett., 119, 031302
  * Joudaki et al. (2020) Joudaki S., et al., 2020, Astron. Astrophys., 638, L1
  * Kaup (1968) Kaup D. J., 1968, Phys. Rev., 172, 1331
  * Khlopov et al. (1985) Khlopov M., Malomed B. A., Zeldovich I. B., 1985, Mon. Not. Roy. Astron. Soc., 215, 575
  * Khmelnitsky & Rubakov (2014) Khmelnitsky A., Rubakov V., 2014, JCAP, 02, 019
  * Kilbinger et al. (2009) Kilbinger M., et al., 2009, Astron. Astrophys., 497, 677
  * Kim & Marsh (2016) Kim J. E., Marsh D. J. E., 2016, Phys. Rev. D, 93, 025027
  * Kirk et al. (2012) Kirk D., Rassat A., Host O., Bridle S., 2012, MNRAS, 424, 1647
  * Kobayashi et al. (2017) Kobayashi T., Murgia R., De Simone A., Iršič V., Viel M., 2017, Phys. Rev. D, 96, 123514
  * Krause et al. (2017) Krause E., et al., 2017
  * Laguë et al. (2022) Laguë A., Bond J. R., Hložek R., Rogers K. K., Marsh D. J. E., Grin D., 2022, JCAP, 01, 049
  * Laureijs et al. (2011) Laureijs R., et al., 2011, arXiv e-prints, p. arXiv:1110.3193
  * Lawrence et al. (2010a) Lawrence E., Heitmann K., White M., Higdon D., Wagner C., Habib S., Williams B., 2010a, ApJ, 713, 1322
  * Lawrence et al. (2010b) Lawrence E., Heitmann K., White M., Higdon D., Wagner C., Habib S., Williams B., 2010b, ApJ, 713, 1322
  * Levkov et al. (2017) Levkov D. G., Panin A. G., Tkachev I. I., 2017, Phys. Rev. Lett., 118, 011301
  * Levkov et al. (2018) Levkov D. G., Panin A. G., Tkachev I. I., 2018, Phys. Rev. Lett., 121, 151301
  * Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, Astrophys. J., 538, 473
  * Li et al. (2014) Li B., Rindler-Daller T., Shapiro P. R., 2014, Phys. Rev. D, 89, 083536
  * Li et al. (2019) Li X., Hui L., Bryan G. L., 2019, Phys. Rev. D, 99, 063509
  * Linares Cedeño et al. (2021) Linares Cedeño F. X., González-Morales A. X., Ureña López L. A., 2021, JCAP, 01, 051
  * Ludlow et al. (2016) Ludlow A. D., Bose S., Angulo R. E., Wang L., Hellwing W. A., Navarro J. F., Cole S., Frenk C. S., 2016, Mon. Not. Roy. Astron. Soc., 460, 1214
  * Lukić et al. (2015) Lukić Z., Stark C. W., Nugent P., White M., Meiksin A. A., Almgren A., 2015, MNRAS, 446, 3697
  * Maleki et al. (2020) Maleki A., Baghram S., Rahvar S., 2020, Phys. Rev. D, 101, 023508
  * Marsh (2016a) Marsh D. J. E., 2016a, arXiv:1605.05973
  * Marsh (2016b) Marsh D. J. E., 2016b, Phys. Rept., 643, 1
  * Marsh & Ferreira (2010) Marsh D. J. E., Ferreira P. G., 2010, Phys. Rev. D, 82, 103528
  * Marsh & Hoof (2021) Marsh D. J. E., Hoof S., 2021, arXiv:2106.08797
  * Marsh & Niemeyer (2019) Marsh D. J. E., Niemeyer J. C., 2019, Phys. Rev. Lett., 123, 051103
  * Marsh & Pop (2015) Marsh D. J. E., Pop A.-R., 2015, Mon. Not. Roy. Astron. Soc., 451, 2479
  * Marsh & Silk (2014) Marsh D. J. E., Silk J., 2014, Mon. Not. Roy. Astron. Soc., 437, 2652
  * Marsh et al. (2012) Marsh D. J. E., Macaulay E., Trebitsch M., Ferreira P. G., 2012, Phys. Rev. D, 85, 103514
  * Massara et al. (2014) Massara E., Villaescusa-Navarro F., Viel M., 2014, JCAP, 12, 053
  * McDonald et al. (2006) McDonald P., et al., 2006, Astrophys. J. Suppl., 163, 80
  * Mead et al. (2015) Mead A., Peacock J., Heymans C., Joudaki S., Heavens A., 2015, Mon. Not. Roy. Astron. Soc., 454, 1958
  * Mead et al. (2021) Mead A. J., Brieden S., Tröster T., Heymans C., 2021, MNRAS, 502, 1401
  * Mehta et al. (2021) Mehta V. M., Demirtas M., Long C., Marsh D. J. E., McAllister L., Stott M. J., 2021, JCAP, 07, 033
  * Michel & Moss (2018) Michel F., Moss I. G., 2018, Phys. Lett. B, 785, 9
  * Mina et al. (2022) Mina M., Mota D. F., Winther H. A., 2022, Astron. Astrophys., 662, A29
  * Mo et al. (2010) Mo H., Van den Bosch F., White S., 2010, Galaxy formation and evolution. Cambridge University Press
  * Mocz et al. (2019) Mocz P., et al., 2019, Phys. Rev. Lett., 123, 141301
  * Molaro et al. (2021) Molaro M., et al., 2021, Mon. Not. Roy. Astron. Soc., 509, 6119
  * Morganson et al. (2018) Morganson E., et al., 2018, PASP, 130, 074501
  * Nadler et al. (2019) Nadler E. O., Gluscevic V., Boddy K. K., Wechsler R. H., 2019, ApJ, 878, L32
  * Narayanan et al. (2000) Narayanan V. K., Spergel D. N., Dave R., Ma C.-P., 2000, Astrophys. J. Lett., 543, L103
  * Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D., 1997, The Astrophysical Journal, 490, 493
  * Nori et al. (2019) Nori M., Murgia R., Iršič V., Baldi M., Viel M., 2019, Mon. Not. Roy. Astron. Soc., 482, 3227
  * PPTA Collaboration (2018) PPTA Collaboration 2018, Phys. Rev. D, 98, 102002
  * Pandey et al. (2021) Pandey S., et al., 2021
  * Peccei & Quinn (1977) Peccei R. D., Quinn H. R., 1977, Phys. Rev. Lett., 38, 1440
  * Porredon et al. (2021) Porredon A., et al., 2021
  * Poulin et al. (2018) Poulin V., Smith T. L., Grin D., Karwal T., Kamionkowski M., 2018, Phys. Rev. D, 98, 083525
  * Preskill et al. (1983) Preskill J., Wise M. B., Wilczek F., 1983, Phys. Lett. B, 120, 127
  * Press & Schechter (1978) Press W., Schechter P., 1978, ApJ, 187, 12
  * Robles et al. (2018) Robles V. H., Bullock J. S., Boylan-Kolchin M., 2018, Monthly Notices of the Royal Astronomical Society, 483, 289
  * Rogers & Peiris (2021a) Rogers K. K., Peiris H. V., 2021a, Phys. Rev. D, 103, 043526
  * Rogers & Peiris (2021b) Rogers K. K., Peiris H. V., 2021b, Phys. Rev. Lett., 126, 071302
  * Rogers et al. (2019) Rogers K. K., Peiris H. V., Pontzen A., Bird S., Verde L., Font-Ribera A., 2019, J. Cosmology Astropart. Phys., 2019, 031
  * Ruffini & Bonazzola (1969) Ruffini R., Bonazzola S., 1969, Phys. Rev., 187, 1767
  * Salehian et al. (2021) Salehian B., Zhang H.-Y., Amin M. A., Kaiser D. I., Namjoo M. H., 2021, JHEP, 09, 050
  * Sarkar et al. (2016) Sarkar A., Mondal R., Das S., Sethi S. K., Bharadwaj S., Marsh D. J. E., 2016, JCAP, 04, 012
  * Sarkar et al. (2021) Sarkar Anjan K., Pandey Kanhaiya L., Sethi K S., 2021, arXiv e-prints, p. arXiv:2101.09917
  * Schive et al. (2014a) Schive H.-Y., Chiueh T., Broadhurst T., 2014a, Nature Physics, 10, 496
  * Schive et al. (2014b) Schive H.-Y., Liao M.-H., Woo T.-P., Wong S.-K., Chiueh T., Broadhurst T., Hwang W. Y. P., 2014b, Phys. Rev. Lett., 113, 261302
  * Schive et al. (2016) Schive H.-Y., Chiueh T., Broadhurst T., Huang K.-W., 2016, Astrophys. J., 818, 89
  * Schneider et al. (2012) Schneider A., Smith R. E., Macciò A. V., Moore B., 2012, MNRAS, 424, 684
  * Schneider et al. (2013) Schneider A., Smith R. E., Reed D., 2013, Mon. Not. Roy. Astron. Soc., 433, 1573
  * Schwabe & Niemeyer (2021) Schwabe B., Niemeyer J. C., 2021, arXiv preprint arXiv:2110.09145
  * Sehgal et al. (2019) Sehgal N., et al., 2019, 51, 6
  * Seidel & Suen (1991) Seidel E., Suen W. M., 1991, Phys. Rev. Lett., 66, 1659
  * Seljak et al. (2006) Seljak U., Slosar A., McDonald P., 2006, JCAP, 10, 014
  * Sheth & Tormen (1999) Sheth R. K., Tormen G., 1999, Mon. Not. Roy. Astron. Soc., 308, 119
  * Sheth et al. (2001) Sheth R. K., Mo H. J., Tormen G., 2001, Mon. Not. Roy. Astron. Soc., 323, 1
  * Smith & Markovic (2011) Smith R. E., Markovic K., 2011, Phys. Rev. D, 84, 063507
  * Smith et al. (2003) Smith R. E., et al., 2003, MNRAS, 341, 1311
  * Stott (2020) Stott M. J., 2020, arXiv e-prints, p. arXiv:2009.07206
  * Stott & Marsh (2018) Stott M. J., Marsh D. J. E., 2018, Phys. Rev. D, 98, 083006
  * Sullivan et al. (2021) Sullivan J. M., Seljak U., Singh S., 2021, JCAP, 11, 026
  * Takahashi et al. (2012) Takahashi R., Sato M., Nishimichi T., Taruya A., Oguri M., 2012, Astrophys. J., 761, 152
  * Tegmark & Zaldarriaga (2002) Tegmark M., Zaldarriaga M., 2002, Phys. Rev. D, 66, 103508
  * Tegmark et al. (2006) Tegmark M., Aguirre A., Rees M., Wilczek F., 2006, Phys. Rev. D, 73, 023505
  * Turner (1983) Turner M. S., 1983, Phys. Rev. D, 28, 1243
  * Uhlemann et al. (2014) Uhlemann C., Kopp M., Haugg T., 2014, Phys. Rev. D, 90, 023517
  * Veltmaat & Niemeyer (2016) Veltmaat J., Niemeyer J. C., 2016, Phys. Rev. D, 94, 123523
  * Veltmaat et al. (2020) Veltmaat J., Schwabe B., Niemeyer J. C., 2020, Phys. Rev. D, 101, 083518
  * Viel et al. (2004) Viel M., Haehnelt M. G., Springel V., 2004, Monthly Notices of the Royal Astronomical Society, 354, 684
  * Viel et al. (2013) Viel M., Becker G. D., Bolton J. S., Haehnelt M. G., 2013, Phys. Rev. D, 88, 043502
  * Wang & White (2007) Wang J., White S. D. M., 2007, Mon. Not. Roy. Astron. Soc., 380, 93
  * Wasserman et al. (2019) Wasserman A., et al., 2019, ApJ, 885, 155
  * Weinberg (1978) Weinberg S., 1978, Phys. Rev. Lett., 40, 223
  * Weinberg (1987) Weinberg S., 1987, Phys. Rev. Lett., 59, 2607
  * Widrow & Kaiser (1993) Widrow L. M., Kaiser N., 1993, ApJ, 416, L71
  * Wilczek (1978) Wilczek F., 1978, Phys. Rev. Lett., 40, 279
  * Wu et al. (2019) Wu X., McQuinn M., Kannan R., D’Aloisio A., Bird S., Marinacci F., Davé R., Hernquist L., 2019, MNRAS, 490, 3177
  * Zhang et al. (2018) Zhang J., Kuo J.-L., Liu H., Sming Tsai Y.-L., Cheung K., Chu M.-C., 2018, ApJ, 863, 73
  * Zuntz et al. (2015) Zuntz J., et al., 2015, Astronomy and Computing, 12, 45
  * de Salas & Widmark (2021) de Salas P. F., Widmark A., 2021, Rept. Prog. Phys., 84, 104901

##  Appendix A DES-Y1 3x2pt

Figure 12: Cosmic shear correlation function between all combinations $i\,j$
of source bins, as indicated in the top left corner of each panel. We display
the $\xi_{+}$ observable on the left hand side and the $\xi_{-}$ observable on
the right hand side (note the different order of bin combinations shown on
each side). We show DES-Y1 data points in black, marking the square-roots of
the diagonal entries of the covariance matrix as vertical bars. We compare the
CDM predictions (black, dashed lines) to the ULA predictions taking four
different values for the axion mass parameter $m_{Ax}$ within our prior range
(solid, colour-coded). The grey shaded regions mask angular scales excluded
from the analysis due to modelling uncertainties. The case $m=10^{-22}\text{
eV}$ is indistinguishable from CDM over the scales shown.  Figure 13: The
$\gamma$ and $w$ correlation functions that form the “3x2pt” observable when
combined with the $\xi_{\pm}.$ The results obtained when using the full data
set are compared to our baseline analysis in Fig. 14, and show a modest
improvement in the bound on the axion mass. Figure 14: Posterior density for
the DM particle mass $m$, the mass fluctuation parameter $S_{8}$, and the
matter density $\Omega_{m}$ for different combinations of data sets. DES
shear-only data set (pink) and the complete 3x2pt data set (red). Note that we
changed the range of the axis compared to Fig. 8 for better visibility.

The DES-Y1 data also contains the galaxy power spectrum in the form of the
correlation function, $w$, and the cross correlation between shear and galaxy
number, $\gamma_{t}$, leading to the complete “3x2pt” dataset. These data
(shown in Figure 13, along with the predictions for CDM and FDM given
reference cosmological and bias parameters) are not included in our baseline
analysis, since $w$ and $\gamma_{t}$ depend on the model of galaxy bias, for
which a linear, scale independent model may not suffice for FDM (as discussed