Complete solutions to the metric of spherically collapsing dust in an expanding spacetime with a cosmological constant

We present semi-analytical solutions to the background equations describing the Lema\^itre-Tolman-Bondi (LTB) metric as well as the homogeneous Friedmann equations, in the presence of dust, curvature and a cosmological constant Lambda. For none of the presented solutions any numerical integration has to be performed. All presented solutions are given for expanding and collapsing phases, preserving continuity in time and radius. Hence, these solutions describe the complete space time of a collapsing spherical object in an expanding universe. In the appendix we present for completeness a solution of the Friedmann equations in the additional presence of radiation, only valid for the Robertson-Walker metric.


Introduction
Cosmic structures today have entered the non-linear regime. They can not on all scales be described by a linear perturbation theory on top of the Friedmann-Lemaître-Robertson-Walker (FLRW) metric [1][2][3][4][5]. The simplest step beyond linear perturbation theory is to look at separate patches, describing the evolution in each patch as a closed system, insensitive to other perturbations outside the patch. Over-densities can be studied in the approximation of spherical collapse, where underdensities expand and potentially become spherical voids. Either way, these spherically symmetric configurations, whether matching to a surrounding FLRW metric or not, are described by the Lemaître-Tolman-Bondi (LTB) metric [6][7][8], with S(r, t) = R ′ (r, t) The LTB metric reduces to the FLRW metric if one sets R(r, t) = r a(t) and E(r) = − 1 2 k r 2 . In Ref. [9] it was shown that when peculiar velocities are small, a seemingly non-linear solution to the metric becomes a linear perturbation on the FLRW metric in the Newtonian gauge. Here we focus on solutions to the metric in full generality.
Spherical collapse is studied for example for the formation of black holes as well as for determining, in cosmology, if and when an initially linear over-density produced during inflation collapses and decouples from the background expansion. In simplest approximation, one considers a homogeneous overclosed patch that expands and collapses, matched to an expanding background by a singular shell [10]. Choosing a continuous curvature profile in Eq. (2), allows for an exact solution without singular shells. Spherical collapse in either the approximation or the exact approach, gives insight in clustering of matter, and thereby has been related to presence of for example Dark Energy, amongst other possibilities [11][12][13][14][15][16][17][18][19][20][21][22][23][24]. One should note that initial velocities could be such that an overdensity evolves to become under-dense, but such a decaying mode corresponds to an inhomogeneous Big Bang, which is at tension with the inflationary paradigm in today's favoured model of cosmology.
The formation of voids is studied for two reasons. One reason is the role of voids in the process of structure formation [25][26][27][28], the other is the effect that unusually deep under-densities can have on our perception of Dark Energy. Some studies consider one large local void (size varying from tens to thousands of megaparsecs) , where others consider a distribution of many voids, the so called Swiss-Cheese universe [55][56][57].
In presence of only dust and curvature, analytical solutions to the LTB-equation (which we write down later) for t(R, r), i.e. time as a function of local angular scale factor R and coordinate radius r, are known in terms of hyperbolic functions (which become trigonometric functions in case of complex arguments). However, observations of distant Supernovae, of the Baryonic Acoustic Oscillations and of the distance to the surface of last scattering of Cosmic Microwave Background photons combined with locally observed expansion rate, demand the presence of a cosmological constant [58][59][60][61]. It is crucial to realize, although it is not the topic of this paper, that the presence of the cosmological constant is only then necessary to explain the observed geometrical distances, if one assumes that on large enough scales the universe is still properly described by the FLRW metric and one assumes that the angular diameter distance-redshift relation is correctly described by the background dynamics only.
In the presence of a cosmological constant, dust and curavture, the solution for t(R, r) is an elliptic integral. Therefore, in most works, authors elude numerical integration of the Einstein equations, although see Refs. [62,63]. However, if one for example wants to solve geodesic equations, one typically would perform a numerical integration of the geodesic equations over a numerical solution of the background. The unknown error in the numerical background solution can propagate into the solution to the geodesic equations, possibly leading to unreliable answers, or very slow and sometimes unstable codes.
In absence of a cosmological constant, but in presence of dust and curvature, the solution for t(R, r) with positive κ(r) is, while for κ(r) < 0, the solution becomes a trigonometric in stead of hyperbolic function if one propagates the sign of κ(r) correctly. In this case, one obtains the inverse solution R(r, t) by numerically inverting t(R, r), which can be done quickly and at ultimate accuracy, since dt dR is known by the definition of the Einstein equations. The thus obtained solution is accurate, fast, and allows for reliable and fast integration of geodesic equations [50].
The purpose of this paper is to provide all solutions to the background equations t(R, r) that have an initial singularity (Big Bang), in presence of dust, curvature and a cosmological constant, in terms of Elliptic Integrals in Carlson's symmetric form, which can be numerically evaluated as accurate and fast as any elementary function. Compared to the known exact solution in the case of no cosmological constant, the solutions presented in this paper are exact, since one eventually obtains R(r, t) by quick and reliable numerical inversion of t(R, r). Hence, throughout this paper no numerical integration is performed, and solutions are exact but only semi-analytical. See for comparison elliptic solutions involving the Weierstrass elliptic function in Refs. [64][65][66][67][68] and references therein. Note that these references only give t(R, r), which in the FLRW case is enough to solve for a(t), but which in the LTB case does not suffice: t(R, t) only straightforwardly leads to the angular scale factor R(r, t), but we also present for the first time the radial scale factor S(r, t), which is more involved as we shall see later.
The main improvement in this work in comparison with the existing literature, is the fact that we provide solutions for all functions appearing in the LTB metric, involving spatial derivatives ∂ r R(r, t), necessary for solving for example geodesic equations. We list the solutions for the limits where the local expansion transits to collapse, while a neighbouring shell continues to experience expansion, at the same time preserving continuity of all functions. Moreover, in the appendix we present linear expansions when either the cosmological constant or curvature is small, or both are small.
The solutions presented in this work allow for a plethora of applications. Let us list but a few. For example, the solutions can be used for: • any numerical work involving a general LCDM background expansion, • obtaining the exact metric around a collapsing structure • simulating the universe as observed by an observer in a large and deep void, in presence of a cosmological constant, allowing for a direct face off between the cosmological constant and the void, • studying the evolution of voids in structure formation, in a ΛCDM universe.
We release a numerical module, written in Fortran, that computes exact metric functions and derivatives for a given curvature profile, that can be easily implemented in any code, at http://web.physik.rwth-aachen.de/download/valkenburg/ColLambda/. A brief example of how to invoke this module is given in Appendix D. This work is organized as follows. In Section 2 we first list for reference the used Elliptic Integrals in Carlson's symmetric form. In Section 3 the LTB metric and its Einstein equations are discussed. Then in Section 4 we present one of the main results of this paper, being t(a) in terms of Carlson's elliptic integrals. Next, in Section 5 we solve analytically for the functions in the metric as a function of time t and scale factor a. In Section 6 we provide an example of an application of the solutions presented in this work. Finally we conclude in Section 7. In Appendix A we provide the asymptotic expansions of the solutions and in Appendix B we show the solution of the Friedmann equation in presence of radiation, matter, curvature and a cosmological constant in terms of Carlson's elliptic integrals.
Throughout this work, square roots of real quantities are taken to be positive, and for all fractional powers of complex numbers x we take the principle value of exp(ln x). Extra minus signs due to the possible crossing of branch cuts are written explicitly. We use units in which G N = c = 1. Overdots denote time derivatives, primes denote radial derivatives. Our notation follows mostly the notation used in for example Refs. [50,56,69].

