Null test for cosmic curvature using Gaussian process

The cosmic curvature $\Omega_{K,0}$ is an important parameter that is related to many fundamental problems in modern cosmology. Knowing whether the universe is spatially open ($\Omega_{K,0}>0$), flat ($\Omega_{K,0}=0$), or closed ($\Omega_{K,0}<0$) is crucial for us to understand its evolution and the property of dark energy. A flat universe is strongly favored by some cosmological observations. For instance, the Planck 2018 cosmic microwave background (CMB) observations combined with the baryon acoustic oscillations (BAO) measurements give $\Omega_{K,0}=0.0007\pm 0.0019$, suggesting that our universe is flat to a $1\sigma$ error of $2\times 10^{-3}$ (Aghanim et al., 2020). However, it was found that the Planck TT,TE,EE+lowE power spectra data alone favors a slightly closed universe, $\Omega_{K,0}=-0.044_{-0.015}^{+0.018}$ (Aghanim et al., 2020; Park and Ratra, 2019; Handley, 2021; Di Valentino et al., 2019). This deviation from a flat universe is interpreted as the undetected systematics, statistical fluctuation, or new physics beyond the $\Lambda$ cold dark matter ($\Lambda$CDM) model. Efstathiou and Gratton (Efstathiou and Gratton, 2020) recently revisited the issue and claimed that the Planck data are still consistent with a flat universe. Whether this crisis really exists is still under debate.

It should be pointed out that most of the curvature parameter estimations assume a specific cosmological model. However, there is a strong degeneracy between the curvature parameter and the dark energy equation of state $w(z)$, so it is difficult to constrain them simultaneously, which hinders our understanding of dark energy. Therefore, it is necessary to measure the cosmic curvature in a cosmological model-independent way. For this purpose, many novel and feasible methods have been proposed; see e.g., Refs. (Bernstein, 2006; Shafieloo and Clarkson, 2010; Li et al., 2014, 2016, 2019; Sapone et al., 2014; Räsänen et al., 2015; Scolnic et al., 2018; Yu and Wang, 2016; L’Huillier and Shafieloo, 2017; Liao, 2019; Rana et al., 2017; Wang et al., 2020a; Wei and Wu, 2017; Xia et al., 2017; Denissenya et al., 2018; Wei, 2018; Witzemann et al., 2018; Collett et al., 2019; Qi et al., 2019; Wei and Melia, 2020; Zhou and Li, 2020; Dhawan et al., 2021; Jesus et al., 2020; Vagnozzi et al., 2021; Zhao et al., 2021; Wei et al., 2022; Wei and Melia, 2022; Koksbang, 2022; Liu et al., 2022; Zhang et al., 2022).

In Ref. (Cai et al., 2016), Cai et al. proposed a model-independent method to test whether the cosmic curvature deviates from zero. They adopted the observational data to reconstruct the reduced Hubble parameter $E(z)$ and distance-redshift relation $[D(z),D^{\prime}(z)]$, and then combined the reconstructions to perform the null test of $\Omega_{K,0}$. In their analysis, the $H(z)$ measurements from cosmic chronometer (CC) and BAO observations as well as the Union2.1 type Ia supernovae (SNe Ia) sample are considered, and the result favors a flat universe. Due to the increase of observational data, Yang and Gong (Yang and Gong, 2021) reperformed the null test using the CC $H(z)$ data and SNe Ia Pantheon compilation, and the result is still consistent with a flat universe. In Ref. (Yang and Gong, 2021), the BAO $H(z)$ data are not considered. It should be pointed out that the CC data may not constitute a reliable source of information due to some concerns (Kjerrgren and Mortsell, 2021). In contrast, the $H(z)$ data from BAO observations are much more accurate and reliable. In the present work, we consider both the CC and BAO $H(z)$ measurements, i.e., we shall use the latest CC+BAO $H(z)$ data and SNe Ia Pantheon sample to test the spatial flatness of the universe.

