© ESO, 2016

1. Introduction

The best current data, ranging from cosmological to astrophysical scales, have established that the Universe is on average highly homogeneous and isotropic and that it is expanding in an accelerating fashion.

The observation that, in comparison with nearby sources, distance sources (z> 0.3) appear dimmer than predictable in a universe with matter only was interpreted as the need for a dark energy component in the luminosity distance (Riess et al. 1998; Permutter et al. 1999). This dimming, however, could be caused by a variation of any of the other components that affect the luminosity distance, namely by inhomogeneities in the energy densities or anisotropies in the expansion (Amendola et al. 2013). On large scales, the matter distribution is statistically homogeneous. Locally, however, matter is distributed according to a pattern of alternate overdensed regions (filaments) and underdensed regions (voids; Scrimgeour et al. 2012; Watson et al. 2013), which is thought to reflect quantum fluctuations away from homogeneity that originated in the very early Universe. In this paper, we hypothesize that the resulting inhomogeneity in the matter distribution might be imprinted in an angular distribution of the cosmological parameters, which is measured at the scale of the matter distribution pattern, about the average parameter values.

Inhomogeneities in the cosmological parameters distribution might create a dynamical mechanism that would emulate a cosmic acceleration as a derived, instead of fundamental, effect (see e.g. Kolb et al. 2006; Larena et al. 2009; Clarkson & Umeh 2011; Chiesa et al. 2014; Adamek et al. 2015). Conversely, inhomogeneities in the cosmological parameters might reflect new physics, such as a dynamical dark energy component or a modification to the Einstein’s theory of general relativity (see e.g. Amendola et al. 2013). On the one hand, inhomogeneities in the matter distribution might generate fluctuations in the dynamical dark energy, which would imply angular variations of the Hubble parameter across the sky that could be observed as anisotropies of the luminosity distance as a function of redshift d_L(z); on the other hand, fluctuations in the new gravitational fields could have imprints on the local expansion and, consequently, on d_L(z) (Cooray et al. 2010).

Supernovae establish the relation of the luminosity distance with redshift up to z ≈ 1. The luminosity distance is sensitive to both the geometry and growth of structure (Bonvin et al. 2006), which renders supernovae a sensitive probe of the late-time expansion history of the Universe. The current amount of data on SNe from the combination of different SN surveys, if large and homogeneous enough, could be used to investigate the angular distribution of inhomogeneities and hence to investigate the nature of cosmic acceleration.

Supernovae have previously been used to constrain inhomogeneities in the cosmic expansion, namely by measuring the hemispherical anisotropy of H₀ (Kalus et al. 2013), assessing the dependence of H₀ with the position on the cosmic web (Wotjak et al. 2014), mapping {H₀,q₀} in the sky (Bengaly et al. 2015)¹, and by estimating the statistical significance of rejecting an isotropic magnitude-redshift relation (Javanmardi et al. 2015)². Supernovae have also been used to constrain inhomogeneities in the dark energy, namely by constraining fluctuations of dark energy (Blomqvist et al. 2010; Cooray et al. 2010) and by constraining radial inhomogeneity for a Lemaître-Tolman-Bondi metric with a cosmological constant (Valkenburg et al. 2014).

In this paper, we develop a method to probe inhomogeneities in the form of a local estimation of cosmological parameters. We demonstrate the method on type Ia supernovae, in particular by dividing the supernova sample into different regions of the sky and estimating the angular distribution of the cosmological parameters that affect the luminosity distance, namely {Ω_M,Ω_Λ,H₀}, across the different regions. Each region is assumed to be described by a Friedmann-Lemaître-Roberston-Walker metric so that the full sky is an inhomogeneous ensemble of disjoint, locally homogeneous regions. This ensemble of homogeneous regions is an approximate solution to the Einstein equation, which we assume to be a sufficient approximation to enable us to interpret the cosmological parameters {Ω_M,Ω_Λ,H₀} as a description of the true energy content. We also suggest measuring the anisotropy in the angular distribution of the parameters not by looking for particular asymmetries at particular scales, but by computing the angular power spectrum that includes all multipoles accessible in the regions of the sky. We then use the angular fluctuations in H₀ to estimate the backreaction in the sense of Buchert.

Since the Einstein tensor is non-linear in the metric, the average evolution of an inhomogeneous space-time is not the same as the evolution of a smooth space-time with the same initial conditions (Buchert 2000). The effect of inhomogeneities on the average properties is known as backreaction. When averaging over the homogeneous regions, angular fluctuations in the expansion factor, and consequently in H₀, induce a backreaction term in the form of an extra positive acceleration (Räsänen 2006). We derive the analytical extra positive acceleration for a toy model of an arbitrary number of disjoint, homogeneous regions and compute the overall deceleration parameter assuming (a) no backreaction and (b) backreaction for the measured angular distribution of H₀.