Carlson's symmetric form of Elliptic Integrals
Before discussing solutions to the LTB equation, let us list some definitions for completeness. We take the following definitions of Carlson's symmetric form of elliptic integrals from Ref. [70], which are defined for {x, y, z} ∈ C | {x, y, z} / ∈ −∞, 0]. These can be evaluated using an iterative procedure, up to unlimited accuracy in very few steps, as explained in Ref. [71]. The definitions are valid for complex arguments, and in all these cases at most one argument is allowed to be zero.

Robertson-Walker and Lemaître-Tolman-Bondi metrics
The Einstein equations for the FLRW metric and the LTB metric can be written in the same form, with the difference that in the former case no other coordinate dependence than time dependence is present, and in the latter case the curvature and scale factor are both radius and time dependent and radiation is absent. The FLRW metric is, and the Friedmann equation for the FLRW metric is where as usual we define today by t = t 0 , a 0 = a(t 0 ) = 1, H 0 = H(t 0 ). The different components and their relative abundances are radiation Ω r , dust Ω m , curvature Ω k and the cosmological constant Ω Λ . The LTB metric is given by where S(r, t) = R ′ (r, t) R(r, t) =r a(r, t).
whereM 2 is an arbitrary parameter defining the length and mass scales, combined with the choice of units G N = c = 1. The Einstein equation leads to the LTB equation [8], where L(r) = r 0 dr M ′ (r) 1 + 2r 2 κ(r)M 2 with M (r) = 4π r 0 dr S(r, t)R 2 (r, t)ρ(r, t) being the total mass 1 inside radius r and ρ(r, t) is the local matter density. Now there are three functions of r that specify the problem, M (r), κ(r) and t BB (r), where the latter is the radially dependent Big Bang time, t BB (r) ≡ t(a = 0, r). One of these three can be fixed to an arbitrary function by redefining the coordinate r →r = f (r) for some monotonic function f (r), without changing the physical description. As discussed in Ref. [72], none of the three possible coordinate gauges where one of the functions is fixed to an arbitrary monotonic function captures all possible configurations. In this work we choose the gauge as follows. Demanding a strictly positive matter density, we have L ′ (r) ≥ 0. We choose L(r) = 4πM 2 r 3 /3, such that L ′ (r) > 0 but L ′ (r) = 0. This implies that are are no vacuum regions. Then we have In this coordinate gauge the LTB equation becomes, and the configuration is completely specified by the two functions κ(r) and t BB (r). The shortcoming of this gauge is, as mentioned above, that it does not allow for solutions with true vacuum over a non-zero range in r, for which M ′ (r) = 0 such that κ(r) → ∞.

