Generalized Darmois–Israel Junction Conditions

Abstract

We present a general method to derive the appropriate Darmois–Israel junction conditions for gravitational theories with higher-order derivative terms by integrating the bulk equations of motion across the singular hypersurface. In higher-derivative theories, the field equations can contain terms which are more singular than the Dirac delta distribution. To handle them appropriately, we formulate a regularization procedure based on representing the delta function as the limit of a sequence of classical functions. This procedure involves imposing suitable constraints on the extrinsic curvature such that the field equations are compatible with the singular source being a delta distribution. As explicit examples of our approach, we demonstrate in detail how to obtain the generalized junction conditions for quadratic gravity, $\mathcal{F}\left(R\right)$ theories, a 4D low-energy effective action in string theory, and action terms that are Euler densities. Our results are novel, and refine the accuracy of previously claimed results in $\mathcal{F}\left(R\right)$ theories and quadratic gravity. In particular, when the coupling constants of quadratic gravity are those for the Gauss–Bonnet case, our junction conditions reduce to the known ones for the latter obtained independently by boundary variation of a surface term in the action. Finally, we briefly discuss a couple of applications to thin-shell wormholes and stellar models.

1. Introduction

In general relativity, the Darmois–Israel junction conditions prescribe the appropriate boundary conditions across a singular hypersurface $(\Sigma )$ supported by a localized energy-momentum source that contains a Dirac delta distribution. Apart from metric continuity at $\Sigma$, these “junction conditions” can be expressed simply as the following equation:

$$\left[K\right]h_{ij}-\left[K_{ij}\right]=8\pi S_{ij},K\equiv K_m^m,$$

(1)

where the square bracket indicates the jump discontinuity across $\Sigma$, $h_{ij}$ is the induced metric, $S_{ij}$ is the singular Dirac delta source localized within $\Sigma$, and all indices pertain to coordinates of $\Sigma$. The Darmois–Israel junction conditions in (1) were historically obtained in [1,2,3,4] by integrating the field equations across the infinitesimal width of $\Sigma$, and were later shown to equivalently follow from the boundary variation of the Gibbons–Hawking–York term in the action [5,6].

Since its discovery, (1) has found numerous applications in a plethora of gravitational topics, being primarily used to govern the construction of spacetime geometries obtained by a “cut-and-paste” approach with $\Sigma$ as the locus of identification. These geometries include thin-shell wormholes [7] (and their other variants such as gravastars [8,9]), static and dynamical stellar models [10,11], such as the Oppenheimer–Snyder solution [12] describing gravitational collapse, and domain wall configurations [13] separating true and false vacua in cosmology.

In recent years, there have been interesting explorations of generalizing (1) to the wider context of gravitational theories beyond Einstein’s theory which we can consider to be a long-distance, effective description connected to some presumably UV-complete framework of quantum gravity such as string theory. When interpreted as such, from the effective field theory perspective, it is natural to consider adding higher-derivative terms to the Einstein–Hilbert action. Models that have attracted interest include supergravity, $\mathcal{F}\left(R\right)$-theories, Gauss–Bonnet gravity, etc. Accordingly, there has been a small number of different proposals for how (1) should be modified in several contexts. For example, junction conditions were proposed in [14,15,16] for $\mathcal{F}\left(R\right)$ gravity, and in [17,18] for quadratic gravity. However, we find these results to be incomplete as they miss contributions which can only be seen with a proper and careful treatment of terms involving products of the step function, Dirac delta function, and their derivatives. For example, consider the case that the extrinsic curvature is discontinuous across a hypersurface: $K_{ij}=K_{ij}^{-}+\left[K_{ij}\right]\Theta \left(n\right)$, where n is the normal coordinate and $\Theta$ is the step function [14,19] with $\Theta \left(0^{-}\right)=0$, $\Theta \left(0^{+}\right)=1$ and $K_{ij}^{\pm }\equiv K_{ij}\left(0^{\pm }\right)$. The Gauss–Codazzi relation $R=\widehat{R}-K^2-2{\partial }_{n}K-K_{ab}{K}^{ab}$ implies that

$$R=R^{+}\Theta +R^{-}(1-\Theta )-\left[K_{ab}\right]\left[K^{ab}\right]({\Theta }^2-\Theta ).$$

(2)

Naively, one may take

$${\Theta }^2=\Theta .$$

(3)

However, the field equations in higher-derivative gravity theory imply that generally, these terms are not isolated. Instead, they appear together with delta functions and their derivatives (this will be discussed extensively below in Section 2). As a result, (3) cannot be taken prematurely from the outset as was assumed in [14,19]. In fact, in any sensible regularization of the delta function, we have $\delta \left(n\right)={\Theta }^{\prime }\left(n\right)$ and $\int dn{\Theta }^n\left(n\right)\delta \left(n\right)=\frac{1}{n+1}.$ This implies that in the distributional sense,

$$\delta \left(n\right)\left({\Theta }^2\left(n\right)-\Theta \left(n\right)\right)\ne 0.$$

In Section 4.1.1 below, we show how the correct treatment of the $({\Theta }^2-\Theta )$-term in (2) leads to contributions of the form $\left[K_{ij}\right]\left[K_{ab}\right]\left[K^{ab}\right]$ in the junction conditions of the $R^2$-theory. Such contributions were missed in the literature [14,15,17]. As products of singular distributions are unavoidable in higher-derivative gravity and the handling of them is subtle, one of the key motivations of this work is to provide a consistent mathematical framework for the derivation of the junction conditions.

In this paper, we present a general method to derive the appropriate junction conditions for gravitational theories with higher-order derivative terms by integrating the bulk equations of motion across the infinitesimal width of the singular surface $\Sigma$, in the same spirit as was applied in the foundational paper of Israel in [2]. Fundamentally, asserting an equivalence between both sides of the field equations in the sense of distribution translates to ensuring the convergence of the integral which is now much more complicated by virtue of the higher-order differential operators. Completing the integral then yields the generalized junction conditions compatible with a Dirac $\delta$-distribution carried by the energy-momentum tensor.

We find that a combination of the following techniques furnishes an effective toolbox for solving this problem. First, we choose to work in Gaussian normal coordinates describing a local neighborhood of $\Sigma$ which is now defined by taking $n=0$, where n parametrizes the proper distance along a geodesic orthogonal to $\Sigma$. In this chart, we invoke the Gauss–Codazzi relations extensively to express various curvature terms in the equations of motion in terms of the intrinsic geometric quantities of $\Sigma$, the extrinsic curvature, and their derivatives. The problem is now reduced to examining the convergence of a one-dimensional integral. For this purpose, we find that expressing the $\delta$-function as a limit of a sequence of classical functions gives us a powerful language for organizing various integrand terms according to their algebraic order of singularity1. This allows us to identify the additional conditions on the embedding geometry that naturally ensure the absence of singular terms that are incompatible with a $\delta$-distribution source term. We call these conditions “regularity constraints”, the vanishing of some function of the extrinsic curvature components on each side of $\Sigma$. As we shall elaborate later on, it turns out that they are intimately related to the Hadamard regularization procedure [21] in the theory of distributions (see, e.g., [22]). In the absence of the regularity constraints, the generalized junction conditions are precisely the Hadamard-finite part of the integral. Our method can be readily adopted and applied to a wide variety of gravitational theories, including matter and gauge couplings, although in this paper, we focus almost exclusively on action terms that are invariants constructed from the Riemann tensor as illustrative examples.

Similar to those in ordinary GR, the generalized junction conditions specify how matter/energy should be constrained in the codimension-one hypersurface $\Sigma$ for higher-derivative gravitational theories. They enable the construction of various interesting solutions that are characterized by localized sources, some important examples being wormholes, domain walls, and gravitational collapse geometries. We included a couple of pedagogical examples in the Section 5, where we apply our results to construct a wormhole solution obtained by the “cut-and-paste” of two copies of Schwarzschild–de Sitter exteriors, and derive the boundary conditions for the energy-momentum content of a few stellar models in the context of $R^2$ gravity. We hope that these serve as useful concrete examples for further study of the generalized junction conditions in the wide landscape of gravitational theories beyond GR.

Our exposition is structured as follows. In Section 2, we present the background canvas of concepts underlying our general approach of deriving the junction conditions via a bulk integration across $\Sigma$, including the notion of delta-convergent sequences, regularity constraints, and how our method relates to Hadamard regularization. Section 3 contains a detailed derivation of the junction conditions for quadratic gravity. This is accompanied by Appendix A containing some integral identities that we developed for our purpose, and by Appendix C containing details of how the junction conditions in quadratic gravity reduce consistently to known ones for Gauss–Bonnet gravity. Section 4 presents the junction conditions for three types of theories: $\mathcal{F}\left(R\right)$ theories, a low-energy effective action descending from string theory, and higher-dimensional Euler density terms that appear in Lovelock gravity. In Section 5, we briefly discuss a couple of applications by examining how thin-shell wormholes and stellar models are governed by junction conditions in the context of $\mathcal{F}\left(R\right)=R+\beta R^2$ theory. Finally, we end with some concluding remarks in Section 6. For the reader who is uninterested in the illustrative examples of gravitational theories we have chosen, Section 2 suffices as a standalone explanation of the tapestry of techniques used in our method, but in our opinion, reviewing a couple of concrete examples in Section 3 or Section 4 should still be useful towards learning how to apply our method for other theories.

Conventions: We work in natural units where $\mathcal{G}=c=1$. Greek indices refer to all of spacetime, whereas Latin indices pertain to coordinates of the singular hypersurface $\Sigma$, with ”n” reserved for the direction normal to $\Sigma$. Our convention for the Riemann tensor is ${R^{\rho }}_{\sigma \mu \nu }={\partial }_{\mu }{\Gamma }_{\nu \sigma }^{\rho }+\dots$.

• Comparison with the existing results

Our results for $\mathcal{F}\left(R\right)$ theory and quadratic gravity refine the accuracy of those reported earlier in literature in a number of ways. For clarity, they are collected here:

1. As mentioned above, we explain in Section 4.1.1 explicitly the origin of a missing term in [14,15,17] for junctions conditions in the theory with $R^2$ action term.

2. In the papers [14,16,17,18], the authors argued that the naive appearance of products of delta function and delta derivatives in the equations of motion leads to inconsistency, and that there is no unambiguous manner to sum them up. As a result, it was assumed that $\left[K\right]=0$ so that R is non-singular. On this issue, we find that the use of delta-convergent sequences clearly quantifies the growth property of each integrand term, and provides an organizational principle for identifying regularity constraints and a well-defined procedure to sum up the residual finite terms, eventually leading us to the junction conditions. For us, singularities are handled by a careful procedure of regularization. This is inspired by the common practice of regularization in quantum field theory. We will show below that, rather than insisting on the use of distribution-valued curvature, our approach is more practical and much less restrictive and allows for more general solutions to the theory.

3. In [17], the authors posited that in quadratic gravity, the energy-momentum tensor should contain the distributional derivative of the $\delta$-function that they derived from the higher-order differential operators acting on a discontinuous Riemann tensor2. In this paper, we focus only on the explicit form of the generalized junction conditions, noting that only the delta-distribution component of $T_{\mu \nu }$ appears in final expression, since any regular component and those containing higher-order derivatives of the delta function do not survive the infinitesimal integral across $\Sigma$.

4. We should also mention that our method of derivation passes a stringent consistency test that is noticeably absent in the past literature, namely, that when applied to quadratic gravity, i.e., taking the Lagrangian $\mathcal{L}={\beta }_1{R}^2+{\beta }_2{R}^{\alpha \beta }{R}_{\alpha \beta }+{\beta }_3{R}^{\alpha \beta \mu \nu }{R}_{\alpha \beta \mu \nu }$, our junction conditions reduce to those of Gauss–Bonnet gravity $({\beta }_2=-4{\beta }_3=-4{\beta }_1)$ that were derived some time ago in [25] via a completely different approach—the boundary variation of a surface term [26] that accompanies the gravitational action. As already argued in [26,27] and definitively shown in [28], the equations of motion of a generic higher-derivative gravitational theory do not descend from a well-posed variational principle with Dirichlet conditions, and hence the approach of obtaining their appropriate junction conditions by boundary variation of surface terms is not always applicable3. An exception lies in the family of theories defined by Euler density terms which are the linear combinations of curvature invariants that generate only second-order field equations, and thus, in principle, there exist appropriate surface terms for them. The simplest example would be the Ricci scalar being the 2D Euler density, with the Gibbons–Hawking–York term $S_{GHY}=\frac{1}{8\pi }\int d^{d-1}x\sqrt{h}K$ being its surface term, and of which variation with respect to the induced metric yields (1). The junction condition derived in this manner nicely furnishes a consistency check for the other derivation route obtained by integrating across $\Sigma$. In Appendix C, we will demonstrate how the complicated junction conditions for a generic quadratic gravity theory reduce to the known ones for Gauss–Bonnet gravity. The regularity constraints need not be imposed precisely for ${\beta }_2=-4{\beta }_3=-4{\beta }_1$, and this also serves as a consistency check for the corresponding equations used in determining them. On this point, we note that previously in [17] where the authors proposed junction conditions for quadratic gravity, they imposed $\left[K_{ij}\right]=0$ even in the Gauss–Bonnet case, which contradicts the known ones derived by a boundary variation [25,29].

2. Generalized Junction Conditions—A General Approach

2.1. Some Preliminaries

For a generic gravitational theory with higher-order derivatives of the metric in the action, let us denote the field equations as

$${\tilde{G}}_{\mu \nu }\equiv R_{\mu \nu }-\frac{1}{2}R{g}_{\mu \nu }+\dots =8\pi T_{\mu \nu },$$

(4)

where the ellipses represent the additional terms descending from the modified action. In this work, we will mainly be studying gravitational theories where these terms are products of various contractions of the Riemann tensor and differential operators acting on them. On the RHS, we assume that the energy-momentum tensor harbors a singular source localized on some codimension-one hypersurface $\Sigma$. If we integrate across the infinitesimal width of $\Sigma$, only singular terms on both sides of (4) survive, and the resulting integral yields the appropriate junction conditions for the gravitational theory considered. The main purpose of this paper is to develop a conceptual and computational basis for this integration procedure. Fundamentally, apart from topological action terms which are not reflected in (4), the bulk equations of motion suffice as the starting point for the derivation of the junction conditions which can be regarded as consistency conditions for the bulk dynamics induced on $\Sigma$.

In the full generality, the induced metric on $\Sigma$$h_{ij}$ can be timelike, spacelike, or null. For the rest of the paper, we will mainly focus on the timelike case, which has garnered the most interest in the past literature. When $\Sigma$ is endowed with a singular energy-momentum tensor localized on it, it represents a singular thin shell of matter (for a spacelike $h_{ij}$, this would be a gravitational instanton, whereas a null $h_{ij}$ would represent a singular light cone).

In the most generic setting, $\Sigma$ is a codimension-one hypersurface that is embedded in a bulk manifold ($\mathcal{M}$) constructed from two distinct ones (${\mathcal{M}}_1,{\mathcal{M}}_2$) that are joined by identifying

$${\Sigma }_1\sim {\Sigma }_2\sim \Sigma ,$$

where ${\Sigma }_{1,2}$ are hypersurfaces in ${\mathcal{M}}_{1,2}$ that can be related via a homeomorphism. In this cut-and-paste procedure, the final manifold $\mathcal{M}$ is the union of $\Sigma$ and the interiors of ${\mathcal{M}}_{1,2}$ previously bounded by ${\Sigma }_{1,2}$. We can introduce a set of coordinates $\left({\zeta }^i\right)$ intrinsic to $\Sigma$ which is defined parametrically by

$$f\left(x^{\mu }\left({\zeta }^i\right)\right)=0,$$

with the induced metric

$$h_{ij}=g_{\mu \nu }\frac{\partial x^{\mu }}{\partial {\zeta }^i}\frac{\partial x^{\nu }}{\partial {\zeta }^j}.$$

(5)

Denoting the sign of $n^2$ by $\xi$, the unit normal vectors to $\Sigma$ read

$$n_{\mu }=\xi N\frac{\partial f}{\partial x^{\mu }},N={\left|g^{\alpha \beta }{\partial }_{\alpha }f{\partial }_{\beta }f\right|}^{-1/2},n^2=\xi ,$$

where $\xi =\{+1,-1,0\}$ for the cases of $\Sigma$ being timelike, spacelike, and null, respectively. The second fundamental form or the extrinsic curvature enacts an essential role in our discussion and can be defined as follows.

$$K_{ab}=h_a^m{h}_b^l{\nabla }_m{n}_l,n^a{K}_{ab}=0.$$

(6)

For our purpose, we find it useful to work in the Gaussian normal coordinate system in the local neighborhood of $\Sigma$ constructed from the family of non-intersecting geodesics orthogonal to $\Sigma$.

$$d{s}^2=\xi d{n}^2+h_{ij}(n,x^k)d{x}^{i}d{x}^j,$$

(7)

with $\Sigma$ defined by taking $n=0$, and $x^i$ are the intrinsic coordinates. Having briefly described the geometric characterization of $\Sigma$, let us now return to the field equations in (4). On the RHS, we claim that the energy-momentum tensor can be expressed as

$$T_{\mu \nu }=S_{\mu \nu }\delta \left(n\right)+\dots ,$$

(8)

where we have denoted the delta function source by $S_{\mu \nu }=S_{\mu \nu }\left(\overrightarrow{x}\right)$ and ⋯ denotes the regular part of the energy-momentum tensor as well as more singular parts that involve linear combination of derivatives of the delta function in general. That this is so follows from a known result in the theory of distribution (see, e.g., [30]) which states that any distribution with point support is a linear combination of derivatives of the delta distribution. Hence, the ansatz (8) indeed represents the most general singular source localized on a codimension-one hypersurface expressed in Gaussian normal coordinates4. Working in the local Gaussian normal chart, we can integrate (4) across the infinitesimal width of $\Sigma$ as follows.

$$\underset{ϵ\to 0^{+}}{lim}{\int }_{-ϵ}^{ϵ}dn{\tilde{G}}_{\mu \nu }=8\pi S_{\mu \nu }.$$

(9)

where, in particular, we note that only the delta term $S_{\mu \nu }$ in (8) survives the integral, and not the regular terms or derivatives of the delta function represented by the ellipses. The infinitesimal integral domain picks up only the $\delta$-distribution part of the energy-momentum tensor, since other components vanish in the $ϵ\to 0$ limit. Formally, we can see this by writing

$$\underset{ϵ\to 0^{+}}{lim}{\int }_{-ϵ}^{ϵ}dn={\int }_{-\infty }^{\infty }dn-\mathrm{p}.\mathrm{v}.{\int }_{-\infty }^{\infty }dn,$$

where the second term on the RHS is the Cauchy principal value of the first integral term. In local Gaussian normal coordinates, the Christoffel symbols read

$${\Gamma }_{nij}=-\frac{1}{2}{\partial }_n{h}_{ij}=-K_{ij},{\Gamma }_{nj}^i=K_j^i=\frac{1}{2}{h}^{im}{\partial }_n{h}_{jm}.$$

(10)

The discontinuity in ${\partial }_n{h}_{ij}$ or the extrinsic curvature induces various singular terms in ${\tilde{G}}_{\mu \nu }$, which is then related to the singular source $S_{\mu \nu }$ in (9).

