Towards determining the neutrino mass hierarchy: weak lensing and galaxy clustering forecasts with baryons and intrinsic alignments

TL;DR

This paper evaluates how baryonic physics and intrinsic alignments affect the ability of future surveys to measure neutrino masses and distinguish their hierarchy, finding that combined data can reduce uncertainties but hierarchy distinction remains challenging.

Contribution

It provides the first combined weak lensing and galaxy clustering Fisher analysis including baryons, intrinsic alignments, and neutrinos for hierarchy determination.

Findings

Baryons significantly degrade neutrino mass forecasts, doubling uncertainties.

Combining galaxy clustering and CMB priors reduces forecast uncertainties to about 9-16%.

Forecasts suggest hierarchy distinction remains difficult with current methods.

Abstract

The capacity of Stage IV lensing surveys to measure the neutrino mass sum and differentiate between the normal and inverted mass hierarchies depends on the impact of nuisance parameters describing small-scale baryonic astrophysics and intrinsic alignments. For a Euclid-like survey, we perform the first combined weak lensing and galaxy clustering Fisher analysis with baryons, intrinsic alignments, and massive neutrinos for both hierarchies. We use a matter power spectrum generated from a halo model that captures the impact of baryonic feedback and adiabatic contraction. For weak lensing, we find that baryons cause severe degradation to forecasts of the neutrino mass sum, $\Sigma$, approximately doubling $\sigma_{\Sigma}$. We show that including galaxy clustering constraints from Euclid and BOSS, and cosmic microwave background (CMB) Planck priors, can reduce this degradation to $\sigma_{\Sigma}$ to 9% and 16% for the normal and inverted hierarchies respectively. The combined forecasts, $\sigma_{\Sigma_{\rm{NH}}}=0.034\, \rm{eV}$ and $\sigma_{\Sigma_{\rm{IH}}}=0.034\, \rm{eV}$, preclude a meaningful distinction of the hierarchies but could be improved upon with future CMB priors on $n_s$ and information from neutrinoless double beta decay to achieve a 2$\sigma$ distinction. The effect of intrinsic alignments on forecasts is shown to be minimal, with $\sigma_{\Sigma}$ even experiencing mild improvements due to information from the intrinsic alignment signal. We find that while adiabatic contraction and intrinsic alignments will require careful calibration to prevent significant biasing of $\Sigma$, there is less risk presented by feedback from energetic events like AGN and supernovae.

Tables

Table 1. Table 1: Survey parameters characterizing a Euclid-like space mission. These include the area of sky covered, A sky subscript 𝐴 sky A_{\rm{sky}} , the survey redshift limits, the median redshift value, z med subscript 𝑧 med z_{\rm{med}} , the number of bins, N z , WL subscript 𝑁 𝑧 WL N_{z,\rm{WL}} , the number density of surveyed galaxies, n gal subscript 𝑛 gal n_{\rm{gal}} , the photometric redshift error, the intrinsic ellipticity dispersion and the range of accessible wavenumbers.

Parameter	Euclid-like value
$A_{\mathrm{sky}}$	$15,000{\mathrm{deg}}^2$
$z_{\min }$	$0.$
$z_{\max }$	$2.0$
$z_{\mathrm{med}}$	$0.9$
$N_{z,\mathrm{WL}}$	$10$
$n_{\mathrm{gal}}$	$30$ gal/arcmin2
${\sigma }_z$	$0.05$
${\sigma }_e$	$0.3$
$l_{\min }$	$10$
$l_{\max }$	$5000$

Table 2. Table 2: 1- σ 𝜎 \sigma error forecasts for a Euclid-like survey of the neutrino mass sum in the normal and inverted hierarchies, without and including marginalization over baryon parameters and the addition of priors on the Λ Λ \Lambda CDM cosmological parameters from Planck CMB measurements. We also include the response factors, R B subscript 𝑅 𝐵 R_{B} , to including baryon marginalization.

	${\sigma }_{\mathrm{WL}}/\mathrm{eV}$	${\sigma }_{\mathrm{WL}}/\mathrm{eV}$	$R_{B,\mathrm{WL}}$	${\sigma }_{\mathrm{WL}+\mathrm{CMB}}/\mathrm{eV}$	${\sigma }_{\mathrm{WL}+\mathrm{CMB}}/\mathrm{eV}$	$R_{B,\mathrm{WL}+\mathrm{CMB}}$
	(ex bar.)	(inc. bar.)		(ex bar.)	(inc. bar.)
${\Sigma }_{\mathrm{NH}}$	0.079	0.156	1.97	0.036	0.044	1.21
${\Sigma }_{\mathrm{IH}}$	0.120	0.204	1.70	0.033	0.042	1.28

Table 3. Table 3: 1- σ 𝜎 \sigma error forecasts for a Euclid-like survey of the neutrino mass sum in the normal and inverted hierarchies, without and including marginalization over baryon parameters and the addition of priors from galaxy clustering and Planck CMB measurements. We also include the response factors, R B subscript 𝑅 𝐵 R_{B} , to including baryon marginalization.

	${\sigma }_{\mathrm{WL}}/\mathrm{eV}$	${\sigma }_{\mathrm{WL}}/\mathrm{eV}$	$R_{B,\mathrm{WL}}$	${\sigma }_{\mathrm{WL}+\mathrm{GC}}/\mathrm{eV}$	${\sigma }_{\mathrm{WL}+\mathrm{GC}}/\mathrm{eV}$	$R_{B,\mathrm{WL}+\mathrm{GC}}$	${\sigma }_{\mathrm{WL}+\mathrm{GC}+\mathrm{CMB}}/\mathrm{eV}$	${\sigma }_{\mathrm{WL}+\mathrm{GC}+\mathrm{CMB}}/\mathrm{eV}$	$R_{B,\mathrm{WL}+\mathrm{GC}+\mathrm{CMB}}$
	(ex bar.)	(inc. bar.)		(ex bar.)	(inc. bar.)		(ex bar.)	(inc. bar.)
${\Sigma }_{\mathrm{NH}}$	0.079	0.156	1.97	0.045	0.053	1.19	0.031	0.034	1.09
${\Sigma }_{\mathrm{IH}}$	0.120	0.204	1.70	0.050	0.078	1.57	0.030	0.034	1.16

Equations

$${\delta }_c\propto 1+0.262f_{\nu }$$

$${\Delta }_v\propto 1+0.916f_{\nu },$$

$$f_{\nu }=\frac{{\Omega }_{\nu }}{{\Omega }_m}\simeq \frac{\Sigma }{94.1{\Omega }_m{h}^2\mathrm{eV}}.$$

$$u\left(k|M\right)⟶u\left({\nu }^{\eta }k|M\right),$$

$$\rho \left(r,M\right)=\frac{{\rho }_s}{{\nu }^{3\eta }}\frac{1}{\left(\frac{r}{{\nu }^{\eta }{r}_s}\right){\left[1+\left(\frac{r}{{\nu }^{\eta }{r}_s}\right)\right]}^2},$$

$$\eta ={\eta }_0-0.3{\sigma }_8\left(z\right),$$

$$F_{\mathrm{tot}}=F_{\mathrm{WL}}+F_{\mathrm{CMB}}.$$

$$F_{ii}^{\prime }⟶F_{ii}\left(1+\alpha \right).$$

$$\delta {\theta }_i=-{\left[F^{\theta \theta }\right]}_{ik}^{-1}{F}_{kj}^{\theta \psi }\delta {\psi }_j,$$

$${\Delta }^2\left(k\right)={\Delta }_{1h}^2\left(k\right)+{\Delta }_{2h}^2\left(k\right).$$

$${\Delta }_{1h}^2\left(k\right)=\frac{k^3}{2{\pi }^2}{\int }_0^{\infty }\mathrm{d}M\frac{M^{2}n\left(M\right)}{{\bar{\rho }}^2}{u}^2\left(k|M\right).$$

$$f\left(\nu \right)=\frac{M}{\bar{\rho }}n\left(M\right)\frac{\mathrm{d}M}{\mathrm{d}\nu }$$

$$f\left(\nu \right)=0.2162\left[1+{\left(0.707{\nu }^2\right)}^{-0.3}\right]\exp \left[-\frac{0.707{\nu }^2}{2}\right].$$

$${\sigma }^2\left(R\right)={\int }_0^{\infty }\mathrm{d}\ln k{\Delta }_{\mathrm{lin}}^2\left(k\right)W^2\left(kR\right),$$

$$u\left(k|M\right)=\frac{4\pi }{M}{\int }_0^{r_v}{r}^2\mathrm{d}r\frac{\sin \left(kr\right)}{kr}\rho \left(r,M\right),$$

$$\rho \left(r,M\right)=\frac{{\rho }_s}{\left(\frac{r}{r_s}\right){\left[1+\left(\frac{r}{r_s}\right)\right]}^2}.$$

$$r_v={\left(\frac{3M}{4\pi \bar{\rho }{\Delta }_v}\right)}^{\frac{1}{3}}.$$

$$c\left(M,z\right)=A_B\frac{1+z_f}{1+z},$$

$$n\left(z\right)\propto z^2\exp \left[-{\left(\frac{z}{z_0}\right)}^{\frac{3}{2}}\right],$$

$$C_{\ell ,ij}=\frac{9}{4}{\Omega }_m^2{\left(\frac{H_0}{c}\right)}^4{\int }_0^{{\chi }_{\max }}\mathrm{d}\chi \frac{g_i\left(\chi \right)g_j\left(\chi \right)}{a^2\left(\chi \right)}P\left(k=\frac{\ell }{f_K\left(\chi \right)},\chi \right),$$

$$g_i\left(\chi \right)={\int }_{\chi }^{{\chi }_{\max }}\mathrm{d}{\chi }^{\prime }{n}_i\left({\chi }^{\prime }\right)\frac{f_K\left(\chi -{\chi }^{\prime }\right)}{f_K\left({\chi }^{\prime }\right)}.$$

$$n_i\left(z\right)=\frac{n\left(z\right){\int }_{z_{i,-}}^{z_{i,+}}\mathrm{d}{z}_{ph}{p}_{ph}\left(z_{ph}|z\right)}{{\int }_{z_{\min }}^{z_{\max }}\mathrm{d}{z}^{\prime }n\left(z^{\prime }\right){\int }_{z_{i,-}}^{z_{i,+}}\mathrm{d}{z}_{ph}{p}_{ph}\left(z_{ph}|z^{\prime }\right)},$$

$$p_{ph}\left(z_{ph}|z\right)=\frac{1}{\sqrt{2\pi }{\sigma }_z\left(1+z\right)}\exp \left\{-{\left[\frac{z-z_{ph}}{\sqrt{2}{\sigma }_z\left(1+z\right)}\right]}^2\right\}.$$

$$C_{\ell ,ij}=C_{\ell ,ij}^{\gamma \gamma }+C_{\ell ,ij}^{I\gamma }+C_{\ell ,ij}^{II},$$

$$C_{\ell ,ij}^{I\gamma }=\frac{3}{2}{\Omega }_m{\left(\frac{H_0}{c}\right)}^2{\int }_0^{{\chi }_{\max }}\mathrm{d}\chi \frac{g_i\left(\chi \right)n_j\left(\chi \right)+g_j\left(\chi \right)n_i\left(\chi \right)}{a\left(\chi \right)\chi }$$

$$\times P_{mI}\left(k=\frac{\ell }{f_K\left(\chi \right)},\chi \right),$$

