Spherically-symmetric solutions in general relativity using a tetrad-based approach

We present a tetrad-based method for solving the Einstein field equations for spherically-symmetric systems and compare it with the widely-used Lemaître– Tolman–Bondi (LTB) model. In particular, we focus on the issues of gauge ambiguity and the use of comoving versus ‘physical’ coordinate systems. We also clarify the correspondences between the two approaches, and illustrate their differences by applying them to the classic examples of the Schwarzschild and Friedmann–Lemaître– Robertson–Walker spacetimes. We demonstrate that the tetrad-based method does not suffer from the gauge freedoms inherent to the LTB model, naturally accommodates non-uniform pressure and has a more transparent physical interpretation. We further apply our tetrad-based method to a generalised form of ‘Swiss cheese’ model, which consists of an interior spherical region surrounded by a spherical shell of vacuum that is embedded in an exterior background universe. In general, we allow the fluid in the interior and exterior regions to support pressure, and do not demand that the interior region be compensated. We pay particular attention to the form of the solution in the intervening vacuum region and illustrate the validity of Birkhoff’s theorem at both the metric and tetrad level. We then reconsider critically the original theoretical arguments underlying the so-called Rh = ct cosmological model, which has recently B Do Young Kim dyk25@mrao.cam.ac.uk Anthony N. Lasenby a.n.lasenby@mrao.cam.ac.uk Michael P. Hobson mph@mrao.cam.ac.uk 1 Astrophysics Group, Cavendish Laboratory, JJ Thomson Avenue, Cambridge CB3 0HE, UK 2 Kavli Institute for Cosmology, Madingley Road, Cambridge CB3 0HA, UK


Introduction
Spherically-symmetric solutions in general relativity are of fundamental importance to the study of compact objects, black holes and cosmology. Indeed, two of the oldest and most commonly studied exact solutions of Einstein's field equations are spherically symmetric: the Schwarzschild metric [1] describes the gravitational field outside a static spherical massive body, and the Friedmann-Lemaître-Robertson-Walker (FLRW) metric [2][3][4][5][6][7][8] describes a homogeneous and isotropic universe in terms of the evolution of its scale factor with cosmic time. Moreover, it was not long before McVittie [9,10] combined the Schwarzschild and FLRW metrics to produce a new spherically-symmetric solution that describes a point mass embedded in an expanding universe, although there still remains some debate regarding its physical interpretation [11,12].
Subsequently, there have been numerous studies of the general-relativistic dynamics of self-gravitating spherical systems. For example, Misner et al. [13] describe the spherically-symmetric collapse of a 'ball of dust' having uniform density and zero pressure that is embedded in a static vacuum exterior spacetime, and later generalise their results to incorporate pressure internal to the object. By contrast, 'Swiss cheese' models [14] consider an exterior expanding FLRW universe, albeit pressureless, in which a pressureless spherical object is embedded and surrounded by a 'compensating void' that itself expands into the background and ensures that there is no net gravitational effect on the exterior universe.
A more realistic description than the Swiss cheese models is provided by models based on the Lemaître-Tolman-Bondi (LTB) solution [15][16][17]. Such models can incorporate an arbitrary (usually continuous) density profile for the central object, which is usually not compensated but can be made so by an appropriate choice of initial radial density and velocity profiles. Nonetheless, these models again assume both the interior and exterior regions to be pressureless, although it is possible to accommodate cosmological models with uniform pressure [18,19]. The standard LTB metric can be extended to include non-uniform pressure but only when it is anisotropic [20][21][22]. Lastly, a more generalised version of the LTB solutions has been presented in [23] to describe a central object with pressure embedded in a static vacuum exterior.
A recent resurgence of interest in Swiss cheese and LTB models has been prompted by the possibility that they may provide an explanation for observations of the acceleration of the universal expansion, without invoking dark energy. This might occur if we, as observers, reside in a part of the universe that happens to be expanding faster than the region exterior to it. By observing a source in the exterior region, one would then measure an apparent acceleration of the universe's expansion, but this would be only a local effect. The effects of local inhomogeneities on the apparent acceler-ation of the universe have been widely studied [24][25][26][27][28][29][30][31], and have been linked with the observations of distant Type-Ia supernova. However it has been shown that such models would induce variations in the CMB black-body spectrum through scattering, which disagrees with observations [32,33]. In addition, LTB models have been used to study the effects of inhomogeneities on observed cosmological parameters, such as the Hubble constant [34][35][36], and to calculate effects of a void as a possible explanation for the cold spot in the cosmic microwave background (CMB) [37,38].
The LTB model does, however, have some limitations. In addition to the usual restriction to pressureless systems, the LTB model is typically expressed in comoving coordinates and thus provides a Lagrangian picture of the fluid evolution that can be difficult to interpret. More importantly, the LTB metric contains a residual gauge freedom that necessitates the imposition of arbitrary initial conditions to determine the system evolution.
As a consequence, we have for some time adopted a different, tetrad-based method for solving the Einstein field equations for spherically-symmetric systems. The method was originally presented in [39] in the language of geometric algebra, and was recently re-expressed in more traditional tetrad notation in [40,41]. The advantages of the approach are that it can straightforwardly accommodate non-uniform pressure, has no gauge ambiguities (except in vacuum regions, as we shall discuss later) and is expressed in terms of a 'physical' (non-comoving) radial coordinate. As a result, the method has a clear and intuitive physical interpretation. Indeed, the gauge choices employed result in equations that are essentially Newtonian in form.
In [42,43], we applied the method to modelling the evolution of a finite-size, spherically-symmetric object with continuous radial density and velocity profiles that is embedded in an expanding background universe (either spatially-flat, open or closed) and compensated so that it does not exert any gravitational influence on the exterior universe; the fluid was assumed to be pressureless throughout. In [40], we used the method to obtain solutions describing a point mass residing in either a spatially-flat, open or closed expanding universe containing a cosmological fluid with pressure. In the spatially-flat case, a simple coordinate transformation relates our solution to the corresponding one derived by McVittie, but for spatially-curved cosmologies our metrics are physically distinct from the corresponding McVittie metrics, as shown in Section 3.2 in [40]. Hence, we believe that the latter in fact do not necessarily describe spatially-curved cosmologies with a point mass in the centre, even though they may be solutions of Einstein's equations. In [41], we extended this study by applying the tetrad-based approach to obtain the solution describing the evolution of a finite spherical region of uniform interior density that is embedded in a background of uniform exterior density, where the fluid in both regions can support pressure and the expansion (or contraction) rates of the two regions are expressed in terms of interior and exterior Hubble parameters that are, in general, independent. We also derived a generalised form of the Oppenheimer-Volkov equation, valid for general time-dependent, spherically-symmetric systems.
In this paper, we present a comparison of our tetrad-based methodology with the LTB model for solving the Einstein field equations for spherically-symmetric systems. In particular, we focus on the issues of gauge ambiguity and the use of comoving versus 'physical' coordinate systems. We also clarify the correspondences, where they exist, between the two approaches. In addition, we extend the analysis presented in [41] by applying our tetrad-based method to a generalised form of 'Swiss cheese' model, which consists of an interior spherical region surrounded by a spherical shell of vacuum that is embedded in an exterior background universe. In general, we allow the fluid in the interior and exterior regions to support pressure, and we demand neither that the interior region be compensated, nor that the interior and exterior regions be uniform. Nonetheless, our principal focus is the case in which the fluid in the interior and exterior regions has uniform (although, in general, different) densities. In particular, we pay special attention to the form of the solution in the intervening vacuum region and verify the validity of Birkhoff's theorem, the usual interpretation of which has recently been brought in question [44]. This investigation allows us to reconsider critically the original theoretical arguments underlying the so-called R h = ct cosmological model [45], which has recently received considerable attention. These considerations in turn elucidate the behaviour of a number of 'horizons' during the general-relativistic evolution of a spherically-symmetric self-gravitating matter distribution, which does not appear to have been widely discussed in the literature.
The structure of this paper is as follows. In Sect. 2, we outline our tetrad-based approach to solving the Einstein equations. In Sect. 3 we compare our approach to the more commonly-used LTB model. Note that the comparison serves as a pedagogical outline which presents the physical motivations behind each approach, with the aim of highlighting the advantages of using the tetrad-based method and 'physical' coordinates, by applying both methods to the familiar FLRW and Schwarzschild spacetimes. We then present original results, where we apply our tetrad-based method to describe the evolution of a generalised form of 'Swiss cheese' model in Sect. 4 and investigate the validity of Birkhoff's theorem in its vacuum region in Sect. 5. We discuss the R h = ct cosmology in Sect. 6 and describe the generic evolution of a number of cosmological 'horizons' in Sect. 7. Finally, we present our conclusions in Sect. 8. We adopt natural units c = G = 1 throughout.

