A weak lensing perspective on nonlinear structure formation with fuzzy dark matter

We investigate nonlinear structure formation in the fuzzy dark matter (FDM) model in comparison to cold dark matter (CDM) models from a weak lensing perspective using perturbative methods. We use Eulerian perturbation theory (PT) up to fourth order to compute the tree-level matter trispectrum and the one-loop matter spectrum and bispectrum from consistently chosen initial conditions. In addition, we predict the non-linear matter power spectra using N -body simulations with CDM and FDM initial conditions. We go on to derive the respective lensing spectra, bispectra and trispectra in CDM and FDM in the context of a Euclid-like weak lensing survey. Finally, we compute the attainable cumulative signal-to-noise ratios and an estimate of the attainable χ 2 -functionals for distinguishing FDM from CDM at particle masses m = 10 − 21 eV, m = 10 − 22 eV and m = 10 − 23 eV. We find that PT predictions cannot be used to reliably distinguish the three models in a weak lensing survey. Assuming that N -body simulations overestimate the late-time small-scale power in the FDM model, future weak lensing survey might be used to distinguish between the FDM and CDM cases up to a mass of m = 10 − 23 eV. However, observations probing the local high-z universe are probably more suited to constrain the FDM mass.


INTRODUCTION
The Fuzzy Dark Matter (FDM) model first proposed by Hu et al. (2000) describes dark matter as a bosonic, scalar field composed of very light particles with typical masses m ∼ 10 −22 eV with negligible self interactions and a macroscopic de Broglie wavelength.The dynamics of FDM in the non-relativistic limit are governed by the Schrödinger-Poisson system (SPS) of equations.Because of its small mass, the FDM boson has an astrophysically relevant de Broglie wavelength at the order of a few kiloparsecs.Its wave-like behaviour suppresses structure formation on small scales while one recovers CDM behaviour on large scales; as such FDM constitutes a small-scale modification of CDM, governed by the particle mass.Therefore, the FDM model has the potential to solve the small-scale crisis of CDM.
The FDM bosons are also called axions or axion-like particles.Axion-like particles with exponentially suppressed masses are naturally generated in supersymmetric theories and theories with extra-dimensions including string theory (Marsh 2016).A currently popular FDM model is one where all of dark matter is composed ⋆ e-mail: alexanderkunkel@ntu.edu.tw† e-mail: bjoern.malte.schaefer@uni-heidelberg.de of axions with m ≈ 10 −22 eV (Hu et al. 2000;Marsh 2016).This mass has also been suggested by Schive et al. (2014) who fitted the ground state density of a halo obtained in a numerical simulation of FDM to the mass distribution of the Fornax dwarf spheroidal galaxy.However, the allowed range of FDM masses that could solve the small-scale crisis of CDM is increasingly narrowed down by a number of different observations: Kobayashi et al. (2017); Iršič et al. (2017); Armengaud et al. (2017) constrained the FDM mass using measurements of the Lyman-α forest flux spectrum from the XQ-100, HIRES/MIKE, SDSS/BOSS and XSHOOTER data sets and derive bounds between m > 7.1 × 10 −22 eV and m > 2.3 × 10 −21 eV.The biggest uncertainties in their approach stem from the modelling of the observed fluctuations in the neutral hydrogen and the use of hydrodynamical CDM simulations using FDM initial conditions (IC), thereby neglecting FDM dynamics at late times.Nori et al. (2018) improved on the simulation side by including the effect of the quantum pressure term in their N-body simulations and derive the bound m > 2.1 × 10 −21 eV.Rogers & Peiris (2021) found a stronger bound of 2 × 10 −20 eV at 95% confidence using a new, robust N-bodybased modelling pipeline marginalising over a IGM model which allows for a wide range of heating and ionization histories.The fact that the different groups' results are consistent indicates that the effect of late-time FDM dynamics on the Lyman-α flux power spectrum may be negligible.At the same time, Zhang et al. (2019) carefully examined the Lyman-α forest constraints on FDM and concluded that full FDM simulation are required for reliable constraints based on the Lyman-α forest.They also pointed out that the simulation uncertainties of hydrodynamic simulations may have been underestimated in previous studies.They nonetheless find that FDM with masses smaller than m = 2.5 × 10 −23 eV can be rejected based on the Lyman-α forest.
Order unity density fluctuations as a result of wave interference in FDM halos are a robust and well-understood prediction of the FDM model.They are expected to dynamically heat stellar orbits in ultra-faint dwarf galaxies (UFDs).Marsh & Niemeyer (2019) compute the gravitational heating of a star cluster in the UFD Eridanus II and find that the range 0.8×10 −21 eV < m < 10 −19 eV is disfavoured.Further, they find that the formation of Eridanus II as a subhalo demands m ≳ 0.8×10 −21 eV.Schive et al. (2020) use FDM simulations to demonstrate that the random walk which the central soliton in the FDM halo undergoes represents an even bigger challenge for the survival of the central star cluster in Eridanus II.As a possible solution, they argue that the star cluster could possibly survive the soliton random walk if the FDM subhalo had been tidally stripped by the influence of the Milky way's gravitational potential.Dalal & Kravtsov (2022) derive the bound m > 3×10 −19 eV by showing that the dynamic stellar heating predicted by FDM is not consistent with the observed stellar velocity dispersion in the UFDs Segue 1 and Segue 2. Amorisco & Loeb (2018) study the dynamical heating of thin stellar streams through the density granulation in FDM halos and find the limit m > 1.5 × 10 −22 eV.In a similar spirit, Church et al. (2019) study the dynamical heating of galactic disks in CDM and FDM and find the lower bound m > 0.6 × 10 −22 eV.
The inner density profile of dwarf galaxies provides another test bed for FDM predictions.Chen et al. (2017); Calabrese & Spergel (2016) compare the flat, inner density profile of Milky way dwarf satellites to the inner solitonic profile of FDM halos and find m ≈ 1 − 6 × 10 −22 eV.Safarzadeh & Spergel (2020) point out that while the FDM prediction for m ≳ 10 −21 eV matches the profiles of the known UFD Milky Way satellites, the density profiles of the dwarf spheroidals Fornax and Sculptor originally considered by Schive et al. (2014) disagrees with the FDM prediction.In turn, a mass around 10 −22 eV seems to suggest that the mass of the UFDs is too large.In other words, no single FDM mass seems to be able to accomodate for the density profiles of all known dwarf galaxies.
Another constraint on the FDM mass can be derived from the circular velocity in disk galaxies.Assuming the soliton-host halo mass relation empirically found by Schive et al. (2014), Bar et al. (2018) argue that the peak circular velocity characterising the host halo on large scales should repeat itself in the central region and rule out FDM masses in the range m = 10 −21 −10 −22 eV.The biggest uncertainty in their approach stems from the empirical soliton hosthalo relation that has only been studied numerically for a limited range of FDM and halo masses.
Using the subhalo mass function inferred from the observed luminosity function of Milky Way Satellites, Nadler et al. (2021) constrain the warm dark matter, self-interacting dark matter and fuzzy dark matter models alike and obtain the bound m > 2.9 × 10 −21 eV at 95 % on the FDM mass.Uncertainty comes from the prediction of halo mass function from the linear power spectrum.A further mass constraint based on early halo formation is given by Lidz & Hui (2018) who use the hyperfine 21 cm transition of hydrogen at redshift around z = 15-20 to constrain the FDM mass to m > 5 × 10 −21 eV.
Future surveys of the 21-cm power spectrum such as HERA are expected to be sensitive to FDM masses up to 10 −18 eV if all dark matter is fuzzy as well as to very small FDM fractions around 1 % in the mass range 10 −25 eV ≲ m ≲ 10 −23 eV (Hotinli et al. 2022).
Gravitational lensing provides yet another complementary probe in addition to the above constraints.Powell et al. (2023) study the granular structure in the main dark matter halo of a single gravitational lens system and derive the bound m > 4.4 × 10 −21 eV.For lower masses, Hložek et al. (2018) identify the possible mass range for FDM via CMB lensing based on the full Planck data set and find no evidence for an FDM component in the mass range 10 −33 eV ≤ m ≤ 10 −24 eV.More recently, Rogers et al. (2023) combine the Planck CMB data and galaxy power spectrum and bispectrum data from the BOSS survey to conclude that fuzzy dark matter with m ≤ 10 −26 eV makes up less than 10 % and fuzzy dark matter with 10 −30 eV ≤ m ≤ 10 −28 eV makes up less than 1 % of dark matter today.
To sum up, essentially the entire mass range from m = 10 −26 eV to m = 10 −16 eV is under tension with a number of different observations.The m ≈ 10 −22 model is particularly constrained by the Lyman-α forest, galactic rotation curves, dynamical heating arguments, the subhalo mass function and strong lensing.
In this work, we estimate whether future weak lensing surveys have the constraining power to further complement some of the existing constraints on the FDM mass.Weak lensing surveys such as Euclid measure the lensing shear spectrum and allow to infer the shape and amplitude of the matter spectrum P δ , thereby giving constraints on cosmological parameters such as the axion mass.Gravitational lensing probes the matter distribution without any biasing assumption.This is important because the preferred value of the galaxy bias b and therefore the normalisation of the galaxy spectrum are different for axion cosmologies than for ΛCDM.As a consequence, the constraining power of galaxy surveys on the FDM mass is lower than that of weak lensing.
In previous weak lensing studies on FDM, Marsh et al. ( 2012) analysed whether adding a small fraction of axions of mass in the range m = 10 −29 eV would be detectable via the lensing convergence spectrum.For modelling nonlinearities, they neglected the quantum pressure and employed the CDM halofit model implemented in the CAMB code (Lewis et al. 2000).More recently, Dentler et al. (2022) combined CMB Planck data with shear correlation data from the Dark Energy Survey year 1 to find a 95% C.L. lower limit m > 10 −23 eV.They modelled the nonlinear FDM spectra using the adapted halo model HMCODE.Lensing can also serve as a tool to investigate axion-related isocurvature fluctuations (Feix et al. 2019(Feix et al. , 2020)), with a sensitivity to axion masses at the scale 10 −19 eV.In contrast, we estimate the spectrum, bispectrum and trispectrum using fluid perturbation theory.To this end, we employ a recasting of SPS into a hydrodynamical form using the so-called Madelung transform.The Madelung equations are Euler-Poisson equations with an additional scale-dependent modification, the socalled quantum pressure term.On large scales, this term vanishes and one recovers the ideal fluid equations for CDM.One of the appeals of the Madelung transform lies in being able to apply standard cosmological perturbation theory to the SPS.Its range of validity is higher than that of wave perturbation theory and one can easily contrast it with perturbation theory in CDM (Woo & Chiueh 2009;Li et al. 2019).Whether the inherently higher degree of small-scale anisotropy in FDM (Dome et al. 2022) can differentiate between dark matter models remains to be investigated.
Building on and extending the work of Li et al. (2019) who computed the one-loop matter spectrum for FDM, we use timedependent nonlinear Eulerian perturbation theory up to fourth order to compute the tree-and loop-level bispectrum as well as the tree-level trispectrum.To estimate the range of validity of our results, we compare the perturbation theory predictions with a set of N-body simulations with CDM and FDM initial conditions.We go on to derive the corresponding lensing spectra for a Euclid-like lensing survey.We estimate the attainable signal-to-noise ratios as well as the χ 2 -functional for distinguishing axions of the masses m = 10 −21 eV, m = 10 −22 eV and m = 10 −23 eV from standard CDM.These masses are interesting because they cover the range of masses that could possibly solve the small-scale crisis of CDM.
The outline of this paper is as follows: Section 2 introduces the SPS, presents linear and nonlinear Eulerian perturbation theory for FDM and computes the matter spectra, bispectra and trispectra.Section 3 introduces the basics of weak lensing surveys, computes the respective lensing spectra and estimates the attainable weak lensing signal-to-noise ratios and χ 2 -functionals.Section 4 summarises and discusses the results.Appendices A through C list the technical and numerical details of the FDM perturbation theory.Appendix D lists the details of the N-body simulations.