$$C_{\ell ,ij}^{II}={\int }_0^{{\chi }_{\max }}\mathrm{d}\chi \frac{n_i\left(\chi \right)n_j\left(\chi \right)}{{\chi }^2}{P}_{II}\left(k=\frac{\ell }{f_K\left(\chi \right)},\chi \right).$$

$$P_{mI}\left(k,\chi \right)=-a_{\mathrm{IA}}\frac{c_{\mathrm{IA}}{\rho }_{\mathrm{crit}}{\Omega }_m}{D\left(z\left(\chi \right)\right)}{P}_m\left(k,\chi \right)$$

$$P_{II}\left(k,\chi \right)={\left(a_{\mathrm{IA}}\frac{c_{\mathrm{IA}}{\rho }_{\mathrm{crit}}{\Omega }_m}{D\left(z\left(\chi \right)\right)}\right)}^2{P}_m\left(k,\chi \right).$$

$$P_{gg}\left(k,\mu ;z\right)={\left(b{\sigma }_8\left(z\right)+f{\sigma }_8\left(z\right){\mu }^2\right)}^2\frac{P_m\left(k;z\right)}{{\sigma }_8^2\left(z\right)},$$

$${\alpha }_{\perp }\left(z\right)\equiv \frac{D_A\left(z\right)}{D_{A,\mathrm{fid}}\left(z\right)}$$

Full text

††pubyear: 2019††pagerange: Towards determining the neutrino mass hierarchy: weak lensing and galaxy clustering forecasts with baryons and intrinsic alignments–D

1 Introduction

Neutrino cosmology is an intersection of particle physics and astrophysics that has the potential to greatly enrich both fields. Flavour oscillations measured in solar, atmospheric, reactor and accelerator experiments have constrained the mass-squared differences, $\Delta m_{12}^{2}$ and $|\Delta m_{23}^{2}|$, for the three neutrino mass eigenstates but cannot determine the absolute masses (e.g., Maltoni et al. 2004; Fogli et al. 2006). The sign of the largest splitting defines two possible orderings, known as the normal (NH, $\Delta m_{23}^{2}>0$) and inverted (IH, $\Delta m_{23}^{2}<0$) hierarchies. Lower bounds on the sum of the masses for each hierarchy, $\Sigma_{\rm{NH}}\gtrsim 0.06\,\rm{eV}$ and $\Sigma_{\rm{IH}}\gtrsim 0.1\,\rm{eV}$, can be arrived at numerically from the oscillation constraints (Lesgourgues & Pastor 2006). Cosmology has been providing increasingly tight upper bounds. Data from the Lyman-$\alpha$ forest in combination with Planck Collaboration et al. (2016) cosmic microwave background (CMB) constraints (Palanque-Delabrouille et al. 2015; Yèche et al. 2017) and the more recent works by Vagnozzi et al. (2017) and Planck Collaboration et al. (2018) analysing CMB and baryon acoustic oscillation (BAO) data all find $\Sigma\lesssim 0.12\,\rm{eV}$ (95% confidence level). Though this starts to confront the limits of the IH, it remains an open question as to whether future cosmological data will be able to distinguish between the hierarchies by constraining the mass sum. The degeneracy between hierarchies is partially broken by different sets of free-streaming scales associated with the mass eigenstates, but these differences are sufficiently small for a fixed mass sum that there is no realistic prospect of achieving a detection in the near future (Hall & Challinor 2012).

Forthcoming Stage IV surveys, such as Euclid (Laureijs et al. 2011) and LSST (LSST Science Collaboration et al. 2009), are aiming to achieve $\Sigma$ errors small enough to determine the hierarchy by using large-scale structure (LSS) probes like tomographic weak gravitational lensing (Bartelmann & Schneider 2001) and spectroscopic galaxy clustering. The former uses multiple photometric redshift bins to extract the influence of cosmological quantities on the expansion history and the growth of the dark matter distribution (Albrecht et al. 2006; Peacock et al. 2006). Understanding how galaxies trace the underlying matter distribution through clustering measurements provides an additional repository of information on the properties of the Universe (Hauser & Peebles 1973; Davis & Peebles 1983). This is subject to a range of systematics and statistical effects including redshift uncertainty (see e.g., Wang et al. 2013), galaxy bias (e.g., Desjacques et al. 2018; Vagnozzi et al. 2018) and anisotropies due to the Alcock-Paczynski (AP) effect (Alcock & Paczynski 1979; Ballinger et al. 1996; Ruggeri et al. 2017) and redshift space distortions (Kaiser 1987; Hamilton 1998).

These probes are expected to observe the imprint of neutrinos through their suppression of matter clustering. Neutrinos are very light and so have large free-streaming lengths. For scales below these they do not cluster. For a fixed total matter density, this results in matter perturbations being slightly smoothed. However, the cumulative effect on the matter power spectrum, $P\left(k\right)$, by $z=0$ is a decrease on linear scales by a factor, $\sim 1-8f_{\nu}$, in linear theory, where $f_{\nu}\equiv\Omega_{\nu}/\Omega_{m}$ is the ratio of the neutrino and total matter density parameters (Hu et al. 1998). On non-linear scales, neutrinos are modelled through perturbative methods (e.g., Pietroni 2008; Lesgourgues et al. 2009; Levi & Vlah 2016) or accounted for in simulations. Generally, the most accurate approaches for the latter have implemented neutrinos as separate low-mass collisionless particle species alongside cold dark matter (CDM) to fully capture their non-linear evolution (e.g. Bird et al. 2012; Massara et al. 2014; Liu et al. 2018). These indicate that the effect of neutrino suppression on haloes depends on their mass, leading to a greater power reduction than linear theory on intermediate scales of $k\sim 1h\,\mathrm{Mpc}^{-1}$ where larger haloes are significant. On deeper non-linear scales where smaller haloes predominate, the power, while suppressed overall, is boosted relative to the linear prediction.

Depending on the choice of probes, methodology and underlying model, 1-$\sigma$ errors on the neutrino mass have been forecast broadly within the range $0.02-0.04\,\rm{eV}$ for Euclid-like surveys (e.g., Carbone et al. 2011; Audren et al. 2013; Boyle & Komatsu 2018). A major systematic neglected in these forecasts is the redistributing of matter on halo scales by baryonic astrophysics. The mechanisms governing adiabatic contraction via radiative cooling or localized high-energy releases from active galactic nuclei (AGN) and supernovae (see e.g., Rudd et al. 2008; van Daalen et al. 2014) are not well-understood. AGN feedback can decrease the fraction of baryons in inner halo regions by several factors and bloat haloes outwards as heated gas expands through the matter distribution (e.g., Duffy et al. 2010). This has been well-demonstrated by high-resolution hydrodynamic simulations such as the OverWhelmingly Large Simulations (OWLS; Schaye et al. 2010). A degree of counterbalance is provided on non-linear scales by the infall of dark matter into gravitational potential wells that have been deepened by the small-scale clustering of cooling gas (Gnedin et al. 2004; Jing et al. 2006; Rudd et al. 2008; Duffy et al. 2010). Overall, $P\left(k\right)$ experiences a $\sim$30% suppression in the mildly non-linear regime, which gives way to a net increase of power through adiabatic contraction on deeply non-linear scales of $k\sim 10h\,\mathrm{Mpc}^{-1}$ (e.g., Semboloni et al. 2011). As neutrinos and baryons both suppress matter power on weakly non-linear scales it is possible that a partial degeneracy exists between the two effects. It is therefore necessary to simultaneously constrain them both. Such degeneracies could potentially be broken by the scale-dependence of the baryonic effect and the redshift-dependence of neutrino suppression.

Marginalizing over the uncertainty in baryon processes can significantly impact parameter forecasts. For example, Copeland et al. (2018) predict 40% degradations to constraints of the dark energy parameter space and up to 80% for other cosmological parameters for a Euclid-like weak lensing survey (see also Semboloni et al. 2011; Zentner et al. 2013; Mohammed et al. 2014). Significant biases in neutrino mass measurements were found by Natarajan et al. (2014) when failing to account for baryons, using simple arguments to model both components in $P\left(k\right)$, although the impact of marginalizing over baryon parameters was not explored. More recently, Parimbelli et al. (2018) found that the neutrino mass is unlikely to be significantly biased by degeneracies with baryons, based on a Markov Chain Monte Carlo (MCMC) method using weak lensing and galaxy clustering that is limited by excluding freedom in the cosmological parameters. This paper explores similar territory but addresses different key questions. We are primarily interested in the degradations to error forecasts from marginalizing over baryonic effects, although the separate issue of bias in the neutrino mass itself is examined in the context of miscalibrations in baryon modelling. The specific focus of this paper is how degradation impacts the capacity to measure $\Sigma_{\rm{NH}}$ and $\Sigma_{\rm{IH}}$, and whether it will be feasible to distinguish between the mass hierarchies. To do this we perform a Fisher analysis of the full range of cosmological parameters probed by Euclid-like weak lensing and galaxy clustering surveys111It is important to emphasize that while this work uses Euclid survey parameters and specifications from e.g., Laureijs et al. (2011) and Amendola et al. (2018), this paper is not the product of a Euclid Consortium collaboration, and the results presented here should not be taken as official Euclid forecasts..

We employ analytic modifications of the halo model to capture the effects of neutrinos and baryons on $P\left(k\right)$ and the weak lensing convergence power spectrum, $C_{\ell}$. Mead et al. (2015) (hereafter M15) provide an empirically driven baryon prescription by calibrating internal halo structure relations to match several implementations of multiple sources of baryonic physics in OWLS. We use this halo model variant, HMCODE, in our work. Improving upon the accuracy of HALOFIT (Smith et al. 2003; Takahashi et al. 2012) by up to a factor of 2 in the non-linear regime, it achieves $\simeq 5$ percent accuracy in $P\left(k\right)$ for scales $k\leq 10\,h\,\mathrm{Mpc}^{-1}$ and redshifts $z\leq 2$. For a wider discussion of the value of this approach compared to treatments that fit directly the stellar and gas physics within the halo (e.g., Semboloni et al. 2011; Mohammed et al. 2014; Fedeli 2014), we refer the reader to Copeland et al. (2018). Mead et al. (2016) (hereafter M16) extend their model to account for non-linear neutrino influences by modifying the parameters governing spherical collapse, as this will be directly impacted by the reduction in matter clustering. Fitting to the simulations of Massara et al. (2014) reproduces $P\left(k\right)$ to within a few percent up to $k=10\,h\,\mathrm{Mpc}^{-1}$ for multiple redshifts. This is a slightly better performance than the Bird et al. (2012) fitting formula.

In our weak lensing analysis, we address one further source of uncertainty. On large scales, gravitational fields cause tidal distortions in galaxies, which contribute to correlations between the intrinsic ellipticities and gravitational shear of galaxies. Failing to account for these intrinsic alignments (IA) can significantly bias lensing parameter estimates (Joachimi et al. 2011; Troxel & Ishak 2015).

Our work is the first to combine baryons, IA and massive neutrinos in a self-consistent Fisher framework in the context of distinguishing the NH and IH in Stage IV surveys. By applying this framework to weak lensing, spectroscopic galaxy clustering and CMB probes, we present a uniquely comprehensive analysis of the current prospects for determining the hierarchy. This paper proceeds as follows. In § 2 we outline the neutrino and baryon modifications to the halo model. We present our results in § 3 for the mass constraints on the normal and inverted hierarchies, analyze the degradation due to baryons and intrinsic alignments, and explore the information gains available from combining weak lensing and galaxy clustering for a Euclid-like survey with priors from BOSS (Alam et al. 2017) and Planck. The significance of model bias is assessed in § 4, before we conclude in § 5. The LSS probes and Fisher formalism we use to make forecasts are explained in Appendix B.

