© ESO 2019

1. Introduction

Diverse sets of cosmological probes have now established that our Universe is undergoing an accelerated expansion: from local type Ia supernovae (SNe Ia; e.g. Riess et al. 1998; Perlmutter et al. 1999; Suzuki et al. 2012) to the combined constraints of the cosmic microwave background (CMB), baryonic acoustic oscillations (BAO), SNe Ia, and lensing provided by the Planck Collaboration VI (2018). The observed accelerated expansion is ascribed to the so-called dark energy (Ω_Λ), which is found to be ∼70% (within the standard concordance model) of the energy budget, thus the dominant component of today’s Universe. In the basic concordance model (i.e. ΛCDM), acceleration is driven by a cosmological constant, Λ, with an additional cold dark matter (Ω_M) component in a flat Universe. Although this model has been overall successful in fitting the data, observational evidence suggesting a deep revision of its underlying assumptions is now mounting. For example, the determination of the Hubble constant (H₀), which represents the normalisation of the Hubble parameter H(z), from local sources (Cepheids, supernovae Ia) is in tension at 4.4σ with respect to the one from the sound horizon observed from the CMB (Riess et al. 2019). If systematic uncertainties are not the main driver of this discrepancy, a significant disagreement would call for fundamental physics beyond the standard model (e.g. time-dependent dark energy equation of state, modified gravity, additional relativistic particles).

Recently, model-independent measurements of the distance modulus–redshift relation (the so-called Hubble diagram) of quasars combined with SNe Ia from the Sloan Digital Sky Survey-II/Supernova Legacy Survey 3 (SDSS-II/SNLS3) Joint Light-curve Analysis (JLA; Betoule et al. 2014) have shown a deviation from the ΛCDM model emerging at high redshifts (z >  1.4) with a statistical significance of ∼4σ (Risaliti & Lusso 2019; hereafter RL19).

In the present work, we present improved measurements of the Hubble diagram leveraging on the analysis presented by RL19. By combining the quasar sample of these latter authors with the most updated compilation of SNe Ia from the Pantheon survey (Scolnic et al. 2018) and publicly available gamma-ray burst (GRB) data (Demianski et al. 2017a), we are now finding a tension between the best-fit cosmographic parameters and the ΛCDM model at > 4σ.

2. The data set

To compute the cosmological parameters, we built a Hubble diagram by combining three samples:

1. Pantheon. A SNe Ia sample consisting of 1048 objects in the range 0.01 <  z <  2.26, significantly extending the redshift coverage at z >  1 with respect to JLA (i.e. maximum redshift z  ≃  1.3). This sample is a combination of 365 spectroscopically confirmed SNe Ia discovered by the Pan-STARRS1 (PS1) Medium Deep Survey together with the subset of 279 PS1 SN Ia (0.03 <  z <  0.68) with distance estimates from SDSS, SNLS, and various low-redshift and HST samples (see Table 4 in Scolnic et al. 2018). For further details on photometry, astrometry, calibration, and systematic uncertainties on SN Ia distances we refer the interested reader to Scolnic et al. (2018).

2. Quasars. The quasar sample is composed of 1598 sources in the redshift range 0.04 <  z <  5.1 with high-quality UV and X-ray flux measurements. Distance moduli are computed from the non-linear relation between the ultraviolet and the X-ray emission observed in quasars (Risaliti & Lusso 2015; Lusso & Risaliti 2016). All the details on sample selection, X-ray, and UV flux computation and the analysis of the non-linear relation, calibration, and a discussion on systematic errors are provided in RL19.