Tetrad-based solution for spherical systems
In a Riemannian spacetime in which events are labelled with a set of coordinates x μ , each point has the corresponding coordinate basis vectors e μ , related to the metric via e μ · e ν = g μν . At each point we may also define a local Lorentz frame by another set of orthogonal basis vectorsê a (Roman indices), which are not derived from any coordinate system and are related to the Minkowski metric η ab = diag(1, −1, −1, −1) viaê a ·ê b = η ab . The relationship between the two sets of basis vectors is defined in terms of tetrads, or vierbeins e a μ , where the inverse is denoted e a μ : It is not difficult to show that the metric elements are given in terms of the tetrads by g μν = η ab e a μ e b ν . The local Lorentz frames at each point define a family of ideal observers whose worldlines are the integral curves of the timelike unit vector fieldê 0 . Along a given worldline, the three spacelike unit vector fieldsê i (i = 1, 2, 3) specify the spatial triad carried by the corresponding observer.
For our spherically-symmetric system, we work in terms of the tetrad components f 1 = e 0 0 , g 1 = e 1 1 and g 2 = e 0 1 , as described in [40], which define the system via the line-element where dΩ is an element of solid angle and we have adopted a 'physical' (noncomoving) coordinate 1 for which the proper area of a sphere of radius r is 4πr 2 . Finally, the remaining gauge freedom (which leaves the line-element unchanged) can be employed (at least in non-vacuum regions) to choose the timelike unit frame vectorê 0 at each point to coincide with the four-velocity of the fluid at that point. Thus, by construction, the four-velocity v of a fluid particle (or an observer comoving with the fluid) has components [v a ] = [1, 0, 0, 0] in the tetrad frame. Since v μ = e μ av a , the four-velocity may be written in terms of the tetrad components and the coordinate basis vectors as v = f 1 e 0 + g 2 e 1 . Thus, the components of a comoving observer's four-velocity in the coordinate basis are simply where dots denote differentiation with respect to the observer's proper time τ .
As a consequence of this final gauge choice, it is convenient to define the two linear differential operators We may identify L t as the derivative with respect to the proper time of a comoving observer, since L t =ṫ∂ t +ṙ ∂ r = d/dτ , and similarly one may show that L r coincides with the derivative with respect to the radial proper distance of a comoving observer. Moreover, since g 2 is the rate of change of the r coordinate of a fluid particle with respect to its proper time, it can be physically interpreted as the fluid 3-velocity. We will therefore, in general, use g 2 and v interchangeably in our analysis. It is also convenient to introduce explicitly the spin-connection coefficients F ≡ ω 0 11 and G ≡ ω 1 00 , as described in [40], which are both, in general, functions of t and r . Since we are assuming standard general relativity, however, for which torsion vanishes, the spin-connection coefficients can be written entirely in terms of the tetrad components and their derivatives. For the torsion to vanish and for the resulting Riemann tensor to satisfy its Bianchi identity, the spin-connection coefficients F and G and the non-zero tetrad components f 1 , g 1 and g 2 must satisfy the relationships where the explicit solution for f 1 contains no arbitrary function of t, because one can always be absorbed by a further t-dependent rescaling of the time coordinate (which does not change f 2 ). For matter in the form of a perfect fluid with proper density ρ and isotropic restframe pressure p, the Einstein field equations and the contracted Bianchi identities lead to the following system 2 of dynamical and continuity equations [39] L r p = −G(ρ + p), where we have defined the function of t and r (in general) and Λ is the cosmological constant and M is the Misner-Sharp mass. The physical interpretation of the functions F, G and M is straightforward. As shown in [40], for an object in general radial motion (not necessarily co-moving with the fluid) with four-velocity components [û a ] = [û 0 ,û 1 , 0, 0] in the tetrad frame, the corresponding components of the object's four-acceleration arê a 0 =u 0 + Gû 0û1 + F(û 1 ) 2 , and its proper acceleration is α = −â bâ b , which provides a physical interpretation of the functions F and G. In particular, for the special case in which the object is co-moving with the fluid, one has [û b ] = [1, 0, 0, 0] and so [â b ] = [0, G, 0, 0]. Thus the proper acceleration of a fluid particle is α = G in the radial direction. Indeed, the L r p-equation in (5) shows that, in the absence of a pressure gradient, G vanishes and so the motion becomes geodesic. The physical interpretation of the function M can be obtained from the forms of the equations in (5) in which it appears. In particular, the L r M-equation can be written simply as ∂ r M = 4πr 2 ρ, which shows that M plays 2 In [39], two further equations are given, namely L r g 1 = Fg 2 + M r 2 − 1 3 Λr − 4πrρ and L t g 2 = Gg 1 − M r 2 + 1 3 Λr − 4πr p, but these may be derived from the L r M and L t M equations, respectively, in combination with the definition of M given in (6). the role of an intrinsic mass that is determined by the amount of mass-energy in a sphere of radius r , also known as the Misner-Sharp mass [46,47]. It is useful to note that [39] have shown that in spherically symmetric systems, M and r appear explicitly in the eigenvalues of the Weyl tensor, and r is also a measurable quantity; hence the name 'physical' coordinate. As they are both intrinsic (i.e. measurable) quantities, it is useful to construct our equations in terms of these variables.
The Eqs. (4)-(6) thus have clear physical interpretations and contain no residual gauge freedom (in non-vacuum regions). In particular, given an equation of state p = p(ρ), and initial data in the form of the density ρ(r, t 0 ) and the velocity g 2 (r, t 0 ), the future evolution of the system is fully determined. This is because ρ determines p and M on a time slice and the definition of M then determines g 1 . The equations for L r g 2 , L r p and L r f 1 then determine the remaining information, namely F, G and f 1 respectively, on the time slice. Finally, the L t ρ equation and L t M equation (together with the definition of M) enable one to propagate ρ and g 2 , respectively, to the next time slice and the repeat the process. The equations can thus be implemented numerically as a simple set of first-order update equations. This approach was illustrated in [39,42,43].
An alternative way of solving the system of Eqs. (4)- (6), which was employed in [40,41,48], is not to impose an equation of state, but instead specify a form for ρ(r, t) for all t or, equivalently, a form for M(r, t) followed by use of the L r M. This method is used in Sect. 4. We merely note here that one may obtain the fluid pressure p(r, t) by first using the L t M equation to eliminate f 1 from the L r p equation, which then yields the 'generalised Oppenheimer-Volkov' equation [41] This equation is, in fact, valid for any spherically-symmetric perfect fluid system and reduces to the standard Oppenheimer-Volkov equation with a cosmological constant [49,50] for a static spherically-symmetric system. Finally, although the system of Eqs. (4)-(6) accommodates non-zero pressure, it is worth considering briefly the special case of a pressureless fluid. In this case, the L r p equation forces G to vanish, so the motion of the fluid particles becomes geodesic and the L r f 1 equation forces f 1 = 1. Consequently, the components in the coordinate basis of the four-velocity of a fluid particle are where dots denote differentiation with respect to the particle's proper time τ . Sincė t = 1, the coordinate time matches the proper time of all observers comoving with the fluid. Hence the Newtonian gauge is a synchronous one: a global 'Newtonian' time is recovered on which all comoving observers agree (provided all clocks are synchronised initially). 3 Furthermore, combining the L t M equation and the definition of M yields (∂ t +g 2 ∂ r )g 2 = −M/r 2 + 1 3 Λr , which has the form of the Euler equation in Newtonian fluid dynamics (recalling that g 2 is the fluid velocity v). Finally, setting Λ = 0 for a moment, the definition of M can itself be rearranged to give 1 2 g 2 − M/r = 1 2 (g 2 1 − 1), which is the Bernoulli equation for zero pressure and total (non-relativistic) energy per unit mass 1 2 (g 2 1 − 1) (i.e. after subtraction of the rest-mass energy).