2 Neutrinos and Baryons In The Halo Model

For a review of the halo model and how the matter power spectrum before including neutrinos and baryons is modelled, we refer the reader to Appendix A.

2.1 Modelling neutrino effects

We use the CAMB Boltzmann code (Lewis et al. 2000) to generate linear $P\left(k\right)$ for the NH and IH. To capture the effect of neutrinos on non-linear scales we use the modifications to the halo model prescribed by M16. It should be noted that our version of HMCODE has been slightly altered from that of M15 and M16, partly to correct for small instabilities that we originally encountered. Massive neutrinos impact the matter distribution primarily by damping matter perturbations as they free-stream through structure. For fixed $\Omega_{m}$222We define the total matter density parameter, $\Omega_{m}=\Omega_{c}+\Omega_{b}+\Omega_{\nu}$, as the sum of contributions from cold dark matter, baryons and neutrinos., the suppression of matter clustering is well-approximated as depending only on $f_{\nu}$. This affects the spherical collapse overdensity, $\delta_{c}$, and the virial density, $\Delta_{v}$, which are characteristic descriptors of halo structure. The former affects the shape of the mass function and determines the halo concentration by defining the formation redshift at which a certain fraction of mass has undergone spherical collapse. The latter defines limits on the density profile. By introducing a linear correction in terms of $f_{\nu}$, such that

[TABLE]

and

[TABLE]

M16 achieve a percent-level fit to N-body massive neutrino and CDM simulations by Massara et al. (2014) up to $k=10\,h\,\mathrm{Mpc}^{-1}$. These relations can be interpreted as the damping of matter fluctuations due to neutrino free-streaming resulting in a smaller fraction of perturbations collapsing by a given redshift, equivalent to raising the density threshold for collapse. Potential improvements could be attained by implementing more sophisticated methods by Ichiki & Takada (2012) or LoVerde (2014) in which, for example, non-linear neutrino clustering effects are considered. However, the M16 fit is sufficiently accurate for our purposes and permits us to use $f_{\nu}$ as a free parameter in fast power spectrum calculations for forecasts. In practice, we treat the mass sum, $\Sigma$, as the actual free parameter, via the relation,

[TABLE]

The Massara et al. (2014) simulations are designed to test the halo model with neutrinos, specifically by providing accurate halo power spectrum responses to different neutrino masses. This makes the M16 fits to these simulations particularly relevant for our purposes, in which an accurate power spectrum is less important than the sensitivity of the power spectrum to changes in cosmology for making forecasts. Other implementations of neutrinos in the halo model (e.g., Takahashi et al. 2012; Bird et al. 2012) generate slightly different non-linear responses in the power spectrum. It should be noted that none of these treatments is entirely robust, being subject to the limited accuracy of the underlying halo model fit on non-linear scales and the challenge of properly accounting for neutrino non-linear clustering effects. There may be future scope to improve on these approaches by extending the methodology of Mead (2017) and Cataneo et al. (2018), who demonstrate the increased accuracy on non-linear scales available to the halo model by incorporating spherical collapse and growth results into power spectrum responses.

In Figure 1 we show the impact of varying $\Sigma$ on the halo mass function with and without the modification to $\delta_{c}$ in equation (1). Neutrinos have a very small impact on the matter distribution on non-linear scales so we have illustrated it through the ratios of the mass function to the case with the minimal NH neutrino mass sum, $\Sigma_{\rm{NH,min}}=0.06\,\rm{eV}$.

Considering first the case with $\delta_{c}$ unmodified, we find that for fixed $\Omega_{m}$ and $\sigma_{8}$, increasing the neutrino mass increases the mass function for high halo masses in the top panel of Figure 1. This reflects the fact that the suppression of power on small scales from neutrino free streaming requires that linear power is boosted on large scales to keep $\sigma_{8}$ fixed. As a result, the variance of the matter distribution is increased for the corresponding large mass scales. This is the regime where $\sigma\left(M\right)$ is smallest, and hence where $\nu$ is largest. At these scales the mass function is dominated by an exponential tail, as shown in equation (12). This leads to small changes in $\nu$, from varying $\Sigma$, generating large changes in $n\left(M\right)$.

The non-linear regime begins on smaller mass scales where the clustering is damped by higher $\Sigma$, reducing the probability of matter fluctuations being dense enough to collapse and contribute to $n\left(M\right)$. Equation (12) is normalized to the mean density, so in order to satisfy the condition that integrating the mass function over all scales returns $\bar{\rho}$, massive neutrinos amplify $n\left(M\right)$ for very small $M$. However, these are highly non-linear scales that are inaccessible to Stage IV surveys and so of limited interest. The higher threshold for collapse in equation (1) suppresses the number density of high-mass haloes. This will propagate into a reduction of matter power on the weakly non-linear scales corresponding to such objects. The suppression of clustering due to $\Sigma$ increases by a factor of $\sim 5$ with this change to $\delta_{c}$, although it should be emphasized that the overall effect of massive neutrinos on the mass function remains very small.

The impact of massive neutrinos on the density profile has a less significant (but still non-trivial) contribution to the matter power spectrum compared to changes in the mass function. The collapse overdensity partially controls the formation redshift determining the concentration relation in equation (18). At late times, regions are more likely to collapse due to the amplitude of matter fluctuations being greater. Increasing $\delta_{c}$ makes halo formation rarer, favouring later times over earlier times, and therefore lowering $z_{f}$. The resulting reduction in halo concentration translates to an increased scale radius for fixed $r_{v}$. Consequently, for a fixed halo mass, the density is suppressed on small scales, with a milder boost on large scales. However, when the reduction of $r_{v}$ through increasing $\Delta_{v}$ via $\Sigma$ is also accounted for, this picture is complicated. If we consider this effect in isolation, for fixed mass and concentration, the density on small scales must increase to reflect the fact that haloes are smaller when meeting the increased density threshold for virialization. Interestingly, when the sensitivity to both $\delta_{c}$ and $\Delta_{v}$ is included, the two contributions mostly cancel, leaving the halo density profile with minimal sensitivity to the neutrino mass on small scales. The relative insensitivity of the density profile, and more importantly the mass function, highlights, at the outset, the challenges of using the non-linear regime to constrain $\Sigma$.

2.2 Baryonic contributions

We adopt the approach of M15 to capture the effects of baryonic astrophysics on haloes333This subsection is a compressed version of a more comprehensive discussion originally presented in Copeland et al. (2018), to which we refer the interested reader. An important addition is the real-space halo density profile with a baryonic feedback modification in equation (5). This was illustrated in Figure 2 of Copeland et al. (2018) but we present the expression for the first time here.. Internal halo structure relations are modified to include two parameters, $A_{B}$ (referred to as $A$ in M15) and $\eta_{0}$, that mimic adiabatic contraction via radiative cooling, and baryonic feedback (from e.g., AGN and supernova sources) respectively. These are fit to a range of simulations with different implementations of baryons (and a dark matter only case) provided by OWLS. Fiducial values, $A_{B,\rm{fid}}=3.13$ and $\eta_{0,\rm{fid}}=0.603$, are determined by fitting to power spectra produced by the emulator code, COSMIC EMU, designed for the high resolution N-body simulations from the Coyote Universe project (Heitmann et al. 2009, 2010; Lawrence et al. 2010; Heitmann et al. 2014). For this fit, HMCODE achieves $\simeq 5$ percent accuracy for the redshifts and scales of interest to Stage IV surveys.

Other methods for modelling baryons introduce additional density profiles describing specific physical contributions. An example of an alternative method for parameterizing baryons is that of Schneider & Teyssier (2015), whose ‘baryon correction model’ incorporates star and gas components to their density profiles to reproduce specific physical contributions. This model is used by Parimbelli et al. (2018) in their analysis of neutrinos and baryons. We advocate instead the M15 philosophy of making empirical corrections to halo model relations because these are primarily motivated by generating the most accurate power spectrum, which is the underlying probe used for our forecasts. For our purposes we require an accurate understanding of the response of halo structures to baryons rather than a detailed framework for the underlying physics. Similar reasoning has previously been applied by surveys such as KiDS (Hildebrandt et al. 2017) for using the M15 model to account for baryons in their data analyses.

In M15, adiabatic contraction is captured by varying the amplitude, $A_{B}$, of the concentration factor in equation (18). Radiative cooling of baryonic matter causes it to cluster, thereby deepening gravitational wells into which dark matter falls, and so leading to the adiabatic contraction of the total matter distribution. On non-linear scales, simulations show that the effect has an impact of several percent (Duffy et al. 2010; Gnedin et al. 2011). This manifests as a reduced scale radius in the halo profile, suppressing the density on large scales and increasing it at smaller scales for a fixed virial radius and halo mass. This is demonstrated in Figure 1 of Copeland et al. (2018), which shows density profiles for a range of $A_{B}$ values.

Baryonic feedback from sources like AGN and supernovae is far less straightforward to model. Significant energy exchanges with their environments involve isotropic radiative transfer or highly directional jets in the case of the former (Schaye et al. 2010; van Daalen et al. 2011; Martizzi et al. 2014) and multiple sources heating surrounding gas in the latter (Pontzen & Governato 2012; Lagos et al. 2013). In both cases envelopes of gas can expand on scales comparable to the virial radius, so a significant range of scales from subparsecs to megaparsecs are affected. This makes modelling feedback analytically difficult, a challenge that is exacerbated by the effect on haloes depending on mass as well as scale. As simulations have demonstrated (e.g., Pontzen & Governato 2012; Teyssier et al. 2013; Martizzi et al. 2013), similar AGN mechanisms for the expulsion of gas from central regions for both lower and higher mass haloes can remove substantial baryonic matter from the former while merely bloating the latter outwards in a less catastrophic process.

This range of effects is captured in M15 by transforming the scale of the window function,

[TABLE]

by the mass-dependent factor, $\nu=\frac{\delta_{c}}{\sigma\left(M\right)}$. We find that in real space the modified density profile takes the form

[TABLE]

We refer the reader to Figure 2 in Copeland et al. (2018), which shows that for more positive values of $\eta$ higher mass haloes (characterized by $\nu>1$) are increasingly bloated while lower mass haloes ($\nu<1$) are left relatively diminished by the outright removal of gas. M15 achieve their best power spectra fits, including for dark-matter-only simulations, when using a non-zero value of $\eta$ with a redshift dependence. The parameter was decomposed such that

[TABLE]

with the constant, $\eta_{0}$, controlling the strength of the feedback impact on halo structure.

Copeland et al. (2018) extend HMCODE by introducing an inner halo core that could parameterize the effects of baryons or, potentially, other physics on the smallest scales of structure. However, we do not include this degree of freedom in this work, and simply use the NFW profile, with $A_{B}$ and $\eta_{0}$ as our baryon parameters.

3 Neutrino Mass Forecasts