Furthermore, we also explore what role the future gravitational wave (GW) standard sirens and CC+BAO+redshift drift (RD) observations will play in the null test of $\Omega_{K,0}$. GWs can serve as standard sirens, since the GW waveform carries the information of the luminosity distance $D_{L}$ to source (Schutz, 1986; Holz and Hughes, 2005; Zhang, 2019). If the source’s redshift can be determined, for example, by identifying the electromagnetic counterpart of the GW event, we can then establish the $D_{L}$-redshift relation. We simulate the GW $D_{L}(z)$ data based on the planned space-based GW detector, DECihertz Interferometer Gravitational wave Observatory (DECIGO) (Kawamura et al., 2021). In the coming decades, with the advent of some powerful optical and radio telescopes, such as the Euclid (Amendola et al., 2018), Subaru Prime Focus Spectrograph (PFS) (Ellis et al., 2014), Dark Energy Spectroscopic Instrument (DESI) Aghamousa et al. (2016a, b), and Square Kilometre Array (SKA) Bacon et al. (2020), we can better measure the Hubble parameter using the CC and BAO methods. We simulate the CC+BAO $H(z)$ data in the redshift interval $0<z<2.5$ according to the observational data. In the future, another promising way to measure $H(z)$ is the RD method (Loeb, 1998). We simulate the high-$z$ Hubble data based on hypothetical RD observations of the upcoming European Extremely Large Telescope (E-ELT). We shall use the simulated GW and CC+BAO+RD data to perform the null test of $\Omega_{K,0}$ and compare the result with that using the current CC+BAO and SNe Ia data.

In this work, we adopt a machine learning method, the Gaussian process (GP), to reconstruct the cosmological functions. It has been widely used in cosmological researches; see e.g., Refs. (Holsclaw et al., 2010; Shafieloo et al., 2012; Seikel et al., 2012; Yahya et al., 2014; Yang et al., 2015; Cai et al., 2016; Wang et al., 2017a; Elizalde and Khurshudyan, 2019; Elizalde et al., 2018; Liao et al., 2019; Cai et al., 2020; Mukherjee and Banerjee, 2021a, 2022a; Mehrabi and Rezaei, 2021; Aljaf et al., 2021; Vazirnia and Mehrabi, 2021; Mukherjee and Mukherjee, 2021; Mukherjee and Banerjee, 2021b; Ruiz-Zapatero et al., 2022; Hwang et al., 2022; Mukherjee and Banerjee, 2022b; Elizalde et al., 2022; Elizalde and Khurshudyan, 2022). The GP method allows one to reconstruct a function and its derivative from data without assuming any particular parametrization, so it is suitable for our purpose. The artificial neural network (ANN) has recently emerged as a promising tool for reconstructing functions (Wang et al., 2020b, 2021; Liu et al., 2021; Qi et al., 2023), however, as far as we know, ANN has difficulties in reconstructing the derivative of a function. For this reason, our research is based on the GP analysis.

The remainder of this paper is organized as follows. We briefly describe the methodology in Sec. II. Sec. III contains the data we adopted. We present the results and make some discussions in Sec. IV. Finally, we give our conclusions in Sec. V.

In the homogeneous and isotropic universe, the FLRW metric is applied to describe its spacetime:

where $c$ is the speed of light, $a$ is the scale factor, and $K$ is a constant that is related to the cosmic curvature by $\Omega_{K}\equiv-Kc^{2}/(aH)^{2}$, with $H$ being the Hubble parameter. We use $\Omega_{K,0}$ to represent the present value of $\Omega_{K}$, and then $\Omega_{K,0}>0$, $\Omega_{K,0}=0$ and $\Omega_{K,0}<0$ correspond to the open, flat and closed universe, respectively. The luminosity distance can be expressed as

where

and $E(z)\equiv H(z)/H_{0}$ is the reduced Hubble parameter.

Differentiating Eq. (2), we have (Clarkson et al., 2007, 2008)

where $D(z)=(H_{0}/c)D_{L}(z)(1+z)^{-1}$ is the dimensionless comoving distance. Obviously, the curvature parameter can be directly determined by using the Hubble parameter and luminosity distance according to Eq. (7). Thus, we can perform the null test of $\Omega_{K,0}$. Note that $D(0)=0$ will bring a singularity at $z=0$. For simplicity, we transform Eq. (7) to

