© The Authors 2024

1. Introduction

The standard cosmological model is based on the assumption that space-time can be well described by the Friedmann-Lemaître-Robertson-Walker (FLRW) metric, which stems from a generalization of the Copernican Principle, stating that our location is not special in the Universe, and also from observations suggesting that, on sufficiently large scales, on the order of 100 Mpc, the Universe is homogeneous and isotropic (Maddox et al. 1990; Cole et al. 2005; Saadeh et al. 2016).

One of the tests of this principle could be achieved by assuming an alternative metric, such as the inhomogeneous metric of the Lemaître-Tolman-Bondi model (LTB), and confronting it with current observations (García-Bellido & Haugbølle 2008; February et al. 2010; Valkenburg et al. 2014; Redlich et al. 2014; Camarena et al. 2021), studying its phenomenology through simulations (Alonso et al. 2012; Marra et al. 2022), or forecasting on the ability of forthcoming surveys to constrain such a deviation from the standard assumption (Amendola et al. 2018; Camarena et al. 2023). In this work, we undertake this task and attempt to determine the precision achievable on the cosmological parameters of the ΛLTB model (i.e., the LTB model endowed with a cosmological constant) by using galaxy cluster abundance data obtained from forthcoming surveys. This follows from previous studies, which show that cluster counts are capable of placing constraints on the standard cosmological model and its extensions (Rapetti et al. 2010; Ilić et al. 2019).

The study of clusters has had a major role in establishing today’s most accepted cosmological model, the cosmological constant and cold dark matter (ΛCDM) model, by pinning down the values of dark matter and dark energy abundances (Bocquet et al. 2019). Clusters’ progenitors are the ingredients of our Universe, associated with the primordial inflationary stage, after which the growth of the density perturbations was amplified by gravity, creating what we now know as large-scale structures (LSSs; Peebles 1982). As such, the abundance of galaxy clusters is very sensitive to the LSS power spectrum and to the growth rate of their density perturbations. On top of that, galaxy cluster abundances are also sensitive both to the cosmological parameters and to the late evolution of the Universe, as these objects were the last to form according to the bottom-up hierarchical formation model (Blumenthal et al. 1984). We need, however, to properly calibrate the mass observable proxy used to relate the mass function, which describes the mass distribution of dark matter halos, to the observed clusters. Therefore, a miscalibration could be one of the sources of the discrepancy found in the value of the matter amplitude fluctuation parameter, σ₈, when constrained by deep probes in comparison to local ones (Sakr et al. 2022). Thus, we discuss later in this study the modeling of the ΛLTB mass function, as well as the calibration of the mass proxy accessed through a scaling relation proper to each survey.

Forecasts on constraints on the ΛLTB model were previously presented in Amendola et al. (2018), and recently in Camarena et al. (2023) using cosmic microwave background anisotropies (Planck Collaboration VI 2020), mock data of baryonic acoustic oscillations (Euclid Collaboration 2020), the Compton-y distortion (Fixsen et al. 1996), and kinetic Sunyaev-Zeldovich data (Reichardt et al. 2021). In this work, we perform a similar study, but aim to benefit from the nonlinear scales that the galaxy clusters additionally probe, expected to be obtained from constructed ‘Stage-IV’-like mock catalogs, such as Euclid (Laureijs et al. 2011; Euclid Collaboration 2022). The latter is an M-class space mission of the European Space Agency (ESA) which will allow the optical detection of clusters previously unattainable in terms of the depth and area covered. Along with it, we also consider clusters detected from the Legacy Survey of Space and Time (LSST) survey, performed at the Vera C. Rubin Observatory (and hereafter named Rubin) (LSST Science Collaboration 2009) and the eROSITA survey X-ray based cluster counts (Pillepich et al. 2012), and the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration 2016a). The former two, using richness as the same mass observable proxy, are complementary in terms of redshift and sky coverage, while eROSITA and DESI will allow us, each with a different detection method – respectively, X-ray and spectroscopy – to further extend the redshift range as well as the number of clusters detected.

Our paper is organized as follows: in Sect. 2 we briefly review the modeling of galaxy cluster counts in ΛLTB. In Sect. 3 we present the surveys’ specifications and scaling relations modeling. Lastly, our results are presented and discussed in Sect. 4 before we conclude in Sect. 5.

2. Modeling galaxy cluster counts in ΛLTB

2.1. The ΛLTB model

The ΛLTB model is a spherically symmetric, inhomogeneous space-time with a cosmological constant, Λ, and cold dark matter (CDM), as in Valkenburg et al. (2014) or Camarena et al. (2023). For the sake of self-consistency, we now summarize the key equations of the model, from the aforementioned papers, which are necessary to understand the galaxy cluster counts methodology detailed in Sect. 2.2.

The ΛLTB line element is of the form

$$\begin{array}{c}\hfill d{s}^2=-c^2\mathrm{d}{t}^2+\frac{R{\prime }^2(t,r)}{1-K(r)r^2}\mathrm{d}{r}^2+R^2(t,r)(\mathrm{d}{\theta }^2+{sin}^2\theta \mathrm{d}{\varphi }^2),\end{array}$$(1)

where c is the speed of light, K(r) is an arbitrary function of the spatial coordinate, r, and a prime denotes a partial derivative with respect to r. If K is constant and R(t, r) = a(t) r, Eq. (1) becomes the FLRW line element where a(t) is the FLRW scale factor. From Eq. (1) one can define longitudinal a_∥ = R′(t, r) and transverse a_⊥ = R(t, r)/r, scale factors, and longitudinal, H_∥, and transverse, H_⊥, expansion rates given by

