© ESO, 2014

1. Introduction

While the discovery of the (apparently) accelerated expansion of the Universe established a non-vanishing dark energy contribution in the framework of the standard cosmological model (Riess et al. 1998; Perlmutter et al. 1999), these observations also motivated many researchers to question the theoretical foundations the standard model is built upon. One of these foundations is the cosmological principle, which asserts that our Universe is spatially isotropic and homogeneous when averaged over sufficiently large scales (≳100 Mpc). This assumption is truly remarkable because if it is correct, it implies that the large-scale geometry of our Universe – which is notably inhomogeneous on small scales – is accurately described by the simple class of spatially isotropic and homogeneous Friedman-Lemaître-Robertson-Walker (FLRW) models (Robertson 1935; Walker 1935). Statistical isotropy about our position has been confirmed by the remarkable uniformity of the cosmic microwave background (CMB) spectrum (Bennett et al. 2013; Planck Collaboration XVI 2014). In contrast, statistical homogeneity on large scales (~Gpc scales) is hard to confirm, mainly because it is difficult to distinguish a temporal from a spatial evolution on the past light cone (see Maartens 2011; Clarkson 2012, for reviews).

This uncertainty inspired many to study inhomogeneous cosmologies (see Marra & Notari 2011; Bolejko et al. 2011, for comprehensive reviews), including non-Copernican models that explain the apparent accelerated expansion of the Universe by means of radial inhomogeneities without requiring any form of dark energy. The basic idea behind these alternative models is quite simple because we know from observations and numerical simulations that the large-scale structure of the Universe consists of filaments and voids (see e.g. Hogg et al. 2005; Springel et al. 2005; Labini 2010, 2011; Einasto et al. 2011b,a; Scrimgeour et al. 2012; Clowes et al. 2013; Nadathur 2013; Nadathur & Hotchkiss 2014; Sutter et al. 2014; Melchior et al. 2014). Einstein’s General theory of Relativity tells us that the expansion rate in space-time regions with lower matter density should be higher than in regions with a higher matter density. If we were to live in a large-scale under-density, the local expansion rate around us would be higher than the average expansion rate in the background. Light-rays propagating from distant sources to us – the observer – would therefore feel an accelerated expansion rate along their path. In comparison to the standard cosmological model, these scenarios hence replace a cosmic acceleration in time (due to dark energy) by a spatially varying expansion rate.

One particular, exact inhomogeneous cosmological model that has extensively been studied is the Lemaître-Tolman-Bondi (LTB) model, which is a spherically symmetric, but radially inhomogeneous dust solution of Einstein’s field equations (Lemaître 1933; Tolman 1934; Bondi 1947). The spatial hypersections of LTB models are spherically symmetric only about one point, and to conserve the remarkable uniformity of the CMB spectrum, we would have to live very close to the symmetry centre (Alnes & Amarzguioui 2006; Foreman et al. 2010). Interpreted as a faithful representation of the Universe, these void models breach the Copernican principle and require a tremendous fine-tuning of our position in the Universe. The plausibility of such scenarios is therefore more than dubious (see Célérier 2012, however, for interesting thoughts). The standard cosmological model also requires significant fine-tuning, however, which gives rise to controversial philosophical discussions (see e.g. Durrer & Maartens 2008; Amendola & Tsujikawa 2010). In this work, we simply demonstrate that a quite general class of LTB void models without dark energy is inconsistent with current observational data, which allows us to set the philosophical discussion aside.

The vast literature on inhomogeneous cosmologies – in particular LTB models – is summarised in the review articles by Bolejko et al. (2011) and Marra & Notari (2011), which allows us to only focus on the aspects that are particularly relevant for this work. For reasons to be clarified later, we only discuss LTB models with synchronous Big Bang throughout, meaning that the bang time function is constant and the Universe has the same global age everywhere. In the first part of this work, we additionally set the cosmological constant to zero. It has long been known that LTB models without dark energy can easily fit supernovae, explaining the apparent acceleration of the Universe by a Gpc-scale void around us (Célérier 2000). In addition, these models can be tuned to fit the small-angle CMB spectrum (Zibin et al. 2008). However, most recent studies agree that a good fit to the CMB requires an unrealistically low local Hubble rate of H₀ ≲ 60 km s^-1 Mpc^-1, which is strongly inconsistent with the observed value of H₀ = (73.8 ± 2.4) km s^-1 Mpc^-1 measured with Cepheid-calibrated type Ia supernovae (Biswas et al. 2010; Moss et al. 2011; Riess et al. 2011; Bull et al. 2012; Zumalacárregui et al. 2012). We here already anticipate that we shall finally arrive at the very same conclusion.

Nevertheless, we believe that our work meaningfully complements the current literature mainly for the following reasons: for the first time, we compare LTB models with the latest Planck 2013 data (Planck Collaboration XVI 2014). This is interesting because the Planck data favour a lower Hubble rate than previous CMB experiments. Moreover, we advocate the use of a recently developed technique for analysing CMB spectra in a model-independent manner, which is particularly useful for investigating alternative cosmological models. Secondly, most recent studies have assumed certain functional forms for the mass or curvature profile of LTB models that may be considered characteristic for voids. These empirical parametrisations might simply be too restrictive, however, and certainly impose artificial constraints on the models when performing maximum-likelihood estimates. To our knowledge, only Zibin et al. (2008) and Moss et al. (2011) considered more flexible spline interpolations for the mass profile of LTB voids.

We extend these ideas and introduce alternative, flexible parametrisations of the local matter density profile, aiming to impose as little a priori constraints on the detailed form as possible. A Monte Carlo technique in combination with recent observational data allows us to systematically vary the matter density profiles of LTB models and derive statistical constraints on the favoured profile shapes. Moreover, we demonstrate that even the enormous flexibility of radially fine-tuned models does not suffice to simultaneously fit the observed local Hubble rate and the CMB.

After this detailed discussion, we provide simple theoretical arguments that explain why not even heavily fine-tuned LTB models without dark energy can be reconciled with current observational data. We then discuss why it seems most natural to extend our models by a non-zero cosmological constant. The resulting ΛLTB models constitute a powerful framework for conducting general tests of the cosmological principle (Marra & Pääkkönen 2010). Recently, Valkenburg et al. (2014) proposed a new technique, based on ΛLTB models, for placing constraints on violations of the Copernican principle. Furthermore, by marginalising over all possible inhomogeneities, these authors derived first observational constraints on the cosmological constant that are free of the usual homogeneity prior (see also Marra et al. 2013b, for a similar ansatz). Marra et al. (2013a) used ΛLTB models to investigate whether fluctuations in the local matter density profile can alleviate the well-known discrepancy between the high local Hubble rate as measured by Riess et al. (2011) and the lower one derived from the Planck 2013 data. Obviously, our previously developed, flexible parametrisations of the local matter density profile in combination with a Monte Carlo technique constitute an ideal tool for conducting similar studies. Again, by imposing (almost) no a priori constraints on the detailed shape of the density profiles, our approach is more general than previous works on this field.

The structure of this paper is as follows: in Sect. 2, we explain our general ansatz for the metric of a Universe that may radially be inhomogeneous on large scales. Section 3 describes the observational data used to constrain LTB models and also discusses some ambiguities that have to be taken into account when fitting non-standard cosmological models to these data. In Sect. 4, we introduce flexible parametrisations of the local matter density profile and describe the algorithm that allows us to statistically constrain the shape of the profile functions. Section 5 summarises our main results of a long series of tests, comparing numerous LTB models without cosmological constant with different combinations of observational data. In Sect. 6, we provide some simple theoretical arguments that explain our empirical results from Sect. 5. We then extend our models by a non-zero cosmological constant and discuss general probes of the cosmological principle in Sect. 7. Finally, we present our conclusions in Sect. 8.

2. LTB ansatz for the metric

We make the following simplified ansatz for the metric of the observable Universe: we maintain the standard inflationary paradigm and assume that the early Universe was highly homogeneous at least until the time of decoupling. Given the remarkable uniformity of the CMB spectrum, we assume that spherical symmetry about our position was conserved until the present epoch. However, we breach the cosmological principle by allowing a fine-tuned, radially inhomogeneous matter density profile on large scales. As an example, we could envisage to live at the centre of a Gpc-scale void that emerged from an – admittedly extreme – under-density of the primordial matter distribution.

Mathematically, this can be realised by describing the local Universe around us with a spherically symmetric, but radially inhomogeneous LTB model that is asymptotically embedded in a homogeneous FLRW background. The LTB model is a spherically symmetric dust solution of Einstein’s field equations with zero pressure (p = 0), that is, the stress-energy tensor is T^μν = ρu^μu^ν. This should be an excellent approximation because the local Universe we observe is also a Universe at a late evolutionary epoch, and at this stage, pressure should be completely negligible for describing the large-scale matter distribution. Following the findings of previous work (Biswas et al. 2010; Zumalacárregui et al. 2012), we generally embed the LTB model into an FLRW background with non-zero curvature to improve the fit to the CMB.

LTB models have already been covered at length in the literature (see e.g. Plebanski & Krasinski 2006; Enqvist 2008; Bolejko et al. 2011, for many details). We therefore only briefly summarise the most important equations here. In coordinates comoving with dust and in the synchronous time gauge, the LTB metric can be written as (with c = 1) $\mathrm{d}{\mathit{s}}^{\mathrm{2}}\mathrm{=}\mathrm{-}\mathrm{d}{\mathit{t}}^{\mathrm{2}}\mathrm{+}\frac{{\mathit{R}}^{\mathrm{\prime }\mathrm{2}}\mathrm{\left(}\mathit{t,r}\mathrm{\right)}}{\mathrm{1}\mathrm{+}\mathrm{2}\mathit{E}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\hspace{0.17em}\mathrm{d}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{R}}^{\mathrm{2}}\mathrm{\left(}\mathit{t,r}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathrm{\Omega }}^{\mathrm{2}}\mathit{,}$(1)where R(t,r) is the areal radius function, a prime denotes the derivative with respect to the radial coordinate r and E(r) is the curvature function. The time evolution of the areal radius function follows from Einstein’s field equations and reads $\mathit{Ṙ}\mathrm{2}\mathrm{\left(}\mathit{t,r}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\mathit{M}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathit{R}\mathrm{\left(}\mathit{t,r}\mathrm{\right)}}\mathrm{+}\mathrm{2}\mathit{E}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{+}\frac{\mathrm{\Lambda }}{\mathrm{3}}\hspace{0.17em}{\mathit{R}}^{\mathrm{2}}\mathrm{\left(}\mathit{t,r}\mathrm{\right)}\mathit{,}$(2)where a dot indicates the derivative with respect to time, M(r) is the active gravitational mass inside a spherical shell at radius r, and Λ is the cosmological constant. The second Einstein field equation relates the effective gravitational mass M(r) to the local matter density ρ(t,r), $\frac{{\mathit{M}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{R}}^{\mathrm{2}}\mathrm{\left(}\mathit{t,r}\mathrm{\right)}{\mathit{R}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{t,r}\mathrm{\right)}}\mathrm{=}\mathrm{4}\mathit{\pi G}\hspace{0.17em}\mathit{\rho }\mathrm{\left(}\mathit{t,r}\mathrm{\right)}\mathit{.}$(3)Equation (3)will be key for our algorithm in Sect. 4.2, because we parametrise the local matter density at the present time with flexible basis functions and use Eq. (3)to derive the corresponding LTB metric.