3. Gamma ray bursts. This sample is a compilation of 162 GRBs ranging from 0.03 <  z <  6.67. Due to convergency issues of the standard cosmographic approach (Capozziello et al. 2011), we limited the GRB redshift range to z = 6.67, thus excluding the two highest redshift objects at z = 8.1 and z = 9.3. The Hubble diagram is constructed by calibrating the correlation between the peak photon energy, E_p, i, and the isotropic equivalent radiated energy, E_iso (Amati et al. 2008). We refer to Demianski et al. (2017a,b) and Amati & Della Valle (2013) for details on the data set and on the calibration of the distance modulus–redshift relation, as well as for discussion and analysis of possible instrumental and selection effects (see also Amati 2006; Ghirlanda et al. 2008; Nava et al. 2012; Amati et al. 2019).

Figure 1 shows the redshift distribution of the three samples used to build the Hubble diagram. The combination of GRBs and quasars is such that we can probe a redshift range better suited for studying dark energy than the one covered by SNe Ia. Nevertheless, SNe Ia are necessary to cross-calibrate both quasars and GRBs to compute absolute distances.

Fig. 1. Redshift distribution of the Pantheon SNe Ia sample (green, Scolnic et al. 2018), the quasars sample (yellow, Risaliti & Lusso 2019), and the GRB sample (open black, Demianski et al. 2017a).

3. Fitting the Hubble diagram: a cosmographic approach

In this analysis, our goal is to use the Hubble diagram obtained by combining SNe Ia, quasars, and GRBs to test the concordance flat ΛCDM model through a model-independent cosmographic method. The cosmographic approach to cosmology is receiving increasing interest as it can be used to obtain a direct fit of the Hubble diagram without any hypothesis concerning the type of dark energy and dark matter required to satisfy the Einstein equation, and simply assuming the minimal priors of isotropy and homogeneity. Here we describe two cosmographic techniques we applied to the data in order to measure the cosmological parameters and test the concordance model.

3.1. Traditional cosmographic approach

To compute cosmological parameters, we followed the standard cosmographic procedure which is fully described in Sect. 2 of Demianski et al. (2017b). We briefly summarise the main steps. We parametrized the space-time geometry described by the Friedmann–Robertson–Lemaître–Walker (FLRW) metric by a Taylor series of the scale factor a(t) and its higher-order derivatives. This expansion leads to a distance–redshift relation that only relies on the assumption of the FLRW metric and is therefore fully model independent. From the cosmographic functions (Visser 2004 and Eqs. from (2) to (6) in Demianski et al. 2017b), we can compute the series expansion to the fifth order in time of the scale factor

$$\begin{array}{cc}\hfill \frac{a(t)}{a(t_0)}& =1+H_0(t-t_0)-\frac{q_0}{2}{H}_0^2{(t-t_0)}^2+\frac{j_0}{3!}{H}_0^3{(t-t_0)}^3\hfill \\ \hfill & +\frac{s_0}{4!}{H}_0^4{(t-t_0)}^4+\frac{l_0}{5!}{H}_0^5{(t-t_0)}^5+O[{(t-t_0)}^6],\hfill \end{array}$$(1)

where q₀, j₀, s₀, and l₀ are commonly indicated as the deceleration, jerk, snap, and lerk parameters, corresponding to the functions evaluated at the present time t₀ (Sahni et al. 2003). With this metric, it is possible to express all the observational quantities, like the luminosity distance d_L, as a power series in the redshift parameter z. From Eq. (2), we can construct the series for the luminosity distance as

$$\begin{array}{c}\hfill D_{\mathrm{L}}(z)=\frac{\mathit{cz}}{H_0}({\mathcal{D}}_{\mathrm{L}}^0+{\mathcal{D}}_{\mathrm{L}}^{1}z+{\mathrm{D}}_{\mathrm{L}}^2{z}^2+{\mathcal{D}}_{\mathrm{L}}^3{z}^3+{\mathcal{D}}_{\mathrm{L}}^4{z}^4+O(z^5)),\end{array}$$(2)

where