A main theme of this work is to seek the conditions characterizing the embedding of $\Sigma$, under which the LHS of (9) is well-defined. This then leads to an appropriate set of junction conditions on the LHS of (9) upon performing the integral. We take this to be the operational meaning of the delta function source in the energy-momentum tensor, since the field equations should be understood in the sense of distribution when the energy-momentum tensor carries a delta-singular source.

Now in our derivation, the Gauss–Codazzi relations turn out to be crucial as they express various curvature quantities in terms of the intrinsic Riemann tensor and the extrinsic curvature tensor, elucidating the form of various singular terms in ${\tilde{G}}_{\mu \nu }$ easily. Recall that the Gauss–Codazzi equations read

$$\begin{array}{ccc}\hfill {\widehat{R}}^a{}_{bcd}& =& \mathcal{P}\left({R^a}_{bcd}\right)+\xi \left(K_c^a{K}_{bd}-K_d^a{K}_{bc}\right),\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill n^a{R}_{aqrs}{h}_b^q{h}_c^r{h}_d^s& =& \xi \left(D_d{K}_{bc}-D_c{K}_{bd}\right),\hfill \end{array}$$

(12)

where $\mathcal{P}$ refers to the indices being projected with the induced metric $h_{ab}$, $D_a$ is the projected covariant derivative defined with the induced metric, and hatted variables are geometrical quantities intrinsic to $\Sigma$. Using Gaussian normal coordinates and contracting some of the indices, we further obtain the following useful set of equations.

$$\begin{array}{ccc}\hfill {R^n}_{snv}& =& \xi \left(-{\partial }_n{K}_{vs}+K_{av}{K}_s^a\right),R_{nijk}=D_k{K}_{ij}-D_j{K}_{ik},\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill R_{ij}& =& \xi \left(-{\partial }_n{K}_{ij}+2K_j^a{K}_{ai}-K{K}_{ij}\right)+{\widehat{R}}_{ij},\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill R_{nn}& =& -{\partial }_{n}K-K^{ab}{K}_{ab}=-h^{ij}{\partial }_n{K}_{ij}+K^{ab}{K}_{ab},\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill R_{in}& =& D^k{K}_{ik}-D_{i}K,\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill R& =& \widehat{R}-\xi \left(2{\partial }_{n}K+K^2+K^{ab}{K}_{ab}\right).\hfill \end{array}$$

(17)

2.2. On the Standard Darmois–Israel Junction Conditions in GR

It is instructive to first review how one can obtain the standard Darmois–Israel junction conditions in ordinary GR before proceeding to more complicated gravitational theories. In (9), we now simply take ${\tilde{G}}_{\mu \nu }=G_{\mu \nu }$ and integrate

$$\underset{ϵ\to 0^{+}}{lim}{\int }_{-ϵ}^{ϵ}dn{G}_{\mu \nu }=8\pi S_{\mu \nu }.$$

(18)

On the LHS, only singular terms in the Einstein tensor survive in the limit $ϵ\to 0$. The singularity of these terms are of the Dirac delta type and can be traced to a discontinuity in the derivative of the metric in the direction normal to $\Sigma$ which, in the Gaussian normal chart, is effectively the discontinuous extrinsic curvature. We keep track of these terms appearing in the Einstein tensor which can be simplified to read

$$\begin{array}{ccc}\hfill & & \xi \left({\partial }_n{K}_j^i-{\delta }_j^i{\partial }_{n}K\right)+\xi \left(-K{K}_j^i+\dots \right)=8\pi T_j^i,\hfill \\ & & \xi \left(-{\nabla }_{i}K+{\nabla }_j{K}_i^j\right)=8\pi T_i^n,-\frac{1}{2}\widehat{R}+\frac{1}{2}\xi \left(K^2-\mathrm{Tr}\left(K^2\right)\right)=8\pi T_n^n,\hfill \end{array}$$

(19)

where the ellipses are terms which do not involve ${\partial }_n$. Integrating over $\Sigma$ as in (18) on both sides of the field equations then yields the junction conditions

$$\xi \left(-\left[K_{ij}\right]+\left[K\right]g_{ij}\right)=8\pi S_{ij},$$

(20)

which are the standard Darmois–Israel junction conditions in GR, and we have used the bracket $[\cdots ]$ to denote the jump or difference between the limiting values of the bracketed quantity on each side of $\Sigma$, for example, ${\left[K\right]\equiv K|}_{n=0^{+}}{-K|}_{n=0^{-}}$. We also note that the junction conditions in other directions are trivial, i.e., in (9),

$$\underset{ϵ\to 0^{+}}{lim}{\int }_{-ϵ}^{ϵ}dn{G}_{in}=\underset{ϵ\to 0^{+}}{lim}{\int }_{-ϵ}^{ϵ}dn{G}_{nn}=0.$$

(21)

Thus in GR, the localized singular source $S_{\mu \nu }$ cannot have components non-parallel to $\Sigma$. As we shall see later, this is not generally true for other gravitational theories. Another feature characterizing the junction condition (20) that differentiates a generic higher-derivative theory from GR is that the junction equation does not only involve bracketed terms but also acquires terms that represent averaging across $\Sigma$. Finally, we note that in general the energy-momentum tensor for diffeomorphism-invariant matter is conserved

$${\nabla }_{\mu }{T}^{\mu \nu }=0$$

(22)

due to the Noether theorem. This holds on-shell and is independent of which higher-derivative gravity theory it is coupled to. When a hypersurface $\Sigma$ is introduced and a localized energy momentum tensor of the form (8) is considered, the conservation of $T^{\mu \nu }$ implies that

$${\nabla }_i{S}^{ij}+K^j{}_i{S}^{in}=0,{\nabla }_i{S}^{in}=0,$$

(23)

for the components $S^{ij}\left(x\right),S^{in}\left(x\right),S^{nn}\left(x\right)$ of the stress tensor that lives on $\Sigma$. As mentioned above, $S^{in}\ne 0$ in general and this leads to the non-conservation ${\nabla }_i{S}^{ij}=G^i\ne 0$ of the stress tensor $S^{ij}$. The non-conservation indicates the presence of an external force $G^j:=-K^j{}_i{S}^{in}$ on the matter localized at the junction. The force is geometric in origin and is sourced by the non-trivial embedding of the surface. Moreover the force satisfies a divergence-free constraint as imposed by the second equation in (23), e.g., ${\nabla }_i{\left(K^{-1}\right)}_j^i{G}^j=0$ when $K_j^i$ is invertible. In GR, $S^{in}=0$ and so $S^{ij}$ is conserved. It is known that in GR, due to the Bianchi identity satisfied by the Einstein tensor, the conservation of the energy momentum is satisfied without requiring the metric to be on-shell. Therefore, the LHS of (20) must also be divergence-free identically. This can be checked explicitly. The same is true for the $\mathcal{F}\left(R\right)$ theory as a result of a generalized Bianchi identity [31] for the corresponding “generalized Einstein tensor”. Therefore, (23) are also satisfied in the $\mathcal{F}\left(R\right)$ theory without requiring the metric to be on-shell, and it follows that a number of identities must be satisfied for the junction terms on the LHS of the junction conditions. e.g., the LHS of (75), (77)–(82) for the case of $R+R^2$ gravity. These identities are automatically satisfied as a result of the generalized Bianchi identity discovered in [31]. We invite interested readers to check them explicitly. For more general higher-derivative gravity theory, e.g., $\mathcal{L}=\mathcal{L}(R_{\mu \nu },R^{\mu }{}_{\nu \alpha \beta })$, it is not clear that such generalized Bianchi identity exists, but the Equation (23) hold on-shell at the very least.

2.3. Delta Sequences and a Double Scaling Limit

A crucial ingredient in our formulation is a procedure that takes into account singular terms more divergent than delta function terms in the LHS of the field Equation (4), that may appear in the equations of motion upon assuming a discontinuous extrinsic curvature (note that this problem is absent in (19)). We demonstrate how to impose suitable regularity constraints on the extrinsic curvature such that the integral in (9) remains convergent, leading to generalized junction conditions which are the appropriate boundary conditions at $\Sigma$. We will find that an essential tool is to regard the delta function as a limit of a sequence of absolutely convergent functions, a delicate subject which we turn to next.

Since the singularities in (4) essentially arise from derivatives of the Heaviside step function parametrizing the discontinuity in the extrinsic curvature, we first consider a suitable representation of it by defining it as a limit $b\to 0$ of a sequence of classical functions $\Theta (n,b)$, with

$$\Theta \left(n\right)=\underset{b\to 0}{lim}\Theta (n,b),\Theta (n,b)=\frac{1}{2}+{\Theta }_X\left(\frac{n}{b}\right),$$

(24)

where we have introduced a function ${\Theta }_X\left(\frac{n}{b}\right)$ such that the delta function can be expressed as

$$\delta \left(n\right)=\underset{b\to 0}{lim}{\partial }_n\Theta (n,b)=\underset{b\to 0}{lim}\frac{1}{b}{\Theta }_X^{\prime }\left(X\right)\equiv \underset{b\to 0}{lim}\frac{1}{b}F\left(X\right),X\equiv \frac{n}{b},$$

(25)

The function $F\left(X\right)$ is sometimes known as a nascent delta function, with ${\Theta }_X\left(X\right)$ being its antiderivative. As a generalized function or distribution, we require that

$$\underset{b\to 0}{lim}{\int }_{-\infty }^{\infty }dn\frac{1}{b}F\left(X\right)f\left(n\right)=f\left(0\right),$$

(26)

for all continuous $f\left(n\right)$ with compact support which implies that $F\left(X\right)$ is normalized as

$${\int }_{-\infty }^{\infty }dXF\left(X\right)=1.$$

(27)

Formally, each nascent delta function gives rise to a distribution of the form

$${\tilde{f}}_{1/b}={\partial }_n\Theta (n,b)=\frac{1}{b}F(n/b),$$

(28)

and the sequence ${\tilde{f}}_1,{\tilde{f}}_2,\dots$ converges to the delta function distribution $\delta \left(n\right)$ if the limit defined on the LHS on (26) exists. Some popular choices of the nascent delta function5 $F\left(X\right)$ appearing in related literature are (i) $\frac{1}{\sqrt{\pi }}{e}^{-X^2}$ (Gaussian distribution), (ii) $\frac{1}{\pi (1+X^2)}$ (Cauchy distribution), (iii) $\frac{sin\left(X\right)}{\pi X}$ (sinc function), etc. For the purpose of our work here, we adopt some choice of $F\left(X\right)$ that is even in n:6

$$F\left(X\right)=F(-X),{\Theta }_X\left(X\right)=-{\Theta }_X\left(-X\right),$$

(29)

with $F\left(X\right)$ being infinitely differentiable, so that it is compatible with our use of its antiderivative ${\Theta }_X\left(X\right)$ to describe functions of varying degrees of smoothness.

Taking $F\left(X\right)$ to be an even function inherits the nature of $\delta \left(n\right)$ being an even distribution, and retains its original physical attribute of being symmetrical on either side of $\Sigma$. If we regard the nascent delta function as a probability density function describing how energy–mass is localized in $\Sigma$ (idealized by the $\delta$-distribution), then (29) translates to $n=0$ being the mean of the energy–mass distribution, which is expected on physical grounds 7. As we shall point out shortly in Section 2.4, the parity assignment (29) also leads to a natural generalization of the derivatives of the $\delta$-function when they are integrated against non-smooth functions 8.

Parametrizing the step function in this manner implies that for a function $g\left(n\right)$ which may not be smooth at (and after) a certain order of its derivative at $n=0$, we can consider it as the limit of a sequence of functions as follows.

$$g\left(n\right)=g_1\left(n\right)+\Theta \left(n\right)(g_2\left(n\right)-g_1\left(n\right))=\overline{g}\left(n\right)+\underset{b\to 0}{lim}{\Theta }_X\left(X\right)\left[g\left(n\right)\right],$$

(30)

where $g_1\left(n\right),g_2\left(n\right)$ are infinitely differentiable functions of n each analytically extending $g\left(n\right)$ beyond $n=0$, and we have defined

$$\overline{g}\left(n\right)=\frac{1}{2}\left(g_1\left(n\right)+g_2\left(n\right)\right),\left[g\left(n\right)\right]=g_2\left(n\right)-g_1\left(n\right),g_2\left(0\right)=g\left(0^{+}\right),g_1\left(0\right)=g\left(0^{-}\right).$$

For example, if $g\left(n\right)$ is of class $C^1$ with its second-order derivative being discontinuous at $n=0$, then $g_1\left(0\right)=g_2\left(0\right),g_1^{\prime }\left(0\right)=g_2^{\prime }\left(0\right),g_1^{″}\left(0\right)\ne g_2^{″}\left(0\right)$, etc. One can perform ordinary derivatives on $g\left(n\right)$ before finally taking the $b\to 0$ limit.

In (9), the integrand is generally a rather complicated function of various non-smooth curvature quantities, and it is crucial that we apply (30) for all terms consistently before evaluating the integral. Formally, this procedure yields a sequence of regular distributions (which is essentially a complicated function of the nascent delta function and other functions of the metric tensor analytic at $\Sigma$) that ensure the convergence of (9).

We will find that this manner of expressing the step function (and hence the delta functions and their derivatives) lends us a powerful language for simplifying and organizing various types of singular terms consistently in the integrand of (9). It also necessitates a more precise prescription of the infinitesimal integration procedure in (9), one which takes into account the relative scale separation between the sequence parameter b and the infinitesimal width $ϵ$ of the surface $\Sigma$. We now address this pivotal point that arises when we apply these considerations to our problem—the relative scaling of the thin-width parameter $ϵ$ and the nascent delta function parameter b.

Recall that in deriving the junction equations, we perform an integral with limits parametrizing the width of the infinitesimally thin $\Sigma$, and hence we send the integral limits to zero after integrating. The integral preserves only terms which accompany a delta function singularity. Consider again (26) with an infinitely differentiable $f\left(n\right)$, but with the integral limits of (9). We still expect to recover $f\left(0\right)$ on the RHS. Expanding $f\left(n\right)$ around $n=0$ to obtain

$$\underset{ϵ\to 0}{lim}\underset{b\to 0}{lim}{\int }_{-ϵ}^{ϵ}dn{\delta }_b\left(n\right)f\left(n\right)=\underset{ϵ\to 0}{lim}\underset{b\to 0}{lim}{\int }_{-ϵ/b}^{ϵ/b}dyF\left(y\right)\left(f\left(0\right)+by{f}^{\prime }\left(0\right)+\frac{1}{2}{b}^2{y}^2{f}^{″}\left(0\right)+\dots \right),$$

(31)

we find that we recover (26) upon taking the double scaling limit

$$ϵ\to 0,b/ϵ\to 0.$$

(32)

We note that on the contrary, if we specify b to vanish such that $b/ϵ\to \mathcal{O}\left(1\right)$ or the converse $ϵ/b\to 0$, then we do not recover $f\left(0\right)$. The physical reasoning behind (32) is simple and intuitive: the thin-width limit must be kept away from/taken after the limit of the sequence of nascent delta functions, or otherwise we are not genuinely integrating a singular delta function across $\Sigma$.

Recall that in the field equations, we seek a solution to

$${\tilde{G}}_{\mu \nu }=R_{\mu \nu }-\frac{1}{2}R{g}_{\mu \nu }+\dots =8\pi \left(S_{\mu \nu }\delta \left(n\right)+\dots \right),$$

(33)

where the ellipses refer to all regular components of the energy-momentum tensor and possibly a linear combination of higher-order derivatives of the delta function, all of which do not survive the integral of (9). A solution to (33) must thus integrate to give $S_{\mu \nu }$ in the following manner.

$$\underset{ϵ\to 0}{lim}\underset{b\to 0}{lim}{\int }_{-ϵ}^{ϵ}dn\left(R_{\mu \nu }-\frac{1}{2}R{g}_{\mu \nu }+\dots \right)=8\pi S_{\mu \nu },$$

(34)

where the limit in b is taken for the sequence of distributions defined by (28). We have to ensure that this sequence of integrals converges, which generically requires additional conditions to be imposed on the extrinsic curvature for the integral in (34) to be well-defined.

This double scaling limit also preserves the fact that we can obtain a well-defined operational meaning for the derivatives of the delta functions when they are integrated against smooth functions, i.e., that we have, in a weak distributional sense,

$${\delta }^{\left(k\right)}\left[g\left(n\right)\right]\sim {(-1)}^k{g}^{\left(k\right)}\left(0\right),$$

(35)

for a $g\left(n\right)$ analytic at $n=0$. This formula needs to be refined generally for a non-smooth $g\left(n\right)$—an important yet tricky point that is relevant for us since, as we shall see, the equations of motion of higher-derivative gravitational theories typically contain terms which are products of non-smooth functions and a number of delta functions each equipped with some order of derivative.

2.4. Generalized Distributions and Regularity Constraints

In this section, we present a well-defined procedure that extends the usual notion of distribution for the delta function (and its derivatives) to one that could be applied when they are integrated against non-smooth functions. Essentially for our purpose here, this will ultimately turn out to furnish a simple and clear method to identify and classify singular terms arising from the junction equations, leading to regularity constraints which we can impose on the extrinsic curvature to eliminate singular terms that render (9) to be ill-defined.

We begin with the basic examples of the $\delta \left(n\right),{\delta }^{\prime }\left(n\right)$.

2.4.1. Warm-Up: More about $\delta \left(n\right),{\delta }^{\prime }\left(n\right)$

Let us consider a discontinuous test function for the delta function which we express as $f\left(n\right)=\overline{f}\left(n\right)+{\Theta }_X\left(X\right)\left[f\left(n\right)\right]$ where $\overline{f}\left(n\right)=\frac{1}{2}\left(f_1\left(n\right)+f_2\left(n\right)\right), \left[f\left(n\right)\right]=f_2\left(n\right)-f_1\left(n\right)$, with both $f_{1,2}\left(n\right)$ being smooth functions that extend $f\left(n\right)$ across $\Sigma$. Within the infinitesimal domain width, we can expand $f\left(n\right)$ about the origin. In the notations introduced in the previous section, we have

$$\underset{ϵ\to 0}{lim}\underset{b\to 0}{lim}{\int }_{-ϵ}^{ϵ}dnf\left(n\right){\delta }_b\left(n\right)=\underset{ϵ\to 0}{lim}\underset{b\to 0}{lim}\sum _{l=0}\frac{b^l}{l!}{\int }_{-ϵ/b}^{ϵ/b}dX{X}^{l}F\left(X\right)\left({\overline{f}}^{\left(l\right)}\left(0\right)+{\Theta }_X\left(X\right)\left[f^{\left(l\right)}\left(0\right)\right]\right)$$

(36)

The only non-vanishing term is $l=0$, and since $F\left(X\right),{\Theta }_X\left(X\right)$ are even and odd in X, respectively, we obtain

$$\underset{ϵ\to 0}{lim}\underset{b\to 0}{lim}{\int }_{-ϵ}^{ϵ}dnf\left(n\right){\delta }_b\left(n\right)=\overline{f}\left(0\right),$$

(37)

a result which is naturally intuitive and reduces correctly to the expected one in the continuous limit. Now we can extend this calculation to derivatives of the delta function, for example, ${\delta }^{\prime }\left(n\right)$. Consider thus the integral

