\section{Introduction \label{sec:intro}}

Modern cosmology has been built on a combination of theoretical considerations, observations, and non-empirical assessments~\citep{2020coce.book.....P}. Illustrative of the last two is the cosmological principle, which combines the observation that there seems to be no preferred direction in the sky (local isotropy), with the non-empirical assessment that we do not occupy a special place in the cosmos (Copernican principle). Together, they imply that the Universe on the largest scales must be spatially homogeneous -- no special place -- and isotropic -- no special direction, -- thereby motivating its description by the  Friedmann-Lema\^{\i}tre-Robertson-Walker (FLRW) geometry.

In addition to symmetries, one must specify dynamical prescriptions in order model the history of cosmic expansion within the FLRW model. This aspect has mostly come from theoretical physics and observation. In the standard lore, the coupling of matter and space-time geometry is described by the hitherto-unchallenged general relativity; at late times, which is the focus of the present article, its content is largely dominated by pressureless matter; the recent acceleration of cosmic expansion is attributed to Einstein's cosmological constant~$\Lambda$, which may also be interpreted as a negative-pressure fluid with equation of state~$w=p/\rho=-1$; and the curvature of homogeneity hyper-surfaces vanishes.

The set of all the aforementioned features is collectively referred to as the standard cosmological model, or $\Lambda$CDM\xspace\,-- which also includes initial conditions set by the inflationary paradigm. To date, $\Lambda$CDM\xspace is practically compatible with all cosmological observations~\citep{Aghanim:2018eyx}, and hence so are its underlying assumptions. Albeit intensely discussed in the community, the so-called $H_0$ tension~\citep{DiValentino:2021izs}, the $\sigma_8$ tension~\citep{Sakr:2018new}, and others~\citep{Perivolaropoulos:2021jda}, have not proved convincing enough to claim the discovery of new physics yet.

Next-generation surveys, such as the \Euclid mission, will represent a leap forward in testing $\Lambda$CDM\xspace and its assumptions, for the amount and quality of their data will be unprecedented. The goal of this article is to assess the performance of a number of \emph{null tests} for $\Lambda$CDM\xspace in the context of \Euclid. These tests are designed to detect any deviation from the history of cosmic expansion, as predicted by the concordance model. They are agnostic in the sense that they do not require any specific alternative to be tested against $\Lambda$CDM\xspace; they merely address the question: Has \Euclid detected new physics? Nevertheless, should one or several tests reviewed here reject $\Lambda$CDM\xspace, their variety would help identifying which underlying assumption of $\Lambda$CDM\xspace would be violated: Is the Universe inhomogeneous or anisotropic on large scales? Is dark energy different from $\Lambda$? Do dark matter and dark energy interact? Is the Universe spatially curved? Etc.




Deviations from spatial homogeneity have been previously constrained with a plethora of different probes~-- see \cite{Clarkson:2012bg} for a review. For example, by using standard candles \citep{Chiang:2017yrq}, the fossil record of galaxies \citep{Heavens:2011mr}, the BOSS DR12 quasar sample \citep{Laurent:2016eqo} and with both spectroscopic \citep{Scrimgeour:2012wt} and  photometric redshift surveys  \citep{Alonso:2014xca}, as well as the kinetic Sunyaev-Zeldovich (kSZ) effect on the Cosmic Microwave Background (CMB)  \citep{Reichardt:2020jrr} and measurements of our peculiar velocity \citep{Nadolny:2021hti}. On the theory side, several tests have been proposed, for example by considering the propagation of light rays in an inhomogeneous universe using numerical relativity \citep{Giblin:2016mjp}, using distances to constrain directly the spatial curvature of the Universe \citep{Clarkson:2007pz,Clarkson:2012bg,Valkenburg:2012td}; using inhomogeneous Lema\^{\i}tre-Tolman-Bondi (LTB) models  \citep{GarciaBellido:2008nz,February:2009pv,Redlich:2014gga}; consistency tests based on dynamical probes, such as the growth rate data \citep{Nesseris:2014mfa,Nesseris:2014qca}; and through a linear model formalism \citep{Marra:2017pst}. 

Future surveys are mainly expected to provide stringent constraints on null tests of $\Lambda$CDM via very accurate and precise distance measurements of baryon acoustic oscillations (BAO), Type Ia supernovae (SNe) and lensing shear correlations. In this work we will explicitly focus on \Euclid, which is an M-class space mission of the European Space Agency (ESA) due for launch in 2022 \citep{Racca:2016qpi}. \Euclid will carry on board the near-infrared spectrophotometric instrument \citep{NISP_paper} and the visible imager \citep{VIS_paper}. Jointly these instruments will perform a spectroscopic and a photometric galaxy survey over 15\,000 deg$^2$ of sky, aiming to map the geometry of the Universe and measure the growth of structures up to $z\sim 2$~\citep{Laureijs:2011gra}.

The main cosmological probes of \Euclid will be: galaxy clustering from the spectroscopic survey, and galaxy clustering and weak lensing from the photometric survey. Spectroscopic galaxy surveys offer much higher radial precision, albeit for fewer objects, while photometric surveys target a larger number of galaxies, but with also larger redshift uncertainties. \Euclid in particular, given its high spectroscopic accuracy, will have very precise galaxy clustering constraints that also include the radial (i.e. along the line of sight) dimension. In our work we will create mock BAO data using \Euclid's spectroscopic survey specifications, following the Fisher matrix approach as in \citet[][]{IST:paper1}, hereafter EC20.

We will also highlight some of the synergies between \Euclid and other large scale structure surveys, namely that of the Dark Energy Spectroscopic Instrument (DESI) \citep{DESI2016} and the Legacy Survey of Space and Time (LSST) survey, performed at the Vera C. Rubin Observatory \citep{Abell:2009aa}. Both will probe the expansion history of the Universe and its large-scale structure (LSS), while being complementary to \Euclid redshift ranges, thus greatly extending the redshift range of our constraints.

In \citet{ReviewDoc}, forecast constraints on deviations from spatial homogeneity and the Copernican principle were presented by performing a joint analysis between the \Euclid galaxy survey \citep{Laureijs:2011gra} and a stage IV supernova mission \citep[][assuming SNAP as a concrete example]{DEtaskforce}. In this work we update these results by relying on more recent \Euclid specifications, \citepalias[see][]{IST:paper1}, we  explore synergies with other surveys and we extend the variety of tests used. Moreover, we also implement machine learning and other model-independent approaches to reconstruct these tests, so as to avoid theoretical biases.

The outline of our paper is as follows: in \Cref{sec:theory} we review the theoretical background behind the fundamental assumptions of the standard cosmological model and various ways in which these can be broken, then in \Cref{sec:tests} we summarize the null tests used in our analysis. In \Cref{sec:cur_data} we describe the analysis methods followed in this work to test the fundamental assumptions of the standard cosmological model, using the currently available data and a machine learning approach, while in \Cref{sec:fore_data} we discuss how we produce mock data based on three different cosmologies: the vanilla 
$\Lambda$CDM\xspace model, used as our null hypothesis, an evolving dark energy equation of state $w(a)$ based on the Chevallier-Polarski-Linder (CPL) parametrization \citep{Chevallier:2000qy,Linder:2002et}, in order to examine the response of our null tests to different expansion histories, and a fiducial cosmology based on the LTB metric. The results of our analysis are then  discussed in \Cref{sec:currentresults} for the currently available data and in \Cref{sec:fore_res} for the mock data and finally, we draw our conclusions in \Cref{sec:conclusions}.

\section{Theoretical background \label{sec:theory}}

The consistency tests to be considered in this work are motivated by various alternatives to the $\Lambda$CDM model that have been proposed in the literature. We consider deviations from this baseline model ranging from simple modifications of the expansion history, to models whose common feature is the relaxation of either the assumption of isotropy or that of spatial homogeneity of the Universe, or both. 

As an example of the first category concerning modifications of the expansion history, we focus in this work on the CPL parametrization of dark energy (DE), where the equation of state parameter $w$ for this component deviates from the constant $w=-1$ value that it takes in $\Lambda$CDM and is allowed to vary with redshift. In this parametrization, the DE equation of state parameter is given by \citep{Chevallier:2000qy,Linder:2002et}
\begin{equation}\label{eq:CPLpar}
w(z) = w_0 + w_a \, \frac{z}{1+z}\, ,
\end{equation}
with $w_0$ and $w_a$ free parameters.

Concerning the second category, which contains models that relax the assumption of isotropy or spatial homogeneity of the Universe, for our purposes in this work the possible alternatives may be grouped under five phenomenological approaches to modify the aforementioned hypotheses.

The first is a simple deviation from the Copernican principle, whose effect is non-perturbative. Its origin could be the presence of a big attractor affecting the observer's position. A typical observer will measure a large  dipole or other multipoles in the cosmic microwave background (CMB) spectrum due to the change in its peculiar velocity with respect to the CMB frame, but also due to changes in the gradients of the expansion rate and the divergence of the shear constructed from the two now-measured directional Hubble parameters.