$$\begin{array}{c}\hfill H_{\Vert }\equiv \frac{\dot{a_{\Vert }}}{a_{\Vert }},H_{\perp }\equiv \frac{{\dot{a}}_{\perp }}{a_{\perp }},\end{array}$$(2)

where the dots represent partial derivatives with respect to cosmic time. Assuming a cosmological constant, the general relativity (GR) equations allow one to derive a FLRW-like equation for the ΛLTB model,

$$\begin{array}{c}\hfill \frac{H_{\perp }^2}{H_{\perp 0}^2}={\mathrm{\Omega }}_{\mathrm{m},0}(\frac{a_{\perp 0}}{a_{\perp }}{)}^3+{\mathrm{\Omega }}_{K,0}(\frac{a_{\perp 0}}{a_{\perp }}{)}^2+{\mathrm{\Omega }}_{\mathrm{\Lambda },0},\end{array}$$(3)

where H_⊥0 ≡ H_⊥(t₀, r), a_⊥0 ≡ a_⊥(t₀, r), and Ω_Λ, 0, Ω_K, 0, Ω_m, 0 are the present-day, t₀, density parameters associated with the cosmological constant, curvature, and matter, respectively. These parameters are functions of the coordinate r and satisfy the closure relation Ω_m, 0(r)+Ω_K, 0(r)+Ω_Λ, 0(r) = 1. Using the GR geodesics equation with the metric Eq. (1) gives

$$\begin{array}{cc}& \frac{\mathrm{d}t}{\mathrm{d}z}=-\frac{1}{(1+z)H_{\Vert }(t,r)},\hfill \end{array}$$(4)

$$\begin{array}{cc}& \frac{\mathrm{d}r}{\mathrm{d}z}=-\frac{c\sqrt{1-K(r)r^2}}{(1+z)a_{\Vert }(t,r)H_{\Vert }(t,r)}.\hfill \end{array}$$(5)

As is seen, Eqs. (4) and (5) are essential in our analysis because their integration provides a relation between r, t, and z, allowing the ΛLTB galaxy cluster’s mass function to be expressed as a function of redshift, which is a key quantity in the modeling of the survey cluster counts described in Sect. 3.

As is shown in Camarena et al. (2023), the ΛLTB model can be parameterized by the arbitrary functions K(r) and m(r) that enter in the definition of

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{m},0}(r)=\frac{2Gm(r)}{H_{\perp 0}^2{a}_{\perp 0}^3{r}^3},{\mathrm{\Omega }}_{K,0}(r)=-\frac{K(r)c^2}{H_{\perp 0}^2{a}_{\perp 0}^2},\end{array}$$(6)

where m(r) is the Euclidean mass profile of the large-scale inhomogeneity

$$\begin{array}{c}\hfill m(r)={\displaystyle {\int }_0^r\mathrm{d}{r}^{\prime }4\pi {\rho }_{\mathrm{m}}(t,r^{\prime })a_{\Vert }{a}_{\perp }^{2}r{\prime }^2,}\end{array}$$(7)

and where ρ_m(t, r) is the matter density. The cosmological constant density parameter is ${\mathrm{\Omega }}_{\mathrm{\Lambda },0}(r)=\mathrm{\Lambda }{c}^2/(3H_{\perp 0}^2)=1-{\mathrm{\Omega }}_{\text{m},\text{0}}(r)-{\mathrm{\Omega }}_{K,0}(r)$. Once the radial coordinate, r, is fixed in such a way that m(r)∝r³, with a normalization described below, the curvature profile, K(r), is the sole free function left in the model, which can be modeled as

$$\begin{array}{c}\hfill K(r)=K_B+(K_C-K_B)P_3(r/r_{\mathsf{B}}),\end{array}$$(8)

where r_B is the comoving radius of the spherical inhomogeneity, K_B is the curvature outside the inhomogeneity, and K_C is the central inhomogeneity curvature. The function P₃ enables us to specify at convenience the radial profile of the central void. We modeled this function, as in Camarena et al. (2023), as