After a separation of variables, Eq. (2)can be integrated in time, ${\mathit{t}}_{\mathrm{0}}\mathrm{-}{\mathit{t}}_{\mathrm{B}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\underset{\mathrm{0}}{\overset{\mathit{R}\mathrm{\left(}{\mathit{t}}_{\mathrm{0}}\mathit{,r}\mathrm{\right)}}{\mathrm{\int }}}\frac{\mathrm{1}}{\sqrt{\frac{\mathrm{2}\mathit{M}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}}\mathrm{+}\mathrm{2}\mathit{E}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{+}\frac{\mathrm{\Lambda }}{\mathrm{3}}\hspace{0.17em}\begin{array}{c}\mathrm{˜}\\ {\mathit{R}}^{\mathrm{2}}\end{array}\mathrm{\left(}\mathit{t,r}\mathrm{\right)}}}\hspace{0.17em}\mathrm{d}\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\mathit{,}$(4)where t₀ denotes the present time and t_B(r) is the so-called bang time function. In general, t_B(r) can be an arbitrary function of r, which means that the Big Bang does not need to occur synchronously, as in FLRW models. It can be shown, however, that fluctuations in the bang time function can be identified with decaying modes in linear perturbation theory (Silk 1977; Zibin 2008). Thus, going back in time, these decaying modes would correspond to inhomogeneities at early times. To conserve the remarkable homogeneity of the CMB spectrum, fluctuations in the age of the Universe must have been smaller than a few hundred years (Bolejko 2009), which is substantially smaller than the present age of the Universe (t_B(r) ≪ t). Complying with our initial assumption of a homogeneous, early Universe (even in regions from which we do not observe CMB photons), we can thus safely neglect the bang time function and assume that it is zero, t_B(r) = const. = 0. This assumption is particularly important for the CMB analysis described in Sect. 3.3. Without cosmological constant (Λ = 0), the integral in Eq. (4)can be solved parametrically and – dependent on the sign of the local curvature E(r) – leads to an elliptic (E< 0), parabolic (E = 0), or hyberbolic (E> 0) evolution. For completeness, we provide these parametric solutions with some more useful expressions in Appendix A. With a non-zero cosmological constant, no general parametric solution exists. Instead, the elliptic integral in Eq. (4)has to be computed numerically, which for certain parameter combinations quickly turns into a difficult numerical problem involving singularities, slow convergence (if at all), and poor precision. However, Valkenburg (2012) showed that the integral in Eq. (4)can be transformed to the so-called Carlson symmetric forms of elliptic integrals. This canonical set of elliptic integrals can be solved using iterative algorithms, which are fast, robust, and quickly converge with machine precision (Carlson 1995). We refer to Valkenburg (2012) for the full derivation, extensive calculations, and many useful details concerning the numerical implementation.

To compare LTB models with observational data, we need to solve for the distance-redshift relation on the past null cone (PNC) of a central observer. Inward radial null geodesics integrated backwards in time are described by the following set of first-order, ordinary differential equations: $\begin{array}{ccc}& & \frac{\mathrm{d}\mathit{t}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{d}\mathit{r}}\mathrm{=}\mathrm{-}\frac{{\mathit{R}}^{\mathrm{\prime }}\left[\mathit{t}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{,r}\right]}{\sqrt{\mathrm{1}\mathrm{+}\mathrm{2}\mathit{E}\mathrm{\left(}\mathit{r}\mathrm{\right)}}}\mathit{,}\\ & & \frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathit{z}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\frac{\mathrm{d}\mathit{z}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{d}\mathit{r}}\mathrm{=}\frac{\mathit{Ṙ}\mathrm{\prime }\left[\mathit{t}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{,r}\right]}{\sqrt{\mathrm{1}\mathrm{+}\mathrm{2}\mathit{E}\mathrm{\left(}\mathit{r}\mathrm{\right)}}}\mathrm{·}\end{array}$These equations can be numerically integrated. After solving for t(r) and z(r) on the PNC, we can numerically invert these relations to arbitrarily transform between t, r and z. For instance, given a redshift z, we can infer the corresponding radius r(z) and compute the time t [ r(z) ] at which an incoming radial null geodesic was at this position. As can be seen from Eq. (1), the angular diameter distance in LTB models is simply given by the areal radius function, ${\mathit{d}}_{\mathrm{A}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathit{R}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathit{R}\left[\mathit{t}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{,r}\mathrm{\left(}\mathit{z}\mathrm{\right)}\right]\mathit{.}$(7)The luminosity distance then follows from the reciprocity theorem (Etherington 1933, 2007; Ellis 2009): ${\mathit{d}}_{\mathrm{L}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{2}}\hspace{0.17em}{\mathit{d}}_{\mathrm{A}}\mathrm{(}\mathit{z}\mathrm{)}\mathit{.}$(8)These equations briefly summarise the most important quantities that are needed for the following sections. As can be seen, LTB models generally depend on three functions: the bang time function t_B(r), the curvature function E(r), and the mass function M(r). However, the metric and all previous formulae are invariant by diffeomorphism symmetry, including coordinate transformations of the form r = f(r′). This gauge freedom can be used to eliminate one function, implying that the physical evolution of LTB models is fully determined by only two free functions. We discussed above that the bang time function can be set to zero by demanding a homogeneous early Universe. This leaves us with only one arbitrary function. In Sect. 4.2, we describe in detail how the curvature and mass functions can be derived from a matter density profile at a given time.

3. Observational data

This section outlines what observational data we used to constrain LTB models. In this context, two complications have to be considered.

Firstly, many cosmological data sets are routinely reduced under the implicit assumption that our Universe is, on large scales, properly described by an FLRW metric. When fitting alternative cosmological models to these data, special care has to be taken that only model-independent observational constraints are used.

Secondly, linear perturbation theory in LTB models is substantially more complicated than in FLRW models. This is mainly because scalar and tensorial perturbations can couple on inhomogeneous backgrounds. Although there has been great progress in developing a gauge-invariant linear perturbation theory for LTB models, the problem has yet to be fully solved (Zibin 2008; Clarkson et al. 2009; February et al. 2014). To be conservative, we can therefore only use observables that do not depend on the details of structure formation.

Perhaps we should emphasise this limitation of the present work. The current state of the theory only allows us to explore constraints on the global structure of the smooth, unperturbed space-time geometry on large scales. We cannot yet realistically model the evolution of linear perturbations in radially inhomogeneous space-times because we lack proper theoretical tools for predicting the statistical properties of the perturbed matter density field, for example. Observations of the local matter distribution, the large-scale structure (including voids and filaments), cluster number counts or galaxy-galaxy correlation functions can therefore not yet be taken into account for constraining the underlying space-time geometry. We also have to exclude constraints from weak-lensing spectra or baryonic acoustic oscillations, for instance, which are widely used to constrain homogeneous and isotropic cosmologies. We are currently working on advancing linear perturbation theory in LTB models, and will include some of the observables listed above in future work.

The lack of linear perturbation theory is not important for Sect. 5, where we show that only constraints from the local Hubble rate and the CMB are sufficient to rule out LTB models without cosmological constant. For the general probe of the cosmological principle presented in Sect. 7, it would certainly be desirable to have more data. On the other hand, Valkenburg et al. (2014) showed that the cosmological observables used in this work are currently the most constraining.

3.1. Local Hubble rate

Riess et al. (2011) used a nearby sample ${\mathrm{(}}^{}\mathrm{0.023}\mathit{<}\mathit{z}\mathit{<}{\mathrm{0.1}}^{\mathrm{)}}$ of Cepheid-calibrated type Ia supernovae to measure the local Hubble rate with a remarkable precision: ${\mathit{H}}_{\mathrm{0}}\mathrm{=}\mathrm{(}\mathrm{73.8}\mathrm{\pm }\mathrm{2.4}\mathrm{)}\hspace{0.17em}\mathrm{km}\hspace{0.17em}{\mathrm{s}}^{-1}\hspace{0.17em}{\mathrm{Mpc}}^{-1}\mathit{.}$(9)Up to now, this is probably the most accurate determination of the local expansion rate of the Universe. Although Efstathiou (2014) recently reanalysed the data from Riess et al. (2011) and proposed a corrected, slightly lower value of H₀ = (72.5 ± 2.5) km s^-1 Mpc^-1; we used the original measurement to constrain our models because the proposed corrections are not significant for our main conclusions.

One essential property of LTB models is their radially dependent expansion rate. More precisely, we can define a longitudinal Hubble rate, H_L(t,r) = Ṙ′(t,r) /R′(t,r), which describes the expansion rate along the radial direction, and a transversal Hubble rate, H_T(t,r) = Ṙ(t,r) /R(t,r), which describes the expansion rate of the individual spherical shells. It is therefore not a priori clear how LTB models should be compared with the above-mentioned measurement.

To mimic the procedure of Riess et al. (2011), we used the following approach: independent of the cosmological model, the (local) luminosity distance can be considered an analytic function of redshift z and hence expanded in a Taylor series, $\begin{array}{ccc}{\mathit{d}}_{\mathrm{L}}\mathrm{\left(}\mathit{z}\mathrm{\right)}& \mathrm{=}& \frac{\mathit{c}}{{\mathit{H}}_{\mathrm{0}}}\left[\mathit{z}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{q}}_{\mathrm{0}}\mathrm{)}\hspace{0.17em}{\mathit{z}}^{\mathrm{2}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{6}}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{q}}_{\mathrm{0}}\mathrm{-}\mathrm{3}{\mathit{q}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{+}{\mathit{j}}_{\mathrm{0}}\mathrm{)}\hspace{0.17em}{\mathit{z}}^{\mathrm{3}}\right]\\ & & \end{array}$(10)where q₀ and j₀ are the deceleration parameter and the jerk, respectively. Like Riess et al. (2011), we fixed these two parameters to q₀ = −0.55 and j₀ = 1. To compute the effective local Hubble rate of a specific LTB model, we first tabulated the luminosity distance d_L(z_i) (cf. Eq. (8)) at N equidistantly spaced steps in the considered redshift range, 0.023 <z₁<...<z_N< 0.1, and then calculated the best-fitting (least-squares) Hubble rate H_LS through these data points using Eq. (10). The deviation from the observed value is then quantified with a simple chi-square, ${\mathit{\chi }}_{{\mathit{H}}_{\mathrm{0}}}^{\mathrm{2}}\mathrm{=}\frac{{\mathrm{(}{\mathit{H}}_{\mathrm{LS}}\mathrm{-}{\mathit{H}}_{\mathrm{0}}\mathrm{)}}^{\mathrm{2}}}{{\mathit{\sigma }}_{{\mathit{H}}_{\mathrm{0}}}^{\mathrm{2}}} \mathrm{(}{\mathit{\sigma }}_{{\mathit{H}}_{\mathrm{0}}}\mathrm{=}\mathrm{2.4}\hspace{0.17em}\mathrm{km}\hspace{0.17em}{\mathrm{s}}^{-1}\hspace{0.17em}{\mathrm{Mpc}}^{-1}\mathrm{)}\mathit{.}$(11)Using this approach, we defined an averaged, effective local Hubble rate for LTB models that closely mimics the one measured with observed type Ia supernovae.

3.2. Supernovae

To constrain the shape of the luminosity distance at even higher redshifts, we used the Union2.1 compilation released by the Supernova Cosmology Project (Suzuki et al. 2012). This catalogue contains 580 uniformly analysed type Ia supernovae and extends out to redshift z ~ 1.5. Currently, the Union2.1 compilation is the largest, and most recent, publicly available sample of standardized type Ia supernovae.