Next we mention the class of LTB models, proposed by \citet{Tomita:1999qn,Celerier:1999hp,Alnes:2005rw} as an alternative to DE \citep[see also][]{Moffat:1994qy,Mustapha:1998jb}.  This is an inhomogeneous universe model with spherical symmetry around a point, thus, exhibiting local isotropy at that point only.
The possibility of a central gigaparsec-scale void  could explain SNe observations without requiring either an exotic type of repulsive content or any modification of gravity, see  \citet[]{GarciaBellido:2008nz}.
The LTB metric is given by 
\begin{equation}
\mathrm{d}s^{2}=-c^2\mathrm{d}t^{2}+\frac{\left[R'(t,r)\right]^{2}}
{1-K(r)}\mathrm{d}r^{2}+R^{2}(t,r)\mathrm{d}\Omega^{2}\,,
\label{eq:LTB}
\end{equation}
where $K(r)$ is an arbitrary function, a prime denotes a partial spatial derivative with respect to the radial dimension $r$.
If $R=a(t)\,r$ and $K(r)=k_0 r^{2}$, we recover the familiar FLRW metric.
For a review of these models see \citet{Marra:2011ct} and references therein.

It is worth stressing that void models as alternatives to dark energy have been ruled out. Models without decaying modes produce a too large kinetic Sunyaev-Zeldovich (kSZ) signal~\citep{GarciaBellido:2008gd,Zhang:2010fa,Moss:2011ze,Bull:2011wi}. On the other hand, models with sizeable decaying modes (which could possibly have a small kSZ signal) are ruled out because of $y$-distortions \citep{Zibin:2011ma}.  The only possibility to save such void models would be by choosing specific (fine-tuned) initial conditions, which would be contrived, or inhomogeneous radiation and baryon fraction~\citep{Clarkson:2010ej}. Consequently, in the following we will consider the so-called $\Lambda$LTB model, that is the LTB model with a cosmological constant. The latter will therefore provide all or most of the dark energy, while we can phenomenologically constrain a small void contribution.

A third class of models that instead abandons the isotropy hypothesis is the Bianchi models, see \citet[]{Ellis:1968vb,Ellis:1998ct}. These are spatially homogeneous models that are endowed with a group of isometries transitively acting on the spatial hypersurfaces. One minimal modification is obtained in the type-I Bianchi models, which are characterized by a vanishing spatial curvature, and contain as a particular case the flat FLRW models. Their line element can be written as
\begin{equation}
\mathrm{d}s^{2} = - c^2\mathrm{d}t^2 + a_x^2(t) \mathrm{d}x^2 + a_y^2(t) \mathrm{d}y^2 + a_z^2(t) \mathrm{d}z^2 \,,
\label{metric}
\end{equation}
where we see that at a given time $t$, the metric is spatial homogeneous (i.e. does not depend on the spatial coordinates). In this metric, a sphere of comoving test particles experiences a change of volume with time and, unlike the FLRW case, distorts in shape into an ellipsoid. The former is characterized by the average expansion
\begin{equation}
    3H = H_x+H_y+H_z\; ,
\end{equation}
where $H_x$, $H_y$ and $H_z$ is the Hubble expansion rate along the $x$, $y$ and $z$ dimensions respectively, while the rate of shape deformation is described by the shear $\sigma$. This can be completed by specifying the components of the spatial metric $h_{\mu\nu}$ at a single point and a unit timelike normal vector to the hypersurface. Then projecting the equations onto the $\vec{n}$ direction gives a generalized Friedmann equation,
\begin{equation}
  3H^2  = 8 \pi G \rho + \sigma^2 -  {}^{(3)}R/2 + \Lambda \, c^2  \,, \label{eq:friedmann}\end{equation}
where $G$ is the gravitational constant, $\rho$ is the
energy density of the matter content, $\sigma^2 \equiv \sigma_{\mu \nu} \sigma^{\mu\nu}/2$ is the shear scalar, $\sigma_{\mu\nu}=\theta_{\mu\nu}-(\theta/3)\,h_{\mu\nu}$ is the shear tensor, while $\theta_{\mu\nu}=h^\alpha{}_\mu\,h^\beta{}_\nu\,X_{(\alpha;\;\beta)}$ is the expansion tensor, $\theta$ is its trace and $h_{\mu\nu}=g_{\mu\nu}+X_\mu\,X_\nu$ for a timelike unit vector $X^\mu$. Finally, ${}^{(3)}R$ is the Ricci scalar of the 3-curvature of the spatial hypersurfaces, and $\Lambda$ is the cosmological constant.

Another phenomenological approach stems from the idea that inhomogeneities in our Universe backreact and affect the average dynamics within our causal horizon such that the observed acceleration might be attributed to their influence. So, if we would like to compute directly the impact of the inhomogeneities, without requiring a highly symmetric solution of Einstein’s equations, one ansatz used in what is usually called ``backreaction theory'' is to construct average quantities that follow equations similar to those of the traditional FLRW model. The averaged scale factor $a_{\cal D}$, with extra contributions, can then be written
\begin{align}
\label{averagedhamilton}
3 \frac{{\dot a}^2_{\cal D}}{a^2_{\cal D}} - 8\pi G \langle{\varrho}\rangle-\Lambda c^2 &= - \frac{\langle{R}_{\cal D}\rangle+{Q}_{\cal D} }{2} \,, \\
3\frac{{\ddot a}_{\cal D}}{a_{\cal D}} + 4\pi G \langle{\varrho}\rangle -\Lambda c^2&= {Q}_{\cal D}\,, \\
\dot{\langle{\varrho}\rangle}_{\cal D} + 3\frac{{\dot a}_{\cal D}}{a_{\cal D}} \langle{\varrho}_{\cal D}\rangle&=0\,,
\end{align}
where the scale factor has been replaced by an average one. The new kinematic back-reaction source term 
\begin{equation}
\label{backreactionterm} 
{Q}_{\cal D} =
\frac{2}{3}\left\langle{\Big(\theta - \langle{\theta}\rangle\Big)^2 }\right\rangle - 
2\left\langle{\sigma^2}\right\rangle \,,
\end{equation}
arises as a result of the expansion rate $\theta$  and shear $\sigma$ fluctuations. For reviews of the underlying mathematical formalism see \citet{Buchert:2007ik,Kolb:2009rp, Andersson:2011za} and references therein. There is, however, a lack of consensus on whether backreaction is a correct framework and that inhomogeneities can generate acceleration, see \citet{Clarkson:2011zq,Green:2014aga, Buchert:2015iva, Kaiser:2017hqn}.

Finally, we briefly mention the Swiss-cheese model universe \citep[]{Einstein:1945id,Marra:2007pm}, where inhomogeneous LTB patches (the holes) are embedded in a flat FLRW background (the cheese), in a way that the dynamics of the holes is scale-independent. In other words, small holes will evolve in the same way as big holes, and isotropy is therefore preserved \citep{Marra:2007pm,Biswas:2007gi}. One simple approach is to assume no shell crossing and to define the initial conditions for each shell at the same time, with FLRW density matched at the boundary of a hole.\\

In what follows, we describe the phenomenological formulations which serve as null tests to detect deviations from the hypothesis of the homogeneity and isotropy of the Universe.
To do that we have to assume a fiducial cosmology to compare with. Even if that could be enough for our `trigger' tests to suggest the need to revise some of the aforementioned assumptions, we decided to also add a more subtle treatment where the fiducial and the null tests are considered in the framework of some of the above models. Therefore, we consider two well studied and widely used models which have been proposed as extensions and alternatives to $\Lambda$CDM: the CPL model described above and the $\Lambda$LTB model we detail later on, which has recently been shown to remain as a viable alternative to $\Lambda$CDM \citep{Camarena:2021mjr}.

 
\section{The consistency tests \label{sec:tests}}
We now present the set of consistency tests of the $\Lambda$CDM\xspace model that we are going to reconstruct in what follows, using mock data from \Euclid and other contemporary surveys. We explicitly use a variety of null tests of the $\Lambda$CDM\xspace model and of the FLRW metric itself, as each one is sensitive to different probes in the redshift range covered by \Euclid. For example, while the Om statistic directly probes the Hubble expansion, the global shear test probes the spatial homogeneity of the Universe and so on. 

\subsection{The Om statistic\label{sec:omtest}}
In a flat $\Lambda$CDM\xspace Universe, the Hubble parameter at late times, when we can safely neglect radiation, is given by 
\begin{equation}
H^2(z)=H^2_0\left[\Omega_{\rm m,0}\,(1+z)^3+1-\Omega_{\rm m,0}\right]\,,\label{eq:Hlcdm}
\end{equation} 
where $H_0$ is the Hubble constant, $\Omega_{\rm m,0}\equiv \rho_\mathrm{m}/\rho_\mathrm{c}$ is the fractional matter density parameter and $\rho_\mathrm{c}$ is the critical density for which the Universe is flat. Then, solving for $\Omega_{\rm m,0}$ we can create the so-called
$\textrm{Om}_\mathrm{H}(z)$ quantity, defined as 
\begin{align}
\textrm{Om}_\mathrm{H}(z) &\equiv \frac{h^2(z)-1}{(1+z)^3-1}
\longrightarrow \Omega_{\rm m,0} \quad \textrm{implies~ $\Lambda$CDM}\,,
\label{eq:omh1} 
\end{align}
which has to be constant and equal to the matter energy density $\Omega_{m,0}$ only when $\Lambda$CDM\xspace is the true model describing the evolution of the Universe and where we have defined the dimensionless Hubble parameter $h(z)\equiv H(z)/H_0$, see \citet[]{Sahni:2008xx,Zunckel:2008ti}. 

Note however that in the zero redshift limit this quantity is ill-defined, which can limit the  precision and accuracy of this test. This might also introduce biases if the numerator and denominator uncertainties do not behave in the same way. In general though, any deviations of \Cref{eq:omh1} from a constant and its $\Lambda$CDM\xspace value given by \Cref{eq:Hlcdm}, imply that $\Lambda$CDM\xspace does not hold, irrespective of the $\Omega_{\rm m,0}$ value. Reconstructions of the $\textrm{Om}_\mathrm{H}(z)$ test with earlier data were performed in \citet[]{Nesseris:2010ep}.

\subsection{Extensions of the Om statistic with curvature}
We also consider the extended $\textrm{Om}_\mathrm{H}(z)$ statistic by \citet{Seikel:2012cs} when curvature is present. In this case we can solve simultaneously the $\Lambda$CDM\xspace Friedmann equations for both the matter and the curvature parameters and doing so, we find the two tests
\begin{align}
\mathcal{O}_{\rm m}(z) &\equiv  2\, \frac{ (1+z)\,[1-h^2(z)]+z\,(2+z)\,h(z)\,h'(z)}{z^2\,(1+z)\,(3+z)}\,, \label{Om-hz}\\
\mathcal{O}_{\rm K}(z) &\equiv \frac{ 3\,(1+z)^2\,[h^2(z) -1]- 2z\, (3+3z+z^2)\,h(z)\,h'(z)}{z^2\,(1+z)\,(3+z)}\,, \label{OK-hz}
\end{align} 
where a prime is a derivative with respect to the redshift $z$ and again we have 
\begin{align}
\mathcal{O}_{\rm m}(z) &\longrightarrow  \Omega_{\rm m,0}
\quad\text{implies $\Lambda$CDM}, \\
\mathcal{O}_{\rm K}(z) &\longrightarrow  \Omega_{\rm k,0}
\quad \text{implies $\Lambda$CDM}. 
\end{align}
The main advantage of these null tests compared to the  $\textrm{Om}_\mathrm{H}(z)$ is that we do not need to make any assumptions on the spatial curvature of the Universe. Other possibilities that also include information from the distances were explored in \citet{Yahya:2013xma}.

\subsection{The $r_0$ test for interactions in the dark sector}
Recently \citet[]{vonMarttens:2018bvz} (see also \citealt{vonMarttens:2020apn}) proposed a null test which is based on the ratio between cold dark matter (CDM) and DE energy densities, that is $r(z) =\rho_\textrm{CDM}(z)/\rho_\textrm{DE}(z)$, which for the $\Lambda$CDM\xspace model  is equal to $r(z)=r_0\,(1 +z)^3$, with $r_0=\Omega_{\rm c,0}/\Omega_{\rm DE,0}$. 
Then the Friedmann equation for $\Lambda$CDM\xspace can be rewritten as
\begin{equation} 
h^2(z) =\Omega_{\rm d,0}\,\frac{1+r_0\,(1+z)^3}{1+r_0}+\Omega_{\rm b,0}(1+z)^3\,, \label{eq:Hz2}
\end{equation} 
where $\Omega_{\rm d,0}=\Omega_{\rm c,0}+\Omega_{\rm DE,0}=1-\Omega_{\rm b,0}$. We assume flatness and we ignore radiation as we only consider low redshift data.

Solving \Cref{eq:Hz2} for $r_0$, we can define a null test for models with interactions in the dark sector as
\begin{equation}
r_0=\frac{1-\Omega_{\rm b,0}+\Omega_{\rm b,0}(1+z)^3-h^2(z)}{h^2(z)-(1+z)^3}, \label{eq:r0null}
\end{equation} 
which for $\Lambda$CDM\xspace has to be constant at all redshifts. Note however, that again this test is ill-defined at low redshifts and also requires an external prior on the baryon density $\Omega_{\rm b,0}$. The natural choice would be to use the value obtained either from CMB data, e.g. \textit{Planck}, or from Big Bang nucleosynthesis (BBN), but currently the two do not agree \citep{Pitrou:2020etk}. Thus, in what follows we will use the \textit{Planck} 2018 best-fit value.

Also note that the tests of \Cref{eq:omh1} and \Cref{eq:r0null} are only degenerate in the case of $\Lambda$CDM, where one may show that indeed
\begin{equation} 
r_0(z)=\frac{\textrm{Om}_\mathrm{H}(z)-\Omega_\mathrm{b,0}}{1-\textrm{Om}_\mathrm{H}(z)}\longrightarrow \text{constant}, 
\end{equation} 
which implies $\Lambda$CDM and that if either of \Cref{eq:omh1} or \Cref{eq:r0null} is a constant, then so is the other. Note, however, that $r_0$ does depend on the value of $\Omega_\mathrm{b,0}$ and in general, the relation between $r_0(z)$ and $\textrm{Om}_\mathrm{H}(z)$ might be more complicated for non-$\Lambda$CDM\xspace models. Furthermore, as discussed in \citet[][]{vonMarttens:2018iav}, a nontrivial $r_0(z)$ can be mapped directly into the DM-DE coupling function. 

Overall, the $r_0(z)$ test probes for interactions in the dark sector, which goes somewhat further than just testing for deviations from the $\Lambda$CDM\xspace model. In our work however, we will not consider fiducial cosmologies with DM-DE interactions.

\subsection{Global Shear}
We can also explore another test of the Copernican principle in the form of the global shear, see \citet{GarciaBellido:2008yq}, given by
\begin{align}
\Sigma(z)
&= \frac{\sqrt{1-kr^2(z)}-h(z)d'(z)}
{3h(z)\,\frac{H_0}{c}\,D_{\rm A}(z)+2\sqrt{1-kr^2(z)}-2h(z)d'(z)}
\nonumber \\[1mm]
&\simeq
\frac{1-h(z)d'(z)}{3h(z)\,\frac{H_0}{c}\,D_{\rm A}(z)+2-2h(z)d'(z)}
\,,
\label{eq:GlobalShear}
\end{align}
where the primes denote a derivative with respect to the redshift $z$ and we have defined the {\it dimensionless} comoving distance:
\begin{equation} 
d(z)=(1+z)\;\frac{H_0}{c}\;D_{\rm A}(z)=\frac{1}{\sqrt{-\Omega_{\rm k,0}}}\sin{\left[
\sqrt{-\Omega_{\rm k,0}}\int_0^z{\frac{\mathrm{d}x}{h(x)}}\right]}\,.\label{eq:distances}
\end{equation}
In the case of a FLRW metric, the global shear given by \Cref{eq:GlobalShear} should always be zero regardless of the curvature \citep[see][]{GarciaBellido:2008yq}.

\subsection{The distance null tests}
Applying the Lagrange inversion theorem to the dimensionless luminosity and angular diameter distances, $d_{\rm L}(z,\Omega_{\rm m,0})\equiv c^{-1}H_0\,D_{\rm L}(z,\Omega_{\rm m,0})$ and  $d_{\rm A}(z,\Omega_{\rm m,0})\equiv c^{-1}H_0\,D_{\rm A}(z,\Omega_{\rm m,0})$ respectively, we can solve for $\Omega_{\rm m,0}$ similarly to the $\textrm{Om}_\mathrm{H}(z)$ test, see \citet[]{Arjona:2019fwb}.

Restricting ourselves to late times, when DE dominates over the other components, we may safely neglect radiation and neutrinos. Thus, the analytical expression of the dimensional luminosity distance for the $\Lambda$CDM\xspace model, assuming a flat Universe, but neglecting radiation and neutrinos, is given by
\begin{multline}
D_{\rm L}(z,\Omega_{\rm m,0})
= c(1+z) \int_0^z\frac{{\rm d} x}{H(x)}
= \frac{c}{H_0}\frac{2(1+z)}{\sqrt{\Omega_{\rm m,0}}}
\\[2mm]
\times
\left\{_2F_1\left(\frac{1}{2},\frac{1}{6};\frac{7}{6};\frac{\Omega_{\rm m,0}-1}{\Omega_{\rm m,0}}\right)-
\frac{_2F_1\left[\frac{1}{2},\frac{1}{6};\frac{7}{6};\frac{\Omega_{\rm m,0}-1}{\Omega_{\rm m,0}(1+z)^{3}}\right]}{\sqrt{1+z}}\right\}
,
\label{eq:dl}
\end{multline}
where $_2F_1\left(a,b;c;x\right)$ is a hypergeometric function.

To invert the previous equation we first do a series expansion on \Cref{eq:dl} around $\Omega_{\rm m,0}=1$ and we keep the first $12$ terms in order to obtain a reliable unbiased estimation and avoid theoretical systematic uncertainties due to having to truncate the series expansion. Alternatively, one may directly Taylor expand the integrand of \Cref{eq:dl} and then perform the integration term by term, so the two approaches are equivalent.

Finally, we apply the Lagrange inversion theorem to actually invert the series and to write the matter density $\Omega_{\rm m,0}$ as a function of the dimensionless luminosity distance $d_{\rm L}$, i.e $\textrm{Om}_\textrm{dL}=\textrm{Om}_\textrm{dL}(z,d_{\rm L})$. For example, the first two terms of the  expansion are
\begin{equation}
\textrm{Om}_\textrm{dL}(a,d_{\rm L})=1-\frac{7a\left(d_{\rm L}-\frac{2-2\sqrt{a}}{a}\right)}{6+\sqrt{a}\left(a^3-7\right)}+O(d^2_{\rm L})\,, \label{eq:OmdLnull}
\end{equation}
where the scale factor $a$ is related to the redshift $z$ as $a=1/(1+z)$ and as mentioned, for the actual calculations we keep the first $12$ terms in the expansion. 

This null test has the main advantage that it does not require taking derivatives of the data as we use the luminosity distance directly. Similarly, assuming the distance duality relation (DDR) holds, as we do in this work, then the angular diameter distance can be calculated via
\begin{equation}
    D_{\rm A}(z)=\frac{D_{\rm L}(z)}{(1+z)^2}\,,
\end{equation}
and assuming the number conservation of photons, we can follow the same approach as for the luminosity distance and find through the Langrange inversion theorem an expression for the matter density  $\Omega_{\rm m,0}$ as a function of the dimensionless angular diameter distance $d_{\rm A}(z)$, i.e  $\textrm{Om}_\textrm{dA}=\textrm{Om}_\textrm{dA}(z,d_{\rm A})$. The first two terms of the expansion in this case are:
\begin{equation}
\textrm{Om}_\textrm{dA}(a,d_{\rm A})=1-\frac{7\left[d_{\rm A}+2a\left(\sqrt{a}-1\right)\right]}{a\left[6+\sqrt{a}\left(a^3-7\right)\right]}+O(d^2_{\rm A})\,. \label{eq:OmdAnull}
\end{equation}
As mentioned earlier, also in this case we keep the first $12$ terms of the expansion for the actual calculations.

\subsection{The curvature test}
A direct test of the Copernican principle is the curvature test of \citet{Clarkson:2007pz}, which considers the distance-redshift relation in a FLRW model of any curvature. Solving \Cref{eq:distances} for the curvature parameter we find
\begin{equation}
\label{eq:OK}
\Omega_{\rm k}(z)=\frac{\left[h(z)d'(z)\right]^2-1}{[d(z)]^2}\longrightarrow\text{constant}\,.
\end{equation}
One way to interpret \Cref{eq:OK} is that it provides the means to measure the current value of the curvature parameter today by using Hubble rate and distance data. In the context of the FLRW metric this parameter will be constant regardless of the dark energy model. Equivalently, we can also write the aforementioned condition as 
\begin{multline}
   \label{C(z)}
\mathcal{C}(z)
= 1+h^2(z)\left\{d(z)\,d''(z)-\left[d'(z)\right]^2\right\}\\
+ h(z)h'(z)d(z)d'(z)
\longrightarrow 0\,. 
\end{multline}
In general for spacetimes other than FLRW metrics (e.g. in LTB models) we have that $\mathcal{C}(z)\neq0$, or $\Omega_{\rm k}(z)\neq$\,const.

\begin{figure}[!t]
\centering
\includegraphics[width = 0.48\textwidth]{figures/flowchart_om_2.pdf}
\caption{A flowchart describing the consistency tests presented in \Cref{sec:tests}. In all cases, we inverted or solved for crucial quantities (a process denoted by an arrow) that describe the expansion history of the Universe based on the $\Lambda$CDM\xspace model. 
As an example, we can solve the Hubble parameter $H(\Omega_\mathrm{m,0},z)$ of $\Lambda$CDM\xspace for $\Omega_\mathrm{m,0}$ and, assuming flatness, we can write the matter density parameter in terms of the dimensionless Hubble parameter $h(z)$, see \Cref{eq:omh1}. Clearly, $\mathrm{Om}_\mathrm{H}(z)\rightarrow \Omega_\mathrm{m,0}$, only if $\Lambda$CDM\xspace is the correct description of the expansion of the Universe. A similar process is also followed for the rest of the tests, which we can write in terms of the dimensionless distances $d_{\rm L}(z)$ and $d_{\rm A}(z)$.}
\label{fig:flow}
\end{figure}

It should be noted however that, as was shown in \citet[]{Sapone:2014nna}, the curvature test $\Omega_{\rm k}(z)$ diverges as $\sim1/z^2$ for small redshifts, hence we will regularize the divergence and in what follows we will instead consider the quantity $z^2 \Omega_{\rm k}(z)$. This is also similar to the complementary test of the Copernican principle, also using the distances and the Hubble rate, which was explored in \citet[]{Arjona:2021zac}. 

Finally, in \Cref{fig:flow} we show a flowchart that systematically organizes the different tests we consider. See \citet{Clarkson:2012bg} for a review of these.

\section{Analysis method of  currently available data\label{sec:cur_data}}

\subsection{Currently available data }
In order to constrain deviations from the $\Lambda$CDM\xspace model and the FLRW metric, we will now analyze a set of data providing information on the Hubble rate $H(z)$, the luminosity $D_{\rm L}(z)$ and the angular diameter distance $D_{\rm A}(z)$. Thus, in what follows, we will describe the currently available BAO and SNe data and present our analysis method.

We first consider the BAO data, given via measurements of the ratio $d_z$, which is defined as
\begin{equation}\label{eq:dz}
d_z\equiv \frac{r_{\rm s}(z_{\rm d})}{D_{\rm V}(z)}\,,
\end{equation}
where $D_{\rm V}$ is the spherically averaged distance
\begin{equation}
D_{\rm V}(z)=\left[(1+z)^2\, D_{\rm A}^2(z)\, \frac{c\,z}{H(z)}\right]^{1/3}\, ,
\end{equation}
while $r_{\rm s}(z_{\rm d})$ is the comoving sound horizon at the drag epoch,
\begin{equation}
r_{\rm s}(z_{\rm d})=\int_{z_{\rm d}}^\infty \frac{c_{\rm s}(z)}{H(z)} \; \text{d}z\, ,
\label{eq:sound-horizon-drag}
\end{equation}
$c_{\rm s}(z)$ is the sound speed of the baryon-photon plasma and $z_{\rm d}$ the redshift at the drag epoch, see Eq.~(4) of \citet{Eisenstein:1997ik} for a fitting function which is on average accurate to about $2\%$ compared to numerical estimates of $r_{\rm s}(z_{\rm d})$ obtained via recombination codes. In \citet[]{BOSS:2013rlg} and \citet[]{Aubourg:2014yra} similar approximations for the sound horizon were obtained  that are accurate to about $0.1\%$, while in \citet[]{Aizpuru:2021vhd} machine learning improved fitting functions, which are accurate to about $0.018\%$, were also presented.


The actual BAO data considered here, described in terms of $d_z$ given by \Cref{eq:dz}, the Hubble distance $D_{\rm H}(z)=c/H(z)$ and the comoving angular diameter distance $D_{\rm M}(z)= (1+z)\, D_\mathrm{A}(z)$, are provided by the surveys 6dFGS \citep{Beutler:2011hx}, WiggleZ \citep{Blake:2012pj}, DES Y3 \citep{DES:2021esc} and SDSS-IV eBOSS \citep{eBOSS:2020yzd}. In what follows, we shall refer to this combination of points as BAO and refer the interested reader to Appendix A of \citet[]{Martinelli:2020hud} and  Table 3 of \citet{eBOSS:2020yzd} for their exact values and their likelihood.

Note that we have taken care of the fact that a fiducial cosmology is usually adopted in order to convert the measured angular scales to distances for the BAO data; we have corrected the BAO results to be consistent with our models where necessary. However, we did not model the change in the nonlinear effects that may modify and damp the position of the BAO. This could introduce systematic uncertainties of a few percent in the reconstructions \citep[see][]{Angulo:2007fw} and may lead to lesser constraining power \citep[see for example][]{Anselmi:2017cuq}. The reason is that the BAO is essentially a feature at linear scales and the damping parameter is marginalized over. Moreover, modelling the nonlinear scales is non-trivial for theories beyond the $\Lambda$CDM\xspace model and the FLRW metric, thus is beyond the scope of this paper. Therefore, we choose not to correct for the aforementioned effects and leave that for future and more focused work on the subject.

We also consider the SNe, which can constrain the luminosity distance $D_{\rm L}(z)$. The key observable is the apparent magnitude $m(z)$
\begin{equation}
    m(z) = M_0+5\log_{10}{\left[\frac{D_{\rm L}(z)}{\rm Mpc}\right]}+25\, ,
    \label{eq:lum_dist}
\end{equation}
where $M_0$ is the intrinsic magnitude of the supernova at a redshift $z$. Clearly, if no external information is provided, the Hubble constant $H_0$ is completely degenerate with $M_0$ as the SNe data alone cannot constrain these two quantities. Here we will consider the updated Pantheon compilation of 1048 points from \citet{Scolnic:2017caz}, while for the likelihood we use the expression given in 
Appendix C of \citet{Conley:2011ku}, which is already marginalized over~$M_0$ and~$H_0$ and takes the covariances of the SNe data into account. 

\subsection{The genetic algorithms\label{sec:GA}} 



We now describe a machine learning approach, called the genetic algorithms (GA). The GAs are a set of stochastic optimization methods, usually applied in the context of non-parametric reconstructions of data. They are inspired by the notion of grammatical evolution and in particular, the genetic operations of mutation, i.e. a random change in an individual, and crossover, i.e. the combination of different individuals to form offspring. These operators behave as an environmental pressure, thus emulating the concept of natural selection in biology.

The ``reproductive'' success of a member of the population or in other words the probability that it will have offspring, is usually assumed to be proportional to its fitness, which is a measure of how accurately every member of the population fits the data. In this work we quantify the fitness of the individuals via a standard $\chi^2$ statistic. Overall, the GA have been used to probe for extensions of the standard model \citep[]{Akrami:2009hp}, deviations from 
$\Lambda$CDM,
both at the background and the perturbations level in linear order \citep[]{Nesseris:2012tt,Arjona:2020kco,Arjona:2019fwb}, to reconstruct a variety of cosmological data, such as SNe data \citep[]{Bogdanos:2009ib, Arjona:2020doi} or to perform reconstructions of null tests, e.g. the so called Om statistic and the curvature test \citep[]{Nesseris:2010ep, Nesseris:2013bia, Sapone:2014nna}.

In this work we perform a simultaneous fit of the SNe and BAO data with the GA, both for the currently available and the mocks presented in later sections, via the following approach. During the initialization phase, a set of functions is randomly chosen among an orthogonal basis that includes polynomials, exp, log and other elementary functions, but also the  operators $+,-,\times,\div, \wedge$, so that every member corresponds to a random guess for the Hubble rate $H(z)$ and the luminosity distance $D_\textrm{L}(z)$. Although we still assume the validity of the distance duality relation between the luminosity and angular diameter distances, $D_\textrm{L}(z)=(1+z)^2\, D_\textrm{A}(z)$, as this case was tested in \citet{Martinelli:2020hud}, we will however relax the assumption of \Cref{eq:distances}, as if it is valid, then the global shear test $\Sigma(z)$ is automatically satisfied regardless of the reconstruction method (GA etc). Thus, the GA will evolve the three functions of redshift that correspond to the Hubble rate $H(z)$, the luminosity distance $D_\textrm{L}(z)$ and the angular diameter distance, given by $D_\textrm{A}(z)=D_\textrm{L}(z)/(1+z)^2$. At this point we also impose some physical priors, such as that the luminosity distance at $z=0$ is zero, but we make no assumption on a DE model. 

After this initial population of function is set up, we then calculate every member's fitness via a $\chi^2$ statistic, using simultaneously as input the SNe and BAO data and their individual covariances, thus assessing the global fitness of the set of the three reconstructed functions. Next, using the  tournament selection approach, see \citet{Bogdanos:2009ib} for more details, a random set of the best-fitting functions in every generation is chosen and the crossover and mutation operators are applied to the selected functions. Finally, we iterate this process thousands of times in order to ensure a good convergence of the algorithm, but we also make runs with different random seeds, in order to avoid biasing the results due to the choice of a specific random seed.

The final output of the GA, after it has converged, is a set of three continuous and differentiable functions of redshift for the Hubble parameter $H(z)$, the luminosity distance $D_\textrm{L}(z)$ and the angular diameter distance $D_\textrm{A}(z)$. In the case of the currently available data we also minimize numerically the $\chi^2$ over the combination $r_{\rm s}(z_{\rm d}) h$, with $r_{\rm s}(z_{\rm d})$ the comoving sound horizon at the drag epoch and $h=H_0/(100\ {\rm km\ s^{-1}\ Mpc^{-1}})$, in order to avoid any assumptions on the Hubble constant for the BAO physics at early times. However, this complication does not exist for the fiducial BAO data, as they have been already marginalized over that quantity and the angular diameter distance is directly available. On a different note, as the BAO data cannot constrain independently the combination $\Omega_{\rm b, 0}h^2$, where necessary we will assume the value $\Omega_{\rm b, 0} h^2=0.02225$ from Planck 2018 \citep{Aghanim:2018eyx}.

The uncertainties of the best-fit Hubble parameter $H(z)$ and the luminosity distance $D_\textrm{L}(z)$ are obtained from the GA via a path integral
approach developed by \citet{Nesseris:2012tt,Nesseris:2013bia}. Specifically, the uncertainties are calculated by functionally integrating the likelihood, i.e. performing a path integral, over the whole functional space covered by the GA. This approach has been extensively tested by \citet{Nesseris:2012tt} and has been found to be in agreement with error estimates obtained via a bootstrap Monte-Carlo. In this work we use the publicly available code \texttt{Genetic Algorithms} for the numerical implementation of the GA.\footnote{\url{https://github.com/snesseris/Genetic-Algorithms}}

\subsection{Binning and parametrized approach\label{sec:parapp}}

We can now describe our parametrized methodology for the null tests presented in the previous section, so as to probe for deviations from $\Lambda$CDM\xspace and spatial homogeneity. To do so, we follow a two-pronged approach.

First, we bin the Hubble rate and distance data in bins with an equal number of points and then we calculate the various tests using the expressions in \Cref{sec:tests}. The first technique to explore the sensitivity of the tests presented in this work consists of evaluating them in redshift bins. Each test is sensitive to a particular combination of the different observables and parameters, hence the binning technique has to be performed individually to take into account the proper dependency and the covariances. To compute the $\mathrm{Om}_{\rm H}$ and $r_0$ tests we simply bin the $H(z)$ data at different redshifts and we propagate their uncertainties to the final tests by taking the diagonal elements of the covariance matrices of the data. We have checked, following the approach of \citet[]{Nesseris:2014vra}, that our choice about the covariance matrices does not affect the results by more than $1\%$.

Furthermore, to compute $\mathrm{Om}_{\rm dA}$ and $\mathrm{Om}_{\rm dL}$ we bin the corresponding distances; $D_{\rm A}(z)$ comes directly from the BAO mocks, whereas for $D_{\rm L}(z)$ we use the distance moduli of SNe by inverting \Cref{eq:lum_dist}. For the latter we further bin the data in 20 equally spaced redshift bins in order to have enough points and help compute the derivatives in a larger redshift range.

The uncertainties on the final $\mathrm{Om}_{\rm dL}$ test are obtained using a standard error propagation technique. Special care has been taken for the $\Sigma(z)$, $\Omega_{\rm k}(z)$, $\mathcal{O}_{\rm m}(z)$ and $\mathcal{O}_{\rm K}(z)$ tests. The global shear test $\Sigma(z)$ depends simultaneously on $H(z)$ and $D_{\rm A}(z)$ with its first derivative with respect to redshift, $D'_{\rm A}(z)$. The derivative of the angular diameter distance was obtained numerically considering two consecutive bins 
\begin{equation}
D'_{\rm A}(z_i) = \frac{D_{\rm A}(z_{i+1})-D_{\rm A}(z_{i})}{z_{i+1}-z_i}\,.
\label{eq:num_der_da}
\end{equation}
The uncertainties on the final test have been obtained constructing the Jacobian matrix which takes into account the derivatives of $\Sigma$ with respect to $H(z)$, $D_{\rm A}(z)$ and $D'_{\rm A}(z)$ and then propagating it with the covariance matrix of the data.

The curvature test $\Omega_{\rm k}(z)$ depends on $H(z)$, $d(z) = (1+z)\,c^{-1}H_0\,D_{\rm A}(z)$ and $d'(z)$, hence we need to recast the $\Omega_{\rm k}$ test as a function of the angular diameter distance in order to propagate the uncertainties directly from the data. Also here we adopt the numerical derivative for two consecutive bins for $d'(z)$, similar to \Cref{eq:num_der_da}, and the final uncertainties have been obtained as for the $\Sigma(z)$ test.

The $\mathcal{O}_{\rm m}(z)$ and $\mathcal{O}_{\rm K}(z)$ terms depend only on $H(z)$ and its first derivative, $H'(z)$. The latter is evaluated numerically in the same way and the uncertainties on the final tests are obtained using the propagation rule for both $H(z)$ and $H'(z)$. Finally, the number of bins for $\Sigma(z)$, $\Omega_{\rm k}(z)$, $\mathcal{O}_{\rm m}(z)$ and $\mathcal{O}_{\rm K}(z)$ is $n_{\rm data}-1$, due to the dependency of the test on the first derivatives of the observables.

Second, for each of the null tests under examination, we fit the binned data, derived in the previous step, with a parametric form that is inspired by the CPL parametrization for the dark energy equation of state, shown in \Cref{eq:CPLpar}, in order to obtain forecast constraints from \Euclid on the deviations from the expected values and on the redshift trends of such null tests.

Specifically, each quantity being tested ($P$) is parametrized with a functional form 
\begin{equation}
P(z)=P_0+P_1 \frac{z}{1+z},\label{eq:CPL}
\end{equation}
where $P_0$ is its value today and $P_1$ the first derivative, both evaluated at $a=1$. This simple parametrization allows us to characterize, in a simple and model-independent fashion, how \Euclid will be able to carry out the tests. The only exception to this form is the one chosen for the test of \Cref{eq:OK}, for which we multiply the parametrization above by a factor $z^2$ in order to ensure there is no divergence at the present time.

Using the binned values, obtained in the first step of this approach, as input data, we fit the predictions given by \Cref{eq:CPL}, sampling the two free parameters through the publicly available sampler \texttt{Cobaya} \citep{Torrado:2020dgo}, which exploits the Metropolis-Hastings algorithm presented in \citet{Lewis:2002ah,Lewis:2013hha}. In each case, we assume flat priors on the two free parameters of the analysis. The results of this parametrized approach will be denoted as (PA) in what follows.


\begin{table}
\begin{center}
\setlength{\tabcolsep}{4pt}
\renewcommand{\arraystretch}{1.35}
\caption{Parameter values for the fiducial models that are used for the mocks. The values used for the $\Lambda$CDM\xspace follow the fiducial of \citetalias{IST:paper1}; in particular, spatial flatness is assumed. $H_0$ is shown in units of km s$^{-1}$~Mpc$^{-1}$. \label{tab:fiducials}}
\begin{tabular}{lcccccccc}
\hline
\hline
model & $M_0$ & $\Omega_{\rm m,0}$ & $\Omega_{\rm b,0}h^2$ & $H_0$ & $w_0$ & $w_a$ & $\delta_0$ & $z_B$\\ \\
 \hline
$\Lambda$CDM\xspace &$-19.3$ & $0.32$ & $0.02225$ & $67$ & $-1$ & $0$ & - & -\\
CPL& $-19.3$ & $0.32$ & $0.02225$ & $67$ & $-0.8$ & $-1$& - & - \\
$\Lambda$LTB\xspace &$-19.3$ & $0.32$ & $0.02225$ & $67$ & $-1$ & $0$&  $-0.65$ & $1.5$ \\
\hline
\hline
\end{tabular}\\
\end{center}
\end{table}

\section{Analysis method of  mock data\label{sec:fore_data}}

\subsection{The fiducial cosmologies}

In order to forecast the ability of \Euclid to improve upon the sensitivity of the null tests mentioned in \Cref{sec:tests} we rely on mock data based on a priori known fiducial cosmologies, which we will now briefly describe. In particular, we will use the $\Lambda$CDM model, the CPL model ($w_0w_a$CDM) and the $\Lambda$LTB\xspace model to create simulated data sets for both the SNe and the BAO measurements based on the \Euclid specifications.

First, for the fiducial cosmologies based on the $\Lambda$CDM\xspace and CPL models we use the values of the parameters shown in \Cref{tab:fiducials}. Both of these fiducial cosmologies assume no violation of the spatial homogeneity and isotropy of the FLRW metric and the $\Lambda$CDM\xspace fiducial was also used in \citetalias{IST:paper1}. In this case, we can jointly write the Hubble expansion rate, at late times when we safely neglect radiation and assuming flatness ($\Omega_{\rm k,0}=0$), for both models as
\begin{equation}
\frac{H^2(a)}{H^2_0}=\Omega_{\rm m,0}\;a^{-3}+(1-\Omega_{\rm m,0})\; \exp\left[-3\int_1^a\frac{1+w(\alpha)}{\alpha}\;
{\rm d}\alpha\right]\,,
\label{eq:HlcdmDE}
\end{equation}  
where the dark energy equation of state for the CPL model follows \Cref{eq:CPLpar}, which includes the $\Lambda$CDM\xspace model for $(w_0,w_{a})=(-1,0)$, while the scale factor is $a= 1/(1+z)$. For the CPL fiducial we chose the parameters $(w_0,w_{a})=(-0.8,-1)$ as an extreme case, which corresponds to the combination  TT,TE,EE+lowE+lensing+SNe+BAO, see Fig.~30 of \citet[]{Aghanim:2018eyx}. Furthermore, we will also assume the values for the SNe absolute magnitude $M_0$, the baryon density $\Omega_{\rm b,0}h^2$ and the Hubble rate today $H_0$ as given in \Cref{tab:fiducials}.

Regarding the $\Lambda$LTB\xspace mock, considering the LTB metric of \Cref{eq:LTB} and using Einstein's equations, we obtain the modified Friedmann equation
\begin{equation}
H_\bot^2(r,t)
\equiv
\left[\frac{\dot{R}(r,t)}{R(r,t)}\right]^2 = \frac{2m(r)}{R^3(r,t)}+\frac{2 r^2k(r)\,M^2}{R^2(r,t)}+\frac{\Lambda\; c^2}{3}, \label{eq:LLTB_fried}
\end{equation}
where the dot is a time derivative and the radial dependent matter density
\begin{equation}
\rho_{\rm m}(r,t)=\frac{m'(r)}{4\pi r^2 R'(r,t)R(r,t)^2},
\end{equation} 
where $M$ is an arbitrary mass scale
and $m(r)$ is the so-called Euclidean mass function.  Note that we have conveniently recast $K(r) \equiv - 2 r^2 k(r) M^2$, where $k(r)$ is the curvature profile of the model.\footnote{Here, $H_\bot$ and $M$ have units of $\mathrm{Mpc}^{-1}$, $m$ has units of $\mathrm{Mpc}$ and that $k(r)$ is dimensionless.} Apart from the mass function, $m(r)$, and the curvature profile, $k(r)$, the LTB solution introduces another free function: the so-called Big Bang function, $t_{\rm BB}(r)$, which appears as the constant of integration of \Cref{eq:LLTB_fried}, defined via
\begin{equation}
t - t_{\rm BB}(r) = \int^{R(r,t)}_0 \frac{d x}{\sqrt{2m(r)x^{-1} + 2 r^2k(r)M^2 + \frac{\Lambda}{3} x^2}} .
\end{equation}
We also assume the compensated curvature profile:
\begin{align}
k(r) &= k_{\rm b}+(k_{\rm c}-k_{\rm b})P_3(r/r_{\rm B})\,,\\
P_{n}(x) &=
\begin{cases}
1-\exp\left[-\frac{1}{x}(1-x)^n\right]
& 0 \le x<1, \\
0 & 1\le x.
\end{cases}
\end{align}
where $k_{\rm b}$ and $k_{\rm c}$ are the curvature outside and at the center of the spherical inhomogeneity, respectively, and $r_{\rm B}$ is the comoving radius of the inhomogeneity. In order to avoid the decaying modes on the matter density field we set the Big Bang function to $t_{\rm BB} (r) = 0$.  The absence of decaying modes of the $\Lambda$LTB model ensures the agreement with the standard scenario of inflation. The remaining arbitrary function, $m(r)$, is effectively a gauge choice, which we fix here by $m(r) = 4\pi M^2 r^3/3$ (see \citealt[][]{Biswas:2010xm} for more details).

Since we assume a compensated model, the $\Lambda$LTB\xspace model becomes exactly $\Lambda$CDM\xspace at scales $r>r_{\rm B}$. It is specified by the six usual $\Lambda$CDM\xspace parameters plus the parameters introduced by the profile $k(r)$, namely $k_{\rm c}$ and $r_{\rm B}$
(from \Cref{eq:LLTB_fried} note that $k_{\rm b} \equiv 4 \pi \Omega_{\rm k,0}/3 \Omega_{\rm m,0}$). Following \citet[]{Camarena:2021mjr} we map $k_{\rm c}$ into the central density contrast $\delta_0 \equiv \rho_{\rm m}(0,t_0)/\rho_{\rm m}(r_{\rm B},t_0) -1$  and $r_{\rm B}$ into its corresponding redshift $z_{\rm B}$. We provide the relevant parameters used for the $\Lambda$LTB\xspace mocks in \Cref{tab:fiducials}, and in
\Cref{fig:LLTB} we show the matter density contrast, defined as $\delta \rho_{\rm m} \equiv \rho_{\rm m}(r,t_0)/\rho_{\rm m}(r_{\rm B},t_0) -1$.

\begin{figure}[!t]
\centering
\includegraphics[width = 0.44\textwidth]{figures/LLTB_Delta-rho.pdf}
\caption{The matter density contrast of the $\Lambda$LTB\xspace model as a function of the redshift $z$ for $\delta_0=-0.65$ and $z_B=1.5$.}
\label{fig:LLTB}
\vspace*{0cm}
\end{figure}

Having established the three cosmologies that we are going to examine, we then calculate the fiducial background quantities for the Hubble rate, the luminosity and angular diameter distances as a function of the redshift $z$. Finally, using the specifications of forthcoming surveys, as described in what follows, we create our mock BAO and SNe data.

Note that while a high precision is expected from forthcoming surveys, such as \Euclid, it is crucial to ensure they are accurate and several efforts have been made to take into account observational systematic uncertainties that will affect the surveys. For example, in \citet{Paykari2020}, the observational systematic effects of the \Euclid VIS instrument were studied, also taking into account the modelling of the point spread function and the charge transfer inefficiency. While some of the specifications for the future instruments could still change prior to their completion, here we will assume that by the time the data arrive, the aforementioned systematic effects will have been well understood. Thus, we will add all the relevant astrophysical systematic effects, such as the galaxy bias etc. as described in 
the following subsections.

Finally, it is worth stressing that the covariance matrices were computed assuming the fiducial $\Lambda$CDM cosmology. Therefore, we are neglecting the error due to the use of a non-$\Lambda$CDM fiducial cosmology. The computation of covariance matrices for alternative cosmologies is indeed an open issue in the exploitation of next-generation survey data \citep[]{Harnois-Deraps:2019rsd,Friedrich:2020dqo}.

\subsection{SNe surveys}

In this work we consider two future SNe surveys, one based on the proposed \Euclid DESIRE survey \citep{Laureijs:2011gra,Astier:2014swa} and the other based on the specifications of the LSST. We will assume that the \Euclid DESIRE survey will observe $1700$ supernovae in the range $z\in[0.7,1.6]$, while the LSST survey will observe $8800$ supernovae in the redshift range $z\in[0.1,1.0]$, for a total of 10500 data points.

In both cases we will assume the redshift distributions described in \citet{Astier:2014swa} and we further assume they are uncorrelated with each other. While the DESIRE survey is not a guaranteed output of \Euclid, we include it here in order to extend the redshift range of LSST, as this is crucial for the GA reconstruction at high $z$. For every SNe event we will assume an observational error of the form
\begin{equation}
    \sigma_{\textrm{tot},i}^2=\delta \mu^2_i+\sigma^2_{\textrm{flux}}+\sigma^2_{\textrm{scat}}+\sigma^2_{\textrm{intr}}\,,
\end{equation}
where the flux, scatter, and intrinsic contributions are the same for each event: $\sigma_{\textrm{flux}} = 0.01$, $\sigma_{\textrm{scat}} = 0.025$, and $\sigma_{\textrm{intr}} = 0.12$, respectively. We also include an error on the distance modulus $\mu=m-M_0$ that evolves linearly with the redshift $z$ as
\begin{equation}
    \delta\mu = e_{\rm M}~z\,,
\end{equation}
where we assume that the parameter $e_{\rm M}$ is drawn from a Gaussian distribution with zero mean and standard deviation $\sigma(e_{\rm M})=0.01$ \citep[see][]{Gong:2009yk,Astier:2014swa}. The aforementioned distance modulus error takes into account the possible redshift evolution of SNe, which has not been accounted for by the distance estimator, see \citet{Astier:2014swa}. Note that  $e_{\rm M}=0.01$ is needed to allow for systematic evolution - but it would just add in quadrature to the effective $e_{\rm M}=0.055$ coming from lensing, to make a single effective $e_{\rm M}$ that is negligibly different from 0.055. 


The uncertainty due to lensing was estimated theoretically to be $\sigma_{\rm lens}= 0.052 z$ \citep[][]{Amendola:2013twa,Quartin:2013moa} and $\sigma_{\rm lens}= 0.056 z$ \citep[][]{Ben-Dayan:2013nkf}, and observationally with Supernova Legacy Survey data to be $\sigma_{\rm lens}= (0.055 \pm 0.04) z$ \citep[][]{Jonsson:2010wx} and $\sigma_{\rm lens}= (0.054 \pm 0.024) z$ \citep[][]{SNLS:2010rmd}.

\subsection{Large-scale structure surveys}

As we are interested in forecasting the sensitivity of the null tests for $\Lambda$CDM\xspace with \Euclid, we will now describe how we simulate BAO data using the Fisher matrix approach. For this we follow the same methodology as was used in \citetalias{IST:paper1} for the spectroscopic survey.

In this work we will mainly focus on the spectroscopic Euclid survey as our goal is to obtain precise measurements of the Hubble parameter $H(z)$ and the angular diameter distance $D_\mathrm{A}(z)$. Overall, the Euclid survey will be able to probe the galaxy power spectrum in the redshift range $z \in [0.9,1.8]$ where, as mentioned in \citetalias{IST:paper1}, the main targets are ${\rm H}_{\alpha}$ emitters. In this case, \Euclid will measure up to $30$ million spectroscopic redshifts with an uncertainty given by $\sigma_z = 0.001(1 + z)$ \citep{Pozzetti:2016cch} and the main observable will be the galaxy power spectrum. The latter carries information on the distortions due to the Alcock-Paczynski effect, the residual shot noise, the redshift uncertainty, the anisotropies due to redshift space distortions and on the galaxy bias. Moreover, nonlinear effects that distort the shape of the power spectrum, e.g. a nonlinear smearing
of the BAO feature, have also been included in the matter power spectrum \citep{Wang:2012bx}. In principle, other effects might also include a nonlinear scale dependent galaxy bias, see for example \citet{delaTorre:2012dg}. 

We use the same binning scheme as in \citet{Martinelli:2020hud}, which is different from that of \citetalias{IST:paper1}; specifically we use nine equally spaced redshift bins of width $\Delta z = 0.1$ instead of four. By rebinning the data as given in \citetalias{IST:paper1}, we find the following specifications for the galaxy number density $n(z)$, in units of Mpc$^{-3}$ and the galaxy bias $b(z)$: 
\begin{align}
n(z)&= \{2.04, 2.08, 1.78, 1.58, 1.39, 1.15, 0.97, 0.7, 0.6\}\times 10^{-4}\,,\nonumber \\
b(z)&= \{1.42, 1.5, 1.57, 1.64, 1.71, 1.78, 1.84, 1.90, 1.96\} \,.\nonumber
\end{align}
This binning scheme allows for more data points and considerably improves the machine learning analysis, as discussed in later sections.
However, we note that we tested this particular choice against that of \citetalias{IST:paper1} in \citet{Martinelli:2020hud}, where we found no statistically significant difference. 

In order to obtain the Fisher matrix for the full set of cosmological parameters, which is used to estimate the parameter covariance matrix and propagate the error estimates, we follow the procedure described in \citetalias{IST:paper1}. The analysis includes the following parameters: the four shape parameters $\{\omega_{\rm m}=\Omega_{\rm m,0}h^2$, $h$, $\omega_{\rm b}=\Omega_{\rm b,0}h^2$, $n_{\rm s}\}$, the two nonlinear parameters $\{\sigma_{\rm p},\,\sigma_{\rm v}\}$ (see \citetalias{IST:paper1}) and the five redshift dependent parameters $\{\ln D_{\rm A},\,\ln H,\,\ln f\sigma_8,\,\ln b\sigma_8,\,P_{\rm s}\}$ evaluated in each redshift bin, where $f\sigma_8\equiv f(z)\sigma_8(z)$ is the linear growth rate times $\sigma_8$, which measures the amplitude of the linear power spectrum at scales of $8 h^{-1} \mathrm{Mpc}$, while $b\sigma_8\equiv b(z)\sigma_8(z)$ and $P_{\rm s}$ characterize the galaxy bias and the shot noise respectively (see \citetalias{IST:paper1}). Doing so we can derive the expected uncertainties from the Euclid survey of both the Hubble parameter $H(z)$ and the angular diameter distance $D_{\rm A}(z)$ in each of the nine redshift bins, while we marginalize over all the other parameters. The final Fisher matrix in principle depends on the particular fiducial cosmology, but here we will assume that this dependence is at best weak.

As the \Euclid spectroscopic survey will only probe the redshift range $z \in [0.9,1.8]$, we will be limited in the range where both the BAO and the SNe data are available. Thus, we will complement our analysis by including the DESI survey, in order to be able to reconstruct the null tests at the full range of the available SNe data. DESI has started survey operations in 2021 and is scheduled to obtain optical spectra for tens of millions of quasars and galaxies up to $z\sim 4$, which will allow for BAO and redshift-space distortion analyses.

We follow the official DESI forecasts for both the Hubble parameter $H(z)$ and the angular diameter distance $D_{\rm A}(z)$ as described in \citet{DESI2016}. These forecasts have also been derived following a Fisher matrix approach, described in \citet{2014JCAP...05..023F}, which is the full anisotropic galaxy power spectrum, i.e. measurements of the matter power spectrum as a function of the angle with respect to the line of sight, but also redshift and wavenumber. Similarly to the \Euclid forecasts, this approach also includes all the available information from the two-point correlation function and not just the position of the BAO peak.

\begin{figure*}[!tbp]
\centering
\includegraphics[width = 0.97\textwidth]{figures/real_final_fig_v2.pdf}
\caption{The GA reconstructions of the null tests mentioned in \Cref{sec:tests} using the currently available BAO and Pantheon SNe data, as described in \Cref{sec:cur_data}. In all cases, the dashed line at zero corresponds to the best-fit $\Lambda$CDM model, described by \Cref{eq:Hlcdm} and the parameters $(\Omega_\mathrm{m,0}=0.303, h=0.659)$, the red line is the GA fit, and the orange shaded region corresponds to the $68.3\%$ GA uncertainties.}
\label{fig:current}
\vspace*{0cm}
\end{figure*}

Specifically, we consider the baseline DESI survey, covering 14\,000\,deg$^2$ and targeting bright galaxies (BGs), luminous red galaxies (LRGs), emission line galaxies (ELGs), and quasars in the redshift range $z\in [0.05,3.55]$, but with a precision that will explicitly depend on the target population. First, the BGs will be in the redshift range $z\in [0.05,0.45]$ given in five equally spaced redshift bins, while the next two targets, the LRGs and ELGs, will be in the range $z\in [0.65,1.85]$ in $13$ equally spaced redshift bins. Moreover, the Ly-$\alpha$ forest quasars will be given in the range $z\in [1.96,3.55]$ in $11$ equally spaced redshift bins. Finally, we assume that these measurements are not correlated with each other.

In our analysis we will only include the DESI data at late times that are not overlapping with the \Euclid points in order to avoid spurious correlations between the two surveys. Furthermore, as the SNe data from LSST + DESIRE will only reach redshift $z=1.6$, we will only include the DESI data from the full BGs survey and the LRGs and ELGs up to $z=0.9$, thus leaving out the Ly-$\alpha$ forest observations.

\section{Results from the current data\label{sec:currentresults}} 

In this section we will now present the reconstructions of the null tests of \Cref{sec:tests} using the currently available SNe and BAO data. In this case, as the BAO data are coming from a plethora of different surveys and in a variety of different forms, i.e. in terms of $d_z(z)$, $1/d_z(z)$ etc, we cannot bin them without introducing new assumptions on their statistical properties. Thus, we will only present the results of the GA analysis, following the methodology presented in \Cref{sec:GA}, and in order to simplify the discussion of the results we have split the tests in groups of related tests.

Specifically, in \Cref{fig:current} we show the GA reconstructions of the null tests where in all cases the dashed line at zero corresponds to the best-fit $\Lambda$CDM model, described by \Cref{eq:Hlcdm} and the parameters $(\Omega_\mathrm{m,0}=0.303, h=0.659)$, the red line is the GA fit, while the orange shaded region corresponds to the $68.3\%$ GA uncertainties. In the first row from the top of \Cref{fig:current} we show the GA reconstructions of the two tests for the matter density, which directly probe the Hubble expansion rate, namely $\textrm{Om}_\mathrm{H}(z)$ given by \Cref{eq:omh1} (left) and $\mathcal{O}_{\rm m}(z)$ given by \Cref{Om-hz} (right). As can be seen, both are in agreement with a constant value and the $\Lambda$CDM\xspace within the $68.3\%$ uncertainties, albeit $\mathcal{O}_{\rm m}(z)$ is closer to the $68.3\%$ confidence level boundary, due to the presence of the derivatives of Hubble parameter which amplify any deviations from the concordance value, see \Cref{Om-hz}.

In the second row from the top we show the reconstructions for the matter density tests that probe directly the cosmological distances, namely $\textrm{Om}_\textrm{dL}(z)$ given by \Cref{eq:OmdLnull} (left panel) and  $\textrm{Om}_\textrm{dA}(z)$ given by \Cref{eq:OmdAnull} (right panel). As can be seen, both of them are as well in good agreement with a constant value and the $\Lambda$CDM\xspace best fit within the errors. Next, in the third row from the top we show the two curvature tests that probe the spatial homogeneity of the Universe, namely $z^2 \Omega_{\rm k}(z)$ given by \Cref{eq:OK}  (left panel) vs $\mathcal{O}_{\rm K}(z)$  given by \Cref{OK-hz} (right panel). In this case, both tests are consistent with a flat universe at the $68.3\%$ confidence level, notice however how the uncertainties seem to increase with redshift for the $z^2 \Omega_{\rm k}(z)$ due to the extra $z^2$ term that suppresses the singularity at late times as explained in \Cref{sec:tests}. On the other hand,  $\mathcal{O}_{\rm K}(z)$ is closer to the $68.3\%$ confidence level boundary, as it is complementary to $\mathcal{O}_{\rm m}(z)$ and is subject to the same impact of the presence of the Hubble derivative in their expression.

Finally, in the bottom row we show the global shear test $\Sigma(z)$ (left panel) given by \Cref{eq:GlobalShear}, which tests the Copernican principle, and the $r_0(z)$ null test (right panel) that probes for dark matter-dark energy interactions and is given by \Cref{eq:r0null}. Again, both tests are consistent within the uncertainties with the assumption of the Copernican principle and no interactions in the dark sector. In some cases, e.g. the $\mathrm{Om}_\mathrm{H}(z)$ and the $r_0(z)$ tests, the errors regions increase at small redshifts due to the fact that these tests are ill-defined in the zero redshift limit, as discussed in \Cref{sec:tests}. Similarly, this also occurs for the $z^2\Omega_\mathrm{k}(z)$ test, but due to the extra $z^2$ which regularizes the singularity at $z=0$, it now occurs at large redshifts instead.

\section{Forecast results\label{sec:fore_res}}
Following the approach described in \Cref{sec:fore_data}, we now present the constraint estimates of our null tests using the SNe and BAO mock data for three different fiducial cosmologies, based on the vanilla flat $\Lambda$CDM\xspace model given by \Cref{eq:Hlcdm}, the CPL parametrization given by 
\Cref{eq:CPLpar} and the $\Lambda$LTB\xspace model described by \Cref{eq:LLTB_fried}, assuming the fiducial values for the parameters shown in \Cref{tab:fiducials}. In our analysis we will also consider two different cases, the first where we only use the \Euclid BAO data, spanning the range $z\in[0.9,1,8]$, in order to quantify the ability of the Euclid survey alone to constrain any deviations from the null tests, and a second case where we also include the DESI BAO data, covering the lower redshifts as discussed in \Cref{sec:fore_data}, in order to highlight the synergies between the two surveys. Note that both cases include all SNe data (\Euclid + LSST) and, similarly to  \Cref{sec:currentresults}, we will group the plots by row, splitting the tests in groups of related tests.

Finally, as in this case the mock BAO data are always given in terms of the angular diameter distance $D_\mathrm{A}(z)$ and the Hubble expansion rate $H(z)$, we are now able to also perform the binning analysis as discussed in \Cref{sec:parapp}, as a complementary model-independent analysis to that of the GA. The binned data produced by this approach are also used to obtain constraint estimates on the considered tests, using a parametric approach based on a CPL-like parametrization given by \Cref{eq:CPL}. Note however, that while we apply the GA simultaneously on the full SNe and BAO mock data, the PA is only applied on the binned data of each test.

\subsection{Results from the $\Lambda$CDM\xspace mocks}
First, we present in \Cref{fig:mockGAEuclid} the results of the GA reconstructions of the null tests described in \Cref{sec:tests} for the \Euclid-only $\Lambda$CDM\xspace mocks and in \Cref{fig:mockGAall} for all mock data (\Euclid + DESI). In all cases, both here and in the plots in later sections, the dashed line corresponds to the fiducial value of the test under consideration for $\Lambda$CDM, the red line is the GA fit, the blue line is the PA fit given by \Cref{eq:CPL} of the binned mock data (black points), while the orange and blue shaded regions correspond to the $68.3\%$ uncertainties of the GA and PA respectively. Note that as the binned mock data include a random realization of error, they are expected to fluctuate around the fiducial value, something which is indeed seen in both \Cref{fig:mockGAEuclid} and \Cref{fig:mockGAall}.

As can be seen in all of the panels of \Cref{fig:mockGAEuclid}, in all cases both the GA and PA reconstructions, along with the binned mock data are able to correctly predict the fiducial model to within the $95.45\%$ uncertainties. Given their agnostic nature, the uncertainties of the GA reconstructions (orange shaded regions) are somewhat larger than those of the PA (blue shaded regions), but are very similar to those of the binning (black points). An exception to this is the $\text{Om}_\mathrm{dA}(z)$ test which in the case of the GA is dominated by the systematic uncertainties of the SNe data, due to the joint reconstruction approach employed in this case, as was also observed in \citet[]{Martinelli:2020hud}.

On the other hand, when we also include the DESI BAO data at small and intermediate redshifts that are not overlapping with the \Euclid points, as seen in \Cref{fig:mockGAall}, we find that the uncertainties in the case of the binning approach dramatically increase due to the lower constraining power of DESI with respect to \Euclid at this redshift range. This affects the PA (blue shaded regions), as in some cases, such as for the $\Sigma(z)$ and the $\mathcal{O}_{\rm K}(z)$ tests, it misses the fiducial model (dashed line) by more than $1\sigma$. The reason for this is that the PA is anchored to the \Euclid points at intermediate redshifts which have smaller error bars than the \Euclid ones, and due to a lack of flexibility of the PA parametrization, discrepancies appear at low redshifts. This also affects the GA reconstruction for $\Sigma(z)$, but in the rest of the cases the GA efficiently predicts the correct cosmology, as it uses the full SNe and BAO data. In either case, these deviations are due to the presence of derivatives in the null tests, which tend to cause instabilities in the reconstructions. 

Overall, after analyzing the $\Lambda$CDM\xspace mocks we find that \Euclid, in combination with other surveys, will be able to improve current constraints by approximately a factor of two or three with the machine learning approach, while the binning and parametric approach will provide a further factor of two improvement over the GA results. 

\subsection{Results from the CPL mocks}
Next, we present in \Cref{fig:mockCPLGAEuclid} the results of the GA and parametric reconstructions of the null tests for the \Euclid-only CPL mocks and in \Cref{fig:mockCPLGAall} for all mock data (\Euclid + DESI). 

This case now is more interesting as, due to the nature of the DE model chosen, many of the tests that probe the expansion history of the Universe, such as the $\mathrm{Om}_\mathrm{H}(z)$, have a distinct evolution with redshift. In the case of the \Euclid-only mocks shown in \Cref{fig:mockCPLGAEuclid} we see that both approaches overall predict the fiducial cosmology within the errors, albeit in some cases there is a small $1\sigma$ deviation. This occurs for example in the case of the $\mathrm{Om}_\mathrm{H}(z)$  test with the GA, as due to the lack of points at small redshifts the GA cannot anchor very well at $z\sim 1$. On the other hand, in the case of the binning, for the $\mathcal{O}_{\rm m}(z)$ and  $\mathcal{O}_{\rm K}(z)$ tests the presence of derivatives in the expressions causes more scatter in the points, which is also propagated to the PA (blue line). This also affects the $r_0(z)$ test, where the effect is further amplified by the need to fix the baryon density parameter $\Omega_\textrm{b,0}$ to its {\it Planck} best-fit value.

When we also include the DESI BAO data, see \Cref{fig:mockCPLGAall}, we can extend the reconstructions of the tests to small redshifts where many of the tests, such as the fiducial value of the $\mathrm{Om}_\mathrm{H}(z)$ test, manifest large deviations from a constant value. In this case we find that while the GA follows the fiducial trend quite well, the PA (blue line) is anchored to the higher redshift \Euclid points, as these have smaller uncertainties compared to the DESI ones, and deviates from the fiducial model by several $\sigma$; a similar effect is visible also for the $r_0(z)$ test.

On the other hand, the GA is not able to reproduce the fiducial model (dashed line) in the case of the $\mathcal{O}_{\rm K}(z)$ test by several $\sigma$, while the PA follows the fiducial much better and only at $68.3\%$ difference. However, in either case of \Euclid-only or \Euclid + DESI BAO data, the global shear test is not affected as expected and the reconstructions are practically identical to that of the $\Lambda$CDM\xspace mocks.

\begin{figure*}[!htbp]
\centering
\includegraphics[width = 0.975\textwidth]{figures/mock_Euclid_LSST_final_fig_LCDM.pdf}
\caption{The GA reconstructions of the null tests mentioned in \Cref{sec:tests} using the \Euclid-only BAO and \Euclid + LSST SNe $\Lambda$CDM\xspace mock data. In all cases, the dashed line corresponds to the fiducial value of the corresponding test for $\Lambda$CDM, the red line is the GA fit, the blue line is the PA fit, the shaded regions correspond to the $68.3\%$ uncertainties, while the black points correspond to the binned \Euclid mock data.}
\label{fig:mockGAEuclid}
\end{figure*}

\begin{figure*}[!htbp]
\centering
\includegraphics[width = 0.975\textwidth]{figures/mock_Euclid_LSST_DESI_final_fig_LCDM.pdf}
\caption{The GA reconstructions of the null tests mentioned in \Cref{sec:tests} using the \Euclid + DESI BAO and \Euclid + LSST SNe $\Lambda$CDM\xspace mock data. In all cases, the dashed line corresponds to the fiducial value of the corresponding test for $\Lambda$CDM, the red line is the GA fit, the blue line is the PA fit, the shaded regions correspond to the $68.3\%$ uncertainties, while the black points correspond to the binned \Euclid and DESI mock data.}
\label{fig:mockGAall}
\end{figure*}


\begin{figure*}[!htbp]
\centering
\includegraphics[width = 0.975\textwidth]{figures/mock_Euclid_LSST_final_fig_CPL.pdf}
\caption{The GA reconstructions of the null tests mentioned in \Cref{sec:tests} using the \Euclid-only BAO and \Euclid + LSST SNe data for the CPL mock for $(w_0,w_a)=(-0.8,-1)$. In all cases, the dashed line corresponds to the fiducial value of the corresponding test for CPL, the red line is the GA fit, the blue line is the PA fit, the shaded regions correspond to the $68.3\%$ uncertainties, while the black points correspond to the binned \Euclid  mock data.}
\label{fig:mockCPLGAEuclid}
\end{figure*}

\begin{figure*}[!htbp]
\centering
\includegraphics[width = 0.975\textwidth]{figures/mock_Euclid_LSST_DESI_final_fig_CPL.pdf}
\caption{The GA reconstructions of the null tests mentioned in \Cref{sec:tests} using the \Euclid + DESI BAO and \Euclid + LSST SNe data for the CPL mock for $(w_0,w_a)=(-0.8,-1)$. In all cases, the dashed line corresponds to the fiducial value of the corresponding test for CPL, the red line is the GA fit, the blue line is the PA fit, the shaded regions correspond to the $68.3\%$ uncertainties, while the black points correspond to the binned \Euclid and DESI mock data.}
\label{fig:mockCPLGAall}
\end{figure*}


\begin{figure*}[!htbp]
\centering
\includegraphics[width = 0.975\textwidth]{figures/mock_Euclid_LSST_final_fig_LTB.pdf}
\caption{The GA reconstructions of the null tests mentioned in \Cref{sec:tests} using the \Euclid-only BAO and \Euclid + LSST SNe data for the $\Lambda$LTB\xspace mock. In all cases, the dashed line corresponds to the fiducial value of the corresponding test for the $\Lambda$LTB\xspace model, the red line is the GA fit, the blue line is the PA fit, the shaded regions correspond to the $68.3\%$ uncertainties, while the black points correspond to the binned \Euclid mock data.}
\label{fig:mockLTBGAEuclid}
\end{figure*}

\begin{figure*}[!htbp]
\centering
\includegraphics[width = 0.975\textwidth]{figures/mock_Euclid_LSST_DESI_final_fig_LTB.pdf}
\caption{The GA reconstructions of the null tests mentioned in \Cref{sec:tests} using the \Euclid + DESI BAO and \Euclid + LSST SNe data for the $\Lambda$LTB\xspace mock. In all cases, the dashed line corresponds to the fiducial value of the corresponding test for the $\Lambda$LTB\xspace model, the red line is the GA fit, the blue line is the PA fit, the shaded regions correspond to the $68.3\%$ uncertainties, while the black points correspond to the binned \Euclid and DESI mock data.}
\label{fig:mockLTBGAall}
\end{figure*}



\subsection{Results from the $\Lambda$LTB\xspace mocks}
Finally, we now present in \Cref{fig:mockLTBGAEuclid} the results of the GA and PA, for the \Euclid-only $\Lambda$LTB\xspace mocks and in \Cref{fig:mockLTBGAall} for all mock data (\Euclid + DESI).
This particular case is even more interesting, as all the tests present large deviations compared to the simple $\Lambda$CDM\xspace case, due to the void profile used. Specifically, when compared to the corresponding plots for the $\Lambda$CDM\xspace mocks, the $\Lambda$LTB\xspace mocks exhibit large deviations from a constant value at small redshifts, due to the shape of the void profile, in particular in the cases of $\mathrm{Om}_\mathrm{H}(z)$, $r_0(z)$, $\mathcal{O}_\mathrm{m}(z)$ and $\mathcal{O}_\mathrm{K}(z)$.

As can be seen in \Cref{fig:mockLTBGAEuclid} in the case of the \Euclid-only BAO data, the PA, given by \Cref{eq:CPL}, is not flexible enough to reproduce the fiducial model, see for example the top left panel of \Cref{fig:mockLTBGAEuclid} for the $\mathrm{Om}_\mathrm{H}(z)$ reconstruction or the $\mathrm{Om}_\mathrm{dL}(z)$ and $\mathrm{Om}_\mathrm{dA}(z)$ tests in the second row. The latter actually give very large values of the matter density, something which would be a smoking gun for the presence of voids. In the case of the $\mathcal{O}_{\rm K}(z)$ test the GA is much closer to the fiducial model, compared to the PA, while for the $r_0(z)$ test, again the deviation is partly due to fixing the value of the baryon density, as before. When we also include the DESI BAO data, the reconstructions of \Cref{fig:mockLTBGAall} clearly highlight the ability of our approach to fully harness the power of the null tests of \Cref{sec:tests}. 

In particular, we see that due to the choice of the $\Lambda$LTB\xspace profile, several of the tests, namely $\mathrm{Om}_\mathrm{H}(z)$, $\mathcal{O}_{\rm K}(z)$, $r_0(z)$, $\Sigma(z)$, have large deviations from a constant value, which the GA is able to capture correctly. Note however that due to the simultaneous fitting of the data with the GA, see \Cref{sec:GA}, there are some correlations between the Hubble parameter and the luminosity distance. As a result, the fits in the redshift range covered by \Euclid might differ between the \Euclid-only and \Euclid + DESI cases, due to the added DESI BAO data at low $z$ in the latter case which may shift the fit. This difference in the fits between the two data sets can be seen for example in the $\mathrm{Om}_\mathrm{dL}(z)$ test shown in \Cref{fig:mockLTBGAEuclid} and  \Cref{fig:mockLTBGAall} respectively for the two sets of mocks.

On the other hand, in these cases the PA again gets anchored at high redshifts due to the small uncertainties of the \Euclid data and, due to its simple redshift evolution, it misses the features at small $z$. The same thing, but in an opposite direction, happens to the reconstructions of the $\mathrm{Om}_\mathrm{dL}(z)$ and $\mathrm{Om}_\mathrm{dA}(z)$ tests, albeit in this case the both the GA and the parametric approach are always within the errors.

\section{Conclusions \label{sec:conclusions}}
In this work we constrained deviations from the Copernican principle, the spatial homogeneity of the Universe on large scales, and from the $\Lambda$CDM\xspace model, using a set of null tests -- identities that have to be equal to a constant value at all redshifts, using current and forecast data. For the latter, we mainly focused on the constraints coming from synergies between \Euclid and other surveys, namely DESI and LSST, using both BAO and SNe data.

In \Cref{sec:theory} we presented several mechanisms where fundamental assumptions of the standard cosmological model, such as spatial homogeneity and the Copernican principle, may break down, as in the case of, for instance, Bianchi or LTB models. For this reason, we employed two theory-agnostic reconstruction approaches, one based on machine learning, called the GA, and another based on binning the data and then fitting a simple CPL-like parametrization, which we called the PA. These approaches allowed us to derive constraints on any deviation from the standard model, making only few assumptions on the redshift trends of these deviations. Avoiding any further assumptions is particularly important as assuming a particular model and then fitting the data, introduces theoretical biases and, as we have shown, makes it prone to misinterpreting deviations from the standard model.

In order to probe for deviations from the Copernican principle, the spatial homogeneity of the Universe on large scales, and from the $\Lambda$CDM\xspace model, we used a set of null tests proposed over the years in the literature, which we summarized in \Cref{sec:tests}. These tests probe different aspects of the aforementioned assumptions and provide a holistic approach to quantifying whether we could detect any deviations from them.

We first considered the currently available SNe and BAO data and using the GA we found that, within the errors, all the tests are compatible with the null hypothesis, that is, no deviations from the assumptions of the standard cosmological model, see \Cref{sec:currentresults} and \Cref{fig:current}. This is to be expected to some level, given the large uncertainties of these data sets and the fact that the BAO points come from a plethora of different surveys, each with their own systematic uncertainties.

The fact that our current constraints are limited by the available BAO data, highlights the importance of future full-sky surveys, such as \Euclid. Thus, we  considered mock data based on the specifications of \Euclid and also allowed for synergies with DESI and LSST, in order to forecast deviations from the assumptions of the standard cosmological model for three different fiducial cosmologies: one based on the $\Lambda$CDM\xspace model, one based on the CPL model for parameters at the $95.45\%$ limit of the current Planck best-fit and finally, another mock based on the $\Lambda$LTB\xspace model, which describes an inhomogeneous universe with a cosmological constant, as described in \Cref{sec:fore_data}.

Performing a simultaneous fit of these mock SNe and BAO data with the GA we found in each case a factor of two to three of improvement over the current constraints for most of the consistency tests, assuming the $\Lambda$CDM\xspace mocks. Alongside the GA we also applied the PA, which relies on the choice of a parameterization for the trend in redshift of the tests. Overall, we find that the GA provides larger errors with respect to the PA and this is for two main reasons. First, the GA is a non-parametric approach, thus is more agnostic, and it spans/explores the much larger functional space, instead of just the usual parameter space as in MCMCs, something which tends to make the errors larger, see \citet[]{Nesseris:2012tt}.

Second, the choice of a parameterization for the PA creates a lack of flexibility, particularly evident in the cases of the CPL and $\Lambda$LTB\xspace mocks, which results in reconstructions that deviate significantly from expectations. Such a lack of flexibility also results in error bars that are tighter than those of the GA, due to the fact that the posterior of PA parameters drops sharply moving away from their best fit. These issues can in principle be avoided increasing the flexibility of the PA, for example considering higher order redshift terms in \Cref{eq:CPL}. However, in realistic situations, one would not know the degree of complexity needed for the reconstruction and we decided here to use a common way to parameterize a redshift trend. This stresses even more strongly the advantage of the GA which does not require any a priori assumption to obtain the required level of flexibility.

Finally, we limited our analysis only to the $\Lambda$CDM, CPL and $\Lambda$LTB\xspace as it is currently not feasible to go beyond these models, e.g. by considering extensions such as the Bianchi type-I or back-reaction models, due to the lack of validated codes that would be used to estimate the necessary functions (distances and expansion rates). Thus, it is quite difficult to test the response of our null tests to these extensions of the standard cosmological model, but with new developments in theoretical model simulation we could explore them in the future.


In summary, this work highlights the benefits of synergies between the \Euclid  BAO survey and external probes in constraining  any deviations from the fundamental assumptions of the standard cosmological model. In particular, we have demonstrated that such a BAO survey will make it possible to constrain deviations from the Copernican principle, the spatial homogeneity of the Universe and the $\Lambda$CDM\xspace model at an unprecedented level in the near future using both non-parametric and parametrized approaches.



\begin{acknowledgements} 
We are grateful to S.~\'Avila for interesting discussions. JGB, MM and SN acknowledge support from the research project  PGC2018-094773-B-C32, and the Centro de Excelencia Severo Ochoa Program SEV-2016-059. SN also acknowledges support from the Ram\'{o}n y Cajal program through Grant No. RYC-2014-15843.  MM also received  support from ``la Caixa'' Foundation (ID 100010434), with fellowship code LCF/BQ/PI19/11690015. DS acknowledges financial support from the Fondecyt Regular project number 1200171. The work of CJM 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 through FCT - Funda\c c\~ao para a Ci\^encia e a Tecnologia in the framework of the project POCI-01-0145-FEDER-028987.
DC thanks CAPES for financial support.
VM thanks CNPq and FAPES for partial financial support.
LL was supported by a Swiss National Science Foundation (SNSF) Professorship grant (No.~170547).
This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 888258.
CC is supported by the UK Science \& Technology Facilities Council (STFC) Consolidated Grant ST/P000592/1.
JPM and AdS acknowledge funding from FCT Projects No.
UID/FIS/04434/2020, No. CERN/FIS-PAR/0037/2019, No. PTDC/FIS-OUT/29048/2017 and IF/01135/2015. ZS acknowledges support from the IRAP and IN2P3 Lyon computing centers.
PF received the support of a fellowship from ``la Caixa'' Foundation (ID 100010434). The fellowship code is LCF/BQ/PI19/11690018.
IT acknowledges support from the Spanish Ministry of Science, Innovation and Universities through grant ESP2017-89838, and the H2020 programme of the European Commission through grant 776247.

\AckEC
\end{acknowledgements}



\bibliographystyle{aa}