$$\begin{array}{c}\hfill P_3(x)=\{\begin{array}{cc}1-exp[-{(1-x)}^3/x]\hfill & \text{for}0\le x<1\hfill \\ 0\hfill & \text{for}1\le x\hfill \end{array},\end{array}$$(9)

with x = r/r_B. From the density parameters in Eq. (6) and the above closure relation, we determined the missing normalization of the mass, $m(r)={\mathrm{\Omega }}_{\mathrm{m},0}^{\mathrm{out}}{(H_0^{\mathrm{out}})}^2{r}^3/(2G)$, as well as $K_B=-{\mathrm{\Omega }}_{K,0}^{\mathrm{out}}{(H_0^{\mathrm{out}}/c)}^2$, which allowed us to determine the volume-averaged integrated matter density contrast,

$$\begin{array}{c}\hfill \delta (r,t_0)=\frac{3m(r)}{4\pi R^3(t_0,r){\rho }_{\mathrm{m}}^{\mathrm{out}}(t_0)}-1=\frac{{\mathrm{\Omega }}_{\mathrm{m},0}(r)}{{\mathrm{\Omega }}_{\mathrm{m},0}^{\mathrm{out}}}{\left[\frac{H_{\perp 0}(r)}{H_{\perp 0}^{\mathrm{out}}}\right]}^2-1.\end{array}$$(10)

In these expressions the quantities denoted by the superscript “out” are taken at r = r_B. We note that the curvature profile outside the LTB void (r ≥ r_B) is set to K(r) = K_B, and the K_C is directly related to the central under-density contrast δ(r = 0, t₀), as is discussed in Camarena et al. (2023).

2.2. Halo mass function in the ΛLTB model

The formalism adopted in this work is based on the mass function derived from LTB N-body simulations in Alonso et al. (2012). Although the latter deals with the original LTB model, we assume that this strategy is still valid for ΛLTB if the dark energy is accounted for. Indeed, simulations in FRLW metrics show that in general the nonlinear LSS formation in ΛCDM is obtained in models that change the background, by replacing the growth in the new theory Baldi (2012), a fix in line with what we adopted next in our modeling.

In the context of the ΛLTB model, we need to account for the density perturbations at cluster scales inside and outside the ΛLTB void. The variance of the density perturbations at very large radii, outside the void, matches that of a FLRW cosmology:

$$\begin{array}{c}\hfill {\sigma }^2(M,z)=\frac{1}{2{\pi }^2}{\displaystyle {\int }_0^{\infty }\mathrm{d}k{k}^2{P}_{\mathrm{m}}(k,z)W^2\left[k{R}_{\mathrm{th}}(M)\right],}\end{array}$$(11)

where M is the mass inside a sphere of radius R_th(M) = [3 M/(4 π ρ_m)]^1/3, ρ_m is the matter density, W(k R_th) is the Fourier transform of a top-hat filter at the scale R_th, and P_m(k, z) is the linear matter power spectrum of a ΛCDM universe at redshift z and wave number k.

Following Alonso et al. (2012), inside the ΛLTB void, σ(M, z) needs to be scaled by a factor, f, which is the ratio between density perturbations at two regimes:

$$\begin{array}{c}\hfill f(t,r)=\frac{{\delta }_{\alpha }(t,r)}{{\delta }_{\alpha }(t,r\to \infty )},\end{array}$$(12)

where δ_α is the linear density perturbation, with δ_α(t, r → ∞) matching the linear density contrast of the ΛCDM model. For small values of shear (see Eq. (14) below), it can be approximated by

$$\begin{array}{c}\hfill {\delta }_{\alpha }(t,r)=\overline{\delta }(t,r)[1+\alpha ϵ(t,r)],\end{array}$$(13)

which depends on the normalized shear parameter, ϵ. The latter is defined as

$$\begin{array}{c}\hfill ϵ\equiv \sqrt{\frac{2}{3}\frac{{\mathrm{\Sigma }}^2}{{\mathrm{\Theta }}^2}}=\frac{H_{\perp }-H_{\Vert }}{2H_{\perp }+H_{\Vert }},\end{array}$$(14)

giving the ratio between the square of the background shear, Σ, and the expansion parameter, Θ, where Σ² = Σ_ij Σ^ij is the square of the background shear, and Θ the expansion parameter of a congruence of comoving geodesics (see García-Bellido & Haugbølle 2009 for further details). Finally, in Eq. (13), α is a correction parameter that can be constrained from simulations or observations that we will discuss in Sect. 4, and $\overline{\delta }(t,r)$ is the local linear density contrast given by

$$\begin{array}{c}\hfill \overline{\delta }(t,r)\propto D({\mathrm{\Omega }}_{\mathrm{m},0}(r),R(t,r)/R_0(r)),\end{array}$$(15)

where

$$\begin{array}{c}\hfill D({\mathrm{\Omega }}_{\mathrm{m},0}(r),a)=a·{}_2{F}_1[\frac{1}{3};1;\frac{11}{6};\frac{{\mathrm{\Omega }}_{\mathrm{m},0}(r)-1}{{\mathrm{\Omega }}_{\mathrm{m},0}(r)}a],\end{array}$$(16)

with ₂F₁(a, b; c; z) the Gauss hypergeometric function and D the growth factor in a ΛCDM universe with matter parameter Ω_m, 0 (Bueno Belloso et al. 2011; Nesseris & Sapone 2015). We note that in our case we observe the redshift, therefore all the above equations were computed in z(r, t) using Eqs. (4) and (5) to then interpolate between the variables. In this framework, Alonso et al. (2012) modeled the mass function in LTB based on a FLRW mass function, with σ replaced by σ(M, z) f(t, r). Here, we adopted the same method and wrote the comoving number density of halos of mass, M, and redshift, z, as

$$\begin{array}{c}\hfill n(M,z)=-\frac{{\rho }_{\mathrm{m},0}}{M}\mathcal{F}[\sigma (M,z)f(t,r)]\frac{\mathrm{d}ln([\sigma (M,z)f(t,r){]}^{-1})}{\mathrm{d}M},\end{array}$$(17)

where ρ_m, 0 = ρ_m(z = 0), and the ℱ function taken from Despali et al. (2016) is

$$\begin{array}{c}\hfill \mathcal{F}(\sigma )=\nu A\sqrt{\frac{2a}{\pi }}[1+{\left(\frac{1}{a{\nu }^2}\right)}^p]exp[-\frac{a{\nu }^2}{2}].\end{array}$$(18)

ν ≡ δ_c/σ(M, z), where δ_c is the linear density threshold at the present time outside the void for nonlinear collapse and (a,  A,  p) = (0.3295,  0.7689,  0.2536). We rescaled following Despali et al. (2016) to match the survey mass definition.

We modified the monteLLTB code used in Camarena et al. (2021) to include a cluster counts module that implements the halo mass function recipe described above¹. The linear mater power spectrum was normalized using the σ_8, 0 parameter, defined as the root mean square of the variance of density fluctuations at R = 8 h⁻¹ Mpc in ΛCDM. In this work, we provide forecasts involving the following cosmological parameters: h = H₀/(100 km s⁻¹ M⁻¹ pc) (where H₀ := H_⊥(r → ∞) = H_∥(r → ∞) is the Hubble constant), Ω_m, 0 := Ω_m, 0(r → ∞), σ_8, 0 (hereafter denoted as σ₈), δ₀ := δ(r = 0, t₀), and z_B is the boundary redshift corresponding to r_B. The last two parameters are specific to the LTB void model, whereas the first three parameters of this list are in common with the ΛCDM cosmological model. Regarding the later, we adopted a flat ΛCDM fiducial model with h = 0.67, ${\mathrm{\Omega }}_{\mathrm{m},0}:={\mathrm{\Omega }}_{\mathrm{m},0}^{\mathrm{out}}=0.32$, ${\mathrm{\Omega }}_{K,0}:={\mathrm{\Omega }}_{K,0}^{\mathrm{out}}=0$ (i.e., K_B = 0), σ₈ = 0.83, and a spectral index of primordial scalar perturbations, n_s = 0.965.

3. Survey specifications and scaling relations modeling

3.1. Euclid and Rubin surveys’ richness-based cluster counts

3.1.1. Euclid-like cluster number counts forecast

To estimate the number of cluster counts, we used the forecast of Sartoris et al. (2016), in which the observable quantity is the observed mass, M_obs, with a sky coverage of Ω_tot = 15 000 deg². The estimated number counts for a redshift bin, z_l, and mass bin, M_obs, m, can be expressed as

$$\begin{array}{cc}\hfill N_{l,m}& ={\displaystyle {\int }_{{\mathrm{\Omega }}_{\mathrm{tot}}}\mathrm{d}\mathrm{\Omega }{\int }_{z_l}^{z_{l+1}}\mathrm{d}z\frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}\mathrm{\Omega }}}\hfill \\ \hfill & \times {\displaystyle {\int }_0^{+\infty }\mathrm{d}Mn(M,z)\frac{1}{2}[\text{erfc}(x_m)-\text{erfc}(x_{m+1})],}\hfill \end{array}$$(19)