Although the shape of type Ia supernovae light-curves is empirically well understood, their absolute magnitude is essentially unknown and needs to be calibrated. Samples like the Union2.1 compilation are therefore reduced by fixing the Hubble rate to an arbitrary value. It is important to remove this artificial constraint from the data when fitting cosmological models. This can either be achieved by analytically marginalising over the assumed Hubble rate – or, equivalently, the absolute magnitude – (see Bridle et al. 2002; or Appendix C.2 of Biswas et al. 2010), or by using the elegant weight matrix formalism described in the appendix of Amanullah et al. (2010). The weight matrix is constructed from the full covariance matrix with systematics (e.g. host mass correction) and incorporates the marginalisation over various nuisance parameters. In particular, the marginalisation over the Hubble rate is included. We use the weight matrix formalism to perform likelihood estimates in the following sections.

Finally, we note that several authors argued that supernova samples reduced with the SALT-II light-curve fitter from Guy et al. (2007) are systematically biased towards the standard cosmological model and tend to disfavour alternative cosmologies (Hicken et al. 2009; Kessler et al. 2009; Smale & Wiltshire 2011). The Union2.1 compilation was reduced with the SALT-II fitter, and hence this potential penalty would also affect the goodness-of-fit of LTB models. However, we can safely neglect this problem for mainly two reasons. Firstly, we did not conduct a detailed, statistical comparison (e.g. Bayesian model comparison) between the standard cosmological model and LTB models here, and therefore the small, systematic corrections are irrelevant for our main conclusions. Secondly, and more importantly, we demonstrate that the supernova data do not impose tight constraints on LTB models anyway. Instead, supernovae can easily be fitted with a variety of different density profiles and certainly do not cause the tension between observations and LTB models that we focus on.

3.3. Cosmic microwave background

The standard approach for analysing CMB spectra is inherently based on the assumption of a spatially isotropic and homogeneous Universe (Bennett et al. 2013; Planck Collaboration XVI 2014). Primary anisotropies are calculated by numerically solving the Boltzmann equations on an FLRW background (Lewis et al. 2000; Lesgourgues 2011). Secondary anisotropies are caused by different forms of interactions between cosmic structures and the CMB photons after the time of decoupling, such as the late integrated Sachs-Wolfe (ISW) effect or rescattering during reionisation. As such, all these processes depend on the details of structure formation and are commonly modelled using perturbation theory in FLRW models. It is obvious that our ansatz for the metric in Sect. 2 strongly violates these usual assumptions. How then can we self-consistently use CMB spectra to constrain LTB models?

We recall that we wish to retain the inflationary paradigm and assume that the early Universe is highly homogeneous. This means that, even in our approach, the early Universe can effectively be described by a FLRW metric. Thus, the modelling of primary CMB anisotropies does not differ from the standard approach at least until the time of decoupling. This was also the motivation for excluding variations of the bang time function t_B(r) from Eq. (4). If we were to drop the assumption of a homogeneous early Universe, we would have to develop a general relativistic formalism for calculating CMB anisotropies on inhomogeneous backgrounds. Clearly, this would go far beyond the scope of this work.

Properly treating secondary CMB anisotropies in LTB models is, however, more complicated because linear perturbation theory is still being developed and cannot yet be used to calculate these secondary effects. We can circumvent this problem by following the elegant work of Vonlanthen et al. (2010), which describes a method for analysing the CMB in a manner that is as independent as possible of late-time cosmology. To this end, the authors begin with identifying the three dominant imprints that the late cosmological model leaves on the observed CMB spectrum. Firstly, CMB photons are lensed as they traverse non-linear structures. The late ISW effect dominates the CMB spectrum on large scales (l ≲ 40). Consequently, low multipoles strongly depend on the detailed properties of the late cosmological model. Secondly, the overall amplitude of the CMB spectrum at l ≫ 40 is reduced by scattering processes during the epoch of reionisation. This suppression is usually parametrised by the factor e^−2τ, where τ is the reionisation optical depth. Thirdly, the angular diameter distance to the last scattering surface (LSS) determines the angular scales of the acoustic peaks. This is a simple projection effect. Variations of the angular diameter distance shift the CMB spectrum in multipole space.

If large scales (l ≲ 40) are excluded from the CMB analysis, the dominant perturbations of the primordial CMB spectrum due to the late cosmological model can therefore be parametrised by a rescaling factor α of the global amplitude and a shift parameter β. Rephrased more mathematically, the expansion coefficients of the CMB power-spectrum in multipole space transform as C_l → αC_βl.

This simple approximation is already astonishingly accurate for large parts of the CMB spectrum. But the crucial point is that this insight allowed Vonlanthen et al. (2010) to encode unknown secondary effects in carefully chosen nuisance parameters (e.g. the global amplitude of the CMB spectrum). These additional parameters can be built into standard parameter estimation codes used for analysing CMB spectra. To derive minimal constraints that do not depend on the detailed properties of the late cosmological model, one simply has to marginalise over these newly introduced nuisance parameters.

Except for excluding high multipoles (l> 800) from the analysis, Vonlanthen et al. (2010) neglected the impact of gravitational lensing on the CMB spectrum. This was well justified given the accuracy of CMB experiments at the time of publishing. Motivated by the improved accuracy of modern CMB experiments, however, Audren et al. (2013a) extended the original method of Vonlanthen et al. (2010) by introducing a new technique for additionally marginalising over the CMB lensing contamination at all multipole orders. Initially, we used this advanced technique in combination with the publicly available parameter estimation code Monte Python (Audren et al. 2013b) to derive model-independent constraints from the latest Planck data (Planck Collaboration XVI 2014). While we prepared this manuscript, Audren (2014) published another slight refinement of their own method and also applied it to the Planck 2013 data. We therefore used the results published in Audren (2014) to estimate likelihoods in the following sections. The relevant model-independent constraints including one-sigma errors are summarised in Table 1.

Finally, we need to explain how LTB models can actually be fitted to these constraints. We begin by noting that the redshift of decoupling z_∗ generally is a function of the physical matter densities w_b and w_c. Assuming standard radiation content, however, this dependence is only weak (Hu & Sugiyama 1996). For simplicity, we can thus fix the decoupling redshift to z_∗ = 1090. Following previous works (see e.g. Zibin et al. 2008; Biswas et al. 2010; Moss et al. 2011, for detailed descriptions), we used an effective FLRW observer approach for computing the angular diameter distance d_A(z_∗) to the LSS: radial null geodesics are numerically integrated in LTB models only up to an embedding redshift z_b. This redshift has to be chosen such that the LTB models are already sufficiently homogeneous (i.e. the gravitational shear ${\mathit{\sigma }}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}{}^{\mathrm{(}}\mathit{H}_{\mathrm{L}}\mathrm{-}{\mathit{H}}_{\mathrm{T}}{}^{\mathrm{)}}$ vanishes) and radiation is still negligible. Typically, z_b ≈ 150 fulfils these criteria.

After reaching this embedding redshift, we continued to compute the distance-redshift relation in an effective FLRW background model up to the decoupling redshift z_∗ = 1090. The effective FLRW model was chosen such that a fictitious observer in this background model would observe the same CMB spectrum. This approach is beneficial because (1) solving the distance-redshift relation in FLRW models is computationally substantially cheaper, and (2) radiation can be included. Lastly, the appropriately scaled matter densities and the calculated angular diameter distance can be compared with the model-independent constraints from Table 1 using a simple chi-square.

3.4. Kinetic Sunyaev-Zel’dovich

Cosmic microwave background photons traversing galaxy clusters interact with the hot intra-cluster gas through inverse Compton scattering. These interactions cause the Sunyaev-Zel’dovich effect, which manifests itself by characteristic distortions of the energy spectrum of the rescattered photons. The dominant contribution in galaxy clusters is due to the thermal Sunyaev-Zel’dovich effect (Sunyaev & Zeldovich 1970, 1972): thermal energy is transferred from the hot intra-cluster gas to the CMB photons, causing a redistribution of photons from lower to higher energy states. A second-order contribution is due to the kinetic Sunyaev-Zel’dovich (kSZ) effect (Sunyaev & Zeldovich 1980): galaxy clusters moving with a non-zero peculiar velocity with respect to the rest frame of the CMB observe an anisotropic CMB spectrum with non-zero dipole moment. This dipole induces a characteristic shift in the spectrum of the reflected CMB photons, which is similar to the relativistic Doppler effect. In principal, the kSZ effect of galaxy clusters can therefore be used to measure their peculiar velocities relative to the CMB.

As was first shown by García-Bellido & Haugbølle (2008b), the kSZ effect can be used as a powerful probe of radially inhomogeneous LTB models. To understand this, we recall that LTB models are isotropic only about the symmetry centre at the coordinate origin. Off-centre observers therefore generally see an anisotropic CMB spectrum (see Alnes & Amarzguioui 2006, for detailed calculations). To first order, the anisotropy can be approximated as a pure dipole, which because of the symmetry of the problem is aligned along the radial direction. This effect is depicted in Fig. 1, which shows an off-centre galaxy cluster that rescatters CMB photons that arrive from the LSS.

The extreme redshifts (as seen from the galaxy cluster) are observed along the radial direction, for CMB photons with redshift z_in that propagated through the centre of the inhomogeneity and, in the opposite direction, for CMB photons with redshift z_out arriving from outside. In the case of a giant void scenario, for instance, photons arriving from inside the void travelled the longest distance through an underdense region. Consequently, they also spent the longest time in a space-time region with a higher expansion rate, so that their redshift z_in reaches the highest possible value. Vice versa, CMB photons arriving from outside the void exhibit the lowest redshift z_out.

Note that the kSZ effect described here is only an apparent kSZ effect; it is not caused by real peculiar motions of the galaxy clusters. Instead, the effect only appears because the space-time geometry around off-centre observers is anisotropic. Note also that the kSZ effect has a distinguished feature in comparison with most other cosmological probes. The galaxy clusters act as mirrors for CMB photons, reflecting radiation from all spatial directions. By analysing the spectrum of the reflected light, we can therefore extract information about space-time regions that would otherwise be inaccessible. In this sense, by measuring the difference between the redshift z_in and z_out, the kSZ effect allows us to (indirectly) look inside our past-null cone (cf. Fig. 1).

As already mentioned, to first order, the anisotropy observed by off-centre galaxy clusters can be approximated as a pure dipole. This approximation is sufficiently accurate as long as the effective size of the LTB inhomogeneity on the sky, as observed by the galaxy cluster, is larger than ~2π (Alnes & Amarzguioui 2006). The amplitude of the dipole is then given by $\mathit{\beta }\mathrm{=}\frac{\mathit{v}}{\mathit{c}}\mathrm{=}\frac{\mathrm{\Delta }\mathit{T}}{\mathit{T}}\mathrm{=}\frac{{\mathit{z}}_{\mathrm{in}}\mathrm{-}{\mathit{z}}_{\mathrm{out}}}{\mathrm{2}\mathrm{+}{\mathit{z}}_{\mathrm{in}}\mathrm{+}{\mathit{z}}_{\mathrm{out}}}\mathrm{·}$(12)Measurements of the kSZ effect of individual galaxy clusters are extremely difficult and suffer from low signal-to-noise ratios. Current data exhibit very large errors and therefore still need to be considered premature. However, even though the uncertainties are huge, García-Bellido & Haugbølle (2008b) and Bull et al. (2012) showed that the currently available measurements already place tight constraints on the depth and radial size of Gpc LTB voids. In this work, we use the kSZ data of the nine galaxy clusters compiled by García-Bellido & Haugbølle (2008b), assuming a conservative scatter of σ_pv = 1600 km s^-1 and zero systematic shift, v_sys = 0; see García-Bellido & Haugbølle (2008b) and references therein for a detailed discussion of the data, sources of errors, and the modelling of the likelihood.

