# Neutrino mass constraints beyond linear order: cosmology dependence and systematic biases

Aoife Boyle, Fabian Schmidt

► To cite this version:

Aoife Boyle, Fabian Schmidt. Neutrino mass constraints beyond linear order: cosmology dependence and systematic biases. Journal of Cosmology and Astroparticle Physics, 2021, 04, pp.022. 10.1088/1475-7516/2021/04/022. [hal-03070688]

HAL Id: hal-03070688

https://hal.science/hal-03070688v1

Submitted on 1 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Neutrino mass constraints beyond linear order: cosmology dependence and systematic biases

To cite this article: Aoife Boyle and Fabian Schmidt JCAP04(2021)022

View the [article online](#) for updates and enhancements.

You may also like

- Detecting Neutrino Mass by Combining Matter Clustering, Halos, and Voids
Adrian E. Bayer, Francisco Villaescusa-Navarro, Elena Massara et al.

- Study of Momentum Diffusion with the Effect of Adiabatic Focusing
J. F. Wang and G. Qin

- WEIGHING NEUTRINOS WITH COSMIC NEUTRAL HYDROGEN
Francisco Villaescusa-Navarro, Philip Bull and Matteo Viel
# Neutrino mass constraints beyond linear order: cosmology dependence and systematic biases

Aoife Boyle$^{a,b}$ and Fabian Schmidt$^a$

$^a$Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching bei München, Germany

$^b$CNRS & Sorbonne Université, UMR 7095, Institut d'Astrophysique de Paris, 98 bis Boulevard Arago, 75014 Paris, France

E-mail: boyle@iap.fr, fabians@mpa-garching.mpg.de

Received December 7, 2020

Accepted February 18, 2021

Published April 9, 2021

**Abstract.** We demonstrate the impact on forecasted neutrino mass constraints of improving galaxy clustering and CMB lensing predictions from linear to next-to-leading-order power spectra. The redshift-space 1-loop power spectrum model we adopt requires an additional four free bias parameters, a velocity bias parameter and two new stochastic parameters. These additional nuisance parameters appreciably weaken the constraints on $M_ν$. CMB lensing plays a significant role in helping to alleviate these degeneracies and tighten the final constraints. The constraint on the optical depth to reionisation $\tau$ has a strong effect on the constraint on $M_ν$, but only when CMB lensing is included in the analysis to keep the degeneracies with the nuisance parameters under control. We also extract constraints when 1) using the BAO signature only as a distance probe, and 2) isolating the scale-dependence of the power spectrum, which, as shown in previous work, provides a cosmology-independent probe of $M_ν$. All constraints except the latter remain strongly sensitive to the assumption of a flat $\Lambda$CDM universe. We perform an analysis of the magnitude of the shift introduced in the inferred $M_ν$ value when neglecting nonlinear corrections, and show that, for a Euclid-like survey, this shift becomes roughly equal to the $1\sigma$ constraint itself even with a conservative cut-off scale of $k_{max} = 0.1h \text{ Mpc}^{-1}$. We also perform a calculation of the approximate expected bias in neutrino mass caused by not including the next, 2-loop order and expect a shift of only about 20% of the $1\sigma$ error for $k_{max} = 0.2h \text{ Mpc}^{-1}$ in a Euclid-like survey.

**Keywords:** cosmological parameters from LSS, neutrino masses from cosmology, power spectrum, redshift surveys

**ArXiv ePrint:** 2011.10594

<figure><img src="image_2.png" /></figure>
# Contents

<table><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 Methodology</td><td>3</td></tr><tr><td>2.1 Implementation of the NLO power spectra</td><td>3</td></tr><tr><td>2.2 Fiducial cosmology, surveys and priors</td><td>5</td></tr><tr><td>2.3 Constraining methods</td><td>6</td></tr><tr><td>3 Results</td><td>7</td></tr><tr><td>3.1 Combined constraints</td><td>7</td></tr><tr><td>3.2 Free-streaming constraints</td><td>11</td></tr><tr><td>3.3 BAO constraints</td><td>12</td></tr><tr><td>4 Discussion</td><td>13</td></tr><tr><td>4.1 Selection effects</td><td>13</td></tr><tr><td>4.2 Varying the fiducial neutrino mass</td><td>14</td></tr><tr><td>4.3 Neutrino-induced bias</td><td>15</td></tr><tr><td>4.4 Degeneracy between bias parameters</td><td>15</td></tr><tr><td>4.5 Comparisons with previous work</td><td>15</td></tr><tr><td>4.6 Biases on M<sub>ν</sub> from the power spectrum model</td><td>16</td></tr><tr><td>4.7 Alternative galaxy surveys</td><td>17</td></tr><tr><td>5 Conclusions</td><td>18</td></tr><tr><td>A Survey parameters</td><td>19</td></tr></table>

# 1 Introduction

The promise of improved constraints on the sum of neutrino masses $M_ν$ is one driving force behind planned large-scale structure surveys [e.g. 1–5] and CMB experiments [e.g. 6, 7]. At present, solar and atmospheric oscillation experiments impose a minimum $M_ν$ (for the normal hierarchy, which is currently favoured at $3σ$ [8, 9]) of just under 0.06 eV. The upper limits provided by cosmology still significantly outstrip all experimental alternatives, with the 2018 Planck data release asserting a $2σ$ constraint<fnref n="1" /> of $M_ν < 0.12$ eV by combining CMB and baryon acoustic oscillation (BAO) data [10] for a flat $Λ$CDM+$M_ν$ cosmology. However, the constraints on the neutrino mass from galaxy clustering and lensing measurements are inherently complex, and are limited by a number of factors including cosmological degeneracies [11–15], baryonic effects [16] and theoretical uncertainties related to the small-scale power spectrum [17–19].

In previous work [13, 14], we have attempted to deconstruct forecasted neutrino mass constraints from future galaxy clustering and CMB surveys. Our aim has been to assess the robustness of expected constraints, considering in particular how they could be weakened by

<footnote n="1">¹Note that the Planck constraints are quoted as a $2σ$ upper limit, but the forecasted constraints we present in this paper are $1σ$ uncertainties.</footnote>
possible alterations to flat $\Lambda$CDM. We have demonstrated that constraints on the total neutrino mass, $M_\nu$, are heavily affected by possible small deviations in the assumed cosmology, even with powerful next-generation surveys. However, we have also shown that isolating the effects of neutrino free-streaming in the power spectrum as a constraining tool can provide cosmology-independent constraints.

In the previous work cited above, we based our forecasts entirely in the linear regime. We assumed the Kaiser formula to describe the tree-level redshift-space galaxy power spectrum (see section 2 for details). This formulation is accurate for quite a limited range of potentially usable scales, and takes a rather crude approach to galaxy bias and redshift-space distortions, and neglects any selection effects. We derived our constraint forecasts for $M_\nu$ assuming a minimum usable comoving scale $k_{max} = 0.2 \ h \ \text{Mpc}^{-1}$, but this is likely to be optimistic. Moving beyond the assumption of linearity provides a natural and important extension to our work. This is particularly important because a substantial contribution to the constraint on $M_\nu$ is provided by the suppression in the power spectrum on small scales, where non-linear corrections become increasingly significant.

Recently, ref. [20] presented the first complete expression for the perturbative next-to-leading-order (NLO) or 1-loop power spectrum in redshift space including all selection effects (see [21–24] for related derivations). Using the approach developed in [25], the formalism presented orders the contributing terms in a consistent way and keeps the contributions appropriately general. It also includes the leading higher-derivative bias and velocity bias terms. This rigorous approach leads to a complex but complete final expression for the redshift-space power spectrum that makes minimal assumptions. We implement this formalism here. Ref. [20] emphasises the importance of accounting for line-of-sight-dependent selection effects in the calculation of the galaxy power spectrum unless physical arguments can be made for neglecting them. We neglect selection effects for most of this article, but justify this and provide some extended calculations in section 4.1. Another recent work, ref. [26], presented a Fisher analysis of the cosmological information available in the 1-loop power spectrum and tree-level galaxy bispectrum, after fully accounting for selection effects, for future galaxy surveys like Euclid [1, 2]. Here we present a similar analysis, but based on the power spectrum only and focused on the neutrino mass, and extending to consider a range of possible cosmologies. We discuss our results in the context of recent related work, such as [19], in section 4.5.

We proceed as follows. In section 2, we discuss the implementation of the NLO galaxy and CMB lensing power spectra in our Fisher matrix code, our fiducial cosmology and surveys and our calculation methods. In section 3 we provide results from the linear and NLO power spectra for three scenarios: using the full broadband power spectrum (we call this the ‘combined’ calculation), using the BAO feature only, and using only the scale-dependent suppression of the power spectrum unique to neutrinos (which we call the ‘free-streaming’ calculation). We also examine how CMB lensing can improve these constraints. We present all of our results in the form of forecasts for Euclid [1, 2], complemented with the expected CMB information from Simons Observatory [6], but discuss the changes in the results for other surveys in section 4.7. In section 4 we provide some discussion and caveats, along with comparisons to previous work in the literature. We conclude in section 5.
# 2 Methodology

## 2.1 Implementation of the NLO power spectra

One of the goals of this work is to compare the results obtained when performing forecasts using the linear and next-to-leading-order (NLO) or 1-loop galaxy power spectrum. The Kaiser formula describes the linear redshift-space galaxy power spectrum at redshift $z$:

$$P_{gg,s}^{\text{lin.}}(k, z) = [b_1(z) + f_{bc}(k, z)\mu^2]^2 P_{bc}(k, z) + P_{\epsilon}^{\{0\}}(z). \quad (2.1)$$