We construct $\left(A_{B},\eta_{0},\Omega_{m},\Omega_{b},\Sigma,h,n_{s},\sigma_{8},w,a_{\rm{IA}}\right)$ Fisher matrices for Euclid-like weak lensing and spectroscopic galaxy surveys. We effectively apply a uniform prior on $\Sigma$ by using the Fisher formalism. This is justified by the fact that different $\Sigma$ values should not significantly change the shape of the multivariate likelihood function. Such a prior can be considered a desirable objective prior that allows the posterior to be determined instead by the likelihood (Heavens & Sellentin 2018). The primary goals of this work are to assess the impacts of IA (a phenomenon exclusive to weak lensing) and baryonic effects on neutrino mass forecasts. As the influence of baryons degrades information significantly on scales $k>1h\,\mathrm{Mpc}^{-1}$ this effect is of far greater concern to weak lensing forecasts than those from galaxy clustering. Hence, our analysis takes the approach of first discussing weak lensing results and the extent of their degradation due to baryons and IA. Information from Planck and galaxy clustering is then added to form a more complete picture of the prospects of a Euclid-like survey to measure the neutrino mass with sufficient accuracy to distinguish the hierarchies.

3.1 Normal and inverted hierarchy results from weak lensing

In Table 1, we state the weak lensing survey parameters specified by the Euclid survey report (Laureijs et al. 2011). $N_{z,\rm{WL}}=10$ redshift bins are chosen in the range $0\leq z\leq 2$ such that each bin contains an equal number density of galaxies. A large range of scales from $\ell_{\mathrm{min}}=10$ to $\ell_{\mathrm{max}}=5000$ are covered so in practice we compute the summation in equation (41) at logarithmic intervals. We include full confidence ellipse plots for the NH in Appendix C. The fiducial values of the cosmological parameters in this work are given by the base $\Lambda$CDM Planck TT,TE,EE+lowP likelihood (see Table 4 in Planck Collaboration et al. 2016). For the neutrino masses we use the minimal bounds as fiducial values i.e., $\Sigma_{\rm{NH,fid}}=0.06\,\rm{eV}$ and $\Sigma_{\rm{IH,fid}}=0.1\,\rm{eV}$.

One can derive significant insight into the results of a Fisher forecast by examining the response of the power spectra to varying parameters across relevant scales. This diagnostic can quickly indicate which scales most Fisher information comes from or is potentially lost due to degenerate power spectrum responses between parameters. In Figures 3 and 4 we show, for the normal hierarchy, responses of $\Delta^{2}\left(k\right)$ at $z=0$ and $C_{\ell}$ respectively to varying each parameter with respect to its fiducial value while fixing the other parameters. In § 3.2 we discuss the responses to the baryon parameters. $\Delta^{2}\left(k\right)$ responses to most parameters exhibit nodes at the scale corresponding to $\sigma_{8}$, at which we normalize power. For example, the response to $\Sigma$ is an enhancement of power on linear scales in order to satisfy power being damped in the non-linear regime due to free streaming. The impact of neutrinos on these scales is small, determined by the limited sensitivity of the mass function and the halo structure parameters, $\delta_{c}$ and $\Delta_{v}$, to $\Sigma$ (discussed in detail in § 2.1). Most of the $C_{\ell}$ responses reflect the corresponding $\Delta^{2}\left(k\right)$ responses, when accounting for redshift-dependent effects along the line of sight. However, parameters that affect the lensing weight function through their impact on cosmological distances can exhibit significantly different $\Delta^{2}\left(k\right)$ and $C_{\ell}$ responses. For example, the competing influences of geometry and growth are apparent in the $w$ responses, which Copeland et al. (2018) examine in depth.

We compare the $\Sigma$ responses for $\Delta^{2}\left(k\right)$ and $C_{\ell}$ for the NH and IH directly in Figure 5. $\Sigma_{\rm{IH,min}}=0.1\,\rm{eV}$ is less than double $\Sigma_{\rm{NH,min}}=0.06\,\rm{eV}$ but the IH exhibits power responses that are more than twice as large as the NH responses. The leading order cause of this is likely a dependence of power on $\Sigma$ that grows with the mass sum from a quadratic minimum which is reached when $\Sigma$ approaches zero.

An additional effect to account for is the influence of the neutrino free-streaming scale, $k_{\rm{fs}}$, which increases with mass. There are different $k_{\rm{fs}}$ for each mass eigenstate in a hierarchy but, within the approximation of the NH as $\left\{m_{1}\approx m_{2}\approx 0\,\rm{eV},m_{3}\approx 0.06\,\rm{eV}\right\}$ and the IH as $\left\{m_{1}\approx m_{2}\approx 0.05\,\rm{eV},m_{3}\approx 0\,\rm{eV}\right\}$, the most relevant free-streaming scales satisfy $k_{\rm{fs,IH}}<k_{\rm{fs,NH}}$. Fixing $\sigma_{8}$ means that, for increases in $\Sigma$, large-scale power is enhanced relative to the fiducial case to satisfy the suppression of power by neutrinos on small scales. There is a similar contribution from the peak of $P\left(k\right)$ shifting with a decrease in the redshift of matter-radiation equality, which depends on $\Omega_{c}$ (which decreases to keep $\Omega_{m}$ constant for increasing $\Sigma$) due to the relativistic nature of neutrinos during this epoch. As the effect of free-streaming neutrinos arises at smaller $k$ for the IH, the relative enhancement of the large-scale IH power is greater than for the NH. Similarly, the responses for the IH are more sensitive to the same fractional change in $\Sigma$ compared to the NH, which we see reflected in Figure 5.

When propagating $\Delta^{2}\left(k\right)$ into $C_{\ell}$, the lowest $\ell$ for a Euclid-like survey correspond to $k$ where the matter power response is weak. As a result, the strongest linear responses for $C_{\ell}$ in Figure 5 are several factors smaller than for $\Delta^{2}\left(k\right)$. Other parameter responses generally experience significantly more modest changes from one probe to another. It should be noted that the background geometry is not sensitive to increasing $\Sigma$ when neutrinos are non-relativistic; the angular distances entering equation (20) depend on $\Omega_{m}$, which is kept constant.

A consequence of the weak lensing probe only being sensitive in the linear regime to scales where the neutrino response is most limited is that the magnitude of the $C_{\ell}$ response is similar between the linear and non-linear regimes. The response for the latter is due mainly to the modest increases in $\delta_{c}$ and $\Delta_{v}$. This highlights the importance of pursuing robust physical modelling for neutrino effects on structure at these scales because they are as sensitive to $\Sigma$ as linear scales.

Weak lensing forecasts of $\Sigma$ are particularly challenging due to this limited sensitivity and the degeneracies apparent between $\Sigma$ and a range of other parameters on non-linear scales. A well-known example is the degeneracy with $\sigma_{8}$, in which both exhibit a ‘spoon’-shaped response (see e.g., Massara et al. 2014). An increase in $\Sigma$ most strongly affects high mass haloes, which dominate over small mass haloes on intermediate scales. Figure 1 shows that in this region higher $\Sigma$ reduces the linear clustering of matter, and therefore the number of high mass haloes that can form, lowering the halo mass function and hence the non-linear matter power. Raising the collapse and virialization density thresholds with $\Sigma$ deepens the resulting dip in the power response and shifts it to larger scales associated with higher masses. At higher $k$ the impact of $\Sigma$ is somewhat decreased and becomes less scale-dependent.

The bump in the $\sigma_{8}$ response is qualitatively similar. By altering the clustering amplitude, a halo mass-dependent and hence scale-dependent power response is induced through the change in the peak height, $\nu=\delta_{c}/\sigma\left(M,z\right)$, a quantity that is also sensitive to $\Sigma$. The response becomes roughly scale-independent on smaller scales. This is similar to the $\Sigma$ response but with the opposite sign, enhancing rather than suppressing power. As most information comes from small scales this degeneracy is significant.

In Figure 6 the contribution per $\ln\ell$ to the Fisher matrix for the mass sum is shown for the NH and IH. This effectively weights the Fisher information by the signal-to-noise in each $\ell$ mode, which increases with $\ell$. The greater number of independent modes in the non-linear regime, and consequently lower noise, results in most of the available information being drawn from these scales. The greater sensitivity of information at high $\ell$ highlights the limitations on achieving strong constraints, given the limited neutrino influence on non-linear power. The IH contributes almost twice as much information as the NH at the peak sensitivity.

The main difference between the hierarchies in their information contribution is from different logarithmic derivatives of $C_{\ell}$. These are partially predicted by the power responses. On non-linear scales the $C_{\ell}$ response for the IH is larger than double the NH for the same fractional change in $\Sigma$. As $\Sigma_{\rm{IH,min}}=0.1\,\rm{eV}$ is less than double $\Sigma_{\rm{NH,min}}=0.06\,\rm{eV}$, this leads to $\partial\ln C_{\ell}/\partial\Sigma_{\rm{IH}}$ being sufficiently larger than $\partial\ln C_{\ell}/\partial\Sigma_{\rm{NH}}$ that the Fisher sensitivity, which depends quadratically on the logarithmic derivative, is up to twice as large.

In Figure 7 we show NH and IH confidence ellipses for $w$-$\Sigma$ for a Euclid-like weak lensing (WL) survey. We do not expect to constrain $\Sigma$ from WL alone, but we are interested in using it as a baseline probe, through which we can examine both the impact of systematics like baryons and IA and the improvements available from adding CMB and galaxy clustering information. Before marginalizing over baryons, we find WL errors of $\sigma_{\Sigma_{\rm{NH}}}=0.079\,\mathrm{eV}$ and $\sigma_{\Sigma_{\rm{IH}}}=0.120\,\mathrm{eV}$ (see by Table 2). We treat the minimal mass for each hierarchy ($\Sigma_{\rm{NH}}=\Sigma_{\rm{NH,min}}=0.06\,\mathrm{eV}$ and $\Sigma_{\rm{IH}}=\Sigma_{\rm{IH,min}}=0.1\,\mathrm{eV}$), as the fiducial cases in this work, so distinguishing the hierarchies requires that error forecasts fall within the threshold, $|\Delta\Sigma_{\rm{min}}|=|\Sigma_{\rm{NH,min}}-\Sigma_{\rm{IH,min}}|=0.04\,\rm{eV}$. It should be noted that any measurement of $\Sigma\geq 0.1\,\rm{eV}$ permits either hierarchy regardless of the strength of the error constraint. The errors we find are therefore too large for a distinction to be made. They would also be too large to permit any positive detection of neutrinos in the scenarios where the mass sum is given by the minimal mass for each hierarchy.

The latter issue can be alleviated by adding priors on cosmology from sources like galaxy clustering or the early Universe. We focus here on the latter through Planck CMB anisotropy measurements, while impact of galaxy clustering is discussed in detail in S 3.4. The CMB provides a wealth of information on the matter energy-density and the geometry of the Universe, which is one of the main contributions to LSS signals. We use the publicly available MCMC chains for the base $\nu\Lambda$CDM combined TT, TE and EE power spectra (see Planck Collaboration et al. 2016) to construct a covariance matrix which is inverted into a Fisher matrix, $F_{\rm{CMB}}$, that encapsulates the information available from the CMB. This is added to the Euclid-like WL information such that the total Fisher information is

[TABLE]