$$\underset{ϵ\to 0}{lim}\underset{b\to 0}{lim}{\int }_{-ϵ}^{ϵ}dnf\left(n\right){\delta }_b^{\prime }\left(n\right)=\underset{ϵ\to 0}{lim}\underset{b\to 0}{lim}\sum _{l=0}\frac{b^{l-1}}{l!}{\int }_{-ϵ/b}^{ϵ/b}dX{X}^l{F}^{\prime }\left(X\right)\left({\overline{f}}^{\left(l\right)}\left(0\right)+{\Theta }_X\left(X\right)\left[f^{\left(l\right)}\left(0\right)\right]\right).$$

(38)

For all $l>1$ this vanishes. For $l=1$, since $F^{\prime }\left(X\right),{\Theta }_X\left(X\right)$ are both odd in X, this term reads

$${\int }_{-\infty }^{\infty }dXX{F}^{\prime }\left(X\right){\overline{f}}^{\prime }\left(0\right)=-{\overline{f}}^{\prime }\left(0\right),$$

(39)

where we have integrated by parts. We are left with the $l=0$ term which is singular, since

$$\underset{b\to 0}{lim}\frac{1}{b}{\int }_{-\infty }^{\infty }dX{F}^{\prime }\left(X\right)\left(\overline{f}\left(0\right)+{\Theta }_X\left(X\right)\left[f\left(0\right)\right]\right)=\underset{b\to 0}{lim}\frac{\left[f\right(0\left)\right]}{b}{\int }_{-\infty }^{\infty }dX{F}^{\prime }\left(X\right){\Theta }_X\left(X\right).$$

(40)

The integral is then well-defined only if we impose from the outset the continuity of f at $n=0$,9

$$f_2\left(0\right)=f_1\left(0\right),$$

(41)

leading to

$$\underset{ϵ\to 0}{lim}{\int }_{-ϵ}^{ϵ}dnf\left(n\right){\delta }^{\prime }\left(n\right)=-{\overline{f}}^{\prime }\left(0\right).$$

(42)

The continuity condition in (41) for ${\delta }^{\prime }\left(n\right)$ is the simplest example of what we shall allude to as regularity constraints that we would need to solve for, and impose when we encounter some integrand in (9) that is typically a product of some nascent delta functions, each possibly equipped with some order of derivative. Subject to (41), we can understand ${\delta }^{\prime }\left(n\right)$ to yield (minus) the average of the first derivative of $f\left(n\right)$ evaluated at $n=0$ should the function be non-differentiable at $n=0$. In a similar vein, one can prove that more generally, we have

$$\underset{ϵ\to 0}{lim}{\int }_{-ϵ}^{ϵ}dnf\left(n\right){\delta }^{\left(k\right)}\left(n\right)={(-1)}^k{\overline{f}}^{\left(k\right)}\left(0\right),k\ge 1,$$

(43)

which is an intuitive generalization of the case where the test functions $f\left(n\right)$ are smooth. We note that the simple form of the RHS of (43) follows by virtue of the symmetry property of the nascent delta function in (29) which, if relaxed, leads to a more complicated relation. For example, if F has no definite parity, then generally, instead of (42), we have $-{\overline{f}}^{\prime }\left(0\right)-\left[f^{\prime }\left(0\right)\right]{\int }_{-\infty }^{\infty }dX{F}^2\left(X\right)X$, which depends on the precise form of F.

Usually, when $\delta \left(n\right),{\delta }^{\prime }\left(n\right)$ are regarded as distributions on some open subset of $\mathbb{R}$, the class of test functions are smooth, compactly supported functions. As shown above, for test functions which are non-smooth, $\delta \left(n\right),{\delta }^{\prime }\left(n\right)$ can still be regarded as distributions with point support at $n=0$, but for ${\delta }^{\prime }\left(n\right)$ to be a sensible distribution, we need to restrict the class of test functions to be at least continuous at $n=0$.

2.4.2. Regularity Constraints for Products of Nascent Delta Functions and Their Derivatives

In the Gaussian normal chart, integrating the bulk equations of motion (of a general higher-derivative gravitational theory) across $\Sigma$ is reduced to evaluating a set of one-dimensional integrals, with the integrand being some complicated product of various derivatives of the extrinsic curvature. Since, generally, we take the extrinsic curvature to be not necessarily continuous at $\Sigma$ in response to a delta-singular energy source, this implies that in the absence of some regularizing constraints, such an integral is generically singular.

In the same vein by which we have studied the basic examples of $\delta \left(n\right),{\delta }^{\prime }\left(n\right)$, we now explain how we can solve for regularity constraints for products of non-smooth functions, delta functions and their derivatives. As we have seen in the previous examples, it is useful to describe a non-smooth function with the ansatz (30) which, in turn, contains smooth functions that can be expanded around $n=0$ within the thin-width domain of integral.

This implies that we can break down the integral in (9) into a linear combination of integrals of the following form.

$$I(l,\overrightarrow{k})=\underset{ϵ\to 0}{lim}{\int }_{-ϵ}^{ϵ}dn{n}^l{\partial }_n^{k_1}\Theta \left(n\right){\partial }_n^{k_2}\Theta \left(n\right)\dots {\partial }_n^{k_j}\Theta \left(n\right),$$

(44)

where l and the indices $\overrightarrow{k}=\{k_1,k_2,\dots ,k_j\}$ are non-negative integers, with each $k_i$ indicating the order of derivative. Together with the scaling limit (32), we can express (44) as

$$I(l,\overrightarrow{k})=\underset{b\to 0}{lim}{b}^{l+1-{\sum }_{m=1}^j{k}_m}{\int }_{-\infty }^{\infty }dX{X}^l{\Theta }_X^{\left(k_1\right)}\left(X\right)\dots {\Theta }_X^{\left(k_j\right)}\left(X\right).$$

(45)

Such a term vanishes for $l\ge {\sum }_m{k}_m$, remains finite for $l+1-{\sum }_m{k}_m=0$, and diverges for

$$l+1-\sum _{m=1}^j{k}_m<0,l=j-\sum _{m=1}^j{k}_m\left(\mathrm{mod}2\right),$$

(46)

with the second condition in (46) being due to the fact that the integral $I(l,\overrightarrow{k})$ vanishes by virtue of ${\Theta }_X\left(X\right)=-{\Theta }_X(-X)$ for $l=1+j-{\sum }_m{k}_m$ (mod 2).

After summing up the linear combination of $I(l,\overrightarrow{k})$ defining (9), for various singular terms which diverge as $b^s$ for some negative index s, we can now sum them up and set the overall coefficient to vanish. This naturally leads to a regularity constraint that we have to impose separately for each order of singularity labeled by each distinct s. A caveat is that the definite integral in (45) is generally representation-dependent. Since we would like to eliminate all singular integrals for any choice of nascent delta function, we should then impose a stronger condition: that for each family of integrals labeled by the same $\overrightarrow{k}$, every singular l satisfying (46) then leads to a regularity constraint. This procedure then renders the integral in (9) to be well-defined for any choice of nascent delta function. In the following, we elaborate on several illustrative examples.

We begin with a basic example by taking ${\overrightarrow{k}}_0=\{1,1,1\}$, with

$$I(l,{\overrightarrow{k}}_0)=\underset{ϵ\to 0}{lim}{\int }_{-ϵ}^{ϵ}dn{n}^l{\left({\Theta }^{\prime }\left(n\right)\right)}^3=\underset{b\to 0}{lim}{b}^{l-2}{\int }_{-\infty }^{\infty }dX{X}^l{F}^3\left(X\right),$$

(47)

which diverges for $l=0$, vanishes for $l=1$ and $l\ge 3$, and for $l=2$, it evaluates to ${\int }_{-\infty }^{\infty }dX{X}^2{F}^3\left(X\right)$. If we replace $n^l$ in the integrand of (47) by a function $\varphi \left(n\right)$ that is analytic at $n=0$, then this implies that the integral diverges unless $\varphi \left(0\right)=0$, in which case we have

$$I(\varphi ,{\overrightarrow{k}}_0)\equiv \underset{ϵ\to 0}{lim}{\int }_{-ϵ}^{ϵ}dn\varphi \left(n\right){\left({\Theta }^{\prime }\left(n\right)\right)}^3=\frac{{\varphi }^{″}\left(0\right)}{2}{\int }_{-\infty }^{\infty }dX{X}^2{F}^3\left(X\right).$$

(48)

In this case, $\varphi \left(0\right)=0$ would be what we call a regularity constraint, which stipulates how fast $\varphi \left(n\right)$ should grow near the origin for the integral to converge. The integral (48) also demonstrates that the naive product of three $\delta$-functions can be understood as a proper distribution, provided we restrict the space of test functions to be those which grow at least as $\varphi \left(n\right)\sim n+\mathcal{O}\left(n^2\right)$ near the origin10. This serves as a simple example of how the use of nascent delta functions resolves the “ambiguity” that may arise in interpreting products of $\delta$-functions as articulated in [17].

As an another example, consider the sum

$$\sum _l{C}_{l}I(l,\overrightarrow{k}),$$

with a fixed $\overrightarrow{k}$, the set of regularity constraints are simply

$$C_l=0,\forall l\mathrm{satisfying}\left(46\right).$$

(49)

In (9), as we shall see through explicit examples in later sections, it involves a linear combination of integrals of the form $I(l,\overrightarrow{k})$ with a finite set of $\overrightarrow{k}$. The entire set of regularity constraints is then the union of all (49) associated with each $\overrightarrow{k}$.

It is also instructive to check how (46) applies to the previous basic examples of $\delta \left(n\right),{\delta }^{\prime }\left(n\right)$ in (36) and (38). We note that for the two terms on the RHS of (36), they correspond to $I(l,\{1\left\}\right)$ and $I(l,\{0,1\left\}\right)$ in the notation of (44). With ${\sum }_m{k}_m=1$ for each term, there is no solution to (46) and thus no regularity constraint is needed. For (38), each of the two terms on the RHS are of the form $I(l,\{2\left\}\right)$ and $I(l,\{0,2\left\}\right)$ respectively, with ${\sum }_m{k}_m=2$ for each term. From (46), we thus see that there is a regularity constraint needed for $I(l,\{0,2\left\}\right)$ (where $j=2$) associated with $l=0$, and this is simply (41).

After imposing the necessary regularity constraints in (9) following the above approach, we are left with a linear combination of well-defined integrals of the form

$$I(l,\overrightarrow{k})=I\left(\sum _{m=1}^j{k}_m-1,\overrightarrow{k}\right)={\int }_{-\infty }^{\infty }dX{X}^{({\sum }_{m=1}^j{k}_m-1)}{\Theta }_X^{\left(k_1\right)}\left(X\right)\dots {\Theta }_X^{\left(k_j\right)}\left(X\right).$$

(50)

Since the integrand of (50) is odd for an even j, such a term only contributes to the junction equations for odd values of j. For example, in (36), the eventual expression corresponds to the first term within the RHS bracket which is of the form (50) with $l=0, j=1, {\sum }_m{k}_m=1$, whereas for (38), we note that (39) is of the form (50) with $l=1, j=1, {\sum }_m{k}_m=2$. As another example, for the sum ${\sum }_l{C}_{l}I(l,\overrightarrow{k})$, the only finite term (after imposing (49)) is

$$C_{r}I\left(r,\overrightarrow{k}\right),r=\sum _{m=1}^j{k}_m-1,$$

(51)

which is non-vanishing only if j is odd. A caveat is that, similar to the definite integral in (46), an expression such as (50) is generally sensitive to the choice of the nascent delta function, apart from terms like

$$\int dX{\Theta }_X^{\prime },\int dX{X}^p{\Theta }_X^{(p+1)},\mathrm{or}\int dX{\left({\Theta }_X\right)}^p{\Theta }_X^{\prime },$$

(52)

which are independent of the choice. For example, consider the integral

$$I(2,\{1,1,1\})={\int }_{-\infty }^{\infty }dX{X}^2{\left({\Theta }_X^{\prime }\left(X\right)\right)}^3.$$

(53)

If we pick the nascent delta function to be the Gaussian $F\left(X\right)=\frac{1}{\sqrt{\pi }}{e}^{-X^2}$, then (53) evaluates to $1/\left(6\pi \sqrt{3}\right)$, whereas a choice of $F\left(X\right)=\frac{sin\left(X\right)}{\pi X}$ yields $1/\left(2{\pi }^2\right)$ instead. It may seem as though such junction terms are “regularization-dependent”. This feature is obviously absent in ordinary GR, and here we see that the higher-order nature of the field equations may probe the form of the nascent delta function. To write down consistent junction conditions that are insensitive to the choice of nascent delta functions, we could additionally set the coefficient $C_r$ in (51) to vanish each time it appears in (9). We are then left with junction terms arising from universal terms such as those in (52).

In most of our working examples for the rest of the paper, these representation-dependent terms turn out not to feature much. The only setting where it arises non-trivially in this work is the case of $R^3$ theory for which we find the terms

$$\int dX{{\Theta }^{\prime }}_X^2{\Theta }_{X}X,\int dX{\Theta }_X^{″}{\Theta }_X^{\prime }{\Theta }_X{X}^2,$$

to be manifest in the junction equations. As we shall demonstrate later, these terms would be absent if we further set $\left[K^{\prime }\right]=0$ as a regularity constraint for the $R^3$ theory’s junction equations. In the general case, our derivation procedure described above allows us to solve for the regularity constraints that will yield the final junction conditions to be representation-independent if desired. Nonetheless, it is noteworthy that the convergence of the integral (9) is compatible with the presence of these terms which, if allowed, implies that the specification of junction conditions is only complete with a choice of nascent delta function.

To summarize, we can now state explicitly how to read off the regularity constraints and junction conditions for (9). After parametrizing each discontinuous geometric quantity and its derivatives by ${\Theta }_X$, expanding the LHS of (9) would yield integrand terms typically of the form

$$I_{\overrightarrow{k}}\equiv \varphi \left(n\right){\Theta }_X^{\left(k_1\right)}{\Theta }_X^{\left(k_2\right)}\dots {\Theta }_X^{\left(k_j\right)},$$

(54)

for some vector index $\overrightarrow{k}$, with $\varphi \left(n\right)$ analytic at $\Sigma$, then the regularity constraints associated with $I_{\overrightarrow{k}}$ are

$${\varphi }^l\left(0\right)=0,\forall l\le \sum _{m=1}^j{k}_m-2,l=j-\sum _{m=1}^j{k}_m\left(\mathrm{mod}2\right),$$

(55)

with the junction term induced by integrating $I_{\overrightarrow{k}}$ across $\Sigma$ being

$$J_{\overrightarrow{k}}=\frac{1}{(\left({\sum }_m{k}_m\right)-1)!}{\varphi }^{(\left({\sum }_m{k}_m\right)-1)}\left(0\right)\left({\int }_{-\infty }^{\infty }dX{X}^{\left({\sum }_m{k}_m\right)-1}{\Theta }_X^{\left(k_1\right)}\dots {\Theta }_X^{\left(k_j\right)}\right).$$

(56)

The final junction condition arising from (9) is then the sum of all the junction terms, each of the form (56), subject to us imposing all the regularity constraints, each of the form (55).

2.5. The Relation to Hadamard Regularization

In this section, we point out a relation between the regularity constraints and Hadamard regularization [21]—a well-studied procedure to regularize divergent integrals typically encountered in the theory of singular integral operators, and also commonly invoked when one handles distributions defined by divergent integrals (for an emphasis in the theory of distributions, see, for example, [22]).

Let us first briefly review the notion of Hadamard regularization with a simple example. Consider the (divergent) integral

$$I_H={\int }_0^{\infty }dx\frac{\varphi \left(x\right)}{x^{3/2}}=\underset{ϵ\to 0}{lim}{\int }_{ϵ}^{\infty }dx\frac{\varphi \left(x\right)}{x^{3/2}},$$

where $\varphi \left(x\right)$ is regular and continuous at the origin. From the mean value theorem, we can write $\varphi \left(x\right)=\varphi \left(0\right)+x{\varphi }^{\prime }\left(tx\right)\equiv \varphi \left(0\right)+x\phi \left(x\right),0<t<1$, and thus we can express the above integral as

$$I_H=\underset{ϵ\to 0}{lim}\frac{2\varphi \left(0\right)}{ϵ}+\underset{ϵ\to 0}{lim}{\int }_{ϵ}^{\infty }dx\frac{\phi \left(x\right)}{\sqrt{x}}.$$

(57)

The second integral in (57) converges and we define it to be Hadamard finite part of $I_H$, writing

$$FP{\int }_0^{\infty }dx\frac{\varphi \left(x\right)}{x^{3/2}}=\underset{ϵ\to 0}{lim}{\int }_{ϵ}^{\infty }dx\frac{\varphi \left(x\right)}{x^{3/2}}-\underset{ϵ\to 0}{lim}\frac{2\varphi \left(0\right)}{\sqrt{ϵ}}.$$

(58)

We have thus extracted the singular piece in the integral, with the first term of RHS of (58) being the residual finite integral. The above example can be generalized to a similar regularization of the divergent integral

$$I_D={\int }_{-\delta }^{\delta }dxf\left(x\right)\varphi \left(x\right),\delta >0,$$

where we take $\varphi \left(x\right)$ to be analytic at the origin, and ${f\left(x\right)\left|x\right|}^m$ has an algebraic singularity of order m at $x=0$, i.e., m is some smallest positive integer such that ${f\left(x\right)\left|x\right|}^m$ is locally integrable.

We can extend this notion to $f\left(x\right)$ being a distribution, the Hadamard regularization of which is then defined as

$$FP\langle f,\varphi \rangle ={\int }_{-\delta }^{\delta }dxf\left(x\right)\left(\varphi \left(x\right)-\left(\varphi \left(0\right)+{\varphi }^{\prime }\left(0\right)x+\dots +\frac{1}{m!}{\varphi }^{(m-1)}\left(0\right)x^{m-1}\right)\right)\Theta \left(1-\frac{\left|x\right|}{ϵ}\right).$$

(59)

We now examine how (59) is relevant for our context. From (44) and (45), we see that the sequence (in parameter b)

$${\Delta }_{\overrightarrow{k}}(n,b)\equiv {\Theta }^{\left(k_1\right)}(n,b)\dots {\Theta }^{\left(k_j\right)}(n,b),$$

converges to a distribution ${\Delta }_{\overrightarrow{k}}\left(n\right)$ that is of order of singularity $\left({\sum }_{m=1}^j{k}_m\right)-1$, and further, upon integrating it against some function $\varphi \left(n\right)$ that is analytic at $n=0$, the Hadamard-regularized distribution reads

$$\begin{array}{ccc}\hfill FP\langle {\Delta }_{\overrightarrow{k}}\left(n\right),\varphi \left(n\right)\rangle & =& \underset{ϵ\to 0}{lim}{\int }_{-ϵ}^{ϵ}dn{\Delta }_{\overrightarrow{k}}\left(n\right)\hfill \\ & & \times \left(\varphi \left(n\right)-\left(\varphi \left(0\right)+{\varphi }^{\prime }\left(0\right)n+\dots +\frac{1}{({\sum }_m{k}_m-2)!}{\varphi }^{({\sum }_m{k}_m-2)}\left(0\right)n^{{\sum }_m{k}_m-2}\right)\right),\hfill \end{array}$$

(60)

where, now, $\delta =ϵ=0^{+}$ for the relevant integral domain. Comparing (60) against (55) and (56), and recalling that $I(l,\overrightarrow{k})$ vanishes for $l=1+j-{\sum }_m{k}_m\left(\mathrm{mod}2\right)$, it is then clear that:

• In the absence of the regularity constraints, the junction condition is nothing but the Hadamard-finite part of the LHS of (9).
• Imposing the regularity constraints then ensures that the integral in (9) converges, being trivially equivalent to its Hadamard regularization.

Although it is nice to have recognized that our derivation of the junction conditions (and regularity constraints) admits a natural interpretation in terms of Hadamard regularization, we stress that our method is independently consistent, and can be understood and implemented without alluding to the latter.

3. Generalized Junction Conditions for Gravitational Theories with Quadratic Terms

In this section, we derive the generalized junction conditions for the following class of gravitational theories with action terms quadratic in (various contractions of) the Riemann tensor.

$${\mathcal{L}}_{quad}=\frac{1}{16\pi }\left(R+{\beta }_1{R}^2+{\beta }_2{R}_{uv}{R}^{uv}+{\beta }_3{R}_{\alpha \beta \mu \nu }{R}^{\alpha \beta \mu \nu }\right).$$

(61)

We integrate the equations of motion ${\tilde{G}}_{\mu \nu }=8\pi T_{\mu \nu }$ across the infinitesimally thin surface ${\int }_{-ϵ}^{ϵ}dn{\tilde{G}}_{\mu \nu }$, identify the regularity constraints, and derive the final explicit covariant form of the generalized junction condition. The equations of motion read

$$\begin{array}{ccc}& & G_{\alpha \beta }+2{\beta }_{1}R{R}_{\alpha \beta }-4{\beta }_3{R}_{\alpha \mu }{R}_{\beta }^{\mu }+2{\beta }_3{R}_{\alpha \rho \mu \nu }{R_{\beta }}^{\rho \mu \nu }+(2{\beta }_2+4{\beta }_3)R_{\alpha \mu \beta \nu }{R}^{\mu \nu }\hfill \\ & & -2({\beta }_1+\frac{1}{2}{\beta }_2+{\beta }_3){\nabla }_{\alpha }{\nabla }_{\beta }R+({\beta }_2+4{\beta }_3)R_{\alpha \beta }-\frac{1}{2}{g}_{\alpha \beta }(-(4{\beta }_1+{\beta }_2)R\hfill \\ & & +{\beta }_1{R}^2+{\beta }_2{R}_{\mu \nu }{R}^{\mu \nu }+{\beta }_3{R}_{\rho \sigma \mu \nu }{R}^{\rho \sigma \mu \nu })=8\pi T_{\alpha \beta }.\hfill \end{array}$$

(62)

To proceed, we will integrate ${\tilde{G}}_{\mu \nu }$ across $\Sigma$. In Appendix A, we collect several useful identities that we developed for evaluating the integrals easily. Using these identities, the results obtained by integrating across $\Sigma$ for the various terms in ${\tilde{G}}_{ij}, {\tilde{G}}_{in}, {\tilde{G}}_{nn}$ are collected in Appendix B. We will find that in contrast to the case in GR, apart from bracketed quantities, the junction terms also involve averaged quantities across $\Sigma$. In particular, the expression $\frac{1}{3}\left(\overline{fg}+2\overline{f}\overline{g}\right)\left[h\right]\equiv \overline{fg}\left[h\right]$ turns out to occur frequently (see (A6) in Appendix A).

3.1. Regularity Constraints

We now examine the conditions under which the integration is well-defined by summing up all singular terms arising from integrating the field equations. As shown in the Appendix B, it is remarkable that every such singular term is associated with the same divergent term

$$\underset{b\to 0}{lim}\underset{ϵ\to 0}{lim}{\int }_{-ϵ}^{ϵ}dn{({\partial }_n\Theta (n,b))}^2=\underset{b\to 0}{lim}\frac{1}{b}{\int }_{-\infty }^{\infty }dX{F}^2\left(X\right).$$

(63)

In the following, we gather the coefficients from each set of components of the field equations which are various functions of jumps of the extrinsic curvature.

From ${\tilde{G}}_{ij}$, the various singular terms sum up to read

$$\begin{array}{ccc}\hfill I_{ij}^{\left(sing\right)}& =& \left(4{\beta }_3-{\beta }_2-8{\beta }_1\right)\left[K\right]\left[K_{ij}\right]-2\left({\beta }_2+4{\beta }_3\right)\left[K_i^l\right]\left[K_{lj}\right]\hfill \\ \hfill & & +g_{ij}\left((2{\beta }_1+\frac{1}{2}{\beta }_2){\left[K\right]}^2-(\frac{1}{2}{\beta }_2+2{\beta }_3)\left[K^{ab}\right]\left[K_{ab}\right]\right),\hfill \end{array}$$

(64)

whereas from ${\tilde{G}}_{nn}$, the various singular terms sum up to be

$$I_{nn}^{\left(sing\right)}=\frac{3}{2}\left({\left[K\right]}^2(4{\beta }_1+{\beta }_2)+\left[K^{ab}\right]\left[K_{ab}\right]({\beta }_2+4{\beta }_3)\right).$$

(65)

Since there are no other singular terms from ${\tilde{G}}_{in}=0$, we proceed to set both (64) and (65) to vanish. From (64), we can take its trace to obtain

$$\left(4{\beta }_3+(2d-8){\beta }_1+(\frac{d}{2}-1){\beta }_2\right){\left[K\right]}^2-({\beta }_2+4{\beta }_3)(2+\frac{d}{2})\left[K^{ab}\right]\left[K_{ab}\right]=0.$$

(66)

where d is spacetime dimensionality. We can read off a few straightforward classes of solutions to $I_{ij}^{\left(sing\right)}=I_{nn}^{\left(sing\right)}=0$ as follows:

(I)

We first search for points in the moduli space such that there is no additional constraints to be imposed on the extrinsic curvature. From (64), we can read off

$${\beta }_2+4{\beta }_3=0,4{\beta }_3-{\beta }_2-8{\beta }_1=0,$$

(67)

which also solves (65) and leads to

$${\beta }_2=-4{\beta }_3=-4{\beta }_1.$$

(68)

This is precisely the combination for the Gauss–Bonnet theory! Recall that the Einstein-Gauss–Bonnet theory is defined with the following addition

$$\alpha \left(R^2-4R_{\alpha \beta }{R}^{\alpha \beta }+R^{\alpha \beta \mu \nu }{R}_{\alpha \beta \mu \nu }\right),$$

(69)

to the Einstein–Hilbert action, with $\alpha$ being some constant parameter (that could be ∼$l_s^2$ where $l_s$ is the string length if there is a string-theoretic origin [26]). In [17], the authors imposed the smoothness condition in which terms naively containing product of delta-functions should be forbidden in the action from the outset. This implies that the extrinsic curvature has to be continuous, and it was argued that the junction conditions are identical to GR for the theory with an added Gauss–Bonnet term in the action. Further, since $\left[K_{ab}\right]=0$, no thin-shell singular sources should be permitted. However, we find that in the related past literature11, a consistent junction condition for Einstein–Gauss–Bonnet theory has been presented in a few papers (which is different from the above conclusion). In the next section, we will show that our results recover those of [25,26,29,38] in the topological limit. This serves as a stringent consistency check for many of our equations, and demonstrates definitively that the Gauss–Bonnet theory does have non-trivial junction conditions.

(II)

Since imposing $\left[K_{ij}\right]=0$ naturally removes all singular terms, we look for less-stringent conditions on the extrinsic curvature. From the form of (64) and (65), we find the following class of solutions

$$\left[K\right]=0,{\beta }_2+4{\beta }_3=0,$$

(70)

corresponding to a family of theories for which we can set just the trace of the extrinsic curvature instead of all its components to vanish. In the space of the couplings $({\beta }_1,{\beta }_2,{\beta }_3)$, this is a plane containing the Gauss–Bonnet “line” (68).

(III)

Finally, we have the trivial solution

$$\left[K_{ij}\right]=0,$$

(71)

which is accompanied by no other constraints on the coupling parameters. There are still non-trivial junction equations to write down even in this case, since, as we have seen earlier, in general, these equations sometimes involve normal derivatives of the extrinsic curvature, i.e., while the first (normal) derivative of the metric has to be continuous, the higher-order ones need not be.

We have focused on the case where $\Sigma$ is timelike for definiteness. Nonetheless, this derivation can be repeated in an identical fashion for a spacelike $\Sigma$ since the difference lies in a few signs to be switched on in the Gauss–Codazzi relations. In particular, we find that the regularity constraints for the spacelike case remain the same, and the classification of junction conditions presented above is also preserved. Explicitly, restoring $\xi$ (the sign of $n^2$), we find that (64) and (65) generalize to

$$\begin{array}{ccc}\hfill I_{ij}^{\left(sing\right)}& =& \left(4{\beta }_3-{\beta }_2-8{\beta }_1\right)\left[K\right]\left[K_{ij}\right]-2\left({\beta }_2+4{\beta }_3\right)\left[K_i^l\right]\left[K_{lj}\right]\hfill \\ \hfill & & +g_{ij}\left((2{\beta }_1+\frac{1}{2}{\beta }_2){\left[K\right]}^2-\left(\frac{1}{2}{\beta }_2+2{\beta }_3\right)\left[K^{ab}\right]\left[K_{ab}\right]\right)\hfill \\ & & +(\xi -1)(4{\beta }_1+{\beta }_2)\left[K\right]\left(\left[K\right]g_{ij}-\left[K_{ij}\right]\right).\hfill \end{array}$$

(72)

There are no singular terms identically from ${\tilde{G}}_{in}$, whereas for ${\tilde{G}}_{nn}$, the various singular terms sum up to be

$$\begin{array}{ccc}\hfill I_{nn}^{\left(sing\right)}& =& \frac{3\xi }{2}\left({\left[K\right]}^2(4{\beta }_1+{\beta }_2)+\left[K^{ab}\right]\left[K_{ab}\right]({\beta }_2+4{\beta }_3)\right)\hfill \\ & & +(1-\xi ){\left[K\right]}^2\left(2{\beta }_1+\frac{3}{2}{\beta }_2+4{\beta }_3\right).\hfill \end{array}$$

(73)

It is then straightforward to show that setting $I_{ij}^{\left(sing\right)}$ and $I_{nn}^{\left(sing\right)}$ to vanish gives the same classification of junction conditions, as described earlier.

In the following section, we simplify and present the final form of the junction conditions for Class II and III theories with their respective regularity constraints implied. It is worthwhile to note that the finite terms multiplied to the representation-dependent factor ${\int }_{-\infty }^{\infty }dX\Theta X{F}^2\left(X\right)$ sum up to zero in the junction equations, since from integrating ${\tilde{G}}_{ij}$ across $\Sigma$, we find that these terms assemble to read

$$2\left((4{\beta }_1+{\beta }_2){\left[K\right]}^2-({\beta }_2+4{\beta }_3)\left[K^{ab}\right]\left[K_{ab}\right]\right)\left[K_{ij}\right]\times {\int }_0^{\infty }dXX{F}^2\left(X\right),$$

(74)

and one can check that for each of the above classes of theories, upon imposing their respective smoothness conditions, (74) vanishes exactly. The final junction equations are thus insensitive to the choice of the nascent delta function for all three families of theories.

3.2. Junction Conditions

Having solved for the appropriate regularity constraints, we can now impose them on the set of general finite junction terms derived earlier, and obtain the generalized Darmois–Israel junction conditions. For a generic choice of coupling parameters, the resulting junction equations can be rather elaborate. In the following, we present the explicit junction equations for each of the three classes of theories defined in the previous section, including the ordinary Darmois–Israel junction terms. On the RHS of the each equation, the boundary energy-momentum tensor is defined as the singular part of the bulk energy-momentum tensor localized within $\Sigma$, denoted by $S_{\alpha \beta }={lim}_{ϵ\to 0}{\int }_{-ϵ}^{ϵ}{T}_{\alpha \beta }$.

For Case (I), we find that our results reduce nicely to those for Gauss–Bonnet obtained independently by the method of boundary variation. This topological theory is defined by the line of couplings ${\beta }_1={\beta }_3=-\frac{1}{4}{\beta }_2$, and we find that the junction conditions simplify to read

$$-2{\beta }_1\left[3J_{ij}-J{h}_{ij}\right]-4{\beta }_1\left[K^{cd}\right]P_{icdj}+\left[K\right]h_{ij}-\left[K_{ij}\right]=8\pi S_{ij},$$

(75)

where

$$\begin{array}{ccc}\hfill J_{ij}& \equiv & \frac{1}{3}\left(2K{K}_{ic}{K}_j^c+K_{cd}{K}^{cd}{K}_{ij}-2K_{ic}{K}^{cd}{K}_{dj}-K^2{K}_{ij}\right),J\equiv J_k^k\hfill \\ \hfill P_{icdj}& \equiv & \frac{1}{2}\widehat{R}{g}_{i[d}{g}_{j]c}+{\widehat{R}}_{icdj}+{\widehat{R}}_{c[d}{g}_{j]i}-{\widehat{R}}_{i[d}{g}_{j]c}.\hfill \end{array}$$

(76)

We note that the bracketed rank-4 curvature tensor in the second line turns out to be the divergence-free component of the intrinsic Riemann tensor. There are no junction terms arising from integrating ${\tilde{G}}_{in}, {\tilde{G}}_{nn}$. We also note that (75) is equivalent to the junction condition derived in earlier literature [25,26,29,38] by taking the boundary variation of the Gauss–Bonnet surface term. Thus, this gives a strong consistency check of our general derivation which is noticeably absent in other previous proposals for junction conditions in quadratic gravity [17,18]. Since this is an important point, in Appendix C, we present a detailed proof of how our equations reduce to (75).

We should mention that a previous work in [17], which proposes junction conditions for quadratic gravity, fails this consistency test. The authors argued that one has to take the extrinsic curvature to be continuous at $\Sigma$, even for such Gauss–Bonnet case, and its junction conditions should read simply as $\left[K_{ij}\right]=0$. This contradicts the result obtained by either the boundary variation of the Gauss–Bonnet surface term, or integrating the equations of motion across $\Sigma$. In contrast, our derivation elucidates how various junction terms involving both averaged and jump quantities assemble nicely in the Gauss–Bonnet limit to yield (75), consistent with what we expect from the action principle.

For Case (II) where we take $\left[K\right]=0$, parametrizing this class of theories by $\{{\beta }_1,{\beta }_3\}$, we have the junction equations

$$\begin{array}{ccc}\hfill & & 2({\beta }_1-{\beta }_3)(-g_{ij}\left[2K^{″}+4K{K}^{\prime }+K{K}^{ab}{K}_{ab}+n^{\alpha }{\nabla }_{\alpha }\left(\mathrm{Tr}\left(K^2\right)\right)\right]+4\overline{K_{ij}{K}^{ab}}\left[K_{ab}\right]\hfill \\ & & -\left[K_{ij}\right]\overline{R}+\overline{K_{ij}}[K^{ab}{K}_{ab}+6K^{\prime }])-2{\beta }_3[J_{ij}-\frac{1}{3}\left(\left(3K{K}_{cd}{K}^{cd}-2K^{cd}{K}_{ac}{K}_d^a\right)h_{ij}\right)\hfill \\ & & +2K^{cd}\left(\frac{1}{2}\widehat{R}{g}_{i[d}{g}_{j]c}+{\widehat{R}}_{icdj}+{\widehat{R}}_{c[d}{g}_{j]i}-{\widehat{R}}_{i[d}{g}_{j]c}\right)-K_{ij}]=8\pi S_{ij},\hfill \end{array}$$

(77)

$$2({\beta }_1-{\beta }_3)h_i^{\alpha }{\nabla }_{\alpha }[2K^{\prime }+K^{ab}{K}_{ab}]=8\pi S_{in},$$

(78)

$$-2({\beta }_1-{\beta }_3)\left[2K{K}^{\prime }+K{K}^{ab}{K}_{ab}\right]=8\pi S_{nn},$$

(79)

where $K^{\prime }=n^{\alpha }{\nabla }_{\alpha }K, K^{″}=n^{\alpha }{n}^{\beta }{\nabla }_{\alpha }{\nabla }_{\beta }K$, restoring covariance in notation. In the ${\beta }_1\to {\beta }_3$ limit, we recover the Gauss–Bonnet junction equations with the additional constraint $\left[K\right]=0$. In the ${\beta }_3=0$ limit, we recover the simplest example of $\mathcal{F}\left(R\right)$-type gravitational theories, of which Lagrangian is an analytic function of the Ricci scalar R.

For Case (III), where $\left[K_{ij}\right]=0$, we have the set of equations

$$\begin{array}{ccc}\hfill -h_{ij}(4{\beta }_1+{\beta }_2)\left(\left[K^{″}\right]+2K\left[K^{\prime }\right]+\frac{1}{2}\left[{\left(\mathrm{Tr}\left(K^2\right)\right)}^{\prime }\right]\right)& +& 4(3{\beta }_1+{\beta }_2+{\beta }_3)K_{ij}\left[K^{\prime }\right]\hfill \\ \hfill +({\beta }_2+4{\beta }_3)\left(K_{(i}^l\left[K_{j)l}^{\prime }\right]-K\left[K_{ij}^{\prime }\right]\right)& =& 8\pi S_{ij},\hfill \end{array}$$

(80)

$$\begin{array}{ccc}\hfill 4({\beta }_1+\frac{1}{2}{\beta }_2+{\beta }_3)h_i^{\alpha }{\nabla }_{\alpha }\left[K^{\prime }\right]+({\beta }_2+4{\beta }_3)n^{\alpha }\left[{\nabla }_{\alpha }{R}_{in}\right]& =& 8\pi S_{in},\hfill \end{array}$$

(81)

$$\begin{array}{ccc}\hfill -(4{\beta }_1+{\beta }_2)\left[K{K}^{\prime }\right]-\frac{1}{2}({\beta }_2+4{\beta }_3)\left[n^{\alpha }{\nabla }_{\alpha }\left(\mathrm{Tr}\left(K^2\right)\right)\right]& =& 8\pi S_{nn}.\hfill \end{array}$$

(82)

As we see here, away from the Gauss–Bonnet limit, it could have non-vanishing components in orthogonal directions. Note that we expect all terms to vanish in the Gauss–Bonnet limit in which setting $\left[K_{ij}\right]$ to vanish implies the absence of junction conditions.

Apart from these generic cases, let us note that, in principle, there could exist other solutions to the vanishing of the singular terms $I_{ij}^{\left(sing\right)}, I_{nn}^{\left(sing\right)}$ in (64) and (65), which would be more fine-tuned than the three classes of solutions we presented above. For example, consider the ansatz

$$\left[K_{ij}\right]=\Omega \left(\overrightarrow{x}\right)h_{ij},$$

(83)

where $\Omega \left(\overrightarrow{x}\right)$ is an arbitrary function of the coordinates of $\Sigma$. Substituting this ansatz into ${\tilde{G}}_{ij},{\tilde{G}}_{nn}$, we find that we can eliminate all singular terms provided we have

$$4{\beta }_3+{\beta }_2+d(4{\beta }_1+{\beta }_2)=0,$$

(84)

where d is the spacetime dimensionality. However, away from the Gauss–Bonnet “line” of couplings, for this family of theories, the finite terms multiplied to the representation-dependent factor ${\int }_{-\infty }^{\infty }dX\Theta X{F}^2\left(X\right)$ do not sum up to vanish in the resulting junction conditions, in contrast to the three regular classes of theories we presented.

