# Shock-Wave Refinement of the Friedmann–Robertson–Walker Metric

Christopher E. Alexander
Department of Mathematics
University College London
London, WC1H 0AY
United Kingdom
<EMAIL_ADDRESS>
& Blake Temple
Department of Mathematics
University of California
Davis, CA 95616
United States
<EMAIL_ADDRESS>
&Joel A. Smoller
Department of Mathematics
University of Michigan
Ann Arbor, MI 48109
United States
Joel Smoller passed away in 2017 and contributed to the original version of
this article, which was supported in part by NSF Applied Mathematics Grant
Number DMS-010-3998.

###### Abstract

The mathematics of general relativistic shock waves is introduced and
considered in a cosmological context. In particular, an expanding
Friedmann–Roberson–Walker metric is matched to a Tolman–Oppenheimer–Volkoff
metric across a spherical shock surface. This is the general relativistic
analogue of a shock-wave explosion within a static singular isothermal fluid
sphere and may be regarded as a model for the Big Bang. These shock waves are
constructed both within and beyond the Hubble radius, which corresponds to a
universe outside and inside its Schwarzschild radius respectively. Certain
self-similar perturbations of the FRW metric lead to an accelerated expansion,
even without a cosmological constant, and thus it is conjectured that such a
mechanism may account for the anomalous acceleration observed today without
recourse to dark energy.

_Keywords_ General Relativity $\cdot$ Einstein-Euler $\cdot$ Self-Similar
$\cdot$ Shock Wave $\cdot$ Cosmology $\cdot$ Dark Energy

###### Contents

  1. 1 Introduction
  2. 2 The FRW Metric
  3. 3 The General Theory of Shock Matching
  4. 4 Shock-Wave Solutions Inside the Hubble Radius - The Case $r_{*}=0$
  5. 5 Shock-Wave Solutions Beyond the Hubble Radius - The Case $r_{*}>0$
  6. 6 Self-Similar Extensions of FRW-TOV Shock-Waves Inside the Hubble Radius
  7. 7 Conclusion

##  1 Introduction

In the Standard Model of Cosmology, the expanding universe of galaxies is
described by a Friedmann–Robertson–Walker (FRW) metric, which in spherical
coordinates has a line element given by,

$ds^{2}=-dt^{2}+R^{2}(t)\left(\frac{dr^{2}}{1-kr^{2}}+r^{2}[d\theta^{2}+\sin^{2}(\theta)d\phi^{2}]\right).$
(1)

In this model, which accounts for physics on the largest length scale, the
Universe is approximated by a space of uniform density and pressure at each
fixed time, and the expansion rate is determined by the cosmological scale
factor $R(t)$ that evolves according to the Einstein field equations.
Astronomical observations indicate that the galaxies are uniform on a scale of
about one billion light-years, and the expansion is critical, that is, $k=0$
in (1). According to (1) with $k=0$, on the largest scale the Universe is
infinitely flat Euclidean space $\mathbb{R}^{3}$ at each fixed time. Matching
the Hubble constant to its observed values and invoking the Einstein field
equations, the FRW model implies that the entire infinite universe
$\mathbb{R}^{3}$ emerged all at once from a singularity ($R=0$) some 13.7
billion years ago, and this event is referred to as the _Big Bang_.

This article, updated in 2023, summarises the work of Smoller and Temple in
[11, 12, 18, 19], and then describes subsequent advances by Smoller, Temple
and Alexander in [22, 20, 1]. The Smoller–Temple solutions describe a two-
parameter family of exact solutions to the perfect fluid Einstein field
equations that refine the FRW metric by a spherical shock-wave cutoff. The
Einstein field equations with a perfect fluid source are also referred to as
the _Einstein-Euler_ equations and we make the distinction that an _exact_
solution is one that is defined by a system of ODE, whereas an _explicit_
solution is one that can be given in terms of elementary mathematical
functions. In the original Smoller–Temple solutions, the flat ($k=0$) FRW
spacetime approximates the expanding wave created behind the general
relativistic version of an explosion into a static singular isothermal sphere.
There are two cases: matching _outside the black hole_ , which means the shock
surface is inside the Hubble radius of the FRW spacetime, and _inside the
black hole_ , that is, beyond the Hubble radius. In the examples given in the
original article, Smoller and Temple recognised that the flat FRW spacetime
did not fully resolve the expanding wave behind the shock, but argued that it
was qualitatively close. After Smoller and Temple introduced a family of
expanding waves in [20], Alexander in [1] accomplished the goal of fully
resolving the expanding wave behind these general relativistic shock waves for
the case of shock waves inside the Hubble radius, that is, outside the black
hole. This is accomplished by incorporating self-similar perturbations of the
FRW spacetime, characterised in [20] as the natural self-similar extension of
the FRW spacetime, that is, spacetimes that are asymptotically FRW at the
origin. Alexander’s solution for pure radiation resolves a long-standing open
problem first proposed by Taub in [3], and later taken up from a different
point of view by Smoller and Temple. Most interestingly, the expanding waves
behind these shocks naturally introduce accelerations relative to the flat FRW
spacetime. Since the flat FRW spacetime is the starting model of Cosmology, it
is intriguing to consider these models as an explanation for cosmic
acceleration without recourse to a cosmological constant. We argue here that
the analysis needs to be extended to a parallel family of solutions beyond the
Hubble radius, as accomplished in [18, 19] for flat FRW shock matching, to
make a plausible case for such an explanation. We first review these
Smoller–Temple solutions and then return in Section 6 to discuss the
refinements made by Smoller, Temple and Alexander.

In order to construct a mathematically simple family of shock-wave refinements
of the FRW metric that satisfy the Einstein-Euler equations exactly, we assume
critical expansion ($k=0$) and restrict to the case that the sound speed in
the fluid on the FRW side of the shock wave is constant, that is, we assume an
FRW equation of state $p=\sigma\rho$ with $0<\sigma\leq c^{2}$, where
$\sigma$, the square of the sound speed $\sqrt{\frac{\partial
p}{\partial\rho}}$, is constant. For $\sigma=c^{2}/3$, this equation of state
describes a state of matter known as _pure radiation_ , as well as the
equation of state of the relativistic limit of free particles, which is
correct during the Radiation Dominated Epoch of the Early Universe [24]. Also,
as $\sigma$ ranges from $0$ to $c^{2}$, we obtain qualitatively correct
approximations to general equations of state. Now by using units such that the
speed of light $c$ and gravitational constant $\mathcal{G}$ are set to unity,
the family of solutions is then determined by two parameters: $0<\sigma\leq 1$
and $r_{*}\geq 0$. The second parameter $r_{*}$ is the FRW radial coordinate
$r$ of the shock in the limit $t\rightarrow 0$, that is, the instant of the
Big Bang.111Since when $k=0$ the FRW metric is invariant under the rescaling
$r\rightarrow\alpha r$ and $R\rightarrow R/\alpha$, we fix the radial
coordinate $r$ by fixing the scale factor $\alpha$ with the condition that
$R(t_{0})=1$ for some time $t_{0}$, say present time. The FRW radial
coordinate $r$ is singular with respect to radial arc-length $\bar{r}=rR$ at
the Big Bang ($R=0$), so setting $r_{*}>0$ does not place the shock wave away
from the origin at time $t=0$. The distance from the FRW centre to the shock
wave tends to zero in the limit $t\rightarrow 0$, even when $r_{*}>0$. In the
limit $r_{*}\rightarrow\infty$, we recover from the family of solutions the
usual (infinite) FRW metric with equation of state $p=\sigma\rho$, that is, we
recover the standard FRW metric in the limit that the shock wave is infinitely
far out. In this sense, our family of solutions of the Einstein-Euler
equations represents a two-parameter refinement of the standard
Friedmann–Robertson–Walker metric.

The explicitly defined solutions in the case $r_{*}=0$ were first constructed
in [12] and are qualitatively different from the exact solutions in the case
$r_{*}>0$, which were constructed later in [19]. The difference is that when
$r_{*}=0$, the shock wave lies closer than one Hubble length from the centre
of the FRW spacetime throughout its motion, but when $r_{*}>0$, the shock wave
emerges at the Big Bang at a distance beyond one Hubble length [17]. The
Hubble length, also referred to as the Hubble radius in certain contexts,
depends on time, and tends to zero as $t\rightarrow 0$. We show in [19] that
one Hubble length, equal to $c/H$ where $H=\dot{R}/R$, is a critical length
scale in a flat FRW metric because the total mass inside one Hubble length has
a Schwarzschild radius equal to exactly one Hubble length.222Since $c/H$ is a
good estimate for the age of the Universe, it follows that the Hubble length
$c/H$ is approximately the distance of light travel starting at the Big Bang
up until present time. In this sense, the Hubble length is a rough estimate
for the distance to the further most objects visible in the Universe. That is,
one Hubble length marks precisely the distance at which the Schwarzschild
radius $\bar{r}_{s}\equiv 2M$ (of the mass $M$ inside a radial shock wave at
distance $\bar{r}$ from the FRW centre) crosses from inside
($\bar{r}_{s}<\bar{r}$) to outside ($\bar{r}_{s}>\bar{r}$) the shock wave. If
the shock wave is at a distance closer than one Hubble length from the FRW
centre, then $\bar{r}>2M$ and we say that the solution lies _outside the black
hole_ , but if the shock wave is at a distance greater than one Hubble length,
then $\bar{r}<2M$ at the shock, and we say the solution lies _inside the black
hole_. Since $M$ increases proportional to $\bar{r}^{3}$, it follows that
$\bar{r}>2M$ for $\bar{r}$ sufficiently small, and $\bar{r}<2M$ for $\bar{r}$
sufficiently large, so there must be a critical radius at which $\bar{r}=2M$.
In Section 2 (taken from [18, 19]), we show that when $k=0$, this critical
radius is exactly the Hubble radius. When the parameter $r_{*}=0$, the family
of solutions for $0<\sigma\leq 1$ starts at the Big Bang, and evolves
thereafter outside the black hole, satisfying $\bar{r}>2M$ everywhere from
$t=0$ onward. But when $r_{*}>0$, the shock wave is further out than one
Hubble length at the instant of the Big Bang, and the solution begins with
$\bar{r}<2M$ at the shock wave. From this time onward, the spacetime expands
until eventually the Hubble radius catches up to the shock wave at
$\bar{r}=2M$ and then passes the shock wave, making $\bar{r}>2M$ thereafter.
Thus when $r_{*}>0$, the whole spacetime begins inside the black hole (with
$\bar{r}<2M$ for sufficiently large $\bar{r}$) but eventually evolves to a
solution outside the black hole. The time when $\bar{r}=2M$ actually marks the
event horizon of a _white hole_ , the time reversal of a black hole, in the
ambient spacetime beyond the shock wave. We show that when $r_{*}>0$, the time
when the Hubble length catches up to the shock wave comes before the time when
the shock wave comes into view at the FRW centre. Furthermore, when
$\bar{r}=2M$, assuming $t$ is so large that we can neglect the pressure from
this time onward, the whole solution emerges from the white hole as a finite
ball of mass expanding into empty space, satisfying $\bar{r}>2M$ everywhere
thereafter. In fact, when $r_{*}>0$, the zero pressure Oppenheimer–Snyder
solution outside the black hole gives the large time asymptotics of the
solution (see [9, 15, 19] and the comments after Theorems 6-8 below).