Here, $P_{\epsilon}^{\{0\}}$ is the stochastic (shot noise) contribution, $P_{bc}(k, z)$ is the linear baryon and cold dark matter power spectrum (obtained using the Boltzmann code CLASS [27]), while $f_{bc}(k, z) \equiv d \ln D(k, z)/d \ln a$ is the scale-dependent linear growth rate. It has been verified by simulations [e.g. 28] that because neutrinos do not cluster on the scales on which galaxies form, the baryon and cold dark matter only power spectrum $P_{bc}(k, z)$ should be used in the calculation of the galaxy power spectrum, along with its corresponding growth rate, $f_{bc}(k, z)$. The presence of massive neutrinos induces a scale-dependence in $f_{bc}(k, z)$ and this contributes to a small degree to constraints on $M_{\nu}$ [13, 29] and is therefore included in equation (2.1). In the following, we will drop the redshift argument for clarity, noting that all bias parameters, power spectra, as well as the growth rate $f_{bc}$ are to be evaluated at the redshift of observation.

The NLO redshift-space power spectrum (which includes the linear component above) provided by [20] is calculated using a total of 28 independent loop integrals, and requires specification of 5 additional free bias parameters, 1 velocity bias parameter and 3 stochastic parameters if selection effects are neglected, as they are throughout most of this work. We note that the effects of massive neutrinos were not explicitly considered in [20], and our results are therefore based on the assumption that the important effects of massive neutrinos are captured by their contribution to the linear power spectrum and growth factor (specifically of the CDM and baryon components only, as described above). This approximation was first introduced in [30, 31]. In section 4.3, we discuss the choice to neglect the specific effect of massive neutrinos on the bias parameters. The NLO redshift-space galaxy power spectrum can be summarised in the following form at a specific redshift:

$$P_{gg,s}^{\text{NLO}}(k, \mu) = P_{gg,s}^{\text{lb+hd}} + P_{gg,s}^{2-2}(k, \mu) + 2P_{gg,s}^{1-3}(k, \mu) \quad (2.2)$$

$$P_{gg,s}^{\text{l+hd}}(k, \mu) = [b_1 + f_{bc}(k)\mu^2]^2 P_{bc}(k) + P_{\epsilon}^{\{0\}} \\ - 2\{b_1 C_{\nabla^2\delta} + \mu^2 f_{bc}(k)[C_{\nabla^2\delta} + b_1\beta_{\nabla^2v}] \\ + \mu^4 f_{bc}^2(k)\beta_{\nabla^2v}\}k^2 P_{bc}(k) \\ + k^2 P_{\epsilon}^{\{2\}} - \mu^2 k^2 P_{\epsilon\epsilon_n}^{\{2\}} \quad (2.3)$$

$$P_{gg,s}^{2-2}(k, \mu) = \sum_{n=0}^{4} \sum_{(m,p)} A_{n(m,p)}(f_{bc}(k), b_{\mathcal{O}}) \mathcal{I}_{mp}(k) \mu^{2n} \quad (2.4)$$

$$P_{gg,s}^{1-3}(k, \mu) = \sum_{l=0,2,4,6} \sum_{n=1}^{5} C_{n}^{1-3,l}(f_{bc}(k), b_{\mathcal{O}}) \mathcal{I}_{n}(k) P_{bc}(k) \mathcal{L}_{l}(\mu) \quad (2.5)$$

In eq. (2.3) we have defined the parameter combination $C_{\nabla^2\delta} = C_s^2 + b_{\nabla^2\delta}$ of the higher-derivative bias proper and the effective sound speed $C_s^2$. The distinction between $C_{\nabla^2\delta}$ and $b_{\nabla^2\delta}$ will become relevant when also including the galaxy-matter cross-power spectrum as probed via the cross-correlation with CMB lensing.
The form of the 23 $\mathcal{I}_{mp}$ and 5 $\mathcal{I}_n$ loop integrals are provided in [20], and the coefficients $A_{n(m,p)}$ and $C_n^{1-3,l}$ and their input bias parameters $b_O$ can be obtained from a MATHEMATICA notebook provided as supplementary material.<fnref n="2" /> These coefficients significantly improve the calculation time by combining contributions to give the minimal number of integrals. We highlight again that we include the scale-dependence in $f_{bc}(k)$ when performing our calculations.

We also analyse the constraining power of CMB lensing. The CMB lensing power spectrum is calculated using the Limber approximation as

$$C_l^{\kappa\kappa} = \left(\frac{4\pi G \rho_{m,0}}{c^2}\right)^2 \int_0^{z_*} dz (1+z)^2 \left(\frac{d_A(z, z_*)}{d_A(z_*)}\right)^2 \frac{P_{mm}^{NLO}\left[k = \frac{l+1/2}{d_A(z)}, z\right]}{H(z)}. \quad (2.6)$$

$\rho_{m,0}$ is the comoving matter density, $z_*$ is the redshift of last scattering, and $d_A(z, z_*)$ represents the comoving angular diameter distance between the two redshifts. $P_{mm}^{NLO}$ is the NLO matter power spectrum (here including the mass contributed by neutrinos). It is calculated as $P_{mm}^{NLO} = P_L(k) + P_{22}(k) + 2P_{13}(k) - 2C_s^2 k^2 P_L(k)$, with $P_{22}(k)$ and $P_{13}(k)$ being calculated using a publicly available code<fnref n="3" /> [32]. Once again, we make the assumption that the effect of massive neutrinos is sufficiently well accounted for by including only the contribution of neutrino perturbations to the linear power spectrum. Calculating the matter power spectrum at the large number of redshifts required by equation (2.6) would be very time-consuming, but it is also improper to scale between different redshifts using the growth factor $D(z)$ because of the redshift-dependent shape of the power spectrum when massive neutrinos are included. To account for this, we calculate the NLO matter power spectra independently for $0 < z < 3$, but scale backwards using the growth factor beyond $z = 3$, where the suppression due to free-streaming is less substantial.

Finally, we can also consider the cross-correlation between galaxy clustering and CMB lensing as a probe. The galaxy-CMB lensing power spectrum is calculated (also using the Limber approximation) as

$$C_l^{g\kappa} = \left(\frac{4\pi G \rho_{m,0}}{c^2}\right) (1+z_g) \frac{d_A(z_g, z_*)}{d_A(z_g) d_A(z_*)} P_{gm} \left[k = \frac{l+1/2}{d_A(z_g)}, z_g\right], \quad (2.7)$$

The NLO galaxy-matter cross-power spectrum required in equation (2.7) can be derived straightforwardly from the NLO galaxy power spectrum presented in [20]. As we use the Limber approximation, we perform the calculation for $\mu = 0$ only. The components can be calculated as follows:

$$P_{gm,s}^{1+\text{hd}}(k) = b_1 P_L(k) + (b_1 C_s^2 - C_{\nabla^2\delta}) k^2 P_L(k) + k^2 P_{\epsilon,gm}^{\{2\}} \quad (2.8)$$

$$P_{gm,s}^{2-2}(k) = \sum_{(m,p)} A_{\text{cross(m,p)}} (f_{bc}(k), b_O) \mathcal{I}_{mp}(k) \quad (2.9)$$

$$P_{gm,s}^{1-3}(k) = \sum_{l=0,2,4,6} \sum_{n=1}^{5} C_{\text{cross,n}}^{1-3,l}(f_{bc}(k), b_O) \mathcal{I}_n(k) P_L(k) \mathcal{L}_l(0) \quad (2.10)$$

$A_{\text{cross(m,p)}}$ has 19 non-zero terms and $C_{\text{cross,n}}^{1-3,l}$ has 12 without selection effects. These coefficients are available in a Mathematica notebook alongside the supplementary material<fnref n="2" />

<footnote n="2">²https://github.com/djeong98/pkgs_supplement.</footnote>

<footnote n="3">³https://wwwmpa.mpa-garching.mpg.de/~komatsu/crl/list-of-routines.html.</footnote>
<table><caption>**Table 1.** Fiducial values assumed for the 6 bias parameters and 3 stochastic parameters required for calculation of the NLO galaxy power spectrum. $R_L$ denotes the Lagrangian radius of halos relevant at each redshift (see appendix A). Note that $C_{\nabla^2\delta}$ also includes the effective sound speed for matter, but this is a subdominant contribution.</caption><thead><tr><th>Bias parameter</th><th>Fiducial value</th></tr></thead><tbody><tr><td>b<sub>1</sub></td><td>b<sub>1</sub>(z)D(z) = 0.76</td></tr><tr><td>C<sub>∇<sup>2</sup>δ</sub></td><td>R<sub>L</sub><sup>2</sup></td></tr><tr><td>β<sub>∇<sup>2</sup>v</sub></td><td>R<sub>L</sub><sup>2</sup></td></tr><tr><td>b<sub>2</sub></td><td>0.412 - 2.143b<sub>1</sub> + 0.929(b<sub>1</sub>)<sup>2</sup> + 0.008(b<sub>1</sub>)<sup>3</sup></td></tr><tr><td>b<sub>K<sup>2</sup></sub></td><td>-<sup>2</sup>&frasl;<sub>7</sub>(b<sub>1</sub> - 1)</td></tr><tr><td>b<sub>td</sub></td><td><sup>23</sup>&frasl;<sub>42</sub>(b<sub>1</sub> - 1)</td></tr><tr><td>P<sub>ε</sub><sup>{0}</sup></td><td>n<sub>g</sub><sup>-1</sup></td></tr><tr><td>P<sub>ε</sub><sup>{2}</sup></td><td>0</td></tr><tr><td>P<sub>εεη</sub><sup>{0}</sup></td><td>0</td></tr></tbody></table>