As can be seen in Figure 7 there is a substantial improvement to $\Sigma$ forecasts for the NH and IH, with errors reduced below the distinction threshold to $\sigma_{\Sigma_{\rm{NH}}}=0.036\,\mathrm{eV}$ and $\sigma_{\Sigma_{\rm{IH}}}=0.033\,\mathrm{eV}$. This is partly a result of parameter degeneracies being broken, for example the CMB provides strong constraints on $\sigma_{8}$ which propagate through to our final $\Sigma$ forecasts. It should be noted that this requires CMB experiments to accurately measure the reionization optical depth, $\tau$, as it is degenerate with the primordial power amplitude. We marginalize over $\tau$, along with a range of other parameters, when constructing $F_{\rm{CMB}}$. Cosmological neutrino measurements stand to benefit significantly from future 21 cm experiments that will probe the reionization epoch with high precision (Allison et al. 2015; Liu et al. 2016). However, this work does not extend to examine this in detail.

The CMB prior provides a significant contribution to our forecasts by breaking parameter degeneracies and, as it is the same for both hierarchies, it leads to less difference between the NH and IH for the WL+CMB forecasts. This is reflected in the weaker WL constraint, $\sigma_{\Sigma_{\rm{IH}}}$, which improves by 72%, benefiting significantly more than $\sigma_{\Sigma_{\rm{NH}}}$, which improves by 54% (see Table 2). Our results are comparable to findings from Audren et al. (2013) who use an MCMC approach to explore the degenerate hierarchy case. However, even these improvements would only achieve a close to 2$\sigma$ detection of $\Sigma$ in the case of the NH (although the IH case fares better with a 3$\sigma$ detection) and a roughly 1$\sigma$ determination of the hierarchy.

3.2 Impact of baryons

When the baryon systematic is accounted for, the prospects of using weak lensing to accurately distinguish hierarchies recede. Figures 3 and 4 include the $\Delta^{2}\left(k\right)$ and $C_{\ell}$ responses to the adiabatic contraction parameter, $A_{B}$, and the baryonic feedback parameter, $\eta_{0}$. Copeland et al. (2018) discuss these in detail while here we highlight key features. The boost to non-linear power from increasing $A_{B}$ is a reflection of the enhanced halo density profiles in this regime (see Figure 1 of Copeland et al. 2018). The response to $\eta_{0}$ exhibits a peak, which is a more subtle consequence of capturing the effects of baryonic feedback over a range of mass and spatial scales. The bloating of higher-mass haloes is the dominant influence compared to the reduction effect experienced by lower-mass haloes, so decreasing $\eta_{0}$ generates a net increase of $\Delta^{2}\left(k\right)$. The peak occurs at deeper non-linear scales for higher redshifts, reflecting the evolution of halo populations.

Baryons and neutrinos both impact the amplitude of the matter power spectrum on small scales, from which most information is drawn, leading to significant degeneracies between both baryon parameters and $\Sigma$. For example, the mass-dependent response for $\eta_{0}$ emerges from its effect on the peak height which, as discussed above, determines non-linear halo properties in a similar manner to the neutrino mass sum. In the M16 model, the peak height depends on $\Sigma$ directly. The total baryon impact on the $\Sigma$ forecast is shown in Figure 7 to be severe. The errors for NH and IH double for the former and increase by 70% for the latter. The greater sensitivity of $\Sigma\Sigma$ Fisher information for the IH relative to the NH on non-linear scales where the baryons limit information results in this comparatively lesser degradation.

In Copeland et al. (2018) we explored several strategies for mitigating baryon impacts on dark energy forecasts. These included increasing the Euclid-like survey limit from $\ell_{\rm{max}}=5000$ to $\ell_{\rm{max}}=10000$ to access information from modes deeper in the non-linear regime. The relative degradation of constraints on the dark energy equation of state increases, although the absolute constraints do improve. For $\Sigma$ we find that baryon degradation is reduced, to 84% and 61% for the NH and IH respectively. Figure 6 shows that there is a repository of potential Fisher information on $\Sigma$ at higher $\ell$. Our response plots suggest that the degeneracies between $A_{B}$, $\eta_{0}$ and $\Sigma$ may be less severe for very high $\ell$ where the neutrino impact on very small haloes approaches scale-independence (Massara et al. 2014) in contrast to the baryons. This may explain why the degradation is reduced, despite the increasing influence of baryons on non-linear scales overall, but ultimately this improvement is limited.

We also explored the gains from adding Fisher information from an external baryon source (e.g., new simulations or improved observations) which in general does mitigate the baryon impact on $\Sigma$, but to remove it as a statistically significant systematic requires more baryon information than could realistically be accessed at present. Therefore, we focus on the most promising source of baryon mitigation from Copeland et al. (2018), the constraints on cosmology from the CMB (Planck Collaboration et al. 2016). These propagate through the Fisher analysis, contributing to the breaking of parameter degeneracies. This limits the baryon degradation, and in the case of neutrinos we find a significant improvement with the impacts of marginalizing over baryons reduced to 21% and 28% for the NH and IH. Although the IH still has a smaller $\Sigma$ error it now experiences greater relative baryon degradation than the NH when CMB priors are included. Despite the improvements, the baryon impact is detrimental for the weak lensing probe, raising the WL+CMB errors beyond the 1$\sigma$ hierarchy distinction threshold, $|\Delta\Sigma_{\rm{min}}|$, to $\sigma_{\Sigma,\rm{NH}}=0.044$ and $\sigma_{\Sigma,\rm{IH}}=0.042$. These would constitute confidence in a positive detection of massive neutrinos at the 1$\sigma$ and 2$\sigma$ level for $\Sigma_{\rm{NH}}$ and $\Sigma_{\rm{IH}}$ respectively. It will therefore be important to include additional information from other sources, such as galaxy clustering.

3.3 Including intrinsic alignments

Biased parameter estimations from failing to account for intrinsic alignments have been extensively studied but it is also important to consider the impact on forecasts of including IA in the $C_{\ell}$ signal and then marginalizing over $a_{\rm{IA}}$. Our results are broadly consistent with the predictions of Krause et al. (2016) for a Euclid-like survey, with the caveat that we do not extend our analysis to marginalizing over nuisance parameters governing a potential luminosity scaling of the IA amplitude, and so we report slightly less significant impacts. We show, in Figure 8, a direct comparison between hierarchies for the $w$-$\Sigma$ forecasts with and without IA.

Our results indicate that the impact is small for most parameters, including $\Sigma$. This can be understood by considering the $C_{\ell}$ responses with and without IA for different redshift bin combinations. The $C_{\ell}^{\rm{tot}}$ responses for the $ij=44$ $z$-bin ($0.79\leq z\leq 0.89$) autocorrelation in Figure 4 are sufficiently similar to the corresponding $C_{\ell}^{\gamma\gamma}$ responses that it is not useful to distinguish between them by explicitly including the latter. This domination of the shear-shear signal is most extreme for auto-correlations in the highest redshift bins, as Figure 2 shows. However, for cross-bin correlations the IA-shear signal becomes important as it arises from the alignment of foreground galaxy ellipticities with background lens potentials. This term is negative and large enough to provide significant cancellation to the shear-shear signal in these cases. In Appendix D, we include Figure 13 to illustrate the impact on the $C_{\ell}^{\rm{tot}}$ and $C_{\ell}^{\gamma\gamma}$ responses for the $ij=09$ cross-correlation. As the IA only significantly affect the cross-power between widely separated bins or the auto-power for low redshift bins, and most of the lensing signal comes from the minimally impacted high redshift bin auto-spectra, the $\Sigma$ Fisher information is similar between cases with and without IA. Any degradation of forecasts would therefore come from degeneracies between $\Sigma$ and $a_{\rm{IA}}$. Figures 4 and 13 show significantly different power responses to these parameters. When combined with the fact that only a small subset of auto- and cross-power spectra are significantly sensitive to the IA, we would therefore expect the impact on constraints to be small.

For both the NH and IH cases, the IA signal adds limited information on all parameters other than $w$ and $\Omega_{m}$. The determinant of the Fisher matrix, which effectively quantifies the total available information, is reduced by the inclusion of IA. However, when the inversion to a parameter covariance matrix is performed this manifests as some individual errors (and parameter correlations) decreasing while others increase. Neutrinos have a stronger effect on the power spectrum for the IH compared to the NH, and we find an even smaller, almost negligible, impact from IA on the $\Sigma_{\rm{IH}}$ forecast than for $\Sigma_{\rm{NH}}$. This suggests there is some sensitivity to the fiducial values of the parameters chosen. In turn this implies non-Gaussianity in the posterior, which would limit the validity of the Fisher approximation. However, the impact of IA for neutrino mass forecasts is overall of little concern, with only small improvements and degradations. The main systematic that must be mitigated remains the baryonic phenomena.

3.4 Galaxy clustering

We combine our WL+CMB forecasts with information from galaxy clustering (GC), which is derived from BOSS data at low redshifts ($0.2<z\leq 0.75$) and from a Euclid-like spectroscopic galaxy survey at higher redshifts ($0.75<z\leq 2.05$). We show, in Figure 9, confidence ellipses comparing the NH and IH with and without baryon marginalization for three cases: weak lensing only; and weak lensing and galaxy clustering; weak lensing, galaxy clustering and CMB priors combined. Table 3 contains the corresponding 1-$\sigma$ errors and the baryon degradation factors in each case. Neutrinos mainly impact the galaxy power spectrum through the matter power spectrum. The power sensitivity to $\Sigma$ is greater for the IH than the NH, leading to more substantial improvements to forecasts for the former. However, both hierarchies experience a significant reduction in errors, with combined WL+GC+CMB probes achieving $\sigma_{\Sigma_{\rm{NH}}}=0.034\,\rm{eV}$ and $\sigma_{\Sigma_{\rm{IH}}}=0.034\,\rm{eV}$, when accounting for baryons. The degradation due to baryons is also approximately halved in the combined case. Again, even though the IH experiences a far greater improvement overall, the relative degradation is almost double that of the NH at 16% and 9% respectively. Interestingly, when examining the impact of adding GC to WL without the inclusion of CMB priors, we see that the NH degradation is reduced by a greater fraction than that of the IH.

The inclusion of BOSS data on low redshift clustering was found to have a significant impact on galaxy clustering constraints compared to the Euclid-like only case. However, the combined WL+GC errors are only mildly improved by $\sim{5}$-10% while there is no appreciable change for the WL+GC+CMB errors. Therefore, in our final analysis the dominant information contribution from galaxy clustering is provided by the Euclid-like forecasts rather than BOSS data. Ultimately, our combined forecast constraints represent a nearly 2$\sigma$ (3$\sigma$) detection of $\Sigma_{\rm{NH}}$ ($\Sigma_{\rm{IH}}$) but are not significantly lower than the minimum threshold required to achieve a distinction between hierarchies.

Our results for a GC+CMB combination are broadly comparable with recent Fisher forecasts by Boyle & Komatsu (2018), although, as Boyle (2018) notes, making comparisons between results across the literature is difficult. Differences in the implementation of e.g., CMB priors, the use of Fisher versus MCMC methodology, the choice of parameters to be varied, and their fiducial values leads to a range of results. For example, Audren et al. (2013) and Sprenger et al. (2018) find stronger constraints in studies restricted to the degenerate mass hierarchy. For galaxy clustering, where neutrino information is coming primarily from the shape of the matter power spectrum on linear and mildly non-linear scales, this choice can have a significant impact. We find that the different free-streaming lengths characterizing the NH and IH (discussed in detail in § 3.1) can have a non-negligible effect on the power sensitivity to $\Sigma$. It follows that addressing whether Stage IV surveys will be able to measure $\Sigma$ or distinguish between hierarchies should not simply rely on predicting $\sigma_{\Sigma}$ for an arbitrary model, but should robustly account for the effects of different hierarchy mass splittings. This represents a sensitivity to fiducial parameters that will change the multivariate Gaussian likelihood function. When computing the Fisher information, one takes expectation values over the sampling distribution of the data given the true parameter values. Therefore, even for a likelihood that is Gaussian in the parameters, the Fisher matrix can still depend on the fiducial parameter values if the posterior covariance on the parameters entering the likelihood depends on the data. However, even if there is no dependence on the data (which is the case in this work), in practice there should still be a dependence on the true parameter values because we assume that the covariance is computed at these values.