Normalization
Normalizing a(r * , t 0 ) = 1 at a chosen {r * , t 0 }, we have which is regular for Λ → 0 and Ω k → 0. Also, we write H * ≡ H 2 (r * , t 0 ). One can choose an arbitrary r * at which to normalize, but one has to fix it once and for all. Clearly, in the FLRW case, the choice of r * is irrelevant since the r-depence of H(r, t) vanishes, and we find H 0 = H * .

Towards solving for t(a, r)
One can define in which case one retrieves the Friedmann equation in absence of radiation when one drops the r-dependences in equation (16) and when one normalizes a(r, t 0 ) = a(t 0 ) = 1. For a generic matter distribution in the LTB metric one however has a(r, t 0 ) = 1. Then, in the LTB metric, these three quantities are the relative content at a(t, r) = 1 in a shell at a given radius: dust, curvature and cosmological constant respectively. At t 0 the relative matter quantities are then depending on the value of a(r, t 0 ), e.g. Ω m a −3 (r, t 0 ) for matter. The reader should note that in this sign convention • Ω k > 0 corresponds to an open universe • and vice versa Ω k < 0 corresponds to a closed universe.
From this point on we will neglect radiation, although the reader is referred to Appendix B for a discussion of the solution in presence of radiation, only valid in the FLRW metric.
Writing A ≡ a(r, t), the general solution to t for the metric (both LTB and FLRW) is given by the integral, where we take the positive square root and the integral is performed at constant r. Note that the Big Bang time t BB (r) acts as an integration constant for the left hand side of this equation.
In the case of existence of at least one positive real root of the polynomial in the denominator (i.e. an over-closed universe (FLRW) or over-closed shell at radius r (LTB)), there are two solutions for t(A), one for the expanding and one for the collapsing phase. For example labeling the smallest positive real root U (r), such that the turning point in the expansion history lies at a(r, t) = U (r), the second (collapsing) solution is where the sign in front of the second integral is determined by whetherȧ changes sign (collapse) or not (continued expansion) at a = U . Throughout the rest of this work, we discard cases wherė a does not change sign, i.e., we only consider cases whereä is non-zero at a = U (r). The case with more than one positive real root then becomes irrelevant: whenȧ changes sign, by symmetry the contraction is identical to the expansion, and we only consider the branch of solutions that experiences an initial singularity (Big Bang). As a shrinks again towards 0, for each value of a the time derivative is exactly −1 times the time derivative for the same value of a during the expansion. The higher roots are never encountered.