provided in [20]. For the effective matter sound speed $C_s^2$, we take the value determined from N-body simulations in [33] at redshift zero of $1.31 h^{-2}$ Mpc$^2$, and rescale with redshift as $C_s^2(z) = [D(z)/D(0)]^4 C_s^2(0)$, which was found to describe the redshift dependence of the effective sound speed well. For the fiducial value of the cross-stochasticity between matter and galaxies, $P_{\epsilon, gm}^{\{2\}}$, for which no reliable simulation measurements exist to our knowledge, we take $C_s^2 P_{\epsilon}^{\{0\}}$. We stress however, that both $C_s^2$ and $P_{\epsilon, hm}^{\{2\}}$ have an entirely negligible effect on our results.

We perform all of our calculations using the Fisher matrix formalism (see [13] for details). Our code has been thoroughly validated through comparisons of results with those of [26]. The covariance for the three-dimensional galaxy power spectrum is calculated as:

$$ \langle \Delta P_{gg,s}(k, \mu)^2 \rangle = \frac{2\pi^2}{Vk^2\Delta k\Delta\mu} 2 [P_{gg,s}(k, \mu)]^2, \quad (2.11) $$

where $V$ is the volume of the redshift bin being observed, and $\Delta k$ and $\Delta \mu$ are the bin sizes for the wavenumber and angle with respect to the line of sight, respectively. The covariances for the angular power spectra are calculated as:

$$ \langle \Delta C_l^{xy} \Delta C_l^{mn} \rangle = \frac{1}{(2l+1)f_{\text{sky}}\Delta l} (C_l^{xm} C_l^{yn} + C_l^{xn} C_l^{ym}). \quad (2.12) $$

Here, $f_{\text{sky}}$ is the fraction of the sky observed. The $C_l$ values on the right-hand side of equation (2.12) must include appropriate noise terms for auto-correlation power spectra.

## 2.2 Fiducial cosmology, surveys and priors

In [13] we used Euclid [1, 2] as a sample survey for which to forecast constraints. In [14], we extended our work to include future CMB surveys and focused on Simons Observatory [6] as our example. Here we continue to focus on these two as representative surveys. The survey parameters for Euclid are given in table 2 in appendix A, and those for Simons Observatory are taken from table 1 of [6]. The Euclid survey area is 15000 deg$^2$ and the total number of galaxies included in the spectroscopic survey is approximately 50 million. The focus of our
work is not to provide quantitative information for a particular survey, but to demonstrate certain principles, and these should be generally extendable to all similar datasets. Results for alternative galaxy surveys are discussed in section 4.7. We assume a maximum usable wavenumber of $k_{max} = 0.2 h \text{ Mpc}^{-1}$ for all galaxy survey calculations for ease of comparison with results in our previous papers, and the effect of this choice is examined in section 4.6. We calculate a consistent $l_{max}$ value in each redshift bin for the galaxy-CMB cross-correlation calculations (eq. (2.7)) as $l_{max}(z) = k_{max} \times d_A(z) - 1/2$. For the CMB lensing auto-correlation (eq. (2.6)), no additional $l_{max}$ is implemented because of the natural limit imposed by the lensing noise. We have checked that repeating our CMB lensing calculations with full Halofit power spectra (as implemented in CLASS) has negligible impact on the results.

Our fiducial cosmology is identical to that in [13, 14]. We have a total of 11 free cosmological parameters. We marginalise our constraints on $M_\nu$ over the six standard Planck convention $\Lambda$CDM parameters ($\theta_s, \omega_b, \omega_{cdm}, \tau, A_s$ and $n_s$) plus $N_{\text{eff}}$. In all cases, we examine the effects of freeing $\Omega_k, w_0$ and $w_a$ on the constraints.

Table 1 shows the 9 free nuisance parameters that remain if selection effects are neglected and how their fiducial values were calculated. $R_L$ represents the Lagrangian radii of haloes, calculated using approximate typical halo masses at each redshift derived from the halo mass function of [34] based on an abundance-matching procedure. Specifically, we use the median mass of halos above a threshold adjusted to yield the expected number density of galaxies $n_g$ in each redshift bin as representative halo mass. The fiducial values of $C_{\nabla^2\delta}$ and $\beta_{\nabla^2\delta}$ are set to $R_L^2$ based on dimensional reasoning. We do not expect the precise values of the fiducial higher-derivative bias parameters to have any impact on our results. The choice for $b_1$ was made for consistency with [13, 14]. The fitting formula for $b_2$ was taken from N-body simulations and provided by [35]. The values of $b_{K2}$ and $b_{td}$ were both derived from the Lagrangian local (in the matter density) assumption [20].

Our basic CMB prior is also the same as that used in [14]. We combine the temperature power spectrum from Planck [10], a forecasted E-mode polarisation power spectrum from Simons Observatory and the cross-correlation between them. These are the unlensed spectra, as we choose to analyse the impact of including CMB lensing separately from that of the primary anisotropies. Note that when this is not of interest, one should perform these calculations using the lensed temperature and polarisation power spectra, as this is what is measured in real experiments. Delensing of measured spectra is a complicated process, and it has been demonstrated previously that the cross-correlation between lensed temperature/polarisation power spectra and the lensing power spectrum is very small [36, 37], so the results will not be misleading. We also include a prior on $\sigma(\tau) = 0.008$ from Planck (as given for TT,TE,EE+lowE in [10]).

## 2.3 Constraining methods

We perform forecasts in this article using three different approaches to the galaxy power spectrum. The most optimistic cosmological parameter constraint forecasts for a particular galaxy survey can be extracted by using the full galaxy power spectrum given in equation (2.2) as the observable with which derivatives for the Fisher matrix are calculated. To include full geometric information, the power spectrum must be adjusted to account for the uncertainty in $H(z)$ and $D_A(z)$ when converting observed scales to physical comoving scales, as follows:

$$P(k_{\parallel}^{\text{obs}}, k_{\perp}^{\text{obs}}) = \frac{H(z)}{H_{\text{fid}}(z)} \left( \frac{D_{A, \text{fid}}(z)}{D_A(z)} \right)^2 P(k_{\parallel}^{\text{com}}, k_{\perp}^{\text{com}}), \qquad (2.13)$$
where $k_{\parallel}^{\text{obs}} = k_{\parallel}^{\text{com}}(H_{\text{fid}}(z)/H(z))$ and $k_{\perp}^{\text{obs}} = k_{\perp}^{\text{com}}(D_A(z)/D_{A,\text{fid}}(z))$. These 'combined' constraints inherently combine information from many different sources within the power spectrum - constraints on distances from characteristic scales (e.g. BAO) and from the Alcock-Paczynski test, constraints on the structure growth rate as a function of redshift from redshift space distortions, and constraints from the overall shape and amplitude of the power spectrum itself. As shown in [13], the constraints extracted in this way can depend heavily on the underlying cosmology, degrading significantly when extensions like non-zero curvature or a varying dark energy equation of state are added.

However, our work in [13, 14] demonstrated that robust, cosmology-independent constraints on the neutrino mass can be extracted from the galaxy power spectrum by isolating a particular probe. Massive neutrinos suppress the growth of structure on small scales relative to that on large scales to an extent that is primarily dependent on the total mass of the neutrino species by free-streaming out of potential wells below a characteristic scale. This process imprints a characteristic scale-dependent signature in the shape of the matter power spectrum and in the structure growth rate. Because changes in other cosmological parameters cannot replicate this effect, the constraints extracted from this scale-dependent suppression in the power spectrum and structure growth rate (we call these 'free-streaming' constraints) are independent of the assumed cosmology, and therefore give much more reliable constraints than those extracted from the full combined power spectrum.

The details of how the free-streaming constraints are achieved in the linear case are given in [13]. In the NLO case, we normalise the input $P_{\text{bc}}(k)$ and $f_{\text{bc}}(k)$ at a fixed scale before calculating the components of equation (2.2) and then calculating numerical derivatives. Therefore only the scale-dependence in the input $P_{\text{bc}}(k)$ and $f_{\text{bc}}(k)$ contributes. Smoothing is also applied to remove any wiggles in the derivatives due to changes in the BAO signal. We marginalise the final Fisher matrix over the uncertainty in the amplitudes of $P_{\text{bc}}(k)$ and $f_{\text{bc}}(k)$. The free-streaming information can be combined with corresponding signatures in the CMB lensing and galaxy-CMB lensing power spectra by marginalising over the large-scale amplitudes of both.

Finally, we present constraints using distance information only extracted from the isolated BAO signal (fitted from the galaxy power spectrum). This proceeds in the same manner as outlined in [13]. A two-step Fisher matrix process is used to first constrain $H(z)$ and $D_A(z)$, then the cosmological parameters. For these constraints (but not for the combined constraints) we also implement an IR resummation procedure in our model to account for BAO damping. In section 3.3 we show that the impact of this on the constraints on $M_\nu$ is small enough to be excluded in the case of the combined constraints.

# 3 Results

## 3.1 Combined constraints