where n(M, z) is the mass function, dΩ is the solid angle element in units of steradian, erfc is the complementary error function, and $\frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}\mathrm{\Omega }}$ is the derivative of the comoving volume with respect to the redshift and solid angle element, defined as

$$\begin{array}{c}\hfill \frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}\mathrm{\Omega }}=\frac{{(1+z)}^2{d}_{\mathrm{A}}{(z)}^{2}c}{H_{\Vert }(z)},\end{array}$$(20)

where d_A(z) = R[t(z),r(z)] is the angular diameter distance and c the speed of light. Finally, x_m = x (M_{obs, m, ℓ}) is defined in each mass bin, m, and redshift bin, ℓ², as

$$\begin{array}{c}\hfill x(M_{\mathrm{obs},m})=\frac{ln(M_{\mathrm{obs},m}/M_0)-ln(M/M_0)-ln(M_{\mathrm{bias}}/M_0)}{\sqrt{2{\sigma }_{M_{\mathrm{obs},m}}^2}},\end{array}$$(21)

where M₀ = 1 h⁻¹ M_⊙. We defined the mass bias and the variance, ${\sigma }_{M_{\text{obs},m}}^2$, to be the same as in S16:

$$\begin{array}{cc}& ln(M_{\mathrm{bias}}/M_0)=B_{M,0}+{\alpha }_{E}ln(1+z),\hfill \end{array}$$(22)

$$\begin{array}{cc}& {\sigma }_{M_{\mathrm{obs},m}}^2={\sigma }_{M_{\mathrm{obs}},0}^2-1+{(1+z)}^{2\beta }.\hfill \end{array}$$(23)

To estimate the number counts, we considered equally spaced redshift bins ranged in z ∈ [0.9, 2.0], with a width of Δz = 0.1. As for the limiting mass, we defined a mass selection function similar to the one in S16. In our analysis, we performed two types of tests. In the first one, we fixed the nuisance parameters above to the best fit of (B_M, 0,  α_E,  σ_Mobs, 0,  β) = (0.0 ± 0.05,  0.0 ± 0.05,  0.2 ± 0.07,  0.125 ± 0.00625), while in the second, we found the best-fit values for the same parameters, as is described in Sect. 4.

3.1.2. Rubin-like cluster number counts forecast

For Rubin, we followed The LSST Dark Energy Science Collaboration 30 (2018) and Murata et al. (2018), which modeled the galaxy cluster number counts as a function of their observed cluster richness. The probability of observing a certain richness, λ, given a mass, M, was modeled as

$$\begin{array}{c}\hfill P(ln\lambda |M)\mathrm{d}ln\lambda \equiv \frac{1}{\sqrt{2\pi }{\sigma }_{ln\lambda |M}}exp[-\frac{x^2(\lambda ,M)}{2{\sigma }_{\lambda |M}^2}]\mathrm{d}ln\lambda ,\end{array}$$(24)

with x defined by

$$\begin{array}{c}\hfill x(\lambda ,M)\equiv ln\lambda -[A+Bln(\frac{M}{M_{\mathrm{pivot}}})+Cln(1+z)],\end{array}$$(25)

with M_pivot = 3 × 10¹⁴ h⁻¹ M_⊙, and A, B, and C the dimensionless constants set below. Thus, the number counts at a given redshift and richness bin can be described by the following expression:

$$\begin{array}{cc}\hfill N& ={\displaystyle {\int }_{{\mathrm{\Omega }}_{\mathrm{tot}}}\mathrm{d}\mathrm{\Omega }{\int }_{z_{\min }}^{z_{\max }}\mathrm{d}z\frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}\mathrm{\Omega }}\int \mathrm{d}Mn(M,z){\int }_{{\lambda }_{\min }}^{{\lambda }_{\max }}\frac{\mathrm{d}\lambda }{\lambda }P(ln\lambda |M)}\hfill \\ \hfill & ={\displaystyle {\int }_{{\mathrm{\Omega }}_{\mathrm{tot}}}\mathrm{d}\mathrm{\Omega }{\int }_{z_{\min }}^{z_{\max }}\mathrm{d}z\frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}\mathrm{\Omega }}\int \mathrm{d}Mn(M,z)S(M|{\lambda }_{\min },{\lambda }_{\max }),}\hfill \end{array}$$(26)