To compute the expected kSZ effect for a given galaxy cluster, we first determined the cluster coordinates (t_cl,r_cl,z_cl) on the past-null cone of the central observer. Starting from this position, we then solved for ingoing and outgoing radial null geodesics by numerically integrating Eqs. (5)and (6)up to the LSS. This procedure yields the two redshifts z_in and z_out, which quantify the CMB dipole as seen by the off-centre galaxy cluster.

4. Monte Carlo approach for constraining the local density profile

Most studies published so far have assumed certain functional forms for either the mass profile or the curvature profile of LTB models that may be considered characteristic for voids (see e.g. Garcia-Bellido & Haugbølle 2008a; Bolejko & Wyithe 2009; Biswas et al. 2010, for typical profiles). Of course, these empirical functional forms impose artificial constraints on the models when performing maximum likelihood estimates.

Our approach here is different because we parametrise the matter density profile of LTB models as flexibly as possible, imposing few a priori constraints on its detailed shape. To derive statistical constraints, we then use a parameter estimation code to systematically vary the density profile. We wish to determine the shape of the favoured profile, and also how tight the constraints on the detailed shape are. We also investigate whether highly flexible profiles allow us to mitigate the reported tension between measurements of the local Hubble rate and the CMB data (Biswas et al. 2010; Moss et al. 2011; Bull et al. 2012; Zumalacárregui et al. 2012). To this end, we proceed by first discussing our choices for the parametrisation of the local density profile, continue with the algorithm that allows us to derive LTB models from these profiles, and conclude this section by explaining some details of the parameter estimation technique.

4.1. Flexible models for the local density profile

Broadly speaking, we have two options to parametrise the local density profile in the most flexible way. We can represent the density profile by a general interpolation scheme, or alternatively decompose the void profile function into a series of appropriately chosen, orthonormal basis functions of the radial coordinate. After some initial testing, we decided to use the following parametrisations:

• 1.

  Linear interpolation: the first free model parameter is the radial size L of the inhomogeneity. For radii r<L, we use linear interpolation, specified by pairs of radial coordinates and corresponding matter densities: [(r₀ = 0,ρ₀), (r₁,ρ₁), ..., (r_n,ρ_n)], where r_n<L and all densities ρ_i are strictly positive. For radii r ≥ L, the density profile is constant and equals the matter density ρ_bg of the homogeneous FLRW background. M′(r) (cf. Eq. (3)), and thereby also radial derivatives of the areal radius function R(t,r) (cf. Eq. (A.10)) are not smooth if the density profile is linearly interpolated. All functions required for this work remain continuous, however, which is sufficient for the computations.

• 2.

  Cubic spline interpolation: this model has the same free parameters as the linear interpolation model: a radial size L, pairs of radial coordinates, and corresponding matter densities [r_i,ρ_i = ρ(r_i)]. Cubic splines are used to interpolate between these nodes. We chose the boundary conditions ρ(r ≥ L) = ρ_bg and ρ′(r = 0) = ρ′(L) = 0 to enforce a smooth void profile at the origin and guarantee a smooth embedding into the homogeneous FLRW background. This is essentially equivalent to the spline model used by Zibin et al. (2008) and Moss et al. (2011).