The left panel of figure 1 shows the forecasted constraints on $M_\nu$ for Euclid using the full galaxy power spectrum as a constraining tool, along with the priors outlined in section 2.2. The vertical dotted tick marks show the constraints predicted if the tree-level galaxy power spectrum is assumed instead. We see that the constraints are significantly weakened on updating our results with the NLO power spectrum. As found for the linear galaxy power spectrum in our previous work [13, 14], these combined constraints are strongly cosmology-dependent, and constraints assuming flatness or a simple cosmological constant should therefore be presented with appropriate caveats both in forecasts and when derived from real data.
<figure><img src="image_3.png" /></figure>

**Figure 1.** Forecasted 1σ constraints on $M_ν$ for various cosmological models, using the full NLO galaxy power spectrum for Euclid. This represents the most optimistic forecast and inherently includes BAO, RSD and shape information. The CMB prior combines the temperature power spectrum from Planck with a forecasted polarisation power spectrum from Simons Observatory (see section 2.2). The left panel shows the constraints with a 1σ prior on $τ$ of 0.008 (from Planck 2018 [10]), and the right panel shows how the constraints would improve if $τ$ were perfectly constrained through some external experiment. In both cases, the dotted tick marks show the constraints if linear theory was assumed for the galaxy power spectrum (with the same $k_{max}$ of 0.2 h Mpc$^{-1}$). Note that CMB lensing is neglected here but included in figure 4.

<figure><img src="image_4.png" /></figure>

**Figure 2.** The log derivatives of the NLO galaxy power spectrum (for $μ = 0.75$) with respect to $M_ν$ and the nine nuisance parameters at $z = 1.35$ (the central redshift of Euclid). A constant factor of 1.0 has been subtracted in the $b_1$ case and the stochastic parameter derivatives have been rescaled by a multiplicative factor (see the legend) for the readability of the plot.

The right panel shows the same constraints as the left panel but with $τ$ held fixed. It is interesting to note that the $M_ν-τ$ degeneracy, a significant focus of much recent forecast work on neutrino mass constraints from large-scale structure [e.g. 11–14, 18, 38, 39] becomes much less significant when considering the NLO galaxy power spectrum (note the change in the linear galaxy power spectrum constraints between the two panels). The $M_ν-τ$ degeneracy
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Contours showing the relationship between the different nuisance parameters and $M_ν$ for the central redshift bin of our Euclid model plus a CMB prior. The blue contours represent the ‘combined’ galaxy clustering case without CMB lensing and the pink contours represent the case when CMB lensing is added, reducing the dependence of the constraints on the nuisance parameters. The parameter $P_{\epsilon\epsilon\eta}^{[0]}$ is excluded because of its negligible impact.</figcaption></figure>

arises as a result of the combination of primary CMB and large-scale structure data. Measurements of the CMB temperature power spectrum constrain the parameter combination $A_s \exp(-2\tau)$. Massive neutrinos affect the growth of structure and therefore $M_ν$ becomes correlated with $A_s$ in large-scale structure measurements. When CMB and large-scale structure data are combined, the result is a strong degeneracy between $M_ν$ and $\tau$. In the right panel of figure 1, however, we see that this plays a much lesser role in the NLO case, as the neutrino mass constraint is supported primarily by other degeneracies.
Figure 2 shows the log derivatives of $P_{gg,s}$ with respect to $M_ν$ and the 9 nuisance parameters. We see that several of these parameters, in particular nonlinear bias parameters, change the power on small scales, and therefore lessen the information on $M_ν$ available from its characteristic suppression of the power spectrum on these scales. Using the NLO power spectrum but fixing the additional free parameters recovers almost the same constraints as in the linear power spectrum case (0.026 eV) for ΛCDM, while the constraint weakens to 0.044 eV when the new parameters are left free to vary. All of the nuisance parameters are somewhat correlated with $M_ν$ (except $P_{\epsilon\epsilon}^{(0)}$, which has little effect), and of particular importance is that the nuisance parameters are quite strongly correlated with each other (see figure 3).

In the linear case, imposing a prior on $b_1$ has the potential to significantly improve constraints on $M_ν$: if $b_1$ is fixed in each redshift bin, the $M_ν$ constraint falls from 0.026 eV to 0.021 eV. In the NLO case, fixing $b_1$ also improves the constraint, from 0.044 eV to 0.034 eV. Without any prior on $b_1$, the constraint in a typical redshift bin is about a factor of 7 weaker in the NLO case compared to in the linear case, because of the degeneracies with other nuisance parameters.

If $b_1$ is well constrained, adding priors on any of the other bias parameters provides further improvements. Fixing the two bias parameters associated with operators derived from the tidal field, $b_{K^2}$ and $b_{td}$, provides a further improvement in the constraints on $M_ν$ to 0.031 eV. $M_ν$ shows some correlation with both parameters. However, these two parameters are very strongly anti-correlated with each other, as is clear from figure 3. We discuss the reason for this strong degeneracy in section 4.4. $M_ν$ shows some correlation with both parameters. Because of the strong degeneracy between the parameters, fixing one without the other has very little impact on the neutrino mass constraint, as the other can compensate for its effect. In this way, these two bias parameters behave like a single effective parameter in this case. $P_{\epsilon}^{(0)}$ is also strongly correlated with $b_{K^2}$.

Fixing $b_1$ and the two tidal bias parameters improves the constraint to within 20% of that in the linear case. Fixing, for example, $b_1$, $C_{\nabla^2\delta}$ and $\beta_{\nabla^2\nu}$ achieves a similar constraint. Because the bias parameters are so correlated, adding priors on a few of them has a significant effect on the constraints on the others and therefore on the constraint on $M_ν$. Robust theoretical priors on the NLO bias parameters could potentially have a significant impact in closing the gap between these results and those derived for the linear power spectrum. Alternatively, measurements of the bispectrum could help reduce the uncertainty on $b_1$.

In [14], we showed that future CMB lensing data will add little to the combined constraints on the neutrino mass from large-scale structure measurements. However, as the NLO constraints are weakened, it is likely that the relative contribution from CMB lensing will be more substantial. The CMB lensing power spectrum is a function of the total matter power spectrum (including neutrino perturbations) and therefore measurements are unencumbered by the new galaxy power spectrum parameters introduced.

Figure 4 shows how the addition of forecast CMB lensing data from Simons Observatory could improve the NLO constraints. CMB lensing contributes primarily by breaking degeneracies between the nuisance parameters and the cosmological parameters (see also figure 3). It also reduces the cosmology dependence of the constraints. The galaxy-CMB lensing cross power spectrum is also included here, but makes very little contribution. Interestingly, once information from CMB lensing is added and the degeneracies between $M_ν$ and the bias parameters are under control, the degeneracy with $\tau$ once again becomes a significant factor in the overall constraints. CMB lensing contributes in particular by improving constraints on $A_s$, which in turn improves constraints on $b_1$, which also impacts the large-scale amplitude of the power spectrum.
<figure><img src="image_6.png" /></figure>

**Figure 4.** The NLO combined constraints with and without CMB lensing. The cross correlation between galaxies and CMB lensing is also included but has very little effect. The same results for the tree-level galaxy power spectrum and linear CMB lensing power spectrum are shown as the vertical tick marks. Comparing to figure 1, one sees that CMB lensing plays a much more substantial role when combined with the 1-loop power spectrum, by helping to break degeneracies between the cosmological parameters and the large number of nuisance parameters. Once these degeneracies are under control, we see in the right panel that the optical depth to reionisation $\tau$ becomes a significant driver of the neutrino mass constraint once again.

<figure><img src="image_7.png" /></figure>

**Figure 5.** As figure 1, but with the constraints derived using only the scale-dependent free-streaming signals in the galaxy, CMB lensing and galaxy-CMB lensing power spectra as a probe (see section 2.3). One can see that the constraints are quite independent of all the cosmological extensions considered in both the linear and NLO cases. These constraints are independent of $\tau$.

## 3.2 Free-streaming constraints

Figure 5 shows the forecasted constraints on the total neutrino mass using only the scale-dependence of the power spectrum as a constraining tool (see section 2.3). It can be seen that although the constraints are weakened somewhat by updating to the NLO power spectrum, they remain independent of changes in curvature and the dark energy equation of state. However, they degrade from around 0.08 eV to approximately 0.11 eV. This is a smaller relative change than in the combined case.

Interestingly, fixing all of the new nuisance parameters in the NLO case improves the constraint only to around 0.1 eV. From figure 6, we can see that the free-streaming signal is actually enhanced in the NLO case (figure 6 shows this for the matter power spectrum), and on investigation of the Fisher matrix, the information content on $M_\nu$ is actually increased. So the cause of this difference is not a change in the strength of signal in the power spectrum itself. The issue lies in the new degeneracy with $b_1$. Note that in the case of the linear power spectrum, extracting only the scale-dependence means treating the amplitude of the power spectrum as a free parameter, and $b_1$ is a simple multiplicative factor, so has no effect on the
<figure><img src="image_8.png" /><figcaption>**Figure 6.** The suppression in the matter power spectrum caused by the addition of a single massive neutrino of mass 0.06 eV relative to the case with no massive neutrinos, for both the linear and 1-loop/NLO cases. The relative suppression is enhanced in the NLO case. This enhanced effect propagates into both the galaxy and CMB lensing power spectra.</figcaption></figure>

final constraint on $M_ν$. However, in the NLO case $b_1$ now has a scale-dependent component and therefore some degeneracy with $M_ν$. Fixing $b_1$ plus all the new nuisance parameters gives a slightly better constraint in the NLO case than in the linear case, as expected.

One can now ask which of the new parameters contribute in particular to the weakening of the constraints on $M_ν$ seen in figure 5. The story is generally very similar to that in the combined case, but with significant changes for $b_1$ and $P_\epsilon^{\{0\}}$. This makes sense, as the other nuisance parameters do not change the large-scale power spectrum, so their effects when only the scale-dependence is considered are identical to those in the combined case.