3.5 Required improvements

To quantify the information required from an arbitrary external source to reduce the errors on $\Sigma_{\rm{NH}}$ and $\Sigma_{\rm{IH}}$ sufficiently for a Euclid-like survey to distinguish the hierarchies beyond the 1$\sigma$ level, we explore the effect of increasing diagonal elements of the (combined WL+GC+CMB) Fisher matrix by a factor, $\left(1+\alpha\right)$,

[TABLE]

In Figure 10 we plot the change in NH and IH errors for increasing $\alpha$ for the parameters $\Sigma$ and $n_{s}$. These provide the strongest improvements, as the impact of marginalizing over other parameters on $\sigma_{\Sigma}$ is already close to fully mitigated by the WL+GC+CMB combination of probes.

For $\Sigma$, while the WL+GC+CMB constraints alone constitute close to a 2$\sigma$ detection for the NH (with the IH being constrained to close to 3$\sigma$), this can be improved to between 3$\sigma$ and 4$\sigma$ by adding priors in the range $0.7\,\sigma_{\Sigma,\rm{con}}\lesssim\sigma_{\Sigma,\rm{prior}}\lesssim 2.0\,\sigma_{\Sigma,\rm{con}}$, where $\sigma_{\Sigma,\rm{con}}$ is the conditional error on the mass sum. It may also be possible to achieve a distinction between hierarchies (assuming that $\Sigma<0.1\,\rm{eV}$) at 2$\sigma$ or even 3$\sigma$ with priors in the range $2.0\,\sigma_{\Sigma,\rm{con}}\lesssim\sigma_{\Sigma,\rm{prior}}\lesssim 2.3\,\sigma_{\Sigma,\rm{con}}$. Oyama et al. (2016) find that combining a precise CMB polarization observation such as Simons Array with a 21 cm line observation from phase 2 of the Square kilometer Array and BAO priors could distinguish between hierarchies at the 95% confidence level provided that $\Sigma$ does not exceed $0.1\,\rm{eV}$. This may indicate scope for similarly significant improvements by combining future CMB polarization data with results from Euclid-like weak lensing surveys. The required information could also be contributed to by e.g., improved data from the Lyman-$\alpha$ forest or possibly further constraints from particle physics. An example of such a particle experiment is measuring the decay rate of neutrinoless double beta decay, as this depends on the sum of mass eigenstates, $m_{\beta\beta}$, when weighted by Majorana phases444Neutrinoless double beta decay is subject to the assumption that at least one neutrino is a Majorana fermion, i.e. that $\nu\leftrightarrow\bar{\nu}$. (see e.g., Capozzi et al. 2017). Although neutrinoless double beta decay has not yet been observed, the KamLAND-Zen experiment derives a lower bound of $m_{\beta\beta}\gtrsim 0.061\rm{-}0.165\,\rm{eV}$ (KamLAND-Zen Collaboration 2016).

There is also considerable scope for improving the $\Sigma$ constraints with additional priors on the spectral index which, if known with absolute certainty, would potentially allow a 2$\sigma$ (5$\sigma$) detection of $\Sigma$ for the NH (IH) and a 2$\sigma$ distinction between the hierarchies. The next generation CMB surveys CMB-S4 (Abazajian et al. 2016) predicts constraints on $n_{s}$ tightening by a factor of 2-3 over Planck, which could contribute to the improvements we see required for $\sigma_{\Sigma}$ in Figure 10. Note that the greater rate of decrease of $\sigma_{\Sigma_{\rm{IH}}}$ with $\alpha$ compared to $\sigma_{\Sigma_{\rm{IH}}}$ is a consequence of stronger degeneracies, e.g. between $n_{s}$-$\Sigma$, being broken for the IH compared to the NH.

4 Model Bias

The parameters controlling the baryon and IA systematics explored in this work are assigned their fiducial values from fits to simulations (Schaye et al. 2010; Mead et al. 2015) and shear observations (see e.g., Hirata & Seljak 2004; Bridle & King 2007) respectively. If the underlying physical models used in simulations is incorrect or if there are inherent limitations in the construction of the simulations, then these fiducial values are biased from the ‘true’ values. It is useful to consider the results of this paper in this context by calculating the degree of calibration bias required to severely bias $\Sigma$ forecasts, for example beyond their 1-$\sigma$ errors.

To first order, Taylor et al. (2007) showed that the bias in a forecast parameter, $\theta$, depends on the bias in a nuisance parameter, $\psi$, such as for baryons or IA, through sub-volumes of the Fisher information according to

[TABLE]

for which there is an implicit sum over $j$ and $k$. In Figure 11 we show the bias in the weak lensing forecasts for $w$-$\Sigma$ when the fiducial values of the baryon and IA parameters have been miscalibrated. We display results for model bias corresponding to the extreme values of the regions $2<A_{B}<4$ and $0.4<\eta_{0}<0.8$ that broadly encompass the range of fits to different OWLS simulations made by M15, and for bias in $a_{\rm{IA}}$ of 10%.

Figure 11 shows that bias in $\Sigma$ from miscalibrating the strength of the IA effect is significant, easily invalidating a measurement for either hierarchy with a 10% $a_{\rm{IA}}$ bias. This is in contrast to the very small impact on $\Sigma$ forecasts from marginalizing over $a_{\rm{IA}}$. However, as we find percent-level constraints on $a_{\rm{IA}}$ we might expect that a 10% miscalibration would indeed produce a very large bias. The problem of IA bias for other parameters such as $w$ is well-documented; we do not address it further here beyond reporting the impact on $\Sigma$. The bias due to baryons is less severe, although an extreme miscalibration of the adiabatic contraction parameter would bias $\Sigma$ beyond its 1-$\sigma$ confidence region. However, biases of $\eta_{0}$, even lying at the extremes of the region of values identified by M15, do not displace the $\Sigma$ estimate beyond the 1-$\sigma$ limit. This difference between $A_{B}$ and $\eta_{0}$ may also be reflected in the fact that marginalizing over $A_{B}$ alone has a significantly greater impact on $\Sigma$ forecasts than marginalizing over $\eta_{0}$ alone, with almost all of the baryon degradation being due to $A_{B}$. This highlights the importance of achieving tight calibrations of the adiabatic contraction parameter and the intrinsic alignment amplitude, while a similarly robust understanding of the strength of the baryonic feedback parameter is less urgent.

A subtlety worth noting is that a bias of $\delta\,\eta_{0}=-0.2$ produces biases of opposite sign in $\Sigma$ of different hierarchies. Because the sensitivity to $\eta_{0}$ is low, it may be that transitioning from one model to the other induces a small change to the bias in $\Sigma$ that is still sufficiently large to change its overall sign. This could complicate attempts to distinguish between hierarchies, but overall the $\eta_{0}$ bias remains subdominant to other effects.

As model bias is a product of subgrid limitations, resolving the issue is challenging without resorting to incorporating external data on baryon phenomenology into future simulations. In the case of IA, there has been substantial exploration into self-calibration methods of cleaning IA signals with independent information on IA correlations from deriving scaling relations that predict (and then subtract out) IA-shear signals from shear observables (Zhang 2010; Troxel & Ishak 2012; Yao et al. 2017, 2019). The IA-IA signal is strongest between close galaxies and so is generally assumed to be more straightforward to treat by purposely using sufficiently high $z$ redshift bins to limit its contribution (Schneider & Bridle 2010; Zhang 2010). A future avenue of research into applying these methods to mitigate the bias impact on parameters like the neutrino mass sum is worth pursuing.

5 Conclusions

In this work we have performed detailed Fisher forecast analyses to assess the potential of a Euclid-like Stage IV survey to measure the neutrino mass sum within sufficient accuracy to distinguish between the normal and inverted mass hierarchies. We applied particular focus to the risk of forecast degradation presented by the need to marginalize over baryonic astrophysics, and further explored the impact of intrinsic alignments. To evaluate the widest possible scope of constraints available at present, we combined forecasts for weak lensing with those of spectroscopic galaxy clustering from Euclid and added prior information from Planck CMB constraints and BOSS measurements of low redshift clustering.

The baryon aspect of this paper was an expansion of the baryon degradation analysis of $w_{0}$-$w_{a}$ constraints in Copeland et al. (2018) to the case of neutrinos. To explore the behaviour of baryons and neutrinos on non-linear scales, we used the halo model presented in Mead et al. (2015) and Mead et al. (2016). This captures baryons by incorporating their effect on adiabatic contraction and baryonic feedback into halo structure relations. M15 achieve accurate fits for the power spectrum to within a few percent by calibrating the associated parameters, $A_{B}$ and $\eta_{0}$, to the OWLS (Schaye et al. 2010) simulations. The impact of neutrinos is modelled by relating the mass sum to the spherical collapse overdensity and the virial density, which exert a significant influence on halo structure by changing the concentration function and the limits of the density profile respectively. Fits to simulations by Massara et al. (2014) reproduce $P\left(k\right)$ to within a few percent up to $k=10\,h\,\mathrm{Mpc}^{-1}$ over different redshifts, improving on the Bird et al. (2012) fitting formula.

We examined both the impact of the $\Sigma$ dependence of the power sensitivity to the mass sum and the effect of different free-streaming scales between hierarchies on the relative sensitivity of the IH compared to the NH. We illustrated that for both hierarchies, there is low sensitivity on scales relevant to our cosmological probes, $\Delta^{2}\left(k\right)$ and $C_{\ell}$. The resulting limited Fisher information available for a Euclid-like weak lensing survey leads to forecasts of $\sigma_{\Sigma_{\rm{NH}}}=0.079\,\mathrm{eV}$ and $\sigma_{\Sigma_{\rm{IH}}}=0.120\,\mathrm{eV}$. These are large enough to rule out making any definitive measurement for either hierarchy.

By incorporating prior CMB information from Planck Collaboration et al. (2016), these constraints were tightened significantly to $\sigma_{\Sigma_{\rm{NH}}}=0.036\,\mathrm{eV}$ and $\sigma_{\Sigma_{\rm{IH}}}=0.033\,\mathrm{eV}$. We discussed the gains provided by this source of cosmological information in breaking degeneracies between $\Sigma$ and other parameters, such as $\sigma_{8}$, which has a qualitatively similar effect as $\Sigma$ on non-linear power. We also illustrated the relative advantages in constraining the IH over the NH (72% and 54% improvements respectively) when including CMB information.

This distinction was found to be more nuanced when we expanded our analysis to include marginalization over the baryon parameters. Studying power responses illustrated specific degeneracies, e.g., between $\Sigma$ and $\eta_{0}$ on non-linear scales, that contributed to the approximate doubling of mass sum error forecasts. While stronger non-linear sensitivity to $\Sigma$ for the IH resulted in a less severe impact than for the NH, when CMB priors were accounted for the degradation for the NH (21%) was appreciably less than that for the IH (28%). Despite this, the IH once again exhibited a much greater overall improvement than for the NH.

