Fuzzy Dark Matter and the Dark Energy Survey Year 1 Data

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 \textsc{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 \emph{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 \emph{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.


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 . Observations constrain the parameter space of PBHs tightly, leaving a single window around M ≈ 10 −12 M 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 Corresponding Author: david.j.marsh@kcl.ac.uk self interactions at early times (Marsh & Hoof 2021). We consider DM composed of a scalar field, φ , with potential V (φ ) = m 2 φ 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 eV m 10 −19 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 ≈ 10 17 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", 1 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 O(few) × 10 −25 eV for the dominant component of DM (Hlozek et al. 2015(Hlozek et al. , 2018Poulin 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 nonlinear scales is required.
One such rigorous approach has been developed for the Lyman-α 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 highestresolution Lyman-α forest probes the matter power spectrum as traced by the intergalactic medium (IGM) for line-of-sight velocity wavenumbers k f < 0.2 s 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 spaces 2 (e.g. Lawrence et al. 2010a). In particular, Bird et al. (2019); Rogers et al. (2019); Rogers & Peiris (2021b,a) built Lyman-α 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 × 10 −20 eV at 95% credibility. The methodology to marginalise uncertainties improves on previous Lyman-α 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-α 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-α 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 dynam- Figure 1. One dimensional posterior on the FDM particle mass, m, assuming a log-flat prior. The DES cosmic shear correlation function, ξ ± 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 ξ ± induced by the FDM Jeans scale, k J,eq ∝ 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.
ics indicate that, for current data, the effect is only on the order of a few percent (Nori et al. 2019). Further, existing Lyman-α forest bounds neglect the large-scale fluctuations (∼ 40 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 (∼ Gpc) and the high-resolution Lyman-α forest (∼ Mpc).
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. 3 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 HM-CODE (Mead et al. 2015) provides such a model, which can be adapted to both warm DM (WDM) and FDM (Marsh 2016a). 4 HM-CODE 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 z 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 concentrationmass 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 Mpc −1 , while the final DES bin covers 1 h Mpc −1 < k < 5 h Mpc −1 . The ULA power spectrum is suppressed below the Jeans scale at matter-radiation equality, given by (Hu et al. 2000): 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.

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 z 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 n i = dz 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(θ ), galaxy-galaxy lensing γ t (θ ), and the two components of the cosmic shear ξ ± (θ ). In detail, w(θ ) measures the distribution of pairs of angular separations of galaxies as compared to what is expected assuming a random distribution. The observables γ t (θ ) and ξ ± (θ ) quantify shape distortions due to lensing by foreground mass distributions. γ t (θ ) measures, for pairs of source and lens galaxies, the ellipticities of the the source galaxies tangential to the line connecting to the lens galaxy. ξ + (θ ) and ξ − (θ ), 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( , χ) along this direction (Abbott et al. 2018a;Krause et al. 2017): where the weighting functions q i g/κ are defined as and The comoving distance is χ, and b is the linear bias. The final observables are calculated by transforming to real-space as follows: where J ν 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(χ)) 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 i j κg ( ) and hence the observable γ t depends linearly on the galaxy bias, while the projection P i j gg ( ) and hence the observable w depends quadratically on the galaxy bias. Thus, measuring both w and γ t allows the linear bias parameters to be measured (Abbott et al. 2018a).
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 ≤ 10 −24 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 nonlinear 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 10 −24 eV consistent with CMB and BOSS lower bounds) the galaxy bias for ULAs will be linear and (almost) scaleindependent, 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 ξ ± ), 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. A1. The full 3x2pt analysis is considered in the Supplementary Material, where we show how the respective observables γ t and w are affected by FDM (Fig. A2). As anticipated in the introduction, for masses 0.5 × 10 −23 eV, the predictions from FDM model differ visibly from the CDM case.
For the cosmic shear observables ξ ± , the angular scales θ measured by DES-Y1 and shown in Fig. A1 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 5 hMpc −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.

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 10 −27 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(Hlozek et al. , 2014. 5 AX-IONCAMB 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 nonlinear 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.

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 HM-CODE (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.

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 ξ (x 1 , x 2 ), the Fourier transform of the power spectrum P(k) Assuming that all matter is contained with halos of density ρ(x), we can factorise ξ (x 1 , x 2 ) as (Cooray & Sheth 2002;Hayashi & White 2008a) ξ  where we define the two-halo term ξ 2h as the correlation function between densities attributed to two different halos h, h , while the one-halo term ξ 1h is defined as the correlation between densities attributed to the same halo. It is commonly assumed that the density distribution within a halo, ρ h , is a function of space with the total mass M of the halo as a single parameter, i.e. ρ h = ρ h (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) 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 ξ hh (r 1 , r 2 |M 1 , M 2 ) encodes the correlation between two halos of mass M 1 at point r 1 and mass M 2 at point r 2 , respectively. If the correlation between halos varies slowly on scales of the order of the halo sizes, the halo density distributions ρ h 1 , ρ h 2 in Eq. (10) can be approximated by Dirac-δ -functions times the respective masses, i.e. ρ h (x − r|M) ∼ M δ D (x − r). The correlation function ξ hh itself can as first approximation be replaced by the correlation function ξ lin obtained from linear theory, times the linear bias b(M) which is justified at large scales |x 1 − x 2 | 1 because at large scales we can again assume ρ h (x − r|M) ∼ M δ D (x − r), and dM n(M) b(M) M =ρ 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 HM-CODE 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 ρ h (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 where the parameters a = 0.707, p = 0.3 and A = 0.2162 are obtained from empirical fits to ΛCDM simulations (Mead et al. 2015). The variable ν is defined by the ratio In the above equation, the variable δ c denotes the critical density defined as the minimum density for which a region collapses under its own gravity. The variable σ 2 , defined by denotes the autocorrelation function measuring the correlation between smoothed overdensities δ s (r, M), which are averaged over a sphere with radius R as follows: where the smoothing radius R = R(M) is chosen such that the corresponding sphere contains on average the mass M = 4π/3ρR 3 (Mo et al. 2010). The conventional model for the second component of the halo model, the halo density profile, ρ h (r|M), is the Navarro-Frenk-White (NFW) profile, defined as (Navarro et al. 1997;Hayashi & White 2008b) ρ h (r|ρ s , r s ) −ρ = ρ s (r/r s )(1 + r/r s ) 2 , which depends on the two parameters scale density ρ s and scale radius r s . In order to obtain a finite halo mass M = ρ 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 and noting we can express the NFW profile in terms of the variables M, R vir and c as follows (Cooray & Sheth 2002) ρ h (r|ρ s , r s ) With the halo density profile and the halo mass function defined, we have constructed a general halo model.

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. 6 For the Halo Model parameters halo mass M, concentration parameter c, and critical density δ c , the authors use the following prescriptions: The halo mass parameter M is defined such that the average density of the haloρ h (M) is greater than the background densitȳ ρ by a factor ∆ v =ρ h /ρ. ∆ 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) where the minimum halo concentration parameter A in the above equation is again fitted to simulations (Mead et al. 2015). The redshift z 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 coll = 0.01 as compared to the current mass at redshift z, estimated according to the relation 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 σ , which we denote explicitly in Eqs. (21)) and (22)). The critical density for collapse, δ 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 δ 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 where f fixed by an empirical fit, and ∆ 2 lin (k) = 4πV (k/2π) 3 P(k), where the power spectrum P(k) is the Fourier transform of the correlation function ξ (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 where ∆ 2 1h,gen = 4πV (k/2π) 3 P 1h,gen (k), with P 1h,gen (k) being the Fourier transform of the generic one-halo term ξ 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: where η 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 ν = 1 does not change, while for halos with ν > 1 (ν < 1) the profile flattens (sharpens), as ρ s decreases (increases) with ν −3η while the scale radius r s and the cut-off increase (decrease) as r s ν η and R vir ν η , 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 with α an empirical fit parameter given in Table 1, and ∆ 2 2h and ∆ 2 1h defined according to Eqs. (23) and (25), respectively. The fit parameter α is determined by the effective spectral index, n eff , given by: where σ 2 is the variance of the linear power spectrum, and the nonlinear scale R nl is defined by σ 2 (R 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 δ c the virialised overdensity ∆ v = 3m/4πR 3 virρ used to define R vir , as well as the parameter A used to define the concentration parameter c(A). See Table 1.

Halo Model and Lensing Observables for FDM
The lensing observables for FDM are shown in Figs. A1 and A2. 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.

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 σ (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 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 ST (P CDM ), introduced in Eq. (13), for a CDM linear power spectrum. In dashed, blue we show the Sheth-Tormen model n ST (P FDM ) for a FDM linear power spectrum. Below some halo mass, the Sheth-Tormen model predicts a suppression of the n ST (P FDM ) with respect to n ST (P 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 ST (P 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 δ c as implemented in HMCODE , we relate the the halo mass M to a virial radius R 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 where n(P CDM ) is the HMF taking into account the additional suppression observed in simulations and n ST (P CDM ) is the original HMF for CDM in the Sheth-Tormen model. The correction function with respect to the CDM case n ST (P CDM ) is defined as with a scaling mass M To implement the additional suppression into HMCODE , we want   Our model of the additional suppression is shown in green. To account for modeling uncertainties we introduce two nuissance parameters α 1 , α 2 (cf. Eq. 31). We show the effect of varying α 1 (α 2 ), keeping α 2 (α 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.
to correct the Sheth-Tormen model for FDM, n ST (P FDM ). We use the ansatz where we choose the correction function to be of the same functional form as Eq. (30), and n ST (P FDM ) is the Sheth-Tormen mass function computed with the FDM linear power spectrum. This is justified because for small halo masses M M (n) 0 , the Sheth-Tormen prediction for both, the CDM and the ULA case are well described by a powerlaw. Re-fitting the slope-parameter α 2 within our fiducial model, we find α 2 = 1.86 for our correction function in Eq. (31). We show our model in solid, green in Fig. 3.
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 B (P CDM ). In thick blue, dashed we show the Bullock et al. (2001) prediction for the FDM case, c B (P CDM ), with the modified collapse redshift, z coll (P FDM ) computed from the FDM linear theory power spectrum from AXIONCAMB . As in the case of   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, γ 1 and γ 2 . We show the effects within each prior, holding the other one fixed, by the light green (γ 1 ) and dark green (γ 2 ) shaded regions. We denote the mass scale M above which formally halos have not yet formed, i.e. z 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. the HMF, c B (P FDM ) is suppressed with respect to c B (P 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 provides a good fit to their simulations. In Eq. (33) the parameter γ 2 controls the slope of the cut-off. The second parameter, γ 1 , shifts the position of the cut-off scale. Note the slightly different form of the suppression function ∆ CDM c as compared to the HMF suppression function ∆ n . As above, M In thin, dashed, we show the prediction using the additional suppression formula Eq. (33), with the best-fit parameters γ 1 = 15 and γ 2 = 0.3. We see that the cut-off scale is moved by a factor ∼ γ 1 = 15 to the left compared to the scale M (c) 0 . To implement this additional suppression within HMCODE , we use the following prescription with 0 , slightly before c B (P CDM ) and c B (P 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 ∆ CDM c , 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 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 ST (P 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, δ 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 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: 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 ∆ FDM n and ∆ 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 σ > 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.

One-to-Two Halo Term Smoothing
HMCODE adopts a smoothing of the transition between the one halo term and the two halo term (Eq. 27), which captures aspects of 1.9 × 10 -24 8.5 × 10 -24 3.7 × 10 -23 Figure 5. Evolution of n eff as function of redshift z as compared to the CDM case, which is calibrated to simulations as described in Mead et al. (2015); Mead et al. (2021). On the far left we show a histogram of the values of n eff determined on the nodes of the calibration emulator. In the top right panel, we show in black the evolution of n eff in the case of CDM for our fiducial cosmology. Colour-coded by the mass m Ax we show the evolution of n 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 src , the normalised redshift distribution of the numbers of galaxies in four redshift bins i = 1 . . . 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 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 Mpc near the quasi-linear regime is smaller than the smallest measured angular bin at 2.8 for our fiducial cosmology. 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 ≈ 10 −24 eV (see Fig. A1, and also in the the Supplementary Material, Fig. A2). 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 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 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 eff simu-lated in the N-body simulations is n eff ≈ −2.6. Fig. 5 also shows the value of n 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 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. 7 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 ≈ 0.5 at which a scale of 1 Mpc is less than 2.8 , the smallest angular bin in the DES-Y1 data. The value of n 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 10 −23 eV, which, as we will see, covers all of the 95% confidence region of our constraints. Furthermore, the parameter n 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 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 eff (z) leaving the calibrated range, we computed the effect on the shear correlation function ξ + in the i-i-bins for two possibilities for n eff (and thus smoothing parameter α), shown in Fig. 6. The first takes α at its value specified by n eff according to the fits in HM-CODE . The second possibility we investigate applies a hard cap to α, 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 ∼ 3 × 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 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 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 eff : 1) leaving n eff free, as implemented in HMCODE , 2) capping n 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 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 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 ∑ m ν or n s , both of which would affect n 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 eff , which can be applied to any input linear power spectrum.

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 COSMO-SIS 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 crosscheck. 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 τ/N where N is the total number of samples (i.e. the number of samples per walker times the number of walkers) and where τ 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 τ using that estimate. An example comparison of the autocorrelation time estimates from the three methods described above can be found in Fig. 7.

DES Data Sets
Our analysis is based on the DES-Y1 3x2pt data set consisting of the observables ξ ± , γ 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 γ 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 8 , finding visual agreement in the twodimensional posteriors.
As discussed by Joudaki et al. (2020), the marginalised limits on the main parameters constrained by weak lensing, (σ 8 ,Ω m ), are prior dependent: the upper limit on σ 8 depends on how A s is sampled, and on the assumptions made about neutrino mass. Furthermore, the limit on Ω 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.

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 φ 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 φ 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 T T angular power spectrum and DES weak lensing arising from the lensing-induced smoothing of smallscale T T 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 10 12 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, ξξ + m = 1. × 10 -25 eV m = 2.7 × 10 -25 eV m = 3.2 × 10 -24 eV Figure 6. Discrepancy in the model predictions for ξ + (solid lines) and ξ − (dashed lines) for different treatments of the n 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 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 eff is capped to remain within the CDM calibrated region. The grey band (grey, dashed line) shows ξ + (ξ − ) 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. . 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 τ. while DES anchors smaller scales and thus constrains the nature of DM further.

Parameters, Priors, and Models
We adopt primary cosmological parameters: where A s , n s are the primordial power spectrum amplitude and spectral index, H 0 is the present day Hubble parameter, Ω b is the baryon density, and Ω m is the total matter density, Ω m = Ω b + Ω d and the DM (Ω d ) is composed of either CDM (for test cases) or FDM, but not both. Massive neutrinos are included with the parameter ∑ m ν for which we assume a single massive neutrino and N eff = 2.04 massless neutrinos. The neutrino mass is related to the neutrino density parameter by Ω ν h 2 = ∑ m ν /94. eV. For detailed discussion of how massive neutrinos are treated in the halo model, we refer to Mead et al. (2015); Mead et al. (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 Ω c h 2 and Ω 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 Ω 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 τ with a uniform prior, and the Planck absolute calibration parameter, A 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.
The axion mass is varied as: 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/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 Mpc −1 < k < 5 h Mpc −1 . Setting k J,eq (m) = 5 h Mpc −1 in Eq. (1) gives m = 1.5 × 10 −23 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-α forest data were sensitive. The disadvantage of a logflat prior is that there is infinite prior volume between any chosen upper limit and the CDM limit m → ∞. The alternative prior: has finite prior volume up to the CDM limit, 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, α, as described in Section 2.4.2.

RESULTS
Our baseline analysis uses the combination of Planck T T , T E, 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 Ω m = Ω d + Ω b , and the "matter fluctuation parameter" S 8 = σ 8 (Ω 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 Ω 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: log 10 (m/eV) ≥ −23.0 (95% C.L.) 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 ≈ 10 −22.5 eV which we attribute to the lower values of S 8 and Ω m measured by this survey in isolation. In the posterior for DES alone, we observe a peak probability at m ≈ 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 MULTI-NEST and EMCEEE and no corresponding peak is present with EM-CEE).
As discussed in Section 2.4.2, the two-halo smoothing term in HMCODE is controlled by a parameter, n 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 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 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 eff . As anticipated in our discussion in sec. 2.4.2, we find that the different models for n 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 eff on the posterior. For DES alone, we found a weak dependence on the n 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, demon-strating that our conclusions are not sensitive to the treatment of the smoothing term outside the range calibrated to simulation.
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 σ 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. 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/eV) and 10 −23 eV/m, respectively. When we calculate the 95% percentile of these distributions, we find a mass bound of m = 9.2 × 10 −24 eV for the log-flat prior and m = 6.7 × 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 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.

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 nonlinear 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/eV) ≥ −24.6 (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/eV) ≥ −23.0 (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-totwo 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 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, ξ ± . 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/eV) ≥ −22.8 (95% C.L.) (see the Supplementary Material).
We also considered varying the statistical aspects of our analy-  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. sis. 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.

Comparison to Other Results
Our model approximates FDM as having vanishing self-interactions, i.e. scalar potential V (φ ) ≈ m 2 φ 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 2 a [1 − cos (φ / 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 λ = 3M 2 pl / f 2 a , 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/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/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-α 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 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 ∼ 5% in a window near m ∼ 10 −27 eV, while the high-Ostriker-Vishniac anisotropies could reach mass fractions of 0.1% near m ∼ 10 −27 eV, and is sensitive to pure FDM up to masses as high as m ∼ 10 −22 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 ≈ 10 −23 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 10 −22 eV (Church et al. 2019), while the survival of the Eridanus-II star cluster excludes a range of masses near 10 −20 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 10 −21 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 σ 8 tensions.

Looking Forward
On small angular scales, θ 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, A1 and A2), 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 ∝ 1/θ for weak lensing correlation functions measured on an angular scale θ . In our analysis, the smallest value of θ was θ ≈ 4 arcmin in ξ + . The smallest scales used in the DES-Y3 analysis (Abbott et al. 2022) are around 1 arcmin. Using again that k J,eq ∝ m 1/2 (Eq. 1), substantial improvement in ULA constraints is possible with existing data, to m ∼ 10 −22 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 ∼ 10 −20 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 COS-MOSIS 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.

DATA AVAILABILITY STATEMENT
The data underlying this article, in the form of chains, will be shared on reasonable request to the corresponding author.  Figure A1. 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 ξ + observable on the left hand side and the ξ − 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 eV is indistinguishable from CDM over the scales shown.
servables. The profiles add to the non-linear power spectrum model in the form of the halo density autocorrelation ρ h ρ h , weighted by the HMF, in the 1h-term given in Eq. (9). Using Eq. (B1), we can factorise the autocorrelation as ρ h ρ h = ρ core ρ core + 2ρ core ρ env + ρ env ρ env . In the bottom row of Fig. B1, we show the halo density autocorrelation. We observe that the core-core correlation dominates the total correlation function up to r r α . The envelope-envelope correlation dominates for r r α , while the core-envelope term is subdominant. We notice that the FDM halo density autocorrelation for M = 5×10 10 M differs significantly from the CDM counterpart on its support r ≤ 2R vir , while the autocorrelations agree up to a few percent for the heavy halo with M = 1 × 10 14 M , except for very small r r α .
However, these differences in the halo auto-correlation between FDM and CDM for m = 10 −23 eV cannot be resolved within the DES survey. This is because the minimum angular bin is θ = 2.8 , which we can translate to the physical distance r in the autocorellation by means of the angular diameter distance D A = r/θ . For the source distribution in our sample, we in general expect the resolution to be highest for the nearest galaxies. We show in Fig. B1 the apparent angular size at redshift z = 0.15, with D A = 558 Mpc for our best fit cosmology (less than 5% of galaxies in the lowest bin of the DES data have a lower redshift). The corresponding region in the autocorrelation is marked as a green band. We observe in Fig. B1 that only the density correlation in the outer regions of large halos can be resolved, where the density profile is dominated by the envelope, which is very similar to the corresponding CDM case. This justifies our approach to use the unmodified CDM halo profile also for FDM halos in our implementation of the halo model, and in our analysis: the effects of solitons appear on small angular scales that are not resolved by DES-Y1.
So far, we have argued that the minimum angular resolution of the DES survey does not allow to measure the halo density autocorrelation of small halos for which the profile deviates significantly from the CDM counterpart. We conclude with the remark that even assuming much higher angular resolution, a detection of FDM features in halos using weak lensing is very challenging, because small FDM halos are suppressed by the additional cut-off in the HMF, as discussed in Section 2.3. We have considered the effects of solitons for m = 10 −23 eV, at the upper edge of our posterior distribution. For larger masses, solitons become smaller, and thus even harder to observe, while for lower masses solitons would lead to stronger deviations from CDM predictions in the observables.

B2 Correlation of the Smooth Component
The halo model, as described in Section 2.3.1, implicitly assumes that all DM is bound in halos. This assumption is valid for models such as CDM with no free streaming scale, where the variance of fluctuations diverges as the length scale in the filter R → 0. For models with a free streaming scale, including mixed DM with CDM and massive neutrinos, with pure WDM, or pure FDM (or any admixture of these models), there is cut-off scale below which some of the DM does not reside in halos. As described in Section 2.4, we implement this effect for FDM by modifying the HMF such that it is significantly suppressed below the cut-off scale. With our choice of  Figure A2. The γ and w correlation functions that form the "3x2pt" observable when combined with the ξ ± . The results obtained when using the full data set are compared to our baseline analysis in Fig. A3, and show a modest improvement in the bound on the axion mass.
the HMF Eq. (31), the fraction of the matter density bound in halos, f halo , can be calculated according to By accounting for the cut-off in the HMF we consistently model the one-halo term in our halo model. However, in addition, there now exists a component of DM unbound to halos, known as the smooth component. This implies that in addition to the two-halo term, we need to consider the auto-correlation of the smooth component as well as the cross-correlation of the smooth component and the matter contained in halos. The halo model can be adapted to account for such a smooth component: e.g. for CDM and neutrinos see Massara et al. (2014), while for pure WDM see Smith & Markovic (2011). The power spectrum with a smooth component is: where s denotes the smooth component, while h denotes the halo component, and we calculate the one-halo term P 1h and the twohalo P 2h term with respect to n(M) including the small scale cut-off. Note that we choose to normalise power spectra P 2h and P 1h as well as the cross-power spectrum P sh with respect to the total mean matter densityρ, not the matter density in halos ρ h , which eliminates the pre-factors of f halo as compared to e.g. Smith & Markovic (2011). As in Section 2.3.1, we assume that halos are biased tracers of the linear matter density. It is furthermore reasonable to assume that also the smooth component traces the linear matter density (Smith & Markovic 2011), and that both the smooth and the halo component can be related to the linear field through the linear bias terms b s (M) and b h (M), respectively.
In the limit that the bias is ignored, the smooth-smooth, smoothhalo, and two-halo terms in the power spectrum all add up to equal the linear power spectrum regardless of the value of f halo (Smith & Markovic 2011). Thus, in the limit that the bias is neglected in the two halo term, one can neglect the smooth component correlations. HMCODE neglects the bias since it is not necessary to compute the non-linear within the desired accuracy: the one-halo term dominates before the bias has a significant effect. Thus it is consistent within the approximations used by HMCODE to neglect the smooth component correlations of FDM, as long as the one-halo term has the correct cut-off. Furthermore, we showed that within the accuracy of DES-Y1 the HMF cut-off does not affect our lower limit to m Ω m 68%, 95% CL DES ξ ± + Planck18 DES 3×2pt + Planck18 Figure A3. Posterior density for the DM particle mass m, the mass fluctuation parameter S 8 , and the matter density Ω 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. (Fig. 10), which further reinforces the claim that smooth component correlations can be neglected within the desired accuracy.

B3 Relativistic Corrections
ULAs are treated in our work at early times fully relativistically in linear perturbation theory within AXIONCAMB . The energy density and effective pressure contains contributions from the kinetic, gradient, and potential energy density of the FDM field, which all enter into the computation of the linear power spectrum.
Our modifications to the halo model, i.e. the cut-offs in n(M) and c(M) are fit to N-body simulations, which take the linear power spectrum as input, but treat the DM entirely non-relativistically. The next level of improved approximation to FDM physics on non-linear scales involves the introduction of the quantum pressure term, and solution of the Schrödinger-Poisson equations. The equations differ from an N-body model, since the equation of motion for the field contains gradient energy, an effective smoothing scale in the Vlasov equation (Widrow & Kaiser 1993;Uhlemann et al. 2014). The quantum pressure term is expected to lead to extra suppression in n(M) and c(M) in addition to what is observed in N-body simulations (although no large enough simulations have been performed including this physics to make a definitive measurement), and we have accounted for this in our model for systematic uncertainty in n(M) and c(M). Quantum pressure also leads to the formation of solitons on small scales in FDM halos, which as we have argued above can be neglected on scales relevant to the DES-Y1 shear correlation.
There are, however, aspects of quantum pressure which are ne-glected even in the most advanced cosmological simulations of the Schrödinger-Poisson equations (Schive et al. 2014a;Veltmaat et al. 2020): namely, relativistic corrections. A full relativistic cosmological simulation of FDM or ULAs is unfeasible, although Numerial Relativity and other semi-relativsitic models can be used to simulate the relativistic collapse of solitons (Helfer et al. 2017;Levkov et al. 2017;Michel & Moss 2018). The present work neglects ULA self-interactions, which is a valid approximation in the mass range of interest, where obtaining the correct relic density requires a decay constant f a 10 17 GeV (Marsh 2016b). In this limit solitons undergo collapse to black holes when they exceed the Kaup (1968) critical mass (Helfer et al. 2017). Solitons formed in DM halos by gravitational hierarchical structure formation and obeying the core-halo mass relation (Schive et al. 2014b) can be shown to remain always below this critical mass for even the most massive observed halos in the Universe, and thus strong gravity effects can safely be neglected.
A relativistic correction to the Schrödinger-Poisson equations that might impact cosmological structure formation is the gravitation of FDM gradient energy, which appears as a correction to the right hand side of the Poisson equation. This correction, and other relativistic corrections, have been studied by Salehian et al. (2021). The leading relativistic corrections are found to lead to smaller solitons than the non-relativistic Schrödinger-Poisson solutions considered in the previous section. Thus we can conclude that our approximation to neglect solitons in the halo model based on the non-relativistic model is a conservative one.  Figure B1. Halo density profiles ρ (top) and auto-correlation functions ρ ρ (bottom) for two different halo masses, for FDM mass m = 10 −23 eV (red) and CDM (blue). We compare a less massive halo M = 5 × 10 10 M (left) to a more massive halo M = 1 × 10 14 M (right). In the case of CDM, we show the NFW profile given in Eq. (17), which corresponds to setting the parameter η = 0 in HMCODE (i.e. no halo bloating). In the case of FDM we show a piece-wise profile with a solitonic core and an envelope parameterised by a modified NFW profile, as given in Eq. (B1). The matching radius r α is set to r α = 2 in the left-hand column, similar to the values found in Robles et al. (2018). In the right-hand column, we chose r α = 1.2 which leads to the same concentration-mass relation for the NFW profile and the profile including the solitonic core. In the top row, we show in dotted (dashed) what the ULA profile would look like if we extrapolated the core (envelope) profile beyond the matching radius. In the bottom row, we depict the full correlation function in solid. For FDM we show in dotted, dashed, and dot-dashed the contributions of the core-core, envelope-envelope and core-envelope correlation functions to the full result (see text for details). To guide the eye, we display the angular scale of the halo density auto-correlation for angular diameter distance D A = 558 Mpc (roughly 5% of the galaxies in the DES-Y1 sample, at z = 0.15). The green shaded region shows the angular resolution of DES. Planck18 DES ξ ± DES ξ ± + Planck18 Figure C1. Lack of degeneracy between massive neutrinos and ULA mass in our analysis.

APPENDIX C: MASSIVE NEUTRINOS AND OTHER SOURCES OF POWER SUPPRESSION
Massive neutrinos, due to free streaming, also suppress P(k), leading to step-like features. For observables only sensitive to nonrelativistic structure formation, this leads to a degeneracy between the effects of massive neutrinos and ULAs (Amendola & Barbieri 2006;Marsh & Ferreira 2010). For neutrinos in the Standard Model of particle physics for the known cross section, neutrinos freeze-out while they are non-relativistic, and the density and mass are not independent parameters, Ω ν h 2 ∝ m ν . Imposing a rough upper limit on the neutrino mass of 1 eV, the power suppression by ULAs only becomes degenerate with neutrinos for m < 10 −28 eV, corresponding to field oscillations beginning in the matter-dominated era ("DE-like ULAs", Hlozek et al. 2015). The degeneracy in P(k) only occurs when ULA fraction is chosen to mimic the P(k) step amplitude of the neutrinos, which demands Ω a 0.12; alternatively the effective number of neutrinos can be varied along with the ULA density to increase the neutrino density parameter while holding m ν fixed (Marsh et al. 2012;Hložek et al. 2017). Once relativistic observables, in particular the CMB anisotropies, are considered, the degeneracy between massive neutrinos and ULAs is broken. This is due to the different behaviour of ULAs and neutrinos in the relativistic regime, where neutrinos have equation of state w = 1/3 (hot relativistic particles) and ULAs have w = −1 (slowly rolling scalar field), and their consequent differing effects on the expansion rate of the Universe, which leads to differing effects in both the Sachs-Wolfe and Silk damping regions of the CMB power spectrum (Hlozek et al. 2015;Hložek et al. 2017).
In the present analysis, we demand that ULAs are all of the DM, which is strongly disfavoured by the data for low masses m < 10 −28 eV where the P(k) degeneracy with neutrinos opens up. This is demonstrated in in Fig. C1, where we show the joint posterior on Ω ν h 2 and log 10 (m/eV). DES data alone does not provide a strong constraint on the neutrino mass (DES-Y1, Abbott et al. 2018a), but demands ULAs have m 10 −23 eV, which cannot be degenerate in P(k) for neutrinos with m ν 1 eV. Including CMB data gives a strong upper bound on the sum of neutrino masses (Aghanim et al. 2020b) and on the ULA mass when all the DM is ULAs. In the allowed parameter space for the CMB only analysis, ULAs are too heavy to be degenerate with neutrinos. This state of affairs can be understood intuitively: the CMB demands that neutrinos are hot, and become non-relativistic at late times. On the other hand the CMB demands that DM be pressureless at matter radiation equality. Consequently the neutrino free-streaming scale and the ULA Jeans scale are separated by many orders of magnitude and the effects on the matter power spectrum are not degenerate.