Application to Schwarzschild spacetime
As an illustration of our approach, we now apply it to the special case in which the matter source is concentrated at the single point r = 0 and the cosmological constant vanishes (see also [39]). For such a solution, ρ = p = 0 everywhere away from the origin and so that M = constant. One can show that the Eqs. (4)-(6) reduce to an under-determined system of equations, such that additional gauge fixing is required to determine an explicit solution. This occurs because in the final part of our gaugefixing procedure described above, one chooses the timelike unit Lorentz frame vector at each point to coincide with the fluid four-velocity at that point, which clearly cannot be performed in a vacuum region. Nonetheless, one may instead choose the timelike unit frame vector to coincide with the four-velocity u of some radially-moving test particle (which need not necessarily be in free-fall), so that its components in the tetrad frame are [û a ] = [1, 0, 0, 0] and hence in the coordinate basis one has [u μ ] ≡ [ṫ,ṙ ,θ,φ] = [ f 1 , g 2 , 0, 0], as previously. This ensures that our previous physical interpretations of the tetrad and spin-connection components still hold. It remains, however, to choose a particular class of radially-moving test particle, and the simplest and most natural choice is a radially free-falling particle that was released from rest at r = ∞. From the definition (6) of M (with Λ ≡ 0), one sees that g 1 corresponds to the total energy per unit rest mass of an infalling particle, and so for a particle released from rest at infinity one should adopt the gauge condition g 1 = 1. It is then a simple matter to obtain expressions for the remaining tetrad components and spin-connection coefficients. The resulting non-zero tetrad components are and the spin-connection coefficients F and G read We note that the condition G = 0 is clearly consistent with the geodesic motion of the test particles. The line-element (2) corresponding to the tetrad components (9) is given by which we recognise as the Schwarzschild spacetime line-element expressed in terms of Painlevé-Gullstrand coordinates [51,52]. This coordinate system has a number of desirable features. For example, the line-element is regular for all positive values of r and the spacelike hypersurfaces t = constant have Euclidean geometry. Moreover, from (9), the non-zero components of the four-velocity of a particle released from rest at infinity are immediatelyṫ and so we recover an essentially Newtonian description of the motion. In particular, we see that t coincides with the proper time of such particles. Note that we can recover the standard form of the Schwarzschild line-element in Schwarzschild coordinates by choosing the preferred class of test particle to have fixed spatial coordinates, which leads to the tetrad components and the spin-connection coefficients We note that G = 0 is consistent with non-geodesic motion of the test particles.

Application to FLRW spacetime
As a second illustration of our approach, we apply it to the special case of a homogeneous and isotropic spacetime, as assumed in cosmology. This corresponds to setting ρ and p to be functions of t only, such that M(r, t) = 4 3 πr 3 ρ. One can show that the non-zero tetrad components are given by where k is an arbitrary constant of integration, and the spin-connection coefficients F and G are We note that the condition G = 0 demonstrates that the fluid particles move geodesically, since there are no pressure gradients. Moreover, one can show One also obtains the further (although clearly not independent) dynamical equation which we recognise as the standard 'acceleration' cosmological field equation expressed in terms of the Hubble parameter H (t). Indeed, from (15), the non-zero components of the four-velocity of a fluid particle are immediatelẏ which both verifies that t coincides with proper time of such particles and recovers Hubble's law. Thus, our approach has led us to work directly with H (t), which is an intrinsic and measurable quantity, rather than the more usual scale factor, which we will denote by S(t). Nonetheless, we can make contact with the latter simply by setting The line-element (2) corresponding to the tetrad components (15) then reads and the dynamical equations (18) and (19) become which we recognise as Friedmann's cosmological field equations in their standard form. Note that on defining the comoving radial coordinater ≡ r/S(t), the lineelement (21) takes the usual FLRW form.

Comparison with LTB model
In contrast to our tetrad-based approach, the LTB model [15][16][17] is based on the use of a comoving radial coordinate, which we denote byr , and the assumption of a diagonal form for the metric. It is also usual to choose the time coordinate, which we denote byt, to coincide with the proper time measured by observers comoving with the fluid, but for the moment we will consider the Lemaître metric [53], a slightly more general version of the LTB metric in which this requirement is not enforced. Thus, we consider a line element of the form where, in general, A, B and R may be arbitrary functions of bothr andt. Note that the LTB metric corresponds to setting A = 1. We may understand the relationship between the line-element (24) and that given in (2), obtained using our tetrad-based approach, by performing a coordinate transformation that expresses the latter in terms of a comoving radial coordinate and brings it into diagonal form. We therefore consider the coordinate transformation Note that, although the time coordinates coincide, we still label the new one ast, since the partial derivatives ∂/∂t and ∂/∂t will, in general, be different because they hold fixed r andr , respectively. One may verify thatr is a comoving radial coordinate by recalling that the four-velocity components of a comoving observer in the Newtonian . The physical nature of the transformation (25) may be further clarified by noting that where L t , defined in (3), is the derivative with respect to the proper time of a comoving observer; indeed this is consistent with our finding above thatṫ = f 1 . Since L t may be considered as a relativistic form of convective derivative, one may interpret the transformation (25) as moving from a Eulerian to a Lagrangian description of the fluid motion. Similarly, one finds that where L r , also defined in (3), is the derivative with respect to the proper radial distance of a comoving observer. Under the transformation (25) to a comoving radial coordinate, the line-element (2) takes the diagonal form One should first note that this has been achieved without having to specify ∂r/∂r; this demonstrates that the Lemaître (and hence LTB) metric (24) possesses a residual gauge freedom, in contrast to the line-element (2) (recall that the final gauge choice made in Sect. 2 leaves the form of (2) unchanged). Comparing (24) and (28), one first identifies that r = R(r ,t) and hence the three non-zero tetrad components used in Sect. 2 are given by where the final result is obtained using (25). The expressions for the spin-connection coefficients F and G are obtained from the relations (4) and read Substituting the expressions (29) and (30) into the dynamical and continuity equations (5), one then obtains where the expression for M in (6) now becomes The assumed form (24) of the line-element and the system of Eqs. (31)-(35) constitute a generalised form of the LTB model that can accommodate pressure and a non-zero cosmological constant. Nonetheless, unlike the tetrad-based approach, this model still possesses a gauge freedom, since ∂r R remains arbitrary. Thus, to determine the evolution of the system, one must first choose a form for the function R(r ,t * ) at some timet * (usually given as an initial condition), which can sometimes be difficult to interpret physically.
To make contact with the standard LTB model, we may now set A = 1 in the line-element (24), so thatt coincides with the proper time measured by observers comoving with the fluid. Hence, like the Newtonian gauge, the standard LTB model employs a synchronous coordinate system. In terms of the tetrad components (29) and spin-connection coefficients (30), setting A = 1 corresponds to setting f 1 = 1 and G = 0, and hence (26) shows that the operators ∂t and L t then coincide. Thus, from the L t g 1 equation in (4) one finds that g 1 = g 1 (r ) is a function only ofr . Using the expression for g 1 in (29) and adopting the standard notation used in LTB models, we therefore define the function E(r ) by which we may choose arbitrarily provided that E(r ) > − 1 2 . It is also immediately clear from (31) that setting A = 1 requires the pressure gradient to vanish. Thus, the standard LTB line-element can at best accommodate a fluid with a non-zero homogeneous pressure. It is usual, however, to assume simply that the fluid is pressureless, in which case (34) shows that M = M(r) is a function only ofr . Given our earlier interpretation of M in the tetrad-based approach, one may thus verify the usual interpretation of M(r ) in the LTB model as the Misner-Sharp mass, which is naturally time-independent in the absence of pressure.