This paper is organized as follows. In Sect. 2 we describe the SN data and the theory for the calculation of the observables. In Sect. 3 we introduce the method by performing both a global and a local estimation, according to the parameter estimation described in Appendix A. From the global estimation, there result parameter values taken for fiducial values. From the local estimation, there result maps of the cosmological parameters. In Sect. 4 we include three tests of the inhomogeneity of the SN surveys. The first test is in the SN distribution in redshift, the second test is in the SN distribution in both redshift and sampling, and the third test is in the removal of a noise bias hypothesized as a measure of the inhomogeneity of the SN sampling. We then proceed to analyse the results. In Sect. 5 we compute the power spectrum of the maps of the parameters. In Sect. 6, for a toy model of backreaction, we compute the values of the parameters averaged over the pixels. These are, to the best of our knowledge, the first estimates of inhomogeneities from the angular distribution of cosmological parameters and the first estimates of backreaction effects. In Sect. 7 we discuss current limitations of the method owing to insufficient data and we indicate future research directions.

2. Supernova data and theory

We use the type Ia supernova sample compiled by the Joint Light-curve Analysis (JLA) SNLS-SDSS collaborative effort (Betoule et al. 2014). The JLA collaboration aimed to include SNe from surveys at different redshifts, improve the accuracy of the calibration, and make the light-curve model uniform across the surveys. The combined sample consists of type Ia SNe from different supernova surveys, namely the Sloan Digital Sky Survey II (SDSS-II) at intermediate redshifts (Sako et al. 2016) and the Supernova Legacy Survey (SNLS) combined with Hubble Space Telescope (HST) surveys at high redshifts and other surveys at low redshifts (Conley et al. 2011), totalling N_SNe = 740 SNe with redshift 0.01 ≤ z ≤ 1.30. The sample includes the SNe in the Hubble flow (z ≥ 0.01) that are not affected by strong Milky Way contamination³. The SNe are distributed on the sky according to Fig. 1.

From the sampled data, we compute the distance modulus for each SN as (1)For each SN we take the peak magnitude m_B inferred from the observed magnitude and corrected for the selection bias, and we define the effective absolute magnitude M_B,eff as the absolute magnitude M_B corrected by the time stretching of the light curve x₁ and the supernova colour at maximum brightness c as follows: (2)The parameters M_B = −19.05 ± 0.02,α = 0.141 ± 0.006, and β = 3.101 ± 0.075 are estimated for all SNe, whereas the parameters { m_B,x₁,c } are estimated for each SN separately by fitting a model of the Ia SN spectral sequence to the photometric data. We compute μ_data for all SNe in the sample and plot it in the left panel of Fig. 2. The error σ_μdata results from the error propagation according to Eqs. (1) and (2) assuming uncorrelated errors.

We now establish the dependence of the distance modulus with the cosmological parameters. The distance modulus is defined in terms of the luminosity distance d_L (measured in Mpc) as (3)For a (3+1)-dimensional line element expressed in spherical coordinates as (4)where κ = { 1,0,−1 } characterizes the spatial curvature and d stands for the two-sphere line element, the luminosity distance is given by (5)where (6)and (7)with the dimensionless energy density of matter, curvature, and dark energy at present time given, respectively, by (8)Assuming that (9)then Hence the distance modulus via the luminosity distance depends on the following set of cosmological parameters { Ω_M,Ω_Λ,H₀ } subject to the condition in Eq. (9).

3. Method

We introduce the method by performing the estimation of the parameters { Ω_M,Ω_Λ,H₀ } as described in Appendix A, first from the complete SN sample and then from SN subsamples grouped by regions in the sky.

3.1. Global parameter estimation

We first use the complete SN sample to fit μ_theo to μ_data and estimate the maximum-likelihood values for the parameters x_i = { Ω_M,Ω_Λ,H₀ } . Whenever unspecified, H₀ is measured in units of km s^-1 Mpc^-1. We ran ten realizations of Markov chains of length 10⁴. The starting point of each realization is randomly generated. We remove the first 20% of entries in the chains to keep the burnt-in phase only. We then thin down the remaining 80% of entries to a half by removing one of each consecutive entry to eliminate correlations within each chain. Averaging over the various realizations of Markov chains, we obtain the following results: If we moreover impose Ω_M + Ω_Λ = 1, thus enforcing Ω_κ = 0, then we obtain the following results: (See Table 1.) These are the mean values and that we use as fiducial values in the subsequent calculations. The errors contain the dispersion in each chain and the dispersion among the averages of the different chains added in quadrature. By fixing H₀ = 70.0 km s^-1 Mpc^-1 and imposing Ω_M + Ω_Λ = 1, we obtain Ω_{M,flat,H0 = 70.0} = 0.322 ± 0.012.