NONLINEAR STRUCTURE FORMATION WITH FUZZY DARK MATTER
Consider a (pseudo-)scalar field φ minimally coupled to gravity: where we follow the convention in Hui et al. (2017).m is the axion mass, which naturally appears in the Compton-scale λ = ℏ/(mc).This action is invariant under parity-and time-inversion because it is quadratic in φ.For QCD axions, this action is applicable after symmetry breaking and after non-perturbative effects have been switched on.Further, we neglect possible non-gravitational selfinteractions of the axion in the form of an axion potential V(φ).In our discussion, it will not be important that the FDM particle is actually an axion or an axion-like particle.We only assume that the particle is bosonic, non-relativistic, has negligible self-interaction and makes up the entirety of dark matter.We can study solutions of the relativistic axion equation of motion for a perturbed Friedmann-Lemaître-Robertson-Walker (FLRW) metric including a cosmological constant Λ, all in Newtonian gauge.This is appropriate for studying structure formation because the virial velocity in a typical galaxy υ vir ∼ 100 km s ≪ c and galaxies are much smaller than the Hubble horizon.On the lower scale end, we are concerned with scales above the axion Compton wave length which would correspond to relativistic scales in the Klein-Gordon equation.Except in the vicinity of black holes, the Newtonian potential Φ obeys |Φ| /c 2 ≪ 1.To study the clustering of axions on nonlinear scales, we further take the WKB-approximation of the form where ψ is a complex scalar field because axions that cluster on galactic scales began oscillating in the very early universe.We apply the previous considerations by taking Φ ∼ ϵ 2 , k/m ∼ ϵ and H/m ∼ ϵ and work to order O(ϵ 2 ).Further, we assume the nonrelativistic approximation | ψ| ≪ mc 2 ℏ |ψ|.The assumption ∂ t ≪ m is non-relativistic because we have ∂ t ∼ ∆/m ∼ k 2 /m and therefore k 2 /m ≪ m.With these simplifications, we obtain the comoving Schrödinger-Poisson system of equations where ρ b (t) is the background density and |ψ| 2 measures the density in a proper volume.We supplement the Schrödinger-Poisson equations with the normalisation condition which fixes the density ρ = |ψ| 2 for N axions of mass m.In Eq. ( 3), positions of particles are described using comoving coordinates x.