Before we move on to more general higher-derivatives theory, let us give some comments on the nature of the regularity constraints in general.

(1)

In the above, we have seen that apart from the Gauss–Bonnet case, where there are no additional constraints to be imposed on $\left[K_{ij}\right]$, for Cases (II) and (III), one needs to impose additional constraints on $\left[K_{ij}\right]$ (i.e., $\left[K\right]=0$ for (II) and $\left[K_{ij}\right]=0$ for (III)) to preserve regularity of the integrated field equations. The presence of singular terms in the integrated equation of motion is, in fact, a generic feature of higher-derivatives gravity theories, and their removal requires the imposition of regularity constraints that are conditions on the extrinsic curvature of the junction $\Sigma$.

(2)

The regularity constraints are not perturbative in the coupling parameters, since the higher-order field equations contain singularities (e.g., (64) and (65) for the quadratic gravity case) that are induced by the discontinuous extrinsic curvature. The singularities are present as long as the couplings are non-vanishing, even when they are small.

(3)

As the regularity constraints impose conditions on the extrinsic curvature of the junction hypersurface, one may wonder whether these constraints may lead to an overdetermined system of equations. The answer will depend on the kind of question asked. Let us analyze this carefully. In general, the regularity constraints are simply boundary conditions on the normal derivatives of the metric at $\Sigma$. In addition, as the order of derivatives is less than that of the equation of motion, the regularity constraints are consistent Cauchy data and do not overdetermine the bulk equation of motion. For example, given an incident field configuration on one side of $\Sigma$, one can consistently invoke the regularity constraints to determine the outgoing configuration on the other side of $\Sigma$; or, one can take $\Sigma$ as a brane, consider the combined gravity–brane matter system, and solve for consistent brane configuration together with its gravitating background. These are some of the popular applications of junction conditions in the literature, and the inclusion of regularity constraints in the higher-derivative generalization is fully consistent. On the other hand, if the problem involves a fixed gravity background with specified asymptotics, for example, as in the construction of wormholes in a fixed background, then the regularity constraints may overconstrain the system. This is so because the hypersurface (co-dimension 1) is specified by a single embedding function, and so if the system of regularity constraints consists of more than one independent condition, then there is no solution in general except for cases where the system of regularity constraints reduces to a single equation due to symmetry. The fact that there is a lack of solution to certain kinds of questions does not imply that the higher-derivative gravity theory is inconsistent. It simply means that some junction configurations are not admissible in the considered higher-derivatives theory. It may also be applied positively. For example, one can take a brane configuration in GR that is believed to exist in the quantum theory (e.g., a certain thin-shell D-brane) as a selection criterion for the possible higher-derivative correction terms to GR.

4. Junction Conditions for Other Examples of Higher-Derivative Theories

4.1. $\mathcal{F}\left(R\right)$ Theories

A class of gravitational theories considered in literature on modified gravity are $\mathcal{F}\left(R\right)$ theories, which refer to a Lagrangian that is an analytic function of the Ricci scalar. The equations of motion read

$${\partial }_R\mathcal{F}\left(R\right)R_{\mu \nu }-\frac{1}{2}\mathcal{F}\left(R\right)g_{\mu \nu }+\left(g_{\mu \nu }-{\nabla }_{\mu }{\nabla }_{\nu }\right){\partial }_R\mathcal{F}\left(R\right)=8\pi T_{\mu \nu }.$$

(85)

4.1.1. More about Junction Conditions of $R^2$ Theory

In the previous section, we covered a simple example of such theories: the $\mathcal{F}\left(R\right)=R+{\beta }_1{R}^2$ Lagrangian which belongs to the class of theories with quadratic curvature invariants that admit the regularity constraint $\left[K\right]=0$. Setting ${\beta }_3=0$ in Equations (77)–(79), we obtain the junction conditions to be

$$\begin{array}{ccc}\hfill & & -\left[K_{ij}\right]+2{\beta }_1(-g_{ij}\left[2K^{″}+4K{K}^{\prime }+K{K}^{ab}{K}_{ab}+{\left(\mathrm{Tr}\left(K^2\right)\right)}^{\prime }\right]\hfill \\ & & +4\overline{K_{ij}{K}^{ab}}\left[K_{ab}\right]-\left[K_{ij}\right]\overline{R}+\overline{K_{ij}}[K^{ab}{K}_{ab}+6K^{\prime }])=8\pi S_{ij},\hfill \end{array}$$

(86)

$$\begin{array}{ccc}\hfill & & 2{\beta }_1{\nabla }_i[2K^{\prime }+K^{ab}{K}_{ab}]=8\pi S_{in},\hfill \end{array}$$

(87)

$$\begin{array}{ccc}\hfill & & -2{\beta }_1\left[2K{K}^{\prime }+K{K}^{ab}{K}_{ab}\right]=8\pi S_{nn}.\hfill \end{array}$$

(88)

Before we briefly comment on some typical features of the junction conditions for a general $\mathcal{F}\left(R\right)$ theory, let us review in detail some aspects of the junction equations in the $R^2$ theory which would serve to highlight certain useful points. The equations of motion for the $R^2$ theory read

$$G_{\mu \nu }+{\beta }_1\left(2R{R}_{\mu \nu }-\frac{1}{2}{R}^2{g}_{\mu \nu }+2(g_{\mu \nu }-{\nabla }_{\mu }{\nabla }_{\nu })R\right)=8\pi T_{\mu \nu }.$$

(89)

Recall that $R=\widehat{R}-2{\partial }_{n}K-K^2-K^{ab}{K}_{ab}$. Setting $\left[K\right]=0$ thus ensures that R is non-singular. From inspection, there is no singular term arising from integrating (89) across $\Sigma$. The $R^2{g}_{\mu \nu }$ term contains no delta-singularity and hence it integrates to zero. Consider the term $2g_{\mu \nu }R$, and take the indices to be parallel to $\Sigma$. The terms which could be non-vanishing after integrating are $g_{ij}{\partial }_n^{2}R+g_{ij}{g}^{kl}{\Gamma }_{kl}^n{\partial }_{n}R=g_{ij}{\partial }_n^{2}R+g_{ij}K{\partial }_{n}R={\partial }_n\left(g_{ij}{\partial }_{n}R\right)-{\partial }_n{g}_{ij}{\partial }_{n}R+g_{ij}K{\partial }_{n}R$. Integrating ${\partial }_n\left(g_{ij}{\partial }_{n}R\right)$ across $\Sigma$ yields $\left[g_{ij}{\partial }_{n}R\right]$. The other two terms require some work. As we discussed in Section 3, we can proceed by expanding all terms using the Gauss–Codazzi relations and then integrating the terms ensuring that the various discontinuities present in each term and its derivatives are handled appropriately (e.g., using the integral identities in Appendix A). We find

$$\begin{array}{ccc}\hfill -\int dn{\partial }_n{g}_{ij}{\partial }_{n}R& =& 2\overline{K_{ij}}\left(2\left[K^{\prime }\right]+\left[K^{ab}{K}_{ab}\right]\right)+\frac{1}{3}\left[K_{ij}\right]\left[K^{ab}\right]\left[K_{ab}\right]=-2\overline{K_{ij}}\left[R\right]+\frac{1}{3}\left[K_{ij}\right]\left[K^{ab}\right]\left[K_{ab}\right],\hfill \\ \hfill \int dn{g}_{ij}K{\partial }_{n}R& =& -g_{ij}K(2\left[K^{\prime }\right]+\left[K^{ab}{K}_{ab}\right])=g_{ij}K\left[R\right].\hfill \end{array}$$

(90)

The remaining terms in (89) integrate to yield

$$\begin{array}{ccc}\hfill \int dnR{R}_{ij}& =& 2\overline{K^{\prime }}\left[K_{ij}\right]-\left(\widehat{R}-K^2-\overline{K_{ab}{K}^{ab}}\right)\left[K_{ij}\right]=-\overline{R}\left[K_{ij}\right]-\frac{1}{6}\left[K^{ab}\right]\left[K_{ab}\right]\left[K_{ij}\right],\hfill \\ \hfill \int dn{\nabla }_i{\nabla }_{j}R& =& \overline{K_{ij}}\left[R\right]-\frac{1}{6}\left[K^{ab}\right]\left[K_{ab}\right]\left[K_{ij}\right].\hfill \end{array}$$

(91)

Summing up all terms, we find the junction condition

$$-\left[K_{ij}\right]+2{\beta }_1\left(g_{ij}\left[n^{\alpha }{\partial }_{\alpha }R\right]-\overline{R}\left[K_{ij}\right]-3\overline{K_{ij}}\left[R\right]+\left[R\right]K{g}_{ij}+\frac{1}{3}\left[K_{ij}\right]\left[K^{ab}\right]\left[K_{ab}\right]\right)=8\pi S_{ij},$$

(92)

which is equivalent to (86) but with some expressions succinctly expressed in terms of the Ricci scalar. We can also express (87) and (88) as

$$2{\beta }_1{\nabla }_i\left[R\right]=-8\pi S_{in},2{\beta }_{1}K\left[R\right]=8\pi S_{nn}.$$

(93)

Incidentally, in [14,15], the authors attempted to study the junction conditions for $\mathcal{F}\left(R\right)$ gravity and obtained different results in general. In [15], the junction conditions were derived under the assumption of a continuous Ricci scalar.

To facilitate comparison of our results against those of [14,15], consider the case where we impose R to be continuous at $\Sigma$. This renders the junction conditions in directions not parallel to $\Sigma$ trivial, whereas (86) or (92) reduce to

$$2{\beta }_1\left(g_{ij}\left[n^{\alpha }{\partial }_{\alpha }R\right]-R\left[K_{ij}\right]+\frac{1}{3}\left[K_{ij}\right]\left[K^{ab}\right]\left[K_{ab}\right]\right)=8\pi S_{ij}.$$

(94)

The first two terms on the LHS of (94) are naively expected from integrating $2g_{ij}{\partial }_n^{2}R$ and $2R{R}_{ij}$, respectively, after using Gauss–Codazzi relations, and were also obtained in [14] (see Equation (A17) in the Appendix B) and [15] (see Equation (3.11)). However, the third term requires a more careful evaluation of the integrals, as we have pinpointed in each of the integrals in (90) and (91).

We note that the appearance of the “anomalous” term $\left[K_{ij}\right]\left[K^{ab}\right]\left[K_{ab}\right]$ in (90) and (91) has been missed in [14,15]. In [15], the authors also employed the technique of integrating the equations of motion across $\Sigma$ but they missed this term due to essentially an incorrect ansatz for the Ricci scalar, as we shall shortly elaborate upon. In [14], a different method for deriving the junction conditions was proposed which, as mentioned earlier, appears to be afflicted with not having a consistent Gauss–Bonnet limit when it is applied to quadratic gravity.

A common feature that characterizes both [14,15] is an improper treatment of ${\Theta }^2\left(n\right)$ term that is present in the Ricci scalar when the extrinsic curvature is discontinuous. As a concrete illustration of its consequence on junction terms, let us revisit (91), this time sparing more details. Acting on the Ricci scalar with covariant derivatives, we have

$${\nabla }_i{\nabla }_{j}R=K_{ij}{\partial }_{n}R+{\partial }_i{\partial }_{j}R-{\Gamma }_{ij}^k{\partial }_{k}R.$$

(95)

Recall (2) that

$$R=R^{+}\Theta +R^{-}(1-\Theta )-\left[K_{ab}\right]\left[K^{ab}\right]({\Theta }^2-\Theta ),$$

(96)

where $\Theta \left(n\right)=\frac{1}{2}+{\Theta }_X\left(n\right)$ is the step function adopted in [14,19]. This implies that only the first term on the RHS of (95) is singular with

$${\partial }_{n}R=\left[R\right]\delta -\left[K_{ab}\right]\left[K^{ab}\right](2\Theta -1)\delta +\dots .$$

(97)

Integrating across the surface, we have

$$\int dn{\nabla }_i{\nabla }_{j}R=\int dn{K}_{ij}\left[R\right]\delta -\int dn\left[K_{ab}\right]\left[K^{ab}\right]K_{ij}(2\Theta -1)\delta .$$

(98)

Substituting

$$K_{ij}=K_{ij}^{-}+\left[K_{ij}\right]\Theta ,\int dn(2\Theta -1)\delta =0,\int dn\Theta (2\Theta -1)\delta =\frac{1}{6},$$

(99)

into (98), we finally have

$$\int dn{\nabla }_i{\nabla }_{j}R=\overline{K_{ij}}\left[R\right]-\frac{1}{6}\left[K_{ij}\right]\left[K_{ab}\right]\left[K^{ab}\right],$$

(100)

which is, of course, (91). It should be clear from this explicit computation that one will miss the term $-\frac{1}{6}\left[K_{ij}\right]\left[K_{ab}\right]\left[K^{ab}\right]$ if we have taken ${\Theta }^2-\Theta =0$ right from the outset erroneously. This error would then lead to the absence of the second term in (98), and only the first term would then remain in (100), as was the case in [14,15]. Fundamentally, our use of nascent delta functions enables us to compute the integral of products of the step function, delta function, and their derivatives in a well-defined and unambiguous manner, as we discussed in a more generic fashion earlier in Section 2.4.

4.1.2. Regularity Constraints and Junction Conditions for $R^3$ Theory

In the following, we derive the junction conditions explicitly for the gravitational theory with additional $\beta R^3$ term in the action with the equation of motion

$${\tilde{G}}_{\mu \nu }\equiv G_{\mu \nu }+\beta \left(3R^2{R}_{\mu \nu }-\frac{1}{2}{R}^3{g}_{\mu \nu }+3(g_{\mu \nu }-{\nabla }_{\mu }{\nabla }_{\nu })R^2\right)=8\pi T_{\mu \nu },$$

(101)

where $g_{\mu \nu }-{\nabla }_{\mu }{\nabla }_{\nu }=g_{\mu \nu }\left({\partial }_n^2+K{\partial }_n+D^2\right)-\left(K_{\mu \nu }{\partial }_n+D_{\mu }{D}_{\nu }\right)$, and $D_{\alpha }$ is the covariant derivative associated with $\Sigma$. We impose the constraint that R be non-singular, thus taking

$$\left[K\right]=0,$$

as our regularity constraint. From inspection, since the highest order n-derivative is a ${\partial }_n^2$ acting on R which would yield at most a delta-like singularity, this suffices to be our only regularity constraint.

The results obtained by integrating across $\Sigma$ for the various terms in ${\tilde{G}}_{ij}, {\tilde{G}}_{in}, {\tilde{G}}_{nn}$ are collected in the Appendix D. We obtain, for the $\mathcal{F}\left(R\right)=R+\beta R^3$ theory with the regularity constraint $\left[K\right]=0$, the junction conditions

$$-\left[K_{ij}\right]+\beta \left(6\left[g_{ij}R{R}^{\prime }\right]-3{\overline{R}}^2\left[K_{ij}\right]+\left[K_{ij}\right]R_{\left(2\right)}(-2\overline{R}+\frac{1}{5}{R}_{\left(2\right)})+6c_1\left[R\right]\overline{R}-\frac{7}{4}\left[K_{ij}\right]{\left[R\right]}^2+{\mathcal{T}}_{rep}\right)=8\pi S_{ij},$$

(102)

and

$$-6\beta \left[R{\partial }_{i}R\right]=8\pi S_{in},6\beta \left[R\right]K\overline{R}=8\pi S_{nn}.$$

(103)

Here in (102), $c_1:=g_{ij}K-3{\overline{K}}_{ij}$, $R_{\left(2\right)}:=-\left[K^{ab}\right]\left[K_{ab}\right]$ and

$${\mathcal{T}}_{rep}\equiv 24\left[K^{\prime }\right]\left[K_{ij}\right]\left(\left[K_{ab}{K}^{ab}\right](\int dn{{\Theta }^{\prime }}^2{\Theta }_{X}n)-2\left[K^{\prime }\right](\int dn{\Theta }^{″}{\Theta }^{\prime }{\Theta }_X{n}^2)\right).$$

(104)

4.1.3. Some General Comments

Having obtained the explicit junction equations for the $R^2$ and $R^3$ theory via a rather intricate integration procedure, let us briefly comments on some features that we can deduce for the general $\mathcal{F}\left(R\right)$ case without alluding to some specific form of $\mathcal{F}$.

• $\left[\mathit{K}\right]=\mathbf{0}$ as the regularity constraint
  We have seen that $\left[K\right]=0$ is an appropriate regularity constraint for the $R^2$ and $R^3$ theories, and it is easy to see that it is valid for a general analytic $\mathcal{F}$. If $\left[K\right]=0$, then R is non-singular, so ${\partial }_R\mathcal{F}{R}_{\mu \nu }$ cannot generate divergent terms upon integration, since there is at most a delta-singularity carried by $R_{\mu \nu }$. The term $\mathcal{F}{g}_{\mu \nu }$ vanishes after integration, leaving us with only $(g_{\mu \nu }-{\nabla }_{\mu }{\nabla }_{\nu }){\partial }_R\mathcal{F}$. Consider the normal derivatives—after an integration by parts, we are left with only ${\partial }_n$ acting on ${\partial }_R\mathcal{F}$ which generate, at most, a delta-singularity that yields finite quantities after integration. This applies to the junction equations in all directions. Hence, $\left[K\right]=0$ is a valid regularity constraint generally.
  When we discussed the case of quadratic gravity earlier, we obtained the general equations of which solutions give all possible regularity constraints. For the $R^2$ theory, we found no regularity constraint apart from $\left[K\right]=0$, and similarly, one can show that this is the case for the $R^3$ theory as well.
• On taking $\mathit{R}$ to be continuous
  Another constraint that we can impose on top of $\left[K\right]=0$ is the continuity of R at $n=0$. This implies that we take

  $$2\left[K^{\prime }\right]+\left[K^{ab}{K}_{ab}\right]=0,\left[K\right]=0.$$

  (105)
  As we observed earlier, the junction equations in the orthogonal directions would be trivial in this limit, leaving only those parallel to $\Sigma$. In the specific cases of $R^2$ and $R^3$ theories, this is evident in Equations (93), (A87) and (A90). Thus, if we further take $\left[R\right]=0$, the discontinuity in the extrinsic curvature can be physically supported purely by a singular source that only has non-vanishing components parallel to $\Sigma$. For the $R^2$ theory, the junction equation in the case of a continuous R can be found in (94). For the $R^3$ theory, we have explicitly

  $$\begin{array}{ccc}\hfill & & \beta (6g_{ij}R\left[R^{\prime }\right]-3R^2\left[K_{ij}\right]+\left[K_{ij}\right]R_{\left(2\right)}(-2R+\frac{1}{5}{R}_{\left(2\right)})\hfill \\ & & +12\left[K^{\prime }\right]\left[K_{ij}\right]R_{\left(2\right)}\int dn{\left({\Theta }^{\prime }\right)}^3{n}^2)-\left[K_{ij}\right]=8\pi S_{ij}\hfill \end{array}$$

  (106)
  In both (94) and (106), we see that although the first two terms on the LHS may be naively expected from integrating $g_{ij}{\partial }_n^2{\partial }_R\mathcal{F}$ and $\left({\partial }_R\mathcal{F}\right)R_{ij}$, respectively, there are non-trivial terms which arise from the intricacies of the integral.
