© The Authors 2024

1. Introduction

Observations of neutrino oscillations (see Gonzalez-Garcia et al. 2016, for a review of the main experimental data) have made it clear that the standard model of particle physics should be extended in order to accommodate massive neutrinos. However, the absolute mass scale of these particles cannot be probed through flavor oscillation experiments, because these are only sensitive to the squared mass differences of the mass eigenstates¹. Complementary efforts are thus needed to obtain additional information, and examples are the kinematic analysis of the β decay (Mertens 2016) and the search for the neutrinoless double-β decay (Dolinski et al. 2019), which probe effective masses involving the mixing angles, Major charge conjugation parity-violating phases, and the physical masses.

In the last few years, cosmology has been recognized as a powerful and independent tool for constraining the neutrino global mass scale. Indeed, the cosmic microwave background (CMB) anisotropy power spectra, galaxy clustering, and weak cosmological lensing are all cosmological observables that depend (at least to first order) on the sum of neutrino masses; hence their relevance for the characterization of these particles and the minimal extensions of the Λ cold dark matter (ΛCDM) cosmological model. In particular, CMB anisotropies reign supreme among all cosmological probes in terms of their power to constrain neutrino mass. For instance, under a ΛCDM model with massive neutrinos, the Planck analysis of the CMB temperature anisotropy power spectrum alone yielded a 95% upper limit of ∑m_ν < 0.54 eV, which was improved to ∑m_ν < 0.24 eV when the polarization signal was considered (Planck Collaboration VI 2020). This upper bound can be further decreased to ∑m_ν < 0.12 eV when baryon acoustic oscillations are additionally taken into account (Planck Collaboration VI 2020).

However, it has been argued that these constraints are potentially subject to the effect of the so-called lensing anomaly (Calabrese & Slosar 2008) present in the CMB power spectra. Indeed, the analysis of Planck data yielded a value for the amplitude of CMB lensing that was higher than expected (Planck Collaboration VI 2020; Valentino et al. 2020), which could bias results on neutrino masses toward smaller values (Capozzi et al. 2021). This motivates either the analysis of CMB measurements that do not suffer from this anomaly (Di Valentino & Melchiorri 2022) or the search for additional and independent cosmological probes that are sufficiently sensitive to neutrino masses.

Recent results from other CMB experiments, such as the South Pole Telescope (Dutcher et al. 2021) or the Atacama Cosmology Telescope (Aiola et al. 2020; Madhavacheril et al. 2024), have yielded relatively stringent upper limits on ∑m_ν, but usually through a combination with Planck or Wilkinson Microwave Anisotropy Probe data, with the exception of the results from the SPT (Balkenhol et al. 2021), where a 95% upper limit of m_ν < 0.30 eV was found from a joint analysis with baryon acoustic oscillation data. In regards to weak lensing, experiments such as the Kilo-Degree Survey (de Jong et al. 2013) or the Dark Energy Survey (Dark Energy Survey Collaboration 2016), which perform cosmic shear measurements along with galaxy clustering and galaxy–galaxy lensing, usually find either very loose constraints on neutrino masses (Tröster et al. 2021) or no constraints at all when the data are not combined with external probes. However, a joint analysis of the Dark Energy Survey Year 3 results with baryon acoustic oscillations, redshift-space distortions, IA supernovae, and Planck CMB data yielded an upper limit of ∑m_ν < 0.20 eV at 95% (Abbott et al. 2023).

The submillimeter galaxy magnification bias has recently been proposed as a novel approach to constrain cosmology via a weak-lensing-induced cross-correlation between a foreground galaxy sample and a background set of submillimeter galaxies (Bonavera et al. 2020, 2021; González-Nuevo et al. 2021). Indeed, the phenomenon of magnification bias (see Bartelmann & Schneider 2001, and references therein) can boost the flux of background sources while increasing the solid angle they subtend. However, imposing a flux threshold effectively creates a mismatch between the two effects, which results in an excess of background sources around those in the foreground with respect to the absence of lensing. Although traditionally deemed inferior to shear analyses for the probing of the galaxy–matter cross-correlation, the realization that submillimeter galaxies provide an optimal background sample for magnification bias studies (as shown by the very significant early detections of this cross-correlation in Wang et al. 2011; González-Nuevo et al. 2014) has turned this observable into an emerging independent cosmological probe. Therefore, the present paper proposes the use of the submillimeter galaxy magnification bias to constrain the sum of neutrino masses. This observable has been shown to be relatively sensitive to cosmological parameters such as Ω_m and σ₈, although by itself it does not currently seem to be able to constrain others, such as the baryon density parameter, Ω_b, or the spectral index of the primordial power spectrum, n_s. However, the realization that choosing a cosmological parametrization based on A_s rather than one in terms of σ₈ can induce very different sensitivities to neutrino masses (while keeping the observable insensitive to the aforementioned unconstrained parameters) motivates us to study its potential to provide a bound on ∑m_ν. To this end, we measured the angular cross-correlation function between a sample of high-redshift background submillimeter galaxies from H-ATLAS (Pilbratt et al. 2010; Eales et al. 2010) and a sample of moderate-redshift foreground galaxies from GAMA II (Driver et al. 2011; Baldry et al. 2010, 2014; Liske et al. 2015). Assuming a flat ΛCDM cosmology with massive neutrinos, we adopted a modified halo model for the nonlinear galaxy-matter cross-power spectrum to derive the posterior probability distribution of the sum of neutrino masses as well as additional halo occupation distribution (HOD) and cosmological parameters.