We can see that the left-hand side of Eq. (8) is non-zero when $z\neq 0$ if $\Omega_{K,0}$ is nonvanishing. Therefore, the null test of $\Omega_{K,0}$ is equivalent to the null test of the left-hand side of Eq. (8). Following Cai et al. (Cai et al., 2016), we define

For a spatially flat universe,

is always true at any redshift, so the deviation from it will imply a nonvanishing cosmic curvature. To carry out the null test of Eq. (10), we need to reconstruct the functions of $E(z)$, $D(z)$, and $D^{\prime}(z)$ using the observational data. In this work, we adopt the GP method to reconstruct the cosmological functions without assuming a specific cosmological model. Therefore, a cosmological model-independent test of whether $\Omega_{K,0}$ deviates from zero is performed in this work; note here that we do not assume a specific dark energy model, but of course an FLRW model of the homogeneous and isotropic universe is adopted.

GP is a non-parametric smoothing method for reconstructing functions (Holsclaw et al., 2010; Shafieloo et al., 2012; Seikel et al., 2012; Yahya et al., 2014; Yang et al., 2015; Cai et al., 2016; Wang et al., 2017a; Elizalde and Khurshudyan, 2019; Elizalde et al., 2018; Liao et al., 2019; Cai et al., 2020; Mukherjee and Banerjee, 2021a, 2022a; Mehrabi and Rezaei, 2021; Aljaf et al., 2021; Vazirnia and Mehrabi, 2021; Mukherjee and Mukherjee, 2021; Mukherjee and Banerjee, 2021b; Ruiz-Zapatero et al., 2022; Hwang et al., 2022; Mukherjee and Banerjee, 2022b; Elizalde et al., 2022; Elizalde and Khurshudyan, 2022), which assumes that at each point $z$, the reconstructed function $f(z)$ is a Gaussian distribution. Furthermore, the functions at different points are related by a covariance function. We can model a data set using a GP as

where $\mu(z)$ is the mean function which provides the mean of random variables at each observational point, and $k(z,\tilde{z})$ is the covariance function which correlates the values of different $f(z)$ at data points $z$ and $\tilde{z}$ separated by $|z-\tilde{z}|$ distance units.

There is a wide range of possible covariance functions, and the choice of covariance function actually affects the reconstruction to some extent (Seikel et al., 2012). Here, we consider the commonly used squared exponential covariance function,

where $\sigma_{f}$ denotes the overall amplitude of the oscillations around the mean and $\ell$ gives a measure of the correlation length between the GP nodes. Both $\sigma_{f}$ and $\ell$ are hyperparameters, which will be optimized by GP with the observational data. Given the GP for $f(z)$, the GP for the first derivative is consequently given by

Therefore, it is convenient to calculate the derivative of reconstructed function. In this work, we use the publicly available GaPP code to implement our analysis (Seikel et al., 2012).

The specific calculation process of GaPP is as follows. For a set of input points, $\bm{Z}=\left\{z_{i}\right\}$, the covariance matrix $K(\bm{Z},\bm{Z})$ is calculated by $[K(\bm{Z},\bm{Z})]_{ij}=k\left(z_{i},z_{j}\right)$. Even in the absence of observations, we can generate a random function $f(z)$ from GP, i.e., we can generate a vector $\bm{f^{*}}$ of function values at $\bm{Z^{*}}=\left\{z_{i}^{*}\right\}$ with $f_{i}^{*}=f\left(z_{i}^{*}\right)$:

where $\bm{\mu}^{*}$ is a prior of the mean of $\bm{f}^{*}$, which is set to zero in this paper. For observations, $\left\{(z_{i},y_{i},\sigma_{i})|_{i=1,\ldots,N}\right\}$, one can also use a GP to describe them. We stress that $y_{i}$ is assumed to be scattered around the underlying function, i.e., $y_{i}=f(z_{i})+\epsilon_{i}$, where Gaussian noise $\epsilon_{i}$ with variance $\sigma_{i}^{2}$ is assumed. Therefore, we need to add the variance to the covariance matrix,

where $\bm{y}$ is the vector of $y_{i}$ values and $C$ is the covariance matrix of the data. For uncorrelated data, we use $C=\text{diag}(\sigma_{i}^{2})$.