• On representation-dependent terms and $\left[{\mathit{K}}^{\prime }\right]$
  In (102), we note the appearance of ${\mathcal{T}}_{rep}$ which is sensitive to the choice of the nascent delta function. We note that such terms always come with at least a factor of $\left[K^{\prime }\right]$. To see this, from (A72) and (A73), we see that the ${\Theta }^{\prime },{\Theta }^{″}$ terms in $R,{\partial }_{n}R$ are each multiplied to a factor of $\left[K\right]$. Now, the representation-dependent terms can always be traced to a product of them. Let m be the total order of derivatives defined as the sum of the order of derivative on each ${\Theta }_X$ in some product. Any such term is generically representation-dependent, with the junction term arising from integrating them against the $(m-1)\mathrm{th}$ derivative of the coefficient. These derivatives must act only on factors of $\left[K\right]$, as otherwise, the term will vanish since we have imposed $\left[K\right]=0$ at $n=0$. Thus, we always have $\left[K^{\prime }\right]$ as part of the overall coefficient of any representation-dependent term. Although the actual form of ${\mathcal{T}}_{rep}$ depends on what $\mathcal{F}\left(R\right)$ is as a function of R, this implies that universally across $\mathcal{F}\left(R\right)$-theories, setting

  $$\left[K^{\prime }\right]=0,$$

  (107)

  implies the absence of these representation-dependent terms. We stress, however, that such a constraint is not necessary for the regularity of the junction equations.
• A restrictive set of regularity constraints for a generic $\mathcal{F}\left(\mathit{R}\right)$ theory
  We consider a set of regularity constraints at $n=0$ which allows us to explicitly derive the appropriate junction equations for a generic analytic $\mathcal{F}\left(R\right)$. This amounts to simply taking all components of R to be continuous at $\Sigma$. Since $R=\widehat{R}-2{\partial }_{n}K-K^2-K^{ab}{K}_{ab}$, this implies that we impose

  $$\left[K_{ij}\right]=0,\left[K^{\prime }\right]=0.$$

  (108)
  From (105) and (107), this set of regularity constraints can also be understood as the one that leads to an absence of representation-dependent terms and a singular energy source that does not have components orthogonal to $\Sigma$. Keeping only terms involving normal derivatives, the equations of motion simplify to read

  $$\begin{array}{ccc}& & -\left({\partial }_R\mathcal{F}\right){\partial }_n{K}_{ij}+g_{ij}{\partial }_n\left({\partial }_{n}R\left({\partial }_R^2\mathcal{F}\right)\right)+g_{ij}{R}^{\prime }K\left({\partial }_R^2\mathcal{F}\right)+K_{ij}{R}^{\prime }\left({\partial }_R^2\mathcal{F}\right)+\left(\dots \right)=8\pi T_{ij},\hfill \\ & & -{\partial }_R^2\mathcal{F}{\partial }_n{\partial }_{i}R-{\partial }_{n}R{\partial }_{i}R{\partial }_R^3\mathcal{F}-2\left({\partial }_R^2\mathcal{F}\right)K_i^m{\partial }_m{\partial }_{n}K+\left(\dots \right)=8\pi T_{in},\hfill \\ & & -\left({\partial }_R\mathcal{F}\right){\partial }_{n}K+K{\partial }_{n}R\left({\partial }_R^2\mathcal{F}\right)+\left(\dots \right)=8\pi T_{nn},\hfill \end{array}$$

  (109)

  where we use the Gaussian normal chart in which ${\Gamma }_{ij}^n=-K_{ij},{\Gamma }_{nj}^i=K_j^i$, arethe Gauss–Codazzi relations and note that the components ${\tilde{G}}_{ni}=0$ identically. Integrating across $\Sigma$ in the double scaling limit then yields the junction conditions to be

  $$g_{ij}\left({\partial }_R^2\mathcal{F}\right)\left(2\left[K^{″}\right]+{\left(K^{ab}{K}_{ab}\right)}^{\prime }\right)=-8\pi S_{ij},$$

  (110)

  with no other junction conditions arising from other components ${\tilde{G}}_{\alpha \beta }$. Restoring covariance, we note that $\left[K^{″}\right]=n^{\alpha }{n}^{\beta }\left[{\nabla }_{\alpha }{\nabla }_{\beta }K\right]$, and similarly ${\left(K^{ab}{K}_{ab}\right)}^{\prime }=n^{\alpha }{\nabla }_{\alpha }\left(K^{ab}{K}_{ab}\right)$.
  Again, we wish to emphasize that the smoothness conditions $\left[K_{ij}\right]=0, \left[{\partial }_{n}K\right]=0$ that arise from a continuous Ricci scalar are not the least restrictive ones. As we have seen in the previous examples of $R^2,R^3$ theories, one could just impose $\left[K\right]=0$, leading to a more complicated set of junction conditions. In particular, there are non-trivial ones which generically require the singular source at $\Sigma$ to have non-vanishing orthogonal components.

4.2. Low-Energy Effective Action from Toroidal Compactification of the Heterotic String

As an illustration of how our method could be applied straightforwardly to the presence of matter couplings in higher-derivative gravitational theories, we consider a simple example motivated by string phenomenology—the low-energy effective action arising from a particular compactification of a ten-dimensional string theory. This effective action involves two scalar fields coupled to the Riemann tensor in a certain manner.

To first-order in ${\alpha }^{\prime }$ expansion, suppressing all gauge fields for simplicity, the 10D effective action of the heterotic superstring reads (see, e.g., Chapter 12 of [39])

$$S=\frac{g_s^2}{{\kappa }^{\left(10\right)}}\int d^{10}x\sqrt{\left|g\right|}{e}^{-2\varphi }\left(R+4{\partial }^{\mu }\varphi {\partial }_{\mu }\varphi -\frac{1}{12}{H}^2+\frac{{\alpha }^{\prime }}{8}{R}_{\mu \nu \rho \sigma }{R}^{\mu \nu \rho \sigma }\right),$$

(111)

where ${\kappa }^{\left(10\right)}=16\pi G^{\left(10\right)}$, $g_s$ is the string coupling and the Riemann-squared term is required for supersymmetry to first order in ${\alpha }^{\prime }$ (essentially this follows from the Chern–Simons terms in the three-form field strength). In various generalizations and applications, one often adds to (111) terms that correspond to solitonic branes of lower dimensions, potentials for the dilaton that arise from spontaneous symmetry-breaking, etc. In the following, for simplicity and definiteness, we furnish a derivation of the junction conditions for this theory following a toroidal compactification to 4D, in the absence of source terms for the scalar fields.

It was shown in [40] that upon compactification on a $T^6$, up to leading order in the string coupling, the effective $4D$ action can be simplified to read

$$S_{eff}=\frac{1}{\kappa }\int d^{4}x\sqrt{\left|g\right|}\left(R-\frac{1}{2}\left({\partial }^{\mu }\varphi {\partial }_{\mu }\varphi +{\partial }^{\mu }\phi {\partial }_{\mu }\phi \right)-\frac{{\alpha }^{\prime }}{8}\varphi \left(R^2-4R^{\mu \nu }{R}_{\mu \nu }+R_{\mu \nu \alpha \beta }{R}^{\mu \nu \alpha \beta }\right)+\frac{{\alpha }^{\prime }}{8}\phi R_{\mu \nu \rho \sigma }{\tilde{R}}^{\mu \nu \rho \sigma }\right),$$

(112)

with $\kappa \sim {\kappa }^{\left(10\right)}/\mathrm{Vol}\left(T^6\right)$, ${\tilde{R}}_{\mu \nu \alpha \beta }=\frac{1}{2}{ϵ}_{\mu \nu \rho \sigma }{R}_{\alpha \beta }^{\rho \sigma }{R}^{\mu \nu \alpha \beta }$, and $\varphi ,\phi$ can be interpreted as the dilaton and axion fields, respectively. We see that the effective action is a sum of the Gauss–Bonnet and Chern–Simons terms, each coupled to a scalar field. Since our method works for any coefficients of the interaction terms, in the following, we derive the junction conditions for the above theory with the factors of $\pm \frac{{\alpha }^{\prime }}{8}$ being replaced by ${\beta }_{2,1}$, respectively, which are arbitrary constants in units of ${\alpha }^{\prime }$.

The equations of motion read

$$\begin{array}{ccc}\hfill & & R_{\mu \nu }-\frac{1}{2}R{g}_{\mu \nu }-{\beta }_1{g}_{\nu \lambda }{\delta }_{\mu \rho \gamma \delta }^{\lambda \sigma \alpha \beta }{R^{\gamma \delta }}_{\alpha \beta }{\nabla }^{\rho }{\nabla }_{\sigma }\varphi +2{\beta }_2{\nabla }^{\rho }{\nabla }^{\sigma }\left({\tilde{R}}_{\rho \left(\mu \nu \right)\sigma }\phi \right)+{\beta }_1\varphi H_{\mu \nu }\hfill \\ & & +\frac{1}{2}\left({\partial }_{\mu }\varphi {\partial }_{\nu }\varphi -\frac{1}{2}{g}_{\mu \nu }{(\partial \varphi )}^2\right)+\frac{1}{2}\left({\partial }_{\mu }\phi {\partial }_{\nu }\phi -\frac{1}{2}{g}_{\mu \nu }{(\partial \phi )}^2\right)=8\pi T_{\mu \nu },\hfill \end{array}$$

(113)

$$\begin{array}{ccc}\hfill & & {\nabla }^2\varphi =-{\beta }_1\left(R^2-4R^{\mu \nu }{R}_{\mu \nu }+R_{\mu \nu \alpha \beta }{R}^{\mu \nu \alpha \beta }\right),\hfill \end{array}$$

(114)

$$\begin{array}{ccc}\hfill & & {\nabla }^2\phi =-{\beta }_2{R}_{\mu \nu \rho \sigma }{\tilde{R}}^{\mu \nu \rho \sigma },\hfill \end{array}$$

(115)

where $T_{\mu \nu }$ represents the external energy-momentum tensor field, and

$$H_{\mu \nu }=2\left(R{R}_{\mu \nu }-2R_{\alpha \mu }{R}_{\nu }^{\alpha }+R_{\mu \alpha \beta \gamma }{R_{\nu }}^{\alpha \beta \gamma }-2R_{\mu \alpha \nu \beta }{R}^{\alpha \beta }\right)-\frac{1}{2}{g}_{\mu \nu }\left(R^2-4R_{\alpha \beta }{R}^{\alpha \beta }+R^{\rho \sigma \alpha \beta }{R}_{\rho \sigma \alpha \beta }\right)$$

is the second Lovelock tensor. For simplicity, we assume that the thin-shell only sources the metric and not the dilaton or axion field. When this assumption is relaxed in more realistic generalizations of (111), we expect the presence of other source terms on the RHS of (114) and (115)12. For the term in ${\beta }_2$, it is useful to invoke Bianchi identities to rewrite it as

$${\beta }_2{\nabla }^{\rho }{\nabla }^{\sigma }\left({\tilde{R}}_{\rho \left(nn\right)\sigma }\phi \right)={\beta }_2{\nabla }_k\left({\partial }^m\phi {ϵ}_{mefn}{R_n}^{kef}\right).$$

(116)

before proceeding to derive the junction equations below.

4.2.1. Regularity Constraints

We first examine the singular terms in (113)–(115). Let us begin with the matter fields’ equations of motion. Since there is no term more singular than the delta function in the Gauss–Bonnet nor the Chern–Simons term, this implies that integrating (114), (115) across $\Sigma$ will yield a finite quantity. In the Gaussian normal chart, ${\nabla }^2={\partial }_n^2+{\nabla }_{\Sigma }^2$. Hence, we have the continuity constraints for the scalar fields:

$$\left[\varphi \right]=0,\left[\phi \right]=0,$$

(117)

since otherwise, the action of ${\partial }_n^2$ on the fields would lead to a singular integral across $\Sigma$. To preserve as much generality as possible, we do not assume, however, that their normal derivatives are continuous at $\Sigma$. We can integrate (114) and (115) across $\Sigma$ and obtain

$$\left[{\varphi }^{\prime }\right]=-4\widehat{R}\left[K\right]+8{\widehat{R}}^{ab}\left[K_{ab}\right]-4\left[K{K}^{ab}{K}_{ab}\right]+\frac{4}{3}\left[K^3\right]+\frac{8}{3}\left[K_l^m{K}^{ln}{K}_{mn}\right],$$

(118)

$$\left[{\phi }^{\prime }\right]=8{\beta }_2{ϵ}_{ijk}{D}^j\overline{K^{kb}}{K}_b^i.$$

(119)

Let us now study if (113) yields any non-trivial regularity constraints. Consider first ${\tilde{G}}_{ij}$, suppressing the Gauss–Bonnet and scalar field terms,

$${\tilde{G}}_{ij}=-{\beta }_1{g}_{j\lambda }{\delta }_{i\rho \gamma \delta }^{\lambda \sigma \alpha \beta }{R}_{\alpha \beta }^{\gamma \delta }{\nabla }^{\rho }{\nabla }_{\sigma }\varphi +{\beta }_2{\nabla }_k\left({\partial }^m\phi {ϵ}_{mef(j}{R_{i)}}^{kef}\right)+\dots ,$$

(120)

where the ellipses refer to finite terms in (9). For a $1/b$-type singularity, we need terms of the form $R_{nb}^{na}$ and ${\partial }_n^2\varphi$ to be present simultaneously or a term that goes as ${\nabla }_n{R}_{nb}^{na}$. By observation, this cannot arise from the dilaton term where, for the axion interaction term, we consider the term where we take the indices $k=f=n$ in (120). After some algebra, we can simplify this term to read $-2{\beta }_2{\partial }^m\phi {ϵ}_{me(j}{\partial }_n^2{K}_{i)}^e$, but integrating this term across $\Sigma$ yields

$$-{{ϵ}^{me}}_j\left[{\partial }_m\phi \right]\left[K_{ei}\right]\underset{b\to 0}{lim}\frac{1}{b}\int dX{F}^2\left(X\right)+\dots$$

where we suppress the manifestly finite quantities. Since $\phi$ and, thus, ${\partial }_m\phi$ are continuous across $\Sigma$, we find no singularity here. Finally, we note that, as shown in the previous section, there is no singular term descending from $H_{\mu \nu }$ (no additional regularity constraint required for pure Einstein–Gauss–Bonnet theory). Thus, from inspection, we can see that there is no singular term from integrating ${\tilde{G}}_{ij}$ across $\Sigma$.

Similar arguments apply for the other components of the field equations, explicitly,

$$\begin{array}{ccc}\hfill {\tilde{G}}_{nn}& =& -{\beta }_1{g}_{nn}{\delta }_{n\rho \gamma \delta }^{n\sigma \alpha \beta }{R}_{\alpha \beta }^{\gamma \delta }{\nabla }^{\rho }{\nabla }_{\sigma }\varphi +2{\beta }_2{\nabla }_k\left({\partial }^m\phi {ϵ}_{mef(n}{R_{n)}}^{kef}\right)+\dots \hfill \\ \hfill \end{array}$$

(121)

$$\begin{array}{ccc}\hfill {\tilde{G}}_{in}& =& -{\beta }_1{g}_{ik}{\delta }_{n\rho \gamma \delta }^{k\sigma \alpha \beta }{R}_{\alpha \beta }^{\gamma \delta }{\nabla }^{\rho }{\nabla }_{\sigma }\varphi +{\beta }_2{\nabla }_k\left({\partial }^m\phi {ϵ}_{mef(i}{R_{n)}}^{kef}\right)+\dots \hfill \end{array}$$

(122)

By inspection, one can again draw the conclusion that no singular terms remain after integrating (121) and (122) across $\Sigma$ by checking the absence of terms of the form ${\nabla }_n{R}_{nb}^{na}$ or products of $R_{nb}^{na}$ and ${\partial }_n^2\varphi$. Altogether, the above considerations reveal that, similar to the pure Gauss–Bonnet theory, there is no additional regularity constraint that we need to impose here. This fits intuitively well with the fact that the interaction terms mixing the otherwise free scalar fields and the graviton happen to be linear in the Gauss–Bonnet and Chern–Simons terms.

4.2.2. Junction Conditions

We now proceed to derive the junction equations. From the form of (121), we find that there are no terms which carry delta-like singularities, so integrating (121) across $\Sigma$ cannot give any junction condition. For (122), the dilaton interaction term could not give rise to any junction term since the generalized Kronecker delta symbol already contains an n-index, and introducing some pair of n-indices (for either $R_{nb}^{na}$ or ${\partial }_n^2\varphi$) would annihilate the term by symmetry of the symbol.

For the axion interaction term in (122), which reads

$$-{\beta }_2\left({\nabla }_a\left({\partial }^{\beta }\phi {ϵ}_{\beta \mu \nu i}{R_n}^{a\mu \nu }\right)+{\nabla }_{\alpha }\left({\partial }^b\phi {ϵ}_{buvn}{R_i}^{\alpha uv}\right)\right).$$

(123)

For the second term in (123), the only non-vanishing term arises from taking the dummy index $\alpha$ to be n, generating the following junction term after integration:

$${\beta }_2{ϵ}_{buvn}{\partial }^b\phi \left[{R_i}^{nuv}\right].$$

For the first term in (123), it is helpful to first expand all terms that come with the covariant derivative.

$$\begin{array}{ccc}& & {\partial }^{\beta }\phi {ϵ}_{\beta \mu \nu i}\left({\partial }_a{R_n}^{a\mu \nu }+{\Gamma }_{al}^a{R_n}^{l\mu \nu }+{\Gamma }_{ak}^{\mu }{R_n}^{ak\nu }+{\Gamma }_{ak}^{\nu }{R_n}^{a\mu k}+{\Gamma }_{an}^{\mu }{R_n}^{an\nu }+{\Gamma }_{an}^{\nu }{R_n}^{a\mu n}\right)\hfill \\ & =& 2{\partial }^{\beta }\phi {ϵ}_{\beta nvi}\left({\partial }_a{R_n}^{anv}+{\Gamma }_{al}^a{R_n}^{lnv}\right)+2{ϵ}_{nuvi}{\phi }^{\prime }{K}_a^u{R_n}^{anv}+2{ϵ}_{buni}{\partial }^b\phi {\Gamma }_{ak}^u{R_n}^{akn}+\left(\dots \right),\hfill \end{array}$$

(124)

where we have only displayed terms that will survive the integral. Invoking the Gauss–Codazzi relation ${R^n}_{snv}=-{\partial }_n{K}_{sv}+K_{bv}{K}_s^b$, and assembling all terms, we integrate to obtain the junction condition for ${\tilde{G}}_{in}$ to be

$$8\pi S_{in}={\beta }_2\left({ϵ}_{buvn}{\partial }^b\phi \left[{R_i}^{nuv}\right]+2{ϵ}_{nmri}\left({\partial }^m\phi \left[D_a{K}^{ar}\right]-\overline{{\phi }^{\prime }{K}_a^m}\left[K^{ar}\right]\right)\right).$$

(125)

We now consider the junction condition for ${\tilde{G}}_{ij}$. Apart from the Gauss–Bonnet term (multiplied to $\varphi$) which integrates to yield $\varphi J_{ij}^{\left(GB\right)}$, where $J_{ij}^{\left(GB\right)}$ is the junction equation of the pure Gauss–Bonnet theory (see (75)), the dilaton and axion interaction terms contribute to the junction condition as well.

Let us begin with the dilaton interaction term which could contribute to the junction equation with the following two terms from (120):