This paper is structured as follows. Section 2 lays out the theoretical background of this work. We discuss the subtleties that arise with the introduction of massive neutrinos in the cosmological setting and describe in detail how the halo model of structure formation should be modified to account for neutrino masses. The methodology we followed is described in Sect. 3, where we present the galaxy samples, the procedure we used to estimate the angular cross-correlation function, and the MCMC setup we used to sample the posterior probability distribution of the parameters involved. Section 4 presents the results we obtained and in Sect. 5 we summarize our conclusions.

2. Theoretical framework

2.1. Neutrino masses in the cosmological setting

Accommodating neutrino masses into the cosmological framework requires an extension of the minimal six-parameter ΛCDM model. The Friedmann equation for a flat ΛCDM universe in the presence of massive neutrinos reads

$$\begin{array}{c}\hfill \frac{H^2(z)}{H_0^2}=[({\mathrm{\Omega }}_{\mathrm{cdm}}+{\mathrm{\Omega }}_{\mathrm{b}}){(1+z)}^3+{\mathrm{\Omega }}_{\gamma }{(1+z)}^4+{\mathrm{\Omega }}_{\mathrm{\Lambda }}+\frac{{\rho }_{\nu }(z)}{{\rho }_{\mathrm{crit},0}}],\end{array}$$

where

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_i\equiv \frac{{\rho }_{i,0}}{{\rho }_{\mathrm{crit},0}}\end{array}$$

is the ratio of the energy density of cosmological species at present i to the critical density, where i runs across cold dark matter (cdm), baryons (b), photons (γ), a cosmological constant (Λ), and neutrinos (ν), and ρ_ν(z)≡∑_jρ_νj(z), where the index j denotes each neutrino mass eigenstate. As neutrinos are the only known particles that undergo a nonrelativistic transition, the redshift scaling of their energy density depends on their mass. Indeed, a neutrino of mass m_νi can be shown to become nonrelativistic at a redshift (Lesgourgues et al. 2013)

$$\begin{array}{c}\hfill 1+z_{{\mathrm{nr}}_i}\approx \frac{m_{{\nu }_i}}{5.28\times {10}^{-4}\mathrm{eV}},\end{array}$$

after which it will contribute to the energy density of matter. Given the experimental data from flavor oscillations (Bilenky 2016), at least two out of the three neutrino states have been nonrelativistic for a long time, which means that the Friedmann equation can be written at our redshifts of interest (0.2 − 0.8) as

$$\begin{array}{c}\hfill \frac{H^2(z)}{H_0^2}=[({\mathrm{\Omega }}_{\mathrm{cdm}}+{\mathrm{\Omega }}_{\mathrm{b}}+{\mathrm{\Omega }}_{\nu }){(1+z)}^3+{\mathrm{\Omega }}_{\gamma }{(1+z)}^4+{\mathrm{\Omega }}_{\mathrm{\Lambda }}].\end{array}$$

Indeed, even if the third neutrino was very light and thus still relativistic at these redshifts, its contribution to Ω_γ would be negligible, and therefore the above equation is a very good approximation.

Relating the neutrino density parameter, Ω_ν, to the sum of their mass, ∑m_ν ≡ ∑_im_νi, requires careful consideration of their thermal history (taking into account that their decoupling from the cosmic plasma is not instantaneous given the proximity in time to e⁺e⁻ annihilation); this yields (Lesgourgues et al. 2013)

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\nu }{h}^2=\frac{\sum m_{{\nu }_i}}{93.14\mathrm{eV}},\end{array}$$(1)

which is a very good approximation for the same reason as in the previous paragraph. As cosmological observables depend on the sum of neutrino masses at first order (Lesgourgues et al. 2013), the minimal extension of the ΛCDM model that includes neutrino masses is thus a mixed dark-matter model with ∑m_ν as an additional parameter.

The cosmological effects that arise from the introduction of neutrino masses – beyond the background evolution – have been extensively studied (Lesgourgues & Pastor 2006; Lesgourgues et al. 2013; Costanzi et al. 2013; Castorina et al. 2014; LoVerde 2014; Massara et al. 2014). The small (but nonzero) mass of neutrinos induces very large thermal velocities, which in turn sets a free streaming scale that, in principle, prevents their clustering within dark matter halos. In reality, neutrinos from the low-velocity tail of the momentum distribution are able to cluster within the potential well of a cold dark matter halo. The strength and scale of the clustering is set by the mass of the neutrinos, the mass of the cold-dark-matter host halo, and redshift (Singh & Ma 2003; Ringwald & Wong 2004; Brandbyge et al. 2010; Villaescusa-Navarro et al. 2013).

In cosmological terms, neutrino-free streaming leads to a slowdown and a suppression of the growth of matter perturbations on scales smaller than the free-streaming scale (Lesgourgues & Pastor 2006). Therefore, power spectra involving the matter field are modified with respect to the absence of massive neutrinos, and the effect, which depends at linear order on the sum of neutrino masses, will transfer onto observables probing the correlation of matter with itself or with other tracers (such as galaxies, as in our case). While the effect can be quantified analytically in the linear regime, N-body simulations are usually needed to describe the matter and neutrino power spectra in the nonlinear regime (Viel et al. 2010; Villaescusa-Navarro et al. 2013; Castorina et al. 2015; Zennaro et al. 2019). Another possibility, first introduced by Massara et al. (2014), is the analytical description (in the presence of massive neutrinos) of the power spectra in the nonlinear regime using the halo model (Cooray & Sheth 2002; Asgari et al. 2023), which is the procedure followed throughout this paper, as detailed in the following subsection.

2.2. The cross-correlation function and the halo model

As explained in Bonavera et al. (2020), González-Nuevo et al. (2021), and Cueli et al. (2022), the phenomenon of magnification bias probes the galaxy–mass correlation via the weak-lensing-induced angular cross-correlation function between two samples of galaxies with nonoverlapping redshift distributions. Under the Limber and flat-sky approximations, this can be expressed as (Cooray & Sheth 2002)