We recall that the JLA collaboration found Ω_M,JLA = 0.295 ± 0.034 for a fixed Hubble parameter H₀ = 70.0 km s^-1 Mpc^-1, which includes our result Ω_{M,flat,H0 = 70.0} within the errors. The results obtained by the JLA collaboration in combination with other data were { Ω_M,Ω_Λ,H₀ } _JLAComb = { 0.305,0.697,68.3 } ± { 0.010,0.010,1.03 } (Betoule et al. 2014). For comparison, in Table 1 we also include the values for these parameters obtained from the ninth-year WMAP data (Hinshaw et al. 2013) and the Planck collaboration (Planck Collaboration XVI 2014). We observe that our results are compatible with the WMAP results, but favour lower values for Ω_M, and simultaneously higher values for Ω_Λ and for H₀, in comparison to the JLAComb and Planck results. In contrast, Riess et al. found H₀ = 73.8 ± 2.4 km s^-1 Mpc^-1 (Riess et al. 2011).

A derived parameter of interest is the curvature energy density Ω_κ, which becomes a dependent variable owing to the flatness condition in Eq. (9). The fiducial values return and The values for Ω_κ returned by the above-mentioned studies are included in Table 1. We observe that our results are compatible with the WMAP result but favour values for Ω_κ that are higher than the JLA Comb results and lower than the Planck results.

Another derived parameter of interest is the cosmic deceleration parameter q₀, which is defined in terms of { Ω_M,Ω_Λ } as (12)The fiducial values return and The values for q₀ returned by the above-mentioned studies are included in Table 1. We observe that our results are compatible with the WMAP result but favour lower values for q₀ in comparison to the JLAComb and Planck results.

We compute μ_theo for the estimated values and plot it as a black line in the left panel of Fig. 2. In the right panel of Fig. 2, we plot the Markov chain of one realization. We observe that (a) the points are aligned approximately parallel to the lines of constant H₀ and (b) the points form approximately an ellipse whose largest axis is not parallel to the lines of constant q₀.

To assess the convergence, we apply the Gelman-Rubin heuristic convergence test. For each parameter x_i and N_chain chains, we compute the mean x_ij and variance of each chain j. We take the mean of all chains and compute the variance of the means of the chains We compute the average of the variances of the chains and compare with forming the ratio (13)If the chains are well mixed and have all sampled the parameter space well, then r_GR ~ 1. We find r_GR ≈ 0.90 for the estimation of both and

3.2. Local parameter estimation: cosmic variance

We proceed to divide the SNe into subsamples grouped by regions in the sky. From each subsample, we aim to estimate the parameters x_i = { Ω_M,Ω_Λ,H₀ } and produce maps such that the values in each region are given by the values of the parameters estimated with the subsample in that region. Hence we must guarantee that the regions have the same surface area and all contain SNe. Ideally each pixel should also have a SN subsample that covers a range of z enough to break the degeneracies among the parameters in Eq. (11). We use the HEALPix pixelation scheme to divide the sphere into pixels of equal surface area (Gorski et al. 2005), hence defining the regions as pixels. In the HEALPix pixelation scheme, we select nside= 1, which produces a map divided into npixel= 12 pixels. The number of SNe in each pixel is indicated in Fig. 3. In this description, the complete sample is equivalent to npixel= 1 pixel.

We then use the SN subsample in each pixel to estimate the maximum-likelihood values for the parameters { Ω_M,Ω_Λ,H₀ } . This estimation is performed, as in the case of the complete sample, on realizations of the original SN distribution in redshift and in space, i.e. for the original dependence of redshift with the position and for the original sampling per pixel. The variance of the estimation in each pixel determines the accuracy with which we can measure the parameters with this sample. For this reason we call this estimation “Cosmic variance”. In each pixel, we ran 30 realizations of Markov chains of length 10⁴. In each pixel, the starting point of each realization is randomly generated as for the case of the complete sample. Also, as in the complete sample, we perform the burnt-in selection and the thinning reduction of the Markov chains in each pixel. Hence from each chain j, at each pixel k, we obtain an estimate for each parameter x_i. The parameter estimation results in maps x_ijk = { Ω_M,Ω_Λ,H₀ } _jk. We average the maps of each parameter over the chains by taking the mean weighted by the number of SNe in each pixel, Since the number of SNe in each pixel is the same for the different chains, this averaging is equivalent to the simple mean of the maps by pixel. We denote the resulting mean maps of the parameters by To assess the convergence, we compute the Gelman–Rubin ratio, finding r_GR ≈ 0.90 for the estimation of all parameters at all pixels.