Ultimately, considering that seven new nuisance parameters have been added to the calculation, all with scale-dependent effects, it could be considered surprising that a degradation in the overall constraints of only about 40% is observed. The cause of this is that the scale-dependence caused by neutrino free-streaming begins to take effect on much larger scales than those on which the scale-dependence of the bias parameters starts to be significant. Increasing the total neutrino mass shifts the free-streaming scale to higher $k$, and could therefore lead to greater degeneracy with the new parameters. This possibility is explored in section 4.2.

### 3.3 BAO constraints

Figure 7 shows the forecasted constraints on $M_ν$ when only constraints on the expansion rate $H(z)$ and angular diameter distance $D_A(z)$ provided by the baryon acoustic oscillations (BAO) in the power spectrum are used as a constraining tool. Importantly, CMB lensing is included, which has a massive impact on breaking cosmological degeneracies in this case. As is well-known, BAO are favoured for their robustness against non-linear effects, and that is reflected here. We see that the constraints degrade little when the 1-loop power spectrum is used instead of the linear power spectrum compared to the case of the full galaxy power spectrum. The constraint on $M_ν$ is degraded by approximately 20% in the ΛCDM case,
<figure><img src="image_9.png" /><figcaption>**Figure 7.** As figure 5, but with the constraints derived from distance information from the BAO signal only.</figcaption></figure>

and much less so for the more complex cosmologies, where the relative contribution of BAO measurements to the overall constraint is diminished. As in the linear case, the constraints remain cosmology-dependent, with a particular sensitivity to the dark energy equation of state (the curvature degeneracy is alleviated by the inclusion of CMB lensing).

We emphasise that in the calculation for figure 7, we include BAO damping using an IR-resummation procedure [23, 40–43], as described below, although BAO damping is not included in our main model of the NLO galaxy power spectrum used in other sections. We add this here to provide a fair comparison of the constraints in the linear and NLO BAO cases. Overall we find it has a minor effect on the $M_\nu$ constraints and we can therefore justify neglecting it in the other calculations. Without IR resummation, the constraints in the NLO case in figure 7 improve by less than 10% in the $\Lambda$CDM case and by an even lesser extent for all of the other cosmologies. For our ‘combined’ calculations, BAO make up only one component of the total constraining power so excluding IR resummation in that case should be reasonable. Finally, we note that, when only the BAO feature is used, the damping of the BAO signal can also be substantially mitigated through BAO reconstruction [44, 45], so the NLO constraints in figure 7 could be viewed as conservative.

We follow [46] to implement anisotropic BAO damping. When calculating the derivatives of the BAO component of the NLO galaxy power spectrum, we multiply the wiggles by a factor $\exp[-k^2\Sigma_{\text{tot}}(\mu)^2]$ where

$$ \Sigma_{\text{tot}}(\mu)^2 = (1 + f\mu^2(2 + f))\Sigma^2 + f^2\mu^2(\mu^2 - 1)\delta\Sigma^2 \quad (3.1) $$

$$ \Sigma^2 = \frac{1}{(2\pi)^3} \frac{4\pi}{3} \int_0^{k_S} dq P_{nw}(q) \left[ 1 - j_0\left(\frac{q}{k_{\text{osc}}}\right) + 2j_2\left(\frac{q}{k_{\text{osc}}}\right) \right] \quad (3.2) $$

$$ \delta\Sigma^2 = \frac{1}{(2\pi)^3} 4\pi \int_0^{k_S} dq P_{nw}(q) j_2\left(\frac{q}{k_{\text{osc}}}\right) \quad (3.3) $$

$P_{nw}(k)$ is the baryon-CDM power spectrum without BAO (‘no-wiggle’), $k_{\text{osc}}$ is the BAO scale $1/(110\ h\ \text{Mpc}^{-1})$, $k_S$ is set to $0.2\ h\ \text{Mpc}^{-1}$, and $j_n$ are spherical Bessel functions of order $n$.

# 4 Discussion

## 4.1 Selection effects

The observed galaxy power spectrum will always have a dependence on galaxy properties that alter the probability of a given galaxy being observed. Biases that arise in the observed
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Contours showing the relationship between $b_\eta$ and a subset of the other parameters in our calculations, for the central redshift bin of the combined galaxy clustering Euclid forecast (without CMB lensing). Apart from $M_\nu$, these parameters include those with the most significant (anti-)correlation with $b_\eta$. Despite some degree of correlation between $b_\eta$, $A_s$ and $\tau$, the effect of varying $b_\eta$ on the constraint on $M_\nu$ is quite mild.</figcaption></figure>

power spectrum as a result of the specific subset of galaxies observed are called selection effects. Many selection effects are included by default in the bias expansion used here, but we have neglected selection effects that depend on the line-of-sight in this analysis. We will now justify this.

One particularly crucial selection effect arises from the dependence of the apparent brightness of a galaxy on its orientation with respect to us, as brighter galaxies are more likely to be detected. This selection effect is degenerate with the Kaiser effect [47, 48], and alters the value of $b_\eta$, which we fix to -1 when neglecting selection effects. An analogous effect arises when galaxies are selected on resonant line flux, such as Lyman-$\alpha$ [49–51]. We can briefly examine the impact of freeing $b_\eta$ to get an impression of its impact.

Without any prior on $b_\eta$, in the combined case for $\Lambda$CDM, the NLO constraint on $M_\nu$ degrades from 0.044 eV to 0.051 eV. Imposing a $1\sigma$ prior on $b_\eta$ of 0.05 improves the constraint to 0.046 eV. Freeing $b_\eta$ does not have any significant effect on the BAO-only or free-streaming constraints. Figure 8 is a contour plot for $b_\eta$ and a small sample of other parameters in our calculations in our combined galaxy clustering forecasts (with marginalisation over the other parameters not shown). There is some clear degeneracy between $b_\eta$, $A_s$ and $\tau$, but the effect of $b_\eta$ on $M_\nu$ seems quite mild despite $M_\nu$ being correlated with the latter two. Of the nuisance parameters, $b_\eta$ shows the most clear degeneracy with $b_2$ and $P_\epsilon^{(0)}$.

We finally note that Euclid will also image the galaxies that they observe [1]. This should allow for the intentional selection or reweighting of a galaxy sample in such a way as to reduce the line-of-sight dependent selection effects.

## 4.2 Varying the fiducial neutrino mass

In section 3.2 we found that the scale-dependent suppression in the power spectrum caused by neutrinos is not particularly degenerate with any of the new parameters. The suppression of the power spectrum by neutrino free-streaming starts at larger scales than those at which the scale-dependent effects of the new bias parameters become significant. However, the free-streaming scale $k_{\text{FS}}$ increases with the neutrino mass. It is therefore worth investigating whether the constraints for a higher fiducial neutrino mass would be further weakened.

We ran the combined results again for an increased $M_\nu$ of 0.27 eV assuming 3 degenerate mass states, and the constraints for all cosmologies were within a few percent of those obtained for the minimal neutrino mass. We conclude that any weakening of the constraints due to greater degeneracy with the non-linear bias parameters is compensated for by the increased amplitude of the suppression that comes with higher mass neutrinos.
## 4.3 Neutrino-induced bias

Something we have neglected in our treatment in this article is the effects of massive neutrinos themselves on the bias parameters. As we have seen, the free-streaming of massive neutrinos introduces a scale-dependence in the growth of perturbations, not just for neutrinos, but also cold dark matter and baryons. This should then result in some scale-dependence in the bias parameters. This was studied in detail by [52] for $b_1$, who modified the spherical collapse and peak-background split derivations of the large-scale bias to account for this. Including this scale-dependent bias should work to counter the suppression caused by massive neutrinos, in principle weakening the signal. The same effect should also cause a scale-dependence in the higher-order bias parameters.

This is a complicated topic and a full model of the next-to-leading-order power spectrum that accounts for this effect has yet to be developed. In addition, there is some disagreement on the precise size of the effect in N-body simulations [28, 53–56]; see [57] for an estimate of the relevance of this effect for future surveys. We therefore leave this for future work.

## 4.4 Degeneracy between bias parameters

We highlighted in section 3.2 that there is a significant anti-correlation between the two parameters associated with the tidal field, $b_{K^2}$ and $b_{td}$. An examination of the contributions of these terms to the power spectrum reveals why.

Both $b_{td}$ and $b_{K^2}$ contribute to the $P_{gg}^{1-3}$ term. In [20], the authors show that the contributions of the associated operators to $P_{gg}^{1-3}$ are linearly proportional. Their contributions are identical apart from a constant factor of 5/2, making the effects of varying the two bias parameters fully degenerate.

The $b_{K^2}$ parameter also contributes to the $P_{gg}^{2-2}$ term, but $b_{td}$ does not, preventing the two parameters from being perfectly anti-correlated. However, on inspection, we see that the relative change in $P_{gg}^{2-2}$ when varying $b_{K^2}$ is significantly smaller than that in $P_{gg}^{1-3}$, so changing $P_{gg}^{1-3}$ is the dominant effect of $b_{K^2}$. The parameters $b_{K^2}$ and $b_{td}$ therefore remain quite strongly anti-correlated, even after marginalising over other parameters.

## 4.5 Comparisons with previous work