The above two GPs for $\bm{f}^{*}$ and $\bm{y}$ can be combined in the joint distribution:

Here $\bm{y}$ is known from observations. To reconstruct $\bm{f}^{*}$, one can consider the conditional distribution,

where

and

are the mean and covariance of $\bm{f}^{*}$, respectively. Eq. (23) is the posterior distribution of Eqs. (14) and (15). To reconstruct $\bm{f}^{*}$ using the above equations, we need to know the hyperparameters $\sigma_{f}$ and $\ell$, which can be determined by maximizing the logarithm marginal likelihood,

The hyperparameters will be fixed after being optimized through Eq. (II), then the reconstructed function $\bm{f}^{*}$ at the chosen points $\bm{Z}^{*}$ can be calculated from Eqs. (24) and (II). Therefore, GP is model-independent and without free parameters.

Here we present the data used for reconstructions. We do not consider to use the CMB data, because (i) we do not use observations to constrain a specific model, (ii) the reconstructions cannot be performed up to the early universe (e.g., the last scattering), and (iii) the inconsistencies in the measurements of the early- and late-universe observations (such as the well-known “Hubble tension”) should also be considered. Thus, we only use the late-universe observations in this work. In the first two subsections III.1 and III.2, we present the CC+BAO and SNe Ia real data, and in the last two subsections III.3 and III.4, we present the GW and CC+BAO+RD mock data.

We first need to normalize the observational Hubble parameter and luminosity distance data to get the $E(z)$ and $D(z)$ data. From Eq. (10), we have

where $\hat{H}_{0}$ is the normalization factor. Since the two factors can cancel out each other, whose value will not influence the null test of the cosmic curvature. In this work, we adopt $\hat{H}_{0}=70.0\ \rm km\ s^{-1}\ Mpc^{-1}$ to normalize the CC+BAO $H(z)$ data as the observational $E(z)$. For consistency, we use the same $\hat{H}_{0}$ to normalize the SNe Ia $D_{L}(z)$ data as the observational $D(z)$. We then adopt the GP method to reconstruct the functions of $E(z)$, $D(z)$, and $D^{\prime}(z)$, and the results are shown in Fig. 1. The black line is the mean of the reconstruction and the shaded grey regions are the $1\sigma$ ($68.3\%$) and $2\sigma$ ($95.4\%$) confidence level (C.L.) of the reconstruction. One may find that the errors of the reconstructed function are significantly smaller than the errors of the data themselves. This is due to the basic assumptions that the distribution of the function at each point is Gaussian and the data points are correlated by the covariance function. As can be seen, the error of $E(z)$ does not increase significantly with redshift due to the relatively accurate BAO data. However, the error of $D(z)$ and $D^{\prime}(z)$ becomes very large at $z>1.5$ because the data in that region are scarce and of poor quality. In addition, all the reconstructed functions in the redshift interval $z<1$ are consistent well with a flat $\Lambda$CDM model with $\Omega_{\rm m}=0.28$, which is adopted for a comparison.

With the GP reconstructions, we can then carry out the null test of the cosmic curvature. In Ref. (Cai et al., 2016), the authors applied the Monte Carlo sampling to determine $\mathcal{O}_{K}(z)$. Different from them, we use the error propagation formula to calculate the error of $\mathcal{O}_{K}(z)$ at each point $z$. The result is shown in the left panel of Fig. 2. The red dashed line refers to the spatially flat universe with $\Omega_{K,0}=0$. It can be seen that the result is consistent with the flat universe within the domain $0<z<2.3$, falling within the $1\sigma$ C.L. We note that the mean of $\mathcal{O}_{K}(z)$ is very close to zero in the redshift interval $0<z<1$, however, it becomes more and more negative at $z>1$. In addition, the universe with $\Omega_{K,0}=0$ almost falls out of the $1\sigma$ C.L. at $z>1.5$. This indicates that there is still a possibility for a spatially closed universe. Notably, the reconstruction at $z>1.5$ prefers a closed universe over an open one. We also plot the $\Lambda$CDM model with negative curvature in Fig. 2. As can be seen, a universe with $\Omega_{K,0}=-0.15$ can fall within the $1\sigma$ C.L. and a universe with $\Omega_{K,0}=-0.3$ can fall within the $2\sigma$ C.L. Note that the curves shown in Fig. 2 are plotted by assuming a $\Lambda$CDM model (with $\Omega_{\rm m}=0.28$) which is adopted here for a comparison.