$$\begin{array}{cc}\hfill {\mathcal{D}}_{\mathrm{L}}^0& =1,\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill {\mathcal{D}}_{\mathrm{L}}^1& =-\frac{1}{2}(-1+q_0),\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill {\mathcal{D}}_{\mathrm{L}}^2& =-\frac{1}{6}(1-q_0-3q_0^2+j_0),\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill {\mathcal{D}}_{\mathrm{L}}^3& =\frac{1}{24}(2-2q_0-15q_0^2-15q_0^3+5j_0+10q_0{j}_0+s_0),\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill {\mathcal{D}}_{\mathrm{L}}^4& =\frac{1}{120}(-6+6q_0+81q_0^2+165q_0^3-105q_0^4-110q_0{j}_0\hfill \\ \hfill & -105q_0^2{j}_0-15q_0{s}_0-27j_0+10j_0^2-11s_0-l_0).\hfill \end{array}$$(7)

To mitigate the usual convergence problems of the cosmographic series expansion at high redshifts, we consider the expansion in the improved parameter y = z/(1 + z) (the so-called y-parameter), which has the great advantage of also holding for z >  1 and is therefore the appropriate tool for handling high-redshift sources such quasars and GRBs (Cattoën & Visser 2007; Vitagliano et al. 2010; Demianski et al. 2012, 2017b).

3.2. Cosmographic parameters from a logarithmic polynomial function

Following a procedure first described in RL19, we fitted the data with a fourth-grade polynomial of log(1 + z) and compare the best fit parameters with the expectations from the standard model. In particular, we define a cosmographic function of the luminosity distance as:

$$\begin{array}{c}\hfill \begin{array}{cc}& D_L(z)=kln(10){\displaystyle \frac{c}{H_0}}\times [log(1+z)+a_2{log}^2(1+z)\hfill \\ \hfill & +a_3{log}^3(1+z)+a_4{log}^4(1+z)]+O[{log}^5(1+z)],\hfill \end{array}\end{array}$$(8)

where k,  a₂,  a₃, and a₄ are free parameters. The parameter k is fitted separately for SNe Ia, quasars, and GRBs in order to cross-calibrate the samples. We note that the Hubble constant H₀ is degenerate with the k parameter, so we can fix it to an arbitrary value (H₀ = 70 km s⁻¹ Mpc⁻¹). The best-fit parameters and their uncertainties were obtained through the Python package emcee (Foreman-Mackey et al. 2013), which is a pure-Python implementation of Goodman & Weare’s affine invariant Markov chain Monte Carlo (MCMC) ensemble sampler.

Within the standard model, the only free parameter beyond H₀ is the matter density, Ω_M. Therefore, we can expand the function D_L(Ω_M, z) in a series of log(1 + z), and obtain the coefficients of each term as a function of Ω_M. This procedure is discussed in detail by RL19 for a third-order expansion, which is shown to reproduce a flat ΛCDM model up to z >  5, with only a few per cent difference at the upper end of the redshift range (see their Supplementary Fig. 5). Since the inclusion of GRBs requires a precise fit up to z ∼ 6.7, here we extend the expansion to the fourth order. In this way, our logarithmic expansion can match the standard model in the whole redshift range with a maximum difference below 2%. After simple algebra we obtain

$$\begin{array}{cc}& a_2=ln(10)\times ({\displaystyle \frac{3}{2}}-{\displaystyle \frac{3}{4}}{\mathrm{\Omega }}_M)\hfill \end{array}$$(9)

$$\begin{array}{cc}& a_3={ln}^2(10)\times ({\displaystyle \frac{9}{8}}{\mathrm{\Omega }}_M^2-2{\mathrm{\Omega }}_M+{\displaystyle \frac{7}{6}})\hfill \end{array}$$(10)

$$\begin{array}{cc}& a_4={ln}^3(10)\times (-{\displaystyle \frac{135}{64}}{\mathrm{\Omega }}_M^3+{\displaystyle \frac{18}{4}}{\mathrm{\Omega }}_M^2-{\displaystyle \frac{47}{16}}{\mathrm{\Omega }}_M+{\displaystyle \frac{5}{8}}).\hfill \end{array}$$(11)

