1 Introduction

The standard cosmological model is based on the assumption that the Universe is homogeneous and isotropic on sufficiently large scales. Nevertheless local observations could be strongly affected by local structure as shown for example in [1], and it is important to study its effects. The analysis of luminosity density data [2] has provided some strong experimental evidence supporting the existence of local inhomogeneities, but it would be important to confirm it using another observable such as the baryonic acoustic oscillations (BAO) measurements [3–8]. The BAO scale allows in fact to determine the expansion rate of the Universe H(z) independently from the luminosity distance and as such provides an important source of information as regards our Universe. If the H(z) estimations obtained from BAO observations data will show deviations from the \(\Lambda CDM\) predictions this could be considered an independent evidence of the existence of local inhomogeneities. This motivates the calculation of a low redshift formula for H(z), able to take into account the effects of inhomogeneities which cannot be fully modeled with perturbation theory, as some of the inhomogeneities found for example in [2]. A low-redshift expansion based on the use of exact solutions of Einstein’s equations is in fact valid also for large values of the density contrast or of the gravitational potential.

The effects of a local inhomogeneity on cosmological observations have been studied already for different cases [1, 9–30] such as the equation of state of dark energy or the luminosity distance [18, 19, 23]. It has been shown for example that the value of the cosmological constant could be affected significantly by the presence of local inhomogeneity seeded by primordial curvature perturbations [9], which could also lead to the wrong conclusion of a varying equation of state for dark energy while only a cosmological constant is present [19]. The origin of these effects is that spatial inhomogeneities can change the energy of the propagating photons, contaminating the cosmological redshift due to the Universe expansion and consequently introducing some errors in the estimation of parameters based on cosmological models which ignore the effects of the inhomogeneities.

In this paper we will focus on the low-redshift effects of inhomogeneities on the Hubble expansion parameter, adopting an analytical approach based on the use of an exact solution of Einstein’s equations to model the local structure. We first derive the redshift expansion of the geodesics equations and use it to obtain the expansion of H(z). We then compute a formula for \(\frac{\delta H}{H}\), the relative difference between the \(\Lambda CDM\) and the inhomogeneous case. Finally we compare the formulas with the numerical calculations based on the integration of the Einstein’s equation and the geodesics equations, finding good agreement. We also check that the low-redshift expansion formulas are more precise than the perturbative calculation, especially when the density contrast is larger.

2 Modeling the local Universe

We use the LTB solution to model the monopole component of the local structure [31–35]