For better visualization, we compute difference maps by subtracting the fiducial values off the corresponding maps { Ω_M,Ω_Λ,H₀ } _{k,cosmic_var}. Averaging over the various realizations of Markov chains at each pixel, we obtain the difference maps shown in Fig. 4. Similarly, using the maps obtained in the “Cosmic variance” estimation, we compute the maps for Ω_κ and q₀. The difference maps ΔΩ_κ and Δq₀, obtained by subtracting and , respectively, from the corresponding maps, are shown Fig. 4. Comparing the difference maps with the fiducial values, we measure fluctuations of order 10–100% for Ω_M, 1–30% for Ω_Λ, and up to 3% for H₀. Similarly, we measure fluctuations of order 0.1–1000% for Ω_κ and of order 10–60% for q₀. Although the fluctuations of Ω_κ are mostly small, they assume very high values in the pixels where { Ω_M,Ω_Λ } fluctuate in the same direction. From the difference maps, we observe that the fluctuations in the Ω_M map are predominantly towards values higher than whereas the fluctuations in the Ω_Λ map are predominantly towards values lower than The fluctuations in the H₀ map are equally distributed from lower to higher values than the fiducial.

To compare our results with previous studies, we define the largest fluctuation in the difference maps as (14)which yields max( { ΔΩ_M,ΔΩ_Λ,ΔH₀ } ) = { 0.352,0.286,3.84 } ,max(ΔΩ_κ) = 0.363, and max(Δq₀) = 0.452. We recall the results obtained in other SN studies. Kalus et al. found 2.0 <max(ΔH₀)<3.4 km s^-1 Mpc^-1 from the Union 2 data for fixed q₀ = −0.601 (Kalus et al. 2013). Bengaly et al. found max(ΔH₀) = { 4.6,5.6 } km s^-1 Mpc^-1 and max(Δq₀) = { 2.56,1.62 } from the Union 2.1 data and JLA data, respectively (Bengaly et al. 2015). Our results for H₀ are intermediate between the previous results, whereas our results for q₀ are more stringent than those of Bengaly et al. Our results for { Ω_M,Ω_Λ } , and subsequently for Ω_κ are, to the best of our knowledge, the first measurements of this kind.

We also observe a mild negative correlation between the number of SNe in each pixel and the fluctuations about the fiducial values, which hints at a poor sampling derived from an inhomogeneous coverage of the sky, both in angular and redshift space, by the SN surveys. This is further supported by the negative correlation between the number of SNe and the standard deviation of the parameters estimated in each pixel. A poor sampling into pixels implies a poor sampling in redshift, which limits the constraint of the degeneracies among the parameters and hence compromises the estimation of the parameters. Figure 5 illustrates these observations. In order to quantify this correlation, we define the correlation between the maps of the parameters and the number of SNe per pixel N_SNk by which takes values in [−1,1] . Similarly, we define the correlation between the sigma maps of the parameters and the number of SNe per pixel N_SNk. We find N_SNk] = {−0.394,0.397,−0.113 } and N_SNk] = {−0.408,−0.452,−0.328 } . The values of show that the angular variations of H₀ are less sensitive to the number of SNe per pixel than the angular variations of Ω_M or Ω_Λ. The comparison of the values of with the values of shows that the angular variations of the parameters are not as correlated with the number of SNe per pixel as the corresponding standard deviations.

The Galactic plane could contribute a spurious inhomogeneity in the local parameter estimation. However, an inspection of the pixels through which the Galactic plane crosses does not consistently show higher or lower values of the fluctuations to suggest a correlation.

These results are based on the assumption that both μ_data(z) and σ_μdata are independent. The existence of correlations would effectively reduce the number of independent SNe, which would aggravate the parameter estimation in the less SN-populated pixels. Assuming independent errors as a first approximation and having observed that the results are highly affected by the inhomogeneity of the SN surveys, the consideration of correlated errors would be a higher order effect and, hence, safely negligible for the current SN sample.

To assess how the local estimation can be improved, we reduce the number of free parameters by fixing one parameter at a time to the corresponding fiducial value estimated from the complete sample. We observe that the estimation of the Hubble parameter is largely unaffected by the size of the pixel subsample. Conversely the estimation of the energy densities depends highly on the size of the pixel subsample.

4. Inhomogeneity tests