The explicitly defined solutions in the case $r_{*}=0$ give a general
relativistic version of an explosion into a static singular isothermal fluid
sphere, qualitatively similar to the corresponding classical explosion outside
the black hole [12]. The main difference physically between the cases
$r_{*}>0$ and $r_{*}=0$ is that when $r_{*}>0$, that is, the case when the
shock wave emerges from the Big Bang at a distance beyond one Hubble length, a
large region of uniform expansion is created behind the shock wave at the
instant of the Big Bang. Thus, when $r_{*}>0$, lightlike (also known as
_null_) information about the shock wave propagates inward from the wave,
rather than outward from the centre, as is the case when $r_{*}=0$ and the
shock lies inside one Hubble length.333One can imagine that when $r_{*}>0$,
the shock wave could get out through a great deal of matter early on when
everything is dense and compressed and still not violate the speed of light
bound. Thus when $r_{*}>0$, the shock wave would _thermalise_ , or more
accurately, _make uniform_ a large region at the centre early on in the
explosion. The authors speculate that such a mechanism might offer a
substitute for inflation in this regard. It follows that when $r_{*}>0$, an
observer positioned in the FRW spacetime inside the shock wave will see
exactly what the Standard Model of Cosmology predicts up until the time when
the shock wave comes into view in the far field. In this sense, the case
$r_{*}>0$ gives a black hole cosmology that refines the standard FRW model of
Cosmology to the case of finite mass. One of the surprising differences
between the case $r_{*}=0$ and the case $r_{*}>0$ is that, when $r_{*}>0$, the
important pure radiation equation of state $p=\frac{1}{3}\rho$ comes out of
the analysis as special at the Big Bang. When $r_{*}>0$, the shock wave
emerges at the instant of the Big Bang at a finite non-zero speed (the speed
of light) only for the special value $\sigma=\frac{1}{3}$. In this case, the
equation of state on both sides of the shock wave tends to the correct
relation $p=\frac{1}{3}\rho$ as $t\rightarrow 0$ and the shock wave
decelerates to a subluminal speed for all positive times thereafter (see [18,
19] and Theorem 8 below).

In all cases $0<\sigma\leq 1$, $r_{*}\geq 0$, the spacetime metric that lies
beyond the shock wave is taken to be a metric of Tolman–Oppenheimer–Volkoff
(TOV) form [8], that is,

$ds^{2}=-B(\bar{r})d\bar{t}^{2}+\frac{1}{A(\bar{r})}d\bar{r}^{2}+\bar{r}^{2}[d\theta^{2}+\sin^{2}(\theta)d\phi^{2}].$
(2)

The metric (2) is in _standard Schwarzschild coordinates_ , that is, diagonal
with a radial coordinate defined by spheres of symmetry. Furthermore, the
metric components depend only on the radial coordinate $\bar{r}$. Barred
coordinates are used to distinguish TOV coordinates from unbarred FRW
coordinates for the purpose of matching the metric at the shock surface (see
Section 3). The mass function $M(\bar{r})$ enters as a metric component
through the relation,

$A=1-\frac{2M(\bar{r})}{\bar{r}}.$ (3)

The TOV metric (2) has a very different character depending on whether $A>0$
or $A<0$, that is, depending on whether the solution lies outside or inside
the black hole. In the case $A>0$, $\bar{r}$ is a spacelike coordinate and the
TOV metric describes a static general relativistic fluid sphere.444The metric
(2) is, for example, the starting point for the stability limits of Buchdahl
and Chandrasekhar for stars [24, 13, 16]. When $A<0$, $\bar{r}$ is the
timelike coordinate and (2) is a dynamic metric that evolves in time. The
explicitly defined shock-wave solutions are obtained by taking $\bar{r}=R(t)r$
to match the spheres of symmetry, and then matching the metrics (1) and (2) at
an interface $\bar{r}=\bar{r}(t)$ across which the metrics are Lipschitz
continuous. This can be done in general. In order for the interface to be a
physically meaningful shock surface, we use the result in Theorem 4 below (see
[11]), so that a single additional conservation constraint is sufficient to
rule out delta function sources at the shock and guarantee that the matched
metric solves the Einstein-Euler equations in the weak sense.555The Einstein
field equations $G=\kappa T$ are second-order in the metric and so delta
function sources will in general be present at a Lipschitz continuous matching
of metrics. The Lipschitz continuous matching of the metrics, together with
the conservation constraint, leads to a system of ordinary differential
equations (ODE) that determine the shock position, together with the TOV
density and pressure at the shock. Since the TOV metric depends only on
$\bar{r}$, the equations thus determine the TOV spacetime beyond the shock
wave. To obtain a physically meaningful outgoing shock wave, we impose the
constraint $\bar{p}\leq\bar{\rho}$ to ensure that the equation of state on the
TOV side of the shock is qualitatively reasonable, and we require that the
shock be compressive as the entropy condition. For an outgoing shock wave,
this corresponds to the conditions $\rho>\bar{\rho}$ and $p>\bar{p}$, that is,
the pressure and density need to be larger on the side of the shock that
receives the mass flux (the FRW side when the shock wave is propagating away
from the FRW centre). This condition breaks the time reversal symmetry of the
equations and is sufficient to rule out rarefaction shocks in classical gas
dynamics [10, 19]. The ODE, together with the equation of state bound and the
conservation and entropy constraints, determine a unique solution of the ODE
for every $0<\sigma\leq 1$ and $\bar{r}_{*}\geq 0$, and this provides the two-
parameter family of solutions discussed here [12, 19].

The Lipschitz matching of the metrics implies that the total mass $M$ is
continuous across the interface, and so when $r_{*}>0$, the total mass of the
entire solution (inside and outside the shock wave) is finite at each time
$t>0$. Both the FRW and TOV spacetimes emerge at the Big Bang. The total mass
$M$ on the FRW side of the shock has the meaning of total mass inside radius
$\bar{r}$ at fixed time, but on the TOV side of the shock, $M$ does not evolve
according to equations that give it the interpretation as a total mass because
the metric is inside the black hole. Nevertheless, after the spacetime emerges
from the black hole, the total mass takes on its usual meaning outside the
black hole and time asymptotically the Big Bang ends with an expansion of
finite total mass in the usual sense. Thus, when $r_{*}>0$, our shock-wave
refinement of the FRW metric leads to a Big Bang of finite total mass.

The Smoller–Temple family of shock-wave solutions are rough models in the
sense that the equation of state on the FRW side has constant $\sigma$ and the
equation of state on the TOV side is determined by the equations and therefore
cannot be imposed. For more accurate equations of state, a more accurate
description of the expanding wave created behind the shock is needed to meet
the conservation constraint and thereby mediate the transition across the
shock wave. At the time of publish of the original article, the authors
thought such expanding waves to be pretty much impossible to model as exact
solutions, however, in the most recent work of Alexander [1], this has been
resolved for shock waves inside the Hubble radius by considering self-similar
perturbations of the flat FRW spacetime. Not only do these modifications
permit enough parameter freedom to impose the same equation of state on each
side of the shock, but the perturbations of the flat FRW spacetime induce an
accelerated expansion without the presence of a cosmological constant. A
rigorous proof for the existence of a general relativistic shock wave inside
the Hubble radius with a pure radiation equation of state each side of the
shock is also demonstrated in [1]. The fact that we can find global solutions
that meet our physical bounds within and beyond the Hubble radius, that are
qualitatively the same for all values of $\sigma\in(0,1]$ and all initial
shock positions, strongly suggests that such a shock wave would be the
dominant wave in a large class of problems.

In Section 2 we derive the FRW solution for constant $\sigma$ and discuss the
Hubble radius as a critical length scale. In Section 3 we state the general
theorems in [11] for matching gravitational metrics across shock waves. In
Section 4 we discuss the construction of the family of solutions in the case
$r_{*}=0$. In Section 5 we discuss the case $r_{*}>0$, and in Section 6 we
discuss the advances made since 2005, in particular, replacing the flat FRW
solution with a family of self-similar perturbations and discussing how such
waves exhibit an accelerated expansion. See [12, 19, 20, 1] for details.

##  2 The FRW Metric

According to Einstein’s Theory of General Relativity, all properties of the
gravitational field are determined by a Lorentzian spacetime metric tensor
$g$, whose line element in a given coordinate system
$\boldsymbol{x}=(x^{0},...,x^{3})$ is given by

$ds^{2}=g_{\mu\nu}dx^{\mu}dx^{\nu}.$ (4)

We use the Einstein summation convention whereby repeated up-down indices are
assumed to be summed from 0 to 3. The components $g_{\mu\nu}$ of the
gravitational metric $g$ satisfy the Einstein field equations,

$\displaystyle G^{\mu\nu}$ $\displaystyle=\kappa T^{\mu\nu},$ $\displaystyle
T^{\mu\nu}$ $\displaystyle=(\rho c^{2}+p)u^{\mu}u^{\nu}+pg^{\mu\nu},$ (5)

where we assume the stress-energy-momentum tensor $T$ is that of a perfect
fluid, making these equations the _Einstein-Euler_ equations. Here $G$ is the
Einstein curvature tensor,

$\kappa=\frac{8\pi\mathcal{G}}{c^{4}}$ (6)

is the coupling constant, $\mathcal{G}$ is Newton’s gravitational constant,
$c$ is the speed of light, $\rho c^{2}$ is the energy density, $p$ is the
pressure and $\boldsymbol{u}=(u^{0},...,u^{3})$ is the fluid four-velocity
[24]. We will also use the convention $c=1$ and $\mathcal{G}=1$ when
convenient.

Putting the metric ansatz (1) into the Einstein-Euler equations (5) gives the
equations for the components of FRW metric [24],

$\displaystyle H^{2}=\left(\frac{\dot{R}}{R}\right)^{2}$
$\displaystyle=\frac{\kappa}{3}\rho-\frac{k}{R^{2}},$ (7)
$\displaystyle\dot{\rho}$ $\displaystyle=-3(p+\rho)H.$ (8)

The unknowns $R$, $\rho$ and $p$ are assumed to be functions of the FRW
coordinate time $t$ alone, with the dot denoting differentiation with respect
to $t$. To verify that the Hubble radius $\bar{r}_{crit}=1/H$ is the limit for
FRW-TOV shock matching outside a black hole, write the FRW metric (1) in
standard Schwarzschild coordinates $\bar{\boldsymbol{x}}=(\bar{t},\bar{r})$
where the metric takes the form

$ds^{2}=-B(\bar{t},\bar{r})d\bar{t}^{2}+\frac{1}{A(\bar{t},\bar{r})}d\bar{r}^{2}+\bar{r}^{2}d\Omega^{2},$
(9)

and the mass function $M(\bar{t},\bar{r})$ is defined through the relation

$A=1-\frac{2M}{\bar{r}}.$ (10)

It is well known that a general spherically symmetric metric can be put in the
form (9) by coordinate transformation [24, 5]. Substituting $\bar{r}=Rr$ into
(1) and diagonalising the resulting metric (see [19] for details), we obtain

$ds^{2}=\frac{1}{1-kr^{2}-H^{2}\bar{r}^{2}}\left(-\frac{1-kr^{2}}{\psi^{2}}d\bar{t}^{2}+d\bar{r}^{2}\right)+\bar{r}^{2}d\Omega^{2},$
(11)

where $\psi$ is an integrating factor that solves the equation

$\frac{\partial}{\partial\bar{r}}\left(\psi\frac{1-kr^{2}-H^{2}\bar{r}^{2}}{1-kr^{2}}\right)-\frac{\partial}{\partial
t}\left(\psi\frac{H\bar{r}}{1-kr^{2}}\right)=0,$ (12)

and the time coordinate $\bar{t}=\bar{t}(t,\bar{r})$ is defined by the exact
differential

$d\bar{t}=\left(\psi\frac{1-kr^{2}-H^{2}\bar{r}^{2}}{1-kr^{2}}\right)dt+\left(\psi\frac{H\bar{r}}{1-kr^{2}}\right)d\bar{r}.$
(13)

Now using (10) in (7), it follows that

$M(t,\bar{r})=\frac{\kappa}{2}\int_{0}^{\bar{r}}\rho(t)s^{2}ds=\frac{1}{3}\frac{\kappa}{2}\rho\bar{r}^{3}.$
(14)