$$\begin{array}{cc}\hfill w_{\mathrm{fb}}(\theta )& =2(\beta -1){\displaystyle {\int }_0^{\infty }}\frac{\mathrm{d}z}{{\chi }^2(z)}\frac{\mathrm{d}{N}_{\mathrm{f}}}{\mathrm{d}z}{W}^{\mathrm{lens}}(z)\hfill \\ \hfill & {\displaystyle {\int }_0^{\infty }}\mathrm{d}l\frac{l}{2\pi }{P}_{\mathrm{g}-\mathrm{m}}(l/\chi (z),z)J_0(l\theta ).\hfill \end{array}$$(2)

In the above equation,

$$\begin{array}{c}\hfill W^{\mathrm{lens}}(z)\equiv \frac{3}{2}\frac{1}{c^2}[\frac{H(z)}{1+z}{]}^2{\displaystyle {\int }_z^{\infty }}\mathrm{d}{z}^{\prime }\frac{\chi (z)\chi (z^{\prime }-z)}{\chi (z^{\prime })}\frac{\mathrm{d}{N}_{\mathrm{b}}}{\mathrm{d}{z}^{\prime }},\end{array}$$(3)

where dN_b/dz (dN_f/dz) is the unit-normalized background (foreground) source distribution, χ(z) is the co-moving distance at redshift z, J₀ is the zeroth-order Bessel function of the first kind, and β is the logarithmic slope of the background source number counts.

The galaxy-matter cross power spectrum, P_g − m(k, z), has been computed within the halo model formalism, as in previous related works. However, in contrast to these works, which assumed a minimal flat ΛCDM cosmology, modifications to the underlying framework are needed for this paper given the existence of a number of subtleties in the massive neutrino setup, as explained below.

As the matter overdensity field is a weighted average of the cold dark matter + baryon and the neutrino fields, the galaxy-matter cross-power spectrum is given by

$$\begin{array}{c}\hfill P_{\mathrm{g}-\mathrm{m}}(k,z)=(1-f_{\nu })P_{\mathrm{g}-\mathrm{c}}(k,z)+f_{\nu }{P}_{\mathrm{g}-\nu }(k,z),\end{array}$$(4)

where f_ν = Ω_ν/Ω_m and “c” denotes the “cold” field (i.e., the weighted average of the cold dark matter and baryon fields). We now turn to the computation of each of the two terms in the above equation within the halo model, although the dominant contribution comes from the first one.

Following the discussion of Massara et al. (2014), both the halo mass function and the linear halo bias – which are crucial ingredients of the halo model – need to be modified to account for neutrino masses. Indeed, it is not clear a priori which field (i.e., the total matter field or only the cold field) should be used to define the number density of cold halos. However, the fact that only the cold field should be used is supported by the works of Castorina et al. (2014) and Villaescusa-Navarro et al. (2014), which show that following this prescription implies a more universal halo mass function and a more scale-independent halo bias at large scales, respectively. Furthermore, the total-matter power spectrum was better reproduced in Massara et al. (2014) when using – following their nomenclature – this “cold dark matter prescription” over that involving the full matter field. This is due to the fact that very few neutrinos are bound by cold halos, which explains why including them all via the so-called matter prescription yields different results. The cold halo mass function is given by

$$\begin{array}{c}\hfill n(M_{\mathrm{c}},z)=\frac{{\overline{\rho }}_{\mathrm{c}}}{M_{\mathrm{c}}^2}f(\nu )|\frac{\mathrm{d}log\nu }{\mathrm{d}log{M}_{\mathrm{c}}}|,\end{array}$$

where ${\overline{\rho }}_{\mathrm{c}}$ is the co-moving cold dark matter + baryon density, f(ν) defines the chosen halo mass function model, and

$$\begin{array}{c}\hfill \nu (M,z)\equiv [\frac{{\widehat{\delta }}_{\mathrm{crit}}(z)}{{\sigma }_{\mathrm{c}}(M,z)}{]}^2.\end{array}$$

In the above expression, ${\widehat{\delta }}_{\mathrm{crit}}(z)$ is the linear critical overdensity at redshift z, computed via the fit of Kitayama & Suto (1996), and

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{c}}^2(M,z)\equiv {\displaystyle {\int }_0^{\infty }}\frac{\mathrm{d}k}{2{\pi }^2}{k}^2{\stackrel{\sim }{W}}^2(kR)P_{\mathrm{c}}^{\mathrm{lin}}(k,z)\end{array}$$(5)

is the variance of the smoothed linear cold overdensity field, where $\widehat{W}$ is the Fourier transform of the filter function (taken to be a top hat in real space) and $P_{\mathrm{c}}^{\mathrm{lin}}(k,z)$ is the linear cold matter power spectrum computed via the Boltzmann code CLASS (Blas et al. 2011).

Given the above discussion, the halo model prescription for the cold matter power spectrum for the galaxy reads

$$\begin{array}{c}\hfill P_{\mathrm{g}-\mathrm{c}}(k,z)=P_{\mathrm{g}-\mathrm{c}}^{1\mathrm{h}}(k,z)+P_{\mathrm{g}-\mathrm{c}}^{2\mathrm{h}}(k,z),\end{array}$$

where