In order to further establish the validity of the method, we perform the following tests of the inhomogeneity of the SN surveys. In the first and second test, following Kalus et al. (2013), we randomize the dependence of the redshift with the position by shuffling the SNe in z. In the first test, we keep the inhomogeneous spatial distribution of the SNe (i.e. we keep the positions of the SNe as in the JLA sample) and, in the second test, we simulate a homogeneous spatial distribution (i.e. we choose random positions for the SNe). In the third test, we hypothesize a noise bias as a measure of the inhomogeneity in the SNe sampling and remove it from the local “Cosmic variance” estimation.

4.1. Inhomogeneity test: shuffle the SNe

In this test, we randomize the dependence of redshift with the position for the original sample by shuffling the SNe in z while keeping the positions in the sky. This way, we average out the inhomogeneities in the SN distribution in redshift for the same subsampling per pixel.

Similar to the “Cosmic variance” estimation, we use the SN subsample in each pixel to estimate the maximum-likelihood values for the parameters { Ω_M,Ω_Λ,H₀ } . We call this estimation “Shuffle SNe”. In each pixel, we ran 100 realizations of Markov chains of length 10⁴. For each chain, the SN redshift is randomly shuffled for the same positions. In each pixel, the starting point of each realization is randomly generated as for the case of the complete sample. Also, as in the complete sample, we perform the burnt-in selection and thinning reduction of the Markov chains in each pixel. From each chain, the parameter estimation results in maps x_ijk = { Ω_M,Ω_Λ,H₀ } _{jk,shuffle_sn}. We average the maps for each parameter over the chains. Since the number of SNe in each pixel is the same for the different chains, we average the maps by taking the simple mean by pixel. We denote the resulting mean maps of the parameters by

For better visualization, we compute difference maps by subtracting the fiducial values from the corresponding maps. Averaging over the various realizations of Markov chains at each pixel, we obtain the difference maps shown in the left column of Fig. 6. For an easier comparison, we use the same colour scale for the maps that we used in the “Cosmic variance” estimation.

Comparing the difference maps with we measure fluctuations of order 1–60% for Ω_M, 1–20% for Ω_Λ, and up to 3% for H₀. We observe that the fluctuations in the Ω_M are predominantly towards values higher than whereas the fluctuations in the Ω_Λ and H₀ are predominantly towards values lower than the corresponding fiducial values. We also observe the same negative correlation between the number of SNe in each pixel and the fluctuations about the fiducial values as in the “Cosmic variance” estimation.

These fluctuations measure the inhomogeneity in the SN sampling, as observed in Fig. 1, which is reflected in the inhomogeneity of the subsampling into pixels. As observed, increasing the number of realizations increasingly homogenizes the sampling in each pixel and hence reduces the fluctuations about the fiducial values in pixels with large subsamples of SNe; however, increasing the number of realizations does not affect the estimation in pixels with small subsamples of SNe because the cosmic variance is large. For a sufficiently large number of realizations, such that the average is consistent with a homogeneous sampling into pixels with large subsamples of SNe, we hypothesize that this test measures the noise bias in the parameter estimation because of the inhomogeneity of the SN sampling. This noise bias is measured by the average maps produced in this test. Hence by removing this noise bias from the “Cosmic variance” estimation, we would obtain a parameter estimation that is independent of sampling inhomogeneities albeit with increased error.

4.2. Inhomogeneity test: randomize the pixelation

In this test, we randomize the dependence of the redshift with the position for a homogeneous sample per pixel by distributing the SNe in randomly chosen positions in the sky. This way we are averaging out inhomogeneities in the SN distribution in redshift and we simultaneously average out inhomogeneities in the sampling.

Similar to the “Cosmic variance” estimation, we use the SN subsample in each pixel to estimate the maximum-likelihood values for the parameters { Ω_M,Ω_Λ,H₀ } . We call this estimation “Random pixelation”. In each pixel, we ran 100 realizations of Markov chains of length 10⁴. For each chain, the SNe redshift is randomly distributed in a homogeneous way throughout the sky. In each pixel, the starting point of each realization is randomly generated as for the case of the complete sample. Also, as in the complete sample, we perform the burnt-in selection and thinning reduction of the Markov chains in each pixel. From each chain, the parameter estimation results in maps x_ijk = { Ω_M,Ω_Λ,H₀ } _{jk,random_pix}. We average the maps for each parameter over the chains. Since the number of SNe in each pixel varies among the chains, we average the maps by taking the mean weighted by the number of SNe in each pixel. We denote the resulting mean maps of the parameters by .

