Towards the Raychaudhuri Equation Beyond General Relativity
Daniel J Burger Saurya Das S. Shajidul Haque Nathan Moynihan Bret Underwood



























    burgerj.daan@gmail.comDepartment of Applied Mathematics and Mathematics,
    
    University of Cape Town, 
    
    Cape Town, South Africasaurya.das@uleth.caDepartment of Physics and Astronomy,
    
    University of Lethbridge, 4401 University Drive,
    
    Lethbridge, Alberta T1K 3M4, Canadashajid.haque@uwindsor.caDepartment of Physics,
    
    University of Windsor,
    
    Windsor, Ontario N9A 0C5, Canadanathantmoynihan@gmail.comDepartment of Applied Mathematics and Mathematics,
    
    University of Cape Town, 
    
    Cape Town, South Africaunderwbj@plu.eduDepartment of Physics,
    
    Pacific Lutheran University,
    
    Tacoma, WA 98447
    
    In General Relativity, gravity is universally attractive, 
    a feature embodied by the Raychaudhuri equation which requires that
    the expansion of a congruence of geodesics is always non-increasing,
    as long as matter obeys the strong or weak energy conditions.
    This behavior of geodesics is an important ingredient in general proofs of singularity theorems,
    which show that many spacetimes are singular in the sense of being geodesically incomplete and 
    suggest that General Relativity is itself incomplete.
    It is possible that alternative theories of gravity, which reduce to General Relativity in some limit,
    can resolve these singularities, so it is of interest to consider how the behavior
    of geodesics is modified in these frameworks.
    We compute the leading corrections to the Raychaudhuri equation for the expansion
    due to models in string theory, braneworld gravity, f(R) theories,
    and Loop Quantum Cosmology, for cosmological and black hole backgrounds,
    and show that while in most cases geodesic convergence is reinforced, in a few cases
    terms representing repulsion arise, weakening geodesic convergence and 
    thereby the conclusions of the singularity theorems. 
    Towards the
    Raychaudhuri Equation Beyond General Relativity
    
        Bret Underwood
        December 30, 2023
    ============================================================
    
    
    
    § INTRODUCTION
    
    
    General Relativity (GR) is widely expected to be an incomplete theory of dynamical spacetime.
    One reason for this is that GR famously predicts its own demise through the existence
    of singularities, as demonstrated by the singularity theorems <cit.>.
    An essential physical ingredient of the singularity theorems is that gravity 
    is attractive, so that congruences of convergent timelike and null geodesics develop singularities
    in finite proper (affine) time.
    More specifically, geodesic congruences with timelike u^M(τ) and null n^M(λ) tangent
    vector fields are characterized by their expansion θ≡∇_M u^M, θ̂≡∇_M n^M,
    respectively, which satisfy the Raychaudhuri equations <cit.>dθ/dτ = -θ^2/D-1-R_MN u^M u^N + ...       
    dθ̂/dλ = -θ̂^2/D-2-R_MN n^M n^N + ...        
    for spacetime dimension D, where the additional +... terms are non-positive; see
    Appendix <ref> for details.
    
    For the timelike Raychaudhuri equation, if the so-called “convergence condition”
    R_MN u^M u^N ≥ 0
    
    holds, then the expansion of a congruence of geodesics is non-increasing.
    Specifically, an initially converging congruence θ_i < 0 develops a singularity θ→ -∞
    in finite proper time τ∼ |θ_i|^-1. A similar singularity θ̂→ -∞ develops for null convergences
    in finite affine time λ∼ |θ̂_i|^-1 if the null convergence condition R_MN n^M n^N ≥ 0 is satisfied.
    These singularities in the expansion don't necessarily imply a pathology of spacetime themselves; such caustics appear in Minkowski
    spacetime, for example <cit.>.
    However, the existence of strictly non-positive contributions to the right-hand side of Eq.(<ref>),
    when combined with other global conditions on the spacetime manifold, form the basis for the general
    existence of singularities in cosmological and black hole spacetimes <cit.>.
    
    
    Within GR it is possible to convert the convergence conditions into energy conditions on the types of matter.
    We can use the (trace-reversed) Einstein equation
    
    R_MN = κ_D^2 (T_MN - 1/D-2 g_MN T^M_M)
    
    (where κ_D^2 ≡ 16π G_D, G_D being the D-dimensional Newton's constant)
    to substitute in Eq.(<ref>)
    
        dθ/dτ = -κ_D^2 (T_MN - 1/D-2 g_MN T^M_M ) u^M u^N + ...       
           dθ̂/dλ = -κ_D^2 T_MN n^M n^N + ...
    
    where the +... terms are again strictly non-positive
    terms.
    Thus, all of the terms on the right-hand side of (<ref>) are non-positive provided that matter satisfies the respective energy conditions
    (T_MN - 1/D-2 g_MN T^M_M ) u^M u^N ≥ 0       
    
    T_MN n^M n^N ≥ 0       
    For an isotropic perfect-fluid energy momentum tensor with energy density ρ and pressure p, these conditions translate
    into
    ρ + D-1/D-3 p ≥ 0       
    ρ + p ≥ 0       
    Most known classical matter obeys both the null and strong energy conditions, while vacuum energy p = -ρ violates the strong
    energy condition but still saturates the null energy condition.
    It may possible to violate the null energy condition with exotic forms of matter <cit.>, non-minimal coupling <cit.>,
    or quantum gravity effects <cit.>, though these approaches often face challenges that we will not explore further here.
    
    However, since we expect corrections to GR of some form, we do not expect the Einstein equations (<ref>) to always hold.
    Corrections to Einstein's equations may make it possible to violate the convergence conditions 
    R_MN u^M u^N < 0 and/or R_MN n^M n^N < 0 without violating
    the energy conditions (<ref>).
    In particular, many corrections to GR appear perturbatively in the form
    
    R_MN = κ_D^2(T_MN - 1/D-2 g_MN T^M_M ) + λ H_MN  ,
    
    where λ controls the strength of the corrections and H_MN is a tensor that contains contributions
    from the metric, curvature, energy-momentum tensor, or additional fields.
    The additional term in Eq.(<ref>) in turn shows up as an additional term on the right-hand side
    of the Raychaudhuri equation
    
        dθ/dτ = -κ_D^2 (T_MN - 1/D-2 g_MN T^M_M ) u^M u^N + λ H_MN u^M u^N+ ...       
           dθ̂/dλ = -κ_D^2 T_MN n^M n^N + λ H_MNn^M n^N+ ...
    
    where it may in principle contribute with any sign.
    In this paper, we will examine corrections to the Raychaudhuri equations
    from four frameworks for corrections to GR: string theory, braneworld gravity, f(R) theories, and
    Loop Quantum Cosmology, for cosmological (and in some cases black hole) backgrounds
    (see also <cit.> for a similar analysis in infinite derivative gravity).
    
    While we are motivated by the existence of singularities in these backgrounds, and the promise of these
    alternatives to GR for resolving the singularities, we will not attempt to prove the absence of singularities
    in this paper. Indeed, a significant amount of work has shown that finding realistic singularity-free 
    spacetimes under computational control is quite challenging, and we expect that true singularity resolution
    will require physics beyond the perturbative approach of (<ref>).
    Instead, we aim for a more modest goal, that of finding corrections to GR that lead to potentially positive
    terms on the right-hand side of the Raychaudhuri equation Eq.(<ref>) for timelike and null
    geodesics – a necessary, but far from sufficient, condition for ultimately resolving spacetime singularities.
    
    
    If indeed such positive terms are found (as we shall see for a handful of cases), further 
    analysis would be required to see, whether:
    (a) a generalization of the singularity theorems holds for these cases, or
    (b) another criterion for singularity (in lieu of geodesic incompleteness)
    can be applied. If none of the above holds, one  
    would be forced to conclude that these spacetimes are indeed non-singular. Such an analysis is beyond the scope of this work, however.
    
    
    In Section <ref> we outline a general set of corrections to GR from string theory in D dimensions,
    and compute the corrections to the Raychaudhuri equation for black hole and cosmological spacetimes.
    In Section <ref> we compute the corrections to the 4-dimensional induced Einstein equation for the braneworld scenario
    in a cosmological spacetime.
    In Section <ref> we compute the corrections to the Einstein equation for so-called f(R) theories, and determine the
    form of the corrections in a cosmological background.
    In Section <ref> we compute the corrections to the Raychaudhuri equation in a cosmological background from Loop
    Quantum Gravity.
    In Section <ref> we conclude with some comments on the potential for theories beyond GR to resolve singularities.
    
    
    
    § STRING THEORY AND GAUSS-BONNET CORRECTIONS
    
    
    In string theory, higher-order α' corrections induce corrections to the action and, correspondingly, to the Einstein Equations.
    In particular, the corrected action at leading order in α' for bosonic, heterotic, type IIA/IIB takes the form <cit.>
    S    =   1/κ_D^2∫ d^Dx √(-g) e^-2ϕ[R_D + 4 (∂ϕ)^2 + λ R_MNPQR^MNPQ + 𝒪(α'^2)] 
    
    where we are ignoring the antisymmetric rank two tensor B_MN, ϕ is the dilaton, R_MNPQ is the D-dimensional Riemann tensor, 
    
    and where
    [The leading α' corrections for type II theories occurs first at 𝒪(α'^3)<cit.>.]λ   =   1/2α'    1/4α'     0    
    Unfortunately, (<ref>) contains higher derivative terms in its equation of motion.
    However, it turns out that (<ref>) is ambiguous up to field redefinitions of the fields g_μν,ϕ
    to next order in α'. It is possible to use these field redefinitions to remove the higher order terms in the
    equations of motion[See the discussion in <cit.>.], giving rise to the action
    
    S_mod = 1/κ_D^2∫ d^Dx √(-g) e^-2ϕ      { R_D + 4 (∂ϕ)^2 +1/2λ[ R_GB^2 + 16 (R^MN - g^MN R_D) ∂_M ϕ∂_N ϕ. . 
          .. -16 ∇^2 ϕ (∂ϕ)^2 + 16 (∂ϕ)^4] + 𝒪(α'^2)} + ℒ_m
    
    where we have allowed for additional matter (including potentially a D-dimensional cosmological constant) through
    ℒ_m, and  R_GB^2 is the 2nd order Gauss-Bonnet combination
     
    R_GB^2 = R_MNPQ R^MNPQ - 4 R_MNR^MN + R_D^2   .
    
    Setting the dilaton to a constant, we are left with the usual Ricci curvature term and the quadratic Gauss-Bonnet
    term in our action
    
    S_EGB = 1/κ_D^2∫ d^Dx √(-g)[R_D + 1/2λ R_GB^2] + ℒ_m  .
    
    We will refer to this as the Einstein-Gauss-Bonnet (EGB) gravity action.
    The resulting equations of motion are
    
    R_MN - 1/2 g_MN R_D = κ_D^2 T_MN + 1/2λĤ_MN  ,
    
    where
    Ĥ_MN =       g_MN/2 R_GB^2 - 2 R_D R_MN + 4 R_MA R^A_N + 4 R^ABR_AMBN- 2 R_MABC R_N^ABC  .
    
    In order to put (<ref>) in the appropriate form relevant for use the Raychaudhuri equation[See <cit.> for an alternative
    method of deriving a generalized Raychaudhuri equation for EGB gravity.], we will
    trace-reverse, giving
    
    R_MN = κ_D^2 (T_MN - g_MN/D-2 T^M_M) + λ/2 H_MN  ,
    
    where now
    
    H_MN = g_MN/D-2 R_GB^2 - 2 R_D R_MN + 4 R_MA R^A_N + 4 R^ABR_AMBN- 2 R_MABC R_N^ABC  .
    
    
    More generally, the terms above are just the leading terms of the more generic Lanczos-Lovelock extensions of gravity <cit.> (see <cit.> for a review).
    The Lanczos-Lovelock extensions include the set of additional terms that can be added to the gravitational action and still
    lead to 2nd order equations of motion.
    In 4 spacetime dimensions, it can be shown that the Ricci scalar and cosmological constant are the only non-topological
    terms that can be added.
    In particular, the correction term in (<ref>) is purely topological
    in 4-dimensions (it is the Euler characteristic χ of the spacetime), not contributing to the equation of motion.
    For higher dimensions D>4, there is a finite series of additional terms, increasing in powers of R, R_MN, R^M_NPQ
    (terminating at some order for a given D).
    We will just focus on this leading order term for now, but will keep in mind that there can be additional corrections to consider.
    
    In order to determine the form of these corrections it is necessary to calculate H_MN for specific backgrounds.
    For simplicity, we will restrict ourselves to corrections to black hole and cosmological backgrounds; 
    in the subsections that follow, we will examine the corrections for these backgrounds in more detail.
    
    
    
     §.§ Black Holes
    
    
    The Gauss-Bonnet correction terms (<ref>) 
    appear as part of a pertubative series of corrections, so we expect solutions for some given matter content
    to also be described by a pertubative series in λ: g_MN(λ) = g_MN^(0) + λ g_MN^(1) + ..., where g_MN^(0) is the
    uncorrected General Relativity solution, g_MN^(1) is the first-order correction, and so on.
    It is remarkable that exact solutions g_MN(λ) of Einstein-Gauss-Bonnet gravity for black hole backgrounds
    are known to all orders in λ<cit.>.
    However, since we expect the Gauss-Bonnet corrections to be only the first term in a series of corrections arising from string theory,
    we cannot trust these exact solutions beyond 𝒪(λ).
    One can therefore consider the effects of the Gauss-Bonnet correction terms on the Raychaudhuri equations up to 𝒪(λ) 
    for black hole backgrounds in two ways:
    
      *  Evaluate  R_MN u^M u^N ∼λ/2 H[g^(0)]_MN u^M u^N on the uncorrected Schwarzschild black hole metric.
    
    
    
      *  Evaluate R[g(λ)]_MN u^M u^N on the known exact metric g_MN(λ) for black hole backgrounds,
    then expand the result to 𝒪(λ).
    
    
    
    (with similar expressions for the null Raychaudhuri equation).
    
    
    
    While both approaches should give identical results, up to 𝒪(λ), we will consider and compare both approaches for completeness.
    
    
    
      §.§.§ Perturbative Black Hole Corrections
    
    
    Since we are considering actions with a constant dilaton and zero form fields, we will restrict ourselves to
    pure gravity solutions.
    We thus first begin by considering the perturbative correction terms H_MN to a D-dimensional black hole solution of the vacuum Einstein equations
    
        ds^2 = -f(r) dt^2 + dr^2/f(r) + r^2 dΩ_D-2^2  ,
    
    where dΩ_D-2^2= ĝ_ijdθ^i dθ^j is the metric of a (D-2)-dimensional
    constant curvature manifold with unit radius (such as a sphere).
    The zeroth order black hole solution takes the form
    
        f(r) = 1-μ/r^D-3  ,
    
    where μ is related to the mass M of the black hole by
    
        M = (D-2)A_D-2/2κ_D^2μ
    
    and A_D-2 is the area of a unit (D-2) sphere.
    Note that for D=4 we have μ = κ_D^2 M/4π.
    
    While the Ricci tensor and scalar vanish R_MN = 0 = R_D, as is expected for a D-dimensional 
    Schwarzschild background, the Riemann tensor does not, and gives the only
    non-zero contribution to the 
    Gauss-Bonnet scalar
    
    R_GB^2 = R_ABCDR^ABCD = (D-3)(D-2)^2(D-1) μ^2/r^2D-2 ,
    
    and corrections to the Einstein equation,
    
    H_MN = g_MN/D-2[R_ABCDR^ABCD] - 2 R_MABCR_N^ABC  .
    
    We are primarily concerned with the correction terms 
    in the (tt) and (rr) directions,
    
    H_tt   =   f(r)μ^2/r^2D-2 (D-4)(D-3)(D-2)(D-1) ; 
    
    H_rr   =    -1/f(r)μ^2/r^2D-2 (D-4)(D-3)(D-2)(D-1) .
    
    Note that these correction terms vanish identically H_MN=0 for
    D ≤ 4, which is what we expect since the Gauss-Bonnet correction is purely topological for D≤ 4, and only acts as a dynamical
    correction for D > 4.
    
    We now consider a radial affine null tangent vector
    
    n^M = (1/f(r),± 1,0⃗)  ,
    
    where ± 1 corresponds to radially outgoing/ingoing.
    The null Raychaudhuri equation takes the form
    dθ̂/dλ = -θ̂^2/D-2-|σ̂|^2 - R_MNn^M n^N  .
    
    Our corrections due to the Gauss-Bonnet term appear on the right hand side, as (recall that T_MN = 0)
    
    R_MNn^M n^N = λ/2 H_MNn^M n^N = λ/2 H_tt n^t n^t + λ/2 H_rr n^r n^r = 0  .
    
    Remarkably, the Gauss-Bonnet corrections  vanish identically everywhere for null rays, implying that the null Raychaudhuri equation for black holes is uncorrected to leading order in λ.
    
    Finally, consider a timelike geodesic described by the tangent vector
    
    u^M = (1/f(r),-(μ/r^D-3)^1/2,0⃗) .
    
    The timelike Raychaudhuri equation takes the form
    dθ/dτ = -θ^2/D-1 - |σ|^2 - R_MN u^M u^N  ,
    
    where again the corrections to the Raychaudhuri equation due to the Gauss-Bonnet corrections come from the last term
    
    R_MN u^M u^N = λ/2H_MN u^M u^N = λ/2H_tt u^t u^t + λ/2H_rr u^r u^r = λ(D-4)(D-3)(D-2)(D-1)/2μ^2/r^2D-2  .
    
    There are a few things to notice about this correction term. First, it does not seem to suffer from any pathologies due to the coordinate singularity of our coordinate system at the horizon; thus, we expect it to be valid throughout the entire spacetime (although the coordinate r will require careful interpretation).
    Second, notice that the corrections are  manifestly positive, thus the corrections make the caustic/conjugate point at
    the location of the putative singularity at r=0 worse, not better!
    
    We have shown that the perturbative EGB corrections to
    the Einstein equations do not provide divergence terms
    in either the null or timelike Raychaudhuri equations;
    in fact, convergence is strengthened in the timelike case.
    These results above are not particularly surprising 
    since it is known that black hole solutions in Einstein-Gauss-Bonnet gravity still posses 
    singularities <cit.>,
    as we will now examine in more detail.
    
    
    
    
      §.§.§ Exact Black Hole Solutions
    
    
    In the previous section we considered how the Gauss-Bonnet curvature squared terms
    lead to perturbative corrections to the Raychaudhuri equation for pure-Einstein gravity black hole 
    solutions.
    However, exact black hole solutions for Einstein-Gauss-Bonnet gravity are well-known
    <cit.>
    (see also <cit.>), so it is also possible to calculate the right-hand-side of the Raychaudhuri
    equation for these fully backreacted solutions.
    
    First, let us review the known black hole solutions of D-dimensional Einstein-Gauss-Bonnet gravity
    <cit.>.
    As before, we will work with the general spherically symmetric metric
    
    ds^2_D = -f(r) dt^2 + dr^2/f(r) + r^2 dΩ^2_D-2  .
    
    Solutions to the Einstein-Gauss-Bonnet equations of motion (<ref>)
    are <cit.> (see also <cit.>)
    
    f(r) = 1 + r^2/λ̂ + σr^2/λ̂√(1+2λ̂μ/r^D-1)  ,
    
    where σ = ± 1 labels different branches of solutions, and we have defined
    λ̂= (D-4)(D-3) λ for convenience.
    
    Consider first the σ = -1 branch. To lowest order in λ, the metric is
    simply that of a Schwarzschild black hole, and is asymptotically flat as r→∞,
    reducing to the usual Einstein gravity solution. For this reason, this is usually called the
    “Einstein branch.”
    In contrast, the σ = +1 branch,
    often called the “Gauss-Bonnet branch,”
    has a non-trivial vacuum structure for μ = 0
    f(r)|_μ =0 = 1 + 2r^2/λ̂  ,
    
    corresponding to anti-deSitter space for λ > 0
    with a large negative cosmological constant Λ_eff∼ -λ̂^-1
    (or deSitter space for λ < 0).
    For non-zero mass, 
    the metric in this branch resembles that of Schwarzschild-anti-deSitter space with a negative
    mass
    
    f(r) ≈ 1 +μ/r^D-3 + 2r^2/λ̂ + 𝒪(λ^2)  .
    
    The analysis of the stability of this branch requires some care (see
    <cit.>).
    Note that since the curvature scale of the corresponding
    AdS space is non-perturbative in λ, it is
    doubtful that we can trust these solutions as solutions to perturbative string theory.
    Further, for λ > 0 there is a naked
    singularity at the origin <cit.>.
    For these reasons, we will restrict our analysis to the Einstein branch σ = -1.
    
    Despite the curvature squared term in EGB gravity,
    the Einstein branch solutions (<ref>,<ref>)
    still have a curvature singularity at the origin.
    This singularity is surrounded by a horizon located
    at r_h, given by the roots
    of the polynomial <cit.>λ̂r_h^D-5 + 2 r_h^D-3 = 2μ  .
    
    For D > 5, this horizon always exists; however,
    for D=5, the existence of the horizon is
    guaranteed only for μ > λ̂/2; thus, for
    “microscopic” black holes, even the EGB solutions
    posses a naked singularity.
    Since the exact EGB black hole solutions still posses a curvature singularity at the origin,
    we are not particularly surprised by our perturbative result from
    the previous subsection indicating that the corrections to the null and timelike Raychaudhuri equations do not allow for terms that could prevent the formation of conjugate points.
    
    Even though we have the exact EGB solutions, which do possess a curvature singularity, in hand, let us nevertheless compute the contribution of EGB gravity to the right-hand side of the Raychaudhuri equation to explore
    the way in which the EGB corrections can affect the formation of conjugate points.
    
    Note that since the solutions (<ref>,<ref>)
    are  exact solutions to the corrected Einstein equations
    
    R_MN - 1/2 g_MN R = κ_D^2 T_MN + 1/2λ H_MN  ,
    
    we can compute the curvature term in the null (timelike) Raychaudhuri equation
    R_MN N^M N^N (R_MN u^M u^N respectively) directly,
    without needing to compute the quadratic curvature terms.
    
    For the metric (<ref>), we have
    
    R_tt   =   1/2 f(r) (f”(r) + (D-2)/r f(r))  ; 
    
    R_rr   =    -1/21/f(r)(f”(r) + (D-2)/r f(r))  ,
    
    where a prime ' denotes a derivative with respect to r.
    Note that these
    vanish for a Schwarzschild solution f(r) = 1-μ/r^D-3,
    as expected.
    
    A radial affine null vector has the same form as
    in the Schwarzschild case
    
    n^M = (1/f(r),± 1,0⃗)  ,
    
    where now f(r) refers to the EGB corrected form (<ref>).
    The curvature term in the Raychaudhuri equation is then
    
    R_MN n^M n^N = R_tt(n^t)^2 + R_rr(n^r)^2 = 0  ,
    
    which again vanishes identically, as we found previously in the perturbative case.
    It is important to note that the vanishing result is independent of the precise functional form
    of f(r), and only requires the generic structure of the metric
    (<ref>).
    
    The geodesic timelike null vector for the metric (<ref>)
    takes the form
    
    u^M = (1/f(r),±√(1-f(r)),0⃗)  ,
    
    where we are assuming f(r) < 1, which is the case for the Einstein branch of solutions.
    The curvature term in the Raychaudhuri equation then takes the form
    
    R_MN u^M u^N = R_tt(u^t)^2 + R_rr(u^r)^2 =1/2(f”(r) + (D-2)/r f'(r))  .
    
    The general expression for f(r) given in (<ref>) is not
    particularly illuminating; however, expanding the solution
    in powers of λ, we obtain
    
    R_MN u^M u^N ≈λ (D-4)(D-3)(D-2)(D-1) μ^2/4r^2D-2 + 𝒪(λ^2)  ,
    
    matching
    our perturbative result (<ref>), as expected.
    
    As in the perturbative analysis, 
    we see that the EGB corrections to the
    Schwarzschild black hole background give rise to
    either vanishing or convergent contributions
    to the null and timelike Raychaudhuri equations.
    Thus, the EGB corrections themselves, despite being quadratic in curvature, do not alleviate convergence
    that leads to conjugate points and singularities in these
    spaces.
    
    We have restricted ourself to spherically symmetric black hole solutions
    for simplicity; however, since we have not found an improvement in the convergence behavior, we do not expect that deviations from spherical symmetry are likely to produce qualitatively different results.
    
    Black hole backgrounds are not the only spacetimes of interest for the study of singularities and the Raychaudhuri equation. In the next subsection, we will
    explore the application of EGB gravity to cosmological
    spacetimes.
    
    
    
     §.§ Cosmology
    
    
    Let's examine the EGB corrections (<ref>) for a D-dimensional
    cosmological background
    [In the Appendix, we consider cosmological models in which spacetime is divided into a d-dimensional external spacetime and m-dimensional
    internal space, each of which with their own time-dependent scale factors.]
    
        ds^2 = -dt^2 + a(t)^2 (dr^2+ r^2 d̂Ω̂_D-2^2)  ,
    
    where d̂Ω̂_D-2^2 = ĝ_ijdθ^i dθ^j is the metric of a (D-2)-dimensional sphere.
    Our matter will consist of a D-dimensional perfect fluid
    
    T_MN = (ρ+p)u_M u_N + p g_MN  ,
    
    where ρ,p are the energy density and pressure of the fluid, respectively.
    
    It is straightforward to compute the Gauss-Bonnet scalar
    
    R_GB^2    =    R_ABCDR^ABCD - 4R_MN R^MN + (R_D)^2 
       =    4ä/aȧ^2/a^2(D-3)(D-2)(D-1) + ȧ^4/a^4(D-4)(D-3)(D-2)(D-1)  ,
    
    and the correction terms (<ref>)
    
    H_tt   =    2 ä/aȧ^2/a^2 (D-4)(D-3)(D-1)-(D-4)(D-3)(D-1) ȧ^4/a^4  ; 
    
    H_rr   =     -2 äȧ^2/a (D-4)(D-3) - ȧ^4/a^2 (D-4)(D-3)^2  .
    
    Notice that these corrections vanish identically for D ≤ 4, as expected since for lower dimensions the EGB terms are topological and don't contributed
    to the equations of motion.
    
    From the radial affine null tangent vector
    
    n^M = (1/a(t),±1/a(t)^2,0⃗)  ,
    
    where again ± refers to radially outgoing/ingoing rays,
    we can compute the corrections to the null Raychaudhuri equation
    
    R_MNn^M n^N    =   κ_D^2 T_MNn^M n^N + λ/2 H_MNn^M n^N 
       =   κ_D^2(ρ+p)/a^2+ λäȧ^2/a^5 (D-4)(D-3)(D-2)(D-1) - λȧ^4/a^6(D-4)(D-3)(D-2) 
       =   κ_D^2(ρ+p)/a^2+λḢ_D H_D^2/a^2 (D-4)(D-3)(D-2)(D-1) + λH_D^4/a^2 (D-4)(D-3)(D-2)^2  ,
    
    where we wrote H_D ≡ȧ/a to simplify our result.
    
    Clearly the first term in (<ref>) is positive
    for matter that obeys the null energy condition, as usual.
    However, the second term can be negative if the time variation
    of the Hubble parameter Ḣ_D is large enough.
    
    In particular, let us consider perturbative solutions for the metric
    (<ref>) of the form
    
    a(t) ≈ a_0(t) + λ a_1(t) + ...
    
    where a_0(t) is the zeroth-order solution to the λ= 0 Einstein equation
    ȧ_0^2/a_0^2   =   κ_D^2/6ρ≡(H_D^(0))^2  ; 
    ä_0/a_0   =    -κ_D^2/12 (ρ+3p) ≡Ḣ_D^(0) + (H_D^(0))^2  ,
    
    and H_D^(0) is the corresponding zeroth-order Hubble parameter.
    Inserting H_D^(0) into (<ref>),
    we obtain the expression
    
    R_MN n^M n^N = κ_D^2 (ρ+p)/a_0^2 -λκ_D^4 (D-4)(D-3)(D-2)/24 a_0^2ρ(D+1/3ρ+D-1/p) + 𝒪(λ^2)  .
    
    Assuming an equation of state p = wρ, the second term is
    negative for equations of state w > -D+1/3(D-1).
    The lower bound is always greater than -1/2, so most ordinary matter
    will satisfy this condition.
    In particular, a D-dimensional universe dominated by radiation w=1/3
    can contribute a divergent term in the null Raychaudhuri equation.
    
    Finally, consider a comoving, proper-time parameterized timelike geodesic described by the tangent vector
    
    u^M = (1,0,0,0)  .
    
    The corrections to the timelike Raychaudhuri equation due to the Gauss-Bonnet corrections come from the last term
    
    R_MN u^M u^N    =    κ_D^2 (T_MNu^M u^N+1/2 T^M_M) + λ/2 H_MN u^M u^N 
       =   κ_D^2 (D-3)ρ+(D-1)p/D-2+λ/2 (D-4)(D-3)(D-1) [2Ḣ_D + H_D^2]H_D^2  ,
    
    where again we substituted H_D ≡ȧ/a to simplify our result.
    The first term is the typical term for matter, and
    is positive for matter that obeys the strong energy condition.
    The last term, proportional to λ, is the new contribution;
    we see that for Ḣ_D ∼ -H_D^2 (as is the case for any a(t) ∼ t^n time-dependence), 
    this term can be negative, giving rise
    to a positive divergent term in the timelike Raychaudhuri equation,
    potentially opposing convergence.
    
    We have seen in both the null and timelike Raychaudhuri equations
    the existence of terms that can give rise to a positive (divergent)
    contribution to the divergence.
    This does not necessarily mean that cosmology with EGB gravity can evade the singularity problem, 
    just that the usual singularity theorems do
    not apply in a straightforward way to these backgrounds.
    
    In particular, consider the EGB-corrected equations of motion
    for the metric
    (<ref>) from (<ref>).
    Examining the (tt) component of (<ref>) we have
    1/2 (D-2)(D-1) ȧ^2/a^2 = κ_D^2 ρ - λ/4 (D-4)(D-3)(D-2)(D-1) ȧ^4/a^4  .
    
    For λκ_D^2 ρ≪ 1, solutions to (<ref>)
    take the usual Einstein form; 
    however, for λκ_D^2 ρ≫ 1, as
    we would expect to occur in the early universe near a big bang singularity,
    we have instead
    ȧ^2/a^2∼√(κ_D^2ρ/λ) .
    
    It is clear that a curvature singularity still exists where a(t)→ 0,
    R_D →∞, despite the presence of the higher
    curvature terms. 
    
    
    
    
    § BRANEWORLD GRAVITY
    
    We will consider a 5-dimensional bulk spacetime (with 5-dimensional cosmological constant Λ_5) 
    with coordinates (t,x⃗,y) and metric g_MN given by <cit.>
    ds^2 = a(t,y)^2 b(t,y)^2 (-dt^2 + dy^2) + a^2(t,y)(dr^2 + r^2 dΩ^2)  ,
    
    which is sufficiently general to capture the backreaction of the singular brane on the bulk
    as well as 4-dimensional homogeneous and isotropic cosmological evolution.
    The extra dimension has the range y ∈ (-∞,∞); however, we will additionally impose a ℤ_2 symmetry <cit.>y → -y, so that the covering space is reduced to y ∈ [0,∞), as in the RS I scenario <cit.>. 
    We embed a (singular) 3-brane at y=0 with unit normal n^M = (0,0,0,0, 1/ab)
    and induced metric q_μν = g_μν - n_μ n_ν; we will use capital Latin letters M,N for 5-dimensional
    coordinate indicies and lowercase Greek letters μ,ν for 4-dimensional brane indicies.
    The 5-dimensional energy-momentum tensor
    has the form
    
    ^(5)T_MN = -Λ_5 g_MN + δ^μ_M δ^ν_N S_μνδ(y),  S_μν = -σ q_μν  + τ_μν
    where σ is the brane tension
    and τ_μν is the brane matter energy momentum tensor.
    Note that the usual Minkowski RS I model has the solution
    
    ds^2 = e^-2K(z)η_μν dx^μ dx^ν + dz^2  .
    
    This can be obtained from (<ref>) under the limits and identifications
    
    b(t,y)    →    1  ; 
     a(t,y)    →    a(y) = e^-K(z)  ; 
    
    z    =   ∫_0^∞ a(y) dy  . 
    
    
    Gravity on the 3-brane hypersurface is induced by its embedding in the extra dimensions; in particular,
    the induced 4-dimensional Einstein equations on the brane are <cit.>
    ^(4)R_μν - 1/2 q_μν ^(4)R = - Λ_4 q_μν + κ_4^2 τ_μν + κ_5^4 π_μν - E_μν,
    
    where Λ_4= 1/2κ_5^2  ( Λ_5 + 1/6κ_5^2  σ^2 ) and κ_4^2= κ_5^4 σ/6,
    and the new tensors π_μν, E_μν are defined as
    π_μν   =   lim_y→ 0[-1/4τ_ματ^α_ν +1/12ττ_μν +1/8 q_μντ_αβτ^αβ-1/24 q_μντ^2]  ; 
    
    E_μν   =   lim_y→ 0[ ^(5)C_μανβ n^α n^β] ,
    
    where ^(5)C_MNPQ is the 5-dimensional Weyl tensor; see Appendix <ref> for more details.
    The first two terms on the right-hand side of (<ref>) are just the usual 4-dimensional cosmological constant
    and energy-momentum tensor sources, while the last two terms arise from the extra dimensional embedding of the brane.
    In order to put this into a form suitable for use with the Raychaudhuri equation we trace-reverse (<ref>)
    
    ^(4)R_μν   =   Λ_4 q_μν + 8π G_N (τ_μν - 1/2q_μντ^μ_μ)
    	+ κ_5^4 (π_μν - 1/2 q_μνπ^μ_μ) - (E_μν - 1/2q_μν E^μ_μ) 
       =   Λ_4 q_μν + 8π G_N (τ_μν - 1/2q_μντ^μ_μ) + H_μν  ,
    
    where we rewrote the last two terms as a correction term H_μν so that the corrections take the form (<ref>) 
    as outlined in the Introduction.
    The timelike and null 4-dimensional Raychaudhuri equations thus have the new terms arising from the braneworld
    dθ/dτ   =    - ^(4)R_μν u^μ u^ν +...=  - κ_5^4 (π_μν u^μ u^ν+1/2π^μ_μ) + (E_μν u^μ u^ν+1/2 E^μ_μ) + ...  ;
    dθ̂/dλ   =    - ^(4)R_μν n^μ n^ν + ...= - κ_5^4 π_μν n^μ n^ν + E_μν n^μ n^ν+ ...  ,
    
    where ... on the right-hand sides denote non-positive terms.
    
    Assuming the brane energy-momentum tensor is of the perfect fluid form, with energy density ρ(t) and pressure p(t),
    we can immediately compute the π_μν tensor contribution for the metric (<ref>)
    
        [                                              π_tt = 1/12ρ^2(t) a^2(t,0) b^2(t,0); π_rr = 1/12ρ(t) (2p(t)+ρ(t)) a^2(t,0); π_θθ = 1/12ρ(t) (2p(t)+ρ(t)) a(t,0)^2 r^2;                              π_ϕϕ = 1/12ρ(t) (2p(t)+ρ(t)) ;a^2(t,0) r^2 sin^2θ  , ]
    
    where we evaluated the terms at y = 0, while we compute E_μν directly from the  metric (<ref>)
    
    E_tt    =    lim_y→ 01/2[b”/b - b'^2/b^2 - b̈/b + ḃ^2/b^2] ≡lim_y→ 01/2β(t,y); 
    
    E_rr    =    1/6lim_y→ 0β(t,y) = E_θθ/r^2 = E_ϕϕ/r^2 sin^2θ  ,
    
    where a dot denotes a derivative with respect to t, and a prime denotes a derivative with respect to y.
    All of the components of the E_μν tensor are proportional to the quantity β(t,y), involving a second derivative b”(t,y).
    
    In order to evaluate the common quantity β(t,y) appearing in E_μν in the limit y→ 0, we need to use 
    the 5-dimensional Einstein equations ^(5)R_MN - 1/2 g_MN ^(5)R = κ_5^2 ^(5) T_MN for the
    metric (<ref>) and sources (<ref>), which take the form <cit.>(tt)        3 [2ȧ^2/a^2 + ȧḃ/ab - a”/a + a'b'/ab] = a^2 b^2 κ_5^2[Λ_5 + (ρ+σ) δ(y)]  ; 
    (yy)        3 [ä/a - ȧḃ/ab - 2a'^2/a^2 - a'b'/ab] = a^2 b^2κ_5^2 Λ_5  ; 
    (ty)        3 [-ȧ'/a + 2 ȧ a'/a^2 + ȧ b'/ab + a'ḃ/ab] = 0  ; 
    (rr) & (ij)       [3ä/a + b̈/b - ḃ^2/b^2 - 3 a”/a - b”/b + b'^2/b^2] = a^2 b^2κ_5^2[Λ_5 + (σ- p) δ(y)]  .
    
    Integrating the (tt) and (ij) equations across the brane gives rise to
    .a'/a|_y=0 = 1/6 a(t,0) b(t,0) κ_5^2(ρ+σ), . b'/b|_y=0 = -1/2 a(t,0) b(t,0) κ_5^2(ρ+p)  .
    
    In addition, the restriction of the (ty) component of the Einstein equations to the brane at y=0
    takes the form of the usual energy-momentum conservation equation on the brane
    ρ̇+ 3 .ȧ/a|_y=0 (ρ+p) = 0  .
    
    Finally, restricting the (yy) component of the Einstein equations to the brane at y=0 
    we can obtain a
    “Hubble”-like equation
    (1/abȧ/a)^2 = κ_5^4 (ρ+σ)^2/36 + κ_5^2Λ_5/6 + μ/a^4  ,
    
    Where the last term arises as an integration constant and is known as the “dark radiation” term <cit.>.
    We will assume μ≥ 0.
    
    We can rearrange (<ref>) to write
    β(t,y) = b”/b - b'^2/b^2 - b̈/b + ḃ^2/b^2 = 3 (ä/a - a”/a) - a^2 b^2 κ_5^2(Λ_5 -  p δ(y))  .
    
    Subtracting (<ref>) and (<ref>), this then becomes
    β(t,y)= 6 (ȧ^2/a^2 - a'^2/a^2) - a^2 b^2 κ_5^2(Λ_5 + (ρ+p)δ(y)) .
    
    Using the jump condition (<ref>) and the Hubble equation (<ref>), this simplifies considerably to
    β(t,y) = 6 μb^2/a^2- a^2 b^2 κ_5^2 (ρ+p) δ(y)  .
    
    The E_μν tensor is evaluated in the limit y→ 0; since the delta-function does not contribute in this 
    limit, we have
    
    E_μν∝lim_y→ 0β(t,y) = 6μb^2/a^2  ,
    
    thus the Weyl tensor term E_μν only potentially contributes to the induced Einstein equations on the brane through
    the dark radiation term. However, the induced
    Einstein equations (<ref>) are not closed,
    and must be supplemented by
    additional constraints <cit.>.
    Examining these constraints for our background
    and matter content in Appendix <ref>, we find that the dark radiation term
    μ must vanish, so that the E_μν contributions
    vanish as well.
    
    In order to evaluate the correction terms (<ref>,<ref>), we will take the
    proper-time parameterized timelike vector and affine null vectors on the brane to be
    
        u^μ   =   ( 1/a(t,0)  b(t,0), 0,0,0)  , 
        
        n^μ   =   (1/a(t,0)  b(t,0), ±1/a(t,0), 0, 0  )  ,
    
    where the ± refers to radially outgoing/ingoing null rays.
    We can then compute the additional braneworld terms
    
    H_μν u^μ u^ν   =   κ_5^4 ρ(ρ + 3/2 p); 
    
    H_μν n^μ n^ν   =   κ_5^4 ρ^2/12  .
    
    We see that the braneworld terms only contribute positively to the additional terms for null vectors, while
    the additional contributions for timelike vectors can be negative if p < -2/3 ρ.
    Including these results, we have the following expressions for the timelike and null Raychaudhuri equations on the brane
    
        dθ/dτ   =   Λ_4-1/2 G_N(ρ+3p) - κ_5^4ρ[ρ+3/2p]  + ...   ; 
        dθ̂/dλ   =    -G_N(ρ+p)-κ_5^4ρ^2/12 + ...   ,
    
    where again ... refer to strictly non-positive terms.
    Thus, we see that the corrections to the induced Einstein equations on the brane due to the braneworld embedding of the brane only contribute
    to convergence for the null Raychaudhuri equation. However, the additional contributions do add divergent (positive) terms to the right-hand
    side of the timelike Raychaudhuri equation; we note, though, that the second term in Eq.(<ref>), which is present in un-modified
    4-dimensional General Relativity, is also divergent in the presence of matter that violates the strong energy condition, p < -1/3 ρ.
    
    
    
    Is the presence of the additional positive term in Eq.(<ref>) arising from induced gravity on the brane sufficient to avoid a singularity?
    Rewriting the Friedmann-like equation Eq.(<ref>) using the relations
    Λ_4= 1/2κ_5^2  ( Λ_5 + 1/6κ_5^2  σ^2 ) and κ_4^2= κ_5^4 σ/6
    and by reparameterizing the time variable on the brane into cosmic time, dt' = a(t,0) b(t,0) dt, we have <cit.>(da/dt'/a)^2 = κ_4^2/3ρ( 1+ ρ/2 σ) + Λ_4/3  . 
    
    Notice that at late times, when the energy density is much smaller than the brane tension ρ≪σ, we recover
    the usual Friedmann equation of General Relativity.
    However, at early times ρ≥σ, assuming an equation of state ρ = wp (w≠ 1) implies ρ∝ a^-3(1+w) (recall that local energy
    conservation on the brane Eq.(<ref>) is unchanged)
    so that <cit.>:
    da/dt'/a∼a^-3(1+w)/√(σ)  ; ⇒ a ∼((1+w)t'/√(σ))^1/3(1+w)  ,
    
    
    
    
    which diverges as t'⟶ 0.
    Thus, we expect the Big Bang singularity persists in braneworld induced gravity models as well.
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    § F(R) THEORIES
    
    
    Higher curvature gravity theories, collectively known as f(R) theories, have been studied both for its improved short distance properties,
    although at the expense of introducing ghosts <cit.>, 
    as well as for their potential for avoiding singularities, and as 
    a model of inflation and an alternative model for dark energy. We start with the action for f(R) gravity in the so-called Jordan frame <cit.>
    
        S = 1/κ_4^2∫ d^4 x √(-g) f(R) + S (g_μν, ψ) ,
    
    where the last term represents
    matter action.
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    From this action, variation of the metric gives rise to the
    modified Einstein equations <cit.>
    R_μν f'(R) = κ_4^2 T_μν+1/2 g_μν f(R) + (∇_μ∇_ν - g_μν) f'(R)  .
    
    From this we can write the null and timelike Raychaudhuri equations for f(R) gravity:
    
        dθ̂/dλ   =    -θ̂^2/2 - [κ_4^2/f'(R) T_μν + 1/2g_μνf(R)/f'(R) + 1/f'(R)(∇_μ∇_ν - g_μν)f'(R)] n^μ n^ν  ; 
        dθ/dτ   =    -θ^2/3 - [κ_4^2/f'(R) T_μν + 1/2g_μνf(R)/f'(R) + 1/f'(R)(∇_μ∇_ν - g_μν)f'(R)] u^μ u^ν  ,
    
    in terms of the expansions θ̂, θ 
    for null and timelike geodesics n^μ(λ),u^μ(τ) parameterized
    by affine and proper times λ,τ, respectively.
    
    Since we are interested in singularities with larger curvatures, we
    will consider an f(R) theory consisting of the usual scalar curvature
    plus a perturbative
    correction of higher power in the curvature 
     f(R)=R [ 1 + (ℓ_Pl^2 R)^n-1], where we will take the Planck
     length ℓ_Pl to control the strength of the corrections, and
    n is a positive integer. 
    Since we are interested in the large curvature regime,
    i.e. R≫ 1/ℓ_Pl^2, we will consider the case where
    the second term in f(R) dominates so that f ∼ R^n.
    In this case, the equations of motion are dominated by the higher curvature
    correction terms, which act as an effective energy-momentum tensor.
    
    
    
    Substituting this into the Raychaudhuri equations (<ref>,<ref>),
    we obtain
    
        dθ̂/dλ = -θ̂^2/2 - 
        
        [
         κ_4^2 T_μν/2 n ℓ_Pl^2n-2 R^n-1 
        + R/2n g_μν
        + 1/n  R^n-1( ∇_μ∇_ν  - g_μν) R^n-1]n^μ n^ν .;
        
              dθ/dτ = -θ^2/3 - 
        
        [
         κ_4^2 T_μν/2 n ℓ_Pl^2n-2 R^n-1 
        + R/2n g_μν
        + 1/n  R^n-1( ∇_μ∇_ν  - g_μν) R^n-1]u^μ u^ν .
    
    Note that (i) the sign of the second term  in square brackets
    in each equation 
    is opposite to the T_μν term, and (ii) the sign of the last term 
    in each equation depends on the spacetime under consideration. 
    
    Next, we consider homogeneous and isotropic cosmological backgrounds with the metric
    
    
        ds^2 = -dt^2 + a(t)^2 [ dr^2/1-kr^2 + r^2 ( dθ^2 + sin^2θ dϕ^2 ) ]  ,
    
    and corresponding scalar curvature
    
        R = 6( ä/a  + ȧ^2/a^2 + k/a^2)   .
    
    The null Raychaudhuri equation (<ref>), using n^μ n^ν g_μν = 0, becomes
    
        dθ̂/dλ = -θ̂^2/2 - κ_4^2 T_μνn^μ n^ν/2n ℓ_Pl^2n-2 R^n-1+ 1/n R^n-1(n^t)^2 ∂_t^2 R^n-1  ,
    
    while the timelike Raychaudhuri equation (<ref>) becomes
    
        dθ/dτ =
        -θ^2/3 - 
        [
        κ_4^2 T_μν/2 n ℓ_Pl^2n-2 R^n-1 
        + 3/n( ä/a  + ȧ^2/a^2 + k/a^2)
         g_μν]u^μ u^ν  .
    
    When the cosmological evolution is dominated by the curvature corrections,
    we have an effective perfect fluid
    
    
        p_ eff= w_ effρ_ eff,  w_ eff = - 6n^2-7n-1/6n^2-9n+3 ,
    
    with scale factor (w_ eff≠ 1):
    
        a = a_0 t^2/3(1+w_ eff),
        ȧ/a = 2/3(1+w_ eff)t,
         ä/a= - 2(1+3w_ eff)/9(1+w_ eff)^2t^2, ä/a  + ȧ^2/a^2 + k/a^2  
        =2(1-3w_ eff)/9(1+w_ eff)^2t^2 + k/a_0^2 t^4/3(1+w_ eff) .
    
    We can substitute this solution for the scale factor into the Raychaudhuri equations to obtain
    
        dθ̂/dλ       = -θ̂^2/2 + 1/n {2(1-3w_ eff)/3n(1+w_ eff)^2t^2 + 3k/n a_0^2 t^4/3(1+w_ eff)}^n-1×
              [-κ_4^2 T_μν n^μ n^ν/2ℓ_Pl^2n-2 + (n^t)^2 ∂_t^2 {2(1-3w_ eff)/3n(1+w_ eff)^2t^2 + 3k/n a_0^2 t^4/3(1+w_ eff)}^n-1] ;  
        dθ/dτ      =
        -θ^2/3 - 
        κ_4^2 ρ/2n ℓ_Pl^2n-2{2(1-3w_ eff)/3n(1+w_ eff)^2t^2 + 3k/n a_0^2 t^4/3(1+w_ eff)}^n-1 
        +  {2(1-3w_ eff)/3n(1+w_ eff)^2t^2 + 3k/n a_0^2 t^4/3(1+w_ eff)} ,
    
    where we considered the rest frame of the fluid
    streamlines, for which T_μν u^μ u^ν =ρ, and used
    u^2=-1.
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    Since we are primarily interested in the corrected Raychaudhuri equation close to the putative Big-Bang
    singularity, i.e. pre-inflation with w≠ -1, we ignore the energy-momentum
    term in Eqs.(<ref>,<ref>), since this term becomes sub dominant at
    larger curvatures.
    Furthermore, omitting exotic matter from our discussions, i.e. w≰ -1/3 such that the k term can also
    be ignored, we arrive at the rather simple form for the Raychaudhuri equations
    
        dθ̂/dλ       = -θ̂^2/2 + 2(n-1)(2n-1) (n^t)^2/n  t^2
        dθ/dτ       = -θ^2/3 + 2(1-3w_ eff)/3n(1+w_ eff)^2t^2 
        = -θ^2/3 + 3(4n^2-5)(2n^2-3n+1)/(n-2)^2 t^2
    
    As can be seen, the second term in both Eqs.(<ref>,<ref>) are repulsive and dominate over the first term 
    (due to the factor of 1/t^2 therein) for all n.
    This causes geodesics to not to converge in this epoch, potentially preventing the formation of a singularity. 
    At later epochs, e.g. during inflation and thereafter, that term is sub-dominant as are additional terms 
    originating in an actual perfect fluid described by the equation p=wρ, with -1 ≤ w ≤ 1/3.
    Therefore the standard conclusions from the Raychaudhuri equation 
    and the singularity theorems hold. 
    Note that while in <cit.>,
    the authors consider constant curvature cosmological
    metrics (Einstein spaces), in <cit.> the authors consider theories of the form
    f(R) = R + α/R^n with α <0 and n ∈ℝ. These are relevant to late 
    time acceleration, as an alternative to dark energy. As mentioned above, here we are
    primarily concerned with very early times, to examine effects of repulsive terms (if any) near the 
    initial singularity. Furthermore, similar to the above references we too assume that the standard energy conditions are valid for the cosmological fluids. 
    
    
    
    
    
    § LOOP QUANTUM COSMOLOGY
    
    
    In loop quantum cosmology, one starts with the standard large scale homogeneous and isotropic model of the Universe, 
    described by metric (<ref>), and obtains an effective Hamiltonian, incorporating the discrete quantum nature
    of spacetime at the fundamental level, quantum back-reaction and the behavior of the scale factor at very small length scales
    <cit.>
    
        H_ eff = - 6/κ_4^2 γ^2sin^2 (λβ)/λ^2 V +  H_ matter ,
    
    where λ is a measure of fundamental discreteness, γ the Immirizi parameter, V the volume, 
    β (=γ H) its conjugate, and H_ matter the matter Lagrangian.
    Next, using the Hamiltonian constraint H_ eff≈ 0 and the Hamilton's equation for V,
    namely V̇ = {V,  H_eff}, one obtains the ` Loop Quantum Corrected Raychaudhuri Equation'
    
    
        ä/a = Ḣ + H^2 = 
         - κ_4^2/4[  
        ρ( 1 - ρ/ρ_c)
        + 3 { P (1-2ρ/ρ_c) - ρ^2/ρ_c}]  ,
    
    where: ρ_c ≡3/8πγ^2 λ^2 V^2.
    
    
    By substituting <cit.>
    
        ȧ/a   =   θ/3 ,  
        θ̂   =   2/a^2[ ȧ - 1/r]
    
    for timelike and null geodesics, respectively, we can
    obtain the corresponding Raychaudhuri equations
    
    
        dθ/dτ   = - θ^2/3 
         - κ_4^2/4[  
        ρ( 1 - ρ/ρ_c)
        + 3 { P (1-2ρ/ρ_c) - ρ^2/ρ_c}]     , 
        
        dθ̂/dλ   = 
        - θ̂^2/2 + 2/a^2Ḣ   
        
                = θ̂^2 - 2/(a^2 r)^2 - 2θ/a^2 r   
        - κ_4^2/2a^2[ ρ( 1 - ρ/ρ_c)
        + 3 { P (1-2ρ/ρ_c) - ρ^2/ρ_c}] 
        
        
        
        
        
        
            .
    
    
    Note that there are three positive (repulsive) terms in the RHS of the RE (<ref>), (<ref>),
    effective near ρ≃ρ_c, i.e. in the early Universe. 
    Since the above equations were derived using a  non-perturbative quantization scheme, without assuming a fixed background metric
    and allowing for back-reaction, one expects that they will hold good all the way to t→ 0. At this point, the repulsive terms
    take over (since ρ→∞) and the singularity theorems and their conclusions cease to hold. 
    
    
    
    
    § CONCLUSIONS
    
    In this paper, we have examined the effects of several different alternative theories of gravity on the geodesic convergence properties of the Raychaudhuri equation
    in some simple backgrounds, and found some cases where these corrections provide repulsive terms.
    
    In the case of string theory, we studied the leading order corrections to the Einstein Equations (for constant dilaton and vanishing form fields)
    arising in the form of Einstein-Gauss-Bonnet gravity. We did not find repulsive contributions to the Raychaudhuri equation for
    D-dimensional black hole backgrounds, which is consistent with known exact black hole solutions in Einstein-Gauss-Bonnet gravity, which still contain
    singularities. We did find repulsive terms for D-dimensional cosmological backgrounds, however these terms do not apear to be significant enough
    to prevent the existence of a Big Bang singularity, at least at leading order.
    It is therefore interesting that the string theory corrections to the Raychaudhuri equation do not have a definite sign, but instead the sign of the corrections
    depends on the background under consideration.
    We note that our results here should be taken as preliminary steps towards analyzing the potential for string theory corrections to resolve
    singularities. In particular, we have only examined the leading-order pure gravity corrections, setting the dilaton to a constant and form fields
    to zero. Allowing these fields to be dynamical will give rise to additional terms in the Raychaudhuri equation that could act repulsively.
    Additionally, the leading corrections we considered here vanish for type IIA/IIB string theory; it would be interesting,
    though technically challenging, to include the next order of corrections.
    
    We also considered corrections to the Einstein equations coming from induced gravity on a brane embedded in a warped 5-dimensional bulk.
    In this case, the corrections show up as quadratic in the brane-localized energy-momentum tensor. For a cosmological background we found that
    the corrections to the null Raychaudhuri equation always increase convergence, while the corrections ot the timelike Raychaudhuri equation
    can give rise to repulsive terms when the brane matter is described by a perfect fluid with pressure p < -2/3 ρ.
    Nonetheless, it appears that a Big Bang singularity persists in this case as well. It would be interesting
    to study these corrections for black hole backgrounds, as well as explore the role that bulk fields play in providing additional
    constraints.
    
    We examined corrections for curvature-dominated 
    f(R) theories of the form f(R) = R [1+(ℓ_Pl R)^n-1], where ℓ_Pl is the Planck length and n is a positive integer.
    We found that the f(R) corrections can contribute repulsive terms to the timelike and null Raychaudhuri equations for w_ eff > -1/3, where w_ eff(n) = -(6n^2-7n-1)/(6n^2-9n+3) is the effective equation of state of the curvature corrections.
    It would be interesting to study these corrections further to examine whether they do indeed lead to a resolution of the cosmological singularity.
    
    Finally, we also considered corrections to the Friedmann equation for Loop Quantum Cosmology.
    These corrections give rise to repulsive terms in both the timelike and null Raychaudhuri equations, suggesting a resolution of the Big Bang singularity.
    
    
    Note that in this case the results are non-perturbative in nature. 
    It is important to explore f(R) theories and Loop Quantum Cosmology further to see if indeed all possible geodesics are
    complete in these theories, and if so, whether there exists another criterion for the existence of singularities. If not, the
    corresponding spacetimes should play an important role near the normally singular regions inside black holes and in cosmology. 
    
    
    In all of our examples we have considered simple isotropic and/or homogeneous backgrounds, and it would be interesting to study how robust our results are
    to anisotropic or inhomogeneous deviations.
    
    Since the power of the singularity theorems is in their general applicability, and generic situations likely would not possess such a high degree of symmetry, it is not clear whether the additional terms to the Raychaudhuri equation would continue to be repulsive in more general backgrounds with less symmetry. We leave this general analysis for future work.
    
    In addition, all of our analysis has focused on the corrections to the classical equations of motion from alternative theories of gravity.
    It would be interesting to examine quantum effects
    (which depend on the wavefunction of the fluid) over and above the above classical terms. However, since
    these effects are always repulsive in nature, and 
    prevents the quantal trajectories (quantum counterparts of geodesics) from crossing, such effects should reinforce the above conclusion <cit.>.
    
    
    
    
     Acknowledgment
    This work is supported by the Natural Sciences and Engineering
    Research Council of Canada. SSH would like to thank Perimeter Institute, where part of the work was done. DJB is supported by a PHD fellowship from the South African National Institute for Theoretical Physics. NM is supported by the South African Research Chairs Initiative
    of the Department of Science and Technology and the National Research Foundation of South Africa. Any opinion, finding and conclusion or recommendation expressed in this material is that of the authors and the NRF does not accept any liability in this regard. 
    
    
    
    § RAYCHAUDHURI EQUATIONS IN D DIMENSIONS
    
    
     §.§ Timelike
    
    Consider a congruence of timelike geodesics parameterized
    by proper time τ with tangent vectors u^M(τ) satisfying
    u^M u_M = -1 and u^M ∇_M u^N = 0
    for some D-dimensional metric g_MN.
    We can define the “spatial” (or transverse) metric as
    
    h_MN≡ g_MN + u_M u_N  ,
    
    which is transverse to the tangent vectors u^M h_MN = 0.
    We will define the following quantities:
    θ≡∇_M u^M       
    σ_MN = 1/2(∇_M u_N + ∇_N u_M) - 1/D-1 h_MNθ      
    ω_MN = 1/2(∇_M u_N - ∇_N u_M)       
    Note that σ_MN, ω_MN are purely spatial (or transverse), since
    u^M σ_MN = 0 = u^M ω_MN.
    In addition, ω_MN = 0 for tangent vectors that are hypersurface orthogonal, as will
    be the case for all of the tangent vectors considered in this paper.
    Following <cit.>, the expansion obeys the Raychaudhuri equation
    dθ/dτ = -θ^2/D-1 - σ_MNσ^MN + ω_MNω^MN - R_MN u^M u^N  .
    
    Since σ_MN is purely spatial, the second term is non-positive, and as mentioned above we will
    be considering cases where ω_MN = 0, so that the right-hand side of Eq.(<ref>)
    is non-positive as long as R_MN u^M u^N ≥ 0.
    
    
    
     §.§ Null
    
    The derivation of the null Raychaudhuri equation proceeds in a similar way as the timelike case above,
    with an additional complication in identifying the transverse directions to the null geodesic.
    
    Consider a congruence of null geodesics parameterized by affine parameter λ with
    tangent vectors n^M(λ), satisfying n^M n_M = 0 and n^M ∇_M n^N = 0.
    In order to define the transverse directions to the geodesic we need an “auxiliary"
    null vector k^M such that[k^M need not be a geodesic, and indeed in most cases it is not possible 
    for both k^M and n^M to be geodesics and cross-normalized to a constant.]k^M k_M = 0 and n^M k_M = -1.
    We define the transverse metric to the null geodesics
    ĥ_MN = g_MN + n_M k_N + k_M n_M
    
    which has dimension g^MNĥ_MN = D-2 and is transverse to both n^M and k^M:
    n^M ĥ_MN = 0 = k^M ĥ_MN.
    We then proceed as in the timelike case by defining:
    θ̂≡∇_M n^M       
    σ̂_MN≡1/2(∇_M n_N + ∇_N n_M )- 1/D-2ĥ_MNθ̂      
    ω̂_MN≡1/2(∇_M n_N - ∇_N n_M)       
    As before, σ_MN, ω_MN are transverse to the null geodesics.
    The expansion obeys the Raychaudhuri equation
    dθ̂/dλ = -θ̂^2/D-2 - σ̂_MNσ̂^MN + ω̂_MNω̂^MN - R_MN n^M n^N  .
    
    As in the timelike case, ω_MN = 0 for geodesics that are hypersurface orthogonal; further, since σ_MN is purely
    spatial the second term in Eq.(<ref>) is non-positive, so that the entire right-hand side is non-positive as well
    as long as R_MN n^M n^N ≥ 0.
    
    
    
    § GAUSS-BONNET CORRECTIONS IN ANISOTROPIC COSMOLOGY
    
    
    
    We will examine the correction terms H_MN for a D-dimensional background
    consisting of a d-dimensional flat, homogeneous spacetime, with coordinates {r,θ^i},
    and a m-dimensional homogeneous spacetime, with coordinates {y^m}, each with their own
    scale factor
    
        ds^2 = -dt^2 + a(t)^2 (dr^2+ r^2 d̂Ω̂_d-1^2) + b(t)^2 g̃_mn dy^m dy^n ,
    
    where D = d + m + 1.
    
    The Gauss-Bonnet scalar is
    
    R_GB^2    =    R_ABCDR^ABCD - 4R_MN R^MN + (R_D)^2 
       =    4d(d-1)(d-2)ä/aȧ^2/a^2 + d(d-1)(d-2)(d-3)ȧ^4/a^4
           + 4 m(m-1)(m-2) b̈/bḃ^2/b^2 + m(m-1)(m-2)(m-3) ḃ^4/b^4
           + 4 dm(m-1) ä/aḃ^2/b^2 + 4 dm(d-1) ȧ^2/a^2b̈/b
    		+ 2 dm(d-1)(m-1) ȧ^2/a^2ḃ^2/b^2
           +  R_m/b^2[4 d ä/a + 2 d(d-1) ȧ^2/a^2 + 4 (m-2) b̈/b^3
    		-6m(m-1) ḃ^2/b^2 - 3  R_m/b^2] .
    
    Note that R_GB^2 simplifies considerably when d=3 and m=1,2,
    corresponding to D=5 or 6-dimensional spacetime, respectively, and the internal space is flat R_m = 0
    R_GB^2|_d=3, m=1,2 = 24äȧ^2/a^3 + 24 m ȧ^2/a^2b̈/b  ,
    
    though this doesn't seem to give any particularly useful interpretation.
    
    We can now compute the correction terms
    
    H_tt   =   g_tt/2 R_GB^2 - 2 R_D R_tt + 4R_tt g^tt R_tt + 4 (g^rr)^2 R_rr R_trtr + 4 g^ii'g^jj' R_i'j' R_titj + 4 R^mn R_tmtn - 2  R_tABC R_t^ABC
       =    -1/2 d(d-1)(d-2)(d-3) ȧ^4/a^4 - 1/2 m(m-1)(m-2)(m-3) ḃ^4/b^4 - d m (d-1)(m-1) ȧ^2/a^2ḃ^2/b^2
           +  R_m/b^2[-d(d-1) ȧ^2/a^2 + 3 m(m-1) ḃ^2/b^2 + 3/2 R_m/b^2] ; 
    
    H_rr   =   g_rr/2 R_GB^2 - 2 R_D R_rr + 4 R_rr g^rr R_rr +4 R_tt (g^tt)^2 R_rtrt + 4 R_i'j'g^ii'g^jj' R_rirj- 2 R_rABCR_r^ABC
       =     2(d-1)(d-2)(d-3) äȧ^2/a + 1/2 (d-1)(d-2)(d-3)(d-4)ȧ^4/a^2
           + 2 m(m-1)(m-2) a^2 b̈ḃ^2/b^3+ m(m-1)(m-2)(m-3) a^2 ḃ^4/b^4
           + 2 m(m-1)(d-1) a äḃ^2/b^2 + 2 m(d-1)(d-2) ȧ^2 b̈/b + m(m-1)(d-1)(d-2) ȧ^2 ḃ^2/b^2
           + a^2  R_m/b^2[2dä/a + d(d-1) ȧ^2/a^2 + 2(m-2) b̈/b^3 - 3 m(m-1) ḃ^2/b^2 - 3/2 R_m/b^2] .
    
    Notice that H_tt vanishes identically for all d =3, m=1,2 and a flat internal space R_m = 0, 
    while H_rr is considerably more complex.
    
    
    
      §.§.§ Raychaudhuri Corrections
    
    
    We now consider a d-dimensional radial, null, affine tangent vector
    
    n^M = (1/a(t),±1/a(t)^2,0⃗,0⃗)  .
    
    Note that n^M n_M = 0 (it is null) and n^M satisfies the affine condition n^N ∇_N n^M = 0, even for the anisotropic background considered here.
    
    The null Raychaudhuri equation takes the form
    dθ̂/dλ = -θ̂^2/D-2-|σ̂|^2 - R_MNn^M n^N  .
    
    Our corrections due to the Gauss-Bonnet term appear on the right hand side, as
          8/α'R_MNn^M n^N ∼   H_MNn^M n^N = H_tt n^t n^t + H_rr n^r n^r 
          = (d-1)(d-2)(d-3)(d-5)/2ȧ^2/a^6 + 2 (d-1)(d-2)(d-3) äȧ^2/a^5
           + (m-1)(m-2)(m-3)(m-5)/2ḃ^4/a^2 b^4 + 2 m(m-1)(m-2) b̈ḃ^2/a^2 b^3 + 2 m (d-1)(d-2) ȧ^2/a^4b̈/b
           + 2 m(m-1)(d-1) ḃ^2/b^2[ä/a^3-ȧ^2/a^4] + 2 (m-1)  R_m b̈/a^2 b^5  .
    
    As  before, choosing d =3, m=1,2 and a flat internal space R_m = 0 simplifies the result considerably
    
    H_MNn^M n^N    =    2m(m-1)(d-1) ḃ^2/b^2[ä/a^3 -ȧ^2/a^4] + 2 m(d-1)(d-2) ȧ^2/a^4b̈/b
       =    4m(m-1) ḃ^2/b^2Ḣ_d/a^2 + 4 mȧ^2/a^4b̈/b  ,
    
    where we wrote H_d ≡ȧ/a, and thus Ḣ_d = ä/a - ȧ^2/a^2.
    The overall sign of this correction term is not clear; typically we expect Ḣ_d ≤ 0 for spacetimes that don't have
    NEC violating matter, so this term is potentially negative. However, it is unclear the scenarios for which b̈/b < 0.
    
    Finally, consider a comoving, proper-time parameterized timelike geodesic described by the tangent vector
    
    u^M = (1,0,0⃗,0⃗) .
    
    It is straightforward to check that u^M u_M = -1 and u^N ∇_N u^M = 0.
    
    The timelike Raychaudhuri equation takes the form
    dθ/dτ = -θ^2/D-1 - |σ|^2 - R_MN u^M u^N  ,
    
    where again the corrections to the Raychaudhuri equation due to the Gauss-Bonnet corrections come from the last term
          8/α' R_MN u^M u^N ∼  H_MN u^M u^N = H_tt u^t u^t 
          = -d(d-1)(d-2)(d-3)/2ȧ^4/a^4 - m(m-1)(m-2)(m-3)/2ḃ^4/b^4 -dm(d-1)(m-1) ȧ^2/a^2ḃ^2/b^2
           +  R_m/b^2[-d(d-1)ȧ^2/a^2+3m(m-1) ḃ^2/b^2 + 3/2 R_m/b^2]  .
    
    Notice here that for a flat internal space R_m = 0, these corrections are  manifestly negative, thus they cause
    divergence in the Raychaudhuri equation.
    
    More generally, we see that whether these additional correction terms due to Gauss-Bonnet act 
    to increase or oppose convergence depends on the particular background, and whether one is considering timelike or
    null rays. Thus, there are backgrounds in which the usual singularity theorems no longer hold.
    
    
    
    
    
    § WEYL CURVATURE TENSOR
    
    
    In Section <ref>, we encountered
    corrections to the induced 4-dimensional
    Einstein equations that involved the 5-dimensional
    Weyl curvature through the tensor
    
    E_μν   =    ^(5)C_μανβ n^α n^β  .
    
    
    In this Appendix, we collect some useful properties
    of the Weyl tensor, as well as evaluate this term for the braneworld metric (<ref>).
    
    
    
     §.§ Weyl Tensor
    
    We begin with some facts about the Weyl curvature tensor. First, the definition of the Weyl curvature tensor is
    
    C_αβγδ = R_αβγδ + 1/n-2[g_αδ R_γβ + g_βγ R_δα - g_αγ R_δβ - g_βδ R_γα]
    	+ 1/(n-1)(n-2)[g_αγ g_δβ - g_αδ g_γβ] R
    
    in terms of the Riemann tensor, Ricci tensor, and Ricci scalar. In effect, the Weyl tensor C_αβγδ is the traceless part of the Riemann curvature tensor.
    As such, it inherits the usual symmetry identities from the Riemann tensor, as well as a traceless condition:
          C_αβγδ = - C_αβδγ = -C_βαγδ = C_γδαβ
           C^α_βαδ = 0 
    The Weyl curvature tensor has an interesting feature: for any 2 metrics that can be related by a conformal factor
    ĝ_αβ = Ω^2 g_αβ  ,
    
    the Weyl curvature tensors are  the sameĈ_αβγδ = C_αβγδ  .
    
    This means that if a metric can be written in a conformally flat form, then C_αβγδ = 0 identically.
    
    In particular, for the RS Minkowksi metric (<ref>), this can be written in conformally flat form (with the coordinate redefinition of dχ = e^-K(z) dz).
    Further, the ansatz metric (<ref>), is  almost conformally flat, aside from the factor b(t,y); thus, we expect
    the Weyl curvature tensor for (<ref>) to be proportional to the time- and space-derivatives of b(t,y) (as we subsequently find it to be).
    
    
    
     §.§ Evaluating E_μν
    
    
    In Section <ref> we evaluated the 
    E_μν tensor 
    in the limit y→ 0, 
    finding that E_μν∝μb^2/a^2, consequently the singular terms do not contribute.
    However, we need to examine the consistency of this result.
    In particular, the induced Einstein equations (<ref>) are not closed, and must be supplemented
    by additional constraints <cit.>.
    The Bianchi identity for the induced Einstein equations (<ref>) implies (assuming that the brane EM tensor τ_μν obeys EM conservation (<ref>))
    ∇̂^μ E_μν = κ_5^2 ∇̂^μπ_μν  ,
    
    where ∇̂^μ is the covariant derivative constructed with respect to the 4d q̂_μν induced metric.
    Thus, we will examine the consistency of our limiting procedure by verifying that our result for E_μν is indeed a solution of the Bianchi identity (<ref>).
    
    
    Starting on the right-hand side, we rewrite the Bianchi identity as ∇̂^μπ_μν = ∇̂_μπ^μ_ν, and using the suitably raised versions of the π^μ_ν tensor,
    it is straightforward to compute
    ∇̂^μπ_μ t = -1/6ρρ̇- 1/2ȧ/aρ (p+ρ) = -ρ/6(ρ̇+ 3ȧ/a(ρ+p)) = 0  ,
    
    where we used energy conservation of the brane matter (<ref>) to set this entire term to zero.
    Similar results are obtained for the ν = r,θ,ϕ components of the Bianchi identity (<ref>).
    
    However, for the left-hand side of the ν = t component of the Bianchi identity, we obtain
    ∇̂_μ E^μ_t = 3ȧ/aμ/a^4 (1-b^2)  .
    
    The only way to satisfy the Bianchi identity (<ref>), then, is to require the dark radiation term to vanish, μ = 0.
    Similar results are obtained for the ν = r, θ, ϕ components.
    
    
    utphysmodb