$$\begin{array}{cc}\hfill P_{\mathrm{g}-\mathrm{c}}^{1\mathrm{h}}(k,z)& ={\displaystyle {\int }_0^{\infty }}\mathrm{d}{M}_{\mathrm{c}}n(M_{\mathrm{c}},z)\frac{M_{\mathrm{c}}}{{\overline{\rho }}_{\mathrm{c}}}\frac{{⟨N_{\mathrm{c}}⟩}_{M_{\mathrm{c}}}}{{\overline{n}}_{\mathrm{g}}(z)}|u_{\mathrm{c}}(k|M_{\mathrm{c}},z)|\hfill \\ \hfill & +{\displaystyle {\int }_0^{\infty }}\mathrm{d}{M}_{\mathrm{c}}n(M_{\mathrm{c}},z)\frac{M_{\mathrm{c}}}{{\overline{\rho }}_{\mathrm{c}}}\frac{{⟨N_{\mathrm{s}}⟩}_{M_{\mathrm{c}}}}{{\overline{n}}_{\mathrm{g}}(z)}|u_{\mathrm{g}}(k|M_{\mathrm{c}},z)||u_{\mathrm{c}}(k|M_{\mathrm{c}},z)|,\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill P_{\mathrm{g}-\mathrm{c}}^{2\mathrm{h}}(k,z)& =P_{\mathrm{c}}^{\mathrm{lin}}(k,z)[{\displaystyle {\int }_0^{\infty }}\mathrm{d}{M}_{\mathrm{c}}{M}_{\mathrm{c}}\frac{n(M_{\mathrm{c}},z)}{{\overline{\rho }}_{\mathrm{c}}}{b}_1(M_{\mathrm{c}},z)u_{\mathrm{c}}(k|M_{\mathrm{c}},z)]\hfill \\ \hfill & [{\displaystyle {\int }_0^{\infty }}\mathrm{d}{M}_{\mathrm{c}}\frac{n(M_{\mathrm{c}},z)}{{\overline{n}}_{\mathrm{g}}(z)}{b}_1(M_{\mathrm{c}},z)({⟨N_{\mathrm{c}}⟩}_{M_{\mathrm{c}}}+{⟨N_{\mathrm{s}}⟩}_{M_{\mathrm{c}}}{u}_{\mathrm{g}}(k|M_{\mathrm{c}},z))].\hfill \end{array}$$

In the above expressions, b₁(M_c, z) is the linear halo bias (computed from the halo mass function via the peak-background split), ${\overline{n}}_{\mathrm{g}}(z)$ is the mean galaxy number density, and

$$\begin{array}{c}\hfill {⟨N⟩}_{M_{\mathrm{c}}}={⟨N_{\mathrm{cen}}⟩}_{M_{\mathrm{c}}}+{⟨N_{\mathrm{sat}}⟩}_{M_{\mathrm{c}}}=[1+(\frac{M_{\mathrm{c}}}{M_1}{)}^{\alpha }]\mathrm{\Theta }(M_{\mathrm{c}}-M_{\min })\end{array}$$

is the mean number of galaxies in a cold halo of mass M_c, split into the contribution of central and satellites according to the three-parameter model of Zehavi et al. (2005). Moreover, u_c and u_g are the normalized Fourier transforms of the cold dark matter halo density profile and the satellite galaxy distribution, respectively, which are assumed to follow Navarro–Frenk–White models.

Regarding the computation of the second term in Eq. (4), that is, the galaxy-neutrino cross-power spectrum, it can be further decomposed following Massara et al. (2014) by splitting the neutrino density field into the “unclustered” – or linear theory (${\delta }_{\nu }^{\mathrm{L}}$) – and clustered (${\delta }_{\nu }^{\mathrm{h}}$) neutrino fields, that is,

$$\begin{array}{c}\hfill P_{\mathrm{g}-\mathrm{m}}(k,z)=(1-f_{\nu })P_{\mathrm{g}-\mathrm{c}}(k,z)+f_{\nu }[(1-F_{\mathrm{h}})P_{\mathrm{g}-\nu }^{\mathrm{L}}(k,z)+F_{\mathrm{h}}{P}_{\mathrm{g}-\nu }(k,z)],\end{array}$$

where F_h denotes the fraction of neutrinos that do cluster within cold halos. As both f_ν and F_h are small², it is clear that the galaxy-matter power spectrum is dominated by the cross-power spectrum between galaxies and the cold field, followed by a subdominant contribution from the cross-power spectrum between galaxies and the linearly evolving neutrino field and a negligible term from the neutrino cross-power spectrum for the galaxy. Of the last two, the former can be approximated as

$$\begin{array}{c}\hfill P_{\mathrm{g}-\nu }^{\mathrm{L}}(k,z)\approx \sqrt{P_{\mathrm{gg}}(k,z)P_{\nu \nu }^{\mathrm{L}}(k,z)},\end{array}$$

and computed via the halo model and the Boltzmann code CAMB (Lewis & Challinor 2011). The latter could in principle be written within the massive neutrino halo model formalism, for which several assumptions and a neutrino density profile are needed. As this term was negligible for our purposes, its contribution was ignored.

3. Methodology

3.1. Galaxy samples

The foreground and background galaxy samples were extracted from the GAMA II (Driver et al. 2011; Baldry et al. 2010, 2014; Liske et al. 2015) and H-ATLAS (Pilbratt et al. 2010; Eales et al. 2010) surveys, respectively. Their common area covers three regions on the celestial equator at 9, 12, and 14.5 h (G09, G12, and G15) and part of the south Galactic pole, which amounts to a total of ∼207 deg². These are the same samples used in most previous related works (see Bonavera et al. 2024, and references therein for further details).

The foreground sample is made up of GAMA II sources with spectroscopic redshifts in the range 0.2 < z < 0.8, resulting in ∼130 000 galaxies with a median redshift of 0.28 and surveyed down to an r-band magnitude of 19.8. The associated redshift distribution is shown in purple in Fig. 1. The background sample is made up of ∼37 000 H-ATLAS reliably detected sources obtained via a photometric redshift selection of 1.2 < z < 4.0 to ensure that there is no overlap with the foreground galaxies. The redshift distribution of the background sample, taking random errors into account, is shown in Fig. 1 (in red).