quantum pressure and Jeans scale
The Madelung transform (Madelung 1927) follows by substituting ψ(x, t) =: ρ(x, t) m e iS (x,t) (6) for real fields ρ(x, t) and S (x, t) into Eq.( 3) and defining the velocity field we obtain the Madelung equations The Madelung equations describe the Schrödinger equation via a system of fluid equations for frictionless, compressible flow in an external gravitational potential Φ.The flow gets modified by the quantum pressure The quantum pressure accounts for the underlying wave dynamics in FDM.For a narrowly located source, the quantum pressure is large and reflects the Heisenberg uncertainty principle in quantum mechanics 1 .On large scales, the Madelung equations reduce to the Euler equations of a pressureless fluid and we recover the dynamics of standard cold dark matter.The quantum pressure Q P is equivalent to an anisotropic pressure stress P i j where We can describe the linearised evolution of the density contrast by assuming δ ≪ 1 and |υ| ≪ 1 for the fluctuation fields and neglecting higher-order perturbations of the form O(δ 2 , υδ, υ 2 ) to obtain For each mode k, Eq. ( 11) describes a harmonic oscillator with time-dependent dampening H(t) and frequency ω(k, t) as Unlike in the CDM case, linear growth in FDM is scale-dependent because of the quantum pressure term.The condition ω = 0 defines the comoving quantum Jeans scale k J (see Khlopov et al. 1985): The Jeans scale describes a force balance between gravity and quantum pressure.For k < k J , i.e. scales larger than λ J , ω(k, t) becomes imaginary and we recover a growing and a decaying mode just like for CDM.Perturbations on these scales are unstable and will gravitationally collapse.For k > k J , i.e. scales smaller than λ J , the frequency ω(k, t) becomes real.Perturbations on small scales therefore undergo oscillations and are stabilised against gravitational collapse.The comoving Jeans wavelength decreases with time as λ J ∼ a − 1 4 which is why more small-scale features can develop as the universe expands.We can find an analytical solution of the linear FDM growth equation (11) in the Einsteinde Sitter case for Ω m = 1.To this end, we substitute the ansatz δ(x, a) = δ(x, a = a 0 )D(a, a 0 ) into Eq.( 11) and rewrite the resulting ODE for the growth factor D in terms of the scale factor a Simplifying with Ω m,0 = 1 yields Chavanis ( 2012) and Suárez & Chavanis (2015) study solutions of Eq. ( 15) in terms of Bessel functions in the context of Bose-Einstein condensate dark matter with self-interaction and the weakfield limit of the Klein-Gordon-Einstein equation.They find where J ∓5/2 are cylindrical Bessel functions of fractional order Li et al. (2019) give a solution of Eq. ( 15) with D ± (k, a 0 , a 0 ) = 1 in terms of the solutions Eq. ( 16) as We immediately notice that Eq. ( 19) can exhibit unphysical divergences at the roots of the Bessel function.Laguë et al. (2020) argue that this is simply a matter of choice of normalisation, but this is not quite true.Eq. ( 19) is not a solution of the linear growth equation ( 15) if it diverges in the domain of interest.We are not aware of a general analytical solution to Eq. ( 15), even in the EdS case.Laguë et al. (2020) suggest two ways to remedy this issue: One way is to approximate the denominator of Eq. ( 19) by a fifth-order Taylor expansion of the Bessel function which removes divergences from Eq. ( 19) and gives fewer oscillations.In contrast, a different number of terms in the Taylor expansion leads to fast oscillations and/or divergences.An alternative that Laguë et al. (2020) adopt is to use a model for the mean growth of the FDM growth function.They describe the growing mode D + in terms of a smoothed Heaviside step function where free parameters are determined via a fit to the axion-CAMB transfer function (Hlozek et al. 2014).For large scales, this approach exactly recovers the CDM growth function, but has the disadvantage that oscillations are entirely neglected.Moreover, the smoothed Heaviside step function falls off exponentially for large k which is not a correct description of the asymptotic behaviour of D + .We therefore opt for integrating Eq. ( 15) numerically.As initial conditions, we choose where D ′ CDM (a 0 ) is obtained via numerical integration of the linear growth equation ( 14) for ℏ = 0 corresponding to CDM.This approach has several advantages: We ensure D(k, a 0 , a 0 ) = 1 via the initial conditions and obtain the correct CDM evolution in a general cosmology in the limit k → 0. Further and most importantly, the growth factor obtained in this way is in fact a solution to the linear growth equation (15).It does not exhibit unphysical divergences as shown in Fig. 1.We see that the numerical solution of the growth equation also exhibits the correct asymptotic behaviour lim k→∞ D FDM (k,a) D CDM (a) = 0.In addition, the numerical solution captures the oscillatory behaviour of the analytical solution for high k.At the same time, there are also several disadvantages to numerically integrating the linear growth equation in FDM.Firstly, we do not capture the growing mode in the oscillating regime.This is because we do not know the correct initial conditions for the growing mode.Therefore, the numerical solution in the oscillating regime will in general be a linear combination of the two modes D + and D − .However, since the correct initial conditions for the linear fluctuation fields are unknown as well, this does not add any uncertainty to the linear growth model in terms of initial conditions.In any case, we recover the correct modes for a > a osc .This is because any component proportional to D − in the initial conditions quickly decays away for a > a osc .Fig. 2 shows the growing solutions obtained for a fiducial cosmology with Ω m,0 = 0.3159 and a dark energy equation of state w = −0.9used in the rest of this work.
To sum up, the key difference between CDM and FDM in linear perturbation theory is the existence of a unique length scale in  19), the renormalised expression using the fifth-order Taylor expansion, the numerical solution obtained by integrating the linear growth equation with initial conditions given by Eq. ( 20) as well as a with the fit parameters α = 0.17, β = 6.50 obtained via fitting to the analytical solution.FDM.Whereas all scales are gravitationally unstable in CDM, the perturbations below the quantum Jeans scale are stabilised in FDM.Nonetheless, nonlinear perturbation theory may alter these conclusions since the quantum Jeans scale is a concept only valid within linear perturbation theory.This can be seen by taking the first few terms of the Taylor expansion of the quantum pressure Eq. ( 9): ) .The linear contribution −∆ 2 δ counteracts gravity.However, the quadratic terms acts in the same direction as gravity and could therefore potentially enhance gravitational collapse (Li et al. 2019).This is related to the fact that the interference of waves can lead to structures that are smaller than their wavelength.In order to estimate the effect of nonlinearities on cosmological structure formation, we develop nonlinear perturbation theory for CDM and FDM in the next section.