$${\tilde{G}}_{ij}=-{\beta }_1{g}_{jl}{\delta }_{ingd}^{lnab}{R}_{ab}^{gd}{\nabla }^n{\nabla }_n\varphi -{\beta }_1{g}_{jl}{\delta }_{irnd}^{lsnb}{R}_{nb}^{nd}{\nabla }^r{\nabla }_s\varphi +\left(\dots \right)$$

(126)

For the first term, expanding the generalized Kronecker symbol, we obtain

$${\beta }_1\left(g_{ji}{R}_{gd}^{gd}-g_{jg}{R}_{id}^{gd}+g_{jd}{R}_{ig}^{gd}-g_{jd}{R}_{gi}^{gd}+g_{jg}{R}_{di}^{gd}-g_{ij}{R}_{dg}^{gd}\right){\nabla }^n{\nabla }_n\varphi =-4{\beta }_1{G}_{ij}{\nabla }^n{\nabla }_n\varphi$$

(127)

Integrating across $\Sigma$ then gives

$$-4{\beta }_1\overline{G_{ij}}\left[{\partial }_n\varphi \right]\equiv -4{\beta }_1\left(-\overline{K_{ij}^{\prime }}+2\overline{K_j^a{K}_{ai}}-\overline{K{K}_{ij}}+{\widehat{R}}_{ij}-\frac{1}{2}{g}_{ij}\left(\widehat{R}-2\overline{K^{\prime }}-\overline{K^2}-\overline{K^{ab}{K}_{ab}}\right)\right)\left[{\partial }_n\varphi \right]$$

(128)

For the second term in (126), we first note that

$$g_{jl}{\delta }_{igd}^{lab}{\nabla }^g{\nabla }_a\varphi K_b^d=g_{ji}{\nabla }_{\Sigma }^2\varphi K-{\nabla }_j{\nabla }_i\varphi K+{\nabla }^m{\nabla }_i\varphi K_{jm}-{\nabla }_{\Sigma }^2\varphi K_{ji}+{\nabla }_j{\nabla }_d\varphi K_i^d-g_{ji}{\nabla }^g{\nabla }_d\varphi K_g^d$$

(129)

after expanding the generalized Kronecker symbol. Integrating the second term of (126) across $\Sigma$ then gives, after some algebra,

$${\beta }_1\left({\nabla }_{\Sigma }^2\varphi (g_{ij}\left[K\right]-\left[K_{ij}\right])+(D_m{D}_a\varphi +\overline{K_{ma}{\varphi }^{\prime }}){\mathcal{C}}_{ij}^{mars}\left[K_{rs}\right]\right),$$

(130)

where we have defined

$${\mathcal{C}}_{ij}^{mars}\equiv {\delta }_i^a{\delta }_j^r{g}^{ms}+{\delta }_j^m{\delta }_i^s{g}^{ar}-{\delta }_j^m{\delta }_i^a{g}^{rs}-g^{ar}{g}^{ms}{g}_{ij}.$$

For the axion interaction term, it is helpful to first write out explicitly three types of terms which contribute to the junction equation. From the second term of (120), we have

$${\beta }_2\left(2{\nabla }_n\left({\partial }^m\phi {ϵ}_{men(j}{R_{i)}}^{nen}\right)+{\nabla }_n\left({\partial }^n\phi {ϵ}_{nef(j}{R_{i)}}^{nef}\right)+2{\nabla }_k\left({\partial }^m\phi {ϵ}_{mnf(j}{R_{i)}}^{knf}\right)\right)+\left(\dots \right)$$

(131)

Expanding the first term of (131) and using Gauss–Codazzi relation to express the Christoffel symbols in terms of the extrinsic curvature gives

$${\nabla }_n\left({\partial }^m\phi {ϵ}_{men(j}{R_{i)}}^{nen}\right)={\partial }_n\left({\partial }^m\phi {ϵ}_{men(j}{R_{i)}}^{nen}\right)+K_l^m{\partial }^l\phi {ϵ}_{men(j}{R_{i)}}^{nen}+{\partial }^m\phi {ϵ}_{men(j}\left({R_{i)}}^{nln}{K}_l^e-K_{i)}^a{R_a}^{nen}\right)$$

(132)

After integration, this gives the junction terms

$$\left[{\partial }^m\phi {ϵ}_{men(j}{R_{i)}}^{nen}\right]+\overline{K_l^m}{\partial }^l\phi {ϵ}_{men(j}\left[K_{i)}^e\right]+\overline{K_l^e}{\partial }^m\phi {ϵ}_{men(j}\left[K_{i)}^l\right]-{\partial }^m\phi {ϵ}_{men(j}\overline{K_{i)}^a}\left[K_a^e\right]$$

(133)

The second term of (131) simply gives, after integration,

$$\left[{\partial }^n\phi {ϵ}_{nef(j}{R_{i)}}^{nef}\right],$$

whereas the third term of (131) yields

$${\partial }^m\phi {ϵ}_{mnf(j}\left(-{\overline{K}}_{i)k}\left[K^{kf}\right]+\overline{K}\left[K_{i)}^f\right]\right),$$

after we note that ${\nabla }_k\left({\partial }^m\phi {ϵ}_{mnf(j}{R_{i)}}^{knf}\right)={\partial }^m\phi {ϵ}_{mnf(j}\left(K_{i)k}{R_n}^{knf}-K{R_{i)}}^{nnf}+\dots \right)$. Assembling all terms together, we find the junction condition associated with ${\tilde{G}}_{ij}$ to be

$$\begin{array}{ccc}\hfill 8\pi S_{ij}& =& {\beta }_2\left(\left[{\partial }^n\phi {ϵ}_{nef(j}{R_{i)}}^{nef}\right]+2{\partial }^m\phi {ϵ}_{mnf(j}\left(-{\overline{K}}_{ki)}\left[K^{kf}\right]+\overline{K}\left[K_{i)}^f\right]\right)\right)\hfill \\ & & +2{\beta }_2\left(\left[{\partial }^m\phi {ϵ}_{men(j}{R_{i)}}^{nen}\right]+{\partial }^l\phi {ϵ}_{men(j}\overline{K_l^m}\left[K_{i)}^e\right]+{\partial }^m\phi {ϵ}_{men(j}\overline{K_l^e}\left[K_{i)}^l\right]\right)\hfill \\ & +& {\beta }_1\left({\nabla }_{\Sigma }^2\varphi (g_{ij}\left[K\right]-\left[K_{ij}\right])+(D_m{D}_a\varphi +\overline{K_{ma}{\varphi }^{\prime }}){\mathcal{C}}_{ij}^{mars}\left[K_{rs}\right]-4\overline{G_{ij}}\left[{\varphi }^{\prime }\right]+\varphi \left[J_{ij}^{\left(GB\right)}\right]\right)\hfill \\ & +& \left[K\right]h_{ij}-\left[K_{ij}\right].\hfill \end{array}$$

(134)

Together with (125) and (134), we specify the appropriate junction condition for the low-energy effective theory defined by the action (112). As a simple consistency check, we note that taking the scalar fields to be constant give the junction conditions for the pure Gauss–Bonnet theory (for a continuous axion, the pure Chern–Simons gravity theory has a trivial junction condition). If we assume a stronger constraint $\left[K_{ij}\right]=0$, then (134) simplifies to read

$$8\pi S_{ij}=-2{\beta }_2{n}_{\alpha }{ϵ}_{mln(j}\left[{\nabla }^{\alpha }{K}_{i)}^l\right]{\partial }^m\phi ,$$

(135)

where only the axion-coupling remains in the junction condition. These junction conditions specify how the matter or energy content of the thin-shell should be constrained in the theory (111). Although we have only treated the simple case where the thin-shell source does not couple to the dilaton/axion, our approach can be adopted straightforwardly to various generalizations that realistically augment (111) with couplings among various fields and potentials for the scalar fields. For example, our derivation here would serve as a useful warm-up for generalizing the junction conditions considered in [41,42,43], where the low-energy effective action of the heterotic string is applied to cosmology.

4.3. Higher-Dimensional Euler Densities

The Gauss–Bonnet term is non-topological beyond four dimensions, and is an example of the Euler characteristic that is the most general extension of the Einstein–Hilbert action that yields at most second-order field equations. The appropriate surface terms were derived some time ago, and their variation with respect to the induced metric of $\Sigma$ gives the junction conditions, as was shown explicitly in the Gauss–Bonnet case in [25,26,29].

The topological Euler density term for a $2m$-dimensional manifold is defined as

$${\mathcal{L}}_m={\Omega }^{a_1{b}_1}\wedge \dots \wedge {\Omega }^{a_m{b}_m}\wedge {ϵ}_{a_1{b}_1\dots a_m{b}_m}=\frac{1}{2^m}{\delta }_{a_1{b}_1\dots a_m{b}_m}^{[c_1{d}_1\dots c_m{d}_m]}{R^{a_1{b}_1}}_{c_1{d}_1}\dots {R^{a_m{b}_m}}_{c_m{d}_m},$$

(136)

where the Kronecker $\delta$-function above is totally antisymmetric in both sets of indices, $\Omega$ is the curvature two-form, and we have normalized it such that the anti-symmetrization symbol in (136) has no other normalization factor. Note that ${\mathcal{L}}_1,{\mathcal{L}}_2$ are the Ricci scalar and Gauss–Bonnet term, respectively. One can consider extending (136) to other dimensions apart from $2m$. For dimensions less than $2m$, it simply vanishes whereas for dimensions higher, it will be non-topological.

Now, the Euler-Lagrange equations of motion of (136) read

$${\tilde{G}}_{ij}=-\frac{1}{2^{m+1}}{g}_{i\mu }{\delta }_{j{b}_1\dots b_{2m}}^{[\mu a_1\dots a_{2m}]}{R_{a_1{a}_2}}^{b_1{b}_2}\dots {R_{a_{2m-1}{a}_{2m}}}^{b_{2m-1}{b}_{2m}}.$$

(137)

This follows from a well-posed variational principle if appropriate surface terms can be added at $\Sigma$ such that setting $\delta g_{ij}=0$ at $\Sigma$ is sufficient for the vanishing of the action variation, and no terms of the form ${\nabla }_k\delta g_{ij}$ survive. In the following, we briefly review how such a surface term was derived in [26] in the language of differential forms. One begins by defining a Chern–Simons form $Q_m$ such that

$$\begin{array}{ccc}\hfill {\mathcal{L}}_m\left(\omega \right)-{\mathcal{L}}_m\left({\omega }_0\right)& =& d{Q}_m(\omega ,{\omega }_0),\hfill \end{array}$$

(138)

$$\begin{array}{ccc}\hfill Q_m& =& m{\int }_0^{1}ds{\theta }^{a_1{b}_1}\wedge {\Omega }_s^{a_2{b}_2}\wedge \dots \wedge {\Omega }_s^{a_m{b}_m}\wedge {ϵ}_{a_1{b}_1\dots a_m{b}_m},\hfill \end{array}$$

(139)

where $\omega$ is the connection one-form, ${\omega }_0$ is the connection one-form defined on $\Sigma$, $\theta =\omega -{\omega }_0$ is the extrinsic curvature/second fundamental form, and finally ${\Omega }_s=d{\omega }_s+{\omega }_s\wedge {\omega }_s$ is the curvature two-form defined with ${\omega }_s=\omega -s\theta$. In [26], it was shown that taking the variation of (138) implies ${\delta }_{\omega }{\int }_M{\mathcal{L}}_m={\int }_{M}d\left({\delta }_{\omega }{Q}_m\right),$ which leads naturally to the surface term in the action

$$I_m={\int }_M{\mathcal{L}}_m-{\int }_{\partial M}{Q}_m.$$

(140)

In coordinate form, the Chern–Simons form that enacts the surface term can be written as [26]

$$Q_m={\int }_0^{1}ds{\delta }_{{\nu }_1\dots {\nu }_{2m-1}}^{{\mu }_1\dots {\mu }_{2m-1}}{K}_{{\mu }_1}^{{\nu }_1}\times \left(\frac{1}{2}{R}_{{\mu }_2{\mu }_3}^{{\nu }_2{\nu }_3}-s^2{K}_{{\mu }_2}^{{\nu }_2}{K}_{{\mu }_3}^{{\nu }_3}\right)\dots \left(\frac{1}{2}{R}_{{\mu }_{2m-2}{\mu }_{2m-1}}^{{\nu }_{2m-2}{\nu }_{2m-1}}-s^2{K}_{{\mu }_{2m-2}}^{{\nu }_{2m-2}}{K}_{{\mu }_{2m-1}}^{{\nu }_{2m-1}}\right).$$

(141)

In principle, one could take the variation of $Q_m$ with respect to the induced metric to obtain the junction conditions. Yet, even for the Gauss–Bonnet case, this can be a rather elaborate calculation, as shown in [25,29].

In the following, we will integrate the equations of motion ${\tilde{G}}_{ij}$ across $\Sigma$ and obtain the junction conditions for the Euler density term valid in dimension $>2 \mathrm{m}$. We first check that there are no singular terms and hence no regularity constraints to impose. This is manifest in the form of (137). Consider again the Gauss–Bonnet term ($m=2$) as an example. The singular terms can only arise from a product of two Riemann tensors, each of which carries two ‘n’ indices. This leads to the antisymmetric delta function of the form

$${\delta }_{jn{\beta }_{2}n{\beta }_4}^{\left[\mu n{\alpha }_{2}n{\alpha }_4\right]},$$

which is identically zero. Similarly for $m>2$, the singular terms arise from delta functions of the form

$${\delta }_{jn{\beta }_{2}n{\beta }_{4}n\dots }^{[\mu n{\alpha }_{2}n{\alpha }_{4}n\dots ]}.$$

The permutations among the “n” indices come in pairs of ± signs and hence they sum to zero. A subtle point is that, as mentioned earlier, in general, there are different classes of divergent terms defined by the number of Riemann tensors with a pair of “n” indices. Each class vanishes separately by the same reason.

Using (137) as the new starting point, we now rederive the junction condition for the Gauss–Bonnet theory which gives us some intuition on how this generalizes for the higher Euler densities. The form and symmetry of the antisymmetric Kronecker $\delta$-function implies that we can pick one of the Riemann tensors to have two “n” indices and thus carry the singular delta function (i.e., the term $\sim -{\partial }_n{K}_{ab}$), whereas the other one should carry all indices parallel to $\Sigma$. Up to some degeneracy factor, we have the junction term

$$\sim -g_{ik}{\delta }_{jcdn{\beta }_4}^{\left[kabn{\alpha }_4\right]}\left({\widehat{R}}_{ab}^{cd}-\overline{K_a^c{K}_b^d}+\overline{K_b^c{K}_a^d}\right)\left[K_{{\alpha }_4}^{{\beta }_4}\right],$$

(142)

where we have used the form of antisymmetric Kronecker $\delta$-function to deduce that Riemann tensor with only one “n” index does not contribute. Note the exchange symmetry between the pairs of indices $(a,c), (b,d)$ and $({\alpha }_4,{\beta }_4)$ and that the term $\left[K_{{\alpha }_4}^{{\beta }_4}\right]$ arises from $-{\partial }_n{K}_{{\alpha }_4}^{{\beta }_4}$. Taking into account the m-dependent coefficient of ${\tilde{G}}_{ij}$ and some symmetry factors, we have the junction term

$$\frac{{(-1)}^2(2\times 4)}{2^{2+1}}{g}_{ik}{\delta }_{jcdn{\beta }_4}^{\left[kabn{\alpha }_4\right]}\left[\left({\widehat{R}}_{ab}^{cd}-\frac{2}{3}{K}_a^c{K}_b^d\right)K_{{\alpha }_4}^{{\beta }_4}\right].$$

(143)

One can show that this is identical to the junction condition we derived earlier in (75). It is apparently in a more compact form due to the choice of expressing the Gauss–Bonnet equations of motion using the antisymmetric Kronecker $\delta$-function, but, more importantly, it generalizes to the higher Euler densities fairly straightforwardly. After some algebra, we find the junction condition for the theory defined by $I_{GR}+\beta I_m, m>2$ to be

$$\begin{array}{ccc}& & \beta \frac{2m}{2^m}{g}_{i\mu }{\delta }_{jn{b}_2{b}_3{b}_4\dots b_{2m}}^{[\mu n{a}_2{a}_3{a}_4\dots a_{2m}]}\sum _{l=0}^{m-1}\frac{{(-2)}^l}{2l+1}\left({\widehat{R}}_{a_3{a}_4}^{b_3{b}_4}\dots {\widehat{R}}_{a_{2(m-l)-1}{a}_{2(m-l)}}^{b_{2(m-l)-1}{b}_{2(m-l)}}{K}_{a_2}^{b_2}{K}_{a_{2(m-l)+1}}^{b_{2(m-l)+1}}\dots K_{a_{2m}}^{b_{2m}}\right)\hfill \\ & & =8\pi S_{ij}+\left[K_{ij}\right]-\left[K\right]h_{ij},\hfill \end{array}$$

(144)

where $\beta$ is the coupling constant for $I_m$ added to the ordinary Einstein–Hilbert term in the action. There are no non-trivial junction conditions in directions non-parallel to $\Sigma$.

5. Applications

5.1. Thin-Shell Wormholes in $R^2$ Gravity

As an application, we examine the energy conditions governing thin-shell wormholes constructed by a cut-and-paste method (similar to the way we define $\Sigma$ as an identification between two manifolds). For definiteness, let us again consider the $R+\beta R^2$ theory which has garnered much interest in recent literature since it was proposed in 1979 by Starobinsky and Gurovich in [44] as a natural model for cosmological inflation. We had earlier derived the junction conditions for such a theory in Section 3 and also in Equations (92) and (93), which we reproduce here for reading convenience.

$$\begin{array}{ccc}\hfill & & -\left[K_{ij}\right]+2\beta \left(g_{ij}\left[{\partial }_{n}R\right]-\overline{R}\left[K_{ij}\right]-\overline{K_{ij}}\left[R\right]+\left[R\right](-2\overline{K_{ij}}+K{g}_{ij})+\frac{1}{3}\left[K_{ij}\right]\left[K^{ab}\right]\left[K_{ab}\right]\right)=8\pi S_{ij},\hfill \\ \hfill & & 2\beta {\nabla }_i\left[R\right]=-8\pi S_{in},2\beta K\left[R\right]=8\pi S_{nn},\hfill \end{array}$$

(145)

together with the condition $\left[K\right]=0$. From the equation of motion (89), we see that any Ricci-flat geometry is valid as a vacuum solution. Thin-shell wormholes have been studied in this theory in a few works, such as [45,46]13, but unfortunately assuming an incorrect set of junction conditions. Consider the following spherically symmetric ansatz

$$d{s}^2=-A\left(r\right)d{t}^2+A{\left(r\right)}^{-1}d{r}^2+r^2(d{\theta }^2+{sin}^2\theta d{\varphi }^2),$$

(146)

with $r>0, \theta \in [0,\pi ], \varphi \in [0,2\pi )$. We can construct a simple model of a thin-shell wormhole by picking some radius a and identifying two copies of the region $r\ge a$ with $\Sigma$ as the hypersurface $r=a$. Such a construction leads to geodesic complete wormhole with $r=a$ being the throat of minimal radius (in terms of the proper radial distance $l={\int }_a^{r}dr{A}^{-1/2}$, the throat is defined by $l=0$). We can also conveniently obtain a family of timelike $\Sigma$ by setting