For better visualization, we compute difference maps by subtracting the fiducial values off the corresponding maps. Averaging over the various realizations of Markov chains at each pixel, we obtain the difference maps { ΔΩ_M,ΔΩ_Λ,ΔH₀ } _{k,random_pix} shown in the right column of Fig. 6.

Comparing the difference maps with we measure fluctuations of order 10% for Ω_M and Ω_Λ, and of order less than 1% for H₀. These fluctuations measure the limitation of an equal subsampling per pixel of the original sample. Since these fluctuations are smaller than the error of the global parameter estimation, they are consistent with a homogeneous sampling in all pixels, which implies that the original sample is large enough to be divided into subsamples.

4.3. Inhomogeneity test: remove the noise bias

We now proceed to test the hypothesis suggested in the context of the “Shuffle SNe” estimation. We let x_ijk and denote maps produced by j chains of the “Cosmic variance” estimation and j′ chains of the “Shuffle SNe” estimation, respectively. Then for each chain j, we define the unbiased maps as where is the average over the j′ realizations of the “Shuffle SNe” estimation. The unbiased maps averaged over the chains j are and the corresponding unbiased difference maps are We show in Fig. 7 to be compared with Fig. 4. Comparing the unbiased difference maps with we measure fluctuations of order 5–95% for Ω_M, 1–25% for Ω_Λ, and up to 5% for H₀. Similarly, we measure fluctuations of order 35–900% for Ω_κ and 1–50% for q₀.

For the correlations of the number of SNe per pixels with the difference maps and sigma maps, we find N_SNk] = { −0.117,0.162,−0.343 } and , respectively. For the unbiased maps in comparison with the original maps, the correlation with the difference maps decreases, whereas the correlation with the error maps increases, as expected. These results seem to indicate the validity of as a measure of the noise bias from the inhomogeneity of the SN sampling. Figure 8 illustrates these observations to be compared with Fig. 5.

5. Angular power spectrum of the parameters

In order to probe anisotropies in the parameters, we compute the angular power spectrum of the difference maps that result from each realization normalized to the corresponding fiducial parameter. We consider the cosmic variance estimation. For each estimated parameter i, for each realization j, we define the map of the fluctuations about the fiducial value at each pixel k as The map defined over the full sky can be decomposed into spherical harmonics according to (18)where each pixel k defines a unit direction vector n_k on the sky. The harmonic coefficients are given by (19)which for a pixelated map with N_pix pixels can be estimated by (20)Assuming statistical isotropy, the â_ℓm,ij can be used to estimate the power spectrum Ĉ_ℓ,ij according to (21)which can be inverted to yield the power spectrum of the parameter i from the map realization j(22)The variable ℓ can be written approximately as ℓ = 2π/θ, where θ ranges from the angular size of the map to twice the size of the pixel, which corresponds to the Nyquist frequency. The minimum value of ℓ is determined by the largest scale probed by the map, hence θ = 2π and ℓ_min = 1. Conversely the maximum value of ℓ is determined by the smallest scale probed by the map, hence θ = 2 × (pixelsize) and ℓ_max = 3.

Averaging the power spectra over the different realizations, we obtain the mean power spectrum Ĉ_ℓ,i for the parameters { Ω_M,Ω_Λ,H₀ } _{cosmic_var}(23)The variance of the estimator is given by (24)Hence the total variance of Ĉ_ℓ,i is the sum of the sample variance about the mean value and the variance of the estimator (25)The resulting power spectra are plotted in Fig. 9 as dashed lines.

For the density parameters, as well as for the deceleration parameter, we observe that the power spectrum has a maximum at the quadrupole (ℓ = 2). The Ω_M power spectrum is significantly larger than the Ω_Λ power spectrum over all scales, ranging from 5 to 20 times larger. For the Hubble parameter, the power spectrum always increases with the multipole. The quadrupole is approximately { 2.0,1.6,5.4 } times the dipole (ℓ = 1) in the { Ω_M,Ω_Λ,H₀ } _{cosmic_var} power spectra, respectively.

In Bengaly et al. (2015), the dipole was the dominant scale of the H₀ power spectrum produced from the low-redshift (z ≤ 0.20) SNe of the JLA sample. This effect was attributed either to the bulk motion (i.e. the recessional velocity component due to the global expansion) or to the inhomogeneous distribution of the data set.