Application to Schwarzschild spacetime
As an illustration of the LTB model, and in particular to compare it with the tetradbased approach, we now apply it to the same physical situation as we considered in Sect. 2.1, namely that of a matter source concentrated into a single point and a vanishing cosmological constant. As previously, for such a solution, ρ = p = 0 everywhere away from the point mass, and so (34) implies that M = constant. Once again, the remaining system of equations is under-determined, and so some gauge-fixing is required. First we must choose a form for the arbitrary function E(r ). Similar to the approach adopted in Sect. 2.1, we may base our choice on some class of radially-moving test particle and, once again, the most natural choice is a radially free-falling particle released from rest at infinity. From the definition (36), we see that the choice of g 1 = 1 in our tetrad-based approach is equivalent to setting E(r ) = 0, which corresponds to the particle having zero energy (after subtraction of its rest mass energy).
One can show that the LTB equations can be then integrated to obtain 2 where the "bang-time"t b (r ) may be an arbitrary function ofr , but is usually chosen such that R(r ,t b (r )) = 0, which in this case requirest b (r ) =r . Thus, after this additional gauge-fixing, which was not required in the tetrad based approach, one obtains which is the line-element for the Schwarzschild spacetime expressed in Lemaître coordinates [15]. It is interesting that, although the tetrad-based approach and the LTB model both employ synchronous time coordinates and are based on the trajectories of radially infalling particles released from rest at infinity, the two methods naturally lead to the very different line-elements (11) and (37). This occurs because of the use of a 'physical' radial coordinate in the former, whereas the latter employs a comoving radial coordinate, and also the requirement that the LTB line-element be diagonal. In the authors' opinion, the former line-element, expressed in Painlevé-Gullstrand coordinates, is the more easily interpreted physically.

Application to FLRW spacetime
As an second illustration of the LTB model, we now apply it to the special case of a homogeneous and isotropic spacetime, as considered in Sect. 2.2 using the tetradbased approach. As before, this corresponds to setting ρ to be a function oft only, but for the LTB model we are limited to considering only pressureless fluids and so p = 0. In contrast to the tetrad-based approach, one must begin by making a gauge choice for the form for M(r). This is most naturally achieved by introducing the scale factor Moreover, the remaining LTB equations (36) and (34) then yield the standard Friedmann equation and cosmological continuity equation, respectively, namely where we have defined the Hubble parameter Thus, we see that the LTB model has led us directly to working in terms of the scale factor, in contrast to the tetrad-based approach used in Sect. 2.2, which led naturally to the Hubble parameter, which is a directly measurable quantity. Moreover, considerable gauge-fixing was required in the LTB model to obtain a definite form for the solution, whereas this was unnecessary in the tetrad-based approach.

Generalised Swiss cheese model
We now apply our tetrad-based approach to a generalised form of the Swiss cheese model. In its classic form, the Swiss cheese model consists of an exterior expanding FLRW universe, albeit pressureless, in which a pressureless spherical object is embedded and surrounded by a 'compensating void' that itself expands into the background and ensures that there is no net gravitational effect on the exterior universe. Such models were employed in some of the earliest attempts to describe non-linear cosmological inhomogeneities [54][55][56][57][58], since they have the advantage that analytical calculations can be performed, and compensation ensures observations in the exterior region can be modelled unambiguously. They have also been used in more recent attempts to characterise effects of inhomogeneities on cosmological observations, such as luminosity distance and perceived dark energy [59][60][61][62][63][64][65][66][67][68][69][70]. Nonetheless, the matter and velocity distributions are clearly unrealistic.
As mentioned in the Introduction, more realistic models can be constructed by working with continuous density and velocity profiles, while still restricting to spherical symmetry and ignoring pressure. Previous work using LTB models [71,72] has usually ignored compensation, and this can lead to subtleties in modelling observations in the exterior region, since it is not described by a homogeneous FLRW cosmology, such as a net gravitational effect seen in the ISW effect [38,73]. The initial density and velocity profiles must be carefully chosen so that streamline crossing is avoided [42,43]. Otherwise, shock fronts form and one must include pressure to produce a realistic model. Compensation in tetrad formalism is discussed in [42,43] and in the LTB formalism in [74,75].
In this section, however, our primary focus is not the modelling of realistic cosmological inhomogeneities or the prediction of observational effects in the exterior region. Rather, we wish merely to extend the analysis presented in [41] (hereinafter NLH3) by applying our tetrad-based method to a generalised form of Swiss cheese model, in which we allow the fluid in the interior and exterior regions to support pressure, in general, and do not demand that the interior region be compensated. Aside from intellectual curiosity, the motivations for this study are two-fold: we first wish to verify that Birkhoff's theorem holds in the vacuum region, the usual interpretation of which has recently been brought into question for related systems [44]; and, second, we wish to consider the validity of the theoretical arguments that underpin the R h = ct cosmological model [45,76,77], which has recently received considerable attention [78][79][80]. These investigations are presented, respectively, in Sects. 5 and 6 below.
As discussed in Sect. 2, instead of imposing an equation of state, p = p(ρ), we solve the system of Eqs. (4)-(6) by specify a form for ρ(r, t) for all t or, equivalently, a form for M(r, t) followed by use of the L r M. In general, the remaining equations need to be solved as a set of coupled PDEs. Nonetheless, as shown in [40,41,48], if ρ(r, t) is piecewise uniform in r , then one may combine the L t ρ, L t M and L r M equations to obtain an ODE in r that may be solved to obtain an expression for the fluid velocity g 2 (r, t) and hence F(r, t), albeit with each containing a time-dependent 'constant' of integration, and the definition of M then determines g 1 (r, t).
One may obtain the fluid pressure p(r, t) by using the 'generalised Oppenheimer-Volkov' equation (8), which requires the imposition of a boundary condition on the pressure at some radius. One may then complete the solution either by obtaining f 1 (r, t) from the L t ρ equation and hence G(r, t) from any other equation that contains it, or by obtaining G(r, t) from the L r p equation and then f 1 (r, t) from the L r f 1 equation. Fig. 1 Generalised Swiss cheese model a spherical interior region of uniform density ρ i (t) and radius a(t) is surrounded by a vacuum region of radius b(t), which itself resides in an exterior region with uniform density ρ e (t). The rates of expansion of the interior, vacuum and exterior regions are characterized by the 'Hubble parameters' H i (t), H v (t) and H e (t), respectively. In general, the fluid can support pressure and the interior region need not be compensated

Model specification
The generalised Swiss cheese model is illustrated in Fig. 1 and consists of a spherical interior region of uniform density ρ i (t) and radius a(t) surrounded by a vacuum region of radius b(t), which itself resides in an exterior region with uniform density ρ e (t).
The rates of expansion of the interior, vacuum and exterior regions are characterized by the 'Hubble parameters' H i (t), H v (t) and H e (t), respectively (the definition of H v (t) is discussed below). These functions are free for us to choose, and together with M(r, t) and initial conditions, completely specify the evolution of the system. As we show later on, these 'Hubble parameters' are, in fact, equal to the covariant Hubble scalar in both the interior and exterior regions. In general, the interior region need not be compensated and the fluid in both the interior and exterior regions can support pressure.
From the figure, we may write down an expression for the total mass-energy M(r, t) contained within a sphere of physical radius r at time t. It is clear that where M 0 is a constant and 3 is the mass contained within b(t) at time t, in excess of that which would be present due to the exterior background alone. For a compensated interior region, one thus has m(t) ≡ 0. We note that the system considered in NLH3 corresponds to setting b(t) = a(t), so that there is no vacuum region.
As we will show below, in the general case where the fluid supports pressure, to determine the dynamical evolution of the system completely one must specify the internal and external Hubble parameters H i (t) and H e (t), together with the evolution a(t) and b(t) of the two boundaries (and the density ρ * ≡ ρ i (t * ) of the interior region at some reference time t = t * ). Typically, one should take H e (t) to correspond to some expanding exterior universe of interest, but H i (t), a(t) and b(t) can, in principle, have any form. 4 This follows both from the presence of the vacuum region and from allowing the relationship between the fluid pressure and density to be arbitrary, since then the interplay between pressure and gravity may allow expansion or contraction of the interior and vacuum regions at any rate. This freedom would disappear, however, if one imposed an equation of state on the fluid. In particular, in the special case of a pressureless fluid, a(t) is straightforwardly determined from H i (t) (and the radius a * ≡ a(t * ) at some reference time t = t * ).
We note that by leaving H i (t) and H e (t) free to choose, we are treating the system as being composed of a mathematical fluid, that is intended to mimic the kinematical evolution (if not the physics) of a combination of baryonic gas and dark matter, having a single effective density and a single effective pressure required for stability (for more details, see [48]). This results in an effective 'equation' of state which depends on both r and t. Choosing the Hubble parameters and specifying the mass-energy M(r, t) means that we can solve for p(r, t) using the generalised Oppenheimer-Volkov equation (8). Hence the effective 'equation of state' is then determined. Treating the fluid as a single mathematical fluid avoids the complication of calculating the non-linear evolution of multi-fluid systems, whereby one would separate the fluid in each region into its baryonic and dark matter components.