Since in the FRW metric $\bar{r}=Rr$ measures arc-length along radial
geodesics at fixed time, we see from (14) that $M(t,\bar{r})$ has the physical
interpretation as the total mass inside radius $\bar{r}$ at time $t$ in the
FRW metric. Restricting to the case of critical expansion ($k=0$), we see from
(7), (13) and (14) that $\bar{r}=1/H$ is equivalent to $\bar{r}=2M$, and so at
fixed time $t$, the following equivalences are valid:

$\bar{r}=\frac{1}{H}\iff\frac{2M}{\bar{r}}=1\iff A=0.$ (15)

We conclude that $\bar{r}=1/H$ is the critical length scale for the FRW metric
at fixed time $t$ in the sense that $A$ changes sign at $\bar{r}=1/H$, and so
the universe lies inside a black hole beyond $\bar{r}=1/H$, as claimed above.
It is shown in [16] that the standard TOV metric outside the black hole cannot
be continued into $A=0$ except in the very special case $\rho=0$, as it takes
an infinite pressure to hold up a static configuration at the event horizon of
a black hole. Thus to do shock matching beyond one Hubble length requires a
metric of a different character, and for this purpose, in [18, 19] we
introduce the TOV metric _inside the black hole_ , that is, a metric of TOV
form with $A<0$ and whose fluid is co-moving with the timelike radial
coordinate $\bar{r}$.

The Hubble radius $\bar{r}_{crit}=c/H$ is also the critical distance at which
the outward expansion of the FRW metric exactly cancels the inward advance of
a radial light ray impinging on an observer positioned at the origin of a flat
($k=0$) FRW metric. Indeed, by (1), a light ray travelling radially inward
toward the centre of an FRW coordinate system satisfies,

$c^{2}dt^{2}=R^{2}dr^{2},$ (16)

so that

$\frac{d\bar{r}}{dt}=\dot{R}r+R\dot{r}=H\bar{r}-c=H\left(\bar{r}-\frac{c}{H}\right)>0,$
(17)

if and only if $\bar{r}>c/H$. Thus the arc-length distance from the origin to
an inward moving light ray at fixed time $t$ in a flat FRW metric will
actually increase as long as the light ray lies beyond the Hubble radius. An
inward moving light ray will, however, eventually cross the Hubble radius and
reach the origin in finite proper time due to the increase in the Hubble
length with time.

We now calculate the infinite redshift limit in terms of the Hubble length. It
is well known that light emitted at $(t_{e},r_{e})$ at wavelength
$\lambda_{e}$ in an FRW spacetime will be observed at $(t_{0},r_{0})$ at
wavelength $\lambda_{0}$ if

$\frac{R_{0}}{R_{e}}=\frac{\lambda_{0}}{\lambda_{e}}.$ (18)

Moreover, the redshift factor $z$ is defined by

$z=\frac{\lambda_{0}}{\lambda_{e}}-1.$ (19)

Thus, infinite redshift occurs in the limit $R_{e}\rightarrow 0$, where $R=0$,
$t=0$ is the Big Bang. Consider now a light ray emitted at the instant of the
Big Bang and observed at the FRW origin at present time $t=t_{0}$. Let
$r_{\infty}$ denote the FRW coordinate at time $t\rightarrow 0$ of the
furthermost objects that can be observed at the FRW origin before time
$t=t_{0}$. Then $r_{\infty}$ marks the position of objects at time $t=0$ whose
radiation would be observed as infinitely redshifted (assuming no scattering).
Note then that a shock wave emanating from $\bar{r}=0$ at the instant of the
Big Bang will be observed at the FRW origin before present time $t=t_{0}$ only
if its position $r$ at the instant of the Big Bang satisfies $r<r_{\infty}$.
To estimate $r_{\infty}$, note first that from (16) it follows that an
incoming radial light ray in an FRW metric follows a null trajectory $r=r(t)$
if

$r-r_{e}=-\int_{t_{e}}^{t}\frac{d\tau}{R(\tau)},$ (20)

and thus

$r_{\infty}=\int_{0}^{t_{0}}\frac{d\tau}{R(\tau)}.$ (21)

Using this, the following theorem is proved in [19].

######  Theorem 1.

If the pressure $p$ satisfies the bounds

$0\leq p\leq\frac{1}{3}\rho,$ (22)

then for any equation of state, the age of the Universe $t_{0}$ and the
infinite red shift limit $r_{\infty}$ are bounded in terms of the Hubble
length by:

$\displaystyle\frac{1}{2H_{0}}\leq$ $\displaystyle t_{0}\leq\frac{2}{3H_{0}},$
(23) $\displaystyle\frac{1}{H_{0}}\leq$ $\displaystyle
r_{\infty}\leq\frac{2}{H_{0}}.$ (24)

Where we have assumed that $R=0$ when $t=0$ and $R=1$ when $t=t_{0}$,
$H=H_{0}$.

The next theorem gives closed form solutions of the FRW equations (7), (8) for
constant $\sigma$. As a special case, we recover the bounds in (23) and (24)
from the cases $\sigma=0$ and $\sigma=\frac{1}{3}$.

######  Theorem 2.

Assume $k=0$ and the equation of state $p=\sigma\rho$, where $0\leq\sigma\leq
1$ is constant. Then, for an expanding FRW universe ($\dot{R}>0$), the
solution of system (7), (8) satisfying $R=0$ at $t=0$ and $R=1$ at $t=t_{0}$
is given by,

$\displaystyle\rho$
$\displaystyle=\frac{4}{3\kappa(1+\sigma)^{2}}\frac{1}{t^{2}},$ (25)
$\displaystyle R$
$\displaystyle=\left(\frac{t}{t_{0}}\right)^{\frac{2}{3(1+\sigma)}},$ (26)
$\displaystyle\frac{H}{H_{0}}$ $\displaystyle=\frac{t_{0}}{t}.$ (27)

Moreover, the age of the Universe $t_{0}$ and the infinite redshift limit
$r_{\infty}$ are given explicitly in terms of the Hubble length by:

$\displaystyle t_{0}$ $\displaystyle=\frac{2}{3(1+\sigma)}\frac{1}{H_{0}},$
(28) $\displaystyle r_{\infty}$
$\displaystyle=\frac{2}{1+3\sigma}\frac{1}{H_{0}}.$ (29)

From (29) we conclude that a shock wave will be observed at the FRW origin
before present time $t=t_{0}$ only if its position $r$ at the instant of the
Big Bang satisfies $r<r_{\infty}$. Note that $r_{\infty}$ ranges from one half
to two Hubble lengths as $\sigma$ ranges from $1$ to $0$, taking the
intermediate value of one Hubble length at $\sigma=\frac{1}{3}$. Note also
that by using (25)-(26) in (14), it follows that

$M=\frac{\kappa}{2}\int_{0}^{\bar{r}}\rho(t)s^{2}ds=\frac{2\bar{r}^{3}}{9(1+\sigma)^{2}t_{0}^{\frac{2}{1+\sigma}}}\frac{1}{t^{\frac{2\sigma}{1+\sigma}}},$
(30)

so $\dot{M}<0$ if $\sigma>0$. It follows that if $p=\sigma\rho$, and $\sigma$
is a positive constant, then the total mass inside a radius of constant $r$
decreases in time.

##  3 The General Theory of Shock Matching

The matching of the FRW and TOV metrics in the next two sections is based on
the following theorems that are derived in [11].666Theorems 3 and 4 apply to
non-null shock surfaces.

######  Theorem 3.

Let $\Sigma$ denote a smooth three-dimensional shock surface with spacelike
normal vector $\boldsymbol{n}$ relative to the spacetime metric $g$, let $K$
denote the second fundamental form on $\Sigma$ and let $G$ denote the Einstein
curvature tensor. Assume that the components $g_{\mu\nu}$ of the gravitational
metric $g$ are continuous up to the boundary on either side separately and
Lipschitz continuous across $\Sigma$ in some fixed coordinate system. Then the
following statements are equivalent:

  1. 1.

$[K]=0$ at each point of $\Sigma$.

  2. 2.

The curvature tensors $R^{\mu}_{\nu\sigma\tau}$ and $G_{\mu\nu}$, viewed as
second-order operators on the metric components $g_{\mu\nu}$, produce no delta
function sources on $\Sigma$.

  3. 3.

For each point $P\in\Sigma$ there exists a $C^{1,1}$ coordinate transformation
defined in a neighbourhood of $P$, such that, in the new coordinates (which
can be taken to be the Gaussian normal coordinates for the surface), the
metric components are $C^{1,1}$ functions of these coordinates.

  4. 4.

For each $P\in\Sigma$, there exists a coordinate frame that is locally
Lorentzian at $P$, and can be reached within the class of $C^{1,1}$ coordinate
transformations.

Moreover, if any one of these equivalencies hold, then the Rankine–Hugoniot
jump conditions,

$[G^{\mu\nu}]n_{\mu}=0$ (31)

hold at each point on $\Sigma$.

The Rankine–Hugoniot jump conditions express the weak form of conservation of
energy and momentum across $\Sigma$ when $G=\kappa T$. Here $[f]$ denotes the
jump in the quantity $f$ across $\Sigma$, which is determined by the metric
separately on each side of $\Sigma$ since $g_{\mu\nu}$ is only Lipschitz
continuous across $\Sigma$. The notation $C^{1,1}$ denotes a function whose
first derivatives are Lipschitz continuous. In the case of spherical symmetry,
a stronger result holds. In this case, the jump conditions (31) are implied by
the single condition

$[G^{\mu\nu}]n_{\mu}n_{\nu}=0$ (32)

so long as the shock surface is not null and the areas of the spheres of
symmetry match smoothly at the shock and change monotonically as the shock
evolves. Note that in general, assuming that the angular variables are
identified across the shock, we expect conservation to entail two conditions,
one for the time and one for the radial components. The fact that the smooth
matching of the spheres of symmetry reduces conservation to one condition can
be interpreted as an instance of the general principle that directions of
smoothness in the metric imply directions of conservation of the sources.

######  Theorem 4.

Assume that $g$ and $\bar{g}$ are two spherically symmetric metrics that match
Lipschitz continuously across a three dimensional shock surface $\Sigma$ to
form the matched metric $g\cup\bar{g}$. That is, assume that $g$ and $\bar{g}$
are Lorentzian metrics given respectively by

$\displaystyle ds^{2}$
$\displaystyle=-a(t,r)dt^{2}+b(t,r)dr^{2}+c(t,r)d\Omega^{2},$ (33)
$\displaystyle d\bar{s}^{2}$
$\displaystyle=-\bar{a}(\bar{t},\bar{r})d\bar{t}^{2}+\bar{b}(\bar{t},\bar{r})d\bar{r}^{2}+\bar{c}(\bar{t},\bar{r})d\Omega^{2},$
(34)

and that there exists a smooth coordinate transformation
$\Psi:(t,r)\rightarrow(\bar{t},\bar{r})$, defined in a neighbourhood of a
shock surface $\Sigma$ given by $r=r(t)$, such that the metrics agree on
$\Sigma$ (with the implicit assumption that $\theta$ and $\varphi$ are
identified). Moreover, assume that

$c(t,r)=\bar{c}(\Psi(t,r)),$ (35)

in an open neighbourhood of the shock surface $\Sigma$, so that, in
particular, the areas of the two-spheres of symmetry in the barred and
unbarred metrics agree at the shock surface. Furthermore, assume that the
shock surface $r=r(t)$ in unbarred coordinates is mapped to the surface
$\bar{r}=\bar{r}(\bar{t})$ by $(\bar{t},\bar{r}(\bar{t}))=\Psi(t,r(t))$, that
the normal $\boldsymbol{n}$ to $\Sigma$ is non-null, and that

$\boldsymbol{n}(c)\neq 0$ (36)