• 3.

  Laguerre polynomials: again, the first free model parameter is a radial length scale L. We introduce the scaled radial coordinate x = r/L and define the matter density profile as $\mathit{\rho }\mathrm{\left(}{\mathit{t}}_{\mathrm{0}}\mathit{,r}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{bg}}\left[\mathrm{1}\mathrm{+}\exp \mathrm{(}\mathrm{-}\mathit{x}\mathrm{)}\sum _{\mathit{i}}{\mathit{a}}_{\mathit{i}}{\mathrm{L}}_{\mathit{i}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\right]\mathit{,}$(13)where L_i denote the Laguerre polynomials, which are an orthogonal basis with respect to the inner product $\mathrm{⟨}{\mathrm{L}}_{\mathit{i}}\mathit{,}{\mathrm{L}}_{\mathit{j}}\mathrm{⟩}\mathrm{=}{{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathrm{\infty }}{\mathrm{L}}_{\mathit{i}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\hspace{0.17em}{\mathrm{L}}_{\mathit{j}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\hspace{0.17em}\exp \mathrm{(}\mathrm{-}\mathit{x}\mathrm{)}\hspace{0.17em}\mathrm{d}\mathit{x}$. In practice, we used Laguerre polynomials up to fifth order, because more degrees of freedom were not constrained by the data. In contrast to the previous interpolation schemes, the Laguerre models do not have a fixed size. Instead, their density profiles only asymptote (characteristic length scale L) to the homogeneous background density.

The different approaches are visualised in Fig. 2. We conducted detailed tests as explained in Sect. 5 with each of the three parametrisations, finding qualitatively the same results. Because the meaning of the free model parameters of the interpolation schemes is quite instructive while on the other hand the geometrical interpretation of expansion coefficients of Laguerre polynomials is not straightforward, we focus our discussion in the subsequent sections on the linear and cubic spline interpolations.

4.2. Algorithm

In our approach, LTB models are generally determined by the dimensionless Hubble parameter h = H₀/ (100 km s^-1 Mpc^-1) and the physical densities of baryonic and cold dark matter, ω_b = Ω_bh² and ω_c = Ω_ch², of the homogeneous FLRW background, as well as a list of model parameters (a₁,...,a_n) that parametrise the local matter density profile ρ(t₀,r) at the present time t₀. If not stated otherwise, these are the base parameters that are later sampled by the Monte Carlo method that was used for fitting LTB models to observational data. Given these base parameters, our algorithm for computing observables in the LTB metric can be outlined as follows:

• 1.

  Age of the Universe: as emphasised before, we demand a homogeneous early Universe and thus set the bang time function t_B(r) to zero. Therefore, the LTB patch and the homogeneous background have the same global age t₀, which we compute using the standard FLRW relation, ${\mathit{t}}_{\mathrm{0}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{H}}_{\mathrm{0}}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}\frac{\sqrt{\mathit{a}}}{\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}\mathrm{+}{\mathrm{\Omega }}_{\mathit{k}}\mathit{a}}}\hspace{0.17em}\mathrm{d}\mathit{a,}$(14)where the curvature parameter is given by Ω_k = 1 − Ω_m, that is, we explicitly consider curved FLRW backgrounds.

• 2.

  Gauge freedom: all formulae given in Sect. 2 are invariant under coordinate transformations of the form $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{r}\end{array}\mathrm{=}\mathit{g}\mathrm{\left(}\mathit{r}\mathrm{\right)}$. We exploit this gauge freedom to scale the radial coordinate r such that it equals the areal radius function at the present time: R(t₀,r) = r.

• 3.

  Mass profile: given the parameters (a₁,...,a_n), we construct the matter density profile at the present time ρ(t₀,r). We reject combinations of parameters that induce negative matter densities. This is important for the cubic spline interpolation and the Laguerre polynomials. In our gauge, the effective gravitational mass M(r) then directly follows from integrating Eq. (3), $\mathit{M}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{4}\mathit{\pi G}\underset{\mathrm{0}}{\overset{\mathit{r}}{\mathrm{\int }}}\mathit{\rho }\mathrm{\left(}{\mathit{t}}_{\mathrm{0}}\mathit{,}\begin{array}{c}\mathrm{˜}\\ \mathit{r}\end{array}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\begin{array}{c}\mathrm{˜}\\ \mathit{r}\end{array}\mathit{.}$(15)

• 4.

  Curvature profile: the curvature function E(r) is implicitly defined by Eq. (4)and can only be computed numerically. We use the parametric solutions from Appendix A in combination with the TOMS 748 root-bracketing algorithm from Alefeld et al. (1995) to determine E(r) as a function of t₀, r and M(r), see Krasiński & Hellaby (2002) for details on choosing the initial bracket. We reject mass models that require E(r) ≤ − 1 / 2 because this would cause singular line-elements (cf. Eq. (1)). The radial derivative E′(r) is numerically computed using a standard fourth-order finite differencing scheme.

• 5.

  Distance-redshift relation: we use a fifth order Dormand-Prince method – which is essentially a Runge-Kutta scheme with error control and adaptive step size (Ahnert & Mulansky 2011) – to numerically integrate the ordinary differential equations describing radial null geodesics on the PNC of a central observer (cf. Eqs. (5)and (6)). The resulting relations t(r) and z(r) are then interpolated and numerically inverted using smooth Akima splines (Akima 1970). This allows us to arbitrarily transform between t, r and z and, in particular, to compute angular diameter and luminosity distances (see Sect. 2). We discard LTB models that exhibit shell-crossings or multivalued redshifts on the PNC (Hellaby & Lake 1985).

This algorithm allowed us to compute all observable quantities that are required for performing likelihood estimates of LTB models.

4.3. Efficient statistical sampling

We are confronted with a typical parameter estimation problem. Given observational data, we need to explore the posterior distribution in a high-dimensional parameter space to estimate the most likely values of the free model parameters (h,ω_b,ω_c,a₁,...,a_n). Monte Carlo methods are the standard approach for solving this kind of problem, and most commonly, variants of the simple Metropolis-Hastings algorithm are used (Hastings 1970).

The Metropolis-Hastings algorithm has an important drawback, however: it requires a fine-tuned proposal distribution to efficiently sample the posterior. If the proposal distribution is thought of as a multivariate Gaussian, this means that each entry of the covariance matrix needs to be tuned. In our case, the matter densities at different radial coordinates can be correlated and hence the covariance matrix is non-diagonal. Consequently, we would have to hand-tune N(N + 1) / 2 unknown parameters (where N is the dimension of the parameter space). This is an extremely time-consuming task, in particular because the fine-tuned parameters strongly depend on the precise parametrisation (e.g. number and position of interpolation nodes) of the matter density profile.

After some testing, we decided to use an alternative Monte Carlo sampler: the so-called stretch-move technique, which was first introduced by Goodman & Weare (2010). This technique is affine-invariant, meaning that it performs equally well under all linear transformations of the parameter space. In particular, it is insensitive to covariances between parameters and therefore requires no fine-tuning. Goodman & Weare (2010) demonstrated the excellent performance of their algorithm (as measured by the auto-correlation time) for several pathological posterior distributions.

In addition, the stretch-move sampler simultaneously explores the parameter space with a whole ensemble of Monte Carlo walkers. The time evolution of this ensemble can easily be parallelised, which greatly reduces the required computing time (wall-clock time) on multi-core machines or large computing clusters (Akeret et al. 2013). Foreman-Mackey et al. (2013) provided an excellent discussion of the stretch-move technique and described a parallelised implementation in detail. Following this, we implemented the ensemble sampler in C++ and parallelised it with the message-passing interface (MPI) system. This enabled us to perform likelihood estimates of many different LTB models within a significantly shorter computing time.

All parameter estimations in the remainder of this work were performed with the stretch-move technique, typically using hundreds of walkers. The length of the burn-in period was estimated by means of the exponential auto-correlation time τ_exp. Convergence of the samples was ensured by letting the random walks proceed for multiple integrated auto-correlation times τ_int after removing all samples drawn during the initial burn-in period. This procedure was advocated by Akeret et al. (2013) and Allison & Dunkley (2014), who provided a detailed discussion of convergence diagnostics with the auto-correlation times described above.

5. Comparing LTB models without cosmological constant with observational data

We now compare LTB models with consecutively different combinations of observational data. This stepwise analysis allows us to carefully explain why LTB models with a constant bang time function and zero cosmological constant are inconsistent with current data. The following results are representative. For each scenario discussed, we fitted numerous different LTB models, varying the radial size of voids, changing the numbers and positions of interpolation nodes, or considering different orders of the expansion in terms of the Laguerre polynomials. Our findings with the different approaches agreed qualitatively well, therefore we only discuss simple parametrisations that already show all important characteristics.

5.1. Constraints: H₀ + supernovae

To begin with, we only analysed the constraints imposed by the local Hubble rate and supernovae. We therefore considered LTB models of fixed radial size L = 3 Gpc and parametrised the matter density profiles with linear and cubic spline interpolation schemes with three equidistant nodes at radii r₁ = 0, r₂ = 1 Gpc, r₃ = 2 Gpc. It is convenient to express the matter densities at these nodes with respect to the matter density ρ_bg of the FLRW background model, viz a_i = ρ(r_i) /ρ_bg. The priors for the parameters a_i are only bounded from below (a_i ≥ 0), so that the stretch-move walkers can essentially freely explore the physically relevant parameter space.

Anticipating the following results, we asymptotically embedded the LTB models in Einstein-de-Sitter (EdS) backgrounds for this particular test, meaning that we explicitly set the spatial curvature of the homogeneous background to zero (Ω_k = 0). We did this to demonstrate that good fits to the data can be achieved even with asymptotically flat backgrounds. As can easily be verified, the models considered are fully determined by the physical matter density ω_m of the EdS background and the profile parameters (a₁,a₂,a₃).

Figure 3 illustrates the preferred shape of the local matter density profile given the local Hubble rate and supernova data. More precisely, we show the 68% confidence intervals for the best-fitting linear and cubic spline interpolation models. Note that the constraints for the two different approaches are almost identical. As expected, the data clearly favour large and deep voids whose density profiles gradually decrease towards the origin. This shape is required to account for the apparent acceleration suggested by the supernova data. Typically, the density contrast at the origin is δρ(t₀,r = 0) /ρ_bg ≈ − 0.75. It can also be seen that the constraints on the matter density profile weaken at larger radii. This is naturally caused by the quality of the data. Small radii correspond to low redshifts. In this range, the profiles are simultaneously constrained by the measurement of the local Hubble rate (z< 0.1) and also by many supernovae with comparably small error bars. At higher redshifts, the amount of observed supernovae decreases, and at the same time, their errors increase (mainly because the supernovae become fainter). Therefore, the matter density profile can substantially fluctuate at larger radii without being penalised too strongly. This freedom also demonstrates that the LTB models considered can easily fit the data without a cosmological constant and zero background curvature. The best-fitting cubic spline model has an excellent log-likelihood value of log (ℒ) = −272.27 and fits the data just as well as the best-fitting ΛCDM model with log (ℒ) = −272.56.

Figure 4 reveals more details of the Markov Chain Monte Carlo (MCMC) simulation specifically for the cubic spline model. While the marginalised posterior distributions, the means, and the standard deviations of the density profile parameters a_i confirm the above discussion, the almost perfect linear correlations between the parameters enforce the idea of a gradually decreasing density profile towards the symmetry centre: If a_i is increased, also a_{i + 1} has to grow to generate the gradient of the expansion rate required by the data. The physical matter density ω_m determines the expansion rate of the background EdS model through $\mathit{h}\mathrm{=}\sqrt{{\mathit{\omega }}_{\mathrm{m}}}$. This is why ω_m and the profile parameters a_i are correlated: If ω_m decreases, so does the expansion rate of the background model, which means that the voids need to become even deeper to maintain the required high local Hubble rate at the origin.

The average increase of the expansion rate around the central observer due to the large-scale under-density is the last aspect worth being considered in more detail. We therefore defined δH = H_LS − H_bg, which measures the difference between the effective local Hubble rate H_LS (cf. Eq. 9) at the void centre and the global Hubble rate H_bg of the background model. From the mean value of the physical matter density ω_m in Fig. 4, we extracted that the data favour H_bg ≈ 56 km s^-1 Mpc^-1, which requires an average shift of δH ≈ 18 km s^-1 Mpc^-1 to fulfil the observationally measured Hubble rate of H₀ ≈ 74 km s^-1 Mpc^-1. The magnitude of this shift agrees well with previous works (Enqvist & Mattsson 2007). The marginalised posterior distribution of δH is plotted in Fig. 5. Even though the voids can become extremely deep, the local Hubble rate rarely increases by more than ~ 24 km s^-1 Mpc^-1. This tendency is also important for Sect. 5.3.

5.2. Constraints: Cosmic microwave background

We proceed by analysing the implications of only the model-independent CMB constraints from Table 1. As reference, we first ran an MCMC simulation to fit curved FLRW models with vanishing cosmological constant to the data. These models can fit the CMB constraints just as well as flat FLRW models with non-vanishing cosmological constant. The data favour closed FLRW models with dimensionless Hubble parameter h = 0.33 ± 0.01, matter density parameter Ω_m = 1.26 ± 0.04, and consequently a curvature parameter of Ω_k ≈ − 0.26. It is important to note that although curved FLRW models can easily fit the minimal CMB constraints without a cosmological constant, the required Hubble rate is extremely low and strongly contradicts the expansion rate measured by Riess et al. (2011).

Next, we compared the same LTB models as in the previous section (L = 3 Gpc; three equidistant interpolation nodes) with the CMB data. This time, however, we asymptotically embedded the LTB models into curved FLRW backgrounds. Of course, these cosmological models fit the data perfectly as well. The constraints on the physical matter density ω_m = ω_b + ω_c and the dimensionless Hubble parameter h of the background model are essentially identical to those found without an LTB inhomogeneity around the observer. Again, an extremely low expansion rate of h = 0.33 ± 0.01 is favoured. The shape of the local matter density profile is not constrained at all. The posteriors of the densities ρ(r_i) at the interpolation nodes are constant over the whole prior ranges, indicating that the Monte Carlo walkers can vary the shape of the local matter density profile without notably deteriorating the fit to the CMB. Moreover, the best-fitting LTB models show almost arbitrarily shaped local density profiles. These results are the same for LTB models with different radial sizes, more or less interpolation nodes, or density profiles represented by Laguerre polynomials.

We conclude that the CMB data alone only constrain the global properties of the background FLRW model. In particular, good fits without a cosmological constant require an extremely low Hubble parameter of h< 0.4. The exact functional form of the local curvature and the matter density profiles are essentially irrelevant. These results agree well with the work of Clifton et al. (2009), who found that the small angle fluctuations of the CMB spectrum only constrain spatial curvature near the surface of last scattering.

5.3. Constraints: H₀ + supernovae + cosmic microwave background

Finally, we simultaneously fitted LTB models to the local Hubble rate, supernovae, and the CMB data. Before that, however, we fitted spatially flat FLRW models with non-vanishing cosmological constant to the data. The best-fitting model has the cosmological parameters (h,Ω_m,Ω_Λ) = (0.719,0.265,0.735) and a log-likelihood value of log (ℒ) = −273.45. This model – and in particular its log-likelihood value – serves as a reference point for the following discussion.

We now considered LTB models with fixed radial size L = 3 Gpc and linearly interpolated matter density profile with three equidistant interpolation nodes at the radial coordinates r₁ = 0, r₂ = 1 Gpc and r₃ = 2 Gpc. For these models, the data favour an FLRW background with dimensionless Hubble parameter h = 0.39 ± 0.01 and matter density parameter Ω_m = 1.09 ± 0.03, demonstrating that curved FLRW backgrounds can indeed improve the fit to the CMB. The statistical constraints on the local matter density profile at the present time t₀ are shown in Fig. 6. While the overall shape of the favoured profiles agrees well with our previous results from Sect. 5.1, the voids become noticeably deeper when including the constraints from the CMB. This change becomes most apparent for the matter density at the origin, ${\mathit{a}}_{\mathrm{1}}\mathrm{=}\mathit{\rho }{}^{\mathrm{(}}\mathit{t}_{\mathrm{0}}\mathit{,r}\mathrm{=}{\mathrm{0}}^{\mathrm{)}}\mathit{/}{\mathit{\rho }}_{\mathrm{bg}\mathrm{\left(}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}}$, indicating that the best agreement with the data can be achieved with an almost vacuum solution around the observer at the coordinate centre.

This tendency demonstrates the problem of LTB models with a constant bang time function we described above: A reasonable fit to the CMB – most importantly the angular diameter distance to the LSS – requires an unrealistically low background Hubble parameter of H₀ ≈ 39 km s^-1 Mpc^-1. On the other hand, the LTB models also need to comply with the observationally measured, high local Hubble rate of H₀ ≈ 73.8 km s^-1 Mpc^-1. The Monte Carlo walkers therefore minimise the matter density at the origin to maximise the local expansion rate. However, Fig. 7 shows that even though the voids can become extremely deep, the effective local Hubble rate never exceeds ~ 62 km s^-1 Mpc^-1. The best-fitting model, for instance, has the profile parameters (a₁,a₂,a₃) = (0.04,0.22,0.32), which means that the matter density almost vanishes at the origin. The local Hubble rate still is much too low, H_LS = 56.83 km s^-1 Mpc^-1, however, and the log-likelihood value of log (ℒ) = −318.55 indicates that the model clearly cannot compete with simple, spatially flat FLRW models with cosmological constant.

As discussed in Sect. 3.1, we computed the effective local Hubble rate H_LS by means of a least-squares fit to the Taylor-expanded luminosity distance for small redshifts, z< 0.1. Legitimately, one might therefore ask whether a radially fine-tuned matter density profile in the nearby range (0 <z(r) < 0.1) might be capable of solving the above tension. We refute this idea by adding one additional interpolation node at r = 500 Mpc and re-running the fitting procedure. The best-fitting model improves only marginally, with a slightly higher effective local Hubble rate of H_LS = 57.26 km s^-1 Mpc^-1 and a still inacceptable log-likelihood value of log (ℒ) = −317.30. Consequently, we note that even radial fine-tuning of the matter density profile cannot mitigate the tension between the required high local Hubble rate and a good fit to the CMB.

Lastly, we discuss the impact of the radial void size on the fit to the data. We considered LTB models whose density profiles are represented by cubic spline interpolations with three equidistant nodes at the radial coordinates r_i = (i − 1) × L/ 3. Table 2 summarises the best-fitting models for ascending radial sizes, L = 3 Gpc, L = 4.5 Gpc and L = 7.5 Gpc. Obviously, the fit to the data significantly improves as the voids become larger. The favoured scenarios are extremely deep Gpc-scale under-densities, which are asymptotically embedded into increasingly dense FLRW backgrounds.

This antipodal behaviour is caused by the constraints from the CMB; the background FLRW models need to increase the convergence to decrease the angular diameter distance to the LSS. Again, even though the best-fitting models are essentially empty at the origin, the effective local Hubble rate is clearly too low compared with the observed value. We thus note that even gigantic voids of radial size 7.5 Gpc with radially fine-tuned matter density profile are inconsistent with current data. Clearly, these models cannot compete with the standard cosmological model, which has a log-likelihood value of log (ℒ) = −273.45.

We could continue this procedure and construct ever larger voids with more free interpolation nodes. Indeed, while testing our code, we observed that LTB voids with radial sizes ≳12Gpc slowly converge to the same log-likelihood value as the standard cosmological model. It does not make much sense, however, to seriously consider even larger voids than those presented in Table 2 as faithful representations of the observable Universe for mainly two reasons. Firstly, these models become yet more implausible with growing size. We recall that we envisaged voids that emerged from fluctuations of the primordial matter distribution. For these voids to reach a size of several Gpc today, the primordial under-densities must have been extreme. These scenarios are therefore highly unlikely within the standard inflationary paradigm (see also Marra et al. 2013a). Secondly, as was shown for example by García-Bellido & Haugbølle (2008b), Zibin & Moss (2011) or Bull et al. (2012), large voids generate a pronounced kSZ effect that is strongly inconsistent with current observational data.

The results of this section agree qualitatively well with previous works, which consistently found that LTB models that fit the CMB data exhibit unrealistically low local Hubble rates (e.g. Biswas et al. 2010; Moss et al. 2011; Bull et al. 2012). We now seek a theoretical explanation for this empirical result.

6. Theoretical arguments for considering ΛLTB models

6.1. Why vacuum solutions maximise the local Hubble rate

In the previous section, we found that huge Gpc-scale voids with almost constant vacuum solution around the origin are favoured when the local Hubble rate and the CMB data are to be fitted simultaneously. Intuitively, it is comprehensible that such vacuum solutions indeed maximise the local Hubble rate. However, we can also understand this result using Raychaudhuri’s equation (Raychaudhuri 1955), which in LTB models is $\mathit{Ḣ}\mathrm{L}\mathrm{+}\mathrm{2}\mathit{Ḣ}\mathrm{T}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{3}}{\mathit{\theta }}^{\mathrm{2}}\mathrm{-}{\mathit{\sigma }}^{\mathrm{2}}\mathrm{-}\mathrm{4}\mathit{\pi G}\hspace{0.17em}\mathit{\rho ,}$(16)with expansion ${\mathit{\theta }}^{\mathrm{2}}\mathrm{=}{{}^{\mathrm{(}}\mathit{H}_{\mathrm{L}}\mathrm{+}\mathrm{2}{\mathit{H}}_{\mathrm{T}}{}^{\mathrm{)}}}^{\mathrm{2}}$, shear ${\mathit{\sigma }}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}{{}^{\mathrm{(}}\mathit{H}_{\mathrm{L}}\mathrm{-}{\mathit{H}}_{\mathrm{T}}{}^{\mathrm{)}}}^{\mathrm{2}}$, and a dot indicating the derivative with respect to an appropriately chosen affine parameter. Generally, the shear in LTB models vanishes as r → 0, implying that the longitudinal and transversal Hubble rates become identical close to the origin (Plebanski & Krasinski 2006). For Gpc-scale voids with almost constant density profile, however, the approximation H_L ≈ H_T is accurate even for larger radii, in particular in the redshift interval (z< 0.1) in which we determine the effective local Hubble rate (cf. Sect. 3.1). In the region around a central observer, LTB solutions are then practically indistinguishable from simple FLRW models, and Eq. (16)simplifies to $\mathrm{3}\mathit{Ḣ}\mathrm{=}\mathrm{-}\mathrm{3}{\mathit{H}}^{\mathrm{2}}\mathrm{-}\mathrm{4}\mathit{\pi G}\hspace{0.17em}\mathit{\rho }\mathit{.}$(17)Equation (17)describes the time evolution of the Hubble rate as a function of the matter density. Starting from an initial Hubble rate H_i at an early time t_i (e.g. a time t_i at which the Universe was still homogeneous), we can readily integrate Eq. (17)to compute the expansion rate at a later time t_f$\mathit{H}\mathrm{\left(}{\mathit{t}}_{\mathrm{f}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\frac{\mathrm{1}}{{\mathit{H}}_{\mathrm{i}}}\mathrm{+}{\mathrm{\int }}_{{\mathit{t}}_{\mathrm{i}}}^{{\mathit{t}}_{\mathrm{f}}}\mathrm{(}\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{m}}\mathrm{(}{\mathit{t}}^{\mathrm{\prime }}{\mathrm{\left)}}^{\mathrm{\right)}}\mathrm{d}{\mathit{t}}^{\mathrm{\prime }}}\mathit{,}$(18)where we introduced the common matter density parameter ${\mathrm{\Omega }}_{\mathrm{m}}\mathrm{=}\frac{\mathrm{8}\mathit{\pi G\rho }}{\mathrm{3}{\mathit{H}}^{\mathrm{2}}}$. Clearly, the final Hubble rate is maximised by minimising the denominator of Eq. (18), which is achieved by setting Ω_m to zero. The local expansion rate of Gpc-scale voids with slowly varying density profile is thus indeed maximised by vacuum solutions.