A recent article [19] performed an MCMC forecast for Euclid to 1-loop order in the power spectrum, as well as the tree-level bispectrum. Combining the current Planck likelihood with their Euclid forecast, they quote a constraint of 0.017 eV (excluding the bispectrum). Unlike us, they do not vary $N_{\text{eff}}$. When we fix $N_{\text{eff}}$, we obtain a constraint of 0.034 eV. There are two obvious differences between [19] and our work: 1) they use a much higher scale cut of $k_{\text{max}} = 0.5h \, \text{Mpc}^{-1}$ but include a theoretical error which is a smooth increasing function of $k$; and 2) they use only six instead of our nine nuisance parameters. In their appendices, the authors show what happens if they impose a sharp cut-off $k_{\text{max}}$ of $0.15 \, h \, \text{Mpc}^{-1}$, and their constraints degrade significantly, although they do not include CMB information in this case and do include the bispectrum, which we do not consider here, so it is not possible to perform a direct comparison. We also note that fixing $\tau$ in our own case takes the constraint from 0.034 eV to 0.023 eV, and a different constraint on $\tau$ may also contribute to the difference in results.

We note that in our work using the linear power spectrum only, our results using Fisher forecasts with sharp cut-offs in $k$ produced results that were remarkably consistent with MCMC forecasts with the controlled theoretical error approach (see the comparison to [18] in [14]).
<figure><img src="image_11.png" /><figcaption>**Figure 9.** The forecasted 1σ constraints on $M_ν$ obtained using the linear (teal line) and 1-loop (purple line) galaxy power spectra for various $k_{max}$. The corresponding parameter biases resulting from not going to the next highest order in both cases are shown as points. The biases on $M_ν$ in the linear case are actually negative. In the 1-loop case, the biases are relatively small and fluctuate between positive and negative depending on the cut-off scale. This explains why the bias does not monotonically increase with $k_{max}$. At $k_{max} = 0.2 h \text{ Mpc}^{-1}$, the value we use in our calculations, the 1-loop bias is approximately 20% of the constraint.</figcaption></figure>

The authors of [19] also examined the specific contribution of BAO to their final constraints, and found it to be modest. This emphasises our point in our previous work and here that neutrino mass constraints from BAO alone waste a great deal of usable information.

## 4.6 Biases on $M_ν$ from the power spectrum model

Using an inappropriate model for an observable (for example, using the linear power spectrum) will result in systematic biases in the parameters eventually constrained. Assume an observable $O(x)$ with observed values $\hat{O}$, where $x$ denotes the elements of the data vector, such as wavenumber and redshift bins, and cosmological parameters $\theta$ with true values $\bar{\theta}_i$. It is straightforward to extend the Fisher matrix formalism to estimate the magnitude of the biases produced by an incomplete model:

$$b[\theta_i] = \langle \theta_{\text{sys},i} \rangle - \langle \bar{\theta}_i \rangle = \sum_j (F^{-1})_{ij} \sum_x \frac{O^{\text{sys}}(x)}{\text{Var}[\hat{O}(x)]} \frac{\partial O^{\text{theory}}(x)}{\partial \theta_j}, \quad (4.1)$$

where $\theta_{\text{sys},i}$ are the incorrectly inferred values of the cosmological parameters, $F$ is the Fisher matrix obtained using the incomplete model and $\partial O^{\text{theory}}/\partial \theta_j$ the corresponding derivatives. $O^{\text{sys}}$ is the systematic error in the model, i.e. the difference between the accurate and systematically incorrect fiducial models.
For our example, we combine power spectra from multiple redshift bins to obtain the final result. Each redshift bin can be seen as an independent experiment. We also want to account for a CMB prior $F_{prior,ij}$ that is not affected by this systematic error. The bias for a given parameter as a result of using the linear galaxy power spectrum as opposed to the NLO power spectrum can then be calculated as

$$b[\theta_i] = \langle\theta_{\text{sys},i}\rangle - \langle\bar{\theta}_i\rangle = \sum_j \left(F_{\text{prior}} + \sum_z F_z\right)_{ij}^{-1} \sum_z \sum_k \frac{P_{gg}^{\text{sys}}(z, k)}{\text{Var}[P_{gg}(z, k)]} \frac{\partial P_{gg}(z, k)}{\partial \theta_j}. \quad (4.2)$$

Here, $P_{gg}^{\text{sys}}(z, k) = P_{gg,1-\text{loop}}(z, k) - P_{gg,L}(z, k)$. Figure 9 shows the biases in the constraints on $M_ν$ inferred from using the linear power spectrum, as a function of $k_{\text{max}}$. It is clear that the bias (shown as teal stars) significantly exceeds the linear power spectrum $1-σ$ constraint (shown as a teal line) even at $k_{\text{max}} = 0.1$ h/Mpc. These biases are also negative, as would be expected: the 1-loop power spectrum has enhanced clustering on small scales. Using the linear power spectrum to infer the neutrino mass would therefore make the power spectrum seem less suppressed on small scales, resulting in an inferred value of $M_ν$ that is too small, resulting in the bias in equation (4.2) being negative.

Very conveniently, we require only the fiducial corrected power spectrum to estimate the bias from using the incorrect model in equation (4.2). This means we can also make an estimate of how much our results would be affected by not moving to the next order, 2-loop power spectrum. We make use of the 2-loop matter power spectrum calculated for [26] and model the 2-loop contribution to the galaxy power spectrum as approximately $P_{gg,2-\text{loop}} = (b_1 + f\mu^2)^2 P_{m,2-\text{loop}}(k)$. It is clear that the parameter biases are much less significant and the 1-loop power spectrum is much more robust up to higher $k_{\text{max}}$. The bias is about 20% of the error at $k_{\text{max}} = 0.2$ h Mpc$^{-1}$.

The above result for the bias on the 1-loop power spectrum is intended to be approximate. The two-loop matter power spectrum used is for a slightly different cosmology and is only calculated at a single redshift and rescaled (so does not include the redshift-dependent effects of neutrinos). More importantly, it also neglects higher-order bias terms arising in the two-loop galaxy power spectrum, whose coefficients could be relatively less suppressed at higher redshifts.

## 4.7 Alternative galaxy surveys

One might ask how the conclusions derived for the Euclid-like survey in this work could change with the survey parameters. The Roman Telescope (formerly WFIRST) survey [3] provides a natural complement to the Euclid survey for examination, having a considerably smaller area (2000 deg$^2$ in our model) but probing deeper redshifts (up to $z = 2.8$).

Our model for the Roman telescope is given in table 3 of appendix A. We have updated our model to match that provided as Design Reference Mission (DRM) 1 in [3], so the results should not be directly compared with those for ‘WFIRST’ provided in [13, 14], particularly as the area of the survey has been updated from 2000 to 3400 deg$^2$. The total number of galaxies here is approximately 16.5 million but with a greater redshift depth than Euclid (up to $z = 2.7$). Note that the survey model is also different to that analysed in [26], which focused exclusively on the Hα-emitting galaxies.

At $k_{\text{max}} = 0.2$ h Mpc$^{-1}$, while Euclid gave $\sigma(M_ν) = 0.042$ eV with a bias of 0.008 eV relative to the 2-loop model, the Roman Telescope gives $\sigma(M_ν) = 0.078$ eV with a bias of 0.0028 eV. The approximate factor of two decrease in the constraint is accompanied by an
improvement in the bias, which is probably a result of the increased redshift range, as the relative importance of nonlinear corrections to the power spectrum shrinks with increasing $z$. Increasing $k_{\text{max}}$ to $0.35 \, h \, \text{Mpc}^{-1}$ is required to reduce the error on $M_\nu$ to that of Euclid, achieving a forecasted error for the Roman Telescope of $0.042 \, \text{eV}$ with a bias of $0.0032$.

# 5 Conclusions

The aim of this paper has been to provide an indication of how forecasted neutrino mass constraints from galaxy clustering and CMB surveys are altered on moving to next-to-leading-order power spectra. The results are relatively encouraging considering that the next-to-leading-order galaxy power spectrum implemented here introduces 7 new unconstrained nuisance parameters. When performing Fisher forecasts for galaxy clustering for a Euclid-like survey and CMB lensing from an experiment like Simons Observatory, we find:

* When the full redshift-space galaxy power spectrum is used, the constraints on $M_\nu$ in the $\Lambda$CDM case degrade from $0.026 \, \text{eV}$ in the linear case to $0.044 \, \text{eV}$ in the NLO case. This degrades to $0.12 \, \text{eV}$ if both curvature and the dark energy equation of state are also allowed to vary. The degeneracy between $M_\nu$ and $\tau$ becomes much less significant in the NLO case, as the error on $M_\nu$ is primarily inflated by uncertainties on the new nuisance parameters. Adding CMB lensing information tightens the NLO constraint to $0.034 \, \text{eV}$ by helping reduce the level of degeneracy with the nuisance parameters, and the connection between the errors on $\tau$ and $M_\nu$ then becomes significant once again.

* When only distance information from BAO are used and combined with CMB lensing, the constraints are approximately the same in the linear and NLO cases: $0.048 \, \text{eV}$ in the $\Lambda$CDM case and $0.82 \, \text{eV}$ when curvature and the dark energy equation of state are left free.

* When only measurements of the direct free-streaming effect in the power spectra are used (the scale-dependence in the matter power spectrum and the structure growth rate), the constraints on $M_\nu$ are independent of $\Omega_k$ and $w$ and approximately $0.11 \, \text{eV}$ in all cases (to be compared to $0.08 \, \text{eV}$ for the linear power spectra).

CMB lensing plays a much more important role here than in [14], where we examined the linear case, but the effect of adding the cross-correlation between galaxy clustering and CMB lensing remains largely negligible in the constraints.