Boundary conditions
Any spatial surface at which the density is discontinuous, and which may in general be moving, will trace out a 3-dimensional (timelike) hypersurface Σ in spacetime on which the solution must satisfy the Israel junction conditions [81,82]. Ifn μ are the covariant components of the unit (spacelike) normal to Σ, pointing from the inside to the outside, then the Israel junction conditions require both of the induced metric h μν = g μν +n μnμ and the extrinsic curvature K αβ = h α μ h β ν ∇ μnν to agree on Σ. For the model illustrated in Fig. 1, two such hypersurfaces are defined by Σ(t, r ) ≡ r − x(t) = 0, where x(t) can equal either a(t) or b(t). As discussed in [41], the componentsn μ are given by where ∂ t x ≡ dx(t)/dt; one may readily verify thatn μn μ = −1, as required. One sees immediately that, for the induced metric to be agree across the boundary x(t), one requires all three non-zero tetrad components f 1 , g 1 and g 2 to be continuous there.
Recalling that g 2 is the rate of change of the r coordinate of a fluid particle with respect to its proper time, and can thus be considered as the fluid velocity, the physical interpretation of the continuity of g 2 is that matter does not cross the boundary x(t) in either direction. This is consistent with the boundary x(t) comoving with the fluid, so that the situation depicted in Fig. 1 does indeed hold at all times. It does not imply, however, that M(x(t), t) is constant, since this quantity denotes the total energy within x(t), which may change as x(t) evolves, as is clear from the L t M and L r M equations in (5). Moreover, since L t corresponds to the derivative with respect to the proper time of an observer comoving with the fluid, then one requires where L t is evaluated at the boundary x(t). Thus, on the hypersurface Σ, one has dx(t)/dt = g 2 / f 1 and the expression (42) simplifies to After a long but straightforward calculation, one then finds that the only non-zero components of the extrinsic curvature of Σ are 5 Since the first Israel junction condition requires f 1 , g 1 and g 2 all to be continuous at the boundary x(t), then the second junction condition requires only that, in addition, ∂ r f 1 is continuous there. The above junction conditions have consequences for the continuity of other variables of interest. In particular, from the L r f 1 equation in (4), one has G = −(g 1 / f 1 )∂ r f 1 , which must therefore also be continuous at the boundary. Moreover, the L r p equation in (5) and the continuity of g 1 and G imply that the pressure p is also continuous across the boundary, although its radial derivative, in general, has a step there.
Finally, we also adopt the boundary condition at large r that all physical quantities tend to those of the exterior cosmology. For spatially-flat and open universes, this corresponds to the limit r → ∞, whereas for a closed universe one must instead consider the limit x(t) r < S(t), where S(t) is the universal scale factor which also corresponds to the curvature scale for a closed universe. In each case, we require the line-element (2) to tend at large r to the corresponding FLRW line-element (21) with H (t) = H e (t). Thus, for large r , one requires

Interior and exterior regions
From (41), one sees that the forms for M in the interior and exterior regions are the same as the model considered in NLH3, albeit with a different definition of m(t). Moreover, the same boundary conditions (46) apply at large r . Therefore, many of the equations derived in NLH3 remain valid. Specifically, in the interior region, the non-zero tetrad components are again given by where, following NLH3, a prime denotes differentiation with respect to t. In order to evaluate the above expressions for f 1,i and g 1,i , one requires forms for ρ i (t) and p i . Substituting the above expression for the fluid velocity g 2,i and the enclosed mass M = 4π 3 ρ i (t)r 3 from (41) into the generalised Oppenheimer-Volkov equation (8), and integrating, will yield an (integral) expression for p i in terms of H i (t) and ρ i (t), after imposing the condition that the pressure is continuous across the boundary a(t) and hence vanishes there. Thus, it only remains to determine ρ i (t), which is straightforwardly obtained for a given boundary evolution a(t) by recalling that M 0 ≡ 4π 3 ρ i (t)a 3 (t) is a constant. In the special case of a pressureless fluid, it is worth noting that one immediately has f 1,i = 1 and so (47) can be integrated to obtain ρ i (t), which then determines a(t).
For the exterior region, the non-zero tetrad elements are given by In order to evaluate the above expressions, one requires forms for p e , ρ e (t) and b(t).
Substituting the expressions for M and g 2,e into the generalised Oppenheimer-Volkov equation (8), integrating and imposing the condition that the pressure is continuous across the boundary b(t) and hence vanishes there, will yield an (integral) expression for p e in terms of ρ e (t), H e (t), b(t) and the (in general) time-dependent uniform fluid pressure p ∞ (t) at large r corresponding to the external cosmological model. One is free to specify H e (t), b(t) and p ∞ (t), and the function ρ e (t) may be determined from the following equations from NLH3 which remain valid in the exterior region: We recognise (53) and (54) as the standard cosmological fluid evolution equation and the Friedmann equation, respectively. Moreover, as in NLH3, these can be combined in the usual manner to yield the dynamical cosmological field equation which thus provides an expression for ρ e (t) in terms of H e (t) and p ∞ (t).
We are free to choose the boundary evolution b(t), and it is most convenient to do this by defining the 'vacuum Hubble parameter' H v (t), such that Equating the above expression with (52) evaluated on the boundary b(t), one obtains This then allows one to write (52) in the elegant form In a similar manner, one may write the expression (51) as where we have momentarily suppressed t-dependencies for brevity. It is worth noting that, for the special case of a pressureless fluid, one immediately has p e = 0 and