6.2. Why not even vacuum LTB solutions exhibit a sufficiently high local Hubble rate

Given the result of the previous section, we can now introduce a vast simplification that helps us explain why not even empty LTB models (that simultaneously fulfil the model-independent constraints from the CMB) exhibit an effective local Hubble rate that comes close to the observationally measured value of H₀ ≈ 73.8 km s^-1 Mpc^-1. We recall that we work in the synchronous time gauge and with a constant bang time function, meaning that the Universe has the same global age t₀ everywhere. We can now envisage a simplified top-hat scenario in which the Universe is globally described by a curved FLRW model, but the observer is located at the centre of an empty sphere (Ω_m = 0, i.e. a Milne model) that was carved out of the homogeneous background. For the purpose of this simplified scenario, technical details concerning the embedding or smooth junction conditions are not relevant. The Hubble function in both regions is given by ${\mathit{H}}^{\mathrm{2}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}{\mathit{H}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{(}{\mathrm{\Omega }}_{\mathrm{m}}{\mathit{a}}^{-3}\mathrm{+}{\mathrm{\Omega }}_{\mathrm{k}}{\mathit{a}}^{-2}\mathrm{)}\mathit{,}$(19)where a denotes the usual scale factor (with a(t₀) = 1) and the curvature parameter is determined by Ω_k = 1 − Ω_m. Integrating Eq. (19)in time yields the standard FLRW relation for the age of the Universe ${\mathit{t}}_{\mathrm{0}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{H}}_{\mathrm{0}}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}\frac{\mathrm{1}}{\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}{\mathit{a}}^{-1}\mathrm{+}{\mathrm{\Omega }}_{\mathrm{k}}}}\hspace{0.17em}\mathrm{d}\mathit{a}\mathit{.}$(20)As this age is assumed to be globally the same, we can directly relate the Hubble rate ${\mathit{H}}_{\mathrm{0}}^{\mathrm{in}}$ inside the empty sphere surrounding the observer to the Hubble rate ${\mathit{H}}_{\mathrm{0}}^{\mathrm{out}}$ of the curved background, ${\mathit{H}}_{\mathrm{0}}^{\mathrm{in}}\mathrm{=}\frac{{\mathit{H}}_{\mathrm{0}}^{\mathrm{out}}}{\mathit{f}\left({\mathrm{\Omega }}_{\mathrm{m}}^{\mathrm{out}}\right)}\mathit{,f}\mathrm{(}{{\mathrm{\Omega }}_{\mathrm{m}}^{\mathrm{out}}}^{\mathrm{)}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}\frac{\mathrm{1}}{\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}^{\mathrm{out}}{\mathit{a}}^{-1}\mathrm{+}\left(\mathrm{1}\mathrm{-}{\mathrm{\Omega }}_{\mathrm{m}}^{\mathrm{out}}\right)}}\mathrm{·}$(21)Figure 8 shows the resulting relation between the Hubble rate ${\mathit{H}}_{\mathrm{0}}^{\mathrm{in}}$ and the matter density parameter ${\mathrm{\Omega }}_{\mathrm{m}}^{\mathrm{out}}$ for a fixed background Hubble rate of ${\mathit{H}}_{\mathrm{0}}^{\mathrm{out}}\mathrm{=}\mathrm{40}\hspace{0.17em}\mathrm{km}\hspace{0.17em}{\mathrm{s}}^{-1}\hspace{0.17em}{\mathrm{Mpc}}^{-1}$, which is the value that was favoured by the CMB data in Sect. 5.3. For instance, if we embed the empty sphere into an Einstein-de Sitter (EdS) background, we find $\mathit{f}\left({\mathrm{\Omega }}_{\mathrm{m}}^{\mathrm{out}}\mathrm{=}\mathrm{1}\right)\mathrm{=}\mathrm{2}\mathit{/}\mathrm{3}$, such that the Hubble rate at the observer position becomes ${\mathit{H}}_{\mathrm{0}}^{\mathrm{in}}\mathrm{=}\mathrm{60}\hspace{0.17em}\mathrm{km}\hspace{0.17em}{\mathrm{s}}^{-1}\hspace{0.17em}{\mathrm{Mpc}}^{-1}$. This agrees well with several previous works, which empirically found that large and deep voids that are embedded into EdS backgrounds cannot increase the local Hubble rate by more than ~ 20 km s^-1 Mpc^-1 (Marra & Notari 2011). Our simple approximation also agrees surprisingly well with the results from Table 2. Most importantly, however, Fig. 8 shows that unrealistically high values for the matter density parameter ${\mathrm{\Omega }}_{\mathrm{m}}^{\mathrm{out}}$ are required to reach the observationally measured local Hubble rate. These theoretical solutions are excluded by current observational data (e.g. constraints from the CMB).

We can thus safely conclude that not even giant LTB voids with radial sizes of several Gpc and vanishing matter density at the origin can simultaneously fit the observed local Hubble rate and the CMB.

6.3. Inhomogeneous Big Bang or cosmological constant?

The decisive property that we exploited in the previous section to link the local Hubble rate around the observer to the Hubble rate of the background model was the assumption of a constant global age throughout the whole Universe. The problem of a too low local Hubble rate can easily be solved by dropping this assumption. To see this very clearly, we used some trivial substitutions to rewrite Eq. (4)in the standard form $\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{B}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{H}}_{\mathrm{0}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}\frac{\mathrm{1}}{\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{a}}^{-3}\mathrm{+}{\mathrm{\Omega }}_{\mathrm{k}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{a}}^{-2}\mathrm{+}{\mathrm{\Omega }}_{\mathrm{\Lambda }}}}\mathrm{d}\mathit{a,}$(22)where the density parameters Ω_x(r) can have an explicit radial dependence in LTB models. Obviously, the local Hubble rate can be increased by decreasing the age of the local Universe. Varying bang time functions that increase towards the origin (r = 0) are thus yet another mechanism to mimic the apparent acceleration of the Universe. In fact, Krasiński (2014) showed that the bang time function can be calibrated such that LTB models exactly reproduce the distance-redshift relation of the standard cosmological model (with dark energy).

These models have previously been discussed in the literature (Clifton et al. 2009; Bull et al. 2012). They were shown to solve the tension between the local Hubble rate and the CMB data. However, Bull et al. (2012) demonstrated that models with the required fluctuations in the bang time function produce a pronounced kSZ effect in galaxy clusters, which is strongly inconsistent with current data.

In addition, we have fundamental objections against LTB models with both flexible density profile and fine-tuned bang time function. As we explained in Sect. 5, the tension between the considered LTB models and current observational data arises as soon as constraints from the CMB are taken into account. In Sect. 2, we briefly discussed that variations in the bang time function can be identified with decaying modes in linear perturbation theory, which in turn would imply (strong) inhomogeneities at early times. However, such features strongly contradict the standard CMB analysis, which is inherently based on the assumption of a spatially isotropic and homogeneous early Universe. We do not believe that it makes sense to introduce a new feature – the varying bang time function – only to reconcile our theoretical model with the CMB, while at the same time this new feature violates the basic assumptions the CMB analysis is based upon. Furthermore, from a statistical point of view, LTB models with varying bang time function become yet more complex, with additional degrees of freedom. Given current observational data, Ockham’s razor would clearly favour the standard cosmological model.

These complications can be partially avoided by extending the considered LTB models in a more natural way. Lovelock (1971, 1972) proved two remarkable, but often overlooked theorems concerning the uniqueness of Einstein’s field equations. Among other things, he showed that under very general simplicity conditions and in four dimensions, any metric theory of gravity locally conserving energy-momentum must have two coupling constants, and its metric must satisfy Einstein’s field equations, ${\mathit{G}}_{\mathit{\mu \nu }}\mathrm{=}\mathit{\kappa }{\mathit{T}}_{\mathit{\mu \nu }}\mathrm{+}\mathrm{\Lambda }{\mathit{g}}_{\mathit{\mu \nu }}\mathit{,}$(23)where G_μν, T_μν and g_μν are, respectively, the Einstein tensor, the energy-momentum tensor, and the metric tensor. κ and Λ are the two coupling constants, which we identify with Newton’s gravitational constant and with the cosmological constant.