where $\boldsymbol{n}(c)$ denotes the derivative of $c$ in the direction of
the vector $\boldsymbol{n}$.777That is, we assume that the areas of the two-
spheres of symmetry change monotonically in the direction normal to the
surface. For example, if $c=r^{2}$, then $\frac{\partial c}{\partial t}=0$, so
the assumption $\boldsymbol{n}(c)\neq 0$ is valid except when
$\boldsymbol{n}=\frac{\partial}{\partial t}$, in which case the rays of the
shock surface would be spacelike. Thus the shock speed would be faster than
the speed of light if our assumption $\boldsymbol{n}(c)\neq 0$ failed in the
case $c=r^{2}$. Then the following are equivalent:

  1. 1.

The components of the metric $g\cup\bar{g}$ in any Gaussian normal coordinate
system are $C^{1,1}$ functions of these coordinates across the surface
$\Sigma$.

  2. 2.

$[G^{\mu\nu}]n_{\mu}=0$.

  3. 3.

$[G^{\mu\nu}]n_{\mu}n_{\nu}=0$.

  4. 4.

$[K]=0$.

where $[f]=\bar{f}-f$ denotes the jump in the quantity $f$ across $\Sigma$,
and $K$ is the second fundamental form on the shock surface.

##  4 Shock-Wave Solutions Inside the Hubble Radius - The Case $r_{*}=0$

To construct the family of shock-wave solutions for parameter values
$0<\sigma\leq 1$ and $r_{*}=0$, we match the explicitly defined solution
(25)-(27) of the FRW metric (1) to the explicitly given TOV metric (2) outside
the black hole, that is, with $A>0$. In this case, we can bypass the problem
of deriving and solving the ODE for the shock surface and constraints
discussed above, by actually deriving the explicit solution of the Einstein-
Euler equations of TOV form that meets these equations. This explicitly
defined solution represents the general relativistic version of a static
singular isothermal fluid sphere, that is, singular because it has an inverse
square density profile and isothermal because the relationship between the
density and pressure is $\bar{p}=\bar{\sigma}\bar{\rho}$ with constant
$\sigma$.

Assuming the stress-energy-momentum tensor for a perfect fluid and assuming
that the density and pressure depend only on $\bar{r}$, the Einstein-Euler
equations for the TOV metric (2) outside the black hole are equivalent to the
Oppenheimer–Volkoff system:

$\displaystyle\frac{dM}{d\bar{r}}$ $\displaystyle=4\pi\bar{r}^{2}\bar{\rho},$
(37) $\displaystyle-\bar{r}^{2}\frac{d\bar{p}}{d\bar{r}}$
$\displaystyle=\mathcal{G}M\bar{\rho}\left(1+\frac{\bar{p}}{\bar{\rho}}\right)\left(1+\frac{4\pi\bar{r}^{3}\bar{p}}{M}\right)\left(1-\frac{2\mathcal{G}M}{\bar{r}}\right)^{-1}.$
(38)

Integrating (37) we obtain the usual interpretation of $M$ as the total mass
inside radius $\bar{r}$,

$M(\bar{r})=\int_{0}^{\bar{r}}4\pi\xi^{2}\bar{\rho}(\xi)d\xi.$ (39)

The metric component $B$ is determined from $\bar{\rho}$ and $M$ through the
equation

$\frac{1}{B}\frac{dB}{d\bar{r}}=-\frac{2}{\bar{p}+\bar{\rho}}\frac{d\bar{p}}{d\bar{r}}.$
(40)

Assuming

$\displaystyle\bar{p}$ $\displaystyle=\bar{\sigma}\bar{\rho},$ (41)
$\displaystyle\bar{\rho}$ $\displaystyle=\frac{\gamma}{\bar{r}^{2}},$ (42)

for some constants $\bar{\sigma}$ and $\gamma$, then substituting this into
(39), we obtain

$M(\bar{r})=4\pi\gamma\bar{r}.$ (43)

Putting (41)-(43) into (38) and simplifying yields the identity

$\gamma=\frac{1}{2\pi\mathcal{G}}\left(\frac{\bar{\sigma}}{1+6\bar{\sigma}+\bar{\sigma}^{2}}\right).$
(44)

From (39) we obtain

$A=1-8\pi\mathcal{G}\gamma<1.$ (45)

Applying (40) leads to

$B=B_{0}\left(\frac{\bar{\rho}}{\bar{\rho}_{0}}\right)^{-\frac{2\bar{\sigma}}{1+\bar{\sigma}}}=B_{0}\left(\frac{\bar{r}}{\bar{r}_{0}}\right)^{\frac{4\bar{\sigma}}{1+\bar{\sigma}}}.$
(46)

By rescaling the time coordinate, we can take $B_{0}=1$ at $\bar{r}_{0}=1$, in
which case (46) reduces to

$B=\bar{r}^{\frac{4\bar{\sigma}}{1+\bar{\sigma}}}.$ (47)

We conclude that when (44) holds, (41)-(45) and (46) provide an explicit
solution of the Einstein-Euler equations of TOV type for each
$0<\bar{\sigma}\leq 1$.888In this case, an explicit solution of TOV type was
first found by Tolman [23] and rediscovered in the case
$\bar{\sigma}=\frac{1}{3}$ by Misner and Zapolsky ([24], page 320). By (45),
these solutions are defined outside the black hole since $\bar{r}>2M$. When
$\bar{\sigma}=\frac{1}{3}$, (44) yields $\gamma=\frac{3}{56\pi{\cal G}}$
([24], equation (11.4.13)).

To match the explicitly given FRW solution (25)-(27) with equation of state
$p=\sigma\rho$ to the explicitly given TOV solution (41)-(47) with equation of
state $\bar{p}=\bar{\sigma}\bar{\rho}$ across a shock interface, we first set
$\bar{r}=Rr$ to match the spheres of symmetry and then match the timelike and
spacelike components of the corresponding metrics in standard Schwarzschild
coordinates. The matching of the $d\bar{r}^{2}$ coefficients yields the
conservation of mass condition that implicitly specifies the shock surface
$\bar{r}=\bar{r}(t)$,

$M(\bar{r})=\frac{4\pi}{3}\rho(t)\bar{r}^{3}.$ (48)

Using this together with (42) and (43) gives the following two relations that
hold at the shock surface:

$\displaystyle\bar{r}$ $\displaystyle=\sqrt{\frac{3\gamma}{\rho(t)}},$ (49)
$\displaystyle\rho$
$\displaystyle=\frac{3}{4\pi}\frac{M}{\bar{r}(t)^{3}}=\frac{3\gamma}{\bar{r}(t)^{2}}=3\bar{\rho}.$
(50)

Matching the $d\bar{t}^{2}$ coefficients on the shock surface determines the
integrating factor $\psi$ (see Section 2) in a neighbourhood of the shock
surface by assigning initial conditions for (12). Finally, the conservation
constraint $[T^{\mu\nu}]n_{\mu}n_{\nu}=0$ leads to the single condition

$(1-A)(\rho+\bar{p})(p+\bar{\rho})^{2}+\left(1-\frac{1}{A}\right)(\bar{\rho}+\bar{p})(\rho+p)^{2}+(p-\bar{p})(\rho-\bar{\rho})^{2}=0,$
(51)

which upon using $p=\sigma\rho$ and $\bar{p}=\bar{\sigma}\bar{\rho}$ is
satisfied providing $\sigma$ and $\bar{\sigma}$ are related by

$\bar{\sigma}=\frac{1}{2}\sqrt{9\sigma^{2}+54\sigma+49}-\frac{3}{2}\sigma-\frac{7}{2}=:H(\sigma).$
(52)

Alternatively, we can solve for $\sigma$ in (52) and write this relation as

$\sigma=\frac{\bar{\sigma}(\bar{\sigma}+7)}{3(1-\bar{\sigma})}.$ (53)

This guarantees that conservation holds across the shock surface. It thus
follows from Theorem 4 that all of the equivalencies in Theorem 3 hold across
the shock surface. Note that $H(0)=0$ and to leading order

$\bar{\sigma}=\frac{3}{7}\sigma+O(\sigma^{2})$ (54)

as $\sigma\rightarrow 0$. Within the physical region
$0\leq\sigma,\bar{\sigma}\leq 1$, $H^{\prime}(\sigma)>0$,
$\bar{\sigma}<\sigma$ and

$\displaystyle H\left(\frac{1}{3}\right)$ $\displaystyle=\sqrt{17}-4\approx
0.1231,$ $\displaystyle H(1)$ $\displaystyle=\sqrt{28}-5\approx 0.2915.$

Using the formulas for the FRW metric in (25)-(27) and setting $R_{0}=1$ at
$\rho=\rho_{0}$, $t=t_{0}$, we obtain the following formulas for the shock
position:

$\displaystyle\bar{r}(t)$ $\displaystyle=\alpha t,$ (55) $\displaystyle r(t)$
$\displaystyle=\frac{\bar{r}(t)}{R(t)}=\beta t^{\frac{1+3\sigma}{3+3\sigma}},$
(56)

where

$\displaystyle\alpha$
$\displaystyle=3(1+\sigma)\sqrt{\frac{\bar{\sigma}}{1+6\bar{\sigma}+\bar{\sigma}^{2}}},$
(57) $\displaystyle\beta$
$\displaystyle=\alpha^{\frac{1+3\sigma}{3+3\sigma}}\left(\frac{3\gamma}{\rho_{0}}\right)^{\frac{1}{3+3\sigma}}.$
(58)

It follows from (43) that $A>0$ and from (56) that $r_{*}=\lim_{t\rightarrow
0}r(t)=0$. The entropy condition that the shock wave be compressive follows
from the fact that $\bar{\sigma}=H(\sigma)<\sigma$. Thus we conclude that for
each $0<\sigma\leq 1$, $r_{*}=0$, the solutions constructed in (41)-(58)
define a one-parameter family of shock-wave solutions that evolve everywhere
outside the black hole, which implies that the distance from the shock wave to
the FRW centre is less than one Hubble length for all $t>0$. Using (55) and
(56), one can determine the shock speed and check when the Lax characteristic
condition [7] holds at the shock. The result is the following (see [12, 1] for
details).999Note that even when the shock speed is larger than $c$, only the
wave, and not the sound speeds or any other physical motion, exceeds the speed
of light.

######  Theorem 5.

Let

$\displaystyle\sigma_{1}$ $\displaystyle=\frac{\sqrt{10}+1}{9}\approx 0.462,$
$\displaystyle\sigma_{2}$ $\displaystyle=\frac{\sqrt{5}}{3}\approx 0.745,$

then the Lax characteristic condition holds at the shock if and only if
$0<\sigma\leq\sigma_{1}$ and the shock speed is subluminal (less than the
speed of light) if and only if $0<\sigma<\sigma_{2}$.

The explicitly defined solution in the case $r_{*}=0$ can be interpreted as
the general relativistic version of a shock-wave explosion into a static
singular isothermal fluid sphere, known in the Newtonian case as a simple
model for star formation [16]. As the scenario goes, a star begins as a
diffuse cloud of gas. The cloud slowly contracts under its own gravitational
force by radiating energy out as gravitational potential energy is converted
into kinetic energy. This contraction continues until the gas cloud reaches
the point where the mean free path for transmission of light is small enough
that light is scattered, instead of transmitted, through the cloud. The
scattering of light within the gas cloud has the effect of equalising the
temperature within the cloud, and at this point the gas begins to drift toward
the most compact configuration of the density that balances the pressure when
the equation of state is isothermal. This configuration is a static singular
isothermal sphere, the general relativistic version of which is the explicitly
given TOV solution beyond the shock wave when $r_{*}=0$. This solution in the
Newtonian case is also inverse square in the density and pressure and so the
density tends to infinity at the centre of the sphere. Eventually, the high
density at the centre ignites a thermonuclear reaction. The result is a shock-
wave explosion emanating from the centre of the sphere, with this explosion
signifying the birth of a star. The solutions when $r_{*}=0$ represent the
general relativistic version of such a shock-wave explosion.

##  5 Shock-Wave Solutions Beyond the Hubble Radius - The Case $r_{*}>0$