Hubble scalar and shear in the interior and exterior regions
The fluid velocity covariant derivative can be split into parts with specific symmetry properties [83,84]. The decomposition is given by (from equation (4.17) in [83], adapted to our metric signature) where ω μν is the vorticity tensor, σ μν is the shear tensor, h μν is the projection tensor into the 3D subspace orthogonal to the fluid velocity, and the volume expansion scalar, θ , is defined as Hence, the Hubble scalar is defined as The shear tensor is defined as the traceless component of the 'fully projected' part of the symmetric piece of ∇ ν v μ . Specifically we define such that the shear tensor is given by The final quantity of interest is the relativistic 4-acceleration vectoṙ which represents the degree to which matter moves under forces other than those of gravity (see for example, Section 2.1 in [84]). We now consider the values of these quantities in the interior and exterior regions of our spherically symmetric solutions.
Since the fluid 4-velocity components, in a (t, r, θ, φ) coordinate system, are given by v μ = dx μ dτ = ( f 1 , g 2 , 0, 0) (where g 2 is the fluid 3-velocity), the Hubble scalar (62) in terms of tetrad components is given by In fact, in both the interior and exterior regions, the Hubble scalar simply reduces to the interior and exterior Hubble parameters respectively, despite the 3-velocity field deviating from the Hubble flow. That this occurs in the interior region is unsurprising, since Eq. (49) shows that the fluid 3-velocity here coincides with a pure Hubble flow, r H i (t). In the exterior region, one again has H s (t) = H e (t) despite the fluid 3-velocity g 2 deviating from a pure Hubble flow, r H e (t) (see Eq. (58)). This occurs because the deviation is a function of time divided by r 2 , for which the contribution to the 4divergence ∇ μ v μ vanishes. We can construct the shear tensor using Eq. (64) above. If we denote S = 1 3 ∂g 2 ∂r − g 2 r , one can find using Eq. (64) that the eigenvalues of the shear tensor are 0 in the t direction, 2S in the r direction and −S in the θ and φ directions. We note this is traceless as required for a shear tensor.
This shows a nice complementarity to what we just found for the Hubble scalar, since here the part of g 2 corresponding to the Hubble flow gives zero shear, as expected, whilst for the exterior, the deviation − b 3 (t) r 2 (H e (t) − H v (t)), as given in Eq. (58), gives us the tangential shear eigenvalue −S. Hence the model studied here, at least in the exterior region, is different from those studied by [85], where the fluid is assumed to be wholly shear-free.
The fact a spherically symmetric fluid can have a Hubble scalar which is everywhere the same as the Hubble parameter at infinity (where there is just a cosmological flow), but nevertheless has non-zero shear (tending to zero at infinity), is perhaps worthy of further comment. Specifically, such a fluid would have to obey the irrotational version of one of the constraint equations which can be derived from the Ricci identity as applied to the velocity 4-vector. To discuss this briefly, we will use the notation employed in [86], which generalises the 1+3 covariant approach to fluid dynamics of [84] to the context of fluids with intrinsic spin, but is convenient here since it uses the same (+, −, −, −) metric as the present work.
The Ricci identity from [86] is and by taking an antisymmetric trace-free part, employing the Einstein equations and (for current purposes) setting the spin and vorticity to zero, we can derive the equation where D λ is the fully projected covariant derivative, which in this case would satisfy Evaluating this constraint equation for the shear tensor and volume expansion scalar found above, we find that there is a non-zero entry for μ = r , given (as a general function of the tetrad components) by This looks as though there is a problem, but in fact the RHS here is proportional to the Einstein tensor entry G tr , which due to our choice of fluid velocity vector, has to vanish (this can be confirmed independently via the explicit solutions for these quantities we have given in each region, including the external one). Thus our results for the shear and expansion factor, and in particular the fact that the expansion parameter can correspond to what one would have for a uniform Hubble flow, whilst the fluid still has non-zero shear, are fully consistent with the Ricci identity constraint on the shear divergence. Finally, we note that the 4-acceleration a μ = v ν ∇ ν v μ simply evaluates to a radial vector with magnitude G, as already discussed in Sect. 2.

Vacuum region
In the vacuum region, one has ρ = p = 0 and so L r M = 0 and L t M = 0, which together imply M is a constant, which we have denoted by M 0 . As we found in our discussion of the Schwarzschild spacetime in Sect. 2.1, in a vacuum region the system of Eqs. (4)-(6) reduces just to the relationships (4) between the tetrad and spin-connection components and the definition of M(= M 0 ) in (6), from which one finds that the quantity is a function of r only. No further equations yield new information, and one thus has an under-determined system of equations that requires additional gauge-fixing to determine an explicit solution. This occurs because in a vacuum region one clearly cannot choose the timelike unit Lorentz frame vector at each point to coincide with the fluid four-velocity at that point. As we did for the Schwarzschild spacetime, however, one may instead choose the timelike unit frame vector to coincide with the four-velocity u of some radiallymoving test particle (which need not necessarily be in free-fall), so that its components in the tetrad frame are [û a ] = [1, 0, 0, 0] and hence in the coordinate basis one has , as previously. Moreover, this ensures that our previous physical interpretations of the tetrad and spin-connection components still hold.
Unlike in the Schwarzschild spacetime, however, there is no simplest or most natural choice for the class of radially-moving test particle to use. All one requires is that the boundary conditions discussed in Sect. 4.2 hold at each boundary a(t) and b(t). It is most convenient to begin by choosing g 2,v = v(r, t), where v may be an arbitrary function satisfying the boundary conditions v(a(t), t) = g 2,i (a(t), t) = a(t)H i (t) and v(b(t), t) = g 2,e (b(t), t) = b(t)H v (t). Then g 1,v is easily found from (71) and is also continuous at both boundaries. Finally, eliminating G between the L r f 1 and L t g 1 equations in (4) gives the general relation and using (71) to substitute for g 1 , one then finds which may be straightforwardly solved for f 1,v . Gathering these results together, in the vacuum region one thus has Let is first consider the general case in which the fluid in the interior and exterior regions supports pressure. Suppose one specifies the form for v(r, t * ) at some time t * . One can then see from (73) that one also requires the profile f 1,v (r, t * ) in order to evolve v in time. Thus, both v(r, t * ) and f 1,v (r, t * ) need to be specified to determine the system. One should note, however, that there is no equation determining the time evolution of f 1,v ; hence f 1,v is free to take any value on any time slice, provided it satisfies the boundary conditions that both f 1 and ∂ r f 1 are continuous at each boundary, as shown in Sect. 4.2. Then the time evolution of v is determined.
The situation is somewhat simpler for the case in which the fluid in the interior and exterior regions is pressureless, since one may take f 1 = 1 everywhere and at all times. Hence, if one specifies the form for v(r, t * ) at some time t * , one can use (73) to evolve v in time [42]. In this case, (73) becomes which may also be derived directly by substituting the definition of M in (6) into L t M = 0. In either case, with or without pressure, the exact choice of the initial profile v(r, t * ) or f 1,v in the vacuum region has no physical effects on, for instance, the total redshift of a photon passing through the inhomogeneity [42], showing that the ambiguity in these functions is a gauge freedom and hence has no physical consequences.

Birkhoff's theorem
One of our motivations for considering the generalised Swiss cheese model is to illustrate that Birkhoff's theorem holds in the vacuum region, despite having timeevolving interior and exterior regions. Although its usual interpretation has recently been brought into question for physical systems of this type [44], one can demonstrate that Birkhoff's theorem does indeed hold in the vacuum region of our generalised Swiss cheese model using a traditional wholly metric-based approach, and that it also holds directly at the level of the tetrad components by transforming into Schwarzschild-de-Sitter form [87]. In particular, we note this does not depend on the radial distribution of matter in the interior or exterior regions or its state of motion, provided spherical symmetry holds.

Comparison with previous work
The validity of Birkhoff's theorem, or at least its usual interpretation, has recently been brought into question by [44]. To be clear, Birkhoff's theorem states that there exist coordinates for which the metric in a vacuum region surrounding any sphericallysymmetric matter distribution takes the standard Schwarzschild (de-Sitter) form with parameter M 0 equal to the enclosed interior mass, even when the vacuum region is itself embedded in an exterior spherically-symmetric matter distribution. Indeed, we have verified this. An immediate corollary of Birkhoff's theorem is that, if the enclosed mass M 0 = 0, there exist coordinates for which the metric in the vacuum region takes the standard Minkowski (de-Sitter) form.
As pointed out by [44], however, a common misinterpretation of Birkhoff's theorem is that the gravitational field anywhere inside a spherically-symmetric matter distribution is determined only by the enclosed mass. While this is true in Newtonian gravity, it does not hold in general relativity. This point is illustrated in [44] by considering the metric corresponding to a static thin spherical shell of mass m s at coordinate radius r = r s , surrounding a spherical central object of mass m i centered on the origin.
A further example of the gravitational field at some radius in a sphericallysymmetric matter distribution being determined by material external to that radius is provided by the system analysed in [41]. As mentioned previously, this system corresponds to setting b(t) = a(t) in the generalised Swiss cheese model discussed in Sect. 4, so that there is no vacuum region. In [41], this system is analysed separately using Newtonian gravity and general relativity. The former case, the Newtonian gravitational potential in the interior region is found to be independent of the properties of the exterior region, whereas the general-relativistic calculation shows that some of the tetrad components, and hence metric elements, in the interior do depend on the properties of the exterior region, such as its density ρ e (t) and Hubble parameter H e (t).