tree-and loop level perturbation theory
The goal of cosmological perturbation theory is to describe the departure of matter evolution from the homogeneous Hubble expansion perturbatively.In Eulerian perturbation theory (PT), one describes the nonlinear gravitational dynamics in terms of solutions δ (1) and υ (1) of the linearised fluid equations in a fixed laboratory frame: where θ = ∇ • υ and the n th -order fluctuation fields δ (n) and θ (n) are proportional to the n th power of the linear fluctuation fields: In addition, we have ∂ a δ (1) = θ (1) from the linearised continuity equation.Therefore, both the velocity and the density field are fully determined by the linear density fluctuations.The coupling of linear modes in the nonlinear theory is described by the nonlinear coupling kernels F n and G n : where F n and G n are homogeneous functions of the wave vectors q 1 , . . ., q n with degree zero and q 1,...,n ≡ q 1 + . . .+ q n .A series of papers (Fry 1984;Goroff et al. 1986;Jain & Bertschinger 1994) developed a method for deriving the time-independent CDM mode coupling kernels in terms of algebraic recursion relations.They rely on the Euler-Poisson system being homogeneous in the scale factor a in the EdS case.Including the quantum pressure breaks this homogeneity and just like for a general cosmology in CDM, the solutions at each order become non-separable functions of time and scale.
A method to nonetheless obtain the PT kernels in this case is described in the series of papers (Scoccimarro 1998(Scoccimarro , 2006) that developed time-dependent Eulerian PT with the aid of Feynman diagrams.Effectively, Li et al. (2019) use this method to compute the kernels F 2 and F 3 in FDM.The nonlinear terms in the Madelung equations are treated as inhomogeneity g(η) while solving the linear growth equation.This leads to an integral equation that can be represented as a Dyson series.The Dyson series allows for a diagrammatic representation and can be recursively solved up to a given order in PT.In appendix A, we retrace the steps taken by Li et al. (2019) and extend them to higher order in perturbation theory for consistent tree-and loop-level computation of bi-and trispectra in FDM, which necessitates the kernels F 2 , F 3 and F 4 ; recent developments in perturbation theory recast this straightforward perturbative approach into more powerful formalisms (Pietroni 2008;Bartelmann et al. 2016;Kozlikin et al. 2021) including FDM (Littek 2018).

spectra, bispectra and trispectra at tree-and loop-level
We derive explicit expressions for the kernels F 2 , F 3 and F 4 in the FDM and CDM case in appendix A. These kernels F n can then be used to perturbatively expand n-point-matter correlation functions P n (k 1 , k 2 , . . ., k n ).The latter are defined as Fourier transform of the connected correlation functions: In the following, we are interested in the equal-time matter spectrum P ≡ P 2 , bispectrum B ≡ P 3 and trispectrum T ≡ P 4 .We express the n-point correlations of non-Gaussian perturbations as integrals over higher-order correlations of Gaussian perturbations.
We then apply Wick's theorem to express higher-order correlations of Gaussian perturbations as two-point correlations of Gaussian perturbations.The latter are completely specified by the initial conditions.Expanding the two-point correlation as a perturbation series P(k, a) = P (0) (k, a) + P (1) (k, a) + . .., we find that the time evolution of the linear spectrum P (0) can be expressed as linear scaling of the initial spectrum: In the following, we omit the explicit time dependencies.The first higher-order contribution P (1) to the spectrum comes at loop-level 2 .The two loop-level contributions to the one-loop spectrum P (1) are given by with P 13 (k) = 6 d3 q F (s) 3 (k, q, −q)P (0) (k)P (0) (q).( 29) The bispectrum at tree-level corresponding to second-order PT is given by The one-loop contribution to the bispectrum consists of four distinct diagram involving up to fourth-order PT kernels: 2 Note that a consistent truncation of series in PT is obtained by including terms up to a certain power m in the linear spectrum.This corresponds to grouping the PT contributions in terms of the number of loops in a diagrammatic representation.
whose explicit expressions are given by In the following, we also examine the reduced bispectrum It is independent of time and normalisation.For scale-free initial conditions P (0) , i.e.P (0) ∝ k n with the spectral index n, Q (0) is also independent of overall scale and for equilateral configurations it is also independent of the spectral index.Finally, the diagrams for the trispectrum involve vertices connecting two and three lines and therefore correspond to second-and third-order PT.One can decompose the tree-level trispectrum into the contributions where We now present the different spectra computed with timedependent PT in CDM and FDM for a dark energy cosmology with w = −0.9 and Ω m,0 = 0.3159 at z = 0.For computing the nonlinear PT kernels F n , we use the EdS growth factors in both CDM and FDM 3 .The initial spectra are computed using the CAMB and axionCAMB codes for CDM and FDM respectively (Lewis et al. 2000;Hlozek et al. 2014).Importantly, the nonlinear CDM matter power spectra obtained via N-body simulations and via PT are only expected to agree to percent-level up to k ∼ 0.1 h/Mpc at z = 0 (Carrasco et al. 2014).Therefore, we also run a total of 8 N-body simulations in two boxes of the sizes L = 256 Mpc/h and L = 30 Mpc/h with N = 512 3 particles using the GADGET-2 code (Springel 2005).The L = 256 Mpc/h box provides a crosscheck of the large-scale loop-level PT matter power spectra.The L = 30 Mpc/h simulations are used to obtain a rough estimate of the full nonlinear power spectra on small scales.Since we probe the matter power spectrum at late times and on scales smaller than the Jeans length, we expect the effect of the quantum pressure to be important.Armengaud et al. (2017) note that especially for masses smaller than m ∼ 10 −22 eV, the quantum pressure affects a significant portion of the dark matter particles on scales k ≳ 1 h/Mpc.Therefore, the N-body power spectra for the m = 10 −22 eV and m = 10 −23 eV cases are expected to exhibit large errors on small scales that could only be eliminated with large-scale FDM simulations including wave dynamics which are arguably beyond the reach of current supercomputers.2018) also compare N-body simulations with and without quantum pressure at z = 0 and find that the quantum pressure reduces small-scale power.Therefore, we believe that the N-body simulations serve as an important cross-validation.Technical details can be found in appendix D.
Fig. 3 shows the CDM and FDM spectra at tree-level and with loop-level corrections as well as the spectra from several N-body simulations with FDM and CDM IC at z = 0.At tree-level, power is strongly suppressed below the Jeans scale in the FDM model.Nonlinear corrections at loop-level transfer power to small scales, but suppression is still dominant.As previously noted by Nori & Baldi (2018); Nori et al. (2018), the initial suppression of power on small scales in the N-body simulations is almost completely overcome by late-time nonlinear CDM dynamics.Note that the low-k modes in the N-body simulations suffer from large variations that stem from a low number of samples.Moreover, there is a mismatch between the high-k modes of the L = 256 Mpc/h and the low-k modes of the L = 30 Mpc/h simulations at k = 6 h/Mpc indicated by the vertical gray line due to the fact that the L = 30 Mpc/h box is too small for the large-scale modes to evolve correctly.Fig. 4 shows the respective spectra at z = 6.1 for reference.The looplevel PT and the N-body results agree up to k = 1 h/Mpc, but still start to differ on smaller scales.The different N-body power spectra demonstrate that late-time CDM dynamics have not yet overcome the initial suppression of power on small-scales at z ∼ 6.
Figs. 5 and 6 show the equilateral matter bispectra and equilateral square matter trispectra and Figs.11 and 12 the respective convergence spectra.As in the case of the spectrum, loop-level corrections for the bispectrum have the effect of adding power on small scales in both CDM and FDM.At both tree-and loop-level, however, suppression below the Jeans scale is still the dominant effect in FDM.
The loop-level corrections given by Eq. ( 27) were computed in a form free of infrared divergences using the CUBA-library Hahn (2005).Details of the numerical integration can be found in appendix C. Fig. 7 shows the angular dependence of the reduced matter bispectrum Q (0) at tree-level.The fact that Q (0) is enhanced for θ = 0, π reflects the fact that large scale flows generated by gravitational instability are mostly parallel to density gradients.As discussed in section 2.1, the kernel F (s)  2 includes the first higher-order correction from the quantum pressure term.Therefore, it does not only counteract gravitational collapse but can also enhance it as ex-10 3 10 2 10 1 10 0 10 1 10 2 momentum k in h/Mpc emplified by the graphs for m = 10 −21 eV and m = 10 −22 eV at scales k around or below the Jeans scale for the respective masses.