Our weak lensing analysis also included an investigation into the impact of including (and marginalizing over) intrinsic alignments in our model on parameter forecasts. Unlike the significant bias induced in parameter estimations by neglecting IA, we found that there is a relatively minimal effect for the separate issue of forecasts because $C_{\ell}$ responses, which determine Fisher information, are very similar in both scenarios (especially for correlations between high redshift bins).

As the WL+CMB errors including baryons represent only a 1$\sigma$ (2$\sigma$) confidence of measuring $\Sigma_{\rm{NH}}$ ($\Sigma_{\rm{IH}}$) and are greater than $|\Delta\Sigma_{\rm{min}}|=0.04\,\rm{eV}$, we included independent forecasts from a Euclid-like spectroscopic galaxy clustering survey. We showed that this halved the baryon degradation factor for both hierarchies, with the NH still experiencing approximately half the impact of the IH. The galaxy clustering analysis highlighted the importance of specifying the mass hierarchy when making neutrino forecasts as most information was derived from linear matter power responses. The configuration of free-streaming scales corresponding to each neutrino mass played a non-negligible role on these scales. Ultimately, the combination of WL+GC+CMB reduced NH (IH) errors to close to 2$\sigma$ (3$\sigma$) and below the $|\Delta\Sigma_{\rm{min}}|$ threshold. We noted the challenge of navigating the variation in forecast methodologies, and therefore results, made across the literature (e.g., Carbone et al. 2011; Audren et al. 2013; Boyle & Komatsu 2018; Sprenger et al. 2018). We argue that our approach of simultaneously incorporating the effects of neutrinos and baryons in the weak lensing power spectrum on non-linear scales in a physically well-motivated manner makes our forecasts reliable and robust.

A potentially encouraging finding was that additional information from e.g., neutrinoless double beta decay on the mass sum, or improved constraints on the spectral index from forthcoming Stage IV CMB experiments, could substantially reduce the errors on $\Sigma$. If these priors prove sufficient to increase the Fisher information on their corresponding parameters by 10-50%, it would achieve a 2$\sigma$ (5$\sigma$) detection of $\Sigma$ for the NH (IH) and distinguish between the hierarchies at the 2$\sigma$ or even 3$\sigma$ level. It is important to emphasize that this conclusion depends on the measured value of the mass sum satisfying $\Sigma<0.1\,\rm{eV}$. Otherwise, both the NH and IH would be viable regardless of the strength of the error constraint.

Finally, we explored the issue of model bias. We used first order approximations (Taylor et al. 2007) to determine the bias in $\Sigma$ from miscalibrations of the baryon and intrinsic alignment parameters. We found that large biases in $A_{B}$ or $a_{\rm{IA}}$ induce a bias in $\Sigma$ larger than $\sigma_{\Sigma}$. By contrast, we found that even significant biases in the baryon feedback parameter, $\eta_{0}$, did not change the $\Sigma$ estimate by more than $\sigma_{\Sigma}$. Though an understanding of model bias is important to place our analysis in its proper context, it should be noted that it will ultimately require additional baryon information, improved simulations or self-calibrating methods for IA to successfully mitigate the issue.

An important limitation of our analysis to recognize is that the Fisher formalism underpinning it is not ideal for making forecasts in the case of distinct, competing models, as is the case with the normal and inverted mass hierarchies. Fisher information is determined by derivatives of likelihood functions around most-likely values. In cosmology these usually take the form of multi-variable unimodal Gaussians. In this paper we have made assumptions that this is the case when analyzing each hierarchy separately. However, to make a robust forecast of $\Sigma$, given uncertainty over which hierarchy is correct, our analysis should incorporate the resulting bimodal Gaussian (with peaks at $\Sigma_{\rm{NH,min}}=0.06\,\rm{eV}$ and $\Sigma_{\rm{IH,min}}=0.1\,\rm{eV}$). This would require a non-trivial extension of the Fisher formalism, which is beyond the scope of this work but that we intend to pursue in the near future.

In summary, by using the M15 and M16 prescriptions for baryonic and neutrino effects on haloes, we have been able to perform a full forecast analysis of the neutrino mass sum in the normal and inverted hierarchies. By combining multiple large scale structure probes through Euclid-like weak lensing and spectroscopic galaxy clustering surveys, Planck CMB constraints and BOSS low redshift galaxy clustering data, we have shown that degradation due to baryon feedback can be reduced to between 9-16%. Using these sources of information, and accounting for baryons and intrinsic alignments, Stage IV surveys could be expected to measure $\Sigma$ with 1-$\sigma$ errors of $\sigma_{\Sigma_{\rm{NH}}}=0.034\,\rm{eV}$ and $\sigma_{\Sigma_{\rm{IH}}}=0.034\,\rm{eV}$ for the NH and IH. These approach the confidence level required to meaningfully distinguish the hierarchies, but with additional future priors on $\Sigma$ and $n_{s}$ there is tentative optimism for achieving more definitive results.

Acknowledgements

The authors would like to thank Alexander Mead for several useful discussions about the implementation of massive neutrinos in HMCODE, and John Peacock for valuable input on the approach to galaxy clustering used in this paper. DC acknowledges the support of an STFC studentship. ANT thanks the Royal Society for a Wolfson Research Merit Award, and the STFC for support from a Consolidate Grant. AH is supported by an STFC Consolidated Grant.

Appendix A The halo model

The halo model provides a highly useful framework for interpreting the behaviour of the matter distribution on non-linear scales (Seljak 2000; Peacock & Smith 2000) (note that several fundamentals of the halo model we outline here overlap with a primer in Copeland et al. (2018)). The prescription of spherically-collapsed haloes positioned randomly in the linear density field permits the matter power spectrum, which describes the two-point statistics of density perturbations, to be written as the sum of two separate contributions,

[TABLE]

Here we have used the dimensionless form of the power spectrum, $\Delta^{2}\left(k\right)$. Correlations between different haloes are represented by the 2-halo term, $\Delta_{2h}^{2}\left(k\right)$, which is commonly approximated as the linear matter power spectrum because it is relevant on large scales. Internal halo structure is described by the 1-halo term, $\Delta_{1h}^{2}\left(k\right)$. This is constructed by summing the self-convolutions of haloes with mass $M$ over the full mass range,

[TABLE]

The integrand is weighted by the total number density of halo pairs with mass $M$. The halo mass function for the comoving number density of haloes within $\mathrm{d}M$ is provided by $n\left(M\right)$ while $\bar{\rho}$ is the background matter density. Press & Schechter (1974) showed that the mass function is nearly independent of cosmology if expressed in terms of the density peak height, $\nu\equiv\delta_{c}/\sigma\left(M\right)$, given by the spherical collapse overdensity, $\delta_{c}$, when the density field is smoothed at the mass scale, $M$. Sheth & Tormen (1999) determined a form of the ‘universal’ mass function,

[TABLE]

that fit simulations to be

[TABLE]

This mass function depends on the square root of the variance of the matter distribution,

[TABLE]

at scale, $R$, defining the radius of a sphere enclosing homogeneous mass, $M$. The field is smoothed by the window function of the spherical top-hat profile. The dimensionless linear matter power spectrum is equivalent to the fractional contribution to $\sigma^{2}\left(R\right)$ per logarithmic interval of $k$, and therefore indirectly influences the 1-halo power.

When normalized by its mass, the Fourier space density profile of a spherical halo is given by

[TABLE]

which transforms a real-space density profile, such as the common NFW case (Navarro et al. 1997),

[TABLE]

The inner (linear) and outer (cubic) regions of the halo are separated at the scale radius, $r_{s}$, which also determines the normalization density, $\rho_{s}$. As a halo stops collapsing once it reaches a virial equilibrium of gravitationally interacting matter shells exchanging energy, the integral in equation (15) is truncated at the associated virial radius,

[TABLE]

This quantity and the related virial density, $\Delta_{v}$, which depends on the specific cosmology in which spherical collapse occurs, are useful parameters for characterizing the final, virialized halo. The relationship between the scale radius and virial radius, $r_{s}=r_{v}/c$, is determined through the Bullock et al. (2001) concentration factor relation,

[TABLE]

which is sensitive to both cosmology and halo mass via the formation redshift, $z_{f}$ at which a fraction $f$ of the total matter in a density perturbation collapses. Bullock et al. (2001) find $f=0.01$ produces the most accurate fits of their model to N-body simulations. The amplitude, $A_{B}$, is discussed in detail in § 2.2.

Appendix B Probes of Large Scale Structure

In this work we explore the impact that marginalizing over baryons and the amplitude of the intrinsic alignment signal has on the capacity of Stage IV surveys to measure $\Sigma$ and distinguish neutrino mass hierarchies. To do so we perform Fisher analyses using weak lensing and spectroscopic galaxy clustering as signals. Here we briefly review the theory behind these probes and the Fisher formalism used to make forecasts.

B.1 Weak gravitational lensing

We approximate the source galaxy distribution as

[TABLE]

where $z_{0}=0.636$, which roughly matches that expected for a Euclid-like survey (Laureijs et al. 2011)555The weak lensing formalism presented here has also been detailed in Copeland et al. (2018), except for intrinsic alignments which were not discussed in that work.. Within $\mathrm{d}\chi$ of comoving position, $\chi$, the number of source galaxies is $\mathrm{d}\chi n\left(\chi\right)$. This is an important ingredient in the weak lensing convergence power spectrum, which, in the Limber approximation is

[TABLE]

where $\left(i,j\right)$ are tomographic redshift bin labels, $f_{K}\left(\chi\right)$ is the comoving angular distance, and $g_{i}\left(\chi\right)$ is the total weighting function, determined by the separation between sources and lenses over the distribution up to the survey limit, $\chi_{\rm{\max}}$. This is given by

[TABLE]

It is necessary to address photometric redshift errors. Within each redshift bin the galaxy distribution when normalized over the survey is given by

[TABLE]

where $p_{ph}\left(z_{ph}|z\right)$ is the probability of galaxies at $z$ being measured at redshift $z_{ph}$. The integrals are equivalent to convolving the true distribution with a binning top-hat function, whilst accounting for the probability of galaxies being measured at $z_{ph}$ outside the bin. The denominator is the result of normalizing over the bin. This can be shown to be equivalent to the distribution used in Taylor et al. (2007). We use the frequently employed probability distribution (Ma et al. 2006; Taylor et al. 2007),

[TABLE]

For the photometric redshift error, we use the Laureijs et al. (2011) value of $\sigma_{z}=0.05$.

B.1.1 Intrinsic alignments

The tidal forces that correlate ellipticities, $\epsilon$, and shears, $\gamma$, are an intrinsic contribution to the observed signal. Accounting for this to linear order by writing the observed ellipticity as $\epsilon=\gamma+\epsilon^{I}$, the weak lensing power spectrum becomes

[TABLE]

where $C_{\ell,ij}^{\gamma\gamma}$ is given by equation (20) and the term, $C_{\ell,ij}^{I\gamma}$, encompasses correlations of foreground ellipticities with background shear and of foreground shear with background ellipticities. The latter correlations should be very small under the assumption that redshift bins do not overlap. The IA-shear term is given by

[TABLE]

and accounts for terms arising from the misordering of redshift bins as well as those from a physical signal. The IA-IA term is given by

[TABLE]

The matter-IA and IA-IA source power spectra in the non-linear alignment model (Bridle & King 2007) are

[TABLE]

and

[TABLE]

$D\left(z\right)$ is the linear growth factor.