R h = ct cosmology
The current standard model of cosmology is based on the ΛCDM model, which provides a good fit to a wide range of cosmological observations. As pointed out by [88], however, for the best-fit ΛCDM model, the present-day Hubble distance is broadly consistent with ct 0 to within observational uncertainties, where t 0 is the current cosmic epoch. This suggests that the universe has expanded by an amount similar to what would have occurred had the expansion rate been constant. In the ΛCDM model, this correspondence is a peculiar coincidence, particularly since this situation should occur only once in the history of the universe; that we observe it to hold at the present epoch is thus intriguing.
It was therefore proposed by [45,76,77] that this correspondence is not coincidental, but should be satisfied at all cosmic times t. The resulting cosmological model, known as the 'R h = ct' model, has received considerable attention over the last few years, since it has been claimed to be favoured over the standard ΛCDM (and its variant wCDM with w = −1) by most observational data [89][90][91][92][93][94]. These claims have, however, recently been brought into question [95,96]. In addition, the comoving Hubble radius is constant in time in the R h = ct model, which cannot account for acoustic peak structure observed in the power spectrum of CMB anisotropies or for baryon acoustic oscillations.
In addition to observational objections, the validity of the physical argument underlying the R h = ct model has also been criticised by a number of authors [78][79][80]97]. Additionally, analytical claims for the model presented in [98] have also been demonstrated to be false [99].
Melia's main proposal begins with the introduction of the 'gravitational radius' R h (t), which is defined by the requirement that where, as previously, we adopt natural units c = G = 1 and M(r, t) = 4 3 πρ(t)r 3 is the mass-energy contained with a sphere of physical radius r at time t. Indeed, substituting this form for M(r, t) into (76), one immediately finds that 1/R 2 h (t) = 8 3 πρ(t). The key point in Melia's argument is the assertion that, in order to satisfy Weyl's postulate, the comoving gravitational radiusr h (t) ≡ R H (t)/S(t) should be independent of t, where S(t) is the scale factor defined in Sect. 2.2. In other words, each fluid particle has a fixed value ofr h .
Then, following Melia's assumption that k = 0 = Λ, the standard Friedmann equation (18) immediately allows one to make the identification R h (t) = 1/H (t), which is the Hubble radius. Since H (t) = ∂ t S(t)/S(t), the comoving gravitational (Hubble) radius in this case is given bŷ The requirement thatr h (t) is constant therefore implies S(t) ∝ t and so R h (t) = 1/H (t) = t or, on momentarily abandoning natural units, one obtains the eponymous R h (t) = ct. It is worth noting that substituting H (t) = 1/t into the dynamical cosmological field equation (19) yields Melia's so-called 'zero active mass' condition ρ + 3 p = 0, which is equivalent to demanding that the overall cosmological fluid has the equation-of-state parameter w = − 1 3 at all epochs. The argument given in [45] for supposing thatr H (t) = constant is that, according to Weyl's postulate, any proper distance in an FLRW spacetime must be the product of the scale factor S(t) and some fixed co-moving radial coordinate, and that the definition of R h (t) as a gravitational radius in (76) implies that it must be a proper distance. No real justification is given for this latter assertion.
In the next section, we show that the requirementr h (t) = constant does not follow from the definition (76). In doing so, our investigations below also elucidate the behaviour of a number of 'horizons' during the general-relativistic evolution of a spherically-symmetric self-gravitating matter distribution; although straightforward and interesting, this behaviour does not appear to have been widely discussed in the literature.

Evolution of horizons
In a homogeneous and isotropic cosmological model, let us consider an imaginary spherical boundary of radius a(t) that is comoving with the fluid and centred on some arbitrary origin. From the discussion in Sect. 4, the equation of motion for a(t) is L t a(t) =ȧ = H (t)a(t) (since f 1 = 1 in this case), where H (t) is the Hubble parameter characterising the evolution of the fluid. If M(t) denotes the mass-energy contained within this sphere, we define the corresponding Schwarzschild radius R S (t) as the solution of which clearly reduces to R S (t) ≡ 2M(t) if Λ = 0. We also define the Hubble radius R H (t) ≡ 1/H (t); note that this coincides with Melia's gravitational radius R h (t) in the case k = 0 = Λ, but differs from it in more general cosmological models (although such models were not considered by Melia). We now consider the behaviour of a(t), R S (t) and R H (t) for a selection of analytical spatially-flat (k = 0) expanding cosmological models, with cosmological constant Λ and fluid equation-of-state parameter w, whose evolution is determined by the Hubble parameter H (t). We note that, in the spatially-flat (k = 0) case, a(t) is the proper distance from the origin to the spherical boundary under consideration, and thatȧ is its rate of change with respect to the proper time of an observer comoving with the fluid. 6 The results are presented in Table 1. The values of w considered correspond to dust (w = 0), radiation (w = 1 3 ) and Melia's zero-active-mass condition (w = − 1 3 ). It is worth noting that it is only in the case w = 0 that the mass-energy M(t) contained within the comoving radius a(t) is constant; for other values of w the presence of non-zero pressure means that the fluid does work and hence M(t) changes with time.
For Λ = 0, the behaviour of the quantities listed in the table is shown in Fig. 2. For both the dust and radiation cases, one sees that a(t), R S (t) and R H (t) cross at a single point that we denote by t = t * , which coincides withȧ(t) dropping to 1 (or c in standard units). This behaviour is quite general and holds for any positive value of the parameter M 0 in Table 1. Thus, the comoving radius a(t) initially lies outside the Hubble radius. Indeed, one may consider thatȧ is allowed to be greater than c during this initial period due to the 'superhorizon' (or, more correctly, super-Hubble-radius) nature of a(t). Note that from the Friedmann equation, we know that 6 We further note that the proper distance is measurable in principle (albeit in a highly impractical manner!), by arranging for a 'cosmic conspiracy' in which comoving observers (who may identify themselves as such by the absence of a CMB dipole) along a radial line from r = 0 to r = a lay down rulers at the same instant of their proper time t (which may be determined, for example, by requiring the CMB temperature they measure to have a particular value). Moreover, by arranging for such a conspiracy to be repeated at a later time t + dt, the 'velocity'ȧ can be determined. Note, however, that when a is bigger than the current particle horizon, it would not be possible to arrange for such a conspiracy.