WEAK LENSING PREDICTIONS
The observational quantity of interest in a weak lensing survey is the shear γ or equivalently the convergence κ.The weak lensing convergence is given by a line-of-sight integration over the density contrast: where the weight function W κ (χ) in the integral is also called weak lensing efficiency and can be modelled as G(χ) is the weighted distance distribution of the lensed galaxies: where q(z) is the galaxy-redshift distribution measured in a weak lensing survey.We assume a simple model for the galaxy redshift distribution used in the Euclid survey (Laureijs et al. 2011) with a median redshift z 0 = 0.9 and β = 1.47.If we know the correlation functions of the convergence, this will give us a way to infer correlation functions of the density contrast.In practice, weak lensing surveys do not measure the convergence κ, but the shear γ by fitting models to or measuring the quadrupole moments of the surface brightness of distant galaxies.The convergence κ is harder to measure since we do not know the intrinsic luminosity of the background galaxies.However, we find |γ| 2 = |κ| 2 in Fourier space and conclude that the respective convergence and shear spectra are equivalent.Schematically, the relationship between the density contrast and the effective convergence is represented in Fig. 8.Our lensing survey parameters are those of a Euclid-like survey.We consider a half-sky survey with f sky = 0.5, a standard deviation of intrinsic ellipticities of galaxies of σ ϵ = 0.4 and an average number of galaxies per steradian n = 4.727 × 10 8 .Perhaps most importantly, we list results up to a maximum multipole moment of ℓ max = 10 4 .Such a high multipole moment is not accessible in a weak lensing survey.In order to see the effect of axion DM in the considered mass range, we need to resolve scales at the order of k = 1 h/Mpc which roughly corresponds to multipole orders ℓ ≳ 10 3 in the case of Euclid.In practice, the highest multipole moment measurable in a weak lensing survey is limited by the shape noise.The signal-to-noise ratio of a single multipole moment drops below 1 at ℓ ∼ 3 × 10 3 .However, we may still gain information by summing over these noise-dominated multipole moments as we make a relatively conservative estimate for the shape noise.At even higher multipole moments of around ℓ ∼ 5 × 10 3 , baryonic feedback becomes a important, a process which we neglect entirely in PT.At this point, our results become completely unreliable.We still show the plots for multipole moments ℓ ≳ 5 × 10 3 in order visually compare the different masses.
More importantly, PT disagrees with N-body simulations on the scales considered and significantly underestimates the full nonlinearity for large k.We therefore compare with an estimate of the lensing spectrum signal obtained by integrating the power spectra obtained from the L = 30 Mpc/h N-body simulations from k = 6 h/Mpc to k = 30 h/Mpc well below the Nyquist frequency k Nyq = π N L = 53.6 h/Mpc.In addition, we demonstrate how the lensing PT results for CDM and the m = 10 −23 eV case change when implementing a hard cutoff on the spectra, that is, we set the input matter spectra to zero above a cutoff scale and study how the PT lensing results depend on the cutoff scale.

lensing spectra, bispectra and trispectra
The angular correlation function for a quantity such as the convergence κ(θ) measured on the sky is given by where the expectation value is computed as ensemble average over statistically equivalent realisations of the field κ(θ).We take the Fourier transform to obtain the angular spectrum in the flat-sky approximation: In order to simplify the computation of angular correlation functions we employ Limber's approximation (Limber 1953).It asserts that if the quantity x(θ) defined in two dimensions is a projection of a quantity y(r) defined in three dimensions with a weight function w x (χ), then the angular spectrum of x is given by where P y (k) is the spectrum of y, evaluated at the three-dimensional wave number k = ℓ/χ.This approximation is applicable if y varies on length scales much smaller than the typical length scale of the weight function w x .Intuitively, we divide χ by ℓ such that we can compare different scales for a given angle.From the Limber approximation, it immediately follows that the convergence spectrum C κ (ℓ) is determined by a weighted line-of-sight integral over the spectrum of the density contrast P δ (k).Likewise, we can express the convergence bi-and trispectrum as appropriately weighted lineof-sight integrals over the bi-and trispectra.All in all, we find where we introduced the subscripts δ to denote the matter spectra as opposed to the convergence spectra denoted by the subscript κ and the weight function W κ is the weak lensing efficiency defined in Eq. ( 37).Fig. 9 shows the respective convergence spectra.Nonlinear corrections significantly increase the magnitude of the dimensionless spectra for multipole moments ℓ ≳ 100.In order to highlight where one would naively expect the models to differ, we translate the comoving quantum Jeans scale into a corresponding multipole order by an order-of-magnitude estimate.We define the quantum Jeans multipole order ℓ J (m) via where χ(z = z 0 ) is the comoving distance at the redshift z 0 amd k J is as defined in Eq. ( 13) at z = 99.For the mean redshift z 0 = 0.9 of the redshift distribution defined in Eq. ( 39 Dimensionless convergence spectra at tree-level, with loop-level corrections and from N-body simulations.The N-body results overlap.The vertical, dotted lines correspond to 0.1 × ℓ J , where the quantum Jeans multipole order ℓ J is defined in Eq. ( 47).
ℓ J m = 10 −23 eV ≈ 5.0×10 3 .Since these quantum Jeans multipole orders are too high to be measurable in a weak lensing survey, the vertical lines in Fig. 9 display 0.1 × ℓ J which roughly describes the multipole order where the CDM PT and FDM PT lensing spectra start to differ.We also compute the lensing spectra using CDM PT with FDM IC via the CDM coupling kernels obtained from the EdS recursion relations as well as the CDM growth factors in the fiducial cosmology.CDM dynamics at tree-and loop-level yield lensing spectra, bispectra and trispectra that are, within the numerical errors, indistinguishable from the ones computed using FDM PT.This is why we refrain from showing the corresponding figures.We conclude that the reason for the suppression of the lensing spectra below the quantum Jeans multipole order in PT lies mainly in the initial conditions and is not the result of the approximation of latetime FDM dynamics through FDM PT.
In contrast, the lensing spectra computed from the N-body spectra in the range k = 6 h/Mpc to k = 30 h/Mpc paint a different picture.The late-time difference in power on small angular scales between the FDM and CDM models is significantly smaller than estimated via PT.Moreover, the FDM models exhibit more power on small angular scales than predicted via PT.Note that the N-body lensing spectra underestimate power for small ℓ because of the cutoff of k.In order to assess, the impact of the high k-modes of the power spectra estimated via PT, we implement a high-k cutoff in the PT lensing results for CDM and m = 10 −23 eV.The tree-level and full loop-level input matter power used as input for the lensing integral are set to zero above the cutoff scale.Fig. 10 demonstrates that a significant part of the discerning power of the model survey comes from the modes with k > 3 h/Mpc.For a smaller cutoff scale, the difference between the CDM and FDM model is not visible and for a cutoff at k = 0.1 h/Mpc where the loop-level PT corrections for CDM are highly accurate, even the difference between the treeand loop-level PT is not discernible.