Solving for the expansion rate
We rewrite Eq. (23) to for n = 3, where we rewrote the polynomial with y i the solutions to We now turn Eq. (25) into Carlson's symmetric form by subsequently making the change of variables a → c = 1 a and c → b = c − 1 A , such that we get These transformations are valid, since no rotations in the complex plane are involved, and none of the roots is transformed onto the path of integration; no branch points come to lie on the positive real axis. At most one root will be at the zero of the integration domain in the variable b, when we integrate ym 0 d a. This exception is allowed in the definition of the symmetric elliptic integrals in Section 2.
In terms of physics this is straight forward to see: for any real negative y m , we have 1 A − 1 ym > 0, i.e. the branch point lies on the negative real axis of b; for any real positive y m , we also have 1 A − 1 ym ≥ 0, since we integrate at most up to A = y m . The scale factor never grows beyond its smallest maximum for t ∈ R.
We now obtained the solution for the time as a function of scale factor expressed as a symmetric elliptic integral in Eq. (28), that is, One of the first three arguments of R J (x, y, z, p) is allowed to be zero, such that the limit of A → y 1 for positive real y 1 is trivial. In order to keep a connection to the well-known Friedmann equation, we so far used a notation in terms of H * and Ω i (r). However, as H * = H(r, t 0 ) is a function of r, it is more convenient to go back to the original form of Eq. (16) and recast Eq. (23) as where the roots y i and z i are related by a simple rescaling.

The roots z i
The roots z i are trivial, For the FLRW metric the roots y i obey the same expressions, replacing X, Y and Z by the corresponding Ω i .

The metric functions and their time derivatives
In this section we present the main result of this paper, the radial scale factor and radial derivatives of other metric functions. In the previous section we discussed the solution t(a). Since we know by definition the exact dt da , one can solve numerically for a(t) using a simple Newton-Raphson algorithm, obtaining a(t) at machine accuracy level at hardly any computational cost. Therefore, in the following we assume {r, t} as input parameters, and the function a(r, t) as a known function. We aim to express all solutions in terms of those quantities.
The functions appearing in the metric are R(r, t) and R ′ (r, t). The time derivatives of these functions are relevant when one wants to solve geodesic equations in this metric, which is why we discuss them here as well. We have R(r, t) ≡ ra(r, t), such that Comparing Eq. (37) to Eq. (38), we see that we only need the solution for a ′ (r, t) in order to be able to calculate all relevant quantities. However, the term 1 2H(r,t) asks for care to be taken at the transition from expansion to collapse, where H(r, t) = 0. We will show in the following that this limit is in fact regular, and that all metric functions remain properly defined throughout all the expansion and collapse history.

Spatial derivative of the scale factor during expansion
Since t is one of the orthogonal coordinates, we have ∂ r t ≡ 0, even if we solve for t by t = t(a(r, t), r). Therefore, The only non-trivial term in (40) is where we continue to use the notation of the previous section, A = a(r, t), X = 8π 3 , Y = 2κ(r) and Z = Λ 3M 2 . This expression is again an elliptic integral. We spare the reader the detailed steps, but the general procedure is very much like in Eq (28), however one not only substitutes a → c = 1 a and c → b = c − 1 A , but one also splits into partial fractions, to arrive at One can take this equation one step further, using Eq.
Hereby one eliminates one evaluation of the function R D (x, y, z), and more importantly, it reveals the kind of singularity that is encountered when one of the arguments goes to zero. Choosing z 1 to be the smallest positive real root, or as in the previous section U (r) = z 1 , we finally have the solution, with A = a(r, t), which can be recast as M a ′ (r, t) =ȧ(r, t)Q(r, t) + P (r, t), with as before A ≡ a(r, t).

Spatial derivative of the scale factor during collapse
During the collapsing phase, we havẽ We write explicitly the r-dependence in U (r) ≡ z 1 (r), to point out that the transition from expansion to collapse is a priori an r-dependent event. 2 Taking the derivative of this expression, we findM where the absolute value |ȧ| is to remind us that we take the positive root of (a − z 1 )(a − z 2 )(a − z 3 ). When we realize that the r-dependence in z 1 (r) is entirely specified by Y (r) as that is the only r-dependent function in all integrals, such that ∂ r z 1 = Y ′ (r)∂ Y z 1 , the first term becomes, As . Using this relation, together with the solution to the same integral in Section 5.1, and some more simple algebra, we arrive at where t U ≡ t(a = z 1 , r), the subscript 'coll' denotes that this expression is valid during a collapsing phase, Q(r, t) and P (r, t) are defined in Eqs. (45,46), and Q(r, t U ) is evaluated by replacing A → z 1 in Eq. (46). It should be understood now that with expressions (44) and (51), a ′ (r, t) is finite and continuous atȧ → 0.