where the subscripts min and max denote the minimum and maximum in each bin and S(M|λ_min, λ_max) is expressed as

$$\begin{array}{cc}\hfill S(M|{\lambda }_{\min },{\lambda }_{\max })& \equiv {\displaystyle {\int }_{{\lambda }_{\min }}^{{\lambda }_{\max }}\mathrm{d}ln\lambda P(ln\lambda |M)=\frac{1}{2}[\text{erf}(\frac{x({\lambda }_{\max },M)}{\sqrt{2{\sigma }_{\lambda |M}}})}\hfill \\ \hfill & -\text{erf}(\frac{x({\lambda }_{\min },M)}{\sqrt{2{\sigma }_{\lambda |M}}})].\hfill \end{array}$$(27)

erf is the error function, and the scatter, σ_lnΛ|M, is expressed as

$$\begin{array}{c}\hfill {\sigma }_{\lambda |M}={\sigma }_0+q_{M}ln(\frac{M}{M_{\mathrm{pivot}}})+q_{z}ln(1+z).\end{array}$$(28)

Here, we considered nine equally spaced redshift bins in the interval [0.0, 0.9], and 20 richness bins within the interval [20, 220] with a sky coverage of Ω_tot = 18 000 deg². In Sect. 4 we show results for the ΛLTB forecast, obtained by fixing the parameters (σ₀,  A,  B,  C,  q_M,  q_z) to the best fit of (0.456, 25 ± 5,  1 ± 0.8,  0 ± 1.2,  0 ± 0.05,  0 ± 0.2), as well as results when finding the best-fit values for these same parameters.

3.2. eROSITA-like survey X-ray-based cluster counts

Unlike the Euclid and Rubin surveys, eROSITA detects X-ray emission from clusters. Thus, our scaling relation is based on X-ray properties, as is described in Pillepich et al. (2012). Galaxy clusters were sorted in terms of the photon counts that will be detected by the eROSITA telescope. The counts were converted to dark matter haloes by taking into account the properties of the X-ray detector and the integration time of the observations. The X-ray cluster counts can be computed in a given redshift bin and observed temperature, T, as

$$\begin{array}{cc}\hfill N& ={\displaystyle {\int }_{{\mathrm{\Omega }}_{\mathrm{tot}}}\mathrm{d}\mathrm{\Omega }{\int }_{z_{\min }}^{z_{\max }}\mathrm{d}z\frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}\mathrm{\Omega }}\int \mathrm{d}Mn(M,z)}\hfill \\ \hfill & \times {\displaystyle {\int }_{T(M_{\min })}\mathrm{d}ln{T}_{\mathrm{X}}P(ln{T}_{\mathrm{X}}|M),}\hfill \end{array}$$(29)

where the subscripts min and max denote the minimum and maximum in each redshift bin and T_X = T/T₀, where T₀ = 1 keV and T(M_min) is the temperature bin lower limit obtained using Eq. (30). This would correspond to the mass taken from the following range, log₁₀(M_min/M₀)∈[13.0,  13.7,  14.0], at the three redshift bins with edges z ∈ [0.24,  0.3,  0.4,  0.52], and a sky coverage of Ω_tot = 27 500 deg² for this survey. P(lnT_X|M) is the probability distribution of the X-ray temperature, T_X, given the mass,

$$\begin{array}{c}\hfill P(ln{T}_{\mathrm{X}}|M)=\frac{1}{\sqrt{2\pi {\sigma }_{\mathrm{TM}}^2}}exp\{-\frac{{[ln{T}_{\mathrm{X}}-{\mu }_{\mathrm{T}}]}^2}{2{\sigma }_{\mathrm{TM}}^2}\},\end{array}$$(30)

μ_T is given by the scaling relation

$$\begin{array}{c}\hfill {\mu }_{\mathrm{T}}={\alpha }_{\mathrm{TM}}ln(\frac{M_{500}}{{\beta }_{\mathrm{TM}}})+{\beta }_{\mathrm{TM}}lnE(z)+ln(\frac{5\mathrm{keV}}{T_0}),\end{array}$$(31)

and E(z) = H(z)/H₀. The best-fit parameter values for this relation and the scatter, σ_TM, are (see Pillepich et al. 2012)

$$\begin{array}{cc}& {\alpha }_{\mathrm{TM}}=0.65\pm 0.03,\hfill \end{array}$$(32)

$$\begin{array}{cc}& {\beta }_{\mathrm{TM}}=(3.02\pm 0.11)\times {10}^{14}{h}^{-1}{M}_{\odot },\hfill \end{array}$$(33)

$$\begin{array}{cc}& {\sigma }_{\mathrm{TM}}=0.119.\hfill \end{array}$$(34)

3.3. DESI-like survey velocity mass determination-based cluster counts

In DESI, the total luminosity of a given cluster within a 1 Mpc radius (L_1 Mpc) is the mass proxy. This luminosity is related (and calibrated) to the mass, according to the scaling relation Eq. (35), below, as in Zou et al. (2021), where a fast clustering algorithm to identify the clusters was applied to the photometric redshift catalog. The total masses of the galaxy clusters were derived using a calibrated richness–mass relation that is based on the observations of X-ray emission, and the Sunyaev-Zeldovich effect:

$$\begin{array}{c}\hfill {\log }_{10}(L_{1\mathrm{Mpc}}/L^{\ast })=\frac{{\log }_{10}(M_{500}/M_0)}{a}-\frac{b}{a}{\log }_{10}(1+z)-\frac{c}{a},\end{array}$$(35)

where L^* is the characteristic luminosity and M₀ = 1 h⁻¹ M_⊙. The best-fit parameters of the above relation are (Zou et al. 2021)

$$\begin{array}{c}\hfill a=0.81\pm 0.02,b=0.50\pm 0.14,c=12.61\pm 0.04.\end{array}$$(36)

The above parameters were obtained from a photometric catalog containing sources that are used to obtain the spectroscopic measurements of their redshift. Therefore, we rescaled the error of the redshift parameter, b, to reflect the improvements in the error of z with respect to the expected accuracy of the spectroscopic measurements in DESI Collaboration (2016b) to reach more realistic constraints, and not heavily underestimate the power of DESI. Given the observables of this survey, we modeled the number counts in a given redshift and luminosity bin as

$$\begin{array}{c}\hfill N={\displaystyle {\int }_{{\mathrm{\Omega }}_{\mathrm{tot}}}\mathrm{d}\mathrm{\Omega }{\int }_{z_{\min }}^{z_{\max }}\mathrm{d}z\frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}\mathrm{\Omega }}\int \mathrm{d}Mn(M,z){\int }_{y_{\min }}\mathrm{d}yP(y|M),}\end{array}$$(37)

where the subscripts min and max denote the minimum and maximum in each redshift bin, where y = log₁₀(L_1 Mpc/L^*), and where P(y |M) has a probability that follows a Gaussian distribution similar to Eq. (30). For DESI, we took redshift bins with a thickness equal to 0.1 in the range z = [0.0,  1.0], a sky coverage of Ω_tot = 20 000 deg², and lower limits for y as presented below, for the numerical evaluation of the integrals in Eq. (37):

$$\begin{array}{cc}\hfill y_{\min }& ={\log }_{10}{(L_{1\mathrm{Mpc}}/L^{\ast })}_{\min }\in [1.65,2.13,2.35,2.50,2.61,2.69,\hfill \\ \hfill & 2.76,2.83,2.88,2.93].\hfill \end{array}$$(38)

Finally, for illustrative purposes, Fig. 1 shows integrated cluster counts in the redshift bins of the surveys used in this work, for a given ΛLTB model versus our fiducial ΛCDM cosmology. The cluster counts modifications within the ΛLTB model are, as expected, correlated with the void profile chosen in this work. In particular, we see that the Euclid sample change to the cluster counts is not noticeable. This is because its redshift bins fall outside the LTB void in the values adopted here. Though it will still serve to break degeneracies with the FRW cosmological parameters, Euclid is also expected to observe lower redshifts than our sample; however, here we already used those redshifts for the Rubin sample, since we considered both as our baseline unified survey. Other estimators were used in Camarena et al. (2023), as was detailed in the introduction, with similar constraints on the cluster counts, as we shall see later, but we still need to model other observables within ΛLTB, such as galaxy lensing or clustering, in order to be able to further benefit from Euclid or other similar high-redshift surveys for which the aforementioned probes constitute their main observational targets.

Fig. 1. Bars with heights representing cluster number counts in each of the redshift bins considered for the different surveys used in this work for ΛCDM. {δ0, zB}={0.0, 0.0} versus ΛLTB with {δ0, zB}={−0.5, 0.5}, while the other parameters were kept at their fiducial values.

In the following sections, we consider three forecast scenarios that differ by the degree of control one may have on the different systematic effects related to the nuisance parameters involving the scaling relations of the surveys described in this section. In Sect. 1 we list the cosmological and nuisance parameters of our modeling as well as the priors assumed in each scenario. These scenarios are:

• Optimistic settings: only the cosmological parameters are left free to vary with flat priors. In this case the nuisance parameters are fixed to their fiducial values.

• Realistic settings: the cosmological parameters are left free, while we consider Gaussian priors for the nuisance parameters of our probes.

• Pessimistic settings: the cosmological as well as the nuisance parameters of all surveys are left free to vary.

4. Results

We are mainly interested in assessing, as closely as possible, the constraints on the ΛLTB model from a Euclid-like survey when optimally combined with Rubin-like survey data. For that, our baseline approach was to combine cluster observations from Euclid and Rubin, knowing that Euclid is expected to perform better at high redshifts, for z > 0.9, whereas Rubin provides information for z < 0.9. We started with an optimistic case, assuming that we would have good control over the systematic effects, and thus we fixed their associated nuisance parameters to their fiducial values. We also fixed the correction parameter, α, in Eq. (13) to α = 2, the value found by simulations to give the best fit (Alonso et al. 2012). Furthermore, we also explored the case where α is left free to show that the impact on the results is minimal, as is seen in Fig. A.1. Finally, as in Camarena et al. (2023), we chose a flat prior for z_B ∈ [0, 0.5], since for δ₀ ≈ 0 the parameter space is highly degenerate with higher arbitrarily values of z_B. This baseline was then compared to the case in which we added data from eROSITA or DESI, these two being independent of each other. Thus, this led to three additional cases of combining survey information: Euclid + Rubin + DESI, Euclid + Rubin + eROSITA, and finally Euclid + Rubin + DESI + eROSITA. We also show the results of combining all the probes. These combinations are expected to offer tighter constraints than our baseline since we also fixed the nuisance parameters for these two probes in our first case.