The amplitude, $a_{\rm{IA}}$, can be taken as a free parameter with a fiducial value of unity, while fits to SuperCOSMOS galaxy data determine the normalization, $c_{\rm{IA}}\approx 0.0134/\rho_{\rm{crit}}$, in terms of the critical density (Brown et al. 2002; Hirata & Seljak 2004). In Figure 2 we show comparisons of the shear-shear, IA-shear and IA-IA auto- and cross- lensing power spectra for several combinations of redshift bins. The IA-IA contributions to the signal are most significant for low redshift auto-power spectra, while the IA-shear term has its greatest impact on cross-power spectra of widely separated bins because it describes the correlations of foreground galaxy ellipticities with background shear.

B.2 Spectroscopic galaxy clustering

The clustering of galaxies provides an additional probe through the anisotropic galaxy power spectrum, $P_{gg}\left(k,\mu;z\right)$, for which the specific observable is H$\alpha$-emitting galaxies for a Euclid-like survey. We outline an approach to constructing $P_{gg}\left(k,\mu;z\right)$ that closely follows e.g., Ballinger et al. (1996) and Seo & Eisenstein (2003). Sources of apparent anisotropy make it necessary to include the dependence on the cosine of the wave mode and radial vector, $\mu$. These include redshift-space distortions arising from a peculiar velocity component in observed redshifts (Kaiser 1987). As a result, the observed strength of galaxy clustering will be sensitive to the orientation of the pair-separation vector. To account for this, the Kaiser factor is incorporated into $P_{gg}\left(k,\mu;z\right)$. The underlying dark matter distribution is biased by the galaxies that trace it (see e.g., Desjacques et al. 2018), so the galaxy power spectrum can then be constructed from the matter power spectrum according to

[TABLE]

where $b\left(z\right)=\sqrt{1+z}$ is the effective bias of the sample of H$\alpha$-emitters used in Amendola et al. (2018) and $f\left(z\right)$ is the logarithmic growth rate of structure, which relates the peculiar velocity to the density. The redshift-dependence of the (linear) matter power spectrum is removed by factoring out $\sigma_{8}\left(z\right)$, which leads to $b\sigma_{8}\left(z\right)$ and $f\sigma_{8}\left(z\right)$ becoming free parameters in galaxy clustering forecasts (White et al. 2009).

Converting measured angles and redshifts into radial and transverse distances requires assumptions about cosmology. If these are incorrect, distortions arise through the AP effect. The ratios,

[TABLE]

and

[TABLE]

of the assumed and true angular diameter distance and Hubble parameter (Alcock & Paczynski 1979) change wavenumbers and angles from their ‘fiducial’ values to (Ballinger et al. 1996)

[TABLE]

and

[TABLE]

Under these projections the matter power spectrum transforms as

[TABLE]

Hereafter, unless explicitly stated, $k$ and $\mu$ may be assumed to account for the AP effect.

Finally, redshift uncertainties along the line-of-sight are accounted for via the factor (e.g., Wang et al. 2013),

[TABLE]

where

[TABLE]

with the spectroscopic redshift error given by $\sigma_{z=0}=0.001$ (Laureijs et al. 2011). The full observed linear galaxy power spectrum is then

[TABLE]

The final term represents any residual shot noise that is not accounted for by Poisson sampling the underlying CDM density field with H$\alpha$-emitting galaxies. We assign a fiducial value of $P_{\rm{shot,fid}}=0$ for all $z$. Non-linear contributions, such as the Finger-of-God effect (Hamilton 1998) or the damping of the BAO component of $P\left(k;z\right)$ (Eisenstein et al. 2007; Seo & Eisenstein 2007), are beyond the scope of this work and are not included.

B.3 Fisher formalism

The Fisher information matrix is given by the expectation of the curvature of a multi-dimensional likelihood function, $L\left(\mn@boldsymbol{x}\lvert\Theta\right)$, evaluated at the most likely values of a parameter set, $\Theta=\left(\theta_{1},...,\theta_{N}\right)$, for the data set, $\mn@boldsymbol{x}$ (e.g., Tegmark, Taylor & Heavens 1997). These values are taken to be the fiducial ones in our model. The inverse of this matrix provides the parameter covariances that underpin our survey forecasts. Therefore, we write the Fisher matrix as

[TABLE]

where $\mathcal{L}=-\ln L$. Assuming a Gaussian likelihood,

[TABLE]

where $\mn@boldsymbol{C}$ is the data covariance matrix and $\boldsymbol{\mu}=\langle\mn@boldsymbol{x}\rangle$ is the mean data vector, the result,

[TABLE]

can be derived. Here, $M_{\alpha\beta}=\boldsymbol{\mu}_{,\alpha}\boldsymbol{\mu}_{,\beta}^{T}+\boldsymbol{\mu}_{,\beta}\boldsymbol{\mu}_{,\alpha}^{T}$ is the expectation value of the second derivative of the data matrix, $\left(\mn@boldsymbol{x}-\boldsymbol{\mu}\right)\left(\mn@boldsymbol{x}-\boldsymbol{\mu}\right)^{T}$, under the Gaussian approximation. Derivatives with respect to parameters are written via the shorthand $,\alpha\equiv\partial/\partial\Theta_{\alpha}$.

We apply the Fisher formalism for two observables: the weak lensing convergence power spectrum and the galaxy power spectrum, which are both treated as Gaussian. For the former we have (Tegmark et al. 1997; Takada & Jain 2004)

[TABLE]

which sums over the spherical harmonic $\ell$ and $m$ modes. The fraction of the sky within the scope of the survey, $f_{\rm{sky}}$, is included as a limiting factor. The auto- and cross-correlations of observed power in redshift bins, $i,j,p,q=\left(1,...,N_{z,\rm{WL}}\right)$, are used to construct the covariance matrix,

[TABLE]

If we ignore non-linear couplings between modes (see e.g., Kiessling et al. (2011) for explorations of how to account for these additional correlations within the Fisher formalism), then the covariance matrix can be treated as block diagonal in $\ell$. Shape noise, averaged over the number of galaxies, is added to the auto-correlations of power within each bin such that

[TABLE]

where $n_{i}$ is the number density of galaxies in redshift bin, $i$, and $\sigma_{e}=0.3$.

The Fisher matrix for galaxy clustering (Seo & Eisenstein 2003),

[TABLE]

is computed by integrating over the angular and radial contributions of the logarithmic galaxy power derivatives across $N_{z,\rm{GC}}$ redshift bins. An upper bound of $k_{\rm{max}}=0.2\,h^{-1}\mathrm{Mpc}$ is chosen, as we restrict our GC forecasts to linear scales. The choice of lower bound has very limited bearing on results so we use the common value, $k_{\rm{min}}=0.001\,h^{-1}\mathrm{Mpc}$. Cross-correlations between redshift bins will have limited impact so they are not included. The effective volume,

[TABLE]

encodes the covariances and is calculated from the survey volume, $V_{s}$, and the number density, $n\left(z\right)$, of H$\alpha$-emitters, for which we use Euclid survey values provided by the reference case in Table 3 of Amendola et al. (2018).

Note that all scales and angles refer to the ‘fiducial’ case, for which we have dropped the label for convenience. Free parameters governing the AP effect, growth, bias and residual shot noise within each redshift bin independently, $\left\{D_{A}\left(z_{i}\right),H\left(z_{i}\right),f\sigma_{8}\left(z_{i}\right),b\sigma_{8}\left(z_{i}\right),P_{\rm{shot}}\left(z_{i}\right)|i=1,...,N_{z,\rm{GC}}\right\}$, can then be propagated in a Fisher analysis alongside redshift-independent parameters, $\left\{\omega_{m},\omega_{b},\omega_{\nu},h,n_{s}\right\}$, controlling the shape of $P\left(k\right)$. The physical densities are related to our original density parameters via $\left\{\omega_{i}=\Omega_{i}\,h^{2}|i=m,b,\nu\right\}$. We marginalize over the nuisance parameters, $\left\{b\sigma_{8}\left(z_{i}\right),P_{\rm{shot}}\left(z_{i}\right)|i=1,...,N_{z,\rm{GC}}\right\}$, by inverting the Fisher matrix, removing rows and columns corresponding to these parameters, and then inverting the remaining sub-block (Seo & Eisenstein 2003). The resulting Fisher matrix, $\tilde{F}$, treats $\left\{D_{A}\left(z_{i}\right),H\left(z_{i}\right),f\sigma_{8}\left(z_{i}\right)|i=1,...,N_{z,\rm{GC}}\right\}$ as free parameters contributing information alongside $\left\{\omega_{m},\omega_{b},\omega_{\nu},h,n_{s}\right\}$. Collecting all these parameters into the set, $\tilde{\Theta}$, we derive a final Fisher matrix for the $w$CDM parameter set, $\Theta=\left\{\Omega_{m},\Omega_{b},\Omega_{\nu},h,n_{s},\sigma_{8},w\right\}$, by performing the transformation,

[TABLE]

where $J=\partial\tilde{\Theta}/\partial\Theta$ is the Jacobian matrix of parameter derivatives. Note that because $\tilde{\Theta}$ is evaluated in the $w$CDM model it can be completely determined by $\Theta$.

We evaluate the Fisher information available to a Euclid-like spectroscopic galaxy survey in the redshift range $0.75<z\leq 2.05$, with bin width $\Delta z=0.1$. For lower redshifts we refer to BOSS as it provides superior information at these $z$ than Euclid-like surveys are designed to achieve. This entails using covariance information provided by Alam et al. (2017) on $\left\{D_{A}\left(z_{i}\right),H\left(z_{i}\right),f\sigma_{8}\left(z_{i}\right)\right\}$, where the $i$ label runs over three partially overlapping redshift bins, spanning the total range $0.2<z\leq 0.75$. We perform a matrix inversion to derive the corresponding Fisher information matrix. The BOSS bins overlap with the full $z$ range available to the Euclid-like survey, so to ensure that covariance between the surveys is limited we only use redshift bins for the latter that lie beyond the range of BOSS. We then evaluate our final galaxy clustering Fisher matrix from the combination of Euclid-like high $z$ forecasts and BOSS low $z$ data.

When evaluating derivatives, step sizes are determined for each parameter as 1% of the fiducial value for that parameter. We have established that modest changes in the step size, up to doubling it, has no significant impact, suggesting our derivatives are stable. Similarly, the results are insensitive to whether we use 2-point or 5-point derivatives. For weak lensing and galaxy clustering, confidence regions for pairs of parameters are determined by inverting the corresponding Fisher matrix to marginalize over the remaining parameters, and extracting the resulting errors, $\sigma_{\alpha}=\sqrt{\left[F^{-1}\right]_{\alpha\alpha}}$, from the elements. Specifically, we compute the 1-$\sigma$ errors.

Appendix C Complete Fisher Forecasts

In Figure 12 we present the results of the full weak lensing Fisher analysis for every parameter combination in $\Theta=\left(A_{B},\eta_{0},\Omega_{m},\Omega_{b},\Sigma,h,\sigma_{8},n_{s},w,a_{\rm{IA}}\right)$ for the NH. These include plots showing comparisons of cases where baryon parameters are fixed to their fiducial values to those where they are marginalized over.

Appendix D Intrinsic alignment power spectra responses

In Figure 13 we provide a comparison of the $C_{\ell}$ responses to varying parameters, both with and without intrinsic alignments. To display this clearly, we show results for a redshift bin combination, $ij=09$, at which the differences between the two cases are most pronounced.