The above equations describe a curve in the (a₂,  a₃,  a₄) parameter space that we can directly compare with the results of the MCMC technique applied to the SNe+quasars+GRBs Hubble diagram.

4. Statistical analysis

4.1. Cosmological parameters from the standard cosmographic approach

To provide statistical constraints on the cosmological parameters, we followed the same procedure as outlined in Demianski et al. (2012, 2017b). Specifically, we used the MCMC method and ran three parallel chains. We tested the convergence through the Gelman–Rubin diagnostic approach, which uses as a test probe the reduction factor R, i.e. the square root of the ratio of the variance between-chain and the variance within-chain. A large R indicates that the former is substantially greater than the latter, meaning that a longer simulation is needed. We required that R converge to 1 for each parameter. We set R − 1 of the order of 0.05, which is more restrictive than the often used and recommended value R − 1 <  0.1 for standard cosmological investigations. Moreover, since methods like the MCMC are based on an algorithm that moves randomly in the parameter space, we a priori imposed some basic consistency controls in the code to reduce the uncertainties of the cosmographic parameters, requiring that all of the (numerically) evaluated values of d_L(z) be positive. We discarded the first 30% of the point iterations at the beginning of each MCMC run, and thinned the chains that were run many times. We finally extracted the constraints on cosmographic parameters by coadding the thinned chains and inferred the median values and the confidence ranges from the marginalisation of the posterior probability distributions of the parameters.

The Hubble diagram resulting from the combination of SNe Ia (cyan points), quasars (yellow points), and GRBs (black points) is presented in Fig. 2. The dashed magenta line corresponds to a flat ΛCDM model fit with Ω_M = 0.3, whilst the green dotted line represents the regression fit (Eq. (2)) with the inclusion of GRBs and the most updated SNe Ia sample.

Fig. 2. Hubble diagram of SNe Ia (cyan points), quasars (yellow points), and GRBs (black points) with the 1σ uncertainties. Red points are the mean (also with 1σ uncertainties) of the distance modulus in narrow redshift bins for quasars only (shown for visualisation purposes). The dashed magenta line shows a flat ΛCDM model fit with ΩM = 0.3. The green dotted line is the best MCMC fit using the standard cosmographic approach (Sect. 3.1) with the inclusion of GRBs and Pantheon, whilst the blue dot-dashed line represents the new regression fit (Sect. 3.2). The black solid line is the best MCMC regression curve of the third-order expansion of log(1 + z) published by RL19 and obtained by fitting the quasars and the SNe Ia from the JLA survey (Betoule et al. 2014).

A summary of the results is provided in Table 1, whilst in Fig. 3 we present the resulting confidence regions for q₀ and j₀. Our statistical MCMC analysis shows that the jerk parameter, j₀, is significantly different from its ΛCDM value j₀ = 1 at more than 4σ.

Fig. 3. Confidence regions in the (q0 − j0) plane as provided by the standard cosmography. The value j0 = 1 (i.e. ΛCDM) is off at > 4σ, thus strengthening the results of RL19 using an independent approach.

4.2. Our new cosmographic approach

Figure 2 shows the results obtained with the new regression fit (Eq. (8)) considering quasars, GRBs and the Pantheon sample. A discrepancy between the ΛCDM model of low-redshift data and the model-independent cosmographic MCMC regression curve obtained with the full data set is observed, in agreement with the findings presented in Sect. 4.1.

For comparison, we also present with the black solid line the best MCMC regression curve of the third-order expansion in log(1 + z) published by RL19 and obtained by fitting the approximately 1600 quasars and the SNe Ia from the JLA survey (Betoule et al. 2014). Our new analysis that considers the updated SNe Ia sample and the GRBs further supports this tension with increased statistics.