$$\Sigma :\{r=a\left(t\right)\},n_{\alpha }=(-\dot{a},1,0,0)/\sqrt{A-{\dot{A}}^2/A}.$$

For such a thin-shell wormhole with $\Sigma$ as its throat, the components of the extrinsic curvature read

$$K_{\theta }^{\theta }=K_{\varphi }^{\varphi }=\pm \frac{1}{a}\sqrt{A\left(a\right)+{a^{\prime }}^2},K_t^t=\pm \frac{{\partial }_{r}A\left(a\right)+2a^{″}}{2\sqrt{A\left(a\right)+{a^{\prime }}^2}}$$

(147)

where $a^{\prime }\equiv \frac{da}{d\tau }$, $\tau$ being the proper time on $\Sigma$ of which induced metric reads

$$d{s}^2=-\left(A-\frac{{\dot{a}}^2}{A}\right)d{t}^2+a^{2}d{\Omega }^2=-d{\tau }^2+a^2\left(\tau \right)d{\Omega }^2.$$

(148)

For simplicity, let us now pick A to correspond to a solution with constant Ricci scalar $R_0$ (so, for example, if $R_0=0$, then $A=1-\frac{2M}{r}$). This implies that terms such as $\left[R\right], \left[R^{\prime }\right]$ vanish in the junction equations. For the theory with Lagrangian $R+\beta R^2$, the junction conditions imply the following for the singular source $S_j^i=\mathrm{diag}(-\sigma ,P,P)$.

$$\begin{array}{ccc}\hfill \sigma & =& \frac{1}{2\pi a}\left(1+2\beta R_0-\frac{2\beta }{3}\left[K^{ab}\right]\left[K_{ab}\right]\right)\left[K_t^t\right]=2P,\hfill \end{array}$$

(149)

$$\begin{array}{ccc}\hfill a^{″}& =& -\frac{1}{2}{\partial }_{r}A-\frac{2}{a}\left(A+{a^{\prime }}^2\right),\hfill \end{array}$$

(150)

$$\begin{array}{ccc}\hfill \left[K_t^t\right]& =& -\frac{4\sqrt{A+{a^{\prime }}^2}}{a},\left[K^{ab}\right]\left[K_{ab}\right]=\frac{3}{2}{\left[K_t^t\right]}^2,\hfill \end{array}$$

(151)

where (150) and (151) arise from the regularity constraint $\left[K\right]=0$.

As a simple example, let us consider static solutions with $a\left(\tau \right)=a_0$, with $a_0$ being some positive constant radius parameter. From (150), we have

$$a_0{\partial }_{r}A\left(a_0\right)+4A\left(a_0\right)=0.$$

(152)

If we take $R_0=0$, we are inevitably led to the Schwarzschild ansatz $A\left(r\right)=1-\frac{2M}{r}$ for which (152) implies that $a_0=3M/2$. Since this is unfortunately smaller than the Schwarzschild radius, we cannot construct a typical thin-shell wormhole in this manner with $r=a_0$ as the time-like throat hypersurface of minimal area. Suppose we take $R_0>0$ and, in particular, that it arises from a positive cosmological constant, with $R_0=4\Lambda$, then we are led to the Schwarzschild–de Sitter ansatz with $A\left(r\right)=1-\frac{2M}{r}-\frac{\Lambda r^2}{3}$. The regularity constraint in (152) translates to

$$\Lambda a_0^3-2a_0+3M=0.$$

(153)

To obtain some finite region of positive $g_{tt}$ in the line element, we need to restrict $\Lambda$ to the domain $0<\Lambda M^2<1/9$, where we have two horizons. The cosmological horizon $r_c$ and black hole event horizon $r_h$ can be analytically solved to read

$$r_c/M=2\sqrt{\frac{1}{M^2\Lambda }}cos\left(\frac{1}{3}{cos}^{-1}(-3\sqrt{M^2\Lambda })\right),r_h/M=2\sqrt{\frac{1}{M^2\Lambda }}cos\left(\frac{1}{3}{cos}^{-1}(-3\sqrt{M^2\Lambda })-\frac{2\pi }{3}\right).$$

(154)

On the other hand, solving for $a_0$ in (153) in this domain yields

$$a_0/M=2\sqrt{\frac{2}{3M^2\Lambda }}cos\left(\frac{1}{3}{cos}^{-1}\left(-\frac{9}{4}\sqrt{\frac{3M^2\Lambda }{2}}\right)\right).$$

(155)

We find that (155) falls nicely between the horizons, as depicted in Figure 1 below.

Thus, we see that contrary to the Schwarzschild case, the regularity constraint $\left[K\right]=0$ is compatible with the wormhole construction procedure of identifying exteriors of the Schwarzschild–de Sitter spacetime. For this class of thin-shell wormholes, we find that the energy density $\sigma$ is unfortunately negative definite if we also adopt the unitarity constraint $1+2\beta R_0>0$ (see, e.g., [45]). For the weak energy condition to be obeyed, we require the coupling parameter $\beta$ to satisfy

$$1+\beta \left(2R_0-{\left[K_t^t\right]}^2\right)\le 0.$$

(156)

We find that $2R_0-{\left[K_t^t\right]}^2$ is positive definite for $\Lambda M^2<1/9$, and that there is no negative $\beta$ which satisfies both (156) and the unitarity condition $1+2\beta R_0>0$. This implies that the family of Schwarzschild–de Sitter wormholes constructed here has to be supported by exotic matter at the throat. They are also unstable under radial perturbation 14.

5.2. Implications for Stellar Models

In the absence of a singular source, the generalized junction equations reduce to a set of conditions for the geometry induced by a non-singular energy-momentum tensor that is possibly discontinuous at $\Sigma$. In the ordinary Einstein theory, this simply translates to continuity in the extrinsic curvature, but taking $S_{\mu \nu }=0$ in the generalized junction conditions typically implies more complicated smoothness conditions on the extrinsic curvature.

In this section, we briefly discuss the form of the junction conditions when the source is non-singular, and some implications for the $R+\beta R^2$ theory. In this case, it turns out to be convenient to begin by first setting $S_{in}=0$, which implies that

$$\left[R\right]=R_0,$$

(157)

where $R_0$ is some constant. In the following, we classify the junction conditions according to whether $R_0$ is zero.

(I)

$\left[R\right]=0$:

From the vanishing of $S_{ij}$, we have

$$2\beta g_{ij}{n}^{\alpha }\left[{\nabla }_{\alpha }R\right]+\left[K_{ij}\right]\left(\frac{2\beta }{3}\left(\left[K^{ab}\right]\left[K_{ab}\right]-3R\right)-1\right)=0.$$

Taking the trace implies that $n^{\alpha }\left[{\nabla }_{\alpha }R\right]=0$, together with the junction conditions

$$\left[K_{ij}\right]=0,\mathrm{or}\left[K^{ab}\right]\left[K_{ab}\right]-3R=\frac{3}{2\beta }.$$

(158)

(II)

$\left[R\right]=R_0\ne 0$:

In this case, $S_{nn}=0$ implies that $K=0$ and from the vanishing of $S_{ij}$, we have $n^{\alpha }\left[{\nabla }_{\alpha }R\right]=0$ together with the junction conditions

$$\left[K_{ij}\right]\left(\frac{2\beta }{3}\left(\left[K^{ab}\right]\left[K_{ab}\right]-3\overline{R}\right)-1\right)=6\beta R_0{\overline{K}}_{ij}.$$

(159)

Generally, if $\Sigma$ is embedded in the bulk with the extrinsic and intrinsic curvature parametrically independent of $\beta$, then the junction conditions for cases (I) and (II) reduce to $n^{\alpha }\left[{\nabla }_{\alpha }R\right]=0$, and

$$\left[R\right]=\left[K_{ij}\right]=0,\mathrm{or}\left[R\right]=R_0\ne 0,K_{ij}=0.$$

(160)

As an application, let us apply (160) to a well-known family of line elements which model static stars with spherical symmetry. In ordinary GR, this class of solutions is constructed by matching a Schwarzschild exterior to a perfect fluid interior with metric of the form

$$d{s}_{int}^2=-A\left(r\right)d{t}^2+B^{-1}\left(r\right)d{r}^2+r^2(d{\theta }^2+{sin}^2\theta d{\varphi }^2),$$

(161)

where

$$A\left(r\right)=e^{2\alpha \left(r\right)},B\left(r\right)=1-\frac{2m\left(r\right)}{r},{\alpha }^{\prime }\left(r\right)=\frac{m\left(r\right)+4\pi r^{3}p\left(r\right)}{r(r-2m(r\left)\right)},$$

(162)

which solves the field equations with the energy-momentum tensor of a perfect fluid

$$T_{\mu \nu }=(\rho \left(r\right)+p\left(r\right))U_{\mu }{U}_{\nu }+p\left(r\right)g_{\mu \nu },U_{\mu }=(\sqrt{A},0,0,0),$$

(163)

with the Tolman–Oppenheimer–Volkoff equation:

$$p^{\prime }\left(r\right)=-\frac{m\left(r\right)}{r^2}\rho \left(r\right)\left(1+\frac{p\left(r\right)}{\rho \left(r\right)}\right)\left(1+\frac{4\pi r^{3}p\left(r\right)}{m}\right){\left(1-\frac{2m\left(r\right)}{r}\right)}^{-1}.$$

(164)

The Schwarzschild solution is still a vacuum solution in $\mathcal{F}\left(R\right)=R+\beta R^2$ theory, and we also take the same interior solution as defined above but now this is sourced by a different energy-momentum which is the sum of (163) and extra terms arising from the field equations of the $\mathcal{F}\left(R\right)$ theory.

The surface $\Sigma$ is the star’s boundary defined by $r=R_s$ for some constant $R_s$ (which typically satisfies an appropriate Buchdahl bound). At $\Sigma$, imposing metric continuity and $\left[K_{ij}\right]=0$ lead to the boundary conditions

$$m\left(R_s\right)=M,p\left(R_s\right)=0,$$

where M is the mass parameter of the Schwarzschild exterior. Note that for the metric of the form (161), the components of the extrinsic curvature are $K_t^t=\frac{1}{2}\frac{\sqrt{B}{A}^{\prime }}{A}, K_{\theta }^{\theta }=K_{\varphi }^{\varphi }=\frac{\sqrt{B}}{R_s}$. Thus, matching it to a Schwarzschild exterior leads to $p\left(R_s\right)=0$ after using (162). Since $R=-8\pi T_{\mu }^{\mu }=8\pi (\rho -3p)$, the additional junction conditions $\left[R\right]=0, \left[{\partial }_{r}R\right]=0$ further impose the additional boundary conditions

$$\rho \left(R_s\right)=0,p^{\prime }\left(R_s\right)=0,{\rho }^{\prime }\left(R_s\right)=0.$$

(165)

For a polytropic equation of state of the form $p\propto {\rho }^{\gamma }$ for some positive constant $\gamma$, these boundary conditions are only compatible with the case of radiative matter $p=\frac{1}{3}\rho$. These results were similarly presented in [46], albeit through a different set of junction conditions. For the Tolman–Oppenheimer–Volkoff stellar model above, our junction conditions lead to an identical final set of boundary conditions on the interior fluid’s density and pressure.

Another well-studied stellar model that is also a cut-and-paste solution involving a Schwarzschild exterior is the Oppenheimer–Snyder solution where the interior is a closed FRW universe sourced by a pressureless dust. For this model of stellar formation, the matching surface $\Sigma$ is taken to preserve the $SO\left(3\right)$ isometry, and is comoving with the FRW interior of which the metric reads

$$d{s}^2=a^2\left(\tau \right)\left(-d{\tau }^2+d{R}^2+{sin}^{2}R(d{\theta }^2+{sin}^2\theta d{\varphi }^2)\right).$$

(166)

The surface $\Sigma$ is defined as the sphere $R=R_c$ for some constant $R_c$, or in the coordinates of the Schwarzschild exterior, $r\left(\tau \right)=a\left(\tau \right)sin\left(R_c\right)$, with the scale factor $a\left(\tau \right)$ satisfying the Friedmann equations for a pressureless dust. It is straightforward to see that this solution is incompatible with the junction conditions (160) since the Ricci scalar $R=8\pi \rho$.

These simple examples appear to indicate that an embedding of GR solutions with discontinuous but non-singular sources is subject to rather stringent constraints associated with the generalized junction conditions in $F\left(R\right)$ theory. It is, however, important to note that we have examined only the simplest embedding, retaining the full GR metric and the geometry of $\Sigma$. The caveat is that this implies that the interior is sourced by an energy-momentum tensor of the form

$$T_{\mu \nu }=T_{\mu \nu }^{\left(GR\right)}+\frac{\beta }{4\pi }\left((g_{\mu \nu }-{\nabla }_{\nu }{\nabla }_{\mu })R+R{R}_{\mu \nu }-\frac{1}{4}{R}^2{g}_{\mu \nu }\right).$$

It would be interesting to consider stellar models where $T_{\mu \nu }$ inherits a more physically motivated form, as well as other shapes of $\Sigma$ which may be dependent on various theory couplings in the $F\left(R\right)$ theory. This would probe a much wider landscape of solutions compatible with our generalized junction conditions.

6. Concluding Remarks

We have presented a general method to derive the appropriate Darmois–Israel junction conditions for gravitational theories with higher-order derivative terms by integrating the bulk equations of motion across the infinitesimal width of the singular hypersurface $\Sigma$, as defined in (9). A salient feature of our work is the presence of regularity constraints which impose conditions on the extrinsic curvature such that the integral in (9) converges. Geometrically, they specify the conditions under which the embedding of $\Sigma$ into the bulk spacetime is compatible with the delta-singular source localized within $\Sigma$.

Our method fundamentally relies on defining the $\delta$-distribution as the limit of a sequence of classical functions as expressed in (24) and (25). We found that the use of delta-convergent sequences yields a powerful language for organizing various terms with different orders of singularities appearing in the integral (9), and is intimately related to the procedure of Hadamard regularization commonly invoked in the theory of distributions. Upon imposing the regularity constraints, the integral in (9) converges and is well-defined. Our method passes a stringent consistency test (that is noticeably absent in previous literature): that the junction conditions for Gauss–Bonnet gravity can be obtained as a suitable limit of those of quadratic gravity when the coupling constants reduce to those of the 4D Euler density term. This is a rigorous check of validity, since the junction conditions for Gauss–Bonnet gravity can also be independently derived by boundary variation of a suitable surface term in the action.

As explicit examples of our approach, we demonstrated in detail how to obtain the regularity constraints and junction conditions for (i) quadratic gravity; (ii) $\mathcal{F}\left(R\right)$ theories; (iii) a 4D low-energy effective action in string theory; and (iv) Euler density action terms which are higher-dimensional analogues of the 4D Gauss–Bonnet term. We have expressed these generalized junction conditions explicitly as functions of the extrinsic curvature tensor and its derivatives. Generically, they also involve components that are non-parallel to $\Sigma$, in contrast to the case in ordinary GR. To our knowledge, all of these generalized junction conditions are novel results. Although there have been past attempts to derive junction conditions for quadratic gravity [17] and $\mathcal{F}\left(R\right)$ theories [14,15,18], their results or underlying methodologies did not appear to demonstrate consistency with the Gauss–Bonnet case. In this aspect, we hope that our work has also clarified some of the ambiguities encountered in these previous studies.

For $\mathcal{F}\left(R\right)$-gravity, the energy momentum tensor satisfies the conservation law as a consequence of the generalized Bianchi identity established in [31] for setup without any distributional singularities. In our approach, the delta function sources on the hypersurface are regularized and smoothed out, so the conservation law on the hypersurface can be considered as a limit of the general conservation law obtained there and remains true. Nevertheless, this is a coincidence that is neither guaranteed nor needed. In fact, there is not such a generalized Bianchi identity for a general higher-derivative gravity theory, and so the energy momentum tensor on the hypersurface is not always conserved:

$$t_{\nu }^{\mu \nu }=G^{\mu }\ne 0.$$

Nevertheless, one may interpret the non-conservation to be due to an external geometric force $G^{\mu }$ exerted by the background geometry on the hypersurface, and the junction can still be a consistent one.

The details of our derivation procedure presented here should be pedagogically useful towards adopting our methodology to derive junction conditions for other more complicated gravitational theories, including those with matter and gauge couplings. We should also mention that although our method applies rather widely to gravitational actions built out of curvature invariants, by definition, it does not apply to topological boundary terms in the action since they do not manifestly modify the bulk equations of motion. In an upcoming work, we derive and examine the generalized junction conditions for Chern–Pontryagin density terms by boundary variation of suitable surface terms. This class of theories includes, in particular, the (non-dynamical) “Chern–Simons gravity” theory in 4D of which surface term was derived in [37].

We hope that these junction equations will furnish the essential first steps towards exploring a potentially rich and phenomenologically interesting landscape of classical solutions which have singular hypersurfaces as their defining geometric feature. In this work, we have briefly touched upon a couple of applications in the $\mathcal{F}\left(R\right)=R+\beta R^2$ theory, where we found a thin-shell wormhole constructed by identifying the exterior regions of two identical copies of Schwarzschild–de Sitter spacetime. We showed that many stellar models in ordinary GR which are lifted directly to this particular $\mathcal{F}\left(R\right)$ theory violate its regularity constraints. It would be very interesting to carry out a more extensive exploration of thin-shell and stellar geometries in many gravitational theories beyond GR, now that we are freshly equipped with the fundamental junction conditions to work with. In particular, we note that thin-shell wormholes have been recently revisited as black hole mimickers for LIGO events [52].

In this paper, we have considered the general quadratic gravity theory and found that for a special linear combination, i.e., the Gauss–Bonnet gravity theory, the theory is regular and there is no need to impose any additional regularity constraint. We have also considered the simplest cubic gravity with a $R^3$ modification term and found that additional regularity constraints are needed. Similar to the quadratic case, one can expect that for a cubic theory with a certain special combinations of the possible cubic gravity terms (the independent ones are $R^3$, $R{R}_{\mu \nu }^2$, $R{R}_{\mu \nu \alpha \beta }^2$, $R_{\mu \nu }^3$, $R_{\mu \alpha }{R}_{\nu \beta }{R}^{\mu \nu \alpha \beta }$, $R^{\mu \nu }{R}_{\mu \alpha \beta \gamma }{R}_{\nu }{}^{\alpha \beta \gamma }$, $R_{\mu \nu }{}^{\alpha \beta }{R}_{\alpha \beta }{}^{\gamma \delta }{R}_{\gamma \delta }{}^{\mu \nu }$, $R_{\mu }{}^{\alpha \beta }{}_{\nu }{R}^{\gamma }{}_{\alpha \beta }{}^{\delta }{R}^{\mu }{}_{\gamma \delta }{}^{\nu }$), there is no higher order singularity and hence no need to impose any additional regularity constraint. These kinds of cubic gravity theories are somewhat special and deserve further studies.

Another natural avenue for future work lies in extending our approach to cover $\Sigma$, which is lightlike. In this case, there is no unique definition of the extrinsic curvature once the induced metric on the surface becomes degenerate, since the normal vector defined in the setting of timelike/spacelike $\Sigma$ is then tangent to $\Sigma$, and naively we need another notion of a “transverse” vector. This subtle point has been addressed in [53], where a proposal for junction conditions in the case of a null surface was presented. It would be interesting to generalize the results of [53] to higher-derivative gravitational theories, and see if some aspects of our method remain useful.