In section 4.6, we attempted to quantify the biases introduced in predictions of the neutrino mass by assuming an incomplete model for the power spectrum. We showed that the bias resulting from assuming the linear galaxy power spectrum already almost exceeds the predicted constraint on $M_\nu$ at $k_{\text{max}}$ of $0.1 \, h \, \text{Mpc}^{-1}$. An error estimate for the 1-loop power spectrum, furnished by the 2-loop matter power spectrum, shows that the former remains robust to much higher $k_{\text{max}}$, with a bias of approximately 20% of the $M_\nu$ constraint at $k_{\text{max}}$ of $0.2 \, h \, \text{Mpc}^{-1}$.

In terms of future work, we note that combined analysis with the bispectrum, which we do not consider here, could help break degeneracies [58] and better constrain the bias parameters. The effects of massive neutrinos themselves on the bias parameters should also be analysed (see section 4.3).

Finally, we emphasise that constraints on the mass of neutrinos from cosmology are generally highly sensitive to the assumed cosmological model, and should be quoted with the appropriate caveats.
<table><caption>**Table 2.** Survey parameters assumed for Euclid.</caption><thead><tr><th>z</th><th>V (h<sup>-3</sup> Gpc<sup>3</sup>)</th><th>10<sup>3</sup> × n<sub>g</sub> (Mpc<sup>3</sup> h<sup>-3</sup>)</th><th>b<sub>1</sub></th><th>R<sub>L</sub> (Mpc h<sup>-1</sup>)</th></tr></thead><tbody><tr><td>0.65</td><td>2.59</td><td>0.637</td><td>1.07</td><td>2.41</td></tr><tr><td>0.75</td><td>3.07</td><td>1.441</td><td>1.12</td><td>1.88</td></tr><tr><td>0.85</td><td>3.52</td><td>1.619</td><td>1.17</td><td>1.79</td></tr><tr><td>0.95</td><td>3.93</td><td>1.489</td><td>1.23</td><td>1.8</td></tr><tr><td>1.05</td><td>4.29</td><td>1.32</td><td>1.28</td><td>1.83</td></tr><tr><td>1.15</td><td>4.62</td><td>1.144</td><td>1.34</td><td>1.84</td></tr><tr><td>1.25</td><td>4.9</td><td>0.995</td><td>1.39</td><td>1.86</td></tr><tr><td>1.35</td><td>5.14</td><td>0.832</td><td>1.45</td><td>1.9</td></tr><tr><td>1.45</td><td>5.35</td><td>0.659</td><td>1.51</td><td>1.94</td></tr><tr><td>1.55</td><td>5.52</td><td>0.503</td><td>1.57</td><td>2.0</td></tr><tr><td>1.65</td><td>5.66</td><td>0.364</td><td>1.62</td><td>2.08</td></tr><tr><td>1.75</td><td>5.78</td><td>0.253</td><td>1.68</td><td>2.18</td></tr><tr><td>1.85</td><td>5.88</td><td>0.166</td><td>1.74</td><td>2.28</td></tr><tr><td>1.95</td><td>5.95</td><td>0.101</td><td>1.8</td><td>2.43</td></tr><tr><td>2.05</td><td>6.01</td><td>0.037</td><td>1.86</td><td>2.72</td></tr></tbody></table>

## Acknowledgments

We would like to thank Eiichiro Komatsu for providing feedback on the first draft of this paper, which appeared as a chapter in AB's doctoral thesis. We would also like to thank Nishant Agrawal, for helping so generously with code and results comparisons. Finally we would like to thank Donghui Jeong and Vincent Desjacques for useful discussions and assistance. FS acknowledges support from the Starting Grant (ERC-2015-STG 678652) "GrInflaGal" of the European Research Council.

## A Survey parameters

Tables 2 and 3 give the survey parameters for Euclid and the Roman Telescope, respectively. The final columns give the estimated halo Lagrangian radii used to determine the fiducial values of $C_{\nabla^{2}\delta}$ and $\beta_{\nabla^{2}v}$.
<table><caption>**Table 3.** Survey parameters assumed for the Roman Telescope.</caption><thead><tr><th>z</th><th>V (h<sup>-3</sup> Gpc<sup>3</sup>)</th><th>10<sup>3</sup> × n<sub>g</sub> (Mpc<sup>3</sup> h<sup>-3</sup>)</th><th>b<sub>1</sub></th><th>R<sub>L</sub> (Mpc h<sup>-1</sup>)</th></tr></thead><tbody><tr><td>1.35</td><td>1.08</td><td>1.253</td><td>1.44</td><td>1.72</td></tr><tr><td>1.45</td><td>1.13</td><td>1.24</td><td>1.48</td><td>1.68</td></tr><tr><td>1.55</td><td>1.17</td><td>1.18</td><td>1.52</td><td>1.65</td></tr><tr><td>1.65</td><td>1.21</td><td>1.091</td><td>1.56</td><td>1.65</td></tr><tr><td>1.75</td><td>1.24</td><td>1.007</td><td>1.6</td><td>1.63</td></tr><tr><td>1.85</td><td>1.26</td><td>0.944</td><td>1.64</td><td>1.62</td></tr><tr><td>1.95</td><td>1.28</td><td>0.961</td><td>1.68</td><td>1.58</td></tr><tr><td>2.05</td><td>1.3</td><td>0.988</td><td>1.72</td><td>1.52</td></tr><tr><td>2.15</td><td>1.31</td><td>1.021</td><td>1.76</td><td>1.47</td></tr><tr><td>2.25</td><td>1.32</td><td>1.001</td><td>1.8</td><td>1.44</td></tr><tr><td>2.35</td><td>1.33</td><td>0.858</td><td>1.84</td><td>1.44</td></tr><tr><td>2.45</td><td>1.33</td><td>0.718</td><td>1.88</td><td>1.46</td></tr><tr><td>2.55</td><td>1.33</td><td>0.552</td><td>1.92</td><td>1.5</td></tr><tr><td>2.65</td><td>1.33</td><td>0.416</td><td>1.96</td><td>1.53</td></tr></tbody></table>

## References

[1] L. Amendola et al., *Cosmology and fundamental physics with the Euclid satellite*, Living Rev. Rel. 21 (2018) 2 [arXiv:1606.00180] [INSPIRE].

[2] EUCLID collaboration, *Euclid definition study report*, arXiv:1110.3193 [INSPIRE].

[3] J. Green et al., *Wide-Field InfraRed Survey Telescope (WFIRST) final report*, arXiv:1208.4012 [INSPIRE].

[4] N. Tamura et al., *Prime Focus Spectrograph (PFS) for the Subaru telescope: overview, recent progress, and future perspectives*, Proc. SPIE Int. Soc. Opt. Eng. 9908 (2016) 99081M [arXiv:1608.01075] [INSPIRE].

[5] DESI collaboration, *The DESI experiment part I: science, targeting, and survey design*, arXiv:1611.00036 [INSPIRE].

[6] SIMONS OBSERVATORY collaboration, *The Simons Observatory: science goals and forecasts*, JCAP 02 (2019) 056 [arXiv:1808.07445] [INSPIRE].

[7] CMB-S4 collaboration, *CMB-S4 science book, first edition*, arXiv:1610.02743 [INSPIRE].

[8] *Valencia-Globalfit webpage*, http://globalfit.astroparticles.es/, (2020).

[9] P.F. de Salas et al., 2020 global reassessment of the neutrino oscillation picture, JHEP 02 (2021) 071 [arXiv:2006.11237] [INSPIRE].

[10] PLANCK collaboration, *Planck 2018 results. VI. Cosmological parameters*, Astron. Astrophys. 641 (2020) A6 [arXiv:1807.06209] [INSPIRE].

[11] R. Allison, P. Caucal, E. Calabrese, J. Dunkley and T. Louis, *Towards a cosmological neutrino mass detection*, Phys. Rev. D 92 (2015) 123535 [arXiv:1509.07471] [INSPIRE].

[12] M. Archidiacono, T. Brinckmann, J. Lesgourgues and V. Poulin, *Physical effects involved in the measurements of neutrino masses with future cosmological data*, JCAP 02 (2017) 052 [arXiv:1610.09852] [INSPIRE].
[13] A. Boyle and E. Komatsu, Deconstructing the neutrino mass constraint from galaxy redshift surveys, JCAP 03 (2018) 035 [arXiv:1712.01857] [INSPIRE].

[14] A. Boyle, Understanding the neutrino mass constraints achievable by combining CMB lensing and spectroscopic galaxy surveys, JCAP 04 (2019) 038 [Erratum ibid. 10 (2019) E01] [arXiv:1811.07636] [INSPIRE].

[15] S. Roy Choudhury and S. Hannestad, Updated results on neutrino mass and mass hierarchy from cosmology with Planck 2018 likelihoods, JCAP 07 (2020) 037 [arXiv:1907.12598] [INSPIRE].

[16] E. Chung, S. Foreman and A. van Engelen, Baryonic effects on CMB lensing and neutrino mass constraints, Phys. Rev. D 101 (2020) 063534 [Erratum ibid. 102 (2020) 109903] [arXiv:1910.09565] [INSPIRE].

[17] T. Baldauf, M. Mirbabayi, M. Simonović and M. Zaldarriaga, LSS constraints with controlled theoretical uncertainties, arXiv:1602.00674 [INSPIRE].