Fig. 1. Normalized redshift distribution of the background (in purple) and foreground (in red) samples of galaxies.

3.2. Cross-correlation measurement

The measurement procedure is explained in detail in Cueli et al. (2023), but we briefly summarize the idea below. We performed a sole measurement of the angular cross-correlation function over the entire available area via the natural modification of the Landy–Szalay estimator (Landy & Szalay 1993):

$$\begin{array}{c}\hfill \stackrel{\sim }{w}(\theta )=\frac{{\mathrm{D}}_{\mathrm{f}}{\mathrm{D}}_{\mathrm{b}}(\theta )-{\mathrm{D}}_{\mathrm{f}}{\mathrm{R}}_{\mathrm{b}}(\theta )-{\mathrm{D}}_{\mathrm{b}}{\mathrm{R}}_{\mathrm{f}}(\theta )+{\mathrm{R}}_{\mathrm{f}}^i{\mathrm{R}}_{\mathrm{b}}(\theta )}{{\mathrm{R}}_{\mathrm{f}}{\mathrm{R}}_{\mathrm{b}}(\theta )},\end{array}$$(6)

where X_fY_b(θ) is the number of foreground-background galaxy pairs at an angular distance θ; when X ≡ D, the galaxies are chosen from the data, whereas X ≡ R implies that the sources are selected from a randomly (unclustered) generated catalog. The associated covariance matrix was computed internally through a Bootstrap resampling method with an oversampling factor of 3 (Norberg et al. 2009), that is,