In Fig. 4 we show the 1σ (dark colours) and 4σ (light colours) error contours in the (a₂,  a₃,  a₄) parameter space resulting from fitting the data using Eq. (8) for the samples Pantheon+quasars and Pantheon+quasars+GRBs, whilst we plot the 1σ and 2σ error contours for the Pantheon+GRBs. The black solid line represents a flat ΛCDM model with Ω_M in the interval 0.1–0.9, where the black point marks the value Ω_M = 0.3 (i.e. Eqs. (9)–(11)). For visualisation purposes, the orientation is chosen to better visualise the discrepancy between the contours and the flat ΛCDM model. From this figure, it is apparent that the tension of our cosmographic fit with the standard model is at 4σ for Pantheon+quasars, and remains significant even with the inclusion of the GRBs (bottom panel of Fig. 4). A summary of the results obtained from the MCMC regression analysis of the fourth order is provided in Table 2.

Fig. 4. 1σ (dark colours) and 4σ (light colours) error contours in the (a2,  a3,  a4) parameter space resulting from fitting the data using Eq. (8) for the samples Pantheon+quasars and Pantheon+quasars+GRBs. We plot the 1σ and 2σ error contours for the Pantheon+GRBs. The black solid line represents a flat ΛCDM model with ΩM in the interval 0.1–0.9, where the black point marks the value ΩM = 0.3 (i.e. Eqs. (9)–(11)). From top to bottom: result from the cosmographic technique described in Sects. 3.2 and 4.2.

As these results reflect an average deviation from the flat ΛCDM over the redshift probed by the different data sets, we can investigate what redshift range drives this tension.

Figure 5 shows how the deviation from a ΛCDM (with Ω_M = 0.3 and normalized to the Pantheon SNe sample up to z = 1) evolves with redshift for the three samples taken singularly. In the top panel we plot the average difference in distance modulus (DM) between the observed DM and the fit assuming a flat ΛCDM model (Δ(DM) = DM − DM_ΛCDM) in log z intervals, whilst the bottom panel represents the same deviation in σ units. Interestingly, all the three samples show that the deviation from the flat ΛCDM becomes increasingly statistically significant at redshifts higher than 1, with with 23 SNe Ia at z >  1 at more than 3σ tension with the flat ΛCDM model.

Fig. 5. Deviation of the distance modulus from the ΛCDM with ΩM = 0.3 and normalized to the Pantheon SNe sample up to z = 1 (Δ(DM) = DM−DMΛCDM) for SNe Ia (blue), quasars (yellow), and GRBs (red). Top panel: Points represent the averages of Δ(DM) in logarithmic redshift intervals. Bottom panel: Deviations in σ. We note that the statistical significance of the discrepancy from the ΛCDM in Fig. 4 is not directly comparable to the one plotted here since the MCMC fits in the former are computed over the entire redshift range.

As a final comment, we mention that the tension with the flat ΛCDM has a peak at z ≃ 3 (where we have excellent pointed observations from our XMM-Newton large program), but it seems to decrease at z >  3 where the data quality reduces. Pointed observations of quasars at z ≃ 4 are required to verify this trend.

5. Conclusions

The spatially flat ΛCDM cosmology is still the most commonly assumed representation of our Universe within the scientific community, yet there is mounting evidence that the standard ΛCDM model (or its simplest extension, e.g. wCDM) is in tension at more than 4σ with local direct measurements of the Hubble constant (e.g. Riess et al. 2019 and references therein).

Recently, RL19 also reported on a ∼4σ tension with the ΛCDM model through a model-independent parametrization of a Hubble Diagram of SNe Ia and quasars. Since this tension could be subject to some level of systematic error, it is paramount to test the ΛCDM model by making use of model-independent approaches and independent samples over a wide redshift range.

To this aim, we present a joint statistical analysis of the Hubble diagram from SNe Ia, quasars, and GRBs, thus using three completely independent samples of standardisable candles spanning a wide redshift interval, from local objects up to z ≃ 6.7. The statistical analysis was performed in two ways: (1) assuming a traditional cosmographic approach and (2) through a cosmographic fourth-order polynomial fit in log(1 + z) of the Hubble diagram, following RL19.