[18] T. Brinckmann, D.C. Hooper, M. Archidiacono, J. Lesgourgues and T. Sprenger, The promising future of a robust cosmological neutrino mass measurement, JCAP 01 (2019) 059 [arXiv:1808.05955] [INSPIRE].

[19] A. Chudaykin and M.M. Ivanov, Measuring neutrino masses with large-scale structure: Euclid forecast with controlled theoretical error, JCAP 11 (2019) 034 [arXiv:1907.06666] [INSPIRE].

[20] V. Desjacques, D. Jeong and F. Schmidt, The galaxy power spectrum and bispectrum in redshift space, JCAP 12 (2018) 035 [arXiv:1806.04015] [INSPIRE].

[21] R. Scoccimarro, Redshift-space distortions, pairwise velocities and nonlinearities, Phys. Rev. D 70 (2004) 083007 [astro-ph/0407214] [INSPIRE].

[22] U. Seljak and P. McDonald, Distribution function approach to redshift space distortions, JCAP 11 (2011) 039 [arXiv:1109.1888] [INSPIRE].

[23] A. Perko, L. Senatore, E. Jennings and R.H. Wechsler, Biased tracers in redshift space in the EFT of large-scale structure, arXiv:1610.09321 [INSPIRE].

[24] L. Fonseca de la Bella, D. Regan, D. Seery and D. Parkinson, Impact of bias and redshift-space modelling for the halo power spectrum: testing the effective field theory of large-scale structure, JCAP 07 (2020) 011 [arXiv:1805.12394] [INSPIRE].

[25] M. Mirbabayi, F. Schmidt and M. Zaldarriaga, Biased tracers and time evolution, JCAP 07 (2015) 030 [arXiv:1412.5169] [INSPIRE].

[26] N. Agarwal, V. Desjacques, D. Jeong and F. Schmidt, Information content in the redshift-space galaxy power spectrum and bispectrum, JCAP 03 (2021) 021 [arXiv:2007.04340] [INSPIRE].

[27] D. Blas, J. Lesgourgues and T. Tram, The Cosmic Linear Anisotropy Solving System (CLASS) II: approximation schemes, JCAP 07 (2011) 034 [arXiv:1104.2933] [INSPIRE].

[28] E. Castorina, C. Carbone, J. Bel, E. Sefusatti and K. Dolag, DEMNUni: the clustering of large-scale structures in the presence of massive neutrinos, JCAP 07 (2015) 043 [arXiv:1505.07148] [INSPIRE].

[29] O.F. Hernández, Neutrino masses, scale-dependent growth, and redshift-space distortions, JCAP 06 (2017) 018 [arXiv:1608.08298] [INSPIRE].

[30] S. Saito, M. Takada and A. Taruya, Impact of massive neutrinos on nonlinear matter power spectrum, Phys. Rev. Lett. 100 (2008) 191301 [arXiv:0801.0607] [INSPIRE].

[31] S. Saito, M. Takada and A. Taruya, Nonlinear power spectrum in the presence of massive neutrinos: perturbation theory approach, galaxy bias and parameter forecasts, Phys. Rev. D 80 (2009) 083528 [arXiv:0907.2922] [INSPIRE].
[32] D. Jeong and E. Komatsu, *Perturbation theory reloaded: analytical calculation of non-linearity in baryonic oscillations in the real space matter power spectrum*, Astrophys. J. 651 (2006) 619 [astro-ph/0604075] [INSPIRE].

[33] T. Lazeyras and F. Schmidt, *A robust measurement of the first higher-derivative bias of dark matter halos*, JCAP 11 (2019) 041 [arXiv:1904.11294] [INSPIRE].

[34] J.L. Tinker et al., *Toward a halo mass function for precision cosmology: the limits of universality*, Astrophys. J. 688 (2008) 709 [arXiv:0803.2706] [INSPIRE].

[35] T. Lazeyras, C. Wagner, T. Baldauf and F. Schmidt, *Precision measurement of the local bias of dark matter halos*, JCAP 02 (2016) 018 [arXiv:1511.01096] [INSPIRE].

[36] J. Peloton, M. Schmittfull, A. Lewis, J. Carron and O. Zahn, *Full covariance of CMB and lensing reconstruction power spectra*, Phys. Rev. D 95 (2017) 043508 [arXiv:1611.01446] [INSPIRE].

[37] J. Rocher, K. Benabed and F. Bouchet, *Probing inflation with CMB polarization: weak lensing effect on the covariance of CMB spectra*, JCAP 05 (2007) 013 [astro-ph/0612662] [INSPIRE].

[38] S. Mishra-Sharma, D. Alonso and J. Dunkley, *Neutrino masses and beyond-ΛCDM cosmology with LSST and future CMB experiments*, Phys. Rev. D 97 (2018) 123544 [arXiv:1803.07561] [INSPIRE].

[39] B. Yu, R.Z. Knight, B.D. Sherwin, S. Ferraro, L. Knox and M. Schmittfull, *Towards neutrino mass from cosmology without optical depth information*, arXiv:1809.02120 [INSPIRE].

[40] L. Senatore and M. Zaldarriaga, *The IR-resummed effective field theory of large scale structures*, JCAP 02 (2015) 013 [arXiv:1404.5954] [INSPIRE].

[41] T. Baldauf, M. Mirbabayi, M. Simonović and M. Zaldarriaga, *Equivalence principle and the baryon acoustic peak*, Phys. Rev. D 92 (2015) 043514 [arXiv:1504.04366] [INSPIRE].

[42] D. Blas, M. Garny, M.M. Ivanov and S. Sibiryakov, *Time-sliced perturbation theory II: baryon acoustic oscillations and infrared resummation*, JCAP 07 (2016) 028 [arXiv:1605.02149] [INSPIRE].

[43] L. Senatore and G. Trevisan, *On the IR-resummation in the EFTofLSS*, JCAP 05 (2018) 019 [arXiv:1710.02178] [INSPIRE].

[44] D.J. Eisenstein, H.-J. Seo, E. Sirko and D. Spergel, *Improving cosmological distance measurements by reconstruction of the baryon acoustic peak*, Astrophys. J. 664 (2007) 675 [astro-ph/0604362] [INSPIRE].

[45] N. Padmanabhan, M. White and J.D. Cohn, *Reconstructing baryon oscillations: a Lagrangian theory perspective*, Phys. Rev. D 79 (2009) 063523 [arXiv:0812.2905] [INSPIRE].

[46] M.M. Ivanov and S. Sibiryakov, *Infrared resummation for biased tracers in redshift space*, JCAP 07 (2018) 053 [arXiv:1804.05080] [INSPIRE].

[47] C.M. Hirata, *Tidal alignments as a contaminant of redshift space distortions*, Mon. Not. Roy. Astron. Soc. 399 (2009) 1074.

[48] E. Krause and C. Hirata, *Tidal alignments as a contaminant of the galaxy bispectrum*, Mon. Not. Roy. Astron. Soc. 410 (2011) 2730 [arXiv:1004.3611] [INSPIRE].

[49] Z. Zheng, R. Cen, H. Trac and J. Miralda-Escudé, *Radiative transfer modeling of Lyman-α emitters. II. New effects in galaxy clustering*, Astrophys. J. 726 (2010) 38 [arXiv:1003.4990] [INSPIRE].

[50] S. Wyithe and M. Dijkstra, *Non-gravitational contributions to the clustering of Ly-α selected galaxies: implications for cosmological surveys*, Mon. Not. Roy. Astron. Soc. 415 (2011) 3929 [arXiv:1104.0712] [INSPIRE].
[51] B. Greig, E. Komatsu and J.S.B. Wyithe, *Cosmology from clustering of Lyman-α galaxies: breaking non-gravitational Lyman-α radiative transfer degeneracies using the bispectrum*, Mon. Not. Roy. Astron. Soc. **431** (2013) 1777 [arXiv:1212.0977] [INSPIRE].

[52] M. LoVerde, *Halo bias in mixed dark matter cosmologies*, Phys. Rev. D **90** (2014) 083530 [arXiv:1405.4855] [INSPIRE].

[53] E. Castorina, E. Sefusatti, R.K. Sheth, F. Villaescusa-Navarro and M. Viel, *Cosmology with massive neutrinos II: on the universality of the halo mass function and bias*, JCAP **02** (2014) 049 [arXiv:1311.1212] [INSPIRE].

[54] F. Villaescusa-Navarro et al., *The imprint of neutrinos on clustering in redshift-space*, Astrophys. J. **861** (2018) 53 [arXiv:1708.01154] [INSPIRE].

[55] C.-T. Chiang, W. Hu, Y. Li and M. Loverde, *Scale-dependent bias and bispectrum in neutrino separate universe simulations*, Phys. Rev. D **97** (2018) 123526 [arXiv:1710.01310] [INSPIRE].

[56] C.-T. Chiang, M. LoVerde and F. Villaescusa-Navarro, *First detection of scale-dependent linear halo bias in N-body simulations with massive neutrinos*, Phys. Rev. Lett. **122** (2019) 041302 [arXiv:1811.12412] [INSPIRE].

[57] A. Raccanelli, L. Verde and F. Villaescusa-Navarro, *Biases from neutrino bias: to worry or not to worry?*, Mon. Not. Roy. Astron. Soc. **483** (2019) 734 [arXiv:1704.07837] [INSPIRE].

[58] C. Hahn, F. Villaescusa-Navarro, E. Castorina and R. Scoccimarro, *Constraining $M_ν$ with the bispectrum*. Part I. *Breaking parameter degeneracies*, JCAP **03** (2020) 040 [arXiv:1909.11107] [INSPIRE].