Table 1
Evolution of the comoving radius a(t), its rate of changeȧ(t), the mass-energy M(t) enclosed within a(t), the Schwarzschild radius (or radii) R S (t) defined in (78) and the Hubble radius R H (t) ≡ 1/H (t) for a selection of analytical spatially-flat (k = 0) expanding cosmological models, with cosmological constant Λ and fluid equation-of-state parameter w, whose evolution is determined by the Hubble parameter H (t) The values of w considered correspond to dust where we consider the general case with Λ = 0. Multiplying the above by a 2 gives uṡ Hence we see that whenȧ = 1, a = R H and also a = R S . This is precisely what we see in Fig. 2. a(t) enters the Hubble radius R H (t) at precisely the same moment as it exits the Schwarzschild radius R S (t); it is allowed to do the latter, since the fluid at the boundary is moving at speed c at this instant. Thus, one has two 'horizon crossings' taking place simultaneously and in opposite directions, which is not usually pointed out in the literature. Moreover, one sees that there is no reason for the Hubble radius to be comoving, which contradicts the central assumption of the R h = ct model. It is also worth noting that in the case of dust, for which the pressure vanishes everywhere, one can consistently 'cut-off' the fluid at the boundary a(t), and thereby consider a finite fluid ball surrounded by vacuum; the results for this case are precisely those given above. In fact, we can consider the case with a finite, pressureless ball in order to verify that the physics is correct. Consider a sphere of radius a(t) filled with dust, embedded in a vacuum. Its evolution must still follow that described in Fig. 2, such that it expands with time and eventually crosses R S . Hence it is the time reversed equivalent of a collapsing ball of dust. In the latter case, the physical event horizon exists only once the boundary of the sphere crosses the point r = 2M(t) (for the Λ = 0 case). Therefore, we can view the expanding sphere as a white hole, which has a horizon at R S , where objects are only allowed to move out initially. Once the boundary of the sphere crosses the horizon, the physical characteristics of the horizon cease to exist, such that it becomes physically possible for objects to fall in through the horizon. To verify that the crossing of the fluid boundary over R S is possible, we can calculate the speed at which a hypothetical observer on R S would move as measured by a comoving observer on the boundary of the sphere, as it crosses R S . One can show that the observed speed is given by dt , which is valid for any Λ and equation of state. As discussed above, when the crossing occurs, the fluid moves at speed c; hence the comoving observer sees R S move at the speed of light when he/she crosses it, for the dust case. This is not problematic, as it is impossible to have an observer hover at R S up until the point in time at which the ball crosses the horizon. In addition, the horizon vanishes when the edge of the fluid crosses it; hence the fluid is able to cross the horizon.
Note that when we have non-zero pressure, we can no longer 'cut-off' the fluid at the boundary a(t); hence we cannot consider a case where the various lines have physical characteristics.
The bottom panel in Fig. 2 shows the behaviour for the case w = − 1 3 , which corresponds to the zero active mass condition required by the R h = ct model. Clearly, in this case, the radii a(t), R S (t) and R H (t) all depend linearly on t. If one chooses M 0 = 1 8 , then one obtains the special case in which a(t) = R S (t) = R H (t) anḋ a(t) = 1 at all times. When M 0 > 1 8 , one hasȧ(t) > 1 and a(t) lies inside the Schwarzschild radius and outside the Hubble radius at all times (this is illustrated in Fig. 2, for which M 0 = 1 2 ). Conversely, if M 0 < 1 8 , one hasȧ(t) < 1 and a(t) lies outside the Schwarzschild radius and inside the Hubble radius for all t.
The behaviour of the quantities listed in Table 1 for Λ = 0 are illustrated in Fig. 3. In both cases, the values of M 0 and Λ have been chosen so that, at least for some values of t, the condition 1 − 9M 2 (t) √ Λ > 0 is satisfied and hence there exist two positive solutions for R S (t), which correspond to the Schwarzschild radius and the de Sitter radius, respectively.
In the case of dust (top panel), the mass-energy M(t) enclosed within the spherical boundary a(t) is constant, as expected, and thus so too are the two solutions for R S (t). One again sees that each solution for R S (t) intersects with a(t) and R H (t) at a single point, and that at both these intersections one hasȧ(t) = 1, as we expect from equation (80). Thus, as in the case Λ = 0, the comoving radius a(t) initially lies outside the Hubble radius and enters it at precisely the same moment as it exits the Schwarzschild radius, at which point the fluid at the boundary a(t) is moving at speed c. In the presence of non-zero Λ, however, the boundary a(t) later exits the Hubble radius again, at precisely the same moment that it enters the de Sitter radius, at which point the fluid at the boundary is again moving with speed c. It is again worth noting in this dust case that the absence of pressure allows one consistently to 'cut-off' the fluid at the boundary a(t), and thereby consider a finite fluid ball surrounded by vacuum with Λ = 0, and obtain identical results. In this case we again have a white hole, with a physical horizon at the Schwarzschild radius inside which objects are swept in, which ceases to exist at the moment in time the edge of the fluid crosses it. Lemaître-Tolman-Bondi (LTB) model. Although the LTB model is widely used, it has a number of limitations. In particular, in its usual form it is restricted to pressureless systems. Moreover, the LTB model is typically expressed in comoving coordinates and thus provides a Lagrangian picture of the fluid evolution that can be difficult to interpret. Perhaps most importantly, however, even in the absence of vacuum regions the LTB metric contains a residual gauge freedom that necessitates the imposition of arbitrary initial conditions to determine the system evolution. As a consequence, we have for some time adopted a different, tetrad-based method for solving the Einstein field equations for spherically-symmetric systems. The method was originally presented in [39] in the language of geometric algebra, and was recently re-expressed in the more traditional tetrad notation in [40,41]. The advantages of the tetrad-based approach are that it can straightforwardly accommodate pressure, has no gauge ambiguities (except in vacuum regions) and is expressed in terms of a 'physical' (non-comoving) radial coordinate. As a result, in contrast to the LTB model, the method has a clear and intuitive physical interpretation. Indeed, the gauge choices employed result in equations that are essentially Newtonian in form.
In comparing our tetrad-based methodology with the LTB model, we have focussed particularly on the issues of gauge ambiguity and the use of comoving versus 'physical' coordinate systems. We have also clarified the correspondences, where they exist, between the two approaches. As an illustration, we applied both methods to the classic examples of the Schwarzschild and Friedmann-Lemaître-Robertson-Walker (FLRW) spacetimes. In the former, we demonstrate that although the tetrad-based and LTB approaches both employ synchronous time coordinates and are based on trajectories of radially-infalling particles released from rest at infinity, the two methods lead to very different results corresponding to the use of Painlevé-Gullstrand and Lemaître coordinates, respectively. For the FLRW spacetime, we find that the LTB approach leads one to work directly in terms of the scale factor, whereas the tetrad-based method leads naturally to a description in terms of the Hubble parameter, which is a measurable quantity. Moreover, considerable gauge-fixing was required in the LTB model to obtain a definite solution, but this was unnecessary in the tetrad-based approach.
We have previously applied our tetrad-based method to modelling the evolution of a finite-size, compensated, spherically-symmetric object with continuous radial density and velocity profiles that is embedded in an expanding background universe, assuming zero pressure throughout [42,43]. We have also previously used the approach to obtain solutions describing a point mass residing in an expanding universe containing a cosmological fluid with pressure [40], and later a finite spherical region of uniform interior density embedded in a background of uniform exterior density, where the pressure may be non-zero in both regions [41]. To illustrate further the use of our tetrad-based approach, we here extended the analysis in [41] to a generalised form of 'Swiss cheese' model, which consists of an interior spherical region surrounded by a spherical shell of vacuum that is embedded in an exterior background universe. In general, we allow the fluid in the interior and exterior regions to support pressure, and we do not demand that the interior region be compensated. We find that much of the analysis in [41], including the specification of boundary conditions, can be applied with little modification, but additional care is needed in determining the solution in the vacuum region, which requires some gauge-fixing, as might be expected.
We paid particular attention to the form of the solution in the vacuum region and verified the validity of Birkhoff's theorem, the usual interpretation of which has recently been brought in question [44]. As pointed out in [44], Birkhoff's theorem is often used to support the common belief that in a vacuum, the gravitational field is only affected by the interior mass. We note that that this is incorrect.
The above investigations allowed us to re-examine critically the original theoretical arguments set out in [45,76,77] for the so-called R h = ct cosmological model, which has recently received considerable attention. After pointing out a number of objections to the R h = ct based on recent observational data, we consider in particular the central assumption underlying the original theoretical argument for the model, namely that the comoving Hubble distance should be constant. We demonstrate that this is not required, and so find no reliable theoretical basis for the R h = ct model.
These considerations in turn elucidated the behaviour of a number of 'horizons' during the general-relativistic evolution of homogeneous and isotropic cosmological models. In particular, we considered the evolution of an imaginary spherical boundary of radius a(t) that is comoving with the fluid and centred on some arbitrary origin. For a selection of analytical spatially-flat cosmological models, we compared a(t) to the Schwarzschild and Hubble radii. In the case of vanishing cosmological constant, we find the generic behaviour (both for dust and radiation models) that the comoving radius a(t) initially lies outside the Hubble radius, but eventually enters it at precisely the same moment as it exits the Schwarzschild radius; it is allowed to do the latter, since the fluid at the boundary is moving at speed c at this instant. Thus, one has two 'horizon crossings' taking place simultaneously and in opposite directions. In the case Λ = 0, one can obtain two positive solutions for R S (t), which correspond to the Schwarzschild radius and the de Sitter radius, respectively. One again finds that the comoving radius a(t) initially lies outside the Hubble radius and enters it at precisely the same moment as it exits the Schwarzschild radius, at which point the fluid at the boundary a(t) is moving at speed c. In the presence of non-zero Λ, however, the boundary a(t) later exits the Hubble radius again, at precisely the same moment that it enters the de Sitter radius, at which point the fluid at the boundary is again moving with speed c. This interesting behaviour is not usually pointed out in the literature.