When the shock wave is beyond one Hubble length from the FRW centre, we obtain
a family of shock-wave solutions for each $0<\sigma\leq 1$ and $r_{*}>0$ by
matching the FRW metric (1) to a TOV metric of the form (2) to form a shock
wave under the assumption that

$\displaystyle A(\bar{r})=1-\frac{2M(\bar{r})}{\bar{r}}=:1-N(\bar{r})<0.$ (59)

In this case, $\bar{r}$ is the timelike variable. Assuming the stress-energy-
momentum tensor $T$ is taken to be that of a perfect fluid co-moving with the
TOV metric, the Einstein equations $G=\kappa T$ (inside the black hole) take
the form [19],

$\displaystyle\bar{p}^{\prime}$
$\displaystyle=\frac{\bar{p}+\bar{\rho}}{2}\frac{N^{\prime}}{N-1},$ (60)
$\displaystyle N^{\prime}$
$\displaystyle=-\left(\frac{N}{\bar{r}}+\kappa\bar{p}\bar{r}\right),$ (61)
$\displaystyle\frac{B^{\prime}}{B}$
$\displaystyle=-\frac{1}{N-1}\left(\frac{N}{\bar{r}}+\kappa\bar{\rho}\right).$
(62)

The system (60)-(62) defines the simplest class of gravitational metrics that
contain matter, evolve inside the black hole and such that the mass function
$M(\bar{r})<\infty$ at each fixed time $\bar{r}$. System (60)-(62) for $A<0$
differs substantially from the TOV equations for $A>0$ because, for example,
the energy density $T^{00}$ is equated with the timelike component $G^{rr}$
when $A<0$ but with $G^{tt}$ when $A>0$. In particular, this implies that
inside the black hole the mass function $M(\bar{r})$ does not have the
interpretation as the total mass inside radius $\bar{r}$ as it does outside
the black hole.

Equations (61), (62) do not have the same character as (37), (38), and the
relation $\bar{p}=\bar{\sigma}\bar{\rho}$ with constant $\sigma$ is
inconsistent with (61), (62) together with the conservation constraint and the
FRW assumption $p=\sigma\rho$ for shock-wave matching. Thus, instead of
looking for an explicit solution of (61), (62) ahead of time, as in the case
$r_{*}=0$, we assume the FRW solution (25)-(27) and derive the ODE that
describe the TOV metrics that match this FRW metric Lipschitz continuously
across a shock surface and then impose the conservation, entropy and equation
of state constraints at the end. Matching a given flat ($k=0$) FRW metric to a
TOV metric inside the black hole across a shock interface leads to the system
of ODE [19],

$\displaystyle\frac{du}{dN}$
$\displaystyle=-\left(\frac{(1+u)}{2(1+3u)N}\right)\left(\frac{(3u-1)(\sigma-u)N+6u(1+u)}{(\sigma-u)N+(1+u)}\right),$
(63) $\displaystyle\frac{d\bar{r}}{dN}$
$\displaystyle=-\frac{1}{1+3u}\frac{\bar{r}}{N},$ (64)

with conservation constraint

$\displaystyle w=\frac{(\sigma-u)N-\sigma(1+u)}{(\sigma-u)N+(1+u)},$ (65)

where

$\displaystyle u$ $\displaystyle=\frac{\bar{p}}{\rho},$ $\displaystyle w$
$\displaystyle=\frac{\bar{\rho}}{\rho},$ $\displaystyle\sigma$
$\displaystyle=\frac{p}{\rho}.$ (66)

Here $\rho$ and $p$ denote the (known) FRW density and pressure and all
variables are evaluated at the shock. Solutions of (63)-(65) determine the
(unknown) TOV metrics that match the given FRW metric Lipschitz continuously
across a shock interface such that conservation of energy and momentum hold
across the shock and such that there are no delta function sources at the
shock [6, 14]. Note that the dependence of (63)-(65) on the FRW metric is only
through the variable $\sigma$, and so the advantage of taking constant
$\sigma$ is that the whole solution is determined by the inhomogeneous scalar
equation (63) when $\sigma$ is constant. We take as the entropy constraint the
condition that

$\displaystyle 0$ $\displaystyle<\bar{p}<p,$ $\displaystyle 0$
$\displaystyle<\bar{\rho}<\rho,$ (67)

and to insure a physically reasonable solution, we impose the equation of
state constraint on the TOV side of the shock,101010This is equivalent to the
dominant energy condition [2].

$0<\bar{p}<\bar{\rho}.$ (68)

Condition (67) implies that outgoing shock waves are compressive. Inequalities
(67) and (68) are both implied by the single condition [19],

$\frac{1}{N}<\left(\frac{1-u}{1+u}\right)\left(\frac{\sigma-u}{\sigma+u}\right).$
(69)

Since $\sigma$ is constant, equation (63) uncouples from (64), and thus
solutions of system (63)-(65) are determined by the scalar non-autonomous
equation (63). Making the change of variable $S=\frac{1}{N}$, which transforms
the Big Bang $N\rightarrow\infty$ over to a rest point at $S\rightarrow 0$, we
obtain,

$\frac{du}{dS}=\left(\frac{(1+u)}{2(1+3u)S}\right)\left(\frac{(3u-1)(\sigma-u)+6u(1+u)S}{(\sigma-u)+(1+u)S}\right).$
(70)

Note that the conditions $N>1$ and $0<\bar{p}<p$ restrict the domain of (70)
to the region $0<u<\sigma<1$, $0<S<1$. The next theorem gives the existence of
solutions for $0<\sigma\leq 1$, $r_{*}>0$ inside the black hole [18].

######  Theorem 6.

For every $0<\sigma<1$ there exists a unique solution $u_{\sigma}(S)$ of (70)
such that (69) holds for all $0<S<1$. Moreover,

$\displaystyle 0$ $\displaystyle<u_{\sigma}(S)<\bar{u},$
$\displaystyle\lim_{S\rightarrow 0}u_{\sigma}(S)$ $\displaystyle=\bar{u},$
$\displaystyle\lim_{S\rightarrow 1}\bar{p}$ $\displaystyle=0,$
$\displaystyle\lim_{S\rightarrow 1}\bar{\rho}$ $\displaystyle=0,$ (71)

where

$\bar{u}=\min\left\\{\frac{1}{3},\sigma\right\\}.$ (72)

Furthermore, for each of the solutions $u_{\sigma}(S)$, the shock position is
determined by the solution of (64), which in turn is determined uniquely by an
initial condition which can be taken to be the FRW radial position of the
shock wave at the instant of the Big Bang,

$r_{*}=\lim_{S\rightarrow 0}r(S)>0.$ (73)

Concerning the the shock speed, we have the following theorem.

######  Theorem 7.

Let $0<\sigma<1$. Then the shock speed $s_{\sigma}(S):=s(u_{\sigma}(S))<1$ for
all $0<S\leq 1$ if and only if $\sigma<\frac{1}{3}$, that is, the shock speed
is subluminal if and only if $\sigma<\frac{1}{3}$.

For the shock speed near the Big Bang ($S=0$), we have the following theorem:

######  Theorem 8.

The shock speed at the Big Bang ($S=0$) is given by:

$\displaystyle\lim_{S\rightarrow 0}s_{\sigma}(S)=0,$ $\displaystyle\sigma$
$\displaystyle<\frac{1}{3},$ (74) $\displaystyle\lim_{S\rightarrow
0}s_{\sigma}(S)=1,$ $\displaystyle\sigma$ $\displaystyle=\frac{1}{3},$ (75)
$\displaystyle\lim_{S\rightarrow 0}s_{\sigma}(S)=\infty,$
$\displaystyle\sigma$ $\displaystyle>\frac{1}{3}.$ (76)

Theorem 8 shows that the equation of state $p=\frac{1}{3}\rho$ plays a special
role in the analysis when $r_{*}>0$, and only for this equation of state does
the shock wave emerge at the Big Bang at a finite non-zero speed, the speed of
light. Moreover, (72) implies that in this case, the correct relation
$\bar{p}=\bar{\sigma}\bar{\rho}$ is also achieved in the limit $S\rightarrow
0$. The result (71) implies that (neglecting the pressure $p$ at this time
onward) the solution continues to a $k=0$ Oppenheimer–Snyder solution outside
the black hole for $S>1$.

It follows that the shock wave will first become visible at the FRW centre
$\bar{r}=0$ at the moment $t=t_{0}$ (where $R(t_{0})=1$) when the Hubble
length $1/H_{0}:=1/H(t_{0})$ satisfies

$\frac{1}{H_{0}}=\frac{1+3\sigma}{2}r_{*},$ (77)

where $r_{*}$ is the FRW position of the shock at the instant of the Big Bang.
At this time, the number of Hubble lengths $\sqrt{N_{0}}$ from the FRW centre
to the shock wave at time $t=t_{0}$ can be estimated by

$1\leq\frac{2}{1+3\sigma}\leq\sqrt{N}_{0}\leq\frac{2}{1+3\sigma}e^{\sqrt{3\sigma}\left(\frac{1+3\sigma}{1+\sigma}\right)}.$
(78)

Thus, in particular, the shock wave will still lie beyond the Hubble radius
$1/H_{0}$ at the FRW time $t_{0}$ when it first becomes visible. Furthermore,
the time $t_{crit}>t_{0}$ at which the shock wave will emerge from the white
hole, given that $t_{0}$ is the first instant at which the shock becomes
visible at the FRW centre, can be estimated by

$\displaystyle\frac{2}{1+3\sigma}e^{\frac{1}{4}\sigma}$
$\displaystyle\leq\frac{t_{crit}}{t_{0}}\leq\frac{2}{1+3\sigma}e^{\frac{\sqrt{12\sigma}}{1+\sigma}},$
$\displaystyle 0$ $\displaystyle<\sigma<\frac{1}{3},$ (79) $\displaystyle
e^{\frac{\sqrt{6}}{4}}$ $\displaystyle\leq\frac{t_{crit}}{t_{0}}\leq
e^{\frac{3}{2}},$ $\displaystyle\sigma$ $\displaystyle=\frac{1}{3}.$ (80)

Inequalities (79), (80) imply, for example, that at the Oppenheimer–Snyder
limit $\sigma=0$:

$\displaystyle\sqrt{N_{0}}$ $\displaystyle=2,$
$\displaystyle\frac{t_{crit}}{t_{0}}$ $\displaystyle=2$

and in the limit $\sigma=\frac{1}{3}$:

$\displaystyle 1$ $\displaystyle<\sqrt{N_{0}}\leq 4.5,$ $\displaystyle 1.8$
$\displaystyle\leq\frac{t_{crit}}{t_{0}}\leq 4.5.$

We can conclude that the moment $t=t_{0}$, when the shock wave first becomes
visible at the FRW centre, the shock wave must lie within 4.5 Hubble lengths
of the FRW centre. Throughout the expansion up until this time, the expanding
universe must lie entirely within a white hole, that is, the Universe will
eventually emerge from this white hole, but not until some later time
$t_{crit}$, where $t_{crit}$ does not exceed $4.5t_{0}$.

##  6 Self-Similar Extensions of FRW-TOV Shock-Waves Inside the Hubble Radius

The previous two sections focused on constructing general relativistic shock
waves with the explicitly known flat ($k=0$) FRW metric forming the interior
expanding wave. In this section we consider a broader family of expanding
waves and demonstrate how to impose conservation across a shock surface with a
TOV spacetime on the exterior. This achieves three things: The first is the
determination of all expanding waves that are regular at the radial centre
($\bar{r}=0$) that can be matched with conservation to a TOV spacetime inside
the Hubble radius (outside the black hole), the second is the determination of
the unique expanding wave with a pure radiation equation of state
($p=\frac{1}{3}\rho$) that can be matched to the pure radiation TOV spacetime,
and the third is the determination of the accelerated expansion that results
from a general relativistic explosion with a TOV exterior.