attainable signal strength
The convergence spectrum is measured by analysing the ellipticities of an ensemble of galaxy images.Averaging over N faint galaxy images, the scatter of the intrinsic ellipticity is reduced to  The vertical, dotted lines correspond to 0.1×ℓ J , where the quantum Jeans multipole order ℓ J is defined in Eq. ( 47).
where σ ϵ is the standard deviation of the intrinsic ellipticity.The angular resolution of this measurement is limited by where n is the average number of source galaxies per squared arc minute.As a result, the observed convergence spectrum C (obs) κ (ℓ) can be modelled as the true spectrum with an additional shot noise contribution: Assuming that the estimates for the spectra can be approximated by a Gaussian distribution, we can estimate the covariance matrices for the lensing spectra, bispectra and trispectra.The covariance of the lensing spectrum is given by where f sky denotes the fraction of the observed sky and we neglect a contribution proportional to the lensing trispectrum due to the non-Gaussianity of the weak lensing field (Kaiser 1998;Scoccimarro et al. 1999).Takada & Jain (2004) provide an expression for the covariance of the weak lensing bispectrum: where ℓ i ≤ ℓ i+1 to count every triangle/quadrilateral configuration only once.∆(ℓ 1 , l 2 , l 3 ) counts the multiplicity of triangle configurations and is defined as otherwise. (53) Similarly, we have for the covariance of the weak lensing trispectrum where ∆(ℓ 1 , ℓ 2 , ℓ 3 , ℓ 4 ) counts the multiplicity of quadrilateral configurations.These covariance matrices allow us to understand the statistical uncertainties on the spectrum measurement.With their help, we can calculate the expected cumulative signal-to-noise ratio Σ(ℓ) for weak lensing measurements of the different spectra up to multipole order ℓ4 :  .Attainable signal to noise-ratio for convergence spectrum at treeand loop-level, convergence bispectrum at tree-and loop-level and convergence trispectrum at tree-level (CDM and FDM identical).
Fig. 13 shows the signal-to-noise ratios obtained in CDM according to Eqs. ( 55), ( 56) and ( 57).We compute the respective covariance matrices using the convergence spectrum with loop-level corrections.This is because the nonvanishing bi-and trispectrum themselves are generated by nonlinear dynamics.Using the treelevel convergence spectrum would therefore underestimate the covariance and overestimate the attainable signal-to-noise ratio.Since the bulk of the cumulative signal comes from the modes with low ℓ, there are no significant differences for the attainable signal-to-noise ratios in CDM and FDM weak lensing surveys for the considered masses.This is also why we do not show the respective N-body results: They would severely underestimate the attainable S2N ratios because of the lack of large-scale power.The sums in Eqs. ( 55), ( 56) and (57) were expressed as integrals and integrated using the CUBA-library (Hahn 2005).We could not compute the respective signal-to-noise ratios for the weak lensing bispectra at loop-level for FDM since the integrals involved proved computationally intractable.
Nevertheless, we compared against the signal-to-noise ratios of the loop-level lensing bispectra computed with CDM PT for FDM IC.Yet, these results are also also subject to substantial numerical uncertainty since the Monte Carlo-integration routine fails to provide error estimates.Fig. 14 visualises the angular dependence of the lensing bispectrum at tree-level for m = 10 −23 eV.The bottom plots reflect that the small-angular scales where the CDM and FDM models actually differ only have a comparatively small signal-to-noise ratio in a weak lensing survey.In contrast, Fig. 15 shows the corresponding loop-level results approximated by CDM with FDM IC.We observe that the loop-level corrections significantly enhance the signal-to-noise ratio for multipole orders where CDM and FDM at m = 10 −23 eV can be distinguished.

differentiability between fuzzy dark matter and standard cold dark matter
We can also give an estimate of whether we can distinguish CDM and FDM experimentally via a weak lensing survey.We assume that the true spectra are given by the CDM spectra and compute the χ 2 -functionals for measuring the noise-weighted mismatch be- tween the true CDM and the wrongly assumed FDM spectra: Fig. 16 shows the χ 2 -functionals for distinguishing CDM and FDM computed according to Eqs. ( 58), ( 59) where all sums are again expressed as integrals and calculated using the CUBA-library (Hahn 2005).
We observe that the χ 2 -values obtained from the N-body simulation in the L = 30 Mpc/h box from k-modes between 6 h/Mpc and 30 h/Mpc are significantly lower than the PT predictions at loop-and even at tree-level.While we expected linear PT to provide a conservative estimate of the distinguishing power of the model lensing survey, it seems that the restoration of small scale power through late-time nonlinear CDM dynamics might cause the real signal to be smaller predicted by PT.Further contrasting FDM PT and CDM PT in Fig. 17 Left column: Dimensionless lensing bispectrum (ℓ 1 ℓ 2 ℓ 3 ) .Differentiability between cold dark matter and fuzzy dark matter in terms of the χ 2 -functional for weak lensing spectra and bispectra as a function of the maximum multipole order ℓ according to Eqs. ( 58) and ( 59).The vertical, dotted lines correspond to 0.1 × ℓ J , where the quantum Jeans multipole order ℓ J is defined in Eq. ( 47).The horizontal blue line corresponds to χ 2 = 1. .Differentiability between cold dark matter and fuzzy dark matter in terms of the χ 2 -functional for weak lensing spectra and bispectra as a function of the maximum multipole order ℓ according to Eqs. ( 58), ( 59) and (60).Both FDM PT and FDM dynamics approximated by CDM PT with FDM IC are shown.The vertical, dotted lines correspond to 0.1 × ℓ J , where the quantum Jeans multipole order ℓ J is defined in Eq. ( 47).
dynamics between FDM and CDM is negligible for the χ 2 -values computed via PT.The only obvious difference at high ℓ arises for the m = 10 −21 eV loop-level PT prediction.Since the other spectra agree well and the FDM PT loop lensing prediction is computed using splines, we conclude that the respective FDM PT χ 2 functional is dominated by noise up to high ℓ.Fig. 18 displaying the tree-level χ 2 -functionals for FDM PT and CDM PT with FDM IC underscores this conclusion.