From the standard approach, we found that j₀ is significantly different from its ΛCDM value, j₀ = 1, at more than 4σ.

Pursuing the cosmographic approach, we again confirmed the tension between the best-fit cosmographic parameters and the ΛCDM model with the joint SNe Ia, quasar, and GRB data set at more than the 4σ statistical level. We also performed the same analysis by considering the SNe Ia+quasar and the SNe Ia+GRB samples separately, finding that the tension is still statistically significant at the ∼4σ and ∼2σ level, respectively. Moreover, we also confirm that this tension becomes statistically significant (above 3σ) only at high redshifts (z >  1) for SNe Ia and quasars taken independently, whilst this is at ∼2σ for GRBs alone.

Summarising, these two complementary (and independent) statistical analyses both confirm a tension with the flat ΛCDM model at 4σ emerging at high redshifts. Moreover, as the completely independent high-redshift Hubble diagrams of quasars and GRBs are fully consistent with each other, this strongly suggests that the deviation from the standard model is not due to unknown systematic effects but to new physics. Possible extensions of the standard cosmological model will be explored in future works.

Acknowledgments

We acknowledge financial contribution from the agreement ASI-INAF n.2017-14-H.O. SB is supported by NASA through the Chandra award no. AR7-18013X (NAS8-03060) and by the grant HST-AR-13240.009. EN acknowledges funding from the EU Horizon 2020 Marie Skłodowska-Curie grant no. 664931.

References

All Tables

All Figures

	Fig. 1. Redshift distribution of the Pantheon SNe Ia sample (green, Scolnic et al. 2018), the quasars sample (yellow, Risaliti & Lusso 2019), and the GRB sample (open black, Demianski et al. 2017a).
In the text

Fig. 2. Hubble diagram of SNe Ia (cyan points), quasars (yellow points), and GRBs (black points) with the 1σ uncertainties. Red points are the mean (also with 1σ uncertainties) of the distance modulus in narrow redshift bins for quasars only (shown for visualisation purposes). The dashed magenta line shows a flat ΛCDM model fit with ΩM = 0.3. The green dotted line is the best MCMC fit using the standard cosmographic approach (Sect. 3.1) with the inclusion of GRBs and Pantheon, whilst the blue dot-dashed line represents the new regression fit (Sect. 3.2). The black solid line is the best MCMC regression curve of the third-order expansion of log(1 + z) published by RL19 and obtained by fitting the quasars and the SNe Ia from the JLA survey (Betoule et al. 2014).

In the text

	Fig. 3. Confidence regions in the (q0 − j0) plane as provided by the standard cosmography. The value j0 = 1 (i.e. ΛCDM) is off at > 4σ, thus strengthening the results of RL19 using an independent approach.
In the text

Fig. 4. 1σ (dark colours) and 4σ (light colours) error contours in the (a2,  a3,  a4) parameter space resulting from fitting the data using Eq. (8) for the samples Pantheon+quasars and Pantheon+quasars+GRBs. We plot the 1σ and 2σ error contours for the Pantheon+GRBs. The black solid line represents a flat ΛCDM model with ΩM in the interval 0.1–0.9, where the black point marks the value ΩM = 0.3 (i.e. Eqs. (9)–(11)). From top to bottom: result from the cosmographic technique described in Sects. 3.2 and 4.2.

In the text

Fig. 5. Deviation of the distance modulus from the ΛCDM with ΩM = 0.3 and normalized to the Pantheon SNe sample up to z = 1 (Δ(DM) = DM−DMΛCDM) for SNe Ia (blue), quasars (yellow), and GRBs (red). Top panel: Points represent the averages of Δ(DM) in logarithmic redshift intervals. Bottom panel: Deviations in σ. We note that the statistical significance of the discrepancy from the ΛCDM in Fig. 4 is not directly comparable to the one plotted here since the MCMC fits in the former are computed over the entire redshift range.

In the text