Extending the flat FRW spacetime to a family of expanding waves is motivated
by the fact that both the FRW and TOV spacetimes are self-similar in the
variable $\xi=r/t$, that is, the metric components and hydrodynamic variables
depend only on the single variable $\xi$. The self-similarity of the TOV
spacetime means that to match an expanding wave to this spacetime with
conservation across the shock and with an equation of state of the form
$p=\sigma\rho$, it must be the case that the expanding wave also be self-
similar in $\xi$ [3]. Furthermore, requiring that the expanding waves have a
_regular centre_ , leaves only one family of self-similar expanding waves,
with these being the one-parameter family of _asymptotically FRW_ spacetimes
described independently in [4] and [20]. Thus it is the case that the family
of spherically symmetric self-similar perturbations of the FRW metric account
for all physically admissible expanding waves that can be matched to a self-
similar TOV metric with conservation across the shock.

Any self-similar metric may be written, without loss of generality, in the
self-similar Schwarzschild coordinate form

$\displaystyle
ds^{2}=-B(\xi)dt^{2}+\frac{1}{A(\xi)}dr^{2}+r^{2}[d\theta^{2}+\sin^{2}(\theta)d\phi^{2}].$
(81)

Under the assumption of spherical symmetry, the fluid four-velocity may also
be written, without loss of generality, as $\boldsymbol{u}=(u^{0},u^{1},0,0)$.
Under the normalisation condition $g(\boldsymbol{u},\boldsymbol{u})=-1$, the
fluid four-velocity has only one independent component and can thus be fully
specified through the _Schwarzschild coordinate velocity_ , defined by

$\displaystyle v=\frac{1}{\sqrt{AB}}\frac{u^{1}}{u^{0}}.$ (82)

Together with $A$, $B$, $\rho$ and $p$, the Schwarzschild coordinate velocity
$v$ is one of five unknown variables that completely specify a solution to the
self-similar Einstein-Euler equations. As there are only four independent
components of the spherically symmetric Einstein-Euler equations, a barotropic
equation of state of the form $p=p(\rho)$ is used to close the system.
However, spherical symmetry and self-similarity in the variable $\xi$ restrict
this equation of state to the form $p=\sigma\rho$ for constant $\sigma$ [3].
Following the development of Smoller and Temple in [20], substituting the
metric ansatz and equation of state into the Einstein-Euler equations yields
the system of nonlinear ODE:

$\displaystyle\xi\frac{dA}{d\xi}$
$\displaystyle=-\frac{(3+3\sigma)(1-A)v}{\\{\cdot\\}_{S}},$ (83)
$\displaystyle\xi\frac{dG}{d\xi}$
$\displaystyle=-G\left[\left(\frac{1-A}{A}\right)\frac{(3+3\sigma)[(1+v^{2})G-2v]}{2\\{\cdot\\}_{S}}-1\right],$
(84) $\displaystyle\xi\frac{dv}{d\xi}$
$\displaystyle=-\left(\frac{1-v^{2}}{2\\{\cdot\\}_{D}}\right)\left[3\sigma\\{\cdot\\}_{S}+\left(\frac{1-A}{A}\right)\frac{(3+3\sigma)^{2}\\{\cdot\\}_{N}}{4\\{\cdot\\}_{S}}\right],$
(85)

in addition to the constraint

$\displaystyle\rho=\frac{3(1-v^{2})(1-A)G}{\kappa r^{2}\\{\cdot\\}_{S}}.$ (86)

The variable $G$, not to be confused with the Einstein curvature tensor, is
defined by

$\displaystyle G=\frac{\xi}{\sqrt{AB}}$ (87)

and the bracketed terms are given (using the form found in [1]) by:

$\displaystyle\\{\cdot\\}_{S}$ $\displaystyle=3(G-v)-3\sigma v(1-Gv),$
$\displaystyle\\{\cdot\\}_{N}$ $\displaystyle=-3(G-v)^{2}+3\sigma
v^{2}(1-Gv)^{2},$ $\displaystyle\\{\cdot\\}_{D}$
$\displaystyle=\frac{3}{4}(3+3\sigma)\left[(G-v)^{2}-\sigma(1-Gv)^{2}\right].$

Thus a solution is specified fully by the three variables $A$, $G$ and $v$.

The asymptotically FRW spacetimes are described by the family of solutions to
(83)-(86) with the leading order form as $\xi\to 0$ [20]:

$\displaystyle A(\xi)$ $\displaystyle\approx
1-\frac{1}{4}a^{2}\xi^{2}+O(\xi^{4}),$ $\displaystyle G(\xi)$
$\displaystyle\approx\frac{1}{4}(3+3\sigma)\xi+O(\xi^{3}),$ $\displaystyle
v(\xi)$ $\displaystyle\approx\frac{1}{2}\xi+O(\xi^{3}).$

The parameter $a$ is referred to as the _acceleration parameter_ , as changing
this parameter value away from $a=1$, which corresponds to the unperturbed FRW
spacetime, changes the accelerated expansion of the spacetime as measured by
an observer at the radial centre. This change is specified through the
modified red-shift versus luminosity relation

$\displaystyle
d_{l}=2ct_{0}\left(z+\frac{a^{2}-1}{2}z^{2}+\frac{(a^{2}-1)(5a^{2}+4)}{10}z^{3}+|a-1|O(z^{4})\right)$
(88)

where $d_{l}$ is the luminosity distance, $t_{0}$ is the time of observation
of the radiation and $z$ is the redshift factor [20, 1].111111Note that (88)
was originally derived in [20], but due to an error in the original paper this
expression was corrected and reproduced in [1]. It is also important to note
that $a\neq 1$ also breaks the spacial homogeneity of the spacetime, with this
change becoming more apparent for larger $|a-1|$ and less apparent closer to
$\xi=0$.

The asymptotically FRW spacetimes are not known explicitly away from $\xi=0$,
instead these solutions are described by the system of ODE (83)-(85) and
constraint (86). In order to match an asymptotically FRW spacetime to a TOV
spacetime with conservation holding across the shock, the following lemma is
required [1].

######  Lemma 1.

Let $(A,G,v)$ denote a solution to (83)-(85) and suppose there exists a
$\xi_{0}>0$ such that

$\displaystyle A(\xi_{0})=1-2M(\bar{\sigma}).$ (89)

Then $(A,G,v)$ can be matched to the TOV spacetime (with equation of state
$\bar{p}=\bar{\sigma}\bar{\rho}$) on the surface $\xi=\xi_{0}$ and the
Rankine–Hugoniot jump condition is given by

$\displaystyle\frac{[\sigma+v^{2}(\xi_{0})]G(\xi_{0})-(1+\sigma)G^{2}(\xi_{0})v(\xi_{0})}{[1+\sigma
v^{2}(\xi_{0})]G(\xi_{0})-(1+\sigma)v(\xi_{0})}=\bar{\sigma}.$ (90)

Thus if a solution satisfies (89) and (90), then by Lemma 1 this solution is a
shock-wave solution. However, for this shock wave to be stable in the gas-
dynamical sense, that is, for the fluid characteristics to impinge on the
shock surface from both sides, the following theorem is needed [1]. This
theorem establishes the _Lax entropy conditions_ , also known as the _Lax
characteristic conditions_.

######  Theorem 9.

Let $(A,G,v)$ denote a solution to (83)-(85). If there exists a $\xi_{0}>0$
such that $(A,G,v)$ can be matched to the TOV spacetime (with equation of
state $\bar{p}=\bar{\sigma}\bar{\rho}$) to form a shock-wave solution with a
subluminal shock speed $(G(\xi_{0})<1)$, then the Lax characteristic
conditions are satisfied if:

  1. 1.

$\sigma=\bar{\sigma}$, or

  2. 2.

$\sigma<\bar{\sigma}$ and $G(\xi_{0})>\sqrt{\bar{\sigma}}$, or

  3. 3.

$\sigma>\bar{\sigma}$ and $\\{\cdot\\}_{D}(\xi_{0})<0$.

The addition of the acceleration parameter $a$ means it is possible to specify
both the FRW and TOV equation of state independently, with the
Rankine–Hugoniot jump condition then determining $a$. We see from Theorem 9
that if we were to require both equations of state to model pure radiation
($p=\frac{1}{3}\rho$), that is, a uniform equation of state across the shock
surface, then the entropy condition is automatically implied from the
Rankine–Hugoniot jump condition. Indeed, numerical approximations of such a
solution imply $a\approx 2.58$ [1]. Such a solution, if taken as a
cosmological model, would yield an accelerated expansion many orders of
magnitude larger than what is currently observed, since the expected value of
$a$ in such a model would be $a\approx 1$ [20]. However, because this model
uses the TOV spacetime outside the black hole, the shock surface would be
within the Hubble radius and thus expected to be visible at the present time.
It remains an active area of research for Alexander and Temple to construct
shock-wave solutions with shock surfaces beyond the Hubble radius and
determine the accelerated expansion exhibited by these spacetimes. For the
solutions constructed beyond the Hubble radius in Section 5, it is intriguing
that both equations of state tend to $p=\frac{1}{3}\rho$ in the limit of the
Big Bang, further reinforcing the expectation that $a\approx 1$. We finish
with the following theorem from [1].

######  Theorem 10.

There exists an $a>1$ such that an asymptotically FRW spacetime can be matched
to a TOV spacetime within the Hubble radius to form a pure radiation general
relativistic shock wave that satisfies the Lax characteristic conditions.

##  7 Conclusion

We have delved into the mathematics behind general relativistic shock waves
that admit asymptotically FRW spacetimes as the expanding wave behind the
shock. By placing a TOV spacetime on the exterior, these shock waves model the
general relativistic analogue of an explosion within a static singular
isothermal fluid sphere. Whether such modifications to the FRW spacetime could
provide an alternative to the Standard Model of Cosmology first depends on
whether the shock surface lies within the Hubble radius. If it does, such a
shock wave should be presently observable, which is not the case. This means
that the TOV spacetime on the exterior must be modified to be within a black
hole. We have seen that it is possible to construct a shock-wave beyond the
Hubble radius through this modification. Such a model has the notable property
of containing finite total mass but does not account for the accelerated
expansion observed in our Universe today. To introduce an accelerated
expansion, it is necessary to extend the flat FRW spacetime to a family of
self-similar perturbations. This extension provides an additional free
parameter $a$, which allows for the equation of state to be specified
independently each side of the shock, and in particular, allows both equations
of state to model pure radiation. Such a shock wave would thus be applicable
in the Radiation Dominated Epoch of the Early Universe.

Constructing general relativistic shock waves within and beyond the Hubble
radius with the explicitly known flat FRW spacetime behind the shock was
considered in Sections 4 and 5 respectively. The extension of Section 4 to
using self-similar perturbations of the flat FRW spacetime was considered in
Section 6 and permitted the construction of general relativistic shock waves
that induce an accelerated expansion. Section 6 resolved the picture within
the Hubble radius, an area which has been in active research since Cahill and
Taub’s seminal paper on self-similar solutions in General Relativity [3]. The
extension of Section 5 to considering self-similar perturbations of the flat
FRW spacetime behind a shock lying beyond the Hubble radius remains an active
area of research for Alexander and Temple.121212In addition, recent studies of
self-similar expanding waves have led to a new characterisation of the
instability of the flat FRW spacetime in [21]. Alexander and Temple are
currently working on a comprehensive extension of this theory. If the
accelerated expansion induced by this shock wave matches what observations
suggest the rate was during the Radiation Dominated Epoch, then such a model
would offer a mathematically independent mechanism for the accelerated
expansion observed today without the need for a cosmological constant, and
thus, without the need for dark energy.