In the GP analysis above, we only consider the squared exponential covariance function. In fact, the choice of covariance function will influence the result. To illustrate the effect, we take the Matern92 covariance function as an example, which is given by

Following the process described above, we perform the null test again, and the updated result is shown in the right panel of Fig. 2. The result is still consistent with the flat universe, falling within the $1\sigma$ C.L. At $z>1$, the mean of $\mathcal{O}_{K}(z)$ also deviates from zero and becomes negative. The difference is that the error of $\mathcal{O}_{K}(z)$ using the Matern92 covariance function is slightly larger than that using the squared exponential covariance function. In general, the difference does exist, but it is not significant enough to change our conclusions. In the following, we adopt only the GP with squared exponential covariance to complete our analysis.

The large error of $\mathcal{O}_{K}(z)$ leaves the big window open for a possible non-flat universe. Therefore, it is necessary to forge new cosmological probes to precisely measure the Hubble parameter and luminosity distance, thus tightening the constraints on the cosmic curvature. In the next decades, the GW standard siren and RD method will be greatly developed. We can adopt them to measure $D_{L}(z)$ and $H(z)$. Then, it is worth further studying the issue of combining the GW and RD observations with the traditional CC and BAO measurements to test the spatial flatness of the universe. We normalize the simulated CC+BAO+RD $H(z)$ data as the $E(z)$ data, and normalize the GW $D_{L}(z)$ data as the $D(z)$ data. Note that the error of RD $H(z)$ is relatively small and decreases with redshift, which is very helpful for us to reconstruct the cosmological function. Having obtained the data sets of $E(z)$ and $D(z)$, we adopt GP to reconstruct $E(z)$, $D(z)$, and $D^{\prime}(z)$, and the results are shown in Fig. 3. As can be seen, the error of $E(z)$ does not increase significantly even up to the redshift of $4.5$. In addition, the error of $D^{\prime}(z)$ reconstructed from the GW data is obviously smaller than that from the SNe Ia data. Here, we wish to note that although the RD method is promising in measuring the Hubble parameter, the actual measurement will be fairly challenging even with the powerful facilities such as the E-ELT.

With the reconstructions of $E(z)$ and $D^{\prime}(z)$, we use the error propagation formula to determine $\mathcal{O}_{K}(z)$, and the result is shown in the left panel of Fig. 4. We see that the mean of $\mathcal{O}_{K}(z)$ is very close to zero in the redshift interval $0<z<4.5$. Some weak deviations are mainly due to the consideration of Gaussian randomness in the simulations. The result strongly favors a flat universe, which is consistent with the assumed flat $\Lambda$CDM model. Note that we only consider $\mathcal{O}_{K}(z)$ at $0<z<4.5$. Even though we have reconstructed the function of $E(z)$ at $0<z<5$, the part of $z>4.5$ is extrapolated, whose accuracy cannot be guaranteed. On the other hand, the GW $D(z)$ data at $z>4.5$ are scarce and of poor quality, so the reconstruction is not convincing. We also plot the $\Lambda$CDM model with negative curvature in the left panel of Fig. 4. As can be seen, the reconstruction at $1.5<z<2.5$ can rule out the universe with $\Omega_{K,0}=-0.07$ at $1\sigma$ C.L. and the universe with $\Omega_{K,0}=-0.14$ at $2\sigma$ C.L. Similarly, we can rule out the universe with positive curvature in this way. We have tested that if the number of CC+BAO $H(z)$ data reaches 300, the reconstruction at $1.5<z<2.5$ can rule out the universe with $\Omega_{K,0}=-0.05$ at $1\sigma$ C.L. and the universe with $\Omega_{K,0}=-0.09$ at $2\sigma$ C.L. We compare the $1\sigma$ error of $\mathcal{O}_{K}(z)$ derived from the different data in the right panel of Fig. 4. It can be seen that the error of $\mathcal{O}_{K}(z)$ derived from the future {CC+BAO+RD, GW} data is significantly smaller than that from the current {CC+BAO, SNe Ia} data. Concretely, for example, the error provided by {CC+BAO+RD, GW} is less than that given by {CC+BAO, SNe Ia} at $z=1$ by 66.6% (and at $z=2$ by $87.3\%$). Moreover, the error of $\mathcal{O}_{K}(z)$ derived from {CC+BAO, SNe Ia} grows rapidly at $z>1.5$, while the error of $\mathcal{O}_{K}(z)$ from {CC+BAO+RD, GW} does not grow rapidly until $z\sim 4$. All these analyses indicate that with the synergy of multiple high-quality observations in the future, we can better determine the spatial topology of the universe.