SUMMARY AND DISCUSSION
In this work, we studied structure formation in the cold dark matter and fuzzy dark matter models and their possible distinction through a Euclid-like weak lensing survey.We extended Eulerian perturbation theory to account for genuine quantum mechanical effects on the de Broglie scale of the dark matter particle.For sufficiently light elementary particles, such as axions and axion-like particles, this scale can be set to be relevant for cosmological structures on the scale of galaxies and below, typically for masses in the range of 10 −23...−21 eV.As a consequence of the Madelung transform of the Schödinger equation, the fluid mechanical equations acquire a quantum-pressure term, counteracting structure formation on scales smaller than the de Broglie scale.We draw the following conclusions: (i) Evolving the density and velocity fields using Eulerian perturbation theory including a quantum pressure term with suitable Gaussian initial conditions leads us to perturbative expressions for the bi-and trispectra, as well as nonlinear corrections to the spectra themselves.We can consistently compute corrections to the spectra and the bispectra at loop-level using FDM perturbation theory and CDM perturbation theory with FDM IC for comparison.Limitations in the numerical evaluation of the resulting integrals via adaptive integration schemes restrict us to tree-level evaluation of the trispectra.In all computations, nonlinear structure formation as predicted by PT does not make up for the lack of initial power on small scales whereas the N-body prediction restores small-scale power on mildly nonlinear scales almost completely.In both CDM .Differentiability between cold dark matter and fuzzy dark matter in terms of the χ 2 -functional for weak lensing spectra, bispectra and trispectra as a function of the maximum multipole order ℓ according to Eqs. ( 58), ( 59) and (60).Both FDM PT and FDM dynamics approximated by CDM PT with FDM IC are shown.The vertical, dotted lines correspond to 0.1 × ℓ J , where the quantum Jeans multipole order ℓ J is defined in Eq. (47).and FDM, non-linear effects leads to increased spectral amplitudes on scales larger than the de Broglie scale.
(ii) Limber projection in the flat-sky approximation with lensing efficiency functions incorporating a Euclid-like source redshift distribution yields lensing spectra, bispectra and trispectra.Again, the lensing quantities are evaluated using adaptive integration.We choose lensing as an observational channel in order to be independent of any biasing assumption typical for galaxy surveys.The respective shear spectra can be found through a correlation analysis on galaxy shapes.As expected, we find a loss in power on the de Broglie scale in the FDM models compared to the CDM model.However, the loss of power predicted by the N-body simulations is significantly smaller than the loss of power predicted by FDM PT.For reference, we convert the de Broglie scale into an angular scale with a typical comoving distance corresponding to the median redshift.Measurements of bi-and even trispectra are clearly within reach of Euclid, with a significance of a few hundred σ: In these estimates we use a full configuration and scale integration with shape noise and a nonlinear covariance.
(iii) With a similar numerical computation, we can estimate whether the spectra, bispectra and trispectra for the CDM and FDM cases are distinguishable: For this purpose, we compute the χ 2functional for a given Gaussian-approximated nonlinear CDM covariance as a function of the FDM particle mass.As in the case of the signal to noise-computations, the χ 2 -functionals follow from a complete configuration space integration up to a limiting multipole, and we consider ∆χ 2 = 1 as a rough criterion of distinguishability: For FDM at tree-level with m = 10 −23 eV this limit is reached at ℓ ∼ 700.The N-body lensing prediction indicates that this limit is reached at ℓ ∼ 1000 − 2000.Yet, for m = 10 −23 eV the latetime impact of the quantum pressure on the scales considered is non-negligible which is why this result suffers from a substantial uncertainty.For masses m = 10 −22 eV and m = 10 −21 eV, the respective signals in our weak lensing survey are too weak to distinguish CDM and FDM below ℓ = 3 × 10 3 ; this is the maximum multipole order accessible in a Euclid-like survey.The lensing bispectrum gives a lower χ 2 -functional than the lensing spectrum at tree-level, but might still serve as an important cross-validation.
(iv) CDM loop-level perturbative results are found to agree reasonably well with the CDM N-body prediction while FDM looplevel PT significantly underestimates the power on scales k ∼ 1 − 100 h/Mpc when compared to the the N-body prediction.Unless, the additional suppression of small-scale power through latetime FDM dynamics significantly enhances the χ 2 -functionals, the values χ 2 = 1 at ℓ ∼ 600 and ℓ ∼ 1500 for m = 10 −23 eV and m = 10 −22 eV by the loop-level lensing spectra are found to be unreliable.In any case, the signal for m = 10 −21 eV is still too weak to be measurable in a realistic weak lensing survey.Considering PT at loop-level, the lensing bispectrum gives slightly higher χ 2 -values than the lensing spectrum.Again, these values are likely unreliable and overestimate the potential of a Euclid-like weak lensing survey.
(v) The main uncertainty in our perturbative approach stems from the modelling of nonlinear structure formation via perturbation theory.In general, the extent to which the FDM and CDM models can be compared depends on the times and scales considered.As the analysis of the impact of a high-k cutoff on the looplevel matter power spectra demonstrates, the bulk of the lensing signal comes from modes with k > 3 h/Mpc for all masses considered.At these scales, loop-level PT is known to underestimate the true power spectra at late times.However, as the comparison with the N-body simulations highlights, full nonlinear CDM dynamics make up for the initial lack of small-scale power while loop-level PT does not.As a consequence, PT seems to significantly overestimate the attainable χ 2 -functionals and is therefore not the right tool to distinguish the considered FDM models from the CDM model in our weak lensing survey.In contrast, PT might be better suited for other observational channels, such as neutral hydrogen surveys for probing the cosmic dawn.They are confined to probing the cosmic state before reionisation at z ∼ 6 where PT still provides a better estimate for the nonlinearity at smaller scales.Technical and numerical details of the FDM perturbation theory are described in detail in appendices A to C, along with the representation of the coupling kernels in terms of Feynman diagrams.
(vi) The estimation of the power spectra from N-body simulations is subject to a number of uncertainties: We only run one simulation per box size and initial condition.Ideally, one averages over a statistically representative sample of simulations with different initial conditions and average the resulting spectra to obtain error estimates.The large-scale power suffers estimated from runs in small simulation volumes is underestimated and suffers from high variance because of a lack of samples.Most importantly, a consistent FDM simulation in a sufficiently large simulation box would be required to reliably estimate the late-time power spectrum on smallscales and as a also consequence the lensing spectra for the FDM masses m = 10 −22 eV and m = 10 −23 eV considered here.There is one free index a that describes the density coupling kernel F 2 (k, k 1 , k 2 ) for a = 1 and the velocity coupling kernel G 2 (k, k 1 , k 2 ) for a = 2.The tree has depth 1 and therefore requires a single time integration.The vertex couplings C (n) owe their name to the diagrammatic representation of Eq. (A13) shown in Fig. A1.In the following, we make use of this diagrammatic language to compute the higherorder coupling kernels.Eq. (A13) can be used to obtain the perturbative solution iteratively: We substitute Ψ = Ψ (1) + O(δ 2 , υ 2 , δυ) into the RHS of Eq. (A13) which gives an equation whose solution is Ψ = Ψ (1) + Ψ (2) + third-order terms.Substituting this expression back into the RHS of Eq. (A13) gives an equation whose solution is Ψ = Ψ (1) + Ψ (2) + Ψ (3) + fourth-order terms and so on.The following section gives expressions for Ψ a (k, τ) up to quartic order in order to compute the kernels F 2 , F 3 and F 4 .The linear velocity divergence fluctuations are related to the density fluctuations via We therefore define the two-component vector f that relates the linear density and velocity divergence fluctuations at time s with the linear density fluctuations at time τ: With this definition in place, Ψ (2) a (k, τ) can be expressed as Comparing this expression to the definition of F n in Eq. ( 23), we can finally read off the coupling kernels.The coupling kernel F 2 can be expressed diagrammatically as shown in Fig. A2 or explic- itly written as Note that all the kernels derived in this section have yet to be symmetrised w.r.t.exchange of momenta.We now proceed by deriving the kernels at third and fourth order in diagrammatic representation.At third order, the diagrammatic representation of the coupling kernels F 3 and G 3 is given by the two types of diagrams shown in Fig. A3.Finally, at fourth order, we find the diagrams shown in Figs.

A1 FDM for EdS cosmology
We now apply the time-dependent PT framework to the FDM case.First, we review the Madelung equations ( 8) in terms of the conformal time τ Next, we Fourier transform the Madelung equations.The Fourier transform of the quantum pressure term is computed by the Taylor expansion up to eighth order using Mathematica (Wolfram Research 2021).We obtain the same continuity equation as in the CDM case as well as the Euler equation with quantum pressure corrections where we omitted the time dependence of all quantities and integrations over momenta with repeated indices are understood.Comparing Eqs. ( A23) and (A24) to Eq. (A4), we can now determine the mode coupling matrices.Following the convention in (Li et al. 2019) we provide them in terms of the time η ≡ 2 √ a = H 0 τ.The linear mode coupling matrix Ω ab reads where the characteristic FDM scale b(k) is defined as Ω 21 vanishes when k equals the comoving quantum Jeans scale k J .We compute the nonlinear mode coupling matrices Γ a,i 1 ,...,in up to n = 4: At second order, the matrices Γ a,i 1 ,i 2 include the nonlinear contribution from the convection term as well as the second-order mode coupling from the quantum pressure term: and all other components vanish.All higher-order contributions to Γ a,i 1 ...in stem solely from the quantum pressure term and therefore represent self-interactions of the density field: and all other contributions at third and fourth order vanish.Next is the Green's function for the linear growth equation ( 15).Using the analytical growth factors give in Eq. ( 19), we find the Green's function Note that the Green's function is free of divergences because the denominators of the growth functions at the initial time cancel.Using equation (A14), the vertex couplings are given by at second order as well as at third and fourth order where all other vertex couplings vanish.This is because above order two there are only self-interactions of the density field stemming from the Taylor expansion of the quantum pressure term.This simplifies the computation and the symmetrisation of the kernels F FDM 3 and F FDM 4 significantly.A more detailed discussion can be found in appendix B. Since 1bc (k, k 1 , k 2 , s, η) inherits this property.As a consequence F FDM 2 as given by Eq. ( A19) is already symmetric under exchange of k 1 and k 2 .