## References

  * [1] C. Alexander (2022). _General Relativistic Shock Waves that Induce Cosmic Acceleration_. Ph.D. Thesis. University of California, Davis.
  * [2] S. K. Blau and A. H. Guth (1987). _Inflationary Cosmology_. In: _Three Hundred Years of Gravitation_ , 524–603, ed. by S. Hawking and W. Israel. Cambridge University Press.
  * [3] M. E. Cahill and A. H. Taub (1971). _Spherically Symmetric Similarity Solutions of the Einstein Field Equations for a Perfect Fluid_. Communications in Mathematical Physics 21, 1–40.
  * [4] B. J. Carr and A. A. Coley (2000). _Complete Classification of Spherically Symmetric Self-Similar Perfect Fluid Solutions_. Physical Review D 62(4), 044023.
  * [5] J. Groah and B. Temple (2004). _Shock-Wave Solutions of the Einstein Equations with Perfect Fluid Sources: Existence and Consistency by a Locally Inertial Glimm Scheme_. Memoirs of the American Mathematical Society 172(813).
  * [6] W. Israel (1966). _Singular Hypersurfaces and Thin Shells in General Relativity_. Il Nuovo Cimento B (1965-1970) 44(1), 1–14.
  * [7] P. D. Lax (1957). _Hyperbolic Systems of Conservation Laws II_. Communications on Pure and Applied Mathematics 10(4), 537–566.
  * [8] J. R. Oppenheimer and G. M. Volkoff (1939). _On Massive Neutron Cores_. Physical Review 55(4), 374–381.
  * [9] J. R. Oppenheimer and H. Snyder (1939). _On Continued Gravitational Contraction_. Physical Review 56(5), 455–459.
  * [10] J. Smoller (1983). _Shock Waves and Reaction–Diffusion Equations_. Springer Science & Business Media.
  * [11] J. Smoller and B. Temple (1994). _Shock-Wave Solutions of the Einstein Equations: The Oppenheimer–Snyder Model of Gravitational Collapse Extended to the Case of Non-Zero Pressure_. Archive of Rational Mechanics and Analysis 128, 249–297.
  * [12] J. Smoller and B. Temple (1995). _Astrophysical Shock-Wave Solutions of the Einstein Equations_. Physical Review D 51(6), 2733–2743.
  * [13] J. Smoller and B. Temple (1997). _Solutions of the Oppenheimer–Volkoff Equations Inside 9/8’ths of the Schwarzschild Radius_. Communications in Mathematical Physics 184, 597–617.
  * [14] J. Smoller and B. Temple (1997). _General Relativistic Shock-Waves that Extend the Oppenheimer–Snyder Model_. Archive of Rational Mechanics and Analysis 138, 239–277.
  * [15] J. Smoller and B. Temple (1998). _Shock-Wave Solutions in Closed Form and the Oppenheimer–Snyder Limit in General Reality_. SIAM Journal on Applied Mathematics 58(1), 15–33.
  * [16] J. Smoller and B. Temple (1998). _On the Oppenheimer–Volkoff Equations in General Relativity_. Archive of Rational Mechanics and Analysis 142(2), 177–191.
  * [17] J. Smoller and B. Temple (2000). _Cosmology with a Shock Wave_. Communications in Mathematical Physics 210, 275–308.
  * [18] J. Smoller and B. Temple (2003). _Shock-Wave Cosmology Inside a Black Hole_. Proceedings of the National Academy of Sciences 100(20), 11216–11218.
  * [19] J. Smoller and B. Temple (2004). _Cosmology, Black Holes and Shock Waves Beyond the Hubble Length_. Methods and Applications of Analysis 11(1), 77–132.
  * [20] J. Smoller and B. Temple (2012). _General Relativistic Self-Similar Waves that Induce an Anomalous Acceleration into the Standard Model of Cosmology_. Memoirs of the American Mathematical Society 218(1025).
  * [21] J. Smoller, B. Temple and Z. Vogler (2017). _An Instability of the Standard Model of Cosmology Creates the Anomalous Acceleration Without Dark Energy_. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 473(2207), 20160887.
  * [22] B. Temple and J. Smoller (2009). _Expanding Wave Solutions of the Einstein Equations that Induce an Anomalous Acceleration into the Standard Model of Cosmology_. Proceedings of the National Academy of Sciences 106(34), 14213–14218.
  * [23] R. C. Tolman (1939). _Static Solutions of Einstein’s Field Equations for Spheres of Fluid_. Physical Review 55(4), 364–374.
  * [24] S. Weinberg (1972). _Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity_. John Wiley & Sons.


# Extrapolating semileptonic form factors using Bayesian-inference fits
regulated by unitarity and analyticity

J.M. Flynn  A. Jüttner  J.T. Tsang

###### Abstract

We discuss our recently proposed [1] model-independent framework for fitting
hadronic form-factor data, which are often only available at discrete
kinematical points, using parameterisations based on unitarity and
analyticity. The accompanying dispersive bound on the form factors (unitarity
constraint) is used to regulate the ill-posed fitting problem and allow model-
independent predictions over the entire physical range. Kinematical
constraints, for example for the vector and scalar form factors in
semileptonic meson decays, can be imposed exactly. The core formulae are
straight-forward to implement with standard math libraries. We demonstrate the
method for the exclusive semileptonic decay $B_{s}\to K\ell\nu$, an example
requiring one to use a generalisation of the original Boyd Grinstein Lebed
(BGL) unitarity constraint. We further present a first application of the
method to $B\to D^{*}\ell\nu$ decays.

##  1 Introduction

In the study of semileptonic $B_{(s)}$-meson decays on the lattice one still
faces the difficulty of reconciling all relevant physical scales at the same
time, while keeping systematic uncertainties at bay.111For a recent review of
$B_{(s)}$-physics from lattice QCD we refer the reader to Ref. [2]. For
instance, the heavy mass of the $b$ quark and the requirement to induce large
spatial final-state momenta both require fine lattice spacings. There is also
the requirement of large physical lattice volumes to limit finite-size
effects. Besides these regulator-dependent limitations, there is also the
problem of the deteriorating signal-to-noise ratio when including data for
larger final-state momenta, data which is required to cover a larger fraction
of the kinematically allowed momentum transfer between the initial and final
state hadrons.

While accommodating all the above constraints remains a long-term goal, a
widely used strategy in the meantime is to create lattice data for
(near-)physical simulation parameters (in particular quark masses), and to
restrict predictions to relatively large momentum transfer $q^{2}$. In order
to make predictions for hadronic form factors over the entire physical
semileptonic region one then requires reliable and model-independent
extrapolation methods constrained by the available lattice data and any
further input quantum-field theory provides. The ideas presented here have
previously been published in Ref. [1] and applied in Ref. [3].

##  2 Fit ansatz

Boyd, Grinstein and Lebed (BGL) proposed one such model-independent ansatz
[4],

$f_{X}(q_{i}^{2})=\frac{1}{B_{X}(q_{i}^{2})\phi_{X}(q_{i}^{2},t_{0})}\sum_{n=0}^{K_{X}-1}a_{X,n}z(q^{2}_{i})^{n}=\sum_{n=0}^{K_{X}-1}Z_{XX,in}a_{X,n},$
(1)

with the unitarity constraint for the coefficients, $|{\bf a}_{X}|^{2}\leq 1$,
derived from dispersion theory. $\phi_{X}(q^{2},t_{0})$ is a known “outer
function” and the Blaschke factor $B_{X}(q^{2})$ is chosen to vanish at the
positions of sub-threshold poles $M^{X}_{i}$. Similar ideas underlie the
approaches in for example [5, 6]. The subscript $X$ specifies the form factor,
e.g. $X=+,0$ for the vector and scalar form factor in tree-level pseudoscalar-
to-pseudoscalar decay $B_{s}\to K\ell\nu$. To the very right of the equation
we introduce a vector-matrix notation for the ansatz. The objective is then to
determine the coefficients $a_{X,n}$ from a finite number $N_{\rm data}$ of
experimental or theory data points.

Within a frequentist fitting strategy the constraint on the number of degrees
of freedom, $N_{\rm dof}=N_{\rm data}-K_{X}\geq 1$, often limits one’s ability
to estimate the truncation error reliably. Moreover, a meaningful
implementation and interpretation of the unitarity constraint within the
Frequentist framework is not straight forward. Here instead we propose to fit
the parameterisation using Bayesian inference. This provides a conceptually
clean way to implement the unitarity constraint and also to determine
coefficients of the parameterisation beyond the Frequentist bound on $N_{\rm
dof}$. In particular, we propose to use the unitarity constraint as a
regulator for the higher-order coefficients. In essence, the unitarity
constraint forces all coefficients to lie within a limited $K_{X}$-dimensional
space. Together with the fact that $|z|\leq 1$ in Eq. (1), this leads to a
suppression of higher-order terms, allowing us to probe from which point on
the statistical error dominates over the truncation error.

##  3 An aside on the unitarity constraint

Figure 1: Illustration of the mapping $q^{2}\to z(q^{2})$, showing the
physical semileptonic region (magenta), poles (crosses) and branch cuts (light
and dark green).

In Fig. 1 we illustrate how the real $q^{2}$ axis maps onto the complex unit
disk where the $z$ parameter in Eq. (1) is defined. The physical semileptonic
region is mapped to a range on the real $z$ axis around zero, and the branch
cut is mapped onto the unit circle. For the $b\to u$ transition in
$B\to\pi\ell\nu$ and $B_{s}\to K\ell\nu$ there are two qualitatively different
situations: for $B\to\pi\ell\nu$ the relevant two-particle production
threshold and lower limit of integration for the corresponding unitarity
constraint is at $q^{2}=(M_{B}+M_{\pi})^{2}$, while the relevant threshold for
$B_{s}\to K\ell\nu$ is at $(M_{B_{s}}+M_{K})^{2}$, which is the part
corresponding to the arc $[-\alpha_{B_{s}K},+\alpha_{B_{s}K}]$ in the
$z$-plane, as indicated in the figure (see also [7, 8, 9]). As a result, the
corresponding unitarity constraints are

$\frac{1}{2\pi
i}\oint_{C}\frac{dz}{z}\theta_{z}|B_{X}(q^{2})\phi_{X}(q^{2},t_{0})f_{X}(q^{2})|^{2}\leq
1\,,$ (2)

where the step function $\theta_{z}=\theta(\alpha_{B_{s}K}-|{\rm arg}[z]|)$ in
the $B_{s}\to K$ case restricts the integration over the unit circle to the
relevant segment, _i.e._ the one corresponding to the branch cut above the
$B_{s}K$ threshold $t_{+}$. Inserting the BGL expansion Eq. (1), the unitarity
constraint takes the compact form [1]

$\sum\limits_{i,j\geq 0}a_{X,i}^{\ast}\langle
z^{i}|z^{j}\rangle_{\alpha_{B_{s}K}}a_{X,j}\equiv|{\bf
a}_{X}|^{2}_{\alpha_{B_{s}K}}\leq 1\,,$ (3)

where the inner product is known analytically,

$\langle
z^{i}|z^{j}\rangle_{\alpha}=\frac{1}{2\pi}\int\limits_{-\alpha}^{\alpha}d\phi(z^{i})^{\ast}z^{j}|_{z=e^{i\phi}}=\begin{cases}\displaystyle\frac{\sin(\alpha(i-j))}{\pi(i-j)}&i\neq
j\,,\\\ \displaystyle\frac{\alpha}{\pi}&i=j\,.\end{cases}$ (4)

##  4 Frequentist fit

We begin with the discussion of fits to lattice results for the $B_{s}\to\
K\ell\nu$ scalar $(f_{0}(q_{i}^{2}))$ and vector $(f_{+}(q_{i}^{2}))$ form
factors, where we assume that we have data for $N_{0}$ and $N_{+}$ data points
respectively. We collect all results in the data vector

${\bf f}^{\,T}=({\bf f}_{+}^{T},{\bf
f}_{0}^{T})=(f_{+}(q^{2}_{0}),f_{+}(q^{2}_{1}),\dots,f_{+}(q^{2}_{N_{+}-1}),f_{0}(q^{2}_{0}),f_{0}(q^{2}_{1}),\dots,f_{0}(q^{2}_{N_{0}-1}))\,.$
(5)