In order to establish the validity of the method for the power spectrum analysis, we use the average “Shuffle SNe” estimation over j′ realizations, denoted by as a measure of the noise bias due to the inhomogeneity of the SN sampling. For each estimated parameter i, we define the average map of the fluctuations about the fiducial value at each pixel k as From the decomposition into spherical harmonics of we compute We define the unbiased power spectrum as (26)with total variance (27)The resulting power spectra are plotted in Fig. 9 with the bias power spectra as dotted lines and the unbiased power spectra as solid lines. The unbiased power spectra of the density parameters and of the Hubble parameter have the same behaviour, albeit less pronounced, as the corresponding power spectra before the noise bias removal. The exception is the deceleration parameter whose subtle maximum at ℓ = 2 is erased with the noise bias removal. After the noise bias removal, the quadrupole is approximately { 1.3,1.4,5.1 } times the dipole in the { Ω_M,Ω_Λ,H₀ } _{cosmic_var} power spectra, respectively. Although these results seem to indicate the validity of as a measure of the noise bias and the dominance of the quadrupole as the presence an inhomogeneity scale, the errors are too large to validate the measurements, thus rendering the results compatible with a flat spectrum.

6. Average over the local estimations

We proceed to compute average values of the parameters from the corresponding locally estimated values. We use the maps obtained in the “Cosmic variance” estimation. Given the mean maps over the chains , we want to average the map of each parameter over the pixel subsamples.

One way of averaging consists in taking the mean weighted by the inverse of the variance of the parameter in each pixel k,(28)For the estimated parameters, we find 0.607,71.06 } ± { 0.099,0.115,0.87 } and 0.724,70.93 } ± { 0.041,0.041,0.64 } . These values are consistent with, but systematically smaller than, the fiducial values estimated from the complete sample (see Table 2). After the noise bias removal, we find { 0.115,0.120,0.94 } and 71.17 } ± { 0.056,0.056,0.90 } . These values are consistent with the fiducial values. We also observe that the flatness condition yields a lower H₀ than in the absence of the flatness condition, which is the opposite trend of that observed in the estimation from the complete sample.

Similarly, we compute the pixel average of the maps for Ω_κ and q₀ derived from the maps obtained in the “Cosmic variance” estimation. The average value of Ω_κ over the pixels is and and after the noise bias removal it is and 0.029 ± 0.079. The average value of q₀ over the pixels is and and after the noise bias removal it is and (see Table 2). After the noise bias removal, the pixel average values of both Ω_κ and q₀ become closer to the corresponding values estimated from the complete sample. We also observe that the flatness condition yields a smaller q₀.

Since we are averaging over inhomogeneous maps, we consider a toy model of an inhomogeneous space-time consisting of disjoint, locally homogeneous regions with different expansion factors. Hence another way of averaging consists in taking the mean weighted by the 3-volume V_k of each pixel k,(29)This averaging requires a choice of time slice. Since all pixels have the same surface area and hence the angular directions expand the same way today, we can assume that the radial direction also expands the same way today, which implies that the three-volume is the same today for all pixels. This amounts to identifying a volume as a pixel. Defining v_k = V_k/ ∑ _kV_k and generalizing the result in Räsänen (2006) for N_pixel disjoint regions, the average of the Hubble parameter is given by (30)Taking the derivative, it follows that (31)which decomposes into a linear term in the pixel average of and a quadratic term in differences of H₀ between pairs of pixels. The quadratic (backreaction) term generates an acceleration that is due not to regions speeding up locally, but instead to the slower regions becoming less represented in the average. In the absence of the quadratic term, the volume average reduces to the pixel average above. Then the volume average of q₀ becomes Using the fluctuations in H₀ measured in the local parameter estimation, we find and and after the noise bias removal we find and (see Table 2). After the noise bias removal, the volume average values become closer to the value for the deceleration estimated from the complete sample. The results from both averaging methods strengthen the validity of as a measure of the noise bias due to the inhomogeneity of the SN sampling.

The quadratic term is of order 10^-3 (or lower) times the linear term, which means that the error of the difference is below the standard deviation. Hence for the angular fluctuations in H₀ measured with this SN sample, the contribution of the quadratic term in Eq. (31) is insignificant, which renders the volume averaging equivalent to the pixel averaging. Hence in the context of this toy model of an inhomogeneous space-time, backreaction is not a viable dynamical mechanism to emulate cosmic acceleration.

7. Conclusions

In this paper, we developed a method to probe the inhomogeneity of the large-scale structure and cosmic acceleration and applied to type Ia SN data. In particular, we subsampled the original SN sample into pixels of equal surface area and in each pixel we estimated the cosmological parameters that affect the luminosity distance, namely { Ω_M,Ω_Λ,H₀ } , assuming a local Friedmann-Lemaître-Roberston-Walker metric in each pixel. This local parameter estimation resulted in maps of the cosmological parameters with the same pixelation as the SN subsamples.