We then repeated the same approach, but this time allowed the nuisance parameters of the baseline probes to vary within large priors in less optimistic configurations. We also combined Euclid and Rubin datasets with the eROSITA and DESI datasets in this scenario in order to break the degeneracies and balance the expected weakening of the constraints in the baseline case with Euclid and Rubin alone. This is expected to tighten the constraints even with wide priors on the primary probes’ nuisance parameters. Finally, we ran a further case, in which we combined all four probes at once, but this time with the most pessimistic settings, leaving all the nuisance parameters free to vary, in order to assess the loss of precision from this full auto-calibration approach.

4.1. Optimistic case assuming good control over the systematic effects

In this subsection, we start by showing, in Fig. 2, the marginalization of the parameters considered for the forecast using our baseline (Euclid + Rubin). We also show, in the same figure, the baseline with the addition of the secondary survey probes described in Sects. 3.2 and 3.3, all with fixed nuisance parameters, as this will be our optimistic case, in which we have control over the systematic effects. We also show in Table 2 the gain in accuracy for all the combinations and cases, which we discuss later with respect to the most pessimistic case in which all the parameters, including the nuisance ones, were left free to vary.

Fig. 2. ΛLTB 68% (darker color) and 95% (lighter color) confidence contours for Euclid and Rubin survey probes (orange), for when we include the DESI survey probe (green), the eROSITA survey probe (blue), and both (yellow) in the optimistic case (see Sect. 4.1).

We first observe that the cosmological parameters, Ω_m, 0, h, and σ₈, are constrained to a percent level while the density of the void, δ₀, does not prefer any value over the null one corresponding to its FLRW value. However, the redshift border, z_B, of the void is constrained to 0.4 and below at the 95% confidence level. Furthermore, as one could expect, when we add additional information from the secondary probes, we better constrain the FLRW-like parameters, especially Ω_m, 0, with the latter showing a gain of about 40%, as is seen in Table 2, while h and σ₈ increase by more than 10%. For the ΛLTB parameters, δ₀ still shows no preference within each optimistic case, while the z_B constraints increase by a factor of three when we combine all of the probes. We also notice that the baseline plus DESI offers similar constraints as the baseline plus eROSITA – due to the fact that eROSITA covers a larger sky area than DESI, while the latter spans a wider range and higher redshifts. Additionally, as was expected, our gain is further improved when we combine all of the probes. Finally, we also compare our results with the work of Camarena et al. (2023), as is shown in Table 3, where we only compare the optimistic case with all of the probes, so as to be on an equal footing with the aforementioned author’s work. Overall, we get better constraints on our parameters than Camarena et al. (2023), getting an improvement of about 50% for the FLRW parameters, 35% for z_B, and almost the same for δ₀.

4.2. Realistic case assuming Gaussian priors on the nuisance parameters

Here we discuss the case in which we consider nuisance parameters with Gaussian priors, with averages and variances taken from those given in Table 1. In this more realistic case, we consider ourselves to have less control over the systematic effects compared to the previous scenario, and therefore we allow some scattering in the nuisance parameters, thereby expecting looser constraints with respect to the optimistic case. On the other hand, when we combine or add more probes we should get tighter constraints on our parameters in any of the settings considered.

In Fig. 3 we show the contours within this scenario for the three cases: the baseline plus DESI, the baseline plus eROSITA, and the baseline plus DESI and eROSITA, while the gains with respect to the most pessimistic case are also in the same Table 2, along with those for the optimistic scenario. As predicted, the optimistic case shows better constraints when compared to its realistic counterpart, close to the same gain we observe when using cluster counts to constrain FLRW cosmological parameters. For the LTB parameters, we observe an overall improvement in the constraints, albeit smaller than that for the FLRW parameters. This is also supported by the general gain seen in the detailed Table 2, where the constraining power is degraded in most cases by more than a factor of two for each analogous probe, with the greatest decrease found when we compare the baseline only, reaching a factor of three for Ω_m, 0 and σ₈. The difference lessens when we choose to add the two additional surveys, DESI or eROSITA, to the baseline, and the lowest value is reached when we simultaneously add these two probes, staying around 20% lower on average. This can also be seen in Fig. 4, where we plot the two different analyses (the optimistic case is plotted in yellow and the realistic case in red) of the chains that include all the considered probes.

Fig. 3. ΛLTB and nuisance 68% (darker color) and 95% (lighter color) confidence contours for our realistic case, in which we considered Gaussian priors for the nuisance parameters in Table 1 in the calibration of the scaling relations. Here, we show the forecasts for the baseline (gray), and we also include secondary probes in the contours (Euclid + Rubin + DESI (black), Euclid + Rubin + eROSITA (pink), and finally Euclid + Rubin + DESi + eROSITA (red)).

4.3. Pessimistic constraints with auto-calibration of the nuisance parameters

Finally, we present the most ‘pessimistic’ case, in which we consider all the nuisance parameters of our probes to be free, by using non-informative flat priors for our Bayesian analysis. In Fig. 4, we show the cosmological and LTB model parameters forecast for this scenario, together with the analogous chains for the realistic and optimistic scenarios. The contours have widened with respect to the more optimistic ones, since most of the gains in Table 2 are positive; however, even when considering non-informative priors for the nuisance parameters, we are still able to put constrains on the LTB parameters at a higher level than in the baseline realistic case. This is seen in Table 2, which shows the power of a combination of probes to break and auto-calibrate the degeneracies between the ΛLTB and the nuisance parameters.

Fig. 4. ΛLTB 68% (darker color) and 95% (lighter color) confidence contours considering all the probes in this paper (Euclid + Rubin + eROSITA + DESI). Here, we compare the optimistic case (yellow) Sect. 4.1 with the realistic one (red) Sect. 4.2 and finally the pessimistic one (fuchsia) Sect. 4.3.

5. Conclusions

