Issue 
A&A
Volume 641, September 2020



Article Number  A130  
Number of page(s)  37  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202038308  
Published online  22 September 2020 
A hydrodynamical halo model for weaklensing cross correlations
^{1}
Institut de Ciències del Cosmos, Universitat de Barcelona, Martí Franquès 1, 08028 Barcelona, Spain
email: alexander.j.mead@googlemail.com
^{2}
Department of Physics and Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada
^{3}
Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK
^{4}
RuhrUniversity Bochum, Astronomical Institute, German Centre for Cosmological Lensing, Universitätsstr. 150, 44801 Bochum, Germany
^{5}
Astrophysics Research Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool L3 5RF, UK
Received:
30
April
2020
Accepted:
30
June
2020
On the scale of galactic haloes, the distribution of matter in the cosmos is affected by energetic, nongravitational processes, the socalled baryonic feedback. A lack of knowledge about the details of how feedback processes redistribute matter is a source of uncertainty for weaklensing surveys, which accurately probe the clustering of matter in the Universe over a wide range of scales. We developed a cosmologydependent model for the matter distribution that simultaneously accounts for the clustering of dark matter, gas, and stars. We informed our model by comparing it to power spectra measured from the BAHAMAS suite of hydrodynamical simulations. In addition to considering matter power spectra, we also considered spectra involving the electronpressure field, which directly relates to the thermal SunyaevZel’dovich (tSZ) effect. We fitted parameters in our model so that it can simultaneously model both matter and pressure data and such that the distribution of gas as inferred from tSZ has an influence on the matter spectrum predicted by our model. We present two variants, one that matches the feedbackinduced suppression seen in the matter–matter power spectrum at the percent level and a second that matches the matter–matter data to a slightly lesser degree (≃2%). However, the latter is able to simultaneously model the matter–electron pressure spectrum at the ≃15% level. We envisage our models being used to simultaneously learn about cosmological parameters and the strength of baryonic feedback using a combination of tSZ and lensing auto and crosscorrelation data.
Key words: cosmology: theory / largescale structure of Universe
© ESO 2020
1. Introduction
Weak gravitational lensing (reviewed by Kilbinger 2015) is the name given to small correlated distortions in the observed shapes of galaxies that are caused by coherent gravitational light deflections between the galaxies and the observer. These shape distortions are of a small amplitude compared to the intrinsic shapes of galaxies and therefore one must have large samples of galaxies in order to tease out signal from the background “shape noise”. The correlation of shape distortions can be used to make maps of a galaxyredshiftweightedversion of the matter distribution, and these can be used to understand the growth and evolution of matter perturbations in the cosmos. Many forthcoming surveys are being designed to measure this effect very precisely, with the ultimate goal of using structure to learn about the accelerated expansion of the Universe.
Weak lensing on its own is a promising tool to infer the distribution of structure. However, in order to extract accurate cosmologicalparameter constraints from data sets, it is necessary to have theoretical models that are significantly more accurate than the survey error bars. On very large scales, density fluctuations are small and linear perturbation theory describes the distribution and evolution or the matter accurately enough for all conceivable future data sets (Lesgourgues 2011). However, the weaklensing signal is predominantly determined by high amplitude, “nonlinear”, fluctuations on small scales (Jain & Seljak 1997) and these have not successfully been described by any abinitio calculation. The standard scheme for modelling the weaklensing signal is therefore to run Nbody simulations over the range of cosmological scenarios under consideration and then to simply measure the fluctuation spectrum as seen in the simulations. One can then either devise fitting functions for the spectra (e.g. Smith et al. 2003; Takahashi et al. 2012; Mead et al. 2015) or build socalled emulation schemes (e.g. Lawrence et al. 2010, 2017; Agarwal et al. 2012; Knabenhans et al. 2019; Angulo et al. 2020) to interpolate between simulation outputs. The lensing signal can then be calculated by integrating over these power spectra along the lineofsight with weightings to account for the projection and lensing.
Most cosmological simulations only consider the gravitational interaction^{1}, partly because it is simpler and partly because nongravitational effects are not understood as well. Early analytic work by White (2004) demonstrated that baryons cooling in a halo may contract and core out the inner halo profile, and it is now known that nongravitational processes can impact the power spectrum of the matter distribution in a way that significantly biases cosmological parameter constraints if not accounted for (e.g. Semboloni et al. 2011; Mohammed et al. 2014; Schneider et al. 2016; Copeland et al. 2018; Huang et al. 2019). The OWLS suite of hydrodynamical simulations of Schaye et al. (2010) were used by van Daalen et al. (2011) to demonstrate that the matter field can be particularly affected by activegalactic nuclei (AGN) feedback and to show that the main effect that reshapes the matter field is the redistribution of gas that is heated as a result of accretion energy that is generated near the central black hole. This heating is extreme enough that gas can even be expelled from the host haloes. Dismayingly, there is a huge range in the possible amplitude of the effects of feedback (Le Brun et al. 2014; McCarthy et al. 2017; Chisari et al. 2018); if this remains unaddressed, then the ability of the next generation of weaklensing surveys to learn about the accelerated expansion will be severely compromised. The effect of baryonic feedback on matter clustering and on the structure of haloes has also been investigated using other hydrodynamic simulations: Martizzi et al. (2012, 2013) look at the effect of feedback on halo profiles. Velliscig et al. (2015) do the same with the COSMOOWLS simulations as well as Mummery et al. (2017), using the more recent “BAryons and HAloes of MAssive Systems” (BAHAMAS^{2}) simulations. All authors agree that the statistics of the matter field on scales relevant for forthcoming weak lensing surveys can be affected at the tens of percent level.
Several techniques have been put forward to mitigate the impact that the unknown feedback strength may have on cosmologicalparameter constraints from weaklensing surveys: Eifler et al. (2015) advocate determining a set of “principal components” from libraries of hydrodynamical simulations and then marginalising over these in a data analysis. Mead et al. (2015) use a modified version of the halo model with parameters that determine the effects of feedback on the halo profiles via a change in concentration and an overall bloating of the halo. The relative merits of these two approaches have been explored by Huang et al. (2019). Several authors (Rudd et al. 2008; Semboloni et al. 2011, 2013; Fedeli et al. 2012, 2014; Fedeli 2014; Debackere et al. 2020) have developed halo models that specifically include a gas component that can be used to model the effect of feedback on the matter–matter spectrum. Recently, van Daalen et al. (2020) have shown that the differences in the effect of feedback on the power spectrum can be mainly attributed to different gas fractions in haloes of ∼10^{14} h^{−1} M_{⊙}, and showed that the feedback amplitude is close to “universal” when expressed in this variable. Another promising technique is the “baryonification” method (Schneider & Teyssier 2015; Schneider et al. 2019; Aricò et al. 2020) where haloes in gravityonly simulations are manually deformed in postprocessing in such a way that they have profiles appropriate for their gas content. The deformation can be informed either from hydrodynamic simulations or from observational data. A variant of this method is investigated by Dai et al. (2018) who use the actual equations of motion governing gas physics to perform the halo deformation.
A more optimistic position to take is to regard the impact of feedback on observables as a positive (e.g. HarnoisDéraps et al. 2015; Foreman et al. 2016; MacCrann et al. 2017); we may be able to learn about energetic processes within galaxies by analysing the weaklensing data. An analysis of CFHTLenS data by MacCrann et al. (2015) included the effect of feedback via a single parameter and a functional form extracted from hydrodynamical simulations. Joudaki et al. (2017a) performed a similar analysis and used the model of Mead et al. (2015) to investigate feedback. Both authors found a modest preference in the data for the presence of feedback. Recent Kilo Degree Survey (KiDS) analyses (Hildebrandt et al. 2017, 2020; Joudaki et al. 2017b) use the model of Mead et al. (2015) and have similar results for the feedback amplitude. Recent results from the Deep Lens Survey (Yoon et al. 2019) using the same model show a stronger preference for strong feedback. However, the matter–matter spectrum may not be the ideal tool if one is interested in energetic processes in galaxies since it is not probing the gas directly.
To measure the gas distribution more directly one can use the thermal SunyaevZel’dovich (tSZ) effect: the inverse Compton scattering of the relatively cold cosmicmicrowave background (CMB) photons by relatively energetic thermal electrons that predominantly exist in galaxy clusters. This scattering results in an increase in energy of the CMB photons, with a characteristic frequencydependence that distorts the standard blackbody CMB spectrum in the direction of the galaxy cluster on the sky. The amplitude of the effect depends on the product of electron number density and temperature along the line of sight, a quantity with the same units as pressure and that is therefore known as “electron pressure”. As a scattering process, the tSZ amplitude does not decay with distance and therefore traces the largescale structure out to high redshift. The characteristic frequency dependence of the induced deviation of the CMB spectrum from blackbody can be exploited to make high fidelity maps of the electron pressure distribution in the Universe (e.g. Planck Collaboration XXII 2016). Due to the fact that tSZ emanates from hot gas in dense clusters, hydrodynamic simulations have been an essential tool for modelling and understanding the effect. Early simulations include those by Suto et al. (1998) and Bond et al. (2005). Trac et al. (2011) developed a technique to paint electron pressure onto a gravityonly simulation and therefore to quickly make mock tSZ maps and power spectrum templates. It was later realised that AGN feedback could have a large impact on the tSZ signal (Dolag et al. 2009; Arnaud et al. 2010; Schaye et al. 2010; McCarthy et al. 2010, 2011; Battaglia et al. 2010), and therefore that these violent feedback episodes needed to be carefully included in simulations, which in turn lead to work to investigate the magnitude of this effect (e.g. Martizzi et al. 2012; Battaglia et al. 2012a,b; Puchwein & Springel 2013; Velliscig et al. 2014; Le Brun et al. 2014, 2015; McCarthy et al. 2014; Sijacki et al. 2015; Dolag et al. 2016).
The lensing–tSZ cross correlation has been measured by Van Waerbeke et al. (2014), Hill & Spergel (2014) and Hojjati et al. (2017). Observational work is often conducted using a stacking approach, where the tSZ map is stacked around the locations of suspected tSZ emitters. Makiya et al. (2018) does the same but with galaxies rather than haloes. Hill et al. (2018) and Tanimura et al. (2020) look at tSZ around SDSS galaxies, see a strong signal and both find evidence for a “twohalo term” – evidence of correlation that extends to scales much beyond the halo boundary. Tanimura et al. (2019a) and de Graaff et al. (2019) looked at residual tSZ emission between close pairs of galaxies and found evidence of the filamentary structure that would be expected to exist given the cosmic web. Tanimura et al. (2019b) mask cluster galaxies and find evidence for a residual “intracluster” tSZ signal. All evidence that supports tSZ emanating from lowdensity regions of the Universe as well as from dense clusters. Hojjati et al. (2015) used the COSMOOWLS simulations to investigate the lensing–tSZ cross correlation measured by Van Waerbeke et al. (2014), similar work was carried out by Battaglia et al. (2015). An advantage of using the cross correlation between lensing and tSZ is that we can ignore many systematic effects that may plague the individual data sets. However, in some cases systematic effects may be correlated between data sets, which may be the case for thermal dust emission from dusty galaxies that will correlate with the lensing signal and may contaminate tSZ maps during the mapmaking process. Yan et al. (2019) showed that this may affect the cross correlation between CMB lensing and tSZ, more than between lowerredshift galaxy lensing and tSZ, due to the higher redshift contribution to CMB lensing. It may be possible to address issues of contamination from thermal dust using the methods of Addison et al. (2012, 2013). More recently, the cross correlation has been measured and investigated by Baxter et al. (2019) and Omori et al. (2019) in the context of a cross correlation of CMB and galaxy lensing with the Dark Energy Survey, and by Osato et al. (2018, 2020) with the Hyper SuprimeCam survey.
To understand the tSZ signal Refregier & Teyssier (2002) construct a halo model and compare this model to a simulation in order to investigate the tSZ autospectrum, demonstrating good agreement with early simulations. Komatsu & Seljak (2002) construct a similar model and demonstrate how the tSZ spectrum scales with cosmological parameters, particularly demonstrating an extreme sensitivity to σ_{8}, the amplitude of the linear power spectrum. These early works suggested that the tSZ autospectrum could contain a wealth of information about cosmological parameters. Holder et al. (2007) examined the relative impacts of feedback (preheating) and changes in cosmology on various tSZ statistics, including the power spectrum. More recent works attempting to utilise tSZ for cosmology are by Shaw et al. (2010), Bolliet et al. (2018) and Hill et al. (2018). In this paper we do not utilise the tSZ autospectrum as we are concerned that it may contain poorly understood internal systematics (Planck Collaboration XXIV 2016; although see Horowitz & Seljak 2017). Hill & Spergel (2014) measured and investigated CMB lensing–tSZ correlation while Ma et al. (2015) investigated the galaxy lensing–tSZ detection from Van Waerbeke et al. (2014). In both cases the cross correlation was modelled using the language of the halo model, with NFW profiles being taken for the dark matter and “universal pressure profiles” (UPP; Arnaud et al. 2010) for the electron pressure.
The purpose of this paper is to develop a halo model that may be used for a joint analysis of lensing and tSZ. We demand that the model be able to reproduce the lensing–lensing auto correlation function and the lensing–tSZ cross correlation. We also demand that the physics underpinning the model be the same so that the two signals depend on each other in a physical way. This is different from previous work in which one profile is used for the matter and a completely separate profile is used for the electron pressure. In this way we hope to be able to learn about cosmological parameters and feedback using both data sets simultaneously. Unless otherwise stated, all figures in this paper are made with cosmological parameters inspired by the WMAP9 (Hinshaw et al. 2013) data analysis: Ω_{m} = 0.2793, Ω_{b} = 0.0463, Ω_{Λ} = 1 − Ω_{m}, σ_{8} = 0.821, n_{s} = 0.972, h = 0.700. These are the same cosmological parameters that are used in the BAHAMAS hydrodynamical simulations.
This paper is set out as follows: in Sect. 2 we show how the halo model can be used to make predictions for the cross spectrum of any field pair, and how we can improve the model by specifically describing a halo in terms of a colddark matter (CDM), gas and stellar component. In Sect. 3 we list our ingredient choices for our halomodel calculations. In Sect. 4 we detail the BAHAMAS simulations, against which we compare our model. In Sect. 5 we tune our halo model by fitting its parameters to simulation data, which is the main results of this paper. Our conclusions are in Sect. 6. We also provide a number of appendices: Appendix A examines a common pitfall regarding computation of the twohalo term. Appendix B discusses how we measure power spectra and shot noise from the BAHAMAS simulations where we have particles with different masses and different properties. In Appendix C we present a series of plots that show how changing our halomodel parameters affects our power spectra. Finally, Appendix D details how we compute projected power spectra of tracers of our underlying threedimensional fields and we show how angular scales of the projected power spectra receive contributions from the underlying threedimensional scales.
2. Halo model
In this section we review the basic halomodel computation of the cross spectrum between two threedimensional cosmological fields (Seljak 2000; Peacock & Smith 2000; Scoccimarro et al. 2001; Cooray & Sheth 2002). We start with a general calculation that can be used to calculate the cross spectrum for any cosmological field pair as long as the halo profiles of the fields are known. Some examples of fields for which this calculation could be applied are “matter overdensity”, “gas overdensity”, “electron pressure”, or “galaxy overdensity”. We then specialise to the details of the calculation appropriate for the lensing–lensing and tSZ–lensing power spectra that are the main focus of this paper, where lensing is generated by the matter overdensity field and tSZ by the electronpressure field^{3}. The modelling presented here draws heavily on work by Refregier & Teyssier (2002), Shaw et al. (2010), Fedeli (2014), Fedeli et al. (2014), Schneider & Teyssier (2015) and Debackere et al. (2020).
2.1. Cross correlations
In the following we omit time dependence from function arguments to make the notation less cluttered. Consider two 3D cosmological fields, u(r) and v(r). The halo model cross spectrum between the two fields at a fixed redshift can be written as the sum of a two and a onehalo term, respectively
$$\begin{array}{cc}& {P}_{2\mathrm{H},uv}(k)={P}_{\mathrm{lin}}(k){\displaystyle \prod _{i=u,v}\left[{\int}_{0}^{\infty}b(M){W}_{i}(M,k)n(M)\phantom{\rule{0.277778em}{0ex}}\mathrm{d}M\right],}\hfill \end{array}$$(1)
$$\begin{array}{cc}& {P}_{1\mathrm{H},uv}(k)={\displaystyle {\int}_{0}^{\infty}{W}_{u}(M,k){W}_{v}(M,k)n(M)\phantom{\rule{0.277778em}{0ex}}\mathrm{d}M,}\hfill \end{array}$$(2)
where P_{lin}(k) is the linear power spectrum of matter fluctuations, M is the halo mass and b(M) is the linear halo bias with respect to the linear matter density (a dimensionless quantity):
$$\begin{array}{c}\hfill {\delta}_{\mathrm{h}}(M,k\to 0)=b(M){\delta}_{\mathrm{m}}(k\to 0).\end{array}$$(3)
n(M) is the halo mass function: the distribution function for the comoving number density of haloes in a mass range (sometimes denoted dn/dM in the literature; units always per mass interval per volume). Equations (1) and (2) contain (spherical) Fourier transforms of the halo profiles of the fields u(r) and v(r) that are being cross correlated:
$$\begin{array}{c}\hfill {W}_{u}(M,k)={\displaystyle {\int}_{0}^{\infty}4\pi {r}^{2}\frac{sin(kr)}{\mathit{kr}}{\theta}_{u}(M,r)\phantom{\rule{0.277778em}{0ex}}\mathrm{d}r,}\end{array}$$(4)
where θ_{u}(M, r) is the averaged radial profile for the field u(r) in a host halo of mass M. For example, if one were interested in the electronpressure field then θ_{u} would be the pressure profile of free electrons, if one were interested in the matter power spectrum, then θ_{u} would be the halo matter overdensity profile. To use these equations to compute the matter–matter power spectrum then one would set ${\theta}_{u}={\theta}_{v}={\rho}_{\mathrm{m}}(M,r)/\overline{\rho}$ where ρ_{m} is the halo matter density profile and $\overline{\rho}$ is the mean matter density. The division by $\overline{\rho}$ ensures that the end result is the overdensity spectrum and we note that ${W}_{\mathrm{m}}(M,k\to 0)=M/\overline{\rho}$. We work using comoving k, so r in the previous equations is also comoving, as is $\overline{\rho}$, which is therefore constant. The units of the halo profile θ_{u} are the units of the field u. The units of the haloFouriertransform functions W_{u} (Eq. (4)) are those of the field u(r) multiplied by volume. We interchange between Δ^{2}(k) and P(k) power spectrum definitions:
$$\begin{array}{c}\hfill {\mathrm{\Delta}}_{\mathit{uv}}^{2}(k)=4\pi {\left(\frac{k}{2\pi}\right)}^{3}{P}_{\mathit{uv}}(k),\end{array}$$(5)
where ${\mathrm{\Delta}}_{uv}^{2}(k)$ has units of the product of fields u and v, while P_{uv}(k) has this and an additional unit of volume.
The dimensionless mass function, g(ν), normalised such that the integral over all ν gives unity, is related to n(M) via
$$\begin{array}{c}\hfill g(\nu )\phantom{\rule{0.277778em}{0ex}}\mathrm{d}\nu =\frac{M}{\overline{\rho}}n(M)\phantom{\rule{0.277778em}{0ex}}\mathrm{d}M,\end{array}$$(6)
and is written in terms of the mass variable ν = δ_{c}/σ(M): δ_{c} is the critical linear density threshold for halo collapse: δ_{c} ≃ 1.686, which has a weak cosmology dependence. σ(M) is the variance in the linear matter field when filtered on a Lagrangian scale R corresponding to mass M = 4πR^{3}/3:
$$\begin{array}{c}\hfill {\sigma}^{2}(R)={\displaystyle {\int}_{0}^{\infty}{\mathrm{\Delta}}_{\mathrm{lin}}^{2}(k){\left[\frac{3}{{(kR)}^{3}}(sin\mathit{kR}kRcos\mathit{kR})\right]}^{2}\phantom{\rule{0.277778em}{0ex}}\mathrm{d}lnk,}\end{array}$$(7)
where the expression in square brackets is the normalised Fourier transform of a realspace tophat filter.
Note that the adopted halo mass function and bias must satisfy the following properties for matter power spectra to have the correct largescale limit^{4}:
$$\begin{array}{cc}& {\displaystyle {\int}_{0}^{\infty}Mn(M)\phantom{\rule{0.277778em}{0ex}}\mathrm{d}M=\overline{\rho},}\hfill \end{array}$$(8)
$$\begin{array}{cc}& {\displaystyle {\int}_{0}^{\infty}b(M)Mn(M)\phantom{\rule{0.277778em}{0ex}}\mathrm{d}M=\overline{\rho}.}\hfill \end{array}$$(9)
In words: these equations enforce that all matter be in a halo of some mass and that, on average, matter is unbiased with respect to itself.
It is worth examining the approximations that lead to the halomodel Eqs. (1) and (2): it has been assumed that haloes trace the underlying linear matter distribution with a linear halo bias, that halo profiles are perfectly spherical with no substructure and that there is no scatter in profile properties at fixed host halo mass. There is also nothing in the modelling to prevent haloes from overlapping. These approximations will break down, and the errors in the eventual power spectrum that they contribute will vary with the fields that are being considered. Unfortunately remedies for these assumptions, where they exist, make the halo model calculation increasingly complex (e.g. Smith et al. 2007; Giocoli et al. 2010; Valageas & Nishimichi 2011; Smith & Markovic 2011; van den Bosch et al. 2013) and are beyond the scope of the work presented here.
2.2. Halo composition and profiles
In this paper we consider each halo of total mass M to be made of separate CDM, gas and stellar components with (possibly timedependent) mass fraction f_{i}(M) in component i, such that these sum to unity
$$\begin{array}{c}\hfill {f}_{\mathrm{c}}(M)+{f}_{\mathrm{g}}(M)+{f}_{\ast}(M)={f}_{\mathrm{m}}(M)=1.\end{array}$$(10)
The gas is further separated into that bound to the halo and that ejected from the halo by some feedback process,
$$\begin{array}{c}\hfill {f}_{\mathrm{g}}(M)={f}_{\mathrm{bnd}}(M)+{f}_{\mathrm{ejc}}(M).\end{array}$$(11)
We use the terms “bound” and “ejected” to differentiate gas that is within the halo virial radius from that outside it. In our formalism it is not necessary that the ejected gas component to be gravitationally bound to the host halo. Indeed, it could be that the ejected gas is far outside the virial radius. All that is required is that the gas component was “associated” with the halo, in that it was part of the initial overdensity from which the halo formed. This means that M should not literally be interpreted as the halo mass, because some mass that was initially associated with the halo may be located outside the virial radius. M is therefore a label, and we take this label to coincide with the measured halo mass in gravityonly simulations, thus justifying our use of a mass function and halo bias calibrated on such simulations. When haloes that are identified in hydrodynamic simulations are compared to nonhydrodynamic simulations that are run with the same initial conditions the halo masses differ objecttoobject. However, there is almost always an objecttoobject correspondence (e.g. van Daalen et al. 2014). Differing halo masses in hydrodynamical simulations arise predominantly because matter has been moved across the (somewhatarbitrary) halo boundary.
We also separate the stellar mass into that in central and satellite galaxies,
$$\begin{array}{c}\hfill {f}_{\ast}(M)={f}_{\mathrm{cen}}(M)+{f}_{\mathrm{sat}}(M),\end{array}$$(12)
where any black hole mass is included in the stellar component.
Each component of the halo is given a spherical density profile, ρ_{i}(M, r), that describes how the component is distributed. These are normalised such that
$$\begin{array}{c}\hfill {f}_{i}(M)M={\displaystyle {\int}_{0}^{{r}_{\mathrm{v}}}4\pi {r}^{2}{\rho}_{i}(M,r)\phantom{\rule{0.277778em}{0ex}}\mathrm{d}r.}\end{array}$$(13)
Note that we can write ρ_{m} = ρ_{c} + ρ_{g} + ρ_{*} and a similar equation for W_{m}.
In this paper we also compute halomodel spectra under the assumption that all matter is CDM, mainly for comparisons with gravityonly simulations. In this case we ignore the gas and stellar components of the halo and set the CDM fraction to unity.
2.3. Largescale limit for matter fields
It is instructive to consider the k → 0 limit of both the two and the onehalo terms in Eqs. (1) and (2) for the specific case of matter fields. Let us illustrate this for the example of the cross spectrum of field u(r) with field v(r), in this case ${\theta}_{u}(M,r)={\rho}_{u}(M,r)/\overline{\rho}$ and Eqs. (4) and (13) give ${W}_{u}(M,k\to 0)={f}_{u}(M)M/\overline{\rho}$, therefore:
$$\begin{array}{cc}& \frac{{P}_{2\mathrm{H},uv}(k\to 0)}{{P}_{\mathrm{lin}}(k\to 0)}={\displaystyle \prod _{i=u,v}\left[\frac{1}{\overline{\rho}}{\int}_{0}^{\infty}b(M){f}_{i}(M)Mn(M)\phantom{\rule{0.277778em}{0ex}}\mathrm{d}M\right],}\hfill \end{array}$$(14)
$$\begin{array}{cc}& {P}_{1\mathrm{H},uv}(k\to 0)=\frac{1}{{\overline{\rho}}^{2}}{\displaystyle {\int}_{0}^{\infty}{f}_{u}(M){f}_{v}(M){M}^{2}n(M)\phantom{\rule{0.277778em}{0ex}}\mathrm{d}M.}\hfill \end{array}$$(15)
If we adopt the notation
$$\begin{array}{c}\hfill \u27e8x\u27e9=\frac{1}{\overline{\rho}}{\displaystyle {\int}_{0}^{\infty}x(M)Mn(M)\phantom{\rule{0.277778em}{0ex}}\mathrm{d}M,}\end{array}$$(16)
for an average of property x over halo mass then Eq. (14) becomes
$$\begin{array}{c}\hfill \frac{{P}_{2\mathrm{H},uv}(k\to 0)}{{P}_{\mathrm{lin}}(k\to 0)}=\u27e8b{f}_{u}\u27e9\u27e8b{f}_{v}\u27e9.\end{array}$$(17)
We see that the amplitude of the twohalo term is governed by the product of the meanhalobiasweighted abundances for fields u(r) and v(r). If we consider the autospectrum of matter overdensity, then u = v = m, f_{m}(M) = 1 and Eq. (17) tells us that the twohalo term is the linear power spectrum because ⟨b⟩ = 1 (Eq. 9). If f_{u}(M) does not depend on halo mass then the average would reduce to Ω_{u}/Ω_{m}, where Ω_{u} is the cosmological density parameter for species u. More generally f_{u}(M) will depend on halo mass, and the twohalo term on large scales is then the linear power spectrum weighted by the field abundance multiplied by the halo bias of the haloes in which the field is to be found.
If instead we focus on the onehalo term, Eq. (15) becomes
$$\begin{array}{c}\hfill {P}_{1\mathrm{H},uv}(k\to 0)=\frac{\u27e8M{f}_{u}{f}_{v}\u27e9}{\overline{\rho}},\end{array}$$(18)
and we see that the amplitude of the onehalo term is governed by the halomassweighted product of the abundances. For the matter–matter case, u = v = m then f_{m} = 1 and the amplitude of the onehalo term is governed by ⟨M⟩, which is the mean halo mass that a unit of matter is to be found in (not the mean halo mass). For our fiducial cosmology at z = 0 this mass is ∼10^{13.5} h^{−1} M_{⊙}, which gives some indication of the typical halo mass responsible for power in the onehalo term of the matter–matter spectrum.
3. Specific ingredients
The method presented in Sect. 2 has been general, but here we list the specific ingredients used in this paper. We are interested in both the lensing auto correlation function and the cross correlation between tSZ and lensing. This led us to compute the 3D matter–matter and matter–electron pressure power spectra. We are also interested in testing how our model for the different components of a halo compares to those components measured in hydrodynamical simulations, so we shall also compute auto and crossspectra of CDM, gas and stars. In this section we take default values for hydrodynamical parameters, but we later free and fit these in Sect. 5.
3.1. Halo demographics
In our calculations we adopt the mass function (g(ν) from Eq. (6)) of Sheth & Tormen (1999)
$$\begin{array}{c}\hfill g(\nu )\phantom{\rule{0.277778em}{0ex}}\mathrm{d}\nu =A[1+\frac{1}{{(q{\nu}^{2})}^{p}}]{\mathrm{e}}^{q{\nu}^{2}/2}\phantom{\rule{0.277778em}{0ex}}\mathrm{d}\nu ,\end{array}$$(19)
with p = 0.3, q = 0.707 and A ≃ 0.216. Note that this function has no explicit redshift dependence. We use the halo bias derived from Eq. (19) using the peakbackground split formalism (Mo & White 1996; Sheth et al. 2001)
$$\begin{array}{c}\hfill b(\nu )=1\frac{1}{{\delta}_{\mathrm{c}}}[1+\nu \frac{\mathrm{d}}{\mathrm{d}\nu}lng(\nu )],\end{array}$$(20)
which then automatically fulfils the meanbias condition in Eq. (9). Explicitly for the mass function in Eq. (19) the peakbackground split gives
$$\begin{array}{c}\hfill b(\nu )=1+\frac{1}{{\delta}_{\mathrm{c}}}[q{\nu}^{2}1+\frac{2p}{1+{(q{\nu}^{2})}^{p}}]\xb7\end{array}$$(21)
The Sheth & Tormen (1999) relation was fitted to haloes that were identified in Nbody simulations with a cosmologydependent virial overdensity criterion taken from the sphericalcollapse model: sphericaloverdense regions that are a factor of Δ_{v} times denser that the background matter density. We take Δ_{v}(z) from the ΛCDM fitting function of Bryan & Norman (1998)
$$\begin{array}{c}\hfill {\mathrm{\Delta}}_{\mathrm{v}}(z)=\frac{1}{{\mathrm{\Omega}}_{\mathrm{m}}(z)}\{18{\pi}^{2}82[1{\mathrm{\Omega}}_{\mathrm{m}}(z)]39{[1{\mathrm{\Omega}}_{\mathrm{m}}(z)]}^{2}\}.\end{array}$$(22)
Δ_{v} is known as the virialcollapse density^{5}. The cosmology dependence of the linearcollapse density, δ_{c}, calculated according to the sphericalcollapse model, was also taken into account in the mass function of Sheth & Tormen (1999). Even though this is a small numerical change it has a larger impact upon the mass function than one might expect due to the fact that it is exponentiated (Courtin et al. 2011; Mead 2017). We use the ΛCDM fitting formula from Nakamura & Suto (1997):
$$\begin{array}{c}\hfill {\delta}_{\mathrm{c}}(z)=\frac{3}{20}{(12\pi )}^{2/3}\{1+0.0123{log}_{10}{\mathrm{\Omega}}_{\mathrm{m}}(z)\}.\end{array}$$(23)
Note that it is important to use halo structural parameters that are consistent with the way haloes are defined in the mass function. We take haloes defined with the virial definition in Eq. (22) and any conversions between halo mass M and ν use the cosmologydependent δ_{c} from Eq. (23).
3.2. Halo composition
We consider haloes to be made from CDM, gas, and stars. We set the halo CDM fraction to be the constant universal fraction
$$\begin{array}{c}\hfill {f}_{\mathrm{c}}(M)=\frac{{\mathrm{\Omega}}_{\mathrm{c}}}{{\mathrm{\Omega}}_{\mathrm{m}}},\end{array}$$(24)
which we consider appropriate because, while feedback may indirectly redistribute CDM within a halo, it will not be able to eject CDM.
Gas is split into bound and ejected gas, with the fraction of bound gas in haloes taken from Schneider & Teyssier (2015):
$$\begin{array}{c}\hfill {f}_{\mathrm{bnd}}(M)=\frac{{\mathrm{\Omega}}_{\mathrm{b}}}{{\mathrm{\Omega}}_{\mathrm{m}}}\frac{{(M/{M}_{0})}^{\beta}}{1+{(M/{M}_{0})}^{\beta}}\xb7\end{array}$$(25)
This function ensures that the bound gas fraction transitions from the universal baryon fraction in highmass haloes (M ≫ M_{0}) to zero in lower mass systems, with β governing the rate of transition; haloes of mass M_{0} having lost half of their initial baryon content. By default we take M_{0} = 10^{14} h^{−1} M_{⊙} and β = 0.6 as in Schneider & Teyssier (2015) as this in in rough accordance with observational results (e.g. Sun et al. 2009; Vikhlinin et al. 2009; Gonzalez et al. 2013). The ejected gas fraction is then fixed such that it accounts for all remaining baryons that were originally associated with the halo assuming that the perturbation was adiabatic
$$\begin{array}{c}\hfill {f}_{\mathrm{ejc}}(M)=\frac{{\mathrm{\Omega}}_{\mathrm{b}}}{{\mathrm{\Omega}}_{\mathrm{m}}}{f}_{\mathrm{bnd}}(M){f}_{\ast}(M).\end{array}$$(26)
With these definitions, the requirement that the total halo mass is accounted for (Eq. (10)) is automatically satisfied^{6}.
We follow Fedeli (2014) and set the stellar fraction to be
$$\begin{array}{c}\hfill {f}_{\ast}(M)={A}_{\ast}exp[\frac{{log}_{10}^{2}(M/{M}_{\ast})}{2{\sigma}_{\ast}^{2}}],\end{array}$$(27)
a form that expresses the fact that star formation efficiency peaks in haloes of mass M_{*} with halo stellar mass fraction A_{*}, while being suppressed for higher and lower halo masses with a logarithmic width of σ_{*}. By default we take A_{*} = 0.03, σ_{*} = 1.2 and M_{*} = 10^{12.5} h^{−1} M_{⊙}, values motivated by Moster et al. (2013) and Kravtsov et al. (2018). We also impose a limit for the highmass end of f_{*}: f_{*}(M > M_{*})≥A_{*}/3, suggested by observational data on the highhalomass saturation of the relation between stellar mass and halo mass (e.g. Leauthaud et al. 2012). The stellar halo mass fraction is further split into central and satellite stars, with the stellar content of lowmass haloes assumed to be dominated by a single central galaxy while the stellar content of highmass haloes is dominated by satellites. For M < M_{*} we take
$$\begin{array}{c}\hfill {f}_{\mathrm{cen}}(M)={f}_{\ast}(M),\phantom{\rule{1em}{0ex}}{f}_{\mathrm{sat}}(M)=0,\end{array}$$(28)
while for M > M_{*} haloes we have
$$\begin{array}{c}\hfill {f}_{\mathrm{cen}}(M)={f}_{\ast}(M){\left(\frac{M}{{M}_{\ast}}\right)}^{\eta},\end{array}$$(29)
$$\begin{array}{c}\hfill {f}_{\mathrm{sat}}(M)={f}_{\ast}(M)[1{\left(\frac{M}{{M}_{\ast}}\right)}^{\eta}],\end{array}$$(30)
with a default of η = −0.3, taken from Moster et al. (2013), such that the stellar mass content of highmass haloes is dominated by satellites while that of lowmass haloes is dominated by centrals.
In Fig. 1 we show the halo mass fractions of the different components as a function of halo mass: The CDM fraction is constant across halo mass. The bound gas fraction is close to the universal baryon fraction for high mass haloes but depletes below half the universal content for M_{0} < 10^{14} h^{−1} M_{⊙}. Lowmass haloes have more baryon content in stars than in gas. The stellar fraction peaks at M_{*} = 10^{12.5} h^{−1} M_{⊙} with a stellar mass fraction A_{*} = 0.03. Stellar mass in satellite galaxies dominates highermass haloes, while stellar mass in central galaxies dominates the lowermass haloes.
Fig. 1. Halo mass fractions f_{u}(M) as a function of halo mass for the different components in our model at z = 0: The halo CDM fraction is constant at Ω_{c}/Ω_{m} for all halo masses. Star formation efficiency peaks around M_{*} = 10^{12.5} h^{−1} M_{⊙} with a mass fraction A_{*} = 0.03, it saturates at high mass at a value of f_{*} = A_{*}/3. Stars are split into central and satellite, with the halo stellar content being dominated by centrals at low masses and by satellites at high mass. The bound gas fraction is near to the universal baryon fraction for high halo masses, but decreases below half the universal fraction for halo masses lower than M_{0} = 10^{14} h^{−1} M_{⊙}, the remaining baryons that are neither bound gas or stars are considered to be in unbound gas that is not in the halo. The total mass fraction drops from unity at lower masses due to the ejected gas. 
3.3. Halo profiles
In gravityonly simulations it has been known since Navarro et al. (1997; NFW) that gravitational collapse causes dark matter to arrange itself into approximately spherical structures following the “NFW” profile
$$\begin{array}{c}\hfill {\rho}_{\mathrm{c}}(M,r)\propto \frac{1}{r/{r}_{\mathrm{s}}{(1+r/{r}_{\mathrm{s}})}^{2}}\xb7\end{array}$$(31)
For our halomodel calculation, the profile is truncated at the virial radius, defined such that this encloses an average density of Δ_{v} times the mean background density:
$$\begin{array}{c}\hfill M=\frac{4\pi}{3}{r}_{\mathrm{v}}^{3}{\mathrm{\Delta}}_{\mathrm{v}}\overline{\rho},\end{array}$$(32)
where Δ_{v}, the virialcollapse density, is taken from Eq. (22). As we always work in comoving units we take $\overline{\rho}$ to be the (constant) comoving matter density and r_{v} is therefore also comoving. The remaining halo scale radius parameter, r_{s} is then usually specified via the concentration c = r_{v}/r_{s} and we use the concentrationmass relation of Duffy et al. (2008) appropriate for the virial definition of halo radius
$$\begin{array}{c}\hfill c(M)=7.85{\left(\frac{M}{2\phantom{\rule{0.166667em}{0ex}}\times \phantom{\rule{0.166667em}{0ex}}{10}^{12}\phantom{\rule{0.166667em}{0ex}}{h}^{1}\phantom{\rule{0.166667em}{0ex}}{M}_{\odot}}\right)}^{0.081}{(1+z)}^{0.71}.\end{array}$$(33)
It has been shown (e.g. Rudd et al. 2008; Velliscig et al. 2014; Mummery et al. 2017) that feedback physics changes the darkmatter halo profiles seen in hydrodynamic simulations when compared to profiles in gravityonly simulations. However, the profiles are still well fit by the NFW form, but have slightly altered concentration parameters. In a similar way to Mead et al. (2015), we introduce a change in concentration parameter caused by the ejection of gas in a way that modifies the standard concentration multiplicatively,
$$\begin{array}{c}\hfill c(M)\to c(M)[1+{\u03f5}_{1}+({\u03f5}_{2}{\u03f5}_{1})\frac{{f}_{\mathrm{bnd}}(M)}{{\mathrm{\Omega}}_{\mathrm{b}}/{\mathrm{\Omega}}_{\mathrm{m}}}]\xb7\end{array}$$(34)
The concentrations of (lowmass) haloes that have lost all of their gas are multiplied by (1 + ϵ_{1}), while those of (highmass) haloes that have not lost any gas are multiplied by (1 + ϵ_{2}). We take ϵ_{1} = ϵ_{2} = 0 as default, which corresponds to no modification. It may be possible to implement an analytical recipe for how gas changes the dark matter profile, such as that presented in White (2004) or Schneider & Teyssier (2015) but we find that this is not necessary for our study.
Gas that is gravitationally bound to a halo is described by the Komatsu & Seljak (2001; KS) density profile, which we take as
$$\begin{array}{c}\hfill {\rho}_{\mathrm{bnd}}(M,r)\propto {\left[\frac{ln(1+r/{r}_{\mathrm{s}})}{r/{r}_{\mathrm{s}}}\right]}^{1/(\mathrm{\Gamma}1)},\end{array}$$(35)
with r_{s} being identical to that in Eq. (31); Γ is the polytropic index for the gas, which we take to be Γ = 1.17 by default. Decreasing Γ makes the gas profile more centrally concentrated. Equation (35) is taken from Martizzi et al. (2013), and is a slightly simplified version of the original KS profile (see Komatsu & Seljak 2002), which we use because it is more convenient to numerically Fourier transform. Rabold & Teyssier (2017) have shown that the KS model provides an extremely good description of haloes that have not undergone violent feedback and Yan et al. (2020) have shown that it also provides a good model for BAHAMAS haloes if Γ is allowed to vary.
We treat the ejected gas in a manner similar to that in Fedeli (2014) and Debackere et al. (2020), which in turn is similar to the way that warm dark matter is treated in Smith & Markovic (2011): we make the approximation that the ejected gas does not contribute to the onehalo term, but only to the twohalo term^{7}. This gives an additive contribution to the twohalo term which is exactly the linear power spectrum multiplied by ⟨bf_{ejc}⟩ per gas field, where
$$\begin{array}{c}\hfill \u27e8b{f}_{\mathrm{ejc}}\u27e9=\frac{1}{\overline{\rho}}{\displaystyle {\int}_{0}^{\infty}b(M){f}_{\mathrm{ejc}}(M)Mn(M)\phantom{\rule{0.277778em}{0ex}}\mathrm{d}M.}\end{array}$$(36)
This is equivalent to taking the ejected gas to be a biased tracer of the linear matter density field. It remains to be determined if this is a reasonable assumption, but it does have the virtue of ensuring that we have accounted for all of the mass and thus preserving the physical relations in Eqs. (8) and (9).
We take the stellar density profile for the central galaxies to be a delta function located at the halo centre
$$\begin{array}{c}\hfill {\rho}_{\mathrm{cen}}(M,r)\propto {\delta}_{\mathrm{D}}(r),\end{array}$$(37)
while we assume that satellite galaxies ρ_{sat}(M, r) follow the NFW profile, given in Eq. (31).
To determine the constant of proportionality for the density profile of each component (Eqs. (31), (35) and (37)) we use Eq. (13): assuring that the component makes up the correct mass fraction of a given halo.
To compute any spectra involving the electron pressure we also need to know the temperature of ionised electrons in haloes. We assume all gas to be ionised, and for the bound gas we can use the KS profile to determine the gas temperature:
$$\begin{array}{c}\hfill {T}_{\mathrm{g}}(M,r)={T}_{\mathrm{v}}(M)\frac{ln(1+r/{r}_{\mathrm{s}})}{r/{r}_{\mathrm{s}}},\end{array}$$(38)
which amounts to assuming hydrostatic equilibrium. We compute the central temperature T_{v}(M) as the halo virial temperature for the gas:
$$\begin{array}{c}\hfill \frac{3}{2}{k}_{\mathrm{B}}{T}_{\mathrm{v}}(M)=\alpha \frac{GM{m}_{\mathrm{p}}{\mu}_{\mathrm{p}}}{a{r}_{\mathrm{v}}},\end{array}$$(39)
where m_{p} is the proton mass, μ_{p} is the mean gas particle mass divided by the proton mass, and the factor of a converts the comoving r_{v} to a physical radius. The total halo electron pressure P_{e} is then given by the ideal gas law
$$\begin{array}{c}\hfill {P}_{\mathrm{e}}(M,r)=\frac{{\rho}_{\mathrm{bnd}}(M,r)}{{m}_{\mathrm{p}}{\mu}_{\mathrm{e}}}{k}_{\mathrm{B}}{T}_{\mathrm{g}}(M,r),\end{array}$$(40)
where ρ_{bnd} is the halo bound gas density, m_{p} is the proton mass and μ_{e} is mean gas particle mass per electron divided by the proton mass^{8}. The first term in this equation is exactly the electron number density. α in Eq. (39) is a parameter that encapsulates deviations from a simple virial relation, for example, unvirialised or unthermalised gas or turbulence, we take the standard α = 1 by default. In our calculations we have neglected relativistic corrections to the effective pressure that contributes to the tSZ spectrum (Itoh et al. 1998; Nozawa et al. 2000), these corrections may be important to consider in the future (Remazeilles et al. 2019).
The ejected gas is treated as if it is distributed with the background linear perturbations, so that it contributes only to the twohalo term. It is given a constant temperature, T_{w}, which is supposed to be the temperature of the warmhot intergalactic medium (WHIM; Cen & Ostriker 1999), but we caution the reader against taking this name too literally. We take T_{w} = 10^{6.5} K as the default temperature of the WHIM (Van Waerbeke et al. 2014).
3.4. Largescale limit for electron pressure
It is informative to calculate the amplitude scaling of the two and onehalo terms for the case of the electron pressure, to compare to matter. If we take the k → 0 limit for the halo Fourier transforms (Eq. (4)) of matter we get W_{m}(M, k → 0) ∝M, but for the electron pressure that originates from bound halo gas this changes to
$$\begin{array}{c}\hfill {W}_{\mathrm{p}}(M,k\to 0)\propto {M}^{5/3}.\end{array}$$(41)
The extra M^{2/3} for electron pressure comes from the fact that pressure is the product of gas density, which is ∝M, and gas temperature (Eq. (39)), which is ∝M^{2/3}. Using the notation from Sect. 2.3 we can write the largescale limit of the twohalo term as
$$\begin{array}{c}\hfill \frac{{P}_{2\mathrm{H},\mathrm{mp}}(k\to 0)}{{P}_{\mathrm{lin}}(k\to 0)}=A\u27e8b{f}_{\mathrm{bnd}}{M}^{2/3}\u27e9+B{T}_{\mathrm{w}}\u27e8b{f}_{\mathrm{ejc}}\u27e9,\end{array}$$(42)
where the second term originates from the pressure of nonhalo WHIM gas (A and B are constants). For the onehalo term we have
$$\begin{array}{c}\hfill {P}_{1\mathrm{H},\mathrm{mp}}(k\to 0)\propto \u27e8{f}_{\mathrm{bnd}}{M}^{5/3}\u27e9.\end{array}$$(43)
The extra M^{2/3} in these equations compared to those for the matter–matter power spectrum cause the matter–electron pressure cross spectrum to be dominated by relatively more massive haloes. A corollary of this is that the transition between the two and onehalo term occurs at larger scales for the matter–electron pressure spectrum compared to the matter–matter spectrum.
3.5. Example power spectra
In Fig. 2 we show the halomodel predictions the spectra of matter–matter, CDM–CDM, gas–gas, stars–stars and matter–electron pressure. We see that at large scales the shape of all curves are identical, a consequence of all components following the linear perturbation distribution (Eq. (14)). For the CDM–CDM spectrum, the ratio at large scales compared to the matter spectrum is exactly (Ω_{c}/Ω_{m})^{2}. For gas–gas and stars–stars it is roughly (Ω_{g}/Ω_{m})^{2} and (Ω_{*}/Ω_{m})^{2}, but this is not exact because the bias of the haloes in which each field lives also contributes to the twohalo term amplitude (see Eq. (17)). At smaller scales we see that the shapes of the onehalo terms are all quite different, which arises from the different halo profiles of each component. Note that for k ≃ 10 h Mpc^{−1} the stellar density distribution eventually has more power than the gas, despite the lower abundance: a consequence of the gas being smoothly distributed on small scales while stars are tightly clustered in halo centres. The CDM–CDM spectrum is very similar to the total matter spectrum, which makes sense since CDM dominates the matter budget; however, there are shape differences at smaller scales where the different scaledependence of contributions of gas–gas and stars–stars are important. The shape of the matter–electron pressure spectrum looks superficially like that of the gas–gas spectrum, but with the two to onehalo transition taking place at a different scale, due to the different scalings of the two and onehalo terms. Note that this spectrum cannot be directly compared with the others since it has units of pressure whereas the others are dimensionless.
Fig. 2. Halomodel predictions at z = 0 for the matter–matter (top), CDM–CDM, gas–gas, stars–stars and the matter–electron pressure (bottom) power spectra. We also show the twohalo (longdashed) and the onehalo (shortdashed) terms for each spectrum apart from matter–matter. At large scales, the shapes of the twohalo terms are all identical, and are that of the linear power spectrum, at small scales the onehalo terms dominates and have significantly different shapes for each spectrum. Note that the transition scale between the two and onehalo term is at much larger scales for the matter–electron pressure cross spectrum, which is a consequence of this spectrum being dominated by contributions from more massive haloes compared to those that contribute to the matter spectra. The downturn at small scales in the gas–gas and matter–electron pressure spectra is partly due to the underrepresentation of lowmass haloes and partly due to their relatively smooth halo profiles. 
In Fig. 3 we show how the z = 0 matter–matter spectrum and the matter–electron pressure spectrum are built up as a function of halo mass. We show cumulative power spectra computed as we vary the upper limit of halo mass in the integrals in Eqs. (1) and (2) from 10^{10} to 10^{16} h^{−1} M_{⊙}. In both cases the integration has converged with an upper limit of 10^{16} h^{−1} M_{⊙}. The matter–matter spectrum builds up more gradually as a function of halo mass compared to the matter–electron pressure spectrum, the latter receiving the majority of the power from haloes between 10^{14} and 10^{15} h^{−1} M_{⊙} and very little power from haloes below 10^{13} h^{−1} M_{⊙}. This, combined with the relatively smooth pressure profiles, contributes to the existence of a peak in the model cross spectrum at k ∼ 7 h Mpc^{−1}, although we caution the reader that this peak may not be physical since we do not model halo substructure. By extension, this implies that the matter–electron pressure cross spectrum is sensitive to higher halo masses than the standard matter–matter spectrum (e.g. Makiya et al. 2018). The reason for this is two fold: first, Eq. (25) tells us that lowmass haloes are deficient in gas and therefore contribute less freeelectron number density to the signal. Second, Eqs. (42) and (43) tell us that the overall amplitude of both the two and the onehalo terms is more sensitive to highmass haloes. This also ensures that the nonlinear scaling of these spectra with σ_{8} is more extreme, as can be seen in Fig. 4, since increasing the power spectrum amplitude a small amount has a comparatively large impact on the highmass tail of the halomass function. We find that roughly ${P}_{\text{mm}}\propto {\sigma}_{8}^{3.5}$, ${P}_{\text{mp}}\propto {\sigma}_{8}^{5.8}$ and ${P}_{\text{pp}}\propto {\sigma}_{8}^{8}$ with the exact exponent depending on wavenumber. This fact has been noticed in the past (e.g. Komatsu & Seljak 2002; Hill & Pajer 2013; Hill & Spergel 2014) and motivates the use of the tSZ autospectrum (Komatsu & Seljak 2002; McCarthy et al. 2014), or the cross correlation of tSZ with lensing (Hill & Spergel 2014; Ma et al. 2015; Hojjati et al. 2017), as a sensitive probe of the powerspectrum amplitude. This statement is equivalent to noting that the highmass tail of the halo mass function is a sensitive probe of the power spectrum amplitude.
Fig. 3. Cumulative halo model power spectra as a function of the upperlimit of halo mass at z = 0. Lefthand panel: matter–matter power spectrum as the upper limit is raised from 10^{10} to 10^{16} h^{−1} M_{⊙} (i.e. the lightest coloured curve shows the contribution from only haloes below 10^{10} h^{−1} M_{⊙}), righthand panel: same but for the matter–electron pressure spectrum. In both cases the converged spectrum from all halo masses is shown in thick black and this has been reached with an upper limit of 10^{16} h^{−1} M_{⊙}. We see that the matter–matter spectrum builds up fairly gradually with halo mass and that all scales receive contributions from a wide range of halo masses. In contrast, the matter–electron pressure spectrum receives very little contribution from haloes less massive than 10^{13} h^{−1} M_{⊙}, with the vast majority of the power coming from haloes between 10^{14} and 10^{15} h^{−1} M_{⊙}. 
Fig. 4. Effect of σ_{8} variations on the matter–matter power spectrum (left) and on the matter–electron pressure power spectrum (right) at z = 0. We show halo model power spectra for σ_{8} = 0.7 (blue) to σ_{8} = 0.9 (red) divided by a central σ_{8} = 0.8 model. For linear theory the plot would simply be horizontal lines. We see that boosting the amplitude of fluctuations increases the amplitude of the matter–electron power spectrum far more than it does for the mattermatter. At a wavenumber of ≃1 h Mpc^{−1} the mattermatter power spectrum scales like ${\sigma}_{8}^{3.5}$ while matter–electron pressure scales like ${\sigma}_{8}^{5.8}$. This is because the electron pressure field is dominated by highmass haloes from the tail of the mass function that are very sensitive to the linear power spectrum amplitude. 
4. Hydrodynamic simulations
In this paper we compare auto and crossspectra from our halo model to those measured from the BAHAMAS simulations (McCarthy et al. 2017). These are a set of smoothparticle hydrodynamical simulations with separate dark matter and gas particles. The dark matter interacts only gravitationally but the gas also experiences a hydrodynamic pressure force. In addition, processes like star formation, gas heating and cooling, supernovae explosions and AGN formation and evolution are also modelled using “subgrid” recipes. A full discussion of the subgrid implementation can be found in Schaye et al. (2010), Le Brun et al. (2014), and McCarthy et al. (2017). Hydrodynamic gas particles have temperature and density as properties, as well as the standard position and velocity, and these can be used to calculate an electron pressure per particle (Appendix B). Star (macro) particles are also created from the gas as the simulation evolves, and these keep track of star formation and death. Black holes are also considered, which can swallow gas particles and grow, although they contribute little to the overall mass.
The BAHAMAS simulation boxes are 400 h^{−1} Mpc cubes with periodic boundary conditions. Each simulation is of 1024^{3} dark matter and 1024^{3} gas particles with cosmological parameters inspired by the WMAP9 (Hinshaw et al. 2013) data analysis. Initial conditions are created at z = 127 using secondorder Lagrangian perturbation theory via the code SGENIC^{9} (Bird et al. 2020) with appropriate, different transfer functions used for dark matter and baryons. We utilise three different hydrodynamic simulations that differ only in the “strength” of their AGN feedback. This strength is determined by a single parameter, the AGN subgrid heating temperature, which is the temperature increase given to gas particles that are targeted for feedback. The default AGN model, for which we have 3 realisations, was calibrated to reproduce the amplitude of the observed gas–halo mass relation of groups and clusters, as inferred from Xray observations (McCarthy et al. 2018). The AGNLO and AGNHI models lowered and raised the subgrid heating temperature so as to approximately bracket the scatter in the observed relation. For comparison purposes there is also a nonhydrodynamic “gravityonly” simulation (DMONLY). This simulation still uses two sets of particles for gas and dark matter, and these different particles start with different initial conditions. The only difference compared to the full hydrodynamic simulations is in the subsequent evolution, where only gravity is considered in DMONLY.
In Fig. 5 we show different fields around a massive halo for the AGN simulation at z = 0. We see that the CDM forms a skeleton around which the other fields are clustered, with the gas density being less tightly clustered than CDM and stellar density being more tightly clustered. Both gas and stars are missing from the lowermass haloes. The electron pressure follows the gas density, but is restricted to emanate from only the highest gasdensity peaks. This is because highmass haloes also have a higher temperature compared to lowmass haloes, which boosts their electronpressure signal relative to the density. Lowmass haloes are therefore severely under represented in the electronpressure distribution.
Fig. 5. Different fields measured from the BAHAMAS AGN simulation at z = 0 around a M ∼ 10^{15} h^{−1} M_{⊙} halo averaged through a square slab of side 20 h^{−1} Mpc and thickness 5 h^{−1} Mpc. We show CDM (top left), gas (top right) and star (bottom left) density contrast as well as electron pressure (bottom right). Colour intensity increases with field value logarithmically and we show a dynamic range of 10^{3} for both density contrast and electron pressure. We see that the gas is less clumped than the CDM and that lowmass sub haloes show a large gas depletion compared to the host. Stars are more tightly clustered in the halo centre than the CDM, but are absent in many of the substructures. The total matter densitycontrast field would be a sum of the CDM, gas and star fields. The electron pressure broadly follows the gas but can only be seen emanating from the highest gasdensity peaks over the dynamic range shown. 
In this paper, we are interested in twopoint statistics that pertain to weak gravitational lensing and to the tSZ effect. We therefore measure the power spectrum of density fluctuations and the electronpressure field from the simulations. We work with the full 3D data, rather than with 2D projections, as these are more directly tied to the modelling. Details of how we measure these power spectra from the particles and how we consider the effects of shot noise can be found in Appendix B. We compute the auto and crossspectra of all combinations of the fields: total matter overdensity δ_{m}, CDM overdensity, δ_{c}, gas overdensity, δ_{g}, stellar overdensity, δ_{*} and electron pressure P_{e}. In our work, all overdensities are defined relative to the total matter density (i.e. $1+{\delta}_{u}={\rho}_{u}/\overline{\rho}$), which ensures
$$\begin{array}{c}\hfill 1+{\delta}_{\mathrm{m}}=(1+{\delta}_{\mathrm{c}})+(1+{\delta}_{\mathrm{g}})+(1+{\delta}_{\ast}).\end{array}$$(44)
The measured power spectra from the BAHAMAS AGN simulation are shown in Fig. 6 together with the halomodel prediction from the model presented in Sect. 3 with the default parameter values. We note that the halo model provides a reasonable model of the data for all spectra shown, but that it is not perfect (the loglog scale can hide some serious defects). At large scales we note that the amplitudes are generally in good agreement, which tells us that we have the overall abundances and halo occupation of each species approximately correct and that we have the mean background electron pressure reasonably well modelled. At smaller scales we note differences that must be due to an incorrect mass function, choice of halo profiles or else due to physical effects that are missing from our simple halo model. For the matter spectra we also note that the standard problem of a power deficit in the transition region between the two and onehalo terms (e.g. Tinker et al. 2005; Valageas & Nishimichi 2011; Mead et al. 2015) is present for all spectra shown. This defect seems to be less of a problem in the spectra involving the electron pressure. With reference to Fig. 2 we can conjecture that this may be because these spectra are more dominated by the onehalo term and because the transition between the two and onehalo terms takes place at comparatively larger scales. The existence of a peak in both the matter–electron pressure cross spectrum and electron pressure auto spectrum is hinted at, but by no means confirmed, by the simulation measurements.
Fig. 6. Comparison of the default halomodel prediction with power spectra from the BAHAMAS AGN simulation at z = 0. Points with errors show the measured power spectra with an erroronthemean calculated from the finite number of modes that contribute to each k. The upper two panels show spectra for “all matter” (purple), CDM (green), gas (light blue), stars (orange) and electron pressure (dark blue). The electron pressure field has units of 100 eV cm^{−3}. The effect of power aliasing can be seen as an upturn in power at the highest k (∼7 h Mpc^{−1}) shown for each spectrum (see Appendix B). Lines show the halomodel predictions using the default model discussed in Sect. 3: solid lines and points show autospectra while broken lines and open points show crossspectra with the total matter field. We see that the trends in the halo model agree reasonably well with the simulations but that there are disagreements in the details. In the lower two panels we show the “response function” calculated with respect to the matter–matter power in the gravityonly model: ${P}_{\mathit{uv}}^{\mathrm{hydro}}(k)/{P}_{\mathrm{mm}}^{\mathrm{gravity}}(k)$. The simulations spectra have been divided by those from the DMONLY simulation and the halo model has been divided by a standard halomodel prediction that assumes all matter to be in NFW haloes. The horizontaldashedblack lines in the lower left panel show the expected large scale Ω_{c}/Ω_{m}, (Ω_{c}/Ω_{m})^{2}, Ω_{b}/Ω_{m}, and (Ω_{b}/Ω_{m})^{2} asymptotes for the CDM and gas cross–matter and autospectra. We see that the response functions are generally smoother than the raw power spectra: a consequence of cosmic variance cancellations at large scales and cancellation of aliasing effects at small scales. 
The lower panels of Fig. 6 show the power spectrum response functions (e.g. Mead 2017; Cataneo et al. 2019), which we define as the ratio of any power spectrum of any field combination to the matter–matter power spectrum measured in a “gravityonly” model. In the case of the BAHAMAS simulations we calculate this response with respect to the DMONLY simulation. The response has the virtue that it cancels the Gaussian cosmic variance at large scales, which leads to a smooth largescale response. The measured response functions for the BAHAMAS spectra that involve the electron pressure are less smooth at large scales, which is probably because of the comparative dominance of the onehalo term at large scales and therefore that the largescale noise is likely to be dominated by Poisson fluctuations in massive halo numbers, rather than the Gaussian variance present in the initial conditions. For the halo model, we calculate this response ratio with respect to a halo model calculation assuming all mass to be in NFW haloes (f_{c} = 1). That the response functions are constant at large scales is indicative that the power on these scales is simply the linear matter spectrum multiplied by some weighting (Eqs. (17) and (42)). For the halo model, the response has the virtue of cancelling some of the standard problems such as the lack of power in the transition region (e.g. Tinker et al. 2005; Mead et al. 2015), as can be seen in Fig. 6, and helps to a lesser extent with the general underestimation of power at smaller scales (Giocoli et al. 2010). Working at the level of the response may also alleviate problems that may arise from not using the most uptodate ingredients for our halo model. For example, we use the mass function of Sheth & Tormen (1999) whereas there exist more recent mass functions such as Tinker et al. (2008, 2010) or Despali et al. (2016). We prefer Sheth & Tormen (1999) because it was calibrated to a wider range of cosmological parameter space than more modern mass functions.
5. Fitting halo model parameters
The level of agreement between the simulations and our model shown in Fig. 6 in encouraging, but demonstrates that the model is not sufficiently accurate to use to draw robust conclusions from forthcoming survey data. Note that the signaltonoise is roughly 3:1 in the lensing–tSZ C(ℓ) measurements of Hojjati et al. (2015) and is expected to be 5:1 for forthcoming KiDS measurements (Tröster et al., in prep.). It is possible that some of this modelling inaccuracy arises from incorrect choices for ingredients, but we think that a substantial amount must also arise from missing features in our basic halomodel calculation. Therefore, following the logic in Mead et al. (2015) we improve our model by fitting parameters that pertain to gas physics to the BAHAMAS simulations at the level of the power spectra. The aim is to find parameters that govern a population of “effective” haloes that describe the power spectra well for the different feedback strengths. We do this at the level of the halomodel “response”, discussed in Sect. 4. Once an accurate model for the response has been developed, we must then multiply this response by an accurate prescription for the matter–matter power spectrum in the gravityonly scenario. In our case we take the accurate prediction to be from HMCODE (Mead et al. 2015, 2016), but we could have used HALOFIT (Smith et al. 2003; Takahashi et al. 2012) or an emulator prediction (e.g. COSMIC EMU; Lawrence et al. 2010, 2017). Note that our model does not generally require us to fit parameters in this way, its utility is more general, fitting just seems to us to be the most obvious way for us to proceed to construct an accurate model.
We present several different models, the utility of each depending on the eventual use case: (1) stars, (2) matter, (3) matter and electron pressure and (4) matter, CDM, gas, and stars. Because the stars are a subdominant component of the total matter budget, we find it necessary to fit parameters that govern the stars separately (1), and then to hold these parameters fixed while we fitted models to the matter (2) and to the matter and electron pressure (3). In case (4), we fitted a model to all the matter fields simultaneously, refitting parameters that govern the star distribution. The free parameters that we have fitted are listed in Table 1 together with descriptions of their physical meaning and their default values. These parameters were chosen by trialanderror to form minimal set that were able to model the data well, without opening the parameter space too widely.
Default hydrodynamical halo model parameters.
To fit we use the Nelder & Mead (1965) simplex algorithm with a large number of initial starting locations so as to avoid local minima. To determine the redshift dependence and relevant parameters, we initially fitted our model to different redshifts separately and we then used this information to parameterise sensible redshift dependencies for some of the parameters. The parameters are fitted to z < 1 BAHAMAS data across redshifts simultaneously with a linear weighting in z and a log weighting in k between 0.015 and 7 h Mpc^{−1}. Our model currently takes ∼0.5 s per z to evaluate all 15 cross spectra we are interested in. The numerical bottleneck is the computation of the Fourier transforms of the gas density and electron pressure profiles, which must be done numerically.
Before presenting our results we caution the reader against taking the details of our modelling as a serious physical description of the underlying haloes. We are fitting parameters in a halo model to power spectra taken from simulations; we are not fitting the simulations at the level of individual halo profiles. Therefore, the parameters of our hydrodynamical halo model should be thought of as pertaining to “effective” haloes that, through the apparatus of the halo model, provide an accurate model of the spectra we are interested in. They should not be confused with the exact parameters that govern the actual internal structure of actual haloes. Although they may relate to these, it would be necessary to check explicitly that this was so. We remind the reader of the shortcomings of the halo model approach: there may be features in the power spectra of the fields we are interested in that will never be accurately described by a linearly biased population of spherical, virialised haloes. Some of our parameter fitting may account for some of these deficiencies, although working with the response also helps to ameliorate some of these problems. Recall Fig. 3, where we saw that haloes of different masses are moreorless important to different power spectra. This implies that we may drastically alter the profiles of haloes that contribute negligibly to a particular power spectrum and leave the spectrum almost unchanged. It would obviously be incorrect to over interpret the physical meaning of fitted parameters in this instance. Also remember that computing a power spectrum via the halo model represents a huge compression of the information that is potentially available from the halo profiles individually, which is smeared out across k by integrations. It is perfectly possible for inaccuracies in the modelling of the halo population to add or cancel in various unintuitive ways during these numerical compressions. Huffenberger & Seljak (2003) demonstrated that halostructural parameters derived directly from haloes extracted from a cosmological simulation were not the same as those derived by fitting the structural parameters in a halomodel power spectrum calculation to the measured power spectrum from the same simulation. Mead et al. (2015) demonstrated that relatively drastic nonphysical additions were required in order to improve halo model predictions for the matter–matter power spectrum from the 30% level to the 5% level, and that these were particularly important in the transition region. Tinker et al. (2005) and van den Bosch et al. (2013) both demonstrated that problems exist when trying to model galaxy clustering and galaxy–galaxy lensing using the halo model. To our knowledge, it has never been demonstrated that power spectra from simplistic halomodel calculations agree with those measured from simulations, even if the ingredients for the halo model calculation are taken directly from the simulation of interest.
5.1. Stars
The first model (1) we present is for the star–star power spectrum in Fig. 7. The free parameters are A_{*}, M_{*} and η with their bestfitting values listed in Table 2. We find that we are able to fit a reasonable model to the different simulations with different AGN feedback strengths (temperatures), and that our model follows the data with an accuracy of ≃5% for the response, nicely demarcating the different feedback scenarios. We note a preference for an inverse correlation between both A_{*} and M_{*} and the feedback temperature. This makes physical sense, since increased AGN activity is suppressing star formation. This could result in both an overall suppression in the number of stars in a halo of a given mass and also in pushing the peak starformation efficiency to lower halo masses. In contrast, the bestfitting η, which governs the split between central and satellite stellar mass, is very similar for each of the different feedback temperatures. Note that there is a scatter in the efficiency of this suppression of star formation by AGN between simulations from different groups. Indeed a correct stellar mass–halo mass relation is one of the targets that hydro simulations strive for, and in practice this can be quite difficult to achieve. We also looked at additionally fitting the parameter σ_{*}, but found that this did not produce a significant improvement in the quality of the fit.
Fig. 7. Halo model (1) for the star–star power spectrum is shown as sold lines while simulation measurements are shown as points with error bars that come from the scatter in power between three different realisations of the AGN BAHAMAS simulation. Note that the power spectra shown here is the measured response multiplied by HMCODE, so cosmic variance is removed at small scales, but we still show the error bar for comparison. Different colours denote the three different AGN feedback temperatures: 10^{7.6} (blue), 10^{7.8} (grey) and 10^{8.0} K (red). Parameters for the model can be found in Table 2. We see that the spectrum is reduced in amplitude as feedback temperature increases, a consequence of AGN feedback suppressing star formation. This is reflected in the modelling as the value of A_{*} and M_{*} decreases as AGN strength increases. 
Bestfitting effective halo model parameters for different combinations of power spectra from the BAHAMAS simulations.
5.2. Matter
In Fig. 8 we present model (2), in which we fitted the parameters ϵ_{1}, Γ and M_{0} to the matter–matter data only. The underlying star model, which contributes to the matter–matter spectrum, is fixed to that shown in Fig. 7 and discussed in the previous subsection. We are able to match the powerspectrum response seen in the simulations at the percent level for each of the feedback models and at all of the redshifts shown. The fitted model has a strong correlation between M_{0} and AGN strength, which makes physical sense as stronger feedback ejects more halo gas. We also see strong trends in the fitted values of ϵ_{1} and Γ, which may be due to different amounts of back reaction on the dark matter and different heating of the residual halo gas. The higher AGN heating temperatures favours higher Γ, which corresponds to a less concentrated gas profile, as might be expected from more violent feedback. The upturn in the response function at small scales originates from the stellar contribution to the matter spectrum. In addition, in Fig. 8 we also show how the model fares when pushed beyond the range of AGN heating temperatures to which it was fitted. To do this we linearly interpolate and extrapolate the parameters of our halo model as a function of log T_{AGN}/K. Indeed, this is the reason to use a physically motivated halo model in the first place, as opposed to a simple fitting function or a more blind interpolation between simulation results. This demonstrates some nice physical features of our model that may not be respected by more simplistic fitting functions (such as that in Mead et al. 2015). The suppression in the matter–matter power spectrum in our model originates because gas has been expelled from haloes, with the threshold mass increasing with the AGN heating temperature. Then, to a lesser extent, there are some additional effects from the nonNFW profile of the remaining gas and some back reaction on the concentration of the darkmatter halo component and then some effect from the very centrally concentrated stellar distribution which shows as an upturn at small scales ∼10 h Mpc^{−1}. Because of this, there are some limits to how our model can behave: At one extreme we could raise the AGN heating temperature to a very high value that would result in almost all halo gas being expelled. The maximum suppression that this could have on the power would be to lower the amplitude of the onehalo term by (Ω_{c}/Ω_{m})^{2} ≃ 0.7. At the other extreme, a low value of the AGN heating temperature would result in almost no gas expulsion and the matter–matter spectrum would be almost unaffected. These limits can be seen to be approached by the extreme temperature values shown in Fig. 8. We also looked at additionally fitting the parameters ϵ_{2} and β, but found that these did not produce a significant improvement in the quality of the fit.
Fig. 8. Similar to Fig. 7, but for the matter–matter spectrum (2). The power spectra are shown in the top row, while the response functions are shown in the bottom row. We see that the model is nearly perfect (errors at the subpercent level) across a wide range of scales and for the different feedback scenarios. We see that the amount of powerspectrum suppression increases with the feedback temperature. The critical parameter that governs this in our modelling is the halo mass M_{0}, which governs the halo mass below which more than half of the cosmic baryon density is missing from the halo. In the fitting parameters we see a strong preference for higher values of M_{0} for stronger feedback temperatures (from blue to red). How our model fares outside the range of AGN heating temperatures to which it was fitted is also shown via the lines that have no corresponding simulation points. To generate these we have interpolated and extrapolated our model parameters as a function of AGN heating temperature, which is shown in 0.1 dex intervals between 10^{7.5} (bluest) and 10^{8.1} K (reddest), indicating that our model is robust to extrapolation. 
We also investigated the predictions of our model at fixed AGN heating temperature as we vary the underlying cosmology. As shown in van Daalen et al. (2020), the impact of AGN feedback on the matter–matter spectrum is quite insensitive to the difference between the Planck and WMAP9 cosmologies. We find a similar insensitivity in our model, which is reassuring, but we do find a small difference in the sense that the matter–matter power suppression in the Planck cosmology is predicted to be slightly less than that in WMAP9 at fixed AGN strength, which was shown by van Daalen et al. (2020). This presumably originates from the slightly different baryon fractions between the two different cosmologies, which means that slightly more halo mass is lost in the WMAP9 case.
5.3. Matter and electron pressure
In Fig. 9 we show the main result of this paper, model (3), which is a result of a joint fit to the matter–matter and matter–electron pressure spectra. This model could be used to predict both shear–shear and shear–y tomographic correlations. We do not consider the electron pressure autospectrum because while testing we discovered that this autospectrum is particularly difficult to fit; we suspect that this is because it receives much of its contributions from only a few of the most massive haloes present in the BAHAMAS simulation, and therefore that the spectra measured from the simulations will have very high Poisson noise contributions. Indeed, we see large variations in this power spectrum between the three different realisations of the BAHAMAS, and this variation dwarfs the differences between the feedback strengths for k < 2 h Mpc^{−1}. The matter–electron pressure spectrum still suffers from this problem somewhat, as can be inferred from the error bars in the second row of Fig. 9, but the wavenumbers for which the variance between the 400 h^{−1} Mpc boxes is greater than the difference between the feedback models are k < 1 h Mpc^{−1}. We see that we are able to simultaneously match the response functions at the few percent level for the wavenumbers at which the feedback scenarios are clearly demarcated. We see a trend that increasing AGN heating temperature causes an increase in α, the parameter that governs departures from simple virial temperature scaling, suggesting that gas that avoids being ejected is nevertheless heated significantly. However, the general trend is for suppressed smallscale power as AGN strength increases, which in our model arises from the fact that more gas is ejected from the more massive haloes, thus decreasing the pressure overall. In the matter–electron pressure case we note a trend that increasing AGN strength suppresses the power spectrum for scales dominated by the onehalo term; however, particularly at the higher z, we see that the power spectrum is relatively enhanced at larger scales. This could be because the heated gas is pushed out to larger scales by the feedback and this in turn gives an excess largescale pressure contribution, either directly, or via shock heating the intergalactic medium. In Fig. 9 we also show how our model for the matter and electron pressure fares when evaluated for AGN heating temperatures outside the range within which the model was constructed. As in Fig. 8 we see that the model behaves sensibly outside of the standard temperature range and thus we suggest that it can be used for cosmological and astrophysical parameter inference. We also looked at additionally fitting the parameters ϵ_{2} and β, but found that these did not produce a significant improvement in the quality of the fit.
Fig. 9. Same as Figs. 7 and 8, but for model (3) for the matter–matter (first and third row; power and response) and matter–electron pressure spectra (second and fourth row; power and response). The match to the matter–matter spectra is slightly degraded compared to the case when we fitted the matter–matter spectra exclusively, but we still match the response at the few percent level while simultaneously we match the matter–electron pressure spectrum well for k > 1 h Mpc^{−1}. If we compare this to the error bars in the power spectra in the second row we see that our fit is good for scales where the three AGN heating temperature strengths are clearly demarcated in k. At larger scales we suspect the larger errors arise because the power is predominantly coming from few massive haloes. Similar to Fig. 8, we show how the model fares outside the range of AGN strengths to which is was fitted. AGN heating temperature is shown in 0.1 dex intervals between 10^{7.5} K and 10^{8.1} K, while the model was fitted between 10^{7.6} K and 10^{8.0} K, thus demonstrating that our model is robust to extrapolation. 
5.4. Matter, CDM, gas and stars
In Fig. 10 we show a model that is the result of a joint fit to all possible combinations of matter, CDM, gas and star spectra. There are 10 distinct power spectra for these 4 fields. In our halo model, the matter field is the sum of the CDM, gas and stars, and so naturally gets upweighted in our fitting when it is included alongside its constituent parts. We see that our model gets reasonable matches to the autospectra shown in Fig. 10 and we have also checked that the crossspectra are well matched (see Fig. 11), with the error always being roughly (but not exactly) the mean of the errors of the autospectra of the two contributing fields. We see that the gas spectrum is particularly well matched, which gives us confidence in our choice of the KS profile for the gas density profile. The match to the CDM–CDM response looks less good in Fig. 10, but note that the deviations from this being a constant are very small, of order a few percent. Physically this implies that the CDM is not deformed very significantly when going from a gravity only to a hydrodynamic model. Indeed, the suppression in the matter–matter power spectrum is driven predominantly by the (lack of) contribution from the gas, and then the upturn at smaller scales mainly derives from the contribution of the stars (see Fig. 6, as well as van Daalen et al. 2020). We also looked at additionally fitting the parameters ϵ_{2}, β and σ_{*}, but found that this did not produce a significant improvement in the quality of the fit.
Fig. 10. Same as Figs. 7–9, but for the bestfitting model (4) to all 10 possible power spectra of matter, CDM, gas and star fields, although we only show the 4 autospectra here. The top half of the figure show power spectra while the bottom half show response functions. Columns show different redshifts. Different colours are different feedback temperatures: T = 10^{7.6} (blue), T = 10^{7.8} (grey) and T = 10^{8.0} K (red). Points (with errors) are from BAHAMAS while lines are from our halo model. 
5.5. Model accuracy
In Fig. 11 we show the full error matrix for the AGN simulation for each of models (2), (3) and (4) but calculated for every possible power spectra of the 5 fields we have available. The matrix is 5 × 5 but only 15 of the 25 elements are independent due to symmetry. Each cell shows the relative difference for the model cross spectrum compared to the simulation, averaged linearly over z between 0 and 1 and logarithmically over k between ≃0.015 and 7 h Mpc^{−1}. Error matrices for AGNLO and AGNHI are not shown, but are similar. We see that when we fitted only the matter–matter spectrum (2) we still have a reasonable match to the CDM, gas and stars spectra and combinations thereof. When we go from (2) to (3) and include the matter–electron pressure spectrum we improve the match to the matter–electron pressure, but at the cost of power spectra that involve the gas. This demonstrates a possible problem in our modelling in how we simultaneously model the gas density and pressure: When we go from case (2) to (4) and try to match matter, CDM, gas and stars power spectra, but ignore electron pressure we match the gas spectrum well, but this success comes at the expense of the pressure, which is ignored in this fit.
Fig. 11. Full error matrix for our fitted models calculated for each field pair for the 10^{7.8} K AGN simulation. The error is averaged linearly over redshift between z = 0 and 1 and over k logarithmically between ≃0.015 and 7 h Mpc^{−1}. This matrix is symmetric because power spectra are symmetric when swapping the field labels. We show the matrix for the three models presented in this work: matter (2, left); matter, electron pressure (3, centre); matter, CDM, gas, stars (4, right). The squares that contain dots are the particular crossspectra that are fitted in the model. In the lefthand plot, we see that in the matter model we get a reasonable fit to the power spectra all of the constituent fields. In the central plot we get a better match to to the matter–electron pressure spectrum than in the lefthand panel, but this is at the expense of any power spectra that directly involves the gas field, which shows some tension in our model between the gas and pressure modelling. In the righthand panel, where we fitted matter, CDM, gas and stars, we see that all the matter spectra are reasonable, and better than in the lefthand panel, but achieving this is at the expense of any spectra that involves the pressure field. 
6. Discussion
6.1. Summary
We have presented a hydrodynamical halo model that can be used to make analytical calculations of power spectra for combinations of the matter, CDM, gas and star density fields, as well as for crossspectra of these fields with the electronpressure field. This model will be used in the future to provide constraints on both cosmology and feedback physics using a combination of lensing and tSZ data. The model has a number of tuneable parameters that pertain to gas physics. We provided four possible sets of values for these parameters that provide good fits to power spectra for: (1) stars, (2) matter, (3) matter and electron pressure and (4) matter, CDM, gas and stars. Parameters of these models are all fitted to power spectra measured from the BAHAMAS simulations for three different AGN feedback strengths, which in this case is parameterised by a single subgrid parameter: the AGN heating temperature. By interpolating between model parameters as a function of AGN heating temperature we suggest that our model can be used in a cosmological analysis and that it is robust to extrapolation. In particular, model (2) will be useful in lensingonly studies as it is able to reproduce the deficit in matter–matter power seen in hydrodynamic simulations, compared to the DMONLY case, at the percent level. Model (3) provides a slightly worse match to the mater–matter power (accurate to 2%), but captures the temperature variations in the matter–electron pressure spectrum that can be directly integrated to calculate the lensing–tSZ cross correlation (accurate to 15%) as well as the lensing–lensing correlation. This accuracy is sufficient for forthcoming lensing and tSZ data sets (Tröster et al., in prep.).
6.2. Identified problems
6.2.1. Electron pressure errors
Our model is less accurate for matter–electron pressure (15%) compared to matter–matter (2%). We suspect that this is partly because we report the error averaged over log k from 0.015 to 7 h Mpc^{−1} and linearly over z from 0 to 1. The matter–electron pressure response is noisier at large scales compared to the matter spectra, and this introduces an unavoidable source of error for our necessarily smooth halo model prediction. We may also be seeing the fact that the relatively small box size of BAHAMAS means that the response (or power) has not converged for k < 1 h Mpc^{−1}, and it is possible that our model would provide a better fit to data from a significantly larger simulation.
6.2.2. Response approach
In this paper, we fitted to the power spectrum response, which we created by dividing a measured hydrodynamical power spectrum by a gravityonly spectrum. This response cancels the Gaussian error in power spectra at large scales for the matter spectra, but it may create complicated correlations at smaller scales. In particular, this response may not be appropriate for the electronpressure spectra, even at large scales, since the largescale variance does not mimic that of the matter spectra. We fitted parameters of our models to the response using a leastsquares approach, with no k or z dependent error. In the future we hope to have access to more and larger simulations that we would need in order investigate this more thoroughly and so that we may derive an appropriate error bar.
6.2.3. Electron pressure auto spectrum
While completing this work we tried to generate a model for the electron pressure autospectrum, together with matter–matter and matter–electron pressure; but this proved to be difficult and was abandoned. We discovered that the realisationtorealisation variance in the electron pressure autospectrum measured from BAHAMAS was huge, and dwarfed the difference between feedback strengths for k < 2 h Mpc^{−1}, which left us with very little signal to fit. This in turn suggests that boxes of much greater volume than the BAHAMAS (400 h^{−1} Mpc)^{3} must be used to ascertain the feedback dependence of this spectrum on large scales. Therefore, our model should not be used to compute tSZ–tSZ power spectrum at this stage. However, with larger hydrodynamical simulations being run in the future, our approach should still be able to be used to match the power spectrum by refitting parameters. We note that it would also be useful to have more accurate calibration data for the simulations to match, as this is one of the limiting factors for BAHAMAS at the moment.
6.2.4. Joint model
We also tried to generate a joint model for matter, CDM, gas, stars and electron pressure, which has not been presented in this paper because it did not work well. We found that it was particularly difficult to simultaneously model spectra that involved the electron pressure and those that involve the gas density with the same underlying physical model. In the current modelling, these two fields are joined in two different ways: first via M_{0}, which determines the overall gas abundance as a function of halo mass, and second via Γ, the polytropic index of the KS gas profile, which jointly determines the density and pressure profiles. The KS profile is the correct model for undisturbed gas that is in hydrostatic equilibrium in an NFW halo, but this state will clearly be violated in the presence of AGN feedback. It was shown in Yan et al. (2020) that the KS model could provide a reasonable match to gasdensity profiles in the BAHAMAS simulations, but it is likely that the polytropic link between density and pressure is not respected in the presence of feedback (Battaglia et al. 2012a,b). It is probable that a better model could be generated by either changing the gas profile, for example by explicitly including heated gas, or by weakening the links between the density and pressure modelling. However, note that some link is necessary if one wants to simultaneously constrain feedback and cosmology with crosscorrelation data. It is also possible that some of the inaccuracy stems from our treatment of nonthermal pressure within clusters. We have a single multiplicative parameter, α, which increases or decreases gas temperature independent of halo mass. In future it may be useful to investigate physical models for nonthermal pressure in clusters, such as those presented in Nelson et al. (2014) or Shi et al. (2015).
6.2.5. Nonlinear halo bias
The halo model provides a good description of the power spectra for twohalo dominated regimes and the onehalo dominated regimes, but it is frustrating that it fails in the transition region between these two regimes. In retrospect, the reason for this is obvious: The assumption made is that haloes are linearly tracing an underlying linear matter field, and this assumption clearly breaks down at scales that correspond to this transition: k ∼ 0.5 h Mpc^{−1} at z = 0 for matter–matter. Linear bias is not valid at this wavenumber (Desjacques et al. 2018) and the twohalo term predicted by the model at these and smaller scales will be, quite simply, wrong. We suggest that a fruitful line of future inquiry would be to improve the treatment of halo bias within the standard halomodel framework (e.g. Smith et al. 2007). We have somewhat sidestepped this issue in this paper by fitting the BAHAMAS response functions with halomodel responses, rather than the power spectra directly. This is certainly of benefit because we cancel the Gaussian noise present at largescales in the simulations. This makes the halo model response an almost perfect tool for understanding matter–matter power response for dark energy models measured in simulations without hydrodynamics (Mead 2017; Cataneo et al. 2019) when the underlying linear spectrum is fixed. This would also be true in hydrodynamic simulations if we could be assured that the troubles in the transition region manifest themselves in the same way, and with the same scale dependence, for other power spectra as they do for matter–matter. Unfortunately, from Fig. 2 we can see that this is not exactly the case, since the transition region occurs at slightly different scales for the different power spectra (particularly the scale for matter–electron pressure is quite different from matter–matter): there simply cannot be some universal kdependent prescription that one could apply to solve all problems in all cases. It is therefore possible that the response is not the optimal quantity to consider when one considers the different power spectra measured in hydrodynamic simulations.
6.2.6. Expelled gas modelling
One further shortcoming in our modelling is in our treatment of the gas expelled from haloes. In this paper we do not consider any sort of onehalo contribution arising from expelled gas, and instead account for the effects of this gas only in the twohalo term. In reality, we know that gas expelled from a halo will still be correlated with that halo, and will be found in an extended shell outside the virial radius. Since this shell is larger than the halo virial radius, it may provide a significant contribution to the power at scales that correspond to the transition between the two and onehalo terms in the classic halo model, and we have so far ignored this. In the early stages of this work we attempted to explicitly include a gas component outside the halo virial radius, but we found that this introduces some problems with the standard halomodel formalism: In the calculation, it is necessary to truncate the density profile of each NFW halo at the virial radius, otherwise they each have an infinite mass. If the gas halo is allowed to extend beyond the virial radius then we introduce an explicit anticorrelation between gas and CDM because we have gas in regions where we explicitly have no CDM. This anticorrelation allows the CDM–gas cross spectrum to be negative in some regions of parameter space. While negative crossspectra are not problematic in general, in cosmology it would be unrealistic due to the general structure of the cosmic web. Indeed, such a negative cross spectrum is not seen in simulations. It may be possible to improve this by using a different CDM halo profile with a finite mass as the radius tends to infinity, but this introduces difficulty in deciding which mass to associate with which halo when identifying haloes in simulations. The fact that we treat the ejected gas in the twohalo term as a simple multiplicative bias also generates another problem: the twohalo term for ejected gas is not suppressed at small scales relative to the shape of linear theory. For other fields, it is always suppressed by a weight that depends on the halo profile and halo bias (Eq. (1)). The details of this twohalo term suppression cannot be taken seriously because they occur at scales where the halo bias is no longer linear, however, this suppression does have the virtue of ensuring that the twohalo term is always much smaller than the onehalo term at small scales. Because the gas profiles themselves lack smallscale structure, the twohalo term for gas (and electron pressure) can be larger than the onehalo term at very small scales (k ∼ 100 h Mpc^{−1}). Clearly this is unphysical and would need to be addressed in further work.
6.2.7. Important halo masses
A further potential problem for this approach is the differing halo masses that are important for the tSZ effect compared to those important for galaxy weak lensing. As shown in Fig. 3, the matter–electron pressure spectrum derives most contribution from highermass haloes than matter–matter. It is therefore these haloes that we learn most about from the lensingtSZ cross correlation. However, it is lowermass haloes that mostly affect the suppression in the power spectrum that is important for the lensinglensing spectrum. A fruitful avenue for future research may be to learn more directly about the gas distribution from the haloes that are most relevant, those of 10^{13} to 10^{14} h^{−1} M_{⊙}. It is difficult to do this using the tSZlensing cross correlation directly, and instead one may have to use some sort of clipping (e.g. Simpson et al. 2011, 2013; Hill & Pajer 2013; Giblin et al. 2018) technique to mask the influence of the highest mass haloes. We note that there is already a great deal known about the impact of baryons on the matter distribution from tSZ and Xray studies of galaxies, groups, and clusters from targeted observational studies. One could imagine using the model presented in this paper to look at both the lensing and tSZ signal (or lensing–tSZ) around haloes of specific masses. This would add information about the signal as it arises from different mass bins, somehow akin to the information in Fig. 3. However, this would necessitate having reliable measurements of individual halo masses and would require modelling effects such as halo miscentring (e.g. Yan et al. 2020). It would also need to be ascertained whether the halo model can be trusted for calculations that are binned in halo mass and how problems with these calculations compare to problems in the standard when integrated over all halo masses.
6.3. Conclusion
We intend to use this model in the near future to analyse the cross correlation between the tSZ map from Planck and weak lensing from the KiDS data, as well as the lensing auto correlation data (Tröster et al., in prep). This work has the potential to simultaneously learn about cosmology and the magnitude of AGN feedback in the Universe and therefore to inform galaxy formation modelling and future hydrodynamical simulations as well as to shrink error bars on cosmological parameters. As alluded to in this paper, the tSZ data itself also has the power to provide tight constraints on σ_{8} and initial tests show that this could break degeneracies that exist in a traditional lensing analysis.
So far, we have only presented the bestfitting model parameters for each scenario and we have not investigated errorson, or degeneracies between, our model parameters, but we hope to be able to do this in future. We note that this would require us to create meaningful error bars on our response measurements.
It is clear that the fitted parameters presented here are very BAHAMAS specific, which may worry some readers who would prefer a model that could encompass a wider range of hydrodynamical simulations. However, we think that our model is general enough that it could be refitted to other simulations in the future. We also note van Daalen et al. (2020) has demonstrated a strong correlation between the matter–matter response function and the gas fraction in haloes of ∼10^{14} h^{−1} M_{⊙} and that the BAHAMAS simulations reproduce observations of gas fractions in haloes of this mass to within the observational error bar. Given that the gas content in haloes of this mass also is important for tSZ suggests that tuning our model to BAHAMAS was a reasonable first choice. In the future, it may be possible to rewrite our results in terms of more physical quantities, such as the halo gas fraction (e.g. Chisari et al. 2019; van Daalen et al. 2020), rather than the unphysical AGN heating temperature. It would also be desirable to confront our model with different hydrodynamical simulations, run with different underlying assumptions and methods, for example those run with the Adaptive Refinement Tree code (Kravtsov et al. 2002; Nagai et al. 2007) or those of ILLUSTRIS (Vogelsberger et al. 2014).
In the future, it may be possible to use this halo model to make predictions for higherorder statistics and these may lead to increased constraints on the model parameters. We also envisage the model being used to make predictions for other “hydrodynamical” observable quantities, such as npoint statistics that involve the Xray field, which is a different function of gas temperature and density compared to the tSZ and provides complimentary information about the state of gas in the universe. The utility of the halo model that we present does not require its parameters to be fitted to simulation data, this was simply one option, and the option that we decided to pursue in this paper. In the future the model can be used to provide quick estimates for the effect of new feedback scenarios that have yet to be simulated, or to explore nonstandard cosmological scenarios that have yet to be hydrodynamical simulated. It would even be possible to investigate the nonlinear coupling of feedback with nonstandard cosmological scenarios.
Often simulations that only consider the gravitational interaction are called “darkmatter only”. We prefer “gravity only” because these simulations do contain baryons, and the initial conditions are set up to account for baryons, but they exclusively consider the gravitational interaction between particles when they are subsequently run.
Achieving these limits is difficult numerically because of the large amount of mass contained in lowmass haloes according to most popular mass functions. Special care must be taken with the twohalo integral in the case of power spectra that involve the matter field. See Appendix A.
The nonstandard factor of Ω_{m}(z) in the denominator of Eq. (22) arises because we work with respect to the matter density, rather than critical density.
If the gas were comprised only of ionised hydrogen and helium, with a hydrogen mass fraction f_{H}, we would have μ_{p} = 4/(3 + 5f_{H}) and μ_{e} = 2/(1 + f_{H}). For the BAHAMAS with f_{H} = 0.752, μ_{e} ≃ 0.59 and μ_{e} ≃ 1.14. However, we adopt the values μ_{p} = 0.61 and μ_{e} = 1.17 to be compatible with gas metallicity in the simulations and the way the electron pressure field was measured in post processesing.
SGENIC: https://github.com/sbird/SGenIC
Acknowledgments
AJM acknowledges support from a CITA National Fellowship and (together with LVW) from NSERC. AJM and TT have received funding from the Horizon 2020 research and innovation programme of the European Union under Marie SkłodowskaCurie grant agreements Nos. 702971 (AJM) and 797794 (TT). TT and CH acknowledge support from the European Research Council under grant number 647112, while IGM acknowledges the same funding under grant number 769130. CH also acknowledges support from the Max Planck Society and the Alexander von Humboldt Foundation in the framework of the Max PlanckHumboldt Research Award endowed by the Federal Ministry of Education and Research. We thank an anonymous referee for useful comments that improved the quality of this paper.
References
 Addison, G. E., Dunkley, J., & Spergel, D. N. 2012, MNRAS, 427, 1741 [NASA ADS] [CrossRef] [Google Scholar]
 Addison, G. E., Dunkley, J., & Bond, J. R. 2013, MNRAS, 436, 1896 [Google Scholar]
 Agarwal, S., Abdalla, F. B., Feldman, H. A., Lahav, O., & Thomas, S. A. 2012, MNRAS, 424, 1409 [CrossRef] [Google Scholar]
 Angulo, R. E., Zennaro, M., Contreras, S., et al. 2020, ArXiv eprints [arXiv:2004.06245] [Google Scholar]
 Aricò, G., Angulo, R. E., HernándezMonteagudo, C., et al. 2020, MNRAS, 495, 4800 [CrossRef] [Google Scholar]
 Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bartelmann, M., & Schneider, P. 2001, Phys. Rep., 340, 291 [NASA ADS] [CrossRef] [Google Scholar]
 Battaglia, N., Bond, J. R., Pfrommer, C., Sievers, J. L., & Sijacki, D. 2010, ApJ, 725, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Battaglia, N., Bond, J. R., Pfrommer, C., & Sievers, J. L. 2012a, ApJ, 758, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Battaglia, N., Bond, J. R., Pfrommer, C., & Sievers, J. L. 2012b, ApJ, 758, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Battaglia, N., Hill, J. C., & Murray, N. 2015, ApJ, 812, 154 [CrossRef] [Google Scholar]
 Baxter, E. J., Omori, Y., Chang, C., et al. 2019, Phys. Rev. D, 99, 023508 [NASA ADS] [CrossRef] [Google Scholar]
 Bird, S., Feng, Y., Pedersen, C., & FontRibera, A. 2020, JCAP, 2020, 002 [CrossRef] [Google Scholar]
 Bolliet, B., Comis, B., Komatsu, E., & MacíasPérez, J. F. 2018, MNRAS, 477, 4957 [NASA ADS] [CrossRef] [Google Scholar]
 Bond, J. R., Contaldi, C. R., Pen, U.L., et al. 2005, ApJ, 626, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Cacciato, M., Lahav, O., van den Bosch, F. C., Hoekstra, H., & Dekel, A. 2012, MNRAS, 426, 566 [NASA ADS] [CrossRef] [Google Scholar]
 Cataneo, M., Lombriser, L., Heymans, C., et al. 2019, MNRAS, 488, 2121 [NASA ADS] [CrossRef] [Google Scholar]
 Cen, R., & Ostriker, J. P. 1999, ApJ, 514, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Chisari, N. E., Richardson, M. L. A., Devriendt, J., et al. 2018, MNRAS, 480, 3962 [NASA ADS] [CrossRef] [Google Scholar]
 Chisari, N. E., Mead, A. J., Joudaki, S., et al. 2019, Open J. Astrophys., 2, 4 [CrossRef] [Google Scholar]
 Cooray, A., & Sheth, R. 2002, Phys. Rep., 372, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Copeland, D., Taylor, A., & Hall, A. 2018, MNRAS, 480, 2247 [CrossRef] [Google Scholar]
 Courtin, J., Rasera, Y., Alimi, J.M., et al. 2011, MNRAS, 410, 1911 [NASA ADS] [Google Scholar]
 Dai, B., Feng, Y., & Seljak, U. 2018, JCAP, 11, 009 [CrossRef] [Google Scholar]
 Debackere, S. N. B., Schaye, J., & Hoekstra, H. 2020, MNRAS, 492, 2285 [CrossRef] [Google Scholar]
 de Graaff, A., Cai, Y.C., Heymans, C., & Peacock, J. A. 2019, A&A, 624, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Desjacques, V., Jeong, D., & Schmidt, F. 2018, Phys. Rep., 733, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Despali, G., Giocoli, C., Angulo, R. E., et al. 2016, MNRAS, 456, 2486 [NASA ADS] [CrossRef] [Google Scholar]
 Dolag, K., Borgani, S., Murante, G., & Springel, V. 2009, MNRAS, 399, 497 [Google Scholar]
 Dolag, K., Komatsu, E., & Sunyaev, R. 2016, MNRAS, 463, 1797 [CrossRef] [Google Scholar]
 Duffy, A. R., Schaye, J., Kay, S. T., & Dalla Vecchia, C. 2008, MNRAS, 390, L64 [NASA ADS] [CrossRef] [Google Scholar]
 Eifler, T., Krause, E., Dodelson, S., et al. 2015, MNRAS, 454, 2451 [NASA ADS] [CrossRef] [Google Scholar]
 Fedeli, C. 2014, JCAP, 4, 28 [Google Scholar]
 Fedeli, C., Dolag, K., & Moscardini, L. 2012, MNRAS, 419, 1588 [NASA ADS] [CrossRef] [Google Scholar]
 Fedeli, C., Semboloni, E., Velliscig, M., et al. 2014, JCAP, 8, 28 [Google Scholar]
 Foreman, S., Becker, M. R., & Wechsler, R. H. 2016, MNRAS, 463, 3326 [NASA ADS] [CrossRef] [Google Scholar]
 Giblin, B., Heymans, C., HarnoisDéraps, J., et al. 2018, MNRAS, 480, 5529 [NASA ADS] [CrossRef] [Google Scholar]
 Giocoli, C., Bartelmann, M., Sheth, R. K., & Cacciato, M. 2010, MNRAS, 408, 300 [NASA ADS] [CrossRef] [Google Scholar]
 Gonzalez, A. H., Sivanandam, S., Zabludoff, A. I., & Zaritsky, D. 2013, ApJ, 778, 14 [NASA ADS] [CrossRef] [Google Scholar]
 HarnoisDéraps, J., van Waerbeke, L., Viola, M., & Heymans, C. 2015, MNRAS, 450, 1212 [NASA ADS] [CrossRef] [Google Scholar]
 Hildebrandt, H., Viola, M., Heymans, C., et al. 2017, MNRAS, 465, 1454 [Google Scholar]
 Hildebrandt, H., Köhlinger, F., van den Busch, J. L., et al. 2020, A&A, 633, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hill, J. C., & Pajer, E. 2013, Phys. Rev. D, 88, 063526 [NASA ADS] [CrossRef] [Google Scholar]
 Hill, J. C., & Spergel, D. N. 2014, JCAP, 2, 030 [NASA ADS] [CrossRef] [Google Scholar]
 Hill, J. C., Baxter, E. J., Lidz, A., Greco, J. P., & Jain, B. 2018, Phys. Rev. D, 97, 083501 [CrossRef] [Google Scholar]
 Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Hojjati, A., McCarthy, I. G., HarnoisDeraps, J., et al. 2015, JCAP, 10, 047 [CrossRef] [Google Scholar]
 Hojjati, A., Tröster, T., HarnoisDéraps, J., et al. 2017, MNRAS, 471, 1565 [CrossRef] [Google Scholar]
 Holder, G. P., McCarthy, I. G., & Babul, A. 2007, MNRAS, 382, 1697 [NASA ADS] [Google Scholar]
 Horowitz, B., & Seljak, U. 2017, MNRAS, 469, 394 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Huang, H.J., Eifler, T., Mandelbaum, R., & Dodelson, S. 2019, MNRAS, 488, 1652 [NASA ADS] [CrossRef] [Google Scholar]
 Huffenberger, K. M., & Seljak, U. 2003, MNRAS, 340, 1199 [NASA ADS] [CrossRef] [Google Scholar]
 Itoh, N., Kohyama, Y., & Nozawa, S. 1998, ApJ, 502, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Jain, B., & Seljak, U. 1997, ApJ, 484, 560 [NASA ADS] [CrossRef] [Google Scholar]
 Joudaki, S., Blake, C., Heymans, C., et al. 2017a, MNRAS, 465, 2033 [NASA ADS] [CrossRef] [Google Scholar]
 Joudaki, S., Mead, A., Blake, C., et al. 2017b, MNRAS, 471, 1259 [NASA ADS] [CrossRef] [Google Scholar]
 Kaiser, N. 1992, ApJ, 388, 272 [NASA ADS] [CrossRef] [Google Scholar]
 Kilbinger, M. 2015, Rep. Prog. Phys., 78, 086901 [Google Scholar]
 Kilbinger, M., Heymans, C., Asgari, M., et al. 2017, MNRAS, 472, 2126 [NASA ADS] [CrossRef] [Google Scholar]
 Knabenhans, M., Stadel, J., Marelli, S., et al. 2019, MNRAS, 484, 5509 [NASA ADS] [CrossRef] [Google Scholar]
 Komatsu, E., & Seljak, U. 2001, MNRAS, 327, 1353 [Google Scholar]
 Komatsu, E., & Seljak, U. 2002, MNRAS, 336, 1256 [NASA ADS] [CrossRef] [Google Scholar]
 Kravtsov, A. V., Klypin, A., & Hoffman, Y. 2002, ApJ, 571, 563 [NASA ADS] [CrossRef] [Google Scholar]
 Kravtsov, A. V., Vikhlinin, A. A., & Meshcheryakov, A. V. 2018, Astron. Lett., 44, 8 [NASA ADS] [CrossRef] [Google Scholar]
 Lawrence, E., Heitmann, K., White, M., et al. 2010, ApJ, 713, 1322 [NASA ADS] [CrossRef] [Google Scholar]
 Lawrence, E., Heitmann, K., Kwan, J., et al. 2017, ApJ, 847, 50 [CrossRef] [Google Scholar]
 Leauthaud, A., Tinker, J., Bundy, K., et al. 2012, ApJ, 744, 159 [NASA ADS] [CrossRef] [Google Scholar]
 Le Brun, A. M. C., McCarthy, I. G., Schaye, J., & Ponman, T. J. 2014, MNRAS, 441, 1270 [NASA ADS] [CrossRef] [Google Scholar]
 Le Brun, A. M. C., McCarthy, I. G., & Melin, J.B. 2015, MNRAS, 451, 3868 [NASA ADS] [CrossRef] [Google Scholar]
 Lemos, P., Challinor, A., & Efstathiou, G. 2017, JCAP, 5, 014 [CrossRef] [Google Scholar]
 Lesgourgues, J. 2011, ArXiv eprints [arXiv:1104.2934] [Google Scholar]
 Lo Verde, M., Miller, A., Shandera, S., & Verde, L. 2008, JCAP, 4, 14 [Google Scholar]
 Ma, Y.Z., Van Waerbeke, L., Hinshaw, G., et al. 2015, JCAP, 9, 046 [NASA ADS] [CrossRef] [Google Scholar]
 MacCrann, N., Zuntz, J., Bridle, S., Jain, B., & Becker, M. R. 2015, MNRAS, 451, 2877 [NASA ADS] [CrossRef] [Google Scholar]
 MacCrann, N., Aleksić, J., Amara, A., et al. 2017, MNRAS, 465, 2567 [CrossRef] [Google Scholar]
 Makiya, R., Ando, S., & Komatsu, E. 2018, MNRAS, 480, 3928 [NASA ADS] [CrossRef] [Google Scholar]
 Martizzi, D., Teyssier, R., Moore, B., & Wentz, T. 2012, MNRAS, 422, 3081 [NASA ADS] [CrossRef] [Google Scholar]
 Martizzi, D., Teyssier, R., & Moore, B. 2013, MNRAS, 432, 1947 [NASA ADS] [CrossRef] [Google Scholar]
 McCarthy, I. G., Schaye, J., Ponman, T. J., et al. 2010, MNRAS, 406, 822 [NASA ADS] [Google Scholar]
 McCarthy, I. G., Schaye, J., Bower, R. G., et al. 2011, MNRAS, 412, 1965 [NASA ADS] [CrossRef] [Google Scholar]
 McCarthy, I. G., Le Brun, A. M. C., Schaye, J., & Holder, G. P. 2014, MNRAS, 440, 3645 [NASA ADS] [CrossRef] [Google Scholar]
 McCarthy, I. G., Schaye, J., Bird, S., & Le Brun, A. M. C. 2017, MNRAS, 465, 2936 [NASA ADS] [CrossRef] [Google Scholar]
 McCarthy, I. G., Bird, S., Schaye, J., et al. 2018, MNRAS, 476, 2999 [NASA ADS] [CrossRef] [Google Scholar]
 Mead, A. J. 2017, MNRAS, 464, 1282 [CrossRef] [Google Scholar]
 Mead, A. J., Peacock, J. A., Heymans, C., Joudaki, S., & Heavens, A. F. 2015, MNRAS, 454, 1958 [NASA ADS] [CrossRef] [Google Scholar]
 Mead, A. J., Heymans, C., Lombriser, L., et al. 2016, MNRAS, 459, 1468 [Google Scholar]
 Mo, H. J., & White, S. D. M. 1996, MNRAS, 282, 347 [NASA ADS] [CrossRef] [Google Scholar]
 Mohammed, I., Martizzi, D., Teyssier, R., & Amara, A. 2014, ArXiv eprints [arXiv:1410.6826] [Google Scholar]
 Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121 [Google Scholar]
 Mummery, B. O., McCarthy, I. G., Bird, S., & Schaye, J. 2017, MNRAS, 471, 227 [CrossRef] [Google Scholar]
 Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007, ApJ, 668, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Nakamura, T. T., & Suto, Y. 1997, Prog. Theor. Phys., 97, 49 [NASA ADS] [CrossRef] [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [NASA ADS] [CrossRef] [Google Scholar]
 Nelder, J. A., & Mead, R. 1965, Comput. J., 7, 308 [Google Scholar]
 Nelson, K., Lau, E. T., & Nagai, D. 2014, ApJ, 792, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Nozawa, S., Itoh, N., Kawana, Y., & Kohyama, Y. 2000, ApJ, 536, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Omori, Y., Giannantonio, T., Porredon, A., et al. 2019, Phys. Rev. D, 100, 043501 [CrossRef] [Google Scholar]
 Osato, K., Flender, S., Nagai, D., Shirasaki, M., & Yoshida, N. 2018, MNRAS, 475, 532 [CrossRef] [Google Scholar]
 Osato, K., Shirasaki, M., Miyatake, H., et al. 2020, MNRAS, 492, 4780 [CrossRef] [Google Scholar]
 Peacock, J. A., & Smith, R. E. 2000, MNRAS, 318, 1144 [NASA ADS] [CrossRef] [Google Scholar]
 Planck Collaboration XXII. 2016, A&A, 594, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIV. 2016, A&A, 594, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puchwein, E., & Springel, V. 2013, MNRAS, 428, 2966 [NASA ADS] [CrossRef] [Google Scholar]
 Rabold, M., & Teyssier, R. 2017, MNRAS, 467, 3188 [NASA ADS] [CrossRef] [Google Scholar]
 Refregier, A., & Teyssier, R. 2002, Phys. Rev. D, 66, 043002 [NASA ADS] [CrossRef] [Google Scholar]
 Remazeilles, M., Bolliet, B., Rotti, A., & Chluba, J. 2019, MNRAS, 483, 3459 [NASA ADS] [CrossRef] [Google Scholar]
 Roncarelli, M., Moscardini, L., Tozzi, P., et al. 2006, MNRAS, 368, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Roncarelli, M., Moscardini, L., Borgani, S., & Dolag, K. 2007, MNRAS, 378, 1259 [NASA ADS] [CrossRef] [Google Scholar]
 Rudd, D. H., Zentner, A. R., & Kravtsov, A. V. 2008, ApJ, 672, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Schaye, J., Dalla Vecchia, C., Booth, C. M., et al. 2010, MNRAS, 402, 1536 [Google Scholar]
 Schmidt, F. 2016, Phys. Rev. D, 93, 063512 [CrossRef] [Google Scholar]
 Schneider, A., & Teyssier, R. 2015, JCAP, 12, 049 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, A., Teyssier, R., Potter, D., et al. 2016, JCAP, 2016, 047 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, A., Teyssier, R., Stadel, J., et al. 2019, JCAP, 2019, 020 [CrossRef] [Google Scholar]
 Scoccimarro, R., Sheth, R. K., Hui, L., & Jain, B. 2001, ApJ, 546, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Seljak, U. 2000, MNRAS, 318, 203 [NASA ADS] [CrossRef] [Google Scholar]
 Semboloni, E., Hoekstra, H., Schaye, J., van Daalen, M. P., & McCarthy, I. G. 2011, MNRAS, 417, 2020 [NASA ADS] [CrossRef] [Google Scholar]
 Semboloni, E., Hoekstra, H., & Schaye, J. 2013, MNRAS, 434, 148 [NASA ADS] [CrossRef] [Google Scholar]
 Shaw, L. D., Nagai, D., Bhattacharya, S., & Lau, E. T. 2010, ApJ, 725, 1452 [NASA ADS] [CrossRef] [Google Scholar]
 Sheth, R. K., & Tormen, G. 1999, MNRAS, 308, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Sheth, R. K., Mo, H. J., & Tormen, G. 2001, MNRAS, 323, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Shi, X., Komatsu, E., Nelson, K., & Nagai, D. 2015, MNRAS, 448, 1020 [NASA ADS] [CrossRef] [Google Scholar]
 Sijacki, D., Vogelsberger, M., Genel, S., et al. 2015, MNRAS, 452, 575 [NASA ADS] [CrossRef] [Google Scholar]
 Simpson, F., James, J. B., Heavens, A. F., & Heymans, C. 2011, Phys. Rev. Let., 107, 271301 [NASA ADS] [CrossRef] [Google Scholar]
 Simpson, F., Heavens, A. F., & Heymans, C. 2013, Phys. Rev. D, 88, 083510 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, R. E., & Markovic, K. 2011, Phys. Rev. D, 84, 063507 [CrossRef] [Google Scholar]
 Smith, R. E., Peacock, J. A., Jenkins, A., et al. 2003, MNRAS, 341, 1311 [Google Scholar]
 Smith, R. E., Scoccimarro, R., & Sheth, R. K. 2007, Phys. Rev. D, 75, 063512 [NASA ADS] [CrossRef] [Google Scholar]
 Sun, M., Voit, G. M., Donahue, M., et al. 2009, ApJ, 693, 1142 [NASA ADS] [CrossRef] [Google Scholar]
 Suto, Y., Sasaki, S., & Makino, N. 1998, ApJ, 509, 544 [NASA ADS] [CrossRef] [Google Scholar]
 Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152 [NASA ADS] [CrossRef] [Google Scholar]
 Tanimura, H., Hinshaw, G., McCarthy, I. G., et al. 2019a, MNRAS, 483, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Tanimura, H., Aghanim, N., Douspis, M., Beelen, A., & Bonjean, V. 2019b, A&A, 625, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tanimura, H., Hinshaw, G., McCarthy, I. G., et al. 2020, MNRAS, 491, 2318 [Google Scholar]
 Tinker, J. L., Weinberg, D. H., Zheng, Z., & Zehavi, I. 2005, ApJ, 631, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709 [NASA ADS] [CrossRef] [Google Scholar]
 Tinker, J. L., Robertson, B. E., Kravtsov, A. V., et al. 2010, ApJ, 724, 878 [NASA ADS] [CrossRef] [Google Scholar]
 Trac, H., Bode, P., & Ostriker, J. P. 2011, ApJ, 727, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Valageas, P., & Nishimichi, T. 2011, A&A, 527, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van Daalen, M. P., Schaye, J., Booth, C. M., & Dalla Vecchia, C. 2011, MNRAS, 415, 3649 [NASA ADS] [CrossRef] [Google Scholar]
 van Daalen, M. P., Schaye, J., McCarthy, I. G., Booth, C. M., & Vecchia, C. D. 2014, MNRAS, 440, 2997 [CrossRef] [Google Scholar]
 van Daalen, M. P., McCarthy, I. G., & Schaye, J. 2020, MNRAS, 491, 2424 [CrossRef] [Google Scholar]
 van den Bosch, F. C., More, S., Cacciato, M., Mo, H., & Yang, X. 2013, MNRAS, 430, 725 [NASA ADS] [CrossRef] [Google Scholar]
 Van Waerbeke, L., Hinshaw, G., & Murray, N. 2014, Phys. Rev. D, 89, 023508 [NASA ADS] [CrossRef] [Google Scholar]
 Velliscig, M., van Daalen, M. P., Schaye, J., et al. 2014, MNRAS, 442, 2641 [NASA ADS] [CrossRef] [Google Scholar]
 Velliscig, M., Cacciato, M., Schaye, J., et al. 2015, MNRAS, 453, 721 [NASA ADS] [CrossRef] [Google Scholar]
 Vikhlinin, A., Burenin, R. A., Ebeling, H., et al. 2009, ApJ, 692, 1033 [NASA ADS] [CrossRef] [Google Scholar]
 Vogelsberger, M., Genel, S., Springel, V., et al. 2014, Nature, 509, 177 [NASA ADS] [CrossRef] [Google Scholar]
 White, M. 2004, Astropart. Phys., 22, 211 [NASA ADS] [CrossRef] [Google Scholar]
 Yan, Z., Hojjati, A., Tröster, T., Hinshaw, G., & van Waerbeke, L. 2019, ApJ, 884, 139 [CrossRef] [Google Scholar]
 Yan, Z., Raza, N., Van Waerbeke, L., et al. 2020, MNRAS, 493, 1120 [CrossRef] [Google Scholar]
 Yoo, J., Tinker, J. L., Weinberg, D. H., et al. 2006, ApJ, 652, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Yoon, M., Jee, M. J., Tyson, J. A., et al. 2019, ApJ, 870, 111 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Evaluation of the twohalo term
There is some confusion in the literature about how exactly to numerically evaluate the integral for the halo model twohalo term:
$$\begin{array}{c}\hfill {P}_{2\mathrm{H},uv}(k)={P}_{\mathrm{lin}}(k){\displaystyle \prod _{i=u,v}\left[{\int}_{0}^{\infty}b(M){W}_{i}(M,k)n(M)\phantom{\rule{0.277778em}{0ex}}\mathrm{d}M\right].}\end{array}$$(A.1)
We shall define the integral appearing in Eq. (A.1) to be
$$\begin{array}{c}\hfill {I}_{u}({M}_{1},{M}_{2},k)={\displaystyle {\int}_{{M}_{1}}^{{M}_{2}}b(M){W}_{u}(M,k)n(M)\phantom{\rule{0.277778em}{0ex}}\mathrm{d}M.}\end{array}$$(A.2)
In this appendix we specialise to the twohalo term for the matter–matter spectrum, u = v = m, as this provides a good illustration of the source of confusion. Equation (A.2) should be evaluated over all halo masses, M_{1} → 0 to M_{2} → ∞. Confusion arises because most numerical implementations evaluate this integral only over a finite range of halo mass, which seems sensible given that usually only a finite range of masses are thought to be relevant. However, note that for the matter spectrum we must have the following largescale limit respected on physical grounds:
$$\begin{array}{c}\hfill {P}_{2\mathrm{H},\mathrm{mm}}(k\to 0)={P}_{\mathrm{lin}}(k\to 0).\end{array}$$(A.3)
So the integral in Eq. (A.2) must tend to unity when u = m in the k → 0 limit. In this case, ${W}_{\mathrm{m}}(M,k\to 0)\to M/\overline{\rho}$ and Eq. (A.2) becomes
$$\begin{array}{c}\hfill {I}_{\mathrm{m}}({M}_{1},{M}_{2},k\to 0)=\frac{1}{\overline{\rho}}{\displaystyle {\int}_{{M}_{1}}^{{M}_{2}}b(M)Mn(M)\phantom{\rule{0.277778em}{0ex}}\mathrm{d}M.}\end{array}$$(A.4)
If we pick sensibleseeming limits of M_{1} = 10^{10} h^{−1} M_{⊙} and M_{2} = 10^{16} h^{−1} M_{⊙}, then we find, for a standard cosmological model and a standard mass function at z = 0, that I_{m}(M_{1} = 10^{10} h^{−1} M_{⊙}, M_{2} = 10^{16} h^{−1} M_{⊙}, k → 0)≃0.67. The problem is that, if taken literally, most popular mass functions have a large amount of the total matter contained in very low mass haloes, and so the convergence of Eq. (A.2) is slow as M_{1} is lowered. The upper limit of M_{2} is not usually a problem as long as it is reasonably high, say 10^{16} h^{−1} M_{⊙} at z = 0, which is effectively infinite here.
To address this issue, one thing that should absolutely not be done is to “renormalise” the mass function by multiplying it by a constant to enforce that I_{m}(M_{1}, ∞, k → 0) = 1 for a particular choice of M_{1}. Physically, this amounts to removing lowmass haloes and putting all their mass into highermass haloes. In an example where I_{m}(M_{1}, ∞, k → 0) = 0.67 this would amount to increasing by ∼50% the numbers of haloes between the integration limits! Using a mass function “renormalised” in this way will drastically change the scale dependence in the two and onehalo terms because it increases the number of haloes above M_{1}.
There are two choices for how to proceed: the first (e.g. Yoo et al. 2006; Cacciato et al. 2012) is to assume that ${W}_{\mathrm{m}}(M<{M}_{1},k)=M/\overline{\rho}$, which is equivalent to taking lowmass haloes to have deltafunction profiles. One can then numerically evaluate Eq. (A.2) above the limit of M_{1} and define a new function
$$\begin{array}{c}\hfill A({M}_{1})=1{I}_{\mathrm{m}}({M}_{1},\infty ,k\to 0),\end{array}$$(A.5)
which is the missing part of the integral from below M_{1} (this only needs to be computed once). One can then add A(M_{1}) every time Eq. (A.2) is evaluated:
$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {I}_{\mathrm{m}}({M}_{1},\infty ,k)\to & {\displaystyle {\int}_{{M}_{1}}^{\infty}b(M){W}_{\mathrm{m}}(M,k)n(M)\phantom{\rule{0.277778em}{0ex}}\mathrm{d}M}\hfill \\ \hfill & +A({M}_{1}).\hfill \end{array}\end{array}$$(A.6)
Note that this correction is additive, not multiplicative. This approximation is sufficient as long as the physical shapes of haloes at and below M_{1} do not influence the calculation; however, it is only the appropriate correction for Eq. (A.2) for matter. The second choice is to alter the mass function so that the missing signal from haloes below M_{1} is considered to be in haloes of mass exactly M_{1} (Schmidt 2016). This amounts to making the substitution
$$\begin{array}{c}\hfill n(M)\to {n}^{\prime}(M)=n(M)+\frac{A({M}_{1}){\delta}_{\mathrm{D}}(M{M}_{1})}{b({M}_{1}){M}_{1}/\overline{\rho}},\end{array}$$(A.7)
in Eq. (A.2). This reduces to a different additive correction to Eq. (A.2) of
$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {I}_{u}({M}_{1},\infty ,k)\to & {\displaystyle {\int}_{{M}_{1}}^{\infty}b(M){W}_{u}(M,k)n(M)\phantom{\rule{0.277778em}{0ex}}\mathrm{d}M}\hfill \\ \hfill & +A({M}_{1})\frac{{W}_{u}({M}_{1},k)}{{M}_{1}/\overline{\rho}}\xb7\hfill \end{array}\end{array}$$(A.8)
This approximation is sufficient provided the physical shapes of haloes below M_{1} do not influence the calculation, but it maintains scaledependence in the twohalo term compared to Eq. (A.6). We prefer the second option because it applies evaluations of Eq. (A.2) for any field u(r), not just for matter (note the “general field” subscript u in Eq. (A.8), compared to the “matter” subscript m in Eq. (A.6)). We apply this massfunction alteration only to the twohalo term in all of our calculations.
Note that the problem of slow convergence with respect to the lower integration limit does not usually affect the onehalo term for matter spectra because lowmass haloes contribute very little to the onehalo term, whose amplitude and shape depends on ⟨M⟩ (notation from Sect. 2.3) which is ∼10^{13.5} h^{−1} M_{⊙} at z = 0 for a standard cosmological model, and therefore is usually included in a sensible integration range for Eq. (A.2).
In cases where the power spectrum is unaffected by lowmass haloes then the discussion in this appendix is not relevant. For example, for the galaxy–galaxy number density power you only need to compute the integrals down to the halo mass below which the haloes contain no galaxies (obviously). In the case of the electronpressure field, so little pressure comes from low mass haloes that this correction can be ignored. Considering the discussion in this appendix will, however, be relevant for any spectrum that involves the matter field because the lowmass haloes contribute significantly to the overall mass.
Appendix B: Power spectrum measurements from hydrodynamical simulations
In this paper we measure the power spectrum of overdensity fluctuations and electron pressure using particle data from hydrodynamic simulations. We do this using a FastFourier Transform (FFT) algorithm applied to Cartesian fields, defined on a regular mesh, that are created from particle data. Generally, we consider a meshdefined field u(r) where u_{i} is the contribution to u(r) per particle located at r_{i}
$$\begin{array}{c}\hfill u(\mathit{r})={\displaystyle \sum _{i=1}^{N}{u}_{i}{\delta}_{\mathrm{D}}(\mathit{r}{\mathit{r}}_{i}),}\end{array}$$(B.1)
with N being the total number of particles contributing to u(r).
To create densitycontrast fields from particles in hydrodynamical simulations it is necessary to assign a different density per particle to the mesh because hydrodynamical particles may have different masses m_{i}. We therefore write
$$\begin{array}{c}\hfill \frac{\rho (\mathit{r})}{\overline{\rho}}={\displaystyle \sum _{i=1}^{N}}\frac{{m}_{i}}{\overline{m}}{\delta}_{\mathrm{D}}(\mathit{r}{\mathit{r}}_{i}),\end{array}$$(B.2)
where $\overline{m}$ is the expected mean mass in a cell. In our definitions of overdensity, this division to create a density contrast is always by the mean mass of all matter, not only of the specific species being considered. The situation is more complicated when creating the electron pressure fields from the particles. Hydrodynamic gas particles (labelled with i) come tagged with an internal energy (k_{B}T_{i}) and a local density (ρ_{i}). To convert the particle gas temperature to a contribution to the total electron pressure per particle, P_{e, i}, a quantity that can be summed to create an electronpressure field on mesh via
$$\begin{array}{c}\hfill {P}_{\mathrm{e}}(\mathit{r})={\displaystyle \sum _{i=1}^{N}{P}_{\mathrm{e},i}{\delta}_{\mathrm{D}}(\mathit{r}{\mathit{r}}_{i}),}\end{array}$$(B.3)
we use the ideal gas law:
$$\begin{array}{c}\hfill {P}_{\mathrm{e},i}V={N}_{\mathrm{e},i}{k}_{\mathrm{B}}{T}_{i}\end{array}$$(B.4)
where N_{e, i} is the total number of free electrons in the gas particle in the meshcell volume V:
$$\begin{array}{c}\hfill {N}_{\mathrm{e},i}=\frac{{m}_{i}}{{m}_{\mathrm{p}}}\frac{1}{{\mu}_{\mathrm{e}}},\end{array}$$(B.5)
where m_{i} is the gas particle mass, m_{p} is the proton mass and μ_{e} is the number of free electrons per proton mass. In creating the electronpressure field from the particles we also calculate the particle hydrogen number density
$$\begin{array}{c}\hfill {n}_{\mathrm{H},i}=\frac{{f}_{\mathrm{H}}{\rho}_{i}}{{m}_{\mathrm{p}}},\end{array}$$(B.6)
where f_{H} ≃ 0.752 is the hydrogen gas mass fraction. We exclude particles with n_{H} > 0.1 cm^{−3} as these particles are presumably cold, neutral and should eventually form stars. We checked that imposing this exclusion does not greatly affect our results.
P_{e, i} from Eq. (B.4) is directly related to the quantity Υ_{i} used in Roncarelli et al. (2006, 2007) and McCarthy et al. (2018):
$$\begin{array}{c}\hfill {\mathrm{{\rm Y}}}_{i}={\sigma}_{\mathrm{T}}\frac{{k}_{\mathrm{B}}{T}_{i}}{{m}_{\mathrm{e}}{c}^{2}}\frac{{m}_{i}}{{\mu}_{\mathrm{e}}{m}_{\mathrm{p}}},\end{array}$$(B.7)
where m_{e} is the electron mass and σ_{T} is the Thompson scattering cross section. Υ_{i} differs from P_{e, i} only by a factor of volume and some constants and has dimensions of area rather than pressure.
In Nbody simulations it is well known that the raw measured power spectrum receive an unphysical additive contribution from the finite number of particles – socalled discreteness or “shot” noise. This contribution arises due to the automatic correlation of particles with themselves at zero separation. In the correlation function this contribution is confined to zero separation, but in P(k) the contribution is spread evenly over all wavenumbers by the Fourier transform (the Fourier transform of a delta function is a constant). In a standard Nbody simulation with equalmass particles the shot noise contribution to the matter–matter P(k) spectrum is:
$$\begin{array}{c}\hfill S=\frac{{L}^{3}}{N},\end{array}$$(B.8)
where N is the total number of particles in the simulation and L is the box size. As expected, the shot noise contribution is larger when there are fewer particles. It is common practice to subtract this contribution from power spectra measured from Nbody simulations since in the limit of an infinite number of particles it would vanish.
From hydrodynamical simulations we are interested in calculating auto and crossspectra between between different fields generated from the particles; for example, density, pressure or temperature fields. In the case that the fields to be correlated are generated from distinct sets of particles then the shot noise contribution is zero (e.g. gas–star density). However, in other cases the fields might be generated using (a subset of) same particles (e.g. matter–star density) or from different properties of the same particles (e.g. gas density–electron pressure, which both come from the same hydrodynamical gas particles). There is also the confounding issue that each particle does not provide the same contribution to the eventual field (e.g. individual stellar and gas particles have different masses, gas particles contribute different pressures).
The additive shotnoise contribution for the cross spectrum, P_{uv}(k), between two fields u(r) and v(r), comprised of particles, where a subset $\stackrel{\sim}{N}$ of those particles contribute to both fields u(r) and v(r), can be shown to be
$$\begin{array}{c}\hfill {S}_{\mathit{uv}}=\frac{{L}^{3}}{{M}^{2}}{\displaystyle \sum _{i=1}^{\stackrel{\sim}{N}}{u}_{i}{v}_{i},}\end{array}$$(B.9)
where L is the simulation box size and M is the total number of mesh cells. Note that the shotnoise contribution is not really a function of M, as this will always cancel with factors of M that come from u_{i} and v_{i}. Equation (B.9) can be derived by taking the Fourier transform of Eq. (B.1) and then considering the part of the power spectrum that would still arise even if the particle positions r_{i} are completely uncorrelated, so that the only contribution is from the selfcorrelation of each particle with itself.
Let us demonstrate that Eq. (B.9) reduces to more standard formulae: if we are interested in shot noise in the matter power spectrum we start from the particle masses from which we create a density field. The densitycontrast field generated from the particles is given by Eq. (B.2). We can write $\overline{m}$ as the total mass in the simulation divided by the total number of mesh cells:
$$\begin{array}{c}\hfill \overline{m}=\frac{1}{M}{\displaystyle \sum _{i=1}^{N}{m}_{i}.}\end{array}$$(B.10)
Using Eq. (B.9) and considering shotnoise in the matter–matter power spectrum m we get:
$$\begin{array}{c}\hfill {S}_{\mathrm{mm}}={L}^{3}\frac{{\sum}_{i=1}^{N}{m}_{i}^{2}}{{\left[{\sum}_{i=1}^{N}{m}_{i}\right]}^{2}}\xb7\end{array}$$(B.11)
If all the particle have equal masses then Eq. (B.11) reduces to Eq. (B.8) as required.
In Fig. B.1 we show the contribution of shot noise to the raw autospectra measured for the different fields from the AGN BAHAMAS simulations. Subtracting shot noise is most important for the gas spectrum, which makes sense because gas is the most diffuse component and therefore the one for which discreteness effects are most obvious. At z = 2 at k = 7 h Mpc^{−1} shot noise can makes up 8% of the raw measured gas–gas spectrum. Shot noise is also important at a similar level in some of the crossspectra that we measure. For all of the results presented in this paper shot noise has been subtracted from the power spectra.
Fig. B.1. Shotnoise contributions to autospectra of the various fields measured from the AGN BAHAMAS simulation. We show the ratio of the shotnoise to the autospectra as a function of scale for matter (purple), CDM (green), gas (blue), stars (orange) and electron pressure (yellow). The dashed line shows the cutoff for a one percent contribution. We see that shot noise is most important for the gas power spectrum at higher redshifts, where it can contribute as much as 8% of the raw measured power for the smallest scales shown. Generally shot noise becomes less important as the simulation evolves and structure develops. 
Appendix C: Variations in halomodel parameters
In this appendix we present figures that demonstrate the effect of varying our fitted parameters on the halo profiles and on the eventual power spectra that we are interested in. We show the effect of varying A_{*} (Fig. C.6), M_{*} (Fig. C.7), η (Fig. C.9), ϵ_{1} (Fig. C.1), Γ (Fig. C.5), M_{0} (Fig. C.3), α (Fig. C.10) and T_{whim} (Fig. C.11) on the power spectra, response, halo density profiles and halo mass fractions. These plots are extremely useful if one wants to gain physical insight on the effect of changing the hydrodynamical parameters on the eventual Fourier Space statistics, which can otherwise be quite difficult to understand.
We see that increasing A_{*} increases the amplitude of the stellar content of haloes independently of their mass, and this has an almost scaleindependent effect on any power spectra that includes the star field. There is a small back reaction on the gas because stars can only be added at the expense of gas, but this is very small because stars are a small fraction of the available halo baryons, even in the most extreme cases. Changing M_{*} changes the halo mass at which star formation peaks, which has a more scaledependent effect on the eventual power spectra as stellar mass is shifted into haloes with different structural properties. Altering η changes the distribution of stellar mass between the extended satellite galaxy part and the central galaxy. This has no effect at the largest scales, but effects the onehalo term. Shifting more stellar mass into central galaxies makes this term more and more shotnoise like.
ϵ_{1} governs the concentration modification of lowermass galaxies that have ejected most of their gas. This has no effect on the power spectra at large scales, but changes them at small scales where there is sensitivity to the halo concentration. Γ governs the gas density and pressure profiles and has large effects on the halo profiles. Decreasing Γ makes the gas and pressure profiles more centrally concentrated and therefore boosts the associated power spectra at small scales. Altering M_{0} changes the mass corresponding to haloes that have lost half of their gas (with lowermass haloes having lost more than half). This has a comparatively large effect on spectra, with a strong scale dependent effect on the matter–gas (no change at large scales) and a more scaleindependent effect on spectra involving the pressure.
Increasing α increases the temperature of gas that is bound to the host halo. This boosts the pressure signal in a perfectly scaleindependent way. There is no effect whatsoever on any matter spectra since these are completely insensitive to gas temperature. Conversely, increasing T_{w} boosts the temperature of the ejected gas, which is correlated only on large scales and thus only affects the largescale, twohalo portions of power spectra that involve the pressure.
Fig. C.1. Effect of varying ϵ_{1}, which governs deviations in halo concentration for lowmass haloes that have ejected their gas, on halo profiles (top right two panels), halo mass fractions (bottom right) and power spectra and responses (left column; top left two panels for matter spectra, bottom left two panels for pressure). Increasing ϵ_{1} increases halo concentration and this generally boosts power. From the top right we see this has the largest impact for lower mass haloes because these have ejected most of their gas. From the lefthand panels, this has comparatively more of an impact on power spectra involving gas than on the other spectra. We can also see similar effects on spectra that involve the electron pressure. 
Fig. C.2. Same as Fig. C.1 but for ϵ_{2}, which governs deviations in halo concentration for highmass haloes that have retained their gas. Increasing ϵ_{2} increases halo concentration and this generally boosts power. From the top right we see this has the largest impact for higher mass haloes because these have retained most of their gas. From the lefthand panels, this has comparatively more of an impact on power spectra involving gas than on the other spectra. We can also see some effects on spectra that involve the electron pressure. 
Fig. C.3. Same as Fig. C.1 but for M_{0}, which governs the halo mass below which haloes have lost most of their gas. As can be seen in the bottomright panel, changing M_{0} changes the split of gas between the bound component, which dominates highmass haloes and the unbound component, which is the case for lowmass haloes that have jettisoned most of their gas. The effects of this on the bound halo profiles can be seen in the topright panel. Decreasing M_{0} adds to the amplitude of the gas density profiles and this boosts power in power spectra involving gas. We see a similar, but more limited effect on spectra involving the electron pressure. This limitation arises because lowmass haloes are lower temperature, so contribute much less to spectra involving the electron pressure than they do to those involving the gas density. There is also an overall amplitude affect on the pressure spectra because the unbound gas is taken to have a fixed T_{w} = 10^{6} K that contributes power on large scales. 
Fig. C.4. Same as Fig. C.1 but for β, the powerlaw index governing the decline in halo gas fraction for lowmass haloes. As can be seen from the bottomright panel, increasing β makes the decline in gas fraction sharper, and this has a comparatively large effect on lowmass haloes as can be seen in the other panels on the righthand side of the figure. High mass haloes are affected in the opposite sense because the gas fraction is forced to be 0.5 at 10^{14} h^{−1} M_{⊙}. The effects on the power spectra can be seen in the lefthand panels, where the matter–gas spectrum can be seen to be most affected. 
Fig. C.5. Same as Fig. C.1 but for Γ, the polytropic index for the gas. Increasing Γ makes the gas profile less concentrated, as can be seen in the haloprofile plots (top right). If Γ is lowered then gas is more concentrated and the power spectrum of gas is boosted at small scales. Effects on power spectra that do not involve gas are much smaller. The same trends are seen in the electronpressure profile and spectra involving the electron pressure, which are also sensitive to this parameter as the gas pressure is partly determined by the gas density. 
Fig. C.6. Same as Fig. C.1, but for A_{*}, which governs the peakstarformation efficiency in haloes. Lowerright panel: A_{*} effects the overall abundance of stellar mass, but not the ratio in central and satellite stellar mass. There is some back reaction on the gas abundance that occurs in order to keep the baryon fraction universal, so that halo gas is slightly depleted as a result of this boost, so there is a small effect on gas profiles. The effect on the halo profiles can be seen in the top right panel. Increasing A_{*} boosts starformation efficiency, which has a fairly uniform effect on the amplitude of any power spectra involving the star density as can be seen in the top two panels in the lefthand column. There is an effect on the matter–matter spectrum at small scales as stars start to dominate the power. The effect on electron pressure is minimal, but not nonexistent due to the gas abundance being affected. 
Fig. C.7. Same as Fig. C.1 but for M_{*}, the halo mass at which star formation efficiency peaks. Changing this parameter shifts the halostellarmass fraction to lower and higher masses, as can be seen in the bottomright panel, there is a small backreaction effect on the gas abundance. Changing M_{*} indirectly affects the halo profiles at different masses, as can be seen in the top right panel, but note that the deltafunction component of the stellar halo profile is not shown in this plot. The power spectra involving stars are affects in a more scaledependent manner as can be seen in the topleft two panels. 
Fig. C.8. Same as Fig. C.1 but for σ_{*}, which governs the width of the stellarmass distribution around the peak (which is at M_{*} = 10^{12.5} h^{−1} M_{⊙} here). As can be seen in the bottomleft panel, decreasing σ_{*} shrinks the stellarmass distribution, resulting in fewer stars in both higher and lower mass haloes than M_{*}. The topright panel shows that this has the biggest impact on 10^{14} h^{−1} M_{⊙} haloes, not highermass haloes due to the saturation we impose at high masses. In the lefthand panels we see that the effects of this are constrained to be on power spectra involving the stars. There is only a tiny effect on gas spectra which is due to baryonic mass conservation. 
Fig. C.9. Same as Fig. C.1 but for η, which determines the split of stellar halo mass between central and satellite galaxies for M > M_{*}, where M_{*} = 10^{12.5} h^{−1} M_{⊙} in this figure. In the lowerright panel we see that increasing η from the fiducial −0.3 means that more mass is put into central galaxies, whereas decreasing it places more mass in satellites. As a result of this, the amplitude of the stellar halo is changed, as can be seen in the topright panel, with decreasing η putting more mass into the NFW part of the profile and thus boosting the halo amplitude. This has a scaledependent effect on the matter–stars power spectrum, as can be seen in the topleft two panels. 
Fig. C.10. Same as Fig. C.1 but for α. α governs deviations of the halo temperature from the virial relation. We see that this affects the pressure profiles (middle right) by uniformly scaling the amplitude, because the pressure is directly proportional to temperature; increasing the temperature increases the pressure. We can see the effect of this on the pressure power spectra in the lower left two panels where the effect is not uniform due to the complicated dependence of the power spectra on the underlying halo populations in the two and onehalo terms, which have different scalings with temperature. 
Fig. C.11. Same as Fig. C.1 but for T_{w}. Varying T_{w} only affects the temperature of the unbound gas that contributes to the twohalo term only. We see that as we increase the gas temperature we boost the twohalo term. Note that there is a floor to the effect that lowering T_{w} can have, which is because the onehalo term contributes at large scales in the pressure spectra too, and lowering the temperature of this gas beyond a certain point means that only the twohalo term becomes dwarfed by the onehalo term at large scales. 
Appendix D: Projected fields
In this paper we are primarily interested in the cross correlation between tSZ and various lensing quantities. We use the metric convention of Bartelmann & Schneider (2001):
$$\begin{array}{c}\hfill \mathrm{d}{s}^{2}=\mathrm{d}{t}^{2}{a}^{2}[\mathrm{d}{r}^{2}+{f}_{k}^{2}(r)(\mathrm{d}{\theta}^{2}+{sin}^{2}\theta \mathrm{d}{\varphi}^{2})],\end{array}$$(D.1)
where a is the scale factor, r is the comoving distance, and f_{k} is the comoving, angulardiameter distance.
The weaklensing signal is produced by the bending of light by the largescale structure of the Universe. With some approximations, gravitational lensing is conveniently summarised by the convergence field κ, which can be written as
$$\begin{array}{c}\hfill \kappa (\mathit{\theta})=\frac{3}{2}{\mathrm{\Omega}}_{\mathrm{m}}{\left(\frac{{H}_{0}}{c}\right)}^{2}{\displaystyle {\int}_{0}^{{r}_{\mathrm{H}}}\frac{{f}_{k}(r)}{a(r)}q(r){\delta}_{\mathrm{m}}(r,\mathit{\theta})\phantom{\rule{0.277778em}{0ex}}\mathrm{d}r.}\end{array}$$(D.2)
Here q(r) is the lensingefficiency kernel, which weights redshifts along the lineofsight according to their contribution to the total lensing
$$\begin{array}{c}\hfill q(r)={\displaystyle {\int}_{z(r)}^{\infty}p({z}^{\prime})\frac{{f}_{k}[{r}^{\prime}({z}^{\prime})r]}{{f}_{k}[{r}^{\prime}({z}^{\prime})]}\phantom{\rule{0.277778em}{0ex}}\mathrm{d}{z}^{\prime},}\end{array}$$(D.3)
where p(z) is the normalised redshift distribution of lensed galaxies. If one is interested in CMB lensing then p(z) = δ_{D}(z − z_{*}) where z_{*} is the redshift of the lastscattering surface. In this case:
$$\begin{array}{c}\hfill q(r,{r}_{\ast})=\frac{{f}_{k}({r}_{\ast}r)}{{f}_{k}({r}_{\ast})},\end{array}$$(D.4)
where r_{*} is the comoving distance to z_{*}.
The tSZ signal in the CMB is produced when high energy, free electrons increase the energy of passing, cool CMB photons via inverse Compton scattering. This results in a boost in the effective temperature of the CMB at the location of the electrons that is in proportion to the electron “pressure”, usually expressed in terms of the dimensionless Compton “y” parameter:
$$\begin{array}{c}\hfill y(\mathit{\theta})=\frac{{\sigma}_{\mathrm{T}}}{{m}_{\mathrm{e}}{c}^{2}}{\displaystyle {\int}_{0}^{{r}_{\mathrm{H}}}\frac{{P}_{\mathrm{e}}(r,\mathit{\theta})}{{a}^{3}(r)}\phantom{\rule{0.277778em}{0ex}}a(r)\mathrm{d}r,}\end{array}$$(D.5)
where σ_{T} is the Thompson scattering cross section, m_{e} is the electron mass, c is the speed of light. The integral is taken over comoving distance between us and the particle horizon, r_{H}. The electron pressure can be written as P_{e} = n_{e}k_{B}T_{e} where n_{e} is the comoving electron number density and T_{e} is the physical electron temperature. The factors of a convert our comoving r and electron pressure to a physical quantities.
Given a general 2D field U(θ) that is the projected version of a 3D field u(r, θ) via some projection kernel X_{U}(r):
$$\begin{array}{c}\hfill U(\mathit{\theta})={\displaystyle {\int}_{0}^{{r}_{\mathrm{H}}}{X}_{U}(r)u(r,\mathit{\theta})\phantom{\rule{0.277778em}{0ex}}\mathrm{d}r}\end{array}$$(D.6)
we can write the angular power spectrum between two such projected fields, U and V, in terms of the 3D power spectrum of u and v using the Limber approximation (Kaiser 1992) as
$$\begin{array}{c}\hfill {C}_{\mathit{UV}}(\ell )={\displaystyle {\int}_{0}^{{r}_{\mathrm{H}}}\frac{{X}_{U}(r){X}_{V}(r)}{{f}_{k}^{2}(r)}{P}_{\mathit{uv}}(k(r),z(r))\phantom{\rule{0.277778em}{0ex}}\mathrm{d}r,}\end{array}$$(D.7)
where k(r) = (ℓ + 1/2)/f_{k}(r) and z(r) corresponds to the redshift at comoving distance r. We have increased the lowℓ accuracy of the approximation to 𝒪(ℓ^{−2}) by including the lowestorder correction ℓ → ℓ + 1/2 (Lo Verde et al. 2008). For the lensing projection kernel we have
$$\begin{array}{c}\hfill {X}_{\kappa}(r)=\frac{3}{2}{\mathrm{\Omega}}_{\mathrm{m}}{\left(\frac{{H}_{0}}{c}\right)}^{2}\frac{{f}_{k}(r)q(r)}{a(r)},\end{array}$$(D.8)
and we must project the 3D matter–matter power spectrum, P_{mm}(k). For y projection we have
$$\begin{array}{c}\hfill {X}_{y}(r)=\frac{{\sigma}_{\mathrm{T}}}{{m}_{\mathrm{e}}{c}^{2}}\frac{1}{{a}^{2}(r)},\end{array}$$(D.9)
and we must project the 3D matter–electron pressure power spectrum, P_{mP}(k), to obtain the lensingtSZ cross correlation.
There has been some recent discussion on the validity of the Limber approximation for the interpretation of cosmological observations. We note that comparisons between this approximation and a full calculation show that the Limber approximation is valid given the scales we are interested in, and that we will be interested in in the near future. We point the interested reader towards Hill & Pajer (2013), Kilbinger et al. (2017) and Lemos et al. (2017). Note that the accuracy of the Limber approximation is generally better for C_{κy}(ℓ) compared to for C_{κκ}(ℓ) due to the broader projection kernel for y compared to that of κ.
The contributions to C_{κκ}(ℓ) and C_{κy}(ℓ) as a function of k, z and r are shown in Fig. D.1 for different ℓ. Note that, for a given ℓ, C_{κy}(ℓ) receives its contributions from lower z and r compared to C_{κκ}(ℓ), and also from a much broader range of k, z and r. To generate this figure we used a p(z) distribution (Eq. (D.3)) taken from KiDS galaxies for redshifts z = 0.1 → 0.9 (Joudaki et al. 2017b).
Fig. D.1. Normalised contribution to each ℓ of C(ℓ) for κ–κ (top) and κ–y (bottom) as a function of scale k, redshift z and (comoving) distance r. The lensing is taken from KiDS galaxies for redshifts z = 0.1 → 0.9. Each of these curves is the normalised integrand in Eq. (D.7) reexpressed as a function of either k, z or r (so dC(ℓ)/dk, dC(ℓ)/dz, dC(ℓ)/dr). In practice, we carry out our integration in r between the observer and z(r) = 9 and we checked that our integrations are robust to this upper limit. The power spectra being integrated are set to zero for k < 10^{−4} h Mpc^{−1} and k > 10^{2} h Mpc^{−1} and again our results are robust to this. Note the generally broader and lower–z contributions for κ–y compared to those for κ–κ. 
All Tables
Bestfitting effective halo model parameters for different combinations of power spectra from the BAHAMAS simulations.
All Figures
Fig. 1. Halo mass fractions f_{u}(M) as a function of halo mass for the different components in our model at z = 0: The halo CDM fraction is constant at Ω_{c}/Ω_{m} for all halo masses. Star formation efficiency peaks around M_{*} = 10^{12.5} h^{−1} M_{⊙} with a mass fraction A_{*} = 0.03, it saturates at high mass at a value of f_{*} = A_{*}/3. Stars are split into central and satellite, with the halo stellar content being dominated by centrals at low masses and by satellites at high mass. The bound gas fraction is near to the universal baryon fraction for high halo masses, but decreases below half the universal fraction for halo masses lower than M_{0} = 10^{14} h^{−1} M_{⊙}, the remaining baryons that are neither bound gas or stars are considered to be in unbound gas that is not in the halo. The total mass fraction drops from unity at lower masses due to the ejected gas. 

In the text 
Fig. 2. Halomodel predictions at z = 0 for the matter–matter (top), CDM–CDM, gas–gas, stars–stars and the matter–electron pressure (bottom) power spectra. We also show the twohalo (longdashed) and the onehalo (shortdashed) terms for each spectrum apart from matter–matter. At large scales, the shapes of the twohalo terms are all identical, and are that of the linear power spectrum, at small scales the onehalo terms dominates and have significantly different shapes for each spectrum. Note that the transition scale between the two and onehalo term is at much larger scales for the matter–electron pressure cross spectrum, which is a consequence of this spectrum being dominated by contributions from more massive haloes compared to those that contribute to the matter spectra. The downturn at small scales in the gas–gas and matter–electron pressure spectra is partly due to the underrepresentation of lowmass haloes and partly due to their relatively smooth halo profiles. 

In the text 
Fig. 3. Cumulative halo model power spectra as a function of the upperlimit of halo mass at z = 0. Lefthand panel: matter–matter power spectrum as the upper limit is raised from 10^{10} to 10^{16} h^{−1} M_{⊙} (i.e. the lightest coloured curve shows the contribution from only haloes below 10^{10} h^{−1} M_{⊙}), righthand panel: same but for the matter–electron pressure spectrum. In both cases the converged spectrum from all halo masses is shown in thick black and this has been reached with an upper limit of 10^{16} h^{−1} M_{⊙}. We see that the matter–matter spectrum builds up fairly gradually with halo mass and that all scales receive contributions from a wide range of halo masses. In contrast, the matter–electron pressure spectrum receives very little contribution from haloes less massive than 10^{13} h^{−1} M_{⊙}, with the vast majority of the power coming from haloes between 10^{14} and 10^{15} h^{−1} M_{⊙}. 

In the text 
Fig. 4. Effect of σ_{8} variations on the matter–matter power spectrum (left) and on the matter–electron pressure power spectrum (right) at z = 0. We show halo model power spectra for σ_{8} = 0.7 (blue) to σ_{8} = 0.9 (red) divided by a central σ_{8} = 0.8 model. For linear theory the plot would simply be horizontal lines. We see that boosting the amplitude of fluctuations increases the amplitude of the matter–electron power spectrum far more than it does for the mattermatter. At a wavenumber of ≃1 h Mpc^{−1} the mattermatter power spectrum scales like ${\sigma}_{8}^{3.5}$ while matter–electron pressure scales like ${\sigma}_{8}^{5.8}$. This is because the electron pressure field is dominated by highmass haloes from the tail of the mass function that are very sensitive to the linear power spectrum amplitude. 

In the text 
Fig. 5. Different fields measured from the BAHAMAS AGN simulation at z = 0 around a M ∼ 10^{15} h^{−1} M_{⊙} halo averaged through a square slab of side 20 h^{−1} Mpc and thickness 5 h^{−1} Mpc. We show CDM (top left), gas (top right) and star (bottom left) density contrast as well as electron pressure (bottom right). Colour intensity increases with field value logarithmically and we show a dynamic range of 10^{3} for both density contrast and electron pressure. We see that the gas is less clumped than the CDM and that lowmass sub haloes show a large gas depletion compared to the host. Stars are more tightly clustered in the halo centre than the CDM, but are absent in many of the substructures. The total matter densitycontrast field would be a sum of the CDM, gas and star fields. The electron pressure broadly follows the gas but can only be seen emanating from the highest gasdensity peaks over the dynamic range shown. 

In the text 
Fig. 6. Comparison of the default halomodel prediction with power spectra from the BAHAMAS AGN simulation at z = 0. Points with errors show the measured power spectra with an erroronthemean calculated from the finite number of modes that contribute to each k. The upper two panels show spectra for “all matter” (purple), CDM (green), gas (light blue), stars (orange) and electron pressure (dark blue). The electron pressure field has units of 100 eV cm^{−3}. The effect of power aliasing can be seen as an upturn in power at the highest k (∼7 h Mpc^{−1}) shown for each spectrum (see Appendix B). Lines show the halomodel predictions using the default model discussed in Sect. 3: solid lines and points show autospectra while broken lines and open points show crossspectra with the total matter field. We see that the trends in the halo model agree reasonably well with the simulations but that there are disagreements in the details. In the lower two panels we show the “response function” calculated with respect to the matter–matter power in the gravityonly model: ${P}_{\mathit{uv}}^{\mathrm{hydro}}(k)/{P}_{\mathrm{mm}}^{\mathrm{gravity}}(k)$. The simulations spectra have been divided by those from the DMONLY simulation and the halo model has been divided by a standard halomodel prediction that assumes all matter to be in NFW haloes. The horizontaldashedblack lines in the lower left panel show the expected large scale Ω_{c}/Ω_{m}, (Ω_{c}/Ω_{m})^{2}, Ω_{b}/Ω_{m}, and (Ω_{b}/Ω_{m})^{2} asymptotes for the CDM and gas cross–matter and autospectra. We see that the response functions are generally smoother than the raw power spectra: a consequence of cosmic variance cancellations at large scales and cancellation of aliasing effects at small scales. 

In the text 
Fig. 7. Halo model (1) for the star–star power spectrum is shown as sold lines while simulation measurements are shown as points with error bars that come from the scatter in power between three different realisations of the AGN BAHAMAS simulation. Note that the power spectra shown here is the measured response multiplied by HMCODE, so cosmic variance is removed at small scales, but we still show the error bar for comparison. Different colours denote the three different AGN feedback temperatures: 10^{7.6} (blue), 10^{7.8} (grey) and 10^{8.0} K (red). Parameters for the model can be found in Table 2. We see that the spectrum is reduced in amplitude as feedback temperature increases, a consequence of AGN feedback suppressing star formation. This is reflected in the modelling as the value of A_{*} and M_{*} decreases as AGN strength increases. 

In the text 
Fig. 8. Similar to Fig. 7, but for the matter–matter spectrum (2). The power spectra are shown in the top row, while the response functions are shown in the bottom row. We see that the model is nearly perfect (errors at the subpercent level) across a wide range of scales and for the different feedback scenarios. We see that the amount of powerspectrum suppression increases with the feedback temperature. The critical parameter that governs this in our modelling is the halo mass M_{0}, which governs the halo mass below which more than half of the cosmic baryon density is missing from the halo. In the fitting parameters we see a strong preference for higher values of M_{0} for stronger feedback temperatures (from blue to red). How our model fares outside the range of AGN heating temperatures to which it was fitted is also shown via the lines that have no corresponding simulation points. To generate these we have interpolated and extrapolated our model parameters as a function of AGN heating temperature, which is shown in 0.1 dex intervals between 10^{7.5} (bluest) and 10^{8.1} K (reddest), indicating that our model is robust to extrapolation. 

In the text 
Fig. 9. Same as Figs. 7 and 8, but for model (3) for the matter–matter (first and third row; power and response) and matter–electron pressure spectra (second and fourth row; power and response). The match to the matter–matter spectra is slightly degraded compared to the case when we fitted the matter–matter spectra exclusively, but we still match the response at the few percent level while simultaneously we match the matter–electron pressure spectrum well for k > 1 h Mpc^{−1}. If we compare this to the error bars in the power spectra in the second row we see that our fit is good for scales where the three AGN heating temperature strengths are clearly demarcated in k. At larger scales we suspect the larger errors arise because the power is predominantly coming from few massive haloes. Similar to Fig. 8, we show how the model fares outside the range of AGN strengths to which is was fitted. AGN heating temperature is shown in 0.1 dex intervals between 10^{7.5} K and 10^{8.1} K, while the model was fitted between 10^{7.6} K and 10^{8.0} K, thus demonstrating that our model is robust to extrapolation. 

In the text 
Fig. 10. Same as Figs. 7–9, but for the bestfitting model (4) to all 10 possible power spectra of matter, CDM, gas and star fields, although we only show the 4 autospectra here. The top half of the figure show power spectra while the bottom half show response functions. Columns show different redshifts. Different colours are different feedback temperatures: T = 10^{7.6} (blue), T = 10^{7.8} (grey) and T = 10^{8.0} K (red). Points (with errors) are from BAHAMAS while lines are from our halo model. 

In the text 
Fig. 11. Full error matrix for our fitted models calculated for each field pair for the 10^{7.8} K AGN simulation. The error is averaged linearly over redshift between z = 0 and 1 and over k logarithmically between ≃0.015 and 7 h Mpc^{−1}. This matrix is symmetric because power spectra are symmetric when swapping the field labels. We show the matrix for the three models presented in this work: matter (2, left); matter, electron pressure (3, centre); matter, CDM, gas, stars (4, right). The squares that contain dots are the particular crossspectra that are fitted in the model. In the lefthand plot, we see that in the matter model we get a reasonable fit to the power spectra all of the constituent fields. In the central plot we get a better match to to the matter–electron pressure spectrum than in the lefthand panel, but this is at the expense of any power spectra that directly involves the gas field, which shows some tension in our model between the gas and pressure modelling. In the righthand panel, where we fitted matter, CDM, gas and stars, we see that all the matter spectra are reasonable, and better than in the lefthand panel, but achieving this is at the expense of any spectra that involves the pressure field. 

In the text 
Fig. B.1. Shotnoise contributions to autospectra of the various fields measured from the AGN BAHAMAS simulation. We show the ratio of the shotnoise to the autospectra as a function of scale for matter (purple), CDM (green), gas (blue), stars (orange) and electron pressure (yellow). The dashed line shows the cutoff for a one percent contribution. We see that shot noise is most important for the gas power spectrum at higher redshifts, where it can contribute as much as 8% of the raw measured power for the smallest scales shown. Generally shot noise becomes less important as the simulation evolves and structure develops. 

In the text 
Fig. C.1. Effect of varying ϵ_{1}, which governs deviations in halo concentration for lowmass haloes that have ejected their gas, on halo profiles (top right two panels), halo mass fractions (bottom right) and power spectra and responses (left column; top left two panels for matter spectra, bottom left two panels for pressure). Increasing ϵ_{1} increases halo concentration and this generally boosts power. From the top right we see this has the largest impact for lower mass haloes because these have ejected most of their gas. From the lefthand panels, this has comparatively more of an impact on power spectra involving gas than on the other spectra. We can also see similar effects on spectra that involve the electron pressure. 

In the text 
Fig. C.2. Same as Fig. C.1 but for ϵ_{2}, which governs deviations in halo concentration for highmass haloes that have retained their gas. Increasing ϵ_{2} increases halo concentration and this generally boosts power. From the top right we see this has the largest impact for higher mass haloes because these have retained most of their gas. From the lefthand panels, this has comparatively more of an impact on power spectra involving gas than on the other spectra. We can also see some effects on spectra that involve the electron pressure. 

In the text 
Fig. C.3. Same as Fig. C.1 but for M_{0}, which governs the halo mass below which haloes have lost most of their gas. As can be seen in the bottomright panel, changing M_{0} changes the split of gas between the bound component, which dominates highmass haloes and the unbound component, which is the case for lowmass haloes that have jettisoned most of their gas. The effects of this on the bound halo profiles can be seen in the topright panel. Decreasing M_{0} adds to the amplitude of the gas density profiles and this boosts power in power spectra involving gas. We see a similar, but more limited effect on spectra involving the electron pressure. This limitation arises because lowmass haloes are lower temperature, so contribute much less to spectra involving the electron pressure than they do to those involving the gas density. There is also an overall amplitude affect on the pressure spectra because the unbound gas is taken to have a fixed T_{w} = 10^{6} K that contributes power on large scales. 

In the text 
Fig. C.4. Same as Fig. C.1 but for β, the powerlaw index governing the decline in halo gas fraction for lowmass haloes. As can be seen from the bottomright panel, increasing β makes the decline in gas fraction sharper, and this has a comparatively large effect on lowmass haloes as can be seen in the other panels on the righthand side of the figure. High mass haloes are affected in the opposite sense because the gas fraction is forced to be 0.5 at 10^{14} h^{−1} M_{⊙}. The effects on the power spectra can be seen in the lefthand panels, where the matter–gas spectrum can be seen to be most affected. 

In the text 
Fig. C.5. Same as Fig. C.1 but for Γ, the polytropic index for the gas. Increasing Γ makes the gas profile less concentrated, as can be seen in the haloprofile plots (top right). If Γ is lowered then gas is more concentrated and the power spectrum of gas is boosted at small scales. Effects on power spectra that do not involve gas are much smaller. The same trends are seen in the electronpressure profile and spectra involving the electron pressure, which are also sensitive to this parameter as the gas pressure is partly determined by the gas density. 

In the text 
Fig. C.6. Same as Fig. C.1, but for A_{*}, which governs the peakstarformation efficiency in haloes. Lowerright panel: A_{*} effects the overall abundance of stellar mass, but not the ratio in central and satellite stellar mass. There is some back reaction on the gas abundance that occurs in order to keep the baryon fraction universal, so that halo gas is slightly depleted as a result of this boost, so there is a small effect on gas profiles. The effect on the halo profiles can be seen in the top right panel. Increasing A_{*} boosts starformation efficiency, which has a fairly uniform effect on the amplitude of any power spectra involving the star density as can be seen in the top two panels in the lefthand column. There is an effect on the matter–matter spectrum at small scales as stars start to dominate the power. The effect on electron pressure is minimal, but not nonexistent due to the gas abundance being affected. 

In the text 
Fig. C.7. Same as Fig. C.1 but for M_{*}, the halo mass at which star formation efficiency peaks. Changing this parameter shifts the halostellarmass fraction to lower and higher masses, as can be seen in the bottomright panel, there is a small backreaction effect on the gas abundance. Changing M_{*} indirectly affects the halo profiles at different masses, as can be seen in the top right panel, but note that the deltafunction component of the stellar halo profile is not shown in this plot. The power spectra involving stars are affects in a more scaledependent manner as can be seen in the topleft two panels. 

In the text 
Fig. C.8. Same as Fig. C.1 but for σ_{*}, which governs the width of the stellarmass distribution around the peak (which is at M_{*} = 10^{12.5} h^{−1} M_{⊙} here). As can be seen in the bottomleft panel, decreasing σ_{*} shrinks the stellarmass distribution, resulting in fewer stars in both higher and lower mass haloes than M_{*}. The topright panel shows that this has the biggest impact on 10^{14} h^{−1} M_{⊙} haloes, not highermass haloes due to the saturation we impose at high masses. In the lefthand panels we see that the effects of this are constrained to be on power spectra involving the stars. There is only a tiny effect on gas spectra which is due to baryonic mass conservation. 

In the text 
Fig. C.9. Same as Fig. C.1 but for η, which determines the split of stellar halo mass between central and satellite galaxies for M > M_{*}, where M_{*} = 10^{12.5} h^{−1} M_{⊙} in this figure. In the lowerright panel we see that increasing η from the fiducial −0.3 means that more mass is put into central galaxies, whereas decreasing it places more mass in satellites. As a result of this, the amplitude of the stellar halo is changed, as can be seen in the topright panel, with decreasing η putting more mass into the NFW part of the profile and thus boosting the halo amplitude. This has a scaledependent effect on the matter–stars power spectrum, as can be seen in the topleft two panels. 

In the text 
Fig. C.10. Same as Fig. C.1 but for α. α governs deviations of the halo temperature from the virial relation. We see that this affects the pressure profiles (middle right) by uniformly scaling the amplitude, because the pressure is directly proportional to temperature; increasing the temperature increases the pressure. We can see the effect of this on the pressure power spectra in the lower left two panels where the effect is not uniform due to the complicated dependence of the power spectra on the underlying halo populations in the two and onehalo terms, which have different scalings with temperature. 

In the text 
Fig. C.11. Same as Fig. C.1 but for T_{w}. Varying T_{w} only affects the temperature of the unbound gas that contributes to the twohalo term only. We see that as we increase the gas temperature we boost the twohalo term. Note that there is a floor to the effect that lowering T_{w} can have, which is because the onehalo term contributes at large scales in the pressure spectra too, and lowering the temperature of this gas beyond a certain point means that only the twohalo term becomes dwarfed by the onehalo term at large scales. 

In the text 
Fig. D.1. Normalised contribution to each ℓ of C(ℓ) for κ–κ (top) and κ–y (bottom) as a function of scale k, redshift z and (comoving) distance r. The lensing is taken from KiDS galaxies for redshifts z = 0.1 → 0.9. Each of these curves is the normalised integrand in Eq. (D.7) reexpressed as a function of either k, z or r (so dC(ℓ)/dk, dC(ℓ)/dz, dC(ℓ)/dr). In practice, we carry out our integration in r between the observer and z(r) = 9 and we checked that our integrations are robust to this upper limit. The power spectra being integrated are set to zero for k < 10^{−4} h Mpc^{−1} and k > 10^{2} h Mpc^{−1} and again our results are robust to this. Note the generally broader and lower–z contributions for κ–y compared to those for κ–κ. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.