The spatial derivative of the Hubble parameter atȧ = 0
During expansion, the expression for H ′ (r, t) in Eq (38) is regular. Atȧ → 0, we can now insert Eq. 44 for a ′ (r, t), to find after some manipulations, which is pefectly regular for the whole domain a(r, t) ∈ {0, z 1 }. One extra minus sign appeared in the last expression, following from a−z1 √ a−z1 = − √ a − z 1 with a ≤ z 1 for all a.

The spatial derivative of the Hubble parameter during collapse
Combining Eq. (51) and Eq. (52), we find during a collapsing phase, preserving continuity atȧ = 0.

Example of application
In Figure 1 we show an example of an application of the solutions presented in this work. The figure shows a comparison between two different density profiles, both parametrized by where κ b is defined in Eq. (17), and with which is the third order of the interpolating function W n (x, α) which interpolates from 1 to 0 in the interval α < x < 1, while remaining C n everywhere. We introduce the interpolating function W n (x, α) in Appendix C. For 0 < α < 1 this function is C 3 on r ∈ [0, ∞ . Hence, all functions in the metric are C 2 everywhere, including at the center (r = 0) and at the matching to FLRW at r = L, guaranteeing a finite and continuous Ricci scalar everywhere. By construction all functions shown in the figure in fact have continuous first derivatives in r, even if by eye it may seem otherwise. The left column in Figure 1 shows quantities for a profile with α = 0, the right column shows the same quantities for α = 3 4 . In the top row we see the curvature profiles as a function of comoving radius, which are time independent. Both profiles share a same matter density at the centre where r = 0 and at the outer radii r > L, however the difference in shape of profile leads to a different total amount of matter inside the over-density. For α = 3 4 , the over-density possesses a larger region with large closed curvature. The second row shows a time-dependent auxiliary radius,  Figure 1: Comparison of the two distinct overdensities. Curvature profile (top row) as a function of comoving radius, the auxiliary FLRW radius r FLRW as a function of radius and time (second row), the matter density as a function of the auxiliary radius r FLRW and time (third row) and finally the time derivative of the metric function S(r, t) as a function of auxiliary radius r FLRW and time (bottom row). Both over-densities are matched to the same homogeneous ΛCDM-universe, at the same comoving radius r = 1 Mpc. The difference between the over-densities is in the value α in Eq. (55), determining the range in r over which the curvature profile falls back to the background value and thereby determining the total mass in the over-density. In all graphs, the time coordinate is represented by the colour of the curves, evolving from red to white to blue (from black to white to black in black-and-white print), indicating time varying from today (t 0 = 13.3 Gyr) to some moment in the past (t 0 = 2.8 Gyr). Additionally, labels inside the graphs indicate times corresponding to different curves. which keeps the length measure along this radius at a given time constant, ds 2 = S(r, t)dr 2 = dr 2 FLRW , but not constant in time. This radius can be interpreted loosely as what an FLRW observer would see. For these scenarios it turns out that roughly r FLRW ≃ R ′ (r, t)/a(∞, t), which is its definition in Refs [55,57].
In the third row we show the relative matter density, normalized to the matter density of the surrounding homogeneous cosmology. While at the centre the time evolution is the same for both profiles, the surrounding under-dense (but still closed curved) shell differs largely between the two cases. Since κ(r) stays large and negative up to higher radius in the α = 3 4 -case, compared to the α = 0-case, a larger range in r is present for which shells have a collapsing solution and actually experience the collapse. Therefore, the range in r in which the shells expand rapidly and become more and more under-dense, is smaller for the α = 3 4 -case. Hence, the resulting surrounding under-density is less dense and more emphasized in the α = 3 4 -case than in the α = 0-case. For a radial null geodesic, the trajectory is defined by, This is why, for illustration, we plotṠ(r, t) in the fourth row in Figure 1.
Obviously, as the effect of the spherical collapse on the outer radii is more violent for the α = 3 4 -case, the red shift that a photon experiences by passing through that region is larger than in the α = 0-case. Even though this observation is not enough to draw conclusions about photon geodesics and collapsing structures, it illustrates how the solutions in this work can be used to further asses the importance of the initial distribution of matter in the line of sight, on distant observations. Note that all quantities remain perfectly smooth and continuous at the transition from expansion to collapse, which in the third and fourth row in Figure 1 occurs roughly where each curve crosses the level of the background, i.e. the level of the same curve at r > L.
The solutions presented in this work allow for practically instantaneous calculation of the quantities presented in Figure 1. We release a module, written in Fortran, which returns all metric functions and derivatives thereof as a function of time for given functions κ(r) and t BB (r), and given cosmological parameters. The module is released at http://web.physik.rwth-aachen.de/download/valkenburg/ColLambda/.

Conclusion
We have presented an as complete as possible overview of the solutions to the Einstein equations governing the Lemaître-Tolman-Bondi metric, including fully continuous solutions for collapsing over-densities surrounded by an expanding universe. The solutions are written in terms of Elliptic Integrals in Carlson's symmetric form, which allow for fast numerical evaluation of the solutions at machine accuracy level. The solutions to all metric functions involve the numerical inversion of one function, t(a), whose derivative is explicitly known a priori, therefore allowing for inversion at machine accuracy level while remaining sufficiently fast.
We finished with a brief example of how these solutions can be applied. These solutions could improve the accuracy and speed of many analyses involving structure formation and inhomogeneous cosmologies.

A Asymptotic expansions
The solutions of the scenario where all three components, dust, curvature and the cosmological constant are non-zero and not asymptotically small, are given in the main body of this paper. For several purposes, a not unimportant one being numerical accuracy, asymptotic expansions are useful. By asymptotic expansions we mean the solutions for a given size of the scale factor, where one or more of the constituents contribute only marginally to the result.
In the following we use a looser definition of Ω i , where at each size of the scale factor Ω i denotes the fractional contribution to H(a). That is, for matter for example we re-define Ω m (a, r) = 8πM 2 3H(a) 2 1 a 3 and so on. We focus only on solutions with a Big Bang, and only solutions with a non-zero matter (dust) content. Then, given Eqs. (16) and (30), which we repeat here, we see that for any choice of κ(r) and Λ, prior to some initial time the equation is dominated by matter. Hence the integral that gives t(a) always has a non-negligible contribution from the matter content, even for a so large that 1 a 3 ≪ Λ 8πM 2 . Therefore we do not consider asymptotic expansions for Ω m ≪ 1.
Similarly, when expanding the integrand in Eq. (60) for small omega's, it only makes sense to expand in Ω k if Ω k ≪ Ω m , regardless of Ω Λ . To summarize, one can use expansions under the conditions, expand in Ω k when 3κ(r)a 4π <ǫ, expand in Ω Λ when Λa 3 8πM 2 + 6κ(r)a <ǫ, guaranteeing that the expansion is valid throughout the whole integration from 0 to a, for t(a) . When one uses a linear expansion, one should set ǫ to ǫ = √ η, with η the desired accuracy, such that the error is O(ǫ 2 ) = O(η). In the case of numerical computation, one sets η to the machine precision, which for double precision (64 bit floating point) is η = 10 −16 , such that ǫ = 10 −8 . That is, in double precision one has approximately 16 significant digits.
For small Λ, the integral in Eq. (60) becomes which reduces to the well-known result as Λ → 0. The turning point from expansion to collapse lies at a = − 4π 3κ(r) if this expression is positive. Care must be taken with the branch cut of the square root here, as for negative κ, 1 For κ − 9 2 the minus signs cancel, and no care has to be taken.
To obtain the derivative of time with respect to radius r in this case, one must first take the derivative of the full expression, and then expand, to arrive at, A.
2 Ω k ≪ 1, Ω m = 0 = Ω Λ For small κ(r), the integral in Eq. (60) becomes The first integral is in principle elliptic, and the full solution in Eq. (31) is applicable when one defines the correct roots z i , however this special case is has a known solution in terms of sinh(x) and can be found in the literature. The second integral is of course the integral that is solved in the main text for t ′ (a, r) for the special case where κ(r) = 0. Hence, where x i are the three solutions of 8π 3 + Λ 3M 2x 3 i = 0. If one allows for a negative Λ, also this scenario can have a postive real x i at which the transition from expansion to collapse occurs. Of course care has to be taken at the limit ofȧ → 0, identical to what is described in the main text concerning a ′ (r, t).
Unfortunately, an expansion to obtain ∂ r [t(A) − t BB (r)] in this asymptotic region is not trivial, where at this moment we do not know of a solution of the last integral in terms of symmetric elliptic integrals, so we leave that integral for future work.
The simplest of the expansions is the scenario with small κ(r) and Λ. The integral becomes, with the spatial derivative,

B Solution in presence of radiation
Writing A ≡ a(r, t), the general solution to t for the Friedmann equations in presence of radiation, matter, curvature and a cosmological constant is given by the integral, which is the same as Eq. (25) for n = 4. We split the integral in two parts, where it is most convenient to take z 1 to be the smallest positive real root, if any root is positive and real (otherwise any root will do), In the second integral we can substitute where if i, j = 2, 3, then k, l = 1, 4). The last equality in Eq. (72) is proven in Ref. [73], and the equality is invariant under the ordering of the roots z m , that is, invariant under the choice for z 4 .
The first integral in Eq. (71) takes another road. Following Ref. [74] we make the change of variables a → A b b+1 ,