(B8)
The contribution J 4 corresponds to the right diagram in Fig. A5: (B10) the line-of-sight integrations into higher-dimensional integrals that we integrate using the CUBA-library.This proves advantageous since the line-of-sight integrations smooth oscillations and make the loop integrations more numerically tractable.The signal-tonoise sums and χ 2 -functionals are also computed as integrals using the CUBA-library.This leads to difficulties during the evaluation of the FDM loop-level bispectrum signal-to-noise sums and χ 2 -functionals.The corresponding PT integrals are six-dimensional and the line-of-sight integral to compute the lensing bispectra makes them seven-dimensional.The computation of a loop-level signal-to-noise ratio therefore requires a three-dimensional integral of a seven-dimensional integral which does not compute in our tests.As an alternative, we combine the signal-to-noise integrals and the perturbative integrals to obtain a seventeen-dimensional integral which again proves computationally intractable.In the CDM case, this approach yields an eleven-dimensional integral that can be evaluated numerically.We use this approach to estimate the loop-level bispectrum signal-to-noise ratios and χ 2 -functionals.

APPENDIX D: N-BODY SIMULATIONS
We run a total of 8 simulations with N = 512 3 particles with the box sizes L = 30 Mpc/h and L = 256 Mpc/h.They serve as an estimate of the full nonlinear evolution of the CDM and FDM models.The initial conditions are created from the respective CAMB and axion-CAMB spectra using the initial condition generator MUSIC (Hahn & Abel 2011).The simulations themselves are run from z = 99 to z = 0 with a total of 99 snapshots and a gravitational softening length of a comoving 14 kpc/h and 1.7 kpc/h respectively.We use the tool GenPK to compute the matter power spectra from the snapshots (Bird 2017).The lensing spectra are integrated using a 2D spline fitted to the matter spectra in time and k-space.Fig. D1 shows the mass projections of the simulations used for the lensing predictions.The m = 10 −23 eV plot in the lower right corner highlights another problem that N-body simulations with a suppression of small-scale initial power suffer from: the formation of spurious halos arranged like beads on a string along the filaments, for instance (Schive et al. 2016).These spurious collapsed objects likely also lead to an overestimation of small-scale power at late-times in the N-body runs.
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.Asymptotic behaviour of growing modes represented via the quotient of FDM and CDM growth factors D +,FDM (k, a)/D +,CDM (a) in a FDMdominated EdS universe for m = 10 −23 eV.The graphs show the analytical, divergent expression (19), the renormalised expression using the fifth-order Taylor expansion, the numerical solution obtained by integrating the linear growth equation with initial conditions given by Eq. (20) as well as a

Figure 2 .
Figure 2. Growing modes D + (a) in CDM and D + (k, a) in FDM in fiducial cosmology obtained by numerical integration of Eq. (15) for ℏ = 0 (CDM) and a mass of m = 10 −22 eV at three different scales k (FDM).Growth factors are normalised to D + (a 0 ) = 1 at a 0 = 0.01.

Figure 3 .Figure 4 .Figure 5 .Figure 6 .Figure 7 .
Figure3.Dark matter spectra P δ (k) at tree-and loop-level as well as from N-body simulation at z = 0.The N-body results overlap.The dashed, grey line indicates the transition between the different simulation boxes.The dotted, vertical lines denote the respective Jeans scales from Eq. (13) at z = 99.

Figure 8 .
Figure 8. Relationship between density contrast, gravitational potential, lensing potential and effective convergence.
Figure9.Dimensionless convergence spectra at tree-level, with loop-level corrections and from N-body simulations.The N-body results overlap.The vertical, dotted lines correspond to 0.1 × ℓ J , where the quantum Jeans multipole order ℓ J is defined in Eq. (47).

Figure 10 .Figure 11 .
Figure10.Dimensionless convergence spectra at tree-level and loop-level corrections with small-scale cutoff in the input power spectra.

Figure 13
Figure13.Attainable signal to noise-ratio for convergence spectrum at treeand loop-level, convergence bispectrum at tree-and loop-level and convergence trispectrum at tree-level (CDM and FDM identical).

Figure 15 .
Fig.16shows the χ 2 -functionals for distinguishing CDM and FDM computed according to Eqs. (58), (59) where all sums are again expressed as integrals and calculated using the CUBA-library(Hahn 2005).We observe that the χ 2 -values obtained from the N-body simulation in the L = 30 Mpc/h box from k-modes between 6 h/Mpc and 30 h/Mpc are significantly lower than the PT predictions at loop-and even at tree-level.While we expected linear PT to provide a conservative estimate of the distinguishing power of the model lensing survey, it seems that the restoration of small scale power through late-time nonlinear CDM dynamics might cause the real signal to be smaller predicted by PT.Further contrasting FDM PT and CDM PT in Fig.17, we conclude that the difference in late-time Figure16.Differentiability between cold dark matter and fuzzy dark matter in terms of the χ 2 -functional for weak lensing spectra and bispectra as a function of the maximum multipole order ℓ according to Eqs. (58) and (59).The vertical, dotted lines correspond to 0.1 × ℓ J , where the quantum Jeans multipole order ℓ J is defined in Eq. (47).The horizontal blue line corresponds to χ 2 = 1.
Figure17.Differentiability between cold dark matter and fuzzy dark matter in terms of the χ 2 -functional for weak lensing spectra and bispectra as a function of the maximum multipole order ℓ according to Eqs. (58), (59) and (60).Both FDM PT and FDM dynamics approximated by CDM PT with FDM IC are shown.The vertical, dotted lines correspond to 0.1 × ℓ J , where the quantum Jeans multipole order ℓ J is defined in Eq. (47).
Figure18.Differentiability between cold dark matter and fuzzy dark matter in terms of the χ 2 -functional for weak lensing spectra, bispectra and trispectra as a function of the maximum multipole order ℓ according to Eqs. (58), (59) and (60).Both FDM PT and FDM dynamics approximated by CDM PT with FDM IC are shown.The vertical, dotted lines correspond to 0.1 × ℓ J , where the quantum Jeans multipole order ℓ J is defined in Eq. (47).

Figure A1 .Figure A2 .
Figure A1.Diagrammatic representation of Eq. (A13) in terms of trees.Outgoing momentum arrows express momentum integrations.The Dirac-Delta functions enforce momentum conservation at the vertices and we employ the Einstein sum convention for repeated indices.The depth of the tree, i.e. 1, denotes the number of time integrations.

Figure A3 .Figure A4 .Figure A5 .
Figure A3.Diagrammatic representation of contributions to PT kernels F 3 and G 3 .The left diagram describes two contributions with the permutation b ↔ c.The depth of the tree is 2. Accordingly, it requires two time integrations.

Figure A6 .
Figure A6.Diagrammatic representation of contributions to PT kernels F 4 and G 4 involving both the second-and third-order vertex couplings C (2) and C (3) .There are one additional permutation b ↔ c for the diagram on the left and two additional permutations b ↔ c ↔ d for the diagram on the right.