According to Lovelock’s theorems, the cosmological constant Λ thus appears naturally not only in Einstein’s field equations, but also in the field equations of any other metric theory of gravity. It could (and perhaps should) therefore be interpreted as a second coupling constant of the metric to matter, on a par with Newton’s gravitational constant. Numerous alternative cosmological models – such as LTB voids – have been constructed to avoid the cosmological constant or dark energy. Nonetheless, from the point of view of Lovelock’s theorem, such models must be considered incomplete unless they justify why the cosmological constant should vanish.

Following this line of reasoning, we believe that considering LTB models with non-zero cosmological constant is the most natural extension of our approach. This is an interesting step because it slightly shifts the research focus away from the dark energy problem to more general tests of the Copernican principle. For example, using flexible parametrisations of the local matter density profile, ΛLTB models allow us to derive statistical constraints on possible deviations from spatial homogeneity. In addition, effects of varying spatial curvature along the past-null cone can be explored. Finally, we can investigate whether fluctuations of the local matter density profile can noticeably influence the statistical inference of the best-fitting parameters of the standard cosmological model (see e.g. Marra et al. 2013b; Valkenburg et al. 2014). We briefly discuss these applications of ΛLTB models in the next section.

7. Probing spatial homogeneity with ΛLTB models

In this section, we constrain the cosmological models considered by all data sets mentioned in Sect. 3, that is, we compute the likelihood given the local Hubble rate, supernovae, model-independent CMB constraints and kSZ data. Again, as a point of reference, we first detail the goodness-of-fit of the standard cosmological model, which is usually assumed to be a spatially flat FLRW model with non-zero cosmological constant. The best-fitting model with cosmological parameters (h,Ω_m,Ω_Λ) = (0.72,0.263,0.737) exhibits a log-likelihood value of log (ℒ) = −275.1, with the individual contributions log [ℒ_H0(1)] ≈ −0.3, log [ℒ_SNe(580)] ≈ −272.9, log [ℒ_CMB(3)] ≈ −0.3 and log [ℒ_kSZ(9)] ≈ −1.6, where we explicitly indicated the number of fitted data points in parenthesis. It is remarkable how well the standard cosmological model agrees with these data. The goodness-of-fit, which can be quantified in terms of the chi-squared per degree of freedom, is truly impressive and hard to improve.

On the other hand, it is well-known that – even within the framework of the standard cosmological model – there seems to be a slight tension concerning the Hubble rate: while measurements with Cepheid-calibrated supernovae yield a local Hubble rate of ${\mathit{H}}_{\mathrm{0}}\mathrm{=}{\mathrm{(}}^{}\mathrm{73.8}\mathrm{\pm }{\mathrm{2.4}}^{\mathrm{)}}\hspace{0.17em}\mathrm{km}\hspace{0.17em}{\mathrm{s}}^{-1}\hspace{0.17em}{\mathrm{Mpc}}^{-1}$, the Planck CMB data alone favour a value of H₀ = (67.3 ± 1.2) km s^-1 Mpc^-1 (Riess et al. 2011; Planck Collaboration XVI 2014). This discrepancy is still debated, but it could be taken as a further motivation to explore more complicated, radially inhomogeneous ΛLTB models with spatially varying Hubble rates (Marra et al. 2013a). Moreover, we can address questions concerning the assumption of spatial homogeneity, such as:

• Can inhomogeneous cosmological models fit the given data even more accurately?

• What limits do current data impose on fluctuations of the local matter density profile? How strongly can we deviate in our assumptions from spatial homogeneity?

• How do the constraints on dark energy or spatial curvature change if we drop the assumption of spatial homogeneity and marginalise over all possible radial inhomogeneities?

As discussed in Sect. 5.2, the CMB data do not constrain the detailed shape of the local matter density profile, but primarily the global properties of the background model and the spatial curvature close to the LSS. The kSZ data constrain the allowed depth of radial inhomogeneities, but since we assumed large errors for the individual measurements, current constraints are rather weak, at least when considering only moderate deviations from homogeneity. We would therefore expect that fine-tuned local matter density profiles can mainly improve the fit to the local Hubble rate and supernovae. As the constraining power of supernovae decreases with increasing redshift (cf. Sect. 5.1), it only makes sense to consider fluctuations well inside the redshift range z< 1.5.

For these reasons, and also for a comparison with the previous sections, we now exemplarily consider ΛLTB models whose density profile is represented by cubic splines with three flexible nodes and a fixed radial size of L = 3 Gpc. While testing, we also tried to treat the radial size L as a free parameter and vary it with the Monte Carlo sampler, but the observational constraints are too weak so that the Monte Carlo chains did not converge. We asymptotically embedded these ΛLTB models into curved FLRW backgrounds.

In summary, the considered ΛLTB models are described by the cosmological parameters h, Ω_m, Ω_Λ, Ω_k = 1 − Ω_m − Ω_Λ of the background, and the three spline nodes ${}^{\mathrm{(}}\mathit{a}_{\mathrm{1}}\mathit{,}{\mathit{a}}_{\mathrm{2}}\mathit{,}{\mathit{a}}_{\mathrm{3}}{\mathrm{)}}^{}$ parametrising possible fluctuations of the matter density profile. In contrast to the flat FLRW model, we hence introduced one new degree of freedom for the cosmological background (Ω_Λ free instead of setting Ω_Λ = 1 − Ω_m) and three additional parameters for radial inhomogeneities. This way, our framework for testing the cosmological principle is most general and able to parametrise a broad class of space-time geometries.

We now sequentially discuss these questions. The best-fitting ΛLTB model has a log-likelihood value of log (ℒ) = −274.7, which is only marginally better than that of the FLRW model. The fits to all observables slightly improve. There is no distinguished data set that is described significantly better.

As an example, we show a comparison of the marginalised posterior distribution of the effective local Hubble rates in Fig. 9. As expected, the more flexible curved ΛLTB models fit the observed value better than the flat FLRW models.

More surprisingly, the best-fitting model is almost perfectly homogeneous, with spline nodes a₁ = 0.99, a₂ = 1.01 and a₃ = 0.99. The cosmological parameters are also very similar to those of the standard model, with h = 0.73, Ω_m = 0.25, Ω_Λ = 0.74 and Ω_k = 0.01.

In summary, we can conclude that, of course, the more flexible models fit the data even better. The improvement is almost negligible (Δlog (ℒ) ~ 0.4) , however, and comes at the high price of introducing four additional free parameters. Occam’s razor penalises more complicated models and can be approximated for instance by the Akaike information criterion or the Bayesian information criterion (Akaike 1974; Schwarz 1978). In our case, these criteria suggest that not even one additional free parameter would be justified. More meaningfully, we can thus conclude that the data used in this work statistically favour the standard cosmological model.

The second question concerns statistical constraints on deviations from spatial homogeneity. The mean values of the spline parameters are given by a₁ = 1.02 ± 0.12, a₂ = 1.02 ± 0.07 and a₃ = 0.96 ± 0.11. The corresponding variations of the density profiles are depicted in Fig. 10. The density profiles were forced to converge to the background value at r = 3 Gpc by construction, so the apparent constraints at r> 2 Gpc are artificial. Clearly, the results agree well with the assumption of spatial homogeneity. However, the scatter of ~ 15% indicates that we need more accurate data to safely confirm the cosmological principle. This demonstrates the importance of fully solving linear perturbation theory in LTB backgrounds (such that it can be applied to predict cosmological observables), since this would allow us to include more cosmological observables.

For instance, it would be natural to expect that Gpc-scale variations of the matter density profile leave a characteristic imprint on the large-scale structure of the local Universe. Such imprints could create non-vanishing amplitudes of the galaxy-galaxy correlation function on Gpc-scales, which would be at odds with data from current spectroscopic surveys (Labini 2011; Scrimgeour et al. 2012). Furthermore, Gpc-variations of the local matter distribution should be well measurable with tomographic weak-lensing methods (Schäfer & Heisenberg 2012). There are other promising observables related to cosmic structures that might help constrain deviations from spatial homogeneity, but we need to advance the numerical algorithms for solving the linear perturbation equations on LTB backgrounds before we are able to reliably calculate these phenomena.

To answer the last question, in Figs. 11 and 12 we show the statistical constraints on the cosmological constant and spatial curvature. It is important to stress that these posterior distributions were marginalised over all possible fluctuations of the matter density profile. In other words, these posteriors show constraints that are independent of the assumption of spatial homogeneity, at least within the limits of our framework. The constraints on the cosmological constant are essentially invariant, with a negligibly larger scatter. This result differs from the findings of Valkenburg et al. (2014), who calculated a shift of ΔΩ_Λ ~ 0.05 when marginalising over inhomogeneities. However, these authors constrained their models by different data, chose another approach for modelling spatial variations of the density profile and additionally included a so-called Copernican prior. Amongst other details, these are important issues that could well explain the difference with respect to our results. Finally, it is remarkable how tightly the spatial curvature of the background model is constrained, even though we considerably relaxed our prior assumptions. The constraints are consistent with a flat background Universe, which agrees very well with the Planck results (Planck Collaboration XVI 2014).

8. Conclusions

What do cosmological observations tell us about the shape of the local matter density profile? Can current data confirm the cosmological principle? These were two central questions that we discussed in this work. To address these problems, we modelled the Universe around us by spherically symmetric, but radially inhomogeneous LTB models with a synchronous Big Bang. In contrast to most previous studies, we did not model the mass profile of LTB models in terms of empirically motivated functional forms, but instead chose to represent the local matter density profile by flexible interpolation schemes or a decomposition into Laguerre polynomials. We chose this alternative, more flexible approach to be able to investigate a broader class of problems (e.g. not restricted to void scenarios) and, more importantly, to not bias our results by prior model assumptions.

In the first part of this work, we reconsidered LTB models without cosmological constant to investigate whether highly flexible, radially fine-tuned mass profiles allow us to simultaneously fit the high local Hubble rate and the CMB data from the Planck satellite (Riess et al. 2011; Planck Collaboration XVI 2014). To this end, we consecutively compared numerous LTB models with different combinations of observational data. The main results of this first part can be summarised as follows:

• 1.

  H₀+ supernovae: as was already well-known, LTB models without cosmological constant can easily fit the observed local Hubble rate and supernovae just as well as the standard cosmological model. These models mimic the apparent acceleration of the Universe by means of large Gpc-scale voids whose matter density profile gradually decreases towards the observer. Having said that, we would like to stress an important issue concerning the interpretation of supernova observations. Interpreted in the framework of a spatially homogeneous FLRW model, supernovae favour a non-zero cosmological constant. Similarly, supernovae favour Gpc-scale LTB voids if and only if we set the cosmological constant to zero and require a synchronous Big Bang. In fact, Célérier et al. (2010) demonstrated that local over-densities (i.e. giant local humps) are favoured if the bang time function is allowed to vary. We mention these ambiguities only to emphasise that the interpretation of supernovae is strongly biased by prior assumptions about the underlying cosmological model.

• 2.

  Cosmic microwave background: the CMB data alone do not constrain the shape of the local matter density profile. Indeed, the density profiles can substantially be varied without notably deteriorating the fit to the data. However, current CMB data impose tight constraints on the geometry and expansion rate of the asymptotic background models. Without a cosmological constant, good fits to the data require asymptotically curved FLRW backgrounds with an unrealistically low Hubble rate of H₀ ≈ 33 km s^-1 Mpc^-1.