We measured fluctuations about the average values of order 5–95% for the matter energy density Ω_M, 1–25% for the dark energy density Ω_Λ, and to 5% for the Hubble parameter H₀, which translate into fluctuations about the average values of order 35–900% for the curvature energy density Ω_κ and 1–50% for the deceleration parameter q₀. The comparison of our results with those obtained in other SN studies indicates convergence of the results and consistency among the measurements. Owing to the spatial inhomogeneity of the SN surveys, as demonstrated by the tests performed, the pixel subsamples are statistically unequal. As a result, well-sampled pixels return robust estimations, whereas poorly sampled pixels return weak constraints. We also assessed how the local estimation can be improved by fixing one parameter. From this test, we observed that the estimation of H₀ is largely unaffected by the size of the pixel subsample, whereas the estimation of the energy densities depends highly on the size of the pixel subsample. For more robust results on { Ω_M,Ω_Λ } , and subsequently on Ω_κ, we would need a SN survey that covered the sky more homogeneously so that all pixels were sufficiently sampled both in number of supernovae and in redshift range.

We also computed the power spectrum of the corresponding maps of the parameters. As a result of the size of each pixel, which was selected so that all pixels contained SNe, the largest multipole accessible is the octopole ℓ = 3. Previous measurements were aimed at the measurement of the dipole ℓ = 1 only.

For an analytical toy model of an inhomogeneous ensemble of homogeneous pixels, we derived the backreaction term in the deceleration parameter due to the fluctuations of H₀ across the sky. The overall deceleration parameter, computed as the average over the pixels assuming (a) no backreaction and (b) backreaction, results in a difference of order 10^-3 between the two averaging methods. This means that, for the toy model considered, the backreaction generates insignificant extra acceleration.

This paper is a proof of concept to measure deviations from homogeneity. Despite the limitations of the SN data available to date, we demonstrated that this method can measure deviations from homogeneity in the cosmological parameters across the sky. Investigating the averaging variability can be relevant to look for a dynamical mechanism that emulates a cosmic acceleration as a derived, instead of fundamental, effect. An idea worth exploring would be to repeat the local estimation for the

dependence of other parameters without an a priori dark energy component, for example using a kinematic parametrization of the luminosity distance expressed in terms of time derivatives of the scale factor. Investigating the angular distribution of the parameters can also be relevant to characterize the pattern of overdensed regions (filaments) and underdensed regions (voids) observed on large scales. Another idea worth exploring would be to correlate the angular distribution of the cosmological parameters with probes of the large-scale structure from which a pattern of matter distribution might emerge. Yet another idea would be to use the σ₈–H₀ degeneracy to explore the implications of the angular variation of the cosmic expansion to neutrino physics (∑ m_ν or N_eff).

Further investigation along these lines can be important for the formulation and validation of inhomogeneous models on the basis of an unprecedented amount of data that we anticipate having in the near future.

---

Acknowledgments

The authors thank I. Hooke, V. Marra, N. Nunes, and S. Räsänen for useful discussions. CSC is funded by Fundação para a Ciência e a Tecnologia (FCT), Grant No. SFRH/BPD/65993/2009.

References

Appendix A: Estimation method

We aim to find the values { Ω_M,Ω_Λ,H₀ } that best fit μ_theo to μ_data. The probability that a given combination of values { Ω_M,Ω_Λ,H₀ } fits the data is known as the posterior probability P_post and is defined as (A.1)The probability of the data given the values for the model parameters (the numerator) is known as the likelihood, and the integral of the likelihood over the parameters’ space (the denominator) is known as the evidence. The prior probability P_prior describes informed constraints that we impose a priori on the parameters.

We assume that the measurements of μ at the different zs are independent and that the noise is Gaussian as a first approximation. Then the likelihood is given by

Each combination of values { Ω_M,Ω_Λ,H₀ } _trial is generated from a Gaussian distribution and defines a trial. We implement the Metropolis-Hastings algorithm to produce Markov chains where the samples are generated from the conditional posterior distribution with the following priors: (a) Θ(Ω_M)Θ(Ω_Λ)Θ(H₀) and (b) Θ(1−Ω_M)Θ(1−Ω_Λ), with Θ standing for the Heaviside function. For each trial we define (A.5)and compare the likelihood between the trial combination and the previous combination by taking the ratio of the corresponding likelihoods times the prior as follows: (A.6)If the ratio is larger than one, then the trial combination is accepted and added to the chain. Otherwise, the ratio is compared with a trial from a uniform distribution defined in [0,1] , and if the ratio is larger, then the trial combination is accepted and added to the chain; otherwise, the trial combination is rejected and the previous combination is added to the chain instead.

All Tables

All Figures