$$\begin{aligned} \mathrm{d}s^2 = -\mathrm{d}t^2 + \frac{(R'(t,r))^2 \mathrm{d}r^2}{1 + 2\,E(r)}+R(t,r)^2 \mathrm{d}\Omega ^2, \end{aligned}$$
(1)

where R is a function of the time coordinate t and the radial coordinate r, E(r) is an arbitrary function of r, and \(R'(t,r)=\partial _rR(t,r)\). The Einstein equations imply that

$$\begin{aligned} \left( {\frac{\dot{R}}{R}}\right) ^2= & {} \frac{2 E(r)}{R^2}+\frac{2M(r)}{R^3}+\frac{\Lambda }{3},\end{aligned}$$
(2)
$$\begin{aligned} \rho (t,r)= & {} \frac{2M'}{R^2 R'}, \end{aligned}$$
(3)

where M(r) is an arbitrary function of r, \(\dot{R}=\partial _tR(t,r)\), and we choose a system of units in which \(c=8\pi G=1\).

To compute H(z) we need to solve the radial null geodesics [36]

$$\begin{aligned} {\mathrm{d}r\over \mathrm{d}z}= & {} {\sqrt{1+2E(r(z))}\over {(1+z)\dot{R'}[t(z),r(z)]}}, \end{aligned}$$
(4)
$$\begin{aligned} {\mathrm{d}t\over \mathrm{d}z}= & {} -\,{R'[t(z),r(z)]\over {(1+z)\dot{R'}[t(z),r(z)]}}, \end{aligned}$$
(5)

and then we substitute in the formula for the Hubble parameter in a LTB space [37, 38]

$$\begin{aligned} H(t,r)= & {} \frac{2}{3} H_{\perp }(t,r) + \frac{1}{3} H_{\Vert }(t,r),\end{aligned}$$
(6)
$$\begin{aligned} H(z)= & {} H(t(z),r(z)), \end{aligned}$$
(7)

where

$$\begin{aligned} H_{\perp }(t,r)\equiv \frac{\dot{R}(t,r)}{R(t,r)},\end{aligned}$$
(8)
$$\begin{aligned} H_{\Vert }(t,r)\equiv \frac{\dot{R}' (t,r)}{R' (t,r)}. \end{aligned}$$
(9)

The analytical solution can be derived [39, 40] introducing a new coordinate \(\eta =\eta (t,r)\), and new functions \(\rho _0(r)\) and k(r) given by

$$\begin{aligned} \frac{\partial \eta }{\partial t}\Big |_r= & {} \frac{r}{R}=\frac{1}{a},\end{aligned}$$
(10)
$$\begin{aligned} \rho _0(r)= & {} \frac{6 M(r)}{r^3},\end{aligned}$$
(11)
$$\begin{aligned} k(r)= & {} -\frac{2E(r)}{r^2}. \end{aligned}$$
(12)

We will adopt, without loss of generality, the coordinate system in which \(\rho _0(r)\) is a constant, the so called FLRW gauge. We can then express Eq. (2) in the form

$$\begin{aligned} \left( \frac{\partial a}{\partial \eta }\right) ^2= -k(r)a^2 + \frac{\rho _0}{3}a + \frac{\Lambda }{3}a^4. \end{aligned}$$
(13)

The coordinate \(\eta \), which can be considered a generalization of the conformal time in a homogeneous FLRW Universe, is defined implicitly by Eq. (10). The relation between t and \(\eta \) is obtained by integrating Eq. (10) and is given by [26]

$$\begin{aligned} t(\eta ,r)=\int _{0}^{\eta } a(x,r) \, \mathrm{d}x +t_b(r), \end{aligned}$$
(14)

where \(t_b(r)\) is a functional constant of integration called the bang function, since it corresponds to the fact that in these models the scalar factor can vanish at different times at different locations. We will consider models with \(t_b(r)=0\). The solution of Eq. (13) can then be written in the form

$$\begin{aligned} a(\eta ,r)=\frac{\rho _0}{k(r)+3 \wp (\frac{\eta }{2};g_2(r),g_3(r))}, \end{aligned}$$
(15)

where \(\wp (x;g_2,g_3)\) is the Weierstrass elliptic function and

$$\begin{aligned} g_2(r)=\frac{4}{3}k(r)^2,\quad g_3(r)=\frac{4}{27} (2k(r)^3 -\Lambda \rho _0^2). \end{aligned}$$
(16)

In terms of \(\eta \) and \(a(\eta ,r)\) the radial null geodesics and the Hubble parameter are given by [12]

$$\begin{aligned} \frac{\mathrm{d}\eta }{\mathrm{d}z}= & {} -\frac{\partial _r t(\eta ,r) + G(\eta ,r)}{(1+z)\partial _\eta G(\eta ,r)},\end{aligned}$$
(17)
$$\begin{aligned} \frac{\mathrm{d}r}{\mathrm{d}z}= & {} \frac{a(\eta ,r)}{(1+z)\partial _{\eta } G(\eta ,r)},\end{aligned}$$
(18)
$$\begin{aligned} H(\eta ,r)= & {} H(t(\eta ,r),r), \end{aligned}$$
(19)

where

$$\begin{aligned} G(\eta ,r)\equiv & {} \frac{R,_r}{\sqrt{1+2E(r)}}\nonumber \\= & {} \frac{[\partial _r (a(\eta ,r)r) - a^{-1} \partial _\eta (a(\eta ,r)r) \partial _r t(\eta ,r)]}{\sqrt{1-k(r)r^2}}.\nonumber \\ \end{aligned}$$
(20)

The function \(G(\eta ,r)\) has an explicit analytical form, making the above geodesics equations particularly suitable for a low-redshift expansion.

3 Low-redshift expansion of the Hubble parameter H(z)

In order to obtain a low-redshift formula for the Hubble parameter we expand the function k(r) as

$$\begin{aligned} k(r)=k_0 + k_1 r + k_2 r^{2} +\cdots , \end{aligned}$$
(21)

We also need an expansion for \(t(\eta ,r)\), which can be obtained from the exact solution for \(a(\eta ,r)\) according to

$$\begin{aligned} t(\eta ,r)= & {} t_0(r)+a(\eta _0,r)(\eta -\eta _0)\nonumber \\&+\,\frac{1}{2} \partial _{\eta }a(\eta _0,r)(\eta -\eta _o)^2+\cdots , \end{aligned}$$
(22)

where we defined the function \(t_0(r)\) by

$$\begin{aligned} t_0(r)\equiv t(\eta _0,r). \end{aligned}$$
(23)

Using the properties of the Weierstrass elliptic functions \(\wp \), \(\zeta \), and \(\sigma \) [41] we can compute the integral in Eq. (14), obtaining

$$\begin{aligned} t(\eta ,r)= & {} \frac{2 \rho _0}{3 \wp '\left( \wp ^{-1} \left( {-\frac{k(r)}{3}} \right) \right) } \left[ \ln \left( \frac{\sigma \left( \frac{\eta }{2}-\wp ^{-1}\left( -\frac{k(r)}{3}\right) \right) }{\sigma \left( \frac{\eta }{2}+\wp ^{-1}\left( -\frac{k(r)}{3}\right) \right) }\right) \right. \nonumber \\&\left. +\,\eta \zeta \left( \wp ^{-1}\left( -\frac{k(r)}{3}\right) \right) \phantom {\left( \frac{\sigma \left( \frac{\eta }{2}-\wp ^{-1}\left( -\frac{k(r)}{3}\right) \right) }{\sigma \left( \frac{\eta }{2}+\wp ^{-1}\left( -\frac{k(r)}{3}\right) \right) }\right) }\right] , \end{aligned}$$
(24)

where \(\wp '\) is the derivative of the Weierstrass elliptic function \(\wp \). \(\wp ^{-1}\), \(\zeta \), and \(\sigma \) are defined by the equations

$$\begin{aligned} \wp ^{-1}\left( \wp \left( x\right) \right)= & {} x, \end{aligned}$$
(25)
$$\begin{aligned} \zeta ' \left( x\right)= & {} - \wp \left( x \right) ,\end{aligned}$$
(26)
$$\begin{aligned} \frac{\sigma ' \left( x \right) }{\sigma \left( x \right) }= & {} \zeta \left( x \right) . \end{aligned}$$
(27)

The use of the exact expression for \(t(\eta ,r)\) improves the accuracy for the expansion for the geodesics with respect to previous calculations [23], which were based on a perturbative expansion of \(t_0(r)\), rather than the use of the exact value.

Now we can find the low redshift Taylor expansion for the geodesic equations [19], and then we calculate the Hubble parameter. We expand the solution of the geodesic equations according to

$$\begin{aligned} r(z)= & {} r_1z+r_2z^2+r_3z^3 +\cdots ,\end{aligned}$$
(28)
$$\begin{aligned} \eta (z)= & {} \eta _0 +\eta _1 z +\eta _2 z^2+\cdots . \end{aligned}$$
(29)

After substituting the above expansion in the geodesic equations we can map the solution of the system of differential equations into the solution of a system of algebraic equations for the coefficients of the expansion. Here we give the formulas for the case in which \(k_{0}=0\), which is enough to understand qualitatively the effects of the inhomogeneity. The term \(k_{0}\) corresponds in fact to the homogeneous component of the curvature function, which in the absence of inhomogeneities is simply the curvature of a FLRW model and as such is not associated to any new physical effect not already known from standard cosmology.

For the geodesics we get

$$\begin{aligned} \eta _1= & {} -\frac{1}{a_0^2 H_0} [a_0+ t_0'(0)],\end{aligned}$$
(30)
$$\begin{aligned} \eta _2= & {} \frac{1}{12 a_0^3 H_0^2 \Omega _M} [ a_0 H_0 t_0'(0) (3 \Omega _M (9 \Omega _M-4)-8 K_1)\nonumber \\&+\,a_0^2 H_0 (9 \Omega _M^2-4 K_1)\nonumber \\&+\,6 \Omega _M (3 H_0 (\Omega _M-1) t_0'(0)^2-t_0''(0))],\end{aligned}$$
(31)
$$\begin{aligned} r_1= & {} \frac{1}{a_0 H_0},\end{aligned}$$
(32)
$$\begin{aligned} r_2= & {} \frac{1}{12 a_0^2 H_0 \Omega _M} [a_0(4 K_1-9 \Omega _M^2)\nonumber \\&+\,6(2-3 \Omega _M)\Omega _M t_0'(0)],\end{aligned}$$
(33)
$$\begin{aligned} r_3= & {} \frac{1}{72 a_0^3 H_0^2 \Omega _{\Lambda } \Omega _M^2} [ a_0^2 H_0(4 K_1^2(2 \Omega _{\Lambda }\nonumber \\&+\,\zeta _0(2-3 \Omega _M)+\Omega _M)-60 K_1 \Omega _{\Lambda } \Omega _M^2\nonumber \\&+\,3 \Omega _{\Lambda } \Omega _M (8 K_2 +3 (9 \Omega _M-4)\Omega _M^2))\nonumber \\&-\,36 a_0 H_0 \Omega _{\Lambda } \Omega _M \text {t0}'(0) (K_1 (4 \Omega _M-2)\nonumber \\&+\,(8-9 \Omega _M)\Omega _M^2)+ 18 \Omega _{\Lambda } \Omega _M^2\nonumber \\&\times \,(6 H_0 (3 \Omega _M^2-4 \Omega _M+1) t_0'(0)^2+(2-3 \Omega _M) t_0''(0))],\nonumber \\ \end{aligned}$$
(34)

where \(\Omega _M\), \(\Omega _{\Lambda }\), \(T_0\), and \(K_n\) are dimensionless quantities given by [16]

$$\begin{aligned} \rho _0= & {} 3 \Omega _M a_0^3 H_0^2, \end{aligned}$$
(35)
$$\begin{aligned} \Lambda= & {} 3 \Omega _{\Lambda } H_0^2,\end{aligned}$$
(36)
$$\begin{aligned} T_0= & {} \eta _0 \left( a_0 H_0 \right) ,\end{aligned}$$
(37)
$$\begin{aligned} K_n= & {} k_n (a_0 H_0)^{n+2}. \end{aligned}$$
(38)

We have also used the following definitions:

$$\begin{aligned} a_0= & {} a (\eta _0,0),\end{aligned}$$
(39)
$$\begin{aligned} H_0= & {} H(\eta _0,0),\end{aligned}$$
(40)
$$\begin{aligned} \zeta _0= & {} \zeta \left( \frac{T_0}{2};\frac{4 K_0^2}{3},\frac{4}{27}\left( 2 K_0^3 - 27 \Omega _{\Lambda } \Omega _M^2\right) \right) , \end{aligned}$$
(41)

where \(\zeta \) is the Weierstrass zeta function. As we can see in the above formulas the effects of the inhomogeneity start to show, respectively, at first order for \(\eta (z)\) and second order for r(z).

Fig. 1
figure 1

We plot \(t_0'(0)\) as a function of \(K_1\). This is the quantity determining the leading order effect for \(\frac{\delta H}{H}\) as shown in Eq. (51)

Full size image

In order to obtain a formula for the Hubble parameter as a function of the redshift we need to substitute Eqs. (17) and (18) in Eq. (19),

$$\begin{aligned} H(z)= & {} H(\eta (z),r(z)). \end{aligned}$$
(42)

After expanding up to second order in z we get

$$\begin{aligned} H(z)= & {} H_0+H_1 z + H_2 z^2,\end{aligned}$$
(43)
$$\begin{aligned} H_1= & {} \frac{1}{2} H_0 \Omega _M \left( \frac{4 t_0'(0)}{a_0}+3\right) ,\end{aligned}$$
(44)
$$\begin{aligned} H_2= & {} \frac{1}{72 a_0^2 \Omega _{\Lambda } \Omega _M} [a_0^2 H_0(20(\zeta _0-1)K_1^2\nonumber \\&+\,48 K_1 \Omega _{\Lambda } \Omega _M+27 \Omega _{\Lambda }(4-3 \Omega _M) \Omega _M^2)\nonumber \\&+\,6 a_0 H_0 \Omega _{\Lambda } \Omega _M t_0'(0)(20 K_1+9 (8-5 \Omega _M)\Omega _M)\nonumber \\&+\,18 \Omega _{\Lambda }\Omega _M^2(H_0 (25+\nonumber \\&-\,12 \Omega _M)t_0'(0)^2+5 t_0''(0))]. \end{aligned}$$
(45)

The procedure to reduce the analytical formula to this form is rather complicated since it involves the expression wherever possible of all the intermediate terms of physically meaningful quantities using the properties of the Weierstrass elliptic functions [41]; we give more details in Appendix A.

As we can see from the first order coefficient \(H_1\), the leading order in \(t_0'(0)\) determines the sign of the correction with respect to the homogeneous case, and for this reason we plot \(t_0'(0)\) as a function of \(K_1\) in Fig. 1. At second order we have a more complicated dependency for \(H_2\), which involves also \(K_2\) and \(t''_0(r)\).

We can easily interpret the linear behavior shown Fig. 1, applying the chain rule for the derivative

$$\begin{aligned} t_0'(0)= & {} \frac{\partial t_0(r) }{\partial k}\frac{\partial k}{\partial r}\Big |_{r=0}=\alpha K_1,\end{aligned}$$
(46)
$$\begin{aligned} \alpha= & {} (a_0 H_0)^{-3} \frac{\partial t_0(r) }{\partial k}\Big |_{k=k_0}\approx -0.57. \end{aligned}$$
(47)
Fig. 2
figure 2

The density contrast \(\delta =\frac{\delta \rho }{{\rho _b}}\) is plotted as a function of the redshift for three different inhomogeneities profiles modeled by LTB solutions

Full size image

4 Relative difference of H(z) with respect to the homogeneous case

For a flat FLRW solution the expansion rate is given by

$$\begin{aligned} H^{\mathrm{FLRW}}(z)= H_0 \sqrt{\Omega _M (1+z)^3 + \Omega _{\Lambda }}. \end{aligned}$$
(48)

Since we want to compare the inhomogeneous case with the flat FLRW case we define the relative difference as

$$\begin{aligned} \frac{\delta H (z)}{H} = \frac{H^{\Lambda LTB}(z)}{ H^{\mathrm{FLRW}}(z)} - 1, \end{aligned}$$
(49)

where we denote by \(H^{\Lambda LTB}(z)\) the expansion rate defined in Eq. (42).

We can now expand the above expression at low redshift to get

$$\begin{aligned} \frac{\delta H (z)}{H}= & {} \frac{\delta H_{1}}{H} z + \frac{\delta H_{2}}{H} z^2 +\cdots ,\end{aligned}$$
(50)
$$\begin{aligned} \frac{\delta H_{1}}{H}= & {} \frac{2 \Omega _M t_{0}'(0)}{a_0}=\frac{2\alpha \Omega _M}{a_0}K_1, \end{aligned}$$
(51)
$$\begin{aligned} \frac{\delta H_{2}}{H}= & {} \frac{1}{36 a_0^2 H_0 \Omega _{\Lambda } \Omega _M} [ 2 a_0^2 H_0 K_1 (5(\zeta _0-1) K_1\nonumber \\&+\,12\Omega _{\Lambda } \Omega _M) +3 a_0 H_0 \Omega _{\Lambda } \Omega _M t_{0}'(0) (20 K_1\nonumber \\&+\,9 (8-9 \Omega _M)\Omega _M)+9 \Omega _{\Lambda } \Omega _M^2\nonumber \\&\times \,(H_0 (25-12 \Omega _M) t_{0}'(0)^2+5 t_{0}''(0))]. \end{aligned}$$
(52)

It is straightforward to verify that in the homogeneous limit, when \(K_0=K_1=K_2=0\), as expected, \(\delta H_1 = \delta H_2 =0\), since \(t_0'(0) = t_0''(0) = 0\). From the expression for \(\frac{\delta H_{1}}{H}\) the crucial role played by \(t_0'(0)\) is clear, which determines the sign of the relative difference at leading order in the redshift, and according to Eqs. (46) and (47) it is proportional to \(K_1\).

Using Eq. (47) we can also see that for a fixed \(K_0\) the sign of \(K_1\) determines the sign of \(\frac{\delta H_{1}}{H}\) at very low redshift when the second order contributions can be neglected. This is in good agreement with Figs. 3 and 4, which show an approximate linear behavior with opposite slope, corresponding to different values of \(K_1\).

Fig. 3
figure 3

The relative difference with respect to the homogeneous case \(\frac{\delta H(z)}{H}\) is plotted as a function of the redshift for three different density contrasts. The colors correspond to the density contrasts in Fig. 2. The solid lines are for the numerical calculation, the dashed lines are for the analytical formulas, and the dot-dashed lines are for the perturbation theory result

Full size image
Fig. 4
figure 4

The relative percentual difference \(\Delta (z) = 100 \frac{\delta H/H-(\delta H/H)^{\mathrm{num}}}{(\delta H /H)^{\mathrm{num}}}\) of \(\delta H/H\) with respect to the numerical results is plotted as a function of the redshift for the analytical formula (dashed) and for perturbation theory (dot-dashed). The colors correspond to the density contrasts in Fig. 2. As can be seen the low-redshift formula is always better, especially for larger density contrasts

Full size image

4.1 Testing the accuracy of the formulas

In order to test the accuracy of the analytical formulas and compare it with perturbation theory we perform numerical calculations using LTB solutions corresponding to the density contrasts shown in Fig. 2. Using large density contrasts we can test the limit of the perturbation theory results and compare them with the low-redshift expansion.

As seen in Figs. 3 and 4 the analytical formula for the relative difference of H(z) with respect to the homogeneous case is in good agreement with \(\frac{\delta H(z)}{H}\) obtained by integrating numerically the geodesics and background equations and is more accurate than the perturbation theory. From Fig. 4 in particular we can see that the agreement with the exact numerical calculation for the redshift expansion is in general better than that of the perturbation theory result. For larger density contrasts, as expected, the perturbative calculation is increasingly inaccurate while the red-shit expansion is still in good agreement with the exact numerical result, since it does not rely on the assumption of a small density contrast. This is due to the fact that perturbation theory is based on the assumption that \(\delta \rho /\rho \ll 1\), while the low-redshift expansion is based on an exact solution of the Einstein equations.

5 Conclusions

We have derived a low-redshift analytical formula for the Hubble parameter for an observer at the center of a spherically symmetric matter distribution, using an exact solution of Einstein’s field equations. Such a formula is in good agreement with numerical calculations and is more accurate than the perturbation theory result, especially for large density contrasts. This is due to the fact that perturbation theory is based on the assumption that \(\delta \rho /\rho \ll 1\), while the low-redshift expansion is based on an exact solution of the Einstein equations.

If the H(z) observations will show deviations from the \(\Lambda CDM\) predictions compatible with the formulas we have derived, this could be considered an independent evidence of the existence of a local inhomogeneity, and the formulas could be used to determine the characteristics of this local structure.

While the expansion for r(z) depends on our choice of radial coordinate, the formulas for the H(z) are independent of it, since both H and the redshift are physically observable quantities and as such are independent of the gauge choice. Since in the derivation of the formulas the inhomogeneity profile is determined by the coefficients of the expansion of the function K(r), there is still some coordinate dependency left in the way we parameterize the inhomogeneity. While our choice of gauge is quite natural, since in the FRW gauge, corresponding to \(\rho _0(r)=\hbox {const}\), the radial coordinates reduces to the radial FRW comoving coordinate in the limit in which the function k(r) goes to zero, a totally coordinate independent formula would still be preferable. For this reason in the future it will be interesting to derive formulas which are completely independent of the choice of the coordinate system, parameterizing the inhomogeneity in terms of the redshift expansion of the density \(\rho (z)\), which is a physical observable, instead of using the expansion of K(r).