$$\begin{array}{c}\hfill \mathrm{Cov}({\theta }_i,{\theta }_j)=\frac{1}{N_{\mathrm{b}}-1}{\displaystyle \sum _{k=1}^{N_{\mathrm{b}}}}[\widehat{w_k}({\theta }_i)-\overline{\widehat{w}}({\theta }_i][\widehat{w_k}({\theta }_j)-\overline{\widehat{w}}({\theta }_j],\end{array}$$(7)

where ${\widehat{\mathit{w}}}_k$ denotes the cross-correlation measurement from the kth Bootstrap sample, $\overline{\widehat{\mathit{w}}}$ is the corresponding average over all Bootstrap samples, and N_b = 10 000. The measurements are plotted in Fig. 2 in black.

Fig. 2. Angular cross-correlation measurements using all four fields (black) and excluding the G15 region (orange).

3.3. Parameter estimation

A Bayesian approach was adopted for the estimation of the probability distribution of the involved parameters. We carried out several MCMC analyses using the open source emcee software package (Foreman-Mackey et al. 2013), which is a Python implementation of the Goodman & Weare so-called affine invariant MCMC ensemble sampler (Goodman & Weare 2010). The log-likelihood functions was assumed to be a multivariate Gaussian distribution, that is,

$$\begin{array}{cc}\hfill log\mathcal{L}({\theta }_1,\dots ,{\theta }_{\mathrm{m}})=-\frac{1}{2}& [mlog(2\pi )+log|C|+{\overrightarrow{\epsilon }}^{\mathrm{T}}{C}^{-1}\overrightarrow{\epsilon }],\hfill \end{array}$$(8)

where

$$\begin{array}{c}\hfill \overrightarrow{ϵ}\equiv [\epsilon ({\theta }_1),\dots ,\epsilon ({\theta }_{\mathrm{m}})]\end{array}$$

and $\epsilon ({\theta }_i)={\mathit{w}}_{\mathrm{fb}}({\theta }_i)-{\widehat{\mathit{w}}}_{\mathrm{mcross}}({\theta }_i)$.

Several different cases were studied in this work according to the parameterization of the model (see the following section), but the estimation procedure involved the same parameters in almost all scenarios, namely the HOD (α, M_min and M₁) and a massive neutrino ΛCDM cosmology with a fixed power spectrum normalization (Ω_m, h and ∑m_ν), along with the number count logarithmic slope, β. The prior distributions used for all cases are summarized in Table 1. They are all uniform with the exception of the logarithmic slope of the background number counts, which was chosen to follow a Gaussian distribution around 2.9 following the analysis of Cueli et al. (2023).

4. Results

4.1. Sensitivity to neutrino masses

It is well known that the minimal ΛCDM cosmological model can be described through the six-parameter family (Ω_m, Ω_b, n_s, σ₈, h, τ), where σ₈ fixes the normalization of the z = 0 linear matter power spectrum. However, this model can also be parameterized by means of the equivalent set (Ω_m, Ω_b, n_s, A_s, h, τ), where A_s is the amplitude of the primordial spectrum of curvature perturbations at a pivot scale k₀, that is³

$$\begin{array}{c}\hfill \frac{k^3}{2{\pi }^2}{P}_{\zeta }(k)=A_{\mathrm{s}}(\frac{k}{k_0}{)}^{n_{\mathrm{s}}-1}.\end{array}$$

These two normalization parameters (σ₈ and A_s) can be linked given the relation between the linear matter spectrum and the primordial spectrum of curvature perturbations (Dodelson & Schmidt 2020):

$$\begin{array}{c}\hfill P_{\mathrm{mm}}^{\mathrm{lin}}(k,z)=\frac{8{\pi }^2}{25}\frac{A_{\mathrm{s}}{k}_0}{{\mathrm{\Omega }}_{\mathrm{m}}^2}{D}^2(z)T^2(k,z)(\frac{k}{k_0}{)}^{n_{\mathrm{s}}}(\frac{c}{H_0}{)}^4.\end{array}$$

In the above equation, the transfer function T(k, z) satisfies T(k → 0, z)→1 and is redshift-dependent in the massive-neutrino setup⁴. In turn, the linear growth factor is given by

$$\begin{array}{c}\hfill D(z)=\frac{5}{2}{\mathrm{\Omega }}_{\mathrm{m}}\frac{H(z)}{H_0}{\displaystyle {\int }_z^{\infty }}\frac{1+z^{\prime }}{{(H(z^{\prime })/H_0)}^3}\mathrm{d}{z}^{\prime }.\end{array}$$(9)

The addition of massive neutrinos does not modify the above equations, but affects the computation of both the growth factor and the transfer function. This entire discussion is related to the fact that, interestingly, the choice of parameter family has a large impact on the constraining power of our observable regarding neutrino masses. Indeed, the left panel of Fig. 3 shows the sensitivity of the theoretical cross-correlation function to changes in ∑m_νi of up to 0.80 eV in the σ₈ parameterization. It is clear that the variation is small and only noticeable at the largest angular scales, where the error bars are still sizeable in comparison. Therefore, an attempt to constrain neutrino masses, even by fixing the σ₈ parameter, is almost certain to produce nonconstraining results, as we show below.

Fig. 3. Sensitivity of the angular cross-correlation function with respect to ∑mν for fixed values of σ8 = 0.79 (left panel) and log(1010 As) = 3.35 (right panel).

The situation in the parameterization that includes A_s is completely different. The right panel of Fig. 3 depicts the sensitivity of the cross-correlation function to the same variations as the previous case, but this time in the A_s parameterization. The response of the observable is now substantial and predominantly takes place at intermediate angular scales, where the error bars are smaller. Although this does not necessarily mean that a full MCMC analysis including A_s instead of σ₈ will compare favorably with the previous case, it does hint at the possibility of restricting neutrino masses when the value of A_s is fixed. This finding is the main motivation behind the present paper.

4.2. Cosmological analysis: A_s vs. σ₈

The first case we studied involved fixing the σ₈ parameter to 0.79, the best fit from the massless neutrino case in Cueli et al. (2023), whereas the rest of parameters (α, M_min, M₁, Ω_m, h, ∑m_ν, β) were allowed to vary. The corresponding statistical results are summarized in the first column of Table 2 and Fig. A.1 (in red). Concerning the HOD, we obtained results in agreement with typical values from the literature (Zehavi et al. 2005, 2011; Abbas et al. 2010) and with the findings of the massless neutrino case from Cueli et al. (2023); we derived mean values of $log{M}_{\min }=11.{64}_{-0.12}^{+0.17}$ and $log{M}_1=13.{39}_{-0.92}^{+0.71}$ for the halo masses and $\alpha =0.{69}_{-0.46}^{+0.30}$ for the satellite logarithmic slope.

Regarding cosmology, we obtained a mean value of ${\mathrm{\Omega }}_{\mathrm{m}}=0.{28}_{-0.08}^{+0.04}$ for the matter density parameter, higher than the previous results by Cueli et al. (2023) and Bonavera et al. (2024), where neutrinos were considered massless. We believe this displacement toward larger values is caused by the additional large-scale freedom provided by neutrino masses, since both parameters affect mainly the high-θ regime. Indeed, as seen on the Ω_m − ∑m_ν plane in Fig. 4 (in red), neutrino masses close to zero imply values of the matter density parameter aligned with ∼0.2, as found by Cueli et al. (2023) and Bonavera et al. (2024). Only an upper limit can be found for the Hubble constant, with h < 0.70 at 68% and, as expected from the analysis in the previous subsection, the neutrino mass sum cannot be constrained in this scenario.

Fig. 4. Marginalized posteriors and probability contours of the cosmological parameters for fixed σ8 (in red) and As (in blue) values according to the best fit from the massless neutrino case of CUE24a.

We next moved to the parametrization of the model where σ₈ is replaced by A_s, which, as in the previous case, is fixed to the best-fit value from the massless neutrino case, which is 2.86 × 10⁻⁹. The statistical results are summarized in the second column of Table 2 and Fig. 4 (in blue). The HOD masses are constrained, with mean values of $log{M}_{\min }=11.{47}_{-0.11}^{+0.14}$ and $log{M}_1=12.{89}_{-0.41}^{+0.56}$, which are in agreement with the previous case within the error bars. The logarithmic slope of the number of satellites, α, is however not as well constrained due to the small-scale influence of neutrino masses, but we can derive a lower limit of α > 0.79 at 68%. It should be noted that, although the fixed parameter in each of the two scenarios has been fixed according to the same best-fit cosmology, A_s and σ₈ are not in a one-to-one relation, as one is a function of all cosmological parameters and not only of the other. Therefore, the quantitative results from the two cases cannot be compared directly, because they are not equivalent given the dependence on the exact value of the normalization parameter.

As far as cosmology is concerned, as seen in Fig. 4 (in blue), we obtained a well-constrained posterior distribution for the matter density parameter, with a mean value of ${\mathrm{\Omega }}_{\mathrm{m}}=0.{18}_{-0.02}^{+0.02}$ and a lower limit of h > 0.77 at 95% for the Hubble constant. These findings, which are aligned with the results from Cueli et al. (2023) and Bonavera et al. (2024) when all fields are considered, are clearly in tension with the results of other cosmological probes (Abbott et al. 2018, 2020; Planck Collaboration VI 2020) and are further discussed in the following subsection.

Regardless, and contrary to the prior scenario, a clear upper limit is now found for the sum of neutrino masses, with ∑m_ν < 0.22 (0.78) eV at 68% (95%) credibility. This result demonstrates the main finding of this paper: by fixing the matter power spectrum normalization, current submillimeter galaxy magnification bias measurements are sensitive to neutrino masses, meaning that a cosmological paramaterization based on A_s as opposed to σ₈ must be adopted.

4.3. Further discussion

In light of the results of the previous subsection, a number of aspects need to be explored regarding their validity. First, one may wonder how the above findings behave with respect to the prior distribution of A_s, as they have been derived assuming the best-fit value to the massless neutrino case. To investigate this matter further, we started by considering a fixed value of A_s according to Planck’s best fit (TT, EE, TE, lowE + lensing), that is, log(10¹⁰A_s) = 3.044 (Planck Collaboration VI 2020).

The statistical results are gathered in the third column of Table 2 and Fig. A.2. Although there are slight variations in the overall HOD, the posterior distributions of M_min and h remain unchanged. Furthermore, as seen in Fig. 5, the sensitivity of the system with respect to neutrino masses stays qualitatively unaltered; indeed, a clear upper limit of ∑m_ν < 0.36 (1.08) eV is still found at 68% (95%). The matter density parameter, however, is displaced toward slightly lower values due to the degeneracy with A_s (larger values of one parameter can be balanced by smaller values of the other).

Fig. 5. Marginalized posteriors and probability contours of the cosmological parameters for fixed As values according to the best fit from CUE24a (in blue) and Planck (in magenta).

Second, and given the current size of the error bars, the question arises as to how wide the prior distribution of A_s can be whilst maintaining the sensitivity of the model to neutrino masses. To investigate this question, we considered a Gaussian prior for A_s around the above Planck value with a standard deviation of 0.1. The results are summarized in the fourth column of Table 2 and Fig. 6 in cyan. Our findings are both qualitatively and quantitatively similar to the previous case, with only slightly larger uncertainties arising from the additional degree of freedom. The posterior distribution of A_s is displaced from its prior to a mean value of $log({10}^{10}{A}_{\mathrm{s}})=3.{17}_{-0.10}^{+0.11}$. This behavior seems to impact the sensitivity to neutrino masses, as seen more clearly in Fig. 6; indeed, although the highest posterior density interval at 68% is ∑m_ν < 0.46 eV, a secondary high-density region is found in the posterior at very high ( ∼ 1.3 eV) neutrino masses. However, this still indicates that neutrino masses are not prior-dominated even in the case of a relatively narrow Gaussian prior for A_s and therefore that sensitivity is still manifest.

Fig. 6. Marginalized posteriors and probability contours of the cosmological parameters for a fixed Planck As value (in magenta) and for a Gaussian prior around that value (in cyan).

All analyses performed up to this point under the A_s parametrization share a common feature: the posterior distributions of Ω_m (and, to a certain extent, that of h) are in tension with common values from the literature. The above tests, however, suggest that the exact prior value of A_s should not be able to account for the entirety of this behavior. Indeed, we believe that the discrepant cosmological results stem from the data themselves.

As already shown in Cueli et al. (2023) and Bonavera et al. (2024), there appears to be an excess of a cross-correlation signal – especially at the largest scales – due the G15 equatorial region. As the galaxy selection criteria are uniform across all four fields (as are the redshift distributions), we hold the view that this stems from the phenomenon of sampling variance and plan to study it further with a larger sample in a follow-up work. In any case, a preliminary analysis performed in Cueli et al. (2023) showed that removing the G15 region induced non-negligible changes in cosmological constraints; in particular, larger values of Ω_m were attained, which are in keeping with standard results from external probes.

Therefore, we decided to assess the influence of leaving out the G15 region in our analysis with neutrino masses. To do so, we fixed A_s once again to the best fit from the massless neutrino case of Cueli et al. (2023). The cross-correlation measurements for this case are shown in Fig. 2 in orange, where the excess of cross-correlation at intermediate and large scales is more visible when compared to all fields.

The results are summarized in the fifth column of Table 2 and Fig. A.4. All three HOD parameters are reasonably well constrained and their distributions are only slightly modified with respect to the four-field case. However, as shown by Fig. 7, the posterior distribution of Ω_m is strongly displaced toward larger values, following the behavior of Cueli et al. (2023) and Bonavera et al. (2024), with a mean of ${\mathrm{\Omega }}_{\mathrm{m}}=0.{34}_{-0.03}^{+0.02}$. Moreover, the Hubble constant is constrained on its two sides, with a clear mean value of $h=0.{72}_{-0.06}^{+0.07}$. However, the sensitivity to neutrino masses is lost, and no constraints can be given in this scenario given the sizeable error bars arising from the exclusion of the G15 region, as confirmed by a simple parameter sensitivity analysis. As expected, the qualitative behavior of this last test remains for a fixed Planck value of A_s.

Fig. 7. Marginalized posteriors and probability contours of the cosmological parameters for a fixed As value according to the best fit from CUE24a. The results using all four fields are shown in blue, while the case where the G15 region was excluded is depicted in orange.

The results without the G15 field are thus in tension with the findings using all four independent regions, and this is even more noticeable than in the massless neutrino scenario. This is still a matter of ongoing investigation and will be addressed further in future work. Nonetheless, as far as neutrino masses are concerned, we find that sacrificing one field has a devastating effect on the sensitivity of cosmic magnification, although it produces cosmological results that are well in agreement with external probes. Therefore, measurement error bars must be at least the size of the four-field case, especially at the smallest scales.

Regarding the validity of the halo model, it is well known (Asgari et al. 2023) that it fails to faithfully reproduce the transition between the one- and two-halo regimes, which can be alleviated by the introduction of nonlinear halo bias, as calibrated by numerical simulations (Nishimichi et al. 2019; Mead & Verde 2021). Moreover, baryonic feedback can induce deviations to the matter and galaxy-matter power spectra with respect to the minimal halo model. Its impact on small scales can be addressed (at least partially) by the assignment of a variable halo concentration amplitude or a separate concentration amplitude for the satellite galaxy distribution (Amon et al. 2023; Dvornik et al. 2023). We will study the introduction of these effects in future work.

5. Summary and conclusions

This paper addresses the sensitivity of the submillimeter galaxy magnification bias to the sum of neutrino masses. We measured the angular cross-correlation function induced by the weak-lensing magnification bias between a sample of high-redshift submillimeter galaxies from H-ATLAS and a sample of moderate-redshift galaxies from GAMA. We describe the nonlinear behavior of the cross-power spectrum of the galaxy matter via the halo model, which allows us to go down to approximately arcminute scales and describe the galaxy–halo connection via a simple three-parameter HOD model. We assume a ΛCDM cosmological model with neutrino masses, for which we implement a modified version of the halo model that entails considering only the cold dark matter + baryon field to describe the nonlinear clustering of matter.

The main goal of this study is to investigate the potential of cosmic magnification with background submillimeter galaxies as a cosmological probe of neutrino masses. We find that a cosmological paramaterization in terms of the amplitude of the primordial power spectrum (A_s) is required, rather than that involving the root mean square of the z = 0 filtered overdensity field (σ₈). The reason is related to the different sensitivity of the angular cross-correlation to the sum of neutrino masses under each paramaterization. Indeed, the variation of the signal is concentrated on the small and intermediate scales, where the error bars are relatively small, for a fixed value of A_s. This is in contrast with the case of a fixed value of σ₈, where the response is located on large angular scales, rendering the current measurements insensitive to neutrino masses.

For a value of A_s fixed to the best fit of the massless neutrino case of Cueli et al. (2023), we obtained clear upper bounds for the sum of neutrino masses: ∑m_ν < 0.22 eV at 68% and < 0.78 at 95%. Although the exact numerical values depend on the chosen value of A_s, we find qualitatively stable results when it is varied. Indeed, for a fixed Planck value of A_s, we obtain ∑m_ν < 0.36 (1.08) eV at 68% (95%). Moreover, assuming a Gaussian prior for A_s of around this Planck value, we find ∑m_ν < 0.46 eV at 68%, with a small high-density region at large neutrino masses due to the additional freedom in A_s. This demonstrates that, by fixing only the normalization of the primordial power spectrum (or by adopting a relatively narrow Gaussian prior distribution for it), cosmic magnification on submillimeter galaxies can be sensitive to the sum of neutrino masses.

However, we find the posterior distributions obtained for the matter density parameter to be in tension with values derived from external probes. This issue was pointed out by Cueli et al. (2023) and Bonavera et al. (2024) and was deemed to be the result of an excess of cross-correlation found within the G15 region. Although the exclusion of this field solves the problem, the larger error bars due to the reduction in galaxy sample size equate to insufficient sensitivity to constrain neutrino masses. Therefore, future work will be directed toward a better understanding of this issue, enabling us to obtain both unbiased and constraining results.

Acknowledgments

LB, JGN, JMC and DC acknowledge the PID2021-125630NB-I00 project funded by MCIN/AEI/10.13039/501100011033/FEDER, UE. LB also acknowledges the CNS2022-135748 project funded by MCIN/AEI/10.13039/501100011033 and by the EU “NextGenerationEU/PRTR”. JMC also acknowledges financial support from the SV-PA-21-AYUD/2021/51301 project. MV is partly supported by the INFN IS INDARK grant.

References

Appendix A: Additional plots

Fig. A.1. Marginalized posterior distributions and probability contours for fixed σ8 (in red) and As (in blue) values according to the best fit from the massless neutrino case of CUE24a.

Fig. A.2. Marginalized posterior distributions and probability contours for a fixed As value according to the best fit from CUE24a (in blue) and to Planck (in magenta).

Fig. A.3. Marginalized posterior distributions and probability contours for a fixed Planck As value (in magenta) and for a Gaussian prior around that value (in cyan).

Fig. A.4. Marginalized posterior distributions and probability contours for a fixed As value according to the best fit from CUE24a. The results using all four fields are shown in blue, while the case where the G15 region was excluded is depicted in orange.

All Tables

All Figures

	Fig. 1. Normalized redshift distribution of the background (in purple) and foreground (in red) samples of galaxies.
In the text

	Fig. 2. Angular cross-correlation measurements using all four fields (black) and excluding the G15 region (orange).
In the text

	Fig. 3. Sensitivity of the angular cross-correlation function with respect to ∑mν for fixed values of σ8 = 0.79 (left panel) and log(1010 As) = 3.35 (right panel).
In the text

	Fig. 4. Marginalized posteriors and probability contours of the cosmological parameters for fixed σ8 (in red) and As (in blue) values according to the best fit from the massless neutrino case of CUE24a.
In the text

	Fig. 5. Marginalized posteriors and probability contours of the cosmological parameters for fixed As values according to the best fit from CUE24a (in blue) and Planck (in magenta).
In the text

	Fig. 6. Marginalized posteriors and probability contours of the cosmological parameters for a fixed Planck As value (in magenta) and for a Gaussian prior around that value (in cyan).
In the text

	Fig. 7. Marginalized posteriors and probability contours of the cosmological parameters for a fixed As value according to the best fit from CUE24a. The results using all four fields are shown in blue, while the case where the G15 region was excluded is depicted in orange.
In the text

	Fig. A.1. Marginalized posterior distributions and probability contours for fixed σ8 (in red) and As (in blue) values according to the best fit from the massless neutrino case of CUE24a.
In the text

	Fig. A.2. Marginalized posterior distributions and probability contours for a fixed As value according to the best fit from CUE24a (in blue) and to Planck (in magenta).
In the text

	Fig. A.3. Marginalized posterior distributions and probability contours for a fixed Planck As value (in magenta) and for a Gaussian prior around that value (in cyan).
In the text

	Fig. A.4. Marginalized posterior distributions and probability contours for a fixed As value according to the best fit from CUE24a. The results using all four fields are shown in blue, while the case where the G15 region was excluded is depicted in orange.
In the text