B.1 Roots z i
In order to write down the roots {z i }, let us use the following definitions, The roots are C The interpolating function W n (x, a) We define the function W n (x, a) as, W n (x, α) ≡ 1 for x < α W 2 (x, α) ≡ 1 2 + 1 2 sin π 1 2 − x−α 1−α W n+1 (x, α) ≡ which is C n everywhere. Note the absolute value inside the integral, which takes care of mirroring two images of the function that interpolates from one to zero, in order to have a function that goes from zero to one to zero. Also, inside the integral the functions are evaluated for α = 0, because the integration limits are always 0 < {x ′ , x ′′ } < 1. One can construct a similar function taking polynomials in stead of sinusoids, by starting with W 0 (x, α) = (1 − x)/(1 − α) for α < x < 1.
We constructed this function by starting with the function f (x) = 1 + sin(x), which equals zero and has a zero derivate at x = −π/2 and at x = 3π/2. Integrating this function over x then leads to a function which has zero first and second derivatives at x = −π/2 and at x = 3π/2. Next one matches this function too a mirror image, normalizes it to zero at both end points of the domain and integrates again, such that one ends up with a function that has zero first, second and third derivatives at the domain borders. One can continue this scheme forever, as explicated in Eq. (88).

D Usage of the numerical module ColLambda
The module ColLambda is written in Fortran90, and can be downloaded from http://web.physik.rwth-aachen.de/download/valkenburg/ColLambda/. Compilation instruc-tions can be found at that website. The module can be included by the statement, As a consequence, any other parameters on which k(r) may depend, such as the maximum size L, maximum curvature κ max , or anything else, must be global variables which the function kofr(r) can access without receiving them as arguments. The values that the subroutine returns have the following notation: Rpd =Ṙ ′ (r, t) Rpdd =R ′ (r, t) S = S(r, t) Sd =Ṡ(r, t) Sdd =S(r, t) a = a(r, t) ap = a ′ (r, t) apd =ȧ ′ (r, t) add =ä(r, t) apdd =ä ′ (r, t) H = H(r, t) Hp = H ′ (r, t) where t U (r) denotes the time at whichȧ(r, t U ) = 0, if it exists. If it exists, the local singularity is reached at t = t BB (r) + 2t U , that is, a(r, t BB (r)) = a(r, t BB (r) + 2t U ) = 0. If this time does not exist, i.e. when the solution is ever expanding, tturn will be set to 10 30 . An example call, where the user has set the normalization variables and has defined the necessary functions κ(r), t BB (r) and their first derivatives, to set myRltb to R(r, t), myS to S(r, t) and myt to t U (r), would look like: call lltb_functions(H0_inf, Lambda, kofr, dkdr, tbbofr, dtbbdr, Mtilde, & r, t, & Rltb=myRltb, S=myS,tturn=myt)