Correspondingly, we collect the BGL parameters into the parameter vector

${\bf a}^{T}=({\bf a}_{+}^{T},{\bf
a}_{0}^{T})=(a_{+,0},a_{+,1},a_{+,2},...,a_{+,K_{+}-1},a_{0,1},...,a_{0,K_{0}-1})\,.$
(6)

With this notation the Frequentist fit is defined as the minimisation of

$\chi^{2}({\bf a},{\bf f})=\left[{\bf f}-Z{\bf a}\right]^{T}C_{\bf
f}^{-1}\left[{\bf f}-Z{\bf a}\right]\,,$ (7)

with respect to the parameters, where $C_{\bf f}$ is the correlation matrix of
the input data ${\bf f}$. Details on how the kinematical constraint
$f_{+}(0)=f_{0}(0)$ can be implemented through the matrix $Z_{XX}$ can be
found in Ref. [1]. Note that we chose to use this constraint to eliminate the
parameter $a_{0,0}$, which is therefore missing from the definition of ${\bf
a}$ in Eq. (6). The solution for the Frequentist fit is

${\bf a}=\left(Z^{T}C_{\bf f}^{-1}Z\right)^{-1}ZC_{\bf f}^{-1}{\bf f}\,,\qquad
C_{\bf a}=\left(Z^{T}C_{\bf f}^{-1}Z\right)^{-1}\,,$ (8)

where $C_{\bf a}$ is the covariance matrix of the fit parameters ${\bf a}$.

##  5 Bayesian fit

Within the Bayesian framework the fit is defined in terms of the expectation
value

$\langle g({\bf a})\rangle=\mathcal{N}\int d{\bf a}\,g({\bf a})\,\pi({\bf
a}|{\bf f},C_{\bf f})\,\pi_{\bf a}\,,$ (9)

where $g({\bf a})$ is a function defined in terms of the BGL parameters. The
probability density for the integral is given as

$\pi_{\bf a}({\bf a}|{\bf f},C_{{\bf f}_{p}})\pi_{\bf a}\propto{\rm
exp}\left(-\frac{1}{2}\chi^{2}({\bf a},{\bf f}_{p})\right)\theta(1-|{\bf
a}_{+}|_{\alpha}^{2})\theta(1-|{\bf a}_{0}|_{\alpha}^{2})\,,$ (10)

where the two Heaviside functions restrict the integration to parameters
compatible with the unitarity constraint. In Ref. [1] we explain how to
estimate the integral in Eq. (9) by means of sampling from a multivariate
normal distribution.

##  6 Results for $B_{s}\to K\ell\nu$

Table 1 shows, as an example, the results for Frequentist fits to the HPQCD 14
data of [10]. The first two columns indicate the order of the fit in a given
row for the vector and scalar form factor, respectively. Acceptable fits (see
$p$-values in the 3rd-last column) are achieved only for $K_{+}\geq 3$ and
$K_{0}\geq 2$. Based on the results in this table, which shows various
possible fits with $N_{\rm dof}\geq 1$, solid conclusions about convergence
cannot be drawn. Tab. 2 shows the results of the Bayesian fit with
$K_{+,0}\leq 8$, where the convergence of the fit parameters is clearly
visible. The higher-order coefficients are constrained by the unitarity
constraint, which acts as regulator – without this the parameters would not be
well constrained.

$K_{+}$ | $K_{0}$ | $a_{0,0}$ | $a_{0,1}$ | $a_{0,2}$ | $a_{+,0}$ | $a_{+,1}$ | $a_{+,2}$ | $p$ | $\chi^{2}_{\rm red}$ | $N_{\rm dof}$
---|---|---|---|---|---|---|---|---|---|---
2 | 2 | 0.0883(44) | -0.250(17) | - | 0.0270(13) | -0.0792(50) | - | 0.03 | 2.93 | 3
2 | 3 | 0.0880(44) | -0.242(19) | 0.053(65) | 0.0273(13) | -0.0760(63) | - | 0.02 | 4.06 | 2
3 | 2 | 0.0906(45) | -0.240(17) | - | 0.0257(14) | -0.0805(50) | 0.068(31) | 0.15 | 1.89 | 2
3 | 3 | 0.0908(46) | -0.215(22) | 0.138(71) | 0.0262(14) | -0.0727(64) | 0.096(34) | 0.97 | 0.00 | 1
Table 1: Frequentist fit results for HPQCD 14 data, where $\chi^{2}_{\rm red}=\chi^{2}/N_{\rm dof}$. $K_{+}$ | $K_{0}$ | $a_{0,0}$ | $a_{0,1}$ | $a_{0,2}$ | $a_{0,3}$ | $a_{0,4}$ | $a_{0,5}$ | $a_{0,6}$ | $a_{0,7}$ |
---|---|---|---|---|---|---|---|---|---|---
2 | 2 | 0.0883(44) | -0.250(17) | - | - | - | - | - | - |
2 | 3 | 0.0880(44) | -0.243(19) | 0.052(65) | - | - | - | - | - |
3 | 2 | 0.0907(46) | -0.240(17) | - | - | - | - | - | - |
3 | 3 | 0.0906(44) | -0.215(22) | 0.137(73) | - | - | - | - | - |
3 | 4 | 0.0907(47) | -0.215(22) | 0.14(11) | -0.01(31) | - | - | - | - |
4 | 3 | 0.0907(45) | -0.214(22) | 0.139(72) | - | - | - | - | - |
4 | 4 | 0.0907(46) | -0.215(25) | 0.12(19) | -0.08(60) | - | - | - | - |
5 | 5 | 0.0909(46) | -0.218(25) | 0.10(19) | -0.12(55) | 0.04(63) | - | - | - |
6 | 6 | 0.0907(45) | -0.217(25) | 0.10(19) | -0.11(53) | 0.06(66) | -0.02(66) | - | - |
7 | 7 | 0.0907(46) | -0.217(26) | 0.11(20) | -0.08(51) | 0.03(73) | 0.03(81) | -0.04(70) | - |
8 | 8 | 0.0908(46) | -0.217(25) | 0.11(20) | -0.08(50) | -0.01(84) | 0.1(1.0) | -0.09(96) | 0.08(74) |
Table 2: Results for BGL coefficients for the scalar form factor (results for
the vector form factor can be found in [1]).

While the Frequentist fit can provide important information about the
compatibility of the fit function with the data, the Bayesian fit allows for a
meaningful imposition of the unitarity constraint and also to study the
convergence of the BGL parameterisation. We think that future analyses of
form-factor data should take advantage of this complementarity.

##  7 $B\to D^{\ast}\ell\nu$ fit

The method can be extended to other decay channels such as the semileptonic
decay of a pseudo-scalar to a vector particle. We consider the case of $B\to
D^{\ast}\ell\nu$ for which lattice data exists from FNAL/MILC [11], HPQCD [12]
and JLQCD [13].

Figure 2: Result for BGL fit to JLQCD data [13] for the $B\to D^{\ast}\ell\nu$
form factor $\mathcal{F}_{1}$. Fit without (left) and with (right) unitarity
constraint.

Fig. 2 shows the result for the form factor $\mathcal{F}_{1}$ from a
simultaneous fit to the JLQCD data over all four form factors
$f,\,g,\,\mathcal{F}_{1},\,\mathcal{F}_{2}$, including kinematical constraints
between $\mathcal{F}_{1}$ and $f$, and $\mathcal{F}_{1}$ and
$\mathcal{F}_{2}$, respectively. The plots show the fit results once without
imposing unitarity and once with. The result of imposing unitarity is a
substantial reduction in the statistical error in the extrapolation of the
lattice data towards zero momentum transfer (towards the right-hand side of
the plots).

##  8 Conclusions

We presented novel ideas for fitting model- and truncation-independent
parametersations to data for semileptonic decay form factors. The combination
of information from Frequentist and Bayesian fits allows for a comprehensive
understanding of fit-quality and truncation dependence. Moreover, the Bayesian
framework allows for a meaningful imposition of the unitarity constraint to
regulate the determination of fit parameters and, as demonstrated for the case
of $B\to D^{\ast}\ell\nu$, to reduce the statistical error in form-factor
extrapolations. The code underlying the results of [1] is publicly available
under [14].

## References

  * [1] J. M. Flynn, A. Jüttner, and J. T. Tsang, Bayesian inference for form-factor fits regulated by unitarity and analyticity, arXiv:2303.11285.
  * [2] J. T. Tsang and M. Della Morte, B-physics from Lattice Gauge Theory, arXiv:2310.02705.
  * [3] RBC/UKQCD Collaboration, J. M. Flynn, R. C. Hill, A. Jüttner, A. Soni, J. T. Tsang, and O. Witzel, Exclusive semileptonic $B_{s}$→K$\ell$$\nu$ decays on the lattice, Phys. Rev. D 107 (2023), no. 11 114512, [arXiv:2303.11280].
  * [4] C. G. Boyd, B. Grinstein, and R. F. Lebed, Constraints on form-factors for exclusive semileptonic heavy to light meson decays, Phys. Rev. Lett. 74 (1995) 4603–4606, [hep-ph/9412324].
  * [5] I. Caprini, L. Lellouch, and M. Neubert, Dispersive bounds on the shape of $\bar{B}\to D^{(*)}\ell\bar{\nu}$ form-factors, Nucl. Phys. B530 (1998) 153–181, [hep-ph/9712417].
  * [6] C. Bourrely, I. Caprini, and L. Lellouch, Model-independent description of $B\to\pi l\nu$ decays and a determination of $|V_{ub}|$, Phys. Rev. D79 (2009) 013008, [arXiv:0807.2722]. erratum: Phys. Rev. D82 (2010) 099902.
  * [7] N. Gubernari, D. van Dyk, and J. Virto, Non-local matrix elements in $B_{(s)}\to\\{K^{(*)},\phi\\}\ell^{+}\ell^{-}$, JHEP 02 (2021) 088, [arXiv:2011.09813].
  * [8] N. Gubernari, M. Reboud, D. van Dyk, and J. Virto, Improved theory predictions and global analysis of exclusive $b\to s{\mu}^{+}{\mu}^{-}$ processes, JHEP 09 (2022) 133, [arXiv:2206.03797].
  * [9] T. Blake, S. Meinel, M. Rahimi, and D. van Dyk, Dispersive bounds for local form factors in $\Lambda_{b}\to\Lambda$ transitions, arXiv:2205.06041.
  * [10] C. Bouchard, G. P. Lepage, C. Monahan, H. Na, and J. Shigemitsu, $B_{s}\to K\ell\nu$ form factors from lattice QCD, Phys.Rev. D90 (2014), no. 5 054506, [arXiv:1406.2279].
  * [11] Fermilab Lattice, MILC, Fermilab Lattice, MILC Collaboration, A. Bazavov et al., Semileptonic form factors for $B\rightarrow D^{*}\ell\nu$ at nonzero recoil from $2+1$-flavor lattice QCD: Fermilab Lattice and MILC Collaborations, Eur. Phys. J. C 82 (2022), no. 12 1141, [arXiv:2105.14019]. [Erratum: Eur.Phys.J.C 83, 21 (2023)].
  * [12] J. Harrison and C. T. H. Davies, $B\rightarrow D^{*}$ vector, axial-vector and tensor form factors for the full $q^{2}$ range from lattice QCD, arXiv:2304.03137.
  * [13] JLQCD Collaboration, Y. Aoki, B. Colquhoun, H. Fukaya, S. Hashimoto, T. Kaneko, R. Kellermann, J. Koponen, and E. Kou, $B\to D^{*}\ell\nu_{\ell}$ semileptonic form factors from lattice QCD with Möbius domain-wall quarks, arXiv:2306.05657.
  * [14] A. Jüttner, BFF – Bayesian Form factor Fit code https://github.com/andreasjuettner/BFF, https://doi.org/10.5281/zenodo.7799543, .