In this work we have forecast the forthcoming constraints for ΛLTB, an inhomogeneous model, using galaxy cluster abundance mock catalogs as our only probe, from ‘Stage-IV’-like experiments such as Euclid and Rubin, as our baseline. We have also studied the impact of taking into account observables from other experiments such as DESI and eROSITA.

We assumed different scenarios, reflecting our degree of confidence in the collected data: an optimistic case in which we fixed the nuisance parameters of the galaxy cluster scaling laws, a realistic case in which we allowed free scaling parameters within the Gaussian prior uncertainties, and finally a pessimistic case in which we assumed non-informative priors on the nuisance parameters.

We find that our combination of probes yields stronger constraints, especially for the non-LTB parameters. The constraints are degraded when the scaling laws are allowed more freedom in their parameter space. However, the combination of all probes is still able to break the degeneracies and auto-calibrate all of the nuisance parameters, even when they are left free to vary.

We also find that the ΛLTB model is still viable, even if the data turn out to be compatible with ΛCDM in the future. This was also the case in Camarena et al. (2023), which combined a plethora of geometrical probes as well as probes relating to the growth of structures, without including clusters. However, since our probe is further tackling the nonlinear scales for structure formation, while retaining some geometrical information from the sky coverage of the different surveys, we are able to get similar or better constraints for some parameters. Moreover, we highlight that we have separated the redshift domain of Euclid and Rubin (for 0 < z < 1 we consider Rubin only, while for Euclid we consider the redshift range of 1 < z < 2). While such an analysis was not performed in Camarena et al. (2023), our conservative approach is still competitive, highlighting the importance of the galaxy cluster counts probe in constraining cosmology within the ΛLTB model.

Acknowledgments

The authors would like to thank Stefano Camera for useful comments and style checks on the draft. Z.S. acknowledges funding from DFG project 456622116 and support from the IRAP and IN2P3 Lyon computing centers. This work was financed by Portuguese funds through FCT – Fundação para a Ciência e a Tecnologia in the framework of the following Research Projects: EXPL/FIS-AST/1368/2021 (A.d.S., A.C., J.P.M.), PTDC/FIS-AST/0054/2021 (N.N., A.d.S., A.C., J.P.M.), UIDB/04434/2020, UIDP/04434/2020, CERN/FIS-PAR/0037/2019 (N.N., A.d.S., A.C., J.P.M.), and 2022.04048.PTDC (C.J.M.). C.J.M. also acknowledges FCT and POCH/FSE (EC) support through Investigador FCT Contract 2021.01214.CEECIND/CP1658/CT0001, and A.C. acknowledges support form the FCT research grant 2020.06644.BD. J.G.B. and S.N. acknowledge support from the Research Project PID2021-123012NB-C43 [MICINN-FEDER], and the Centro de Excelencia Severo Ochoa Program CEX2020-001007-S. D.S. acknowledges financial support from the Fondecyt Regular project number 1200171. D.C. would also like to thank the Robert E. Young Origins of the Universe Chair fund for its generous support. This work was financed by FEDER – Fundo Europeu de Desenvolvimento Regional – funds through the COMPETE 2020 – Operational Programme for Competitiveness and Internationalisation (POCI) – and by Portuguese funds of FCT, under projects PTDC/FIS-AST/0054/2021 and UIDB/04434/2020 and UIDP/04434/2020.

References

Appendix A: Case with a free nonlinear correction prescription parameter

As is mentioned in the main text, here we show the optimistic case of our baseline (Euclid + Rubin) in Fig. A.1, where we compare the constraints obtained from having a fixed α = 2 (Eq. 13) to when we set this parameter as free to vary. We observe that the contours of these two cases are similar, as was expected from the simulations in Alonso et al. (2010) justifying fixing α for the other case’s analysis.

Fig. A.1. ΛLTB 68% (darker color) and 95% (lighter color) confidence contours for the baseline of our analysis: the Euclid and Rubin survey probes. In yellow we present the case in which we model α to a fixed value, whereas in red we display the case where our modeling of the problem includes a free α with a wide prior, ∈[0, 4].

All Tables

All Figures

	Fig. 1. Bars with heights representing cluster number counts in each of the redshift bins considered for the different surveys used in this work for ΛCDM. {δ0, zB}={0.0, 0.0} versus ΛLTB with {δ0, zB}={−0.5, 0.5}, while the other parameters were kept at their fiducial values.
In the text

	Fig. 2. ΛLTB 68% (darker color) and 95% (lighter color) confidence contours for Euclid and Rubin survey probes (orange), for when we include the DESI survey probe (green), the eROSITA survey probe (blue), and both (yellow) in the optimistic case (see Sect. 4.1).
In the text

Fig. 3. ΛLTB and nuisance 68% (darker color) and 95% (lighter color) confidence contours for our realistic case, in which we considered Gaussian priors for the nuisance parameters in Table 1 in the calibration of the scaling relations. Here, we show the forecasts for the baseline (gray), and we also include secondary probes in the contours (Euclid + Rubin + DESI (black), Euclid + Rubin + eROSITA (pink), and finally Euclid + Rubin + DESi + eROSITA (red)).

In the text

	Fig. 4. ΛLTB 68% (darker color) and 95% (lighter color) confidence contours considering all the probes in this paper (Euclid + Rubin + eROSITA + DESI). Here, we compare the optimistic case (yellow) Sect. 4.1 with the realistic one (red) Sect. 4.2 and finally the pessimistic one (fuchsia) Sect. 4.3.
In the text

	Fig. A.1. ΛLTB 68% (darker color) and 95% (lighter color) confidence contours for the baseline of our analysis: the Euclid and Rubin survey probes. In yellow we present the case in which we model α to a fixed value, whereas in red we display the case where our modeling of the problem includes a free α with a wide prior, ∈[0, 4].
In the text