• 3.

  H₀+ supernovae + CMB: the separate CMB analysis clearly highlighted the problem of the considered class of LTB models: a good fit to the CMB requires an extremely low background Hubble rate, which is in strong tension with the observed local Hubble rate of H₀ = (73.8 ± 2.4) km s^-1 Mpc^-1 (Riess et al. 2011). In the last part of our analysis, we thus focused on the question whether or not radially tuned matter density profiles can be found that simultaneously comply with all observational constraints. However, even though we considered a wide variety of huge, extremely deep and heavily fine-tuned void profiles, the effective local Hubble rate remained too low, H₀ ≲ 62 km s^-1 Mpc^-1. To make matters worse, we note that we merely used minimal, model-independent constraints from the CMB data to arrive at these conclusions (see Sect. 3.3, for details). We did not even use the complete information contained in the full CMB spectra, which appear to be even more problematic for void models (Moss et al. 2011).

After this detailed analysis, we presented simple theoretical arguments that explain why not even heavy fine-tuning of the radial matter density profile can solve the tension between the local Hubble rate and the CMB. Various solutions, such as varying bang time functions, dynamical effects of radiation, or modifications of the primordial curvature power-spectrum, were proposed to reconcile giant void scenarios with current data (Clarkson & Regis 2011; Nadathur & Sarkar 2011; Bull et al. 2012). However, all these modifications require deviations from the standard inflationary paradigm or introduce more complexity, resulting in fine-tuned, alternative cosmological models. Ockham’s razor would clearly favour the standard cosmological model. In addition, according to Lovelock’s theorems, the cosmological constant appears naturally in the field equations of general relativity and of any other metric theory of gravity (Lovelock 1971, 1972). We thus argued that considering LTB models with non-zero cosmological constant seems to be the most natural step.

In the final part, we therefore analysed LTB models with non-zero cosmological constant, which can be considered as the simplest, spatially isotropic, but radially inhomogeneous extension of the standard cosmological model. ΛLTB models are a valuable framework to systematically study deviations from spatial homogeneity, verify or falsify the cosmological principle, or simply explore effects of varying curvature along the past null cone. We showed that current data provide no evidence for radial inhomogeneities on Gpc-scales. Instead, spatially flat FLRW models with homogeneous matter distribution are favoured. These results statistically support the cosmological principle. However, we also showed that fluctuations of ~ 15% with respect to a homogeneous matter density profile are still compatible with current data.

As emphasised in the main text (cf. Sects. 3 and 7), our analysis was limited by an important drawback: linear perturbation theory in LTB models is substantially more complicated than in FLRW models, mainly because scalar, vectorial, and tensorial perturbations do not decouple on inhomogeneous backgrounds. Although great progress has been made, substantial additional efforts are required before linear perturbations can reliably be computed in realistic cosmological settings. We therefore had to neglect all cosmological observables that depend on the details of linear structure formation. For instance, we were not yet able to calculate statistical properties of the perturbed matter density in the LTB models considered, meaning that we could not yet compare these models with the statistical properties of the observed large-scale structure (e.g. galaxy-galaxy correlation functions, and cluster number counts). Such tests should certainly help to constrain the underlying space-time geometry. We also had to neglect important cosmological probes such as baryonic acoustic oscillations and weak-lensing spectra, which are widely (and successfully) used to constrain homogeneous and isotropic cosmologies. We will focus our research on advancing linear perturbation theory in LTB models to tighten observational constraints on the shape of the local matter density profile in future works.

Acknowledgments

M.R. thanks the Sydney Institute for Astronomy (SIfA, University of Sydney) for hospitality and the German Academic Exchange Service (DAAD, Doktoranden-Stipendium) for the financial support during a research visit. This work was supported in part by the Deutsche Forschungsgemeinschaft (DFG). Most simulations required for this work were performed on the bwGRiD cluster (http://www.bw-grid.de), member of the German D-Grid initiative, funded by the Ministry for Education and Research (Bundesministerium für Bildung und Forschung) and the Ministry for Science, Research and Arts Baden-Wuerttemberg (Ministerium für Wissenschaft, Forschung und Kunst Baden-Württemberg).

References

Appendix A: Parametric solutions for the LTB model without cosmological constant

In this appendix, we amend some more useful relations for LTB models without cosmological constant that we skipped in the main part to improve the readability of the paper.

To begin with, we repeat Eq. (4)(with Λ = 0) from Sect. 2, which results from integrating one of Einstein’s field equations in time: ${\mathit{t}}_{\mathrm{0}}\mathrm{-}{\mathit{t}}_{\mathrm{B}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\underset{\mathrm{0}}{\overset{\mathit{R}\mathrm{\left(}{\mathit{t}}_{\mathrm{0}}\mathit{,r}\mathrm{\right)}}{\mathrm{\int }}}\frac{\mathrm{1}}{\sqrt{\frac{\mathrm{2}\mathit{M}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}}\mathrm{+}\mathrm{2}\mathit{E}\mathrm{\left(}\mathit{r}\mathrm{\right)}}}\hspace{0.17em}\mathrm{d}\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\mathit{.}$(A.1)Dependent on the sign of the curvature function E(r), this integral has three different parametric solutions:

• Elliptic evolution: E(r) < 0$\begin{array}{ccc}& & \mathit{R}\mathrm{\left(}\mathit{t,r}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathit{M}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{2}\mathit{E}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathrm{(}\mathrm{1}\mathrm{-}\cos \mathit{\eta }\mathrm{)}\mathit{,}\\ & & \mathit{\eta }\mathrm{-}\sin \mathit{\eta }\mathrm{=}\frac{{\left[\mathrm{-}\mathrm{2}\mathit{E}\mathrm{\left(}\mathit{r}\mathrm{\right)}\right]}^{\mathrm{3}\mathit{/}\mathrm{2}}}{\mathit{M}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\left[\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{B}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\right]\mathit{.}\end{array}$

• Parabolic evolution: E(r) = 0$\mathit{R}\mathrm{\left(}\mathit{t,r}\mathrm{\right)}\mathrm{=}{\left\{\frac{\mathrm{9}}{\mathrm{2}}\mathit{M}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\left[\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{B}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\right]}^{\mathrm{2}}\right\}}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{.}$(A.4)

• Hyperbolic evolution: E(r) > 0$\begin{array}{ccc}& & \mathit{R}\mathrm{\left(}\mathit{t,r}\mathrm{\right)}\mathrm{=}\frac{\mathit{M}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{2}\mathit{E}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathrm{(}\cosh \mathit{\eta }\mathrm{-}\mathrm{1}\mathrm{)}\mathit{,}\\ & & \sinh \mathit{\eta }\mathrm{-}\mathit{\eta }\mathrm{=}\frac{{\left[\mathrm{2}\mathit{E}\mathrm{\left(}\mathit{r}\mathrm{\right)}\right]}^{\mathrm{3}\mathit{/}\mathrm{2}}}{\mathit{M}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\left[\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{B}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\right]\mathit{.}\end{array}$

As described in Sect. 4.2, we determine the curvature function E(r) as a function of time t, effective gravitational mass M(r), and areal radius R(t,r). For this purpose, it is useful to rewrite the parametric solutions as $\begin{array}{ccc}& & \mathit{t}\mathrm{=}{\mathit{t}}_{\mathrm{B}}\mathrm{+}\frac{\mathit{M}}{\mathrm{(}\mathrm{-}\mathrm{2}\mathit{E}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\left[\arccos \left(\mathrm{1}\mathrm{+}\frac{\mathrm{2}\mathit{ER}}{\mathit{M}}\right)\mathrm{-}\mathrm{2}\sqrt{\frac{\mathrm{-}\mathit{ER}}{\mathit{M}}\left(\mathrm{1}\mathrm{+}\frac{\mathit{ER}}{\mathit{M}}\right)}\hspace{0.17em}\right]\\ & & \mathrm{0}\mathrm{\le }\mathit{\eta }\mathrm{\le }\mathit{\pi ,}\end{array}$(A.7)$\begin{array}{ccc}& & \mathit{t}\mathrm{=}{\mathit{t}}_{\mathrm{B}}\mathrm{+}\frac{\mathit{M}}{\mathrm{(}\mathrm{-}\mathrm{2}\mathit{E}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\left[\mathit{\pi }\mathrm{+}\arccos \left(\mathrm{-}\mathrm{1}\mathrm{-}\frac{\mathrm{2}\mathit{ER}}{\mathit{M}}\right)\mathrm{+}\mathrm{2}\sqrt{\frac{\mathrm{-}\mathit{ER}}{\mathit{M}}\left(\mathrm{1}\mathrm{+}\frac{\mathit{ER}}{\mathit{M}}\right)}\hspace{0.17em}\right]\mathit{,}\\ & & \mathit{\pi }\mathrm{\le }\mathit{\eta }\mathrm{\le }\mathrm{2}\mathit{\pi ,}\end{array}$(A.8)for the expanding and collapsing elliptic cases, and $\mathit{t}\mathrm{=}{\mathit{t}}_{\mathrm{B}}\mathrm{+}\frac{\mathit{M}}{\mathrm{(}\mathrm{2}\mathit{E}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\left[\mathrm{2}\sqrt{\frac{\mathit{ER}}{\mathit{M}}\left(\mathrm{1}\mathrm{+}\frac{\mathit{ER}}{\mathit{M}}\right)}\mathrm{-}\mathrm{arcosh}\left(\mathrm{1}\mathrm{+}\frac{\mathrm{2}\mathit{ER}}{\mathit{M}}\right)\right]$(A.9)for the hyperbolic case. We refer for example to Plebanski & Krasinski (2006) for more details on these solutions. Equations (A.7)–(A.9)can be numerically unstable in the near parabolic limit (| E | ≪ 1). In this case, we use an inverse series expansion as explained in Appendix B of Hellaby & Krasiński (2006).

Finally, we note that the parametric solutions (A.2)–(A.6)can be combined to derive an analytic expression for the radial derivative of the areal radius function, ${\mathit{R}}^{\mathrm{\prime }}\mathrm{=}\left(\frac{{\mathit{M}}^{\mathrm{\prime }}}{\mathit{M}}\mathrm{-}\frac{{\mathit{E}}^{\mathrm{\prime }}}{\mathit{E}}\right)\mathit{R}\mathrm{+}\left[\left(\frac{\mathrm{3}}{\mathrm{2}}\frac{{\mathit{E}}^{\mathrm{\prime }}}{\mathit{E}}\mathrm{-}\frac{{\mathit{M}}^{\mathrm{\prime }}}{\mathit{M}}\right)\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{B}}\mathrm{)}\mathrm{-}{\mathit{t}}_{\mathrm{B}}^{\mathrm{\prime }}\right]\mathit{Ṙ}\mathit{.}$(A.10)This relation, together with the radial derivative of Eq. (2), $\mathit{Ṙ}\mathrm{\prime }\mathrm{=}\frac{\mathrm{1}}{\mathit{Ṙ}}\hspace{0.17em}\left[\left(\frac{{\mathit{M}}^{\mathrm{\prime }}\mathit{R}\mathrm{-}\mathit{M}{\mathit{R}}^{\mathrm{\prime }}}{{\mathit{R}}^{\mathrm{2}}}\right)\mathrm{+}{\mathit{E}}^{\mathrm{\prime }}\right]\mathit{,}$(A.11)allows us to analytically compute R′ and Ṙ′, which renders the numerical integration of radial null geodesics more efficient (cf. Eqs. (5)and (6)).

All Tables

All Figures