In this work, we simulated only five RD measurements at high redshifts by observing the Ly-$\alpha$ absorption lines of QSOs. It should be pointed out that the SKA Phase I can measure the RD at $0<z<0.3$ by observing the H i emission lines of galaxies (Weltman et al., 2020; Liu et al., 2020; Qi et al., 2021). However, due to the excellent performance of CC+BAO in reconstructing $E(z)$, we did not consider this case. In addition, we did not consider the future observations of supernovae, because the GW data can reconstruct $D(z)$ and $D^{\prime}(z)$ very well, as shown in Fig. 3. We note that the GW standard siren method also has the potential to measure the Hubble parameter (Nishizawa et al., 2011). In principle, the observations of luminosity distances to GW sources across the sky should not be directional. However, mainly due to the local motion of the observer, there are tiny anisotropies in the luminosity distance, which enables us to measure $H(z)$. In this work, we considered a conservative scenario, i.e., $5000$ GW events with determined redshifts. In very optimistic scenarios, DECIGO is expected to detect $10^{5}\sim 10^{6}$ GW events with determined redshifts. If that can be done, we can reconstruct the functions of $D(z)$ and $D^{\prime}(z)$ with breathtaking precision and measure the Hubble parameter at $0<z\lesssim 3$ with a few percent accuracy (Nishizawa et al., 2011). Then we can test the spatial flatness of the universe in a cosmological model-independent way using only the GW data. We plan to explore this possibility in a future work.

In this paper, we adopt a cosmological model-independent method to test whether the cosmic curvature $\Omega_{K,0}$ deviates from zero. We use the Gaussian process method to reconstruct the reduced Hubble parameter $E(z)$ and the distance-redshift relation $[D(z),D^{\prime}(z)]$, independently. In the reconstruction, we do not assume any specific cosmological model. By combining the reconstructions of $E(z)$ and $D^{\prime}(z)$, we can determine $\mathcal{O}_{K}(z)$, which is zero at any redshift for a spatially flat universe with $\Omega_{K,0}=0$. Thus, we can carry out the null test of $\Omega_{K,0}$. We adopt the latest CC Hubble data, radial BAO Hubble data, and Pantheon SNe Ia data to implement our analysis.

Our result is consistent with a universe with $\Omega_{K,0}=0$ within the domain of reconstruction $0<z<2.3$, falling within the $1\sigma$ confidence level. We stress that the reconstruction favors a flat universe at $0<z<1$, however, it tends to favor a closed universe at $z>1$. In this sense, there is still a possibility for a closed universe. The error of the reconstructed function grows rapidly at $z>1.5$ due to the poor-quality observational data in that region, so it is necessary to forge new cosmological probes to precisely measure the luminosity distance and Hubble parameter. The GW standard siren and redshift drift observations that can be used to measure $D_{L}(z)$ and $H(z)$ will be greatly developed in the next decades. We simulated the GW standard siren and RD data based on the hypothetical observations of the upcoming DECIGO and E-ELT, respectively. The traditional methods for measuring the Hubble parameter are also promising, and we simulated the CC+BAO $H(z)$ data for the next decades. Combining these mock data, we performed the flatness test of the universe. We find that with the synergy of multiple high-quality observations in the future, we can tightly constrain the spatial geometry of the universe or exclude the flat universe with $\Omega_{K,0}=0$.