#  YITP-21-142 Gravitational collapse of spherical shells of fluid in the
isotropic homogeneous universe

Shuichi Yokoyama,
∗ Center for Gravitational Physics,
Yukawa Institute for Theoretical Physics, Kyoto University,
Kitashirakawa Oiwake-cho, Sakyo-Ku, Kyoto 606-8502, Japan
shuichi.yokoyama[at]yukawa.kyoto-u.ac.jp

We investigate gravitational collapse of thick shell of fluid in the isotropic
homogeneous universe without radiation described by the Einstein gravity with
cosmological constant. We construct analytic solutions of this kind
interpolating the Friedmann-Lemaitre-Robertson-Walker (FLRW) metric and the
de-Sitter Schwarzschild black hole one by the Gullstrand-Painlevé one. After
determining the scale factor for a perfect fluid coexisting with cosmological
constant, we determine the orbit of the collapsing shells of fluid from the
continuity condition for density and the perpendicular component of pressure
at interface. However the continuity condition cannot fix the orbit when there
does not exist a perfect fluid inside the shells. In this case we determine
the orbit from an equation of state for the fluid consisting a shell near the
surface. Finally we confirm that the total energy defined in arXiv:2005.13233
is independent of the given time evolution in this system.

##  1 Introduction

The observation of a supermassive black hole and its event horizon [1, 2] has
made its existence manifest and endorsed that laws of gravitational physics
are governed by General Relativity. A typical formation of a black hole is due
to gravitational collapse, by which a compact object becomes too heavy to
support itself by the internal repulsive force by the matter constituents.
Therefore the investigation of gravitational collapse is highly motivated and
it has been one of the traditional subjects in gravitational physics from an
early period [3, 4, 5, 6, 7, 8, 9, 10]. (See [11] for a review and references
therein.)

Another telescopic observation of high redshift supernova [12, 13] and cosmic
microwave background [14, 15] indicates the existence of nonzero cosmological
constant, whose value may be much smaller than the scale expected from
particle physics [16]. By taking into account this result, it is motivated to
investigate gravitational collapse in the universe with cosmological constant.
There have been several works on gravitational collapse in the Einstein
gravity with cosmological constant, for instance, [17, 18, 19, 20, 21, 22, 23,
24], which were all studied by the so-called junction method [25, 26, 27, 28].
The junction method provides a general way to construct a metric to describe
gravitational collapse by preparing two coordinate patches for the spacetime
and gluing them smoothly at their overlap, though phenomena described by the
glued metric may not be lucid as a whole.

On the other hand, there is another approach to gravitational collapse by
preparing for one good coordinate system to cover the entire space-time in
which gravitational collapse happens [29, 30, 31]. For convenience we refer to
this approach as the interpolating method. This method may be more refined
than the junction method with regard that one needs to find out a good
coordinate system which will smoothly interpolate the two coordinate systems
glued by the junction method. In order to take the latter approach to studying
gravitational collapse in our universe, one needs to prepare for a good
coordinate system which interpolates the Friedmann-Lemaitre-Robertson-Walker
(FLRW) metric and the de-Sitter Schwarzschild black hole metric. It was
pointed out in [31] that such a coordinate system is given by the Gullstrand-
Painlevé one [32, 33]. The Gullstrand-Painlevé coordinate system describes
space-time like flows of a river flowing through a flat background [34], and
it is useful to describe gravitational collapse of our interest. (See also
[35].)

The aim of this paper is to construct an analytic solution of a catastrophic
gravitational collapse of thick shells of fluid into a black hole as described
in the figure drawn in [6] in the isotropic homogeneous universe governed by
the Einstein gravity with cosmological constant adopting the interpolating
method. This model has a richer structure than the one studied in [36] in the
sense that there exists a perfect fluid inside the collapsing or blowing
shells of fluid. (See also [37, 38].) The construction of such a dynamical
solution with fluid collapsing into a black hole is indeed nontrivial since it
was shown that an emergent curvature singularity is not generally covered by
an event horizon and the cosmic censorship could be violated locally or
globally in a collapsing process [39, 40, 41, 42]. (See [43, 11] for reviews
and further references.) For reader’s convenience, we illustrate the possible
local causal structures of a gravitationally collapsing model of fluid
considered in this paper by Penrose-Carter diagrams in Fig. 1. We consider a
situation where there exists a positive cosmological constant, which generates
the de Sitter horizon and changes the far asymptotic structure compared to the
flat case. Therefore the Penrose-Carter diagrams for the local structures are
modified compared to those given in [39].

(a) Emergent black hole (b) Emergent local naked singularity 1 (c) Emergent
local naked singularity 2 (d) Emergent global naked singularity

Figure 1: The Penrose-Carter diagrams to describe possible local causal
structures of the gravitational collapse of thick shell of fluid are shown.
“HR”, “Future $\infty$” mean “horizon”, “Future infinity”, respectively. The
bold line emanating from the singularity to the event horizon is the infinite
redshift surface or the trapping horizon. The space-time is separated by the
thick shell of fluid, whose region is denoted by II and covered by some color.
The region inside, which is colored by another color and labeled by I, is
filled with a perfect fluid, while in the region outside denoted by III there
exists only cosmological constant or dark energy.

For the purpose of constructing a completely collapsing solution without a
naked singularity, the strategy is to set an ansatz for an interpolating
metric connecting to a static black hole solution as the final state and
determine an unfixed parameter in the ansatz so as to satisfy the Einstein
equation. We determine the orbit of the shells of fluid by the continuity
condition for density and the perpendicular component of pressure. It turns
out that this gravitationally collapsing or blowing model is correlated to the
status of the FLRW universe. That is, in this model, a gravitationally
collapsing solution happens only when the region inside described by the FLRW
metric is under accelerating shrinkage, while a gravitationally blowing one
occurs only when the region inside is under accelerating expansion. What we
call a gravitationally blowing solution is obtained by a time reversal for a
gravitationally collapsing one, and thus such a solution describes the process
of a classical evaporation of a white hole. We call this solution
gravitational evaporation.

As a special case, we can consider a situation where there does not exist a
perfect fluid inside the shell. In this situation, the region inside the shell
expands or shrink exponentially, since there exists only cosmological constant
or dark energy inside. Then the orbit of the shell is not determined only by
the continuity condition. This is natural because the property of the fluid
consisting of shells is not determined in this situation. There is room to
specify the property of the shells, and, we determine the orbit of the shells
of the fluid by imposing an equation of state for fluid at a layer near the
surface.

Finally we compute the energy distribution inside a fixed radius by adopting
the manifest covariant definition of energy presented in [44, 45, 46]. From
this energy distribution one can see that the energy of a perfect fluid is
converted into that of the shells of fluid in the region where the shells
exist. Using this definition we easily check that the total energy of the
system is independent of the given time evolution and matches the quasi-local
mass evaluated at the behavior of the metric in the far asymptotic region.

The rest of this paper is organized as follows. In section 2 we explain that
the Gullstrand-Painlevé coordinate system suits to describing gravitational
collapse of thick shells of fluid in the isotropic homogeneous universe. In
section 3 we determine the scale factor for a perfect fluid together with
cosmological constant and investigate the continuity condition for physical
quantities between the region of a perfect fluid and that of shells of fluid,
from which the orbit of the shells is determined. In section 4 we compute the
energy distribution of the system and confirm the conservation of the total
energy. Section 5 is devoted to summary and discussion. In appendix A we
determine the orbit of shells of fluid in the absence of a perfect fluid
inside by imposing an equation of state near the surface.

##  2 Gullstrand-Painlevé coordinates as interpolating metric

We investigate the spherically symmetric gravitational collapse into a black
hole or the gravitational evaporation from a white hole in the homogeneous
isotropic universe, which is described by the FLRW metric with conformally
flat space as

$\displaystyle ds_{\rm FLRW}^{2}=$
$\displaystyle-(dx^{0})^{2}+a(x^{0})^{2}({d\check{r}^{2}}+\check{r}^{2}\tilde{g}_{ij}dx^{i}dx^{j}),$
(2.1)

where the scale factor $a(x^{0})>0$ is only dependent on the time coordinate
$x^{0}$, and $\tilde{g}_{ij}$ is a metric for the unit $(d-2)$-dimensional
Einstein manifold with $\tilde{R}_{ij}=(d-3)\tilde{g}_{ij}$. The behavior of
the scale factor is determined by a state of a perfect fluid filled in the
region described by the metric.

In order to construct a solution for gravitational collapse or evaporation in
the homogeneous isotropic universe, what we first do is to find a metric
interpolating between the FLRW one (2.1) and the Kottler spacetime [47, 48]
given by

$\displaystyle ds_{\rm Kott}^{2}=$
$\displaystyle-(1+u(r))(dx_{s}^{0})^{2}+{1\over(1+u(r))}dr^{2}+r^{2}\tilde{g}_{ij}dx^{i}dx^{j},$
(2.2)

where $u(r)=-\frac{2\Lambda}{(d-2)(d-1)}r^{2}-\frac{m}{r^{d-3}}$ with $m_{0}$
a positive parameter related to the black hole mass and $\Lambda$ the
cosmological constant. It was pointed out in [31] that such an interpolation
can be nicely done by the Gullstrand-Painlevé coordinates [32, 33]:

$\displaystyle ds_{\rm GP}^{2}=-(1-\psi_{\pm}^{2})(dt_{\pm})^{2}\pm
2\psi_{\pm}dt_{\pm}dr+dr^{2}+r^{2}\tilde{g}_{ij}dx^{i}dx^{j},~{}~{}$ (2.3)

where $\psi_{\pm}$ generally depends on the time and the radial coordinates.
We call the coordinate system with the upper sign chosen the upper patch
denoted by $U_{+}$ and the one with the lower sign the lower patch by $U_{-}$
for convenience. These two local patches are not transformed to each other for
general $\psi_{\pm}$, but they are if the two functions have the same absolute
value, $\psi_{+}=\pm\psi_{-}$, in which case the coordinate transformation is
simply $t_{+}=\mp t_{-}$ with the double signs in the same order. In what
follows, for notational simplicity, we write $t=\pm t_{\pm}$,
$\psi=\psi_{\pm}$, so that

$\displaystyle ds_{\rm GP}^{2}=-(1-\psi^{2})(dt)^{2}+2\psi
dtdr+dr^{2}+r^{2}\tilde{g}_{ij}dx^{i}dx^{j}.$ (2.4)

To see the interpolation, we change from the comoving radial coordinate
$\check{r}$ to the proper one by performing a coordinate transformation for
the FLRW metric (2.1) to absorb the scale factor into the radial coordinate by
$r=a(x^{0})\check{r},t=x^{0}$ with other coordinates intact. Then the above
FLRW metric is expressed by the Gullstrand-Painlevé coordinates (2.4) with
$\psi=-H(t)r$ and $H(t)={d\log a(t)\over dt}$. Remark that $\psi$ has the
opposite sign to $H(t)$, so that $\psi$ is negative for accelerating
expansion, while $\psi$ is positive for accelerating shrinkage. On the other
hand, to see the Gullstrand-Painlevé coordinates to admit the Kottler black
hole metric, we perform a coordinate transformation such that

$t=x^{0}_{s}-g(r),~{}~{}~{}g(r)=\int dr{-\psi(r)\over 1+\psi(r)^{2}},$ (2.5)

where we assume that $\psi$ is a function only of the radial coordinate $r$.
Then the Gullstrand-Painlevé coordinates (2.4) describe the Kottler metric
(2.2) with $\psi=-\sqrt{-u(r)}$. Notice that this transformation is meaningful
as long as the function $u(r)$ as well as $\psi$ are nonpositive for all
$r>0$. Therefore, firstly in order to study gravitational collapse or
evaporation described by the Gullstrand-Painlevé coordinates, we need to
assume that the cosmological constant is nonnegative. This means that there
exists a nonnegative constant $H_{0}$ such that

$H_{0}^{2}={2\Lambda\over(d-1)(d-2)}.$ (2.6)

Secondly, the region of the Kottler black hole and the one described by the
FLRW metric cannot be covered by only one of the coordinate patches $U_{\pm}$,
which will be seen in the next section.

Now that the Gullstrand-Painlevé coordinates (2.4) suit to describe a
gravitational collapse or evaporation in the FLRW universe, we choose it as
the ansatz to solve the Einstein equation. By computing the Einstein tensor
from the Gullstrand-Painlevé coordinates (2.4), we can determine the form of
the stress energy tensor to satisfy the Einstein equation
$T_{\mu\nu}=\frac{1}{8\pi G_{N}}(R_{\mu\nu}-{1\over 2}Rg_{\mu\nu}+\Lambda
g_{\mu\nu})$, where $G_{N}$ is the Newton constant:

$\displaystyle T^{t}\\!_{t}=$ $\displaystyle{d-2\over 16\pi
G_{N}}{(r^{d-3}(-\psi^{2}+H_{0}^{2}r^{2}))^{\prime}\over r^{d-2}},~{}~{}$
$\displaystyle T^{t}\\!_{r}=$ $\displaystyle{d-2\over 16\pi
G_{N}}{2\psi\dot{\psi}\over r},~{}~{}~{}~{}T^{r}\\!_{t}=0,~{}~{}$
$\displaystyle T^{r}\\!_{r}=$ $\displaystyle{d-2\over 16\pi
G_{N}}\left({(r^{d-3}(-\psi^{2}+H_{0}^{2}r^{2}))^{\prime}\over
r^{d-2}}+\frac{2\dot{\psi}}{r}\right),$ $\displaystyle T^{i}\\!_{j}=$
$\displaystyle{\delta^{i}_{j}\over 16\pi
G_{N}}\left(\frac{(r^{d-3}(-\psi^{2}+H_{0}^{2}r^{2}))^{\prime\prime}}{r^{d-3}}+2\frac{(r^{d-3}\dot{\psi})^{\prime}}{r^{d-3}}\right),$
(2.7)

where $\dot{F}:=\partial_{t}F,F^{\prime}:=\partial_{r}F$. We have checked that
this stress energy tensor satisfies the covariant conservation equation at
$d=4$. The stress energy tensor is related to macroscopic observable
quantities for the infalling fluid as

$T^{\mu}\\!_{\nu}=\rho
u^{\mu}u_{\nu}+q^{\mu}u_{\nu}+u^{\mu}q_{\nu}+P^{\mu}\\!_{\nu},$ (2.8)

where $u^{\mu}$ is the fluid velocity relative to a radially moving frame,
$\rho$ is the density, $q^{\mu},P^{\mu}\\!_{\nu}$ are a vector and a tensor
perpendicular to $u^{\mu}$, respectively. Note that due to the
perpendicularity the tensor $P^{\mu}\\!_{\nu}$ can be written as

$P^{a}\\!_{b}=p_{\perp}(\delta^{a}\\!_{b}+u^{a}u_{b}),~{}~{}P^{i}\\!_{j}=p_{\parallel}\delta^{i}\\!_{j},$
(2.9)

where $a,b=t,r$ and $p_{\perp},p_{\parallel}$ are the radial, angular
component of the pressure, respectively. We fix the radial moving frame so
that the fluid velocity be well-defined during the whole dynamical process.
From a radially moving frame the fluid velocity is written as
$(u^{\mu})=(1,u^{r},\vec{0})$, where the time component is normalized to be
unity. Furthermore we consider a perfect fluid, which is generally timelike, a
fluid velocity is normalized so that $u^{\mu}u_{\mu}=-1$. From these two
equations the radial component is uniquely fixed as $u^{r}=-\psi$. Then the
physical quantities read

$\begin{split}\rho=&-{d-2\over 16\pi
G_{N}}{(r^{d-3}(-\psi^{2}+H_{0}^{2}r^{2}))^{\prime}\over r^{d-2}},\\\
p_{\perp}=&{d-2\over 16\pi
G_{N}}({(r^{d-3}(-\psi^{2}+H_{0}^{2}r^{2}))^{\prime}\over
r^{d-2}}+\frac{2\dot{\psi}}{r}),\\\ p_{\parallel}=&{1\over 16\pi
G_{N}}\frac{(r^{d-3}(-\psi^{2}+H_{0}^{2}r^{2}))^{\prime\prime}+2(r^{d-3}\dot{\psi})^{\prime}}{r^{d-3}},\end{split}$
(2.10)

and $q^{\mu}=0$. We shall see that these quantities are indeed well-behaved in
later sections.

Before moving on detailed calculation, let us comment on a radial geodesic on
the Gullstrand-Painlevé coordinates. The particle radial geodesics is
determined by $ds^{2}=-d\tau^{2}$ with $\tau$ the particle proper time, which
reads

$-(1-\psi^{2})({dt\over d\tau})^{2}+2\psi{dt\over d\tau}{dr\over
d\tau}+({dr\over d\tau})^{2}=-1.$ (2.11)

On the other hand, the time component of the radial geodesic equation reads

$\displaystyle{d^{2}t\over d\tau^{2}}-\psi^{2}\psi^{\prime}({dt\over
d\tau})^{2}-2\psi\psi^{\prime}{dt\over d\tau}{dr\over
d\tau}-\psi^{\prime}({dr\over d\tau})^{2}=0,$ (2.12)

where we used
$\Gamma^{t}_{tt}=-\psi^{2}\psi^{\prime},~{}\Gamma^{t}_{tr}=-\psi\psi^{\prime},~{}\Gamma^{t}_{rr}=-\psi^{\prime}$.
Removing ${dr\over d\tau}$ from these two equations we find

$\displaystyle{d^{2}t\over d\tau^{2}}-\psi^{\prime}({dt\over
d\tau})^{2}=-\psi^{\prime}.$ (2.13)

This equation has an obvious solution of $t=\tau+\mathrm{const}$ for any
function $\psi$. This means that the time coordinate in the Gullstrand-
Painlevé metric can be understood as the proper time of a free particle. Then
this reduces the particle radial geodesic equation to

$\displaystyle{dr\over dt}=-\psi(t,r(t)).$ (2.14)

Plugging $\psi=-H(t)r$ into this gives $\dot{r}=H(t)r(t)$, which shows that
the radial coordinate of a particle geodesic $r(t)$ is proportional to the
scale factor $a(t)$. This is also expected from the transformation from the
FLRW metric to the Gullstrand-Painlevé one: $r=a(t)\check{r}$. On the other
hand, the radial null geodesic is determined by $ds^{2}=0$. This reads
${dr\over dt}=\mp 1-\psi$. Thus the infinite redshift surface is determined by
$\psi=\mp 1$. Note that the positive sign corresponds to an outgoing light-
ray, while the minus one an ingoing one in the coordinate system (2.4).

##  3 Collapsing shells of fluid with a perfect fluid inside

###  3.1 Scale factor for a perfect fluid with the cosmological constant

Before studying gravitational collapse or evaporation, we consider spacetime
filled with a perfect fluid and the dark energy or the cosmological constant
in general and determine the behavior of its scale factor. As argued in the
previous section, this situation is realized by setting $\psi=-H(t)r$ in the
Gullstrand-Painlevé metric (2.4). Indeed in this situation the radial
component of the pressure and the angular one determined as (2.10) become
identical:

$\displaystyle p_{\parallel}=p_{\perp}={d-2\over 16\pi
G_{N}}(-(d-1)(H(t)^{2}-H_{0}^{2})-2\dot{H}(t)).$ (3.1)

Let us consider a perfect fluid to satisfy an equation of state such that

$\displaystyle p=w\rho,$ (3.2)

where $p=p_{\parallel}=p_{\perp}$, $w$ is a constant greater than or equal to
$-1$. Using (2.10) we can rewrite this as

$\displaystyle\dot{H}(t)=-{1\over 2}(w+1)(d-1)(H(t)^{2}-H_{0}^{2}).$ (3.3)

Nontrivial solutions are either

$H(t)=H_{0}\tanh\left({H_{0}\over 2}(w+1)(d-1)(t-t_{0})\right)~{}~{}{\rm
or}~{}~{}H_{0}\coth\left({H_{0}\over 2}(w+1)(d-1)(t-t_{0})\right),$ (3.4)

where $t_{0}$ is an integration constant. The former solution goes to
vanishing in the zero cosmological constant limit $H_{0}\to 0$, which is not
the correct behavior of the scale factor. Therefore we choose the latter. Then
we can determine the scale factor $a(t)$ by solving $\dot{a}(t)=H(t)a(t)$ as

$\displaystyle a(t)=C_{1}\left(\sinh(\pm{H_{0}\over
2}(w+1)(d-1)(t-t_{0}))\right)^{2\over(w+1)(d-1)},$ (3.5)

where $C_{1}$ is another integration constant and the double sign is chosen so
that the term inside the bracket becomes positive. This solution correctly
reduces to the well-known scale factor for a perfect fluid with the equation
of state (3.2) in the zero cosmological constant limit:

$\displaystyle a(t)\to a_{0}|H_{0}(t-t_{0})|^{2\over(w+1)(d-1)},$ (3.6)

where we set $a_{0}:=C_{1}({(w+1)(d-1)\over 2})^{2\over(w+1)(d-1)}$.

Without loss of generality, we can set the integration constant $t_{0}$ to
zero by shifting the time coordinate, so that

$H(t)=H_{0}\coth\left({H_{0}\over 2}(w+1)(d-1)t\right).$ (3.7)

Then the sign of $H(t)$ is correlated to that of the time coordinate: for
$t<0$, $H(t)<0$, or the internal region is shrinking with acceleration, while
for $t>0$, $H(t)>0$, the internal region is expanding with acceleration. In
what follows, we use this convention.

Note that neither solutions of (3.4) go to the desired solution of (3.3) in
the limit $w\to-1$, in which case the correct solution would be such that
$H(t)$ is a general constant. Therefore we have to deal with the case of
$w=-1$ or $H(t)=H_{0}$ separately when we construct a solution describing
gravitational collapse or evaporation in this case. In this situation we
cannot determine the orbit of the shell only by the continuity condition. One
of the ways to determine the orbit of the shell is to specify the property of
the fluid by imposing an equation of state near the surface of the shell [36].
The analysis becomes more complicated. For completeness we perform it in
appendix A.

###  3.2 Shell structure and continuity condition

Let us construct a solution for gravitational collapse or evaporation of
shells of fluid in the presence of a perfect fluid inside and the cosmological
constant filled in the spacetime. In order to construct such a solution of our
interest, we first set the ansatz to realize such a configuration that a thick
shell of fluid exists on top of a perfect fluid. Thus at a fixed time slice
the space is generally separated into three regions denoted by I, II, III as
explained in the introduction: the region I is filled with a uniform perfect
fluid, the region II is with thick shell of fluid infalling, and the third one
is only with dark energy or cosmological constant. To describe the regions
precisely we prepare for a function $h(r)$ to describe the radial orbit of the
thick shell of fluid such that the orbits of the two edges of the thick shell
are described by $t=h(r)$ and $t=h(r)+\Delta$, where $\Delta$ is a positive
constant describing a width of the shells.111We flip the sign of the function
$h(r)$ in this paper compared to the one in [36].  More precisely, if $h(r)$
is a monotonically decreasing function, in which case gravitational collapse
happens, each region is described as follows:

$\displaystyle\left\\{\begin{array}[]{rcl}\mbox{region I}&=&\\{(t,r);0\leq
r\leq h^{-1}(t)\\}\\\ \mbox{region II}&=&\\{(t,r);h^{-1}(t)\leq r\leq
h^{-1}(t-\Delta)\\}\\\ \mbox{region III}&=&\\{(t,r);h^{-1}(t-\Delta)\leq
r\\}\\\ \end{array}\right.\quad\mbox{if $h^{\prime}(r)<0$.}$ (3.11)

On the other hand, if $h(r)$ monotonically increases, then gravitational
evaporation happens, and each region is described as

$\displaystyle\left\\{\begin{array}[]{rcl}\mbox{region I}&=&\\{(t,r);0\leq
r\leq h^{-1}(t-\Delta)\\}\\\ \mbox{region II}&=&\\{(t,r);h^{-1}(t-\Delta)\leq
r\leq h^{-1}(t)\\}\\\ \mbox{region III}&=&\\{(t,r);h^{-1}(t)\leq r\\}\\\
\end{array}\right.\quad\mbox{if $h^{\prime}(r)>0$.}$ (3.15)

A goal of the analysis is to determine the function $h(r)$ to satisfy the
Einstein equation with some physical condition such as the continuity at
interfaces and an equation of state of fluid. Note that the orbit of a
particle which consists of the shells satisfies a differential equation such
that

${dr\over dt}=\frac{1}{h^{\prime}(r(t))}.$ (3.16)

Let us take a strategy so as to consider solutions in the region I and III and
connect them to the region II smoothly. Without losing generality we can use
the coordinate system with the upper patch $U_{+}$ to cover the region I. As
discussed in the previous section, the Einstein equation can be solved in the
region I with the function $\psi_{+}$ chosen as

$\psi_{+}=-H(t)r,\quad\mbox{$(t,r)\in$ region I},$ (3.17)

where $H(t)$ is given by (3.7). Therefore if $t<0$ then $\psi_{+}\geq 0$,
while $\psi_{+}\leq 0$ if $t>0$, in this region. On the other hand, since the
metric describes the Schwarzschild black hole in the region III, as argued in
the previous section, $\psi_{+}$ may be written as

$\tilde{\psi}_{+}\overset{?}{=}-\sqrt{{m_{0}\over
r^{d-3}}+H_{0}^{2}r^{2}},\quad\mbox{$(t,r)\in$ region III}.$ (3.18)

Let us try to connect these solutions via the region II. A candidate of a
desired function may have a form such that

$\tilde{\psi}_{+}\overset{?}{=}-\sqrt{{m(t,r)\over
r^{d-3}}+H(t,r)^{2}r^{2}},\quad\mbox{$(t,r)\in$ region II},$ (3.19)

where $m(t,r),H(t,r)$ are smooth functions satisfying

$\begin{split}\lim_{(t,r)\to\rm region\,I}m(t,r)=0,~{}~{}\lim_{(t,r)\to\rm
region\,I}H(t,r)=H(t),\\\ \lim_{(t,r)\to\rm
region\,III}m(t,r)=m_{0},~{}~{}\lim_{(t,r)\to\rm
region\,III}H(t,r)=H_{0}.\end{split}$ (3.20)

However, it is soon noticed that this function cannot be smooth at the
interface between the region I and II, since $H(t)<0$ for $t<0$, so that
$\lim_{(t,r)\to\rm region\,I}\tilde{\psi}_{+}(t,r)=H(t)r\not=\psi_{+}$ for
$t<0$.

This shows that it is not possible to cover all the regions I, II and III by
one local patch $U_{+}$. Thus we use the local upper patch $U_{+}$ to cover
the region I and II and set the function $\psi_{+}$ in the region II as

$\psi_{+}(t_{+},r)=\sqrt{{m(t_{+},r)\over
r^{d-3}}+H(t_{+},r)^{2}r^{2}},\quad\mbox{$(t_{+},r)\in$ region I and II},$
(3.21)

where $m(t,r),H(t,r)$ satisfy (3.20), while we cover the region II and III by
the local patch $U_{-}$, where $\psi_{-}$ is chosen as

$\psi_{-}(-t_{-},r)=\sqrt{{m(-t_{-},r)\over
r^{d-3}}+H(-t_{-},r)^{2}r^{2}},\quad\mbox{$(-t_{-},r)\in$ region II and III}.$
(3.22)

Then $\psi_{+}$ with (3.21) smoothly connects to the one with (3.17) when
there exists $t_{+}<0$ in the region I, and the de-Sitter Schwarzschild black
hole solution is realized in the region III covered by the lower patch
$U_{-}$. These two patches can be glued consistently in the intersection by
the coordinate transformation $t_{+}=-t_{-}$ since $\psi_{+}=\psi_{-}$ in
$U_{+}\cap U_{-}$. In order to argue for the case with $t_{+}>0$ present in
the region I, we consider one global patch $U=U_{+}\cup U_{-}$ and cover all
the regions I, II and III by $U$, on which the metric is given by (2.4) with

$\psi(t,r)=\sqrt{{m(t,r)\over r^{d-3}}+H(t,r)^{2}r^{2}}.$ (3.23)

Then we do not need to care about the issue how to cover each region by local
patches. That is, it is sufficient to construct a smooth function $\psi$ for
all the regions to satisfy the Einstein equation, and have only to choose an
appropriate patch for the physical interpretation.

The function $m(t,r)$ describes a profile of thick shell of fluid and $H(t,r)$
a profile of a perfect fluid inside it. Thus we realize these functions by
preparing the corresponding slope functions $F(x)$, $\Theta(x)$, which are
monotonically increasing or decreasing functions valued between $0$ and $1$ in
the interval $[0,1]$ as

$\displaystyle m(t,r)=$ $\displaystyle m_{0}\theta(r)F({t-h(r)\over\Delta}),$
(3.24) $\displaystyle H(t,r)^{2}=$ $\displaystyle H_{0}^{2}+\delta
H(t)^{2}\Theta({t-h(r)\over\Delta}),$ (3.25)

with $\theta(r)$ the step function with $\theta(0)=0$,

$\delta
H(t)^{2}=H(t)^{2}-H_{0}^{2}=\left\\{\begin{array}[]{lr}\frac{H_{0}^{2}}{\sinh^{2}({H_{0}\over
2}(w+1)(d-1)t)}&w\not=-1\\\ 0&w=-1\end{array}\right..$ (3.26)

More specifically whether slope functions $F(x),\Theta(x)$ are upslope or
downslope depends on a situation such that

$\left\\{\begin{array}[]{lc}F^{\prime}(x)\geq 0,\Theta^{\prime}(x)\leq
0&\mbox{with $h^{\prime}(r)<0$ or gravitational collapse}\\\ F^{\prime}(x)\leq
0,\Theta^{\prime}(x)\geq 0&\mbox{with $h^{\prime}(r)>0$ or gravitational
evaporation}.\end{array}\right.$ (3.27)

The condition (3.20) is satisfied
$F(0)=\Theta(1)=\sigma,F(1)=\Theta(0)=1-\sigma$, where $\sigma$ is defined by

$\sigma:=\left\\{\begin{array}[]{lc}0,&\mbox{with $h^{\prime}(r)<0$ or
gravitational collapse}\\\ 1,&\mbox{with $h^{\prime}(r)>0$ or gravitational
evaporation}\end{array}\right.$ (3.28)

More formally, $\sigma=\theta(h^{\prime}(r))$.

These slope functions are inserted so as for some physical quantities to be
continuous at the interfaces. We request the density and the perpendicular or
radial component of the pressure to be continuous. We first consider the
continuity condition at interface between the region I and II. The continuity
of the density at this interface is mathematized as

$\lim_{x\to\sigma+0}\rho|_{t=h(r)+\Delta
x}=\lim_{x\to\sigma-0}\rho|_{t=h(r)+\Delta x},$ (3.29)

The density is given in each region as

$\rho=\left\\{\begin{array}[]{lr}{(d-2)\delta H(t)^{2}(d-1)\over 16\pi
G_{N}},&\mbox{region I}\\\
{(d-2)(m_{0}F^{\prime}({t-h(r)\over\Delta}){-h^{\prime}(r)\over\Delta}+\delta
H(t)^{2}\\{\Theta^{\prime}({t-h(r)\over\Delta}){-h^{\prime}(r)\over\Delta}r^{d-1}+\Theta({t-h(r)\over\Delta})(d-1)r^{d-2}\\})\over
16\pi G_{N}r^{d-2}},&\mbox{region II}\\\ 0,&\mbox{region III}\\\
\end{array}\right.$ (3.30)

Therefore the continuity condition reads

$\displaystyle\lim_{x\to\sigma}[m_{0}F^{\prime}(x)+\delta
H(t)^{2}\Theta^{\prime}(x)r^{d-1}]=0.$ (3.31)

For the case with $w=-1$, $\delta H(t)$ vanishes. Thus this condition is
equivalent to $F^{\prime}(\sigma)=0$. For the case with $w\not=-1$, the value

$\alpha:=-\lim_{x\to\sigma}\frac{\Theta^{\prime}(x)}{F^{\prime}(x)}$ (3.32)

is nonnegative since the slope functions $F,\Theta$ have the opposite sign of
the gradient. We assume $\alpha$ to be nonzero. Then employing (3.26) we
obtain

$\displaystyle\sinh^{2}({H_{0}\over
2}(w+1)(d-1)(h(r)+\sigma\Delta))={H_{0}^{2}\alpha r^{d-1}\over m_{0}}.$ (3.33)

This can be solved as
$h(r)+\sigma\Delta=\pm\frac{2}{{H_{0}}(w+1)(d-1)}\mathrm{arcsinh}(\sqrt{{\alpha\over
m_{0}}}H_{0}r^{d-1\over 2}),$ which can be rewritten as

$\displaystyle
h(r)=\left\\{\begin{array}[]{lr}-\frac{2}{{H_{0}}(w+1)(d-1)}\mathrm{arcsinh}(\sqrt{{\alpha\over
m_{0}}}H_{0}r^{d-1\over 2}),&h^{\prime}(r)<0\\\
-\Delta+\frac{2}{{H_{0}}(w+1)(d-1)}\mathrm{arcsinh}(\sqrt{{\alpha\over
m_{0}}}H_{0}r^{d-1\over 2}),&h^{\prime}(r)>0\\\ \end{array}\right..$ (3.36)

Comments are in order.222 A similar expression was obtained as (40) in [18]
from the junction condition. Just for comparison, we take the result in [20],
where a collapsing pressureless fluid (or dust) with density finely tuned as
$\rho=\mu/a^{3}$ was studied. The equation of the trajectory of the surface is
determined as the equation (40), which reads ${dt\over
dr}=\frac{-1}{(1-(H_{0}^{2}r^{2}+\frac{m_{0}}{r}))\sqrt{H_{0}^{2}r^{2}+\frac{m_{0}}{r}}}$
in our notation. On the other hand, in our study of a collapsing fluid with a
perfect fluid inside, the corresponding quantity is ${dh\over
dr}=\frac{-1}{(1+w)\sqrt{H_{0}^{2}r^{2}+{m_{0}\over\alpha}\frac{1}{r}}}$.
Firstly, the gravitational collapsing/evaporating solution is realized
accordingly when the FLRW universe shrinks/expands respectively in this model.
Secondly, in the zero cosmological constant limit, this solution reduces to

$\displaystyle
h(r)\to\left\\{\begin{array}[]{lr}-\frac{2}{(w+1)(d-1)}(\sqrt{{\alpha\over
m_{0}}}r^{d-1\over 2}),&h^{\prime}(r)<0\\\
-\Delta+\frac{2}{(w+1)(d-1)}(\sqrt{{\alpha\over m_{0}}}r^{d-1\over
2}),&h^{\prime}(r)>0\\\ \end{array}\right..$ (3.39)

The dependence of the radial coordinate matches the result in [31] at $d=4$.
We shall show soon that the emergent curvature singularity is always covered
by the emergent event horizon. On the other hand, on the continuity of the
radial component of the pressure, which reads $p_{\perp}=-\rho+{d-2\over 16\pi
G_{N}}{2\dot{\psi}(t,r)\over r}$, it is necessary and sufficient to consider
the continuity for $\dot{\psi}$ taking into account the continuity of $\rho$.
$\dot{\psi}$ is given by

$\dot{\psi}=\left\\{\begin{array}[]{lr}\frac{1}{2\psi}(\partial_{t}\delta
H(t)^{2}r^{2}),&\mbox{region I}\\\
\frac{1}{2\psi}\left({m_{0}F^{\prime}({t-h(r)\over\Delta}){1\over\Delta}\over
r^{d-3}}+\\{(\partial_{t}\delta H(t)^{2})\Theta({t-h(r)\over\Delta})+\delta
H(t)^{2}\Theta^{\prime}({t-h(r)\over\Delta})\frac{1}{\Delta}\\}r^{2}\right),&\mbox{region
II}\\\ 0,&\mbox{region III}\\\ \end{array}\right.$ (3.40)

Thus the continuity condition, $\lim_{x\to\sigma+0}\dot{\psi}|_{t=h(r)+\Delta
x}=\lim_{x\to\sigma-0}\dot{\psi}|_{t=h(r)+\Delta x}$, reduces to the same one
for the density $\rho$, (3.31). As a result the orbit of the shell of fluid is
determined and both gravitational collapse and evaporation are possible.

We next consider the continuity condition at the interface between II and III.
That of the density is now written as

$\lim_{x\to\sigma+0}\rho|_{t=h(r)+\Delta-\Delta
x}=\lim_{x\to\sigma-0}\rho|_{t=h(r)+\Delta-\Delta x}.$ (3.41)

From (3.30) the continuity condition for the density reads

$\displaystyle m_{0}F^{\prime}(1-\sigma){-h^{\prime}(r)\over\Delta}+\delta
H(t)^{2}\left(\Theta^{\prime}(1-\sigma){-h^{\prime}(r)\over\Delta}r^{d-1}+\Theta(1-\sigma)(d-1)r^{d-2}\right)=0.$
(3.42)

For the case with $w=-1$, $\delta H(t)=0$, which reduces this condition to
$F^{\prime}(1-\sigma)=0$. For the case with $w\not=-1$, from the continuity
condition between the region I and II, $h(r)$ is determined as (3.36).
Combining this, the above condition (3.42) is equivalent to

$\displaystyle F^{\prime}(1-\sigma)=\Theta^{\prime}(1-\sigma)=0.$ (3.43)

As in the same computation of the continuity condition between the region I
and II, we can show that the continuity of the perpendicular component of the
pressure reduces to the same continuity condition for the density (3.43).

We conclude that the function $h(r)$ is determined as (3.36) from the
continuity condition for $w\not=-1$, in which there exists a non-trivial
perfect fluid inside. Gravitational collapse/evaporation happens when the
internal FLRW universe shrinks/expands, respectively. On the other hand, the
function $h(r)$ cannot be determined only from the continuity condition in the
case with $w=-1$, where the inside of the shell of fluid is filled only with
the dark energy. This implies that one needs to impose a condition on thick
shell of fluid to determine its orbit in this situation, which is physically
reasonable since the property of fluid is not specified in this case. We
determine the function $h(r)$ in this case by imposing an equation of state
near the interface in the appendix.

For instruction we draw a schematic picture for the case of gravitational
collapse of this model in four dimensions, $d=4$, in Fig. 2. The location of
the event horizon of the final black hole is determined by $\psi=1$, where
$\psi$ takes the form in the region III:
$\psi=\sqrt{\frac{m_{0}}{r}+r^{2}H_{0}^{2}}$. This equation is solved by the
Cardano’s formula as $r=r_{k}$, where

$r_{k}=\omega^{k}\sqrt[3]{-\frac{m_{0}}{2H_{0}^{2}}+\sqrt{(\frac{m_{0}}{2H_{0}^{2}})^{2}-(\frac{1}{3H_{0}^{2}})^{3}}}+\omega^{3-k}\sqrt[3]{-\frac{m_{0}}{2H_{0}^{2}}-\sqrt{(\frac{m_{0}}{2H_{0}^{2}})^{2}-(\frac{1}{3H_{0}^{2}})^{3}}},~{}~{}~{}~{}~{}$
(3.44)

with $k=0,1,2$, and $\omega={-1+i\sqrt{3}\over 2}$ the cube root of unity. As
long as the $H_{0}$ is smaller than the critical value
$\frac{2}{3\sqrt{3}m_{0}}$, there exist 2 positive solutions out of 3 and the
smaller one is the radius of the location of the final black hole horizon,
$r_{\rm BH}=r_{1}$, while the bigger one is that of the location of the de-
Sitter horizon, $r_{\rm dS}=r_{2}$. On the other hand, the growing infinite
redshift surface is determined by $\psi=1$ in the region II, where $\psi$ is
given by
$\psi=\sqrt{\frac{m_{0}F(\frac{t-h(r)}{\Delta})}{r}+r^{2}H_{0}^{2}(1+\frac{\Theta(\frac{t-h(r)}{\Delta})}{\sinh^{2}(\frac{H_{0}}{2}(w+1)3t)})}$.
The infinite redshift surface extends up to the location of the final black
hole horizon, which intersect at the time $t_{BH}=\Delta+h(r_{BH})$. For later
convenience we denote the intersecting point $(t_{BH},r_{BH})$ by $P$. Note
that the explicit expression of the infinite redshift surface depends on a
choice of the slope functions $F,\Theta$.

(a) Blackhole formation (b) Whitehole evaporation

Figure 2: (a) depicts the blackhole formation from infalling shells of fluid.
The bottom one colored by pink is the region I, where a perfect fluid and the
dark energy coexist, the middle one colored by orange is the region II, where
the shells of fluid are infalling, and the upper one is the region III, where
there exists only the dark energy. The orbit of the thick shells of fluid is
drawn by choosing the parameters $w=0.1,\Delta=10,m_{0}=4,\Lambda=0.025$,
$F(x)=\sin^{2}{\pi x\over 2},\Theta(x)=\cos^{2}{\pi x\over 2}$. The innermost
shell reaches the origin at $t=0$, at which the curvature singularity starts
to emerge. The growing infinite redshift surface is given by $t=h(r)+\Delta
F^{-1}(\frac{r}{m}(1-\frac{\Lambda}{3}r^{2}))$. Other curves such as the last
ray-out, the location of the cosmological horizon, which are described by some
differential equations explained in the main text, are schematically drawn and
not precise. (b) depicts the whitehole evaporation into blowing shells of
fluid. We chose the same values for the parameters.

Now let us show that the emergent curvature singularity is covered by the
emergent event horizon. This is sufficient to show that there exists a light-
ray which is emitted outward from the origin before the singularity emerges
but does not go far away to infinity. A light-ray emitted outward from the
origin is described by the radial null geodesic ${dt\over dr}=1-\psi$. For a
while after the emittance, the light-ray moves in the perfect fluid or the
region I. Therefore during this period $\psi=-H(t)r$, where $H(t)$ is given by
(3.7). Note that during this period, $t<0$, and thus $H(t)<0$. Then after some
time the right-ray reaches the interface between the perfect fluid and the
thick shell, and starts to go through the thick shell or the region II. During
this period, $\psi=\sqrt{{m(t,r)\over r^{d-3}}+H(t,r)^{2}r^{2}}$. This holds
for all the right-rays before they get out of the region II. When the right-
lay comes out the region II, whether it goes outward or inward depends on when
it is emitted. The light-ray which passes through the point $P$ does not come
out, because it is trapped at the event horizon of the final static black hole
geometry. This is the last light-ray and its orbit yields that of the event
horizon of the growing black hole. Any light-ray emitted before the last
light-ray goes outward to infinity, while any one emitted after the last
light-ray goes inward and is trapped inside. Thus the emergent curvature
singularity in this model is always covered by the event horizon.

##  4 Energy distribution and conservation

We compute the energy distribution of the system using the definition given in
[44, 45, 46]. We denote the energy inside the radius $R$ at time $t$ by
$E(t,R)$:

$\displaystyle E(t,R)=$ $\displaystyle\int_{\sqrt{x^{2}}\leq
R}d^{d-1}x\sqrt{|g|}T^{t}{}_{\mu}n^{\mu},$ (4.1)

where $n^{\mu}=-\delta^{\mu}_{t}$. Using (2.7) we can compute this as follows.

$\displaystyle E(t,R)=$ $\displaystyle{d-2\over 16\pi
G_{N}}V_{d-2}\int_{0}^{R}dr{\partial_{r}(m(t,r)+(H(t,r)^{2}-H_{0}^{2})r^{d-1})}$
$\displaystyle=:$ $\displaystyle E_{1}(t,R)+E_{2}(t,R)+E_{3}(t,R)$ (4.2)

where $V_{d-2}$ is the volume of the internal Einstein manifold,

$\displaystyle E_{1}(t,R)=$ $\displaystyle M\int_{0}^{R}dr({d\theta(r)\over
dr}\widetilde{F}({t-h(r)\over\Delta}))=M\widetilde{F}({t-h(0)\over\Delta}),$
(4.3) $\displaystyle E_{2}(t,R)=$ $\displaystyle
M\int_{0}^{R}dr(\partial_{r}\widetilde{F}({t-h(r)\over\Delta}))=M(\widetilde{F}({t-h(R)\over\Delta})-\widetilde{F}(\frac{t-h(0)}{\Delta})),$
(4.4) $\displaystyle E_{3}(t,R)=$ $\displaystyle
M\frac{(H(t)^{2}-H_{0}^{2})}{m_{0}}\widetilde{\Theta}({t-h(R)\over\Delta})R^{d-1},$
(4.5)

where $M={d-2\over 16\pi G_{N}}V_{d-2}m_{0}$ and
$\widetilde{F},\widetilde{\Theta}$ are extended functions of $F,\Theta$ to the
region outside the interval $[0,1]$ with the values unchanged:

$\widetilde{F}(x)=\left\\{\begin{array}[]{lc}F(0),&x\leq 0\\\ F(x),&0\leq
x\leq 1\\\ F(1),&1\leq x\\\
\end{array}\right.,\quad\widetilde{\Theta}(x)=\left\\{\begin{array}[]{lc}\Theta(0),&x\leq
0\\\ \Theta(x),&0\leq x\leq 1\\\ \Theta(1),&1\leq x\\\ \end{array}\right.$
(4.6)

$E_{1}(t,R)$ is the energy localized at the origin or the mass of the black
hole at a fixed time $t$. $E_{2}(t,R)+E_{3}(t,R)$ is the energy of both the
perfect fluid and the shell of fluid inside the radius $R$. Note that the
energy of the perfect fluid inside the shell is converted into that of the
shells of fluid in the region II.

We comment on the conservation of the total energy, which is defined by
$E(t)=\lim_{R\to\infty}E(t,R)$.333The momentum for the radial direction is
defined by
$P^{r}=\int_{\mathbf{R}^{d-1}}d^{d-1}x\sqrt{|g|}T^{t}{}_{\mu}v^{\mu}$, where
$v^{\mu}=\delta^{\mu}_{r}$. Since this system does not radiate anything, the
total energy is expected to be conserved. Whether the total energy is
conserved or not can be tested by the conservation test equation [46]. The
statement is that the charge defined by
$Q=\int_{\mathbf{R}^{d-1}}d^{d-1}x\sqrt{|g|}T^{t}{}_{\mu}v^{\mu},$ is
conserved if the following equation is satisfied:

$T^{\mu}{}_{\nu}\nabla_{\mu}v^{\nu}=0.$ (4.7)

To discuss the conservation of the energy, we choose $v^{\mu}=n^{\mu}$. Since
$\Gamma^{i}_{\mu t}=0$ for all $\mu$, this equation reduces to

$T^{a}{}_{b}\nabla_{a}n^{b}=0$ (4.8)

where $a,b=t,r$. One can check that this conservation test equation holds by
using

$\nabla_{a}n^{b}=\left(\begin{array}[]{cc}\psi^{2}\psi^{\prime}&\psi\psi^{\prime}\\\
-\psi^{\prime}\psi^{3}+\psi^{\prime}\psi-\dot{\psi}&-\psi^{2}\psi^{\prime}\\\
\end{array}\right).$ (4.9)

It is also not difficult to compute the total energy directly by using
$\widetilde{F}({t-h(\infty)\over\Delta})=1,\widetilde{\Theta}({t-h(\infty)\over\Delta})=0$.
The result is $E(t)=M$, where $M$ is the final or initial black hole mass in
the case of gravitational collapse or evaporation, respectively. The total
energy of the system is conserved as expected, and matches the ADM or quasi-
local mass computed in the asymptotic behavior of the metric.

##  5 Discussion

We have investigated gravitational collapse or evaporation of shells of fluid
in the isotropic homogeneous universe without any radiation described by the
Einstein gravity with cosmological constant. To this end we have chosen the
ansatz of a metric written by the Gullstrand-Painlevé coordinates which suits
to describe the isotropic homogeneous space inside the shells of fluid and the
de-Sitter black hole outside them. After determining the scale factor for a
perfect fluid to obey an equation of state together with cosmological
constant, we have determined the orbit of thick shell of fluid from the
continuity conditions for the density and the perpendicular component of the
pressure at interface. When there is no perfect fluid inside the shells, the
orbit of the shells cannot be determined by the continuity condition. We have
determined the orbit of the shells by imposing an equation of state near its
surface as done in [36] in appendix A. We have confirmed that the total energy
of this model is independent of the given time evolution.

In this paper we dealt with the matter constituents of shells as fluid. This
approximation can be justified only when the gravitational coupling constant
is small enough to neglect the back-reaction. If we specified the matter field
and took into account back-reaction, then the analysis would be more
complicated. It would be interesting to study the system specifying the matter
and taking into account the back-reaction perturbatively and what happens when
the coupling constant is not small. In such a situation the shell crossing
singularity would become strong obstruction to constructing a solution as
pointed out in [36].

Another way to make the system closer to a realistic situation is to consider
the time dependence of the radial component of the Gullstrand-Painlevé metric
as considered in [5]. Such a metric can describe not only the radiation of
gravitational wave but also the isotropic homogeneous universe with different
topology. It would be interesting to evaluate the energy of the gravitational
wave by computing the differences of energy between the shells and the final
black hole. It might be amusing to study the possibility to construct a
solution of the Einstein gravity with a change of topology of the universe
employing such a generalized Gullstrand-Painlevé metric.

We hope to come back to these issues in near future.

## Acknowledgement

The author would like to thank Shigeki Sugimoto for valuable comments on the
draft. This work is supported in part by the Grant-in-Aid of the Japanese
Ministry of Education, Sciences and Technology, Sports and Culture (MEXT) for
Scientific Research (No. JP19K03847, 19H01897).

##  Appendix A White hole evaporation during inflation

In this appendix, as asserted in section 3.2, we study the case with $w=-1$ or
$H(t)=H_{0}$, where the internal of the shell is filled only with dark energy,
so that the function $\psi$ in the Gullstrand-Painlevé metric is given by

$\psi=-\sqrt{{m(t,r)\over r^{d-3}}+H_{0}^{2}r^{2}},$ (A.1)

where $m(t,r)$ takes the form of (3.24). In this case only the continuity
condition does not determine the function of the orbit of the shells $h(r)$.
We determine this function by imposing an equation of state near the surface.
An equation of state studied in [36] is the usual one such that pressure is
proportional to density unless the pressure is mixed with its components
$p_{\parallel}$ and $p_{\perp}$:

$P=w_{0}\rho,$ (A.2)

where $w_{0}\geq-1$ and $P$ is some mixture of the components of the pressure
so that $P=p_{\parallel}+v_{0}p_{\perp}$ with $v_{0}$ a constant. Recall that
the surface or the outermost shell is specified by $t=h(r)+\Delta(1-\sigma)$,
where $\sigma$ is given by (3.28) and takes $0$ or $1$ for the case of
gravitational collapse or evaporation, respectively. We impose the equation of
state (A.2) at a layer specified by $t=h(r)+\Delta(1-\sigma-\varepsilon)$,
where $\varepsilon$ is a small positive number.

###  A.1 Transverse pressure

We first study the case where the pressure in (A.2) consists only of the
parallel component or $v_{0}=0$. In this situation
$p_{\parallel}=w^{(0)}_{\parallel}\rho$ with $w^{(0)}_{\parallel}\geq-1$.
Using (2.10) one can rewrite this as

$(-m^{\prime}(t,r)+{\dot{m}(t,r)\over\psi})^{\prime}=w_{\parallel}\frac{1}{r}m^{\prime}(t,r),$
(A.3)

where $w_{\parallel}=(d-2)w^{(0)}_{\parallel}$. This can be further computed
by employing (3.24) as

$\displaystyle-
yh^{\prime}(h^{\prime}+\frac{1}{\bar{\psi}})+(h^{\prime\prime}-\frac{{\bar{\psi}}^{\prime}}{{\bar{\psi}}^{2}})=-\frac{w_{\parallel}}{r}h^{\prime},$
(A.4)

where $y={F^{\prime\prime}(1-\sigma-\varepsilon)\over
F^{\prime}(1-\sigma-\varepsilon)\Delta}$, ${\bar{\psi}}=-\sqrt{{\bar{m}\over
r^{d-3}}+H_{0}^{2}r^{2}}$ with $\bar{m}=m_{0}F(1-\sigma-\varepsilon)$ .

We cannot solve the differential equation (A.4) analytically with general
parameterization, so we solve it simplifying the situation with the help of
mathematica. We first consider the case with $y=0$. In this situation (A.4)
reduces to

$\displaystyle
h^{\prime\prime}-\frac{{\bar{\psi}}^{\prime}}{{\bar{\psi}}^{2}}=-\frac{w_{\parallel}}{r}h^{\prime},$
(A.5)

This can be solved in terms of $h^{\prime}(r)$ as

$\displaystyle h^{\prime}(r)|_{y=0}=$ $\displaystyle
c_{1}r^{-w_{\parallel}}-r^{\frac{d-3}{2}}\left(\frac{2w_{\parallel}\,{}_{2}F_{1}\left(\frac{1}{2},\frac{d+2w_{\parallel}-3}{2(d-1)};\frac{-3d-2w_{\parallel}+5}{2-2d};-\frac{H_{0}^{2}r^{d-1}}{\bar{m}}\right)}{(d+2w_{\parallel}-3)\sqrt{\bar{m}}}-\frac{1}{\sqrt{H_{0}^{2}r^{d-1}+\bar{m}}}\right),$

where $c_{1}$ is an integration constant and ${}_{2}F_{1}$ is the
hypergeometric function. This can be further integrated as

$\displaystyle h(r)|_{y=0}=$
$\displaystyle\frac{c_{1}r^{1-w_{\parallel}}}{1-w_{\parallel}}+\frac{2\mathrm{arcsinh}\left(\frac{H_{0}}{\sqrt{\bar{m}}}r^{\frac{d-1}{2}}\right)}{H_{0}(d-1)}$
$\displaystyle-\frac{2w_{\parallel}r^{\frac{d-1}{2}}\Gamma\left(\frac{d+2w_{\parallel}-3}{2(d-1)}\right)\,_{3}{F}_{2}\left(\frac{1}{2},\frac{d+2w_{\parallel}-3}{2(d-1)},\frac{1}{2};\frac{-3d-2w_{\parallel}+5}{2-2d},\frac{3}{2};-\frac{H_{0}^{2}r^{d-1}}{\bar{m}}\right)}{(d-1)^{2}\sqrt{\bar{m}}\Gamma\left(\frac{d+2w_{\parallel}-5}{2(d-1)}\right)},$
(A.6)

where ${}_{p}F_{q}$ is the generalized hypergeometric function. The last term
appears due to the effect of pressure which a particle consisting of fluid
receives. We fix the integration constant $c_{1}$ by requesting that the
particle radial geodesic reduces to that of a freely falling one in the
pressureless limit $w_{\parallel}\to 0$. In this limit the solution reduces to

$\displaystyle h(r)|_{y=0}\to$ $\displaystyle
c_{1}r+\frac{2\mathrm{arcsinh}\left(\frac{H_{0}}{\sqrt{\bar{m}}}r^{\frac{d-1}{2}}\right)}{(d-1)H_{0}},$
(A.7)

which satisfies $h^{\prime}|_{y=w_{\parallel}=0}=c_{1}-\frac{1}{\bar{\psi}}$.
Then the particle orbit (3.16) becomes ${dr\over dt}={1\over
c_{1}-\frac{1}{\bar{\psi}}}$. This equation matches the geodesic of a free
particle given in (2.14) if and only if $c_{1}$ vanishes. As a result we
obtain the solution in the case with $y=0$ as

$\displaystyle h(r)|_{y=0}=$
$\displaystyle\frac{2\mathrm{arcsinh}\left(\frac{H_{0}}{\sqrt{\bar{m}}}r^{\frac{d-1}{2}}\right)}{H_{0}(d-1)}-\frac{2w_{\parallel}r^{\frac{d-1}{2}}\Gamma\left(\frac{d+2w_{\parallel}-3}{2(d-1)}\right)\,_{3}{F}_{2}\left(\frac{1}{2},\frac{d+2w_{\parallel}-3}{2(d-1)},\frac{1}{2};\frac{-3d-2w_{\parallel}+5}{2-2d},\frac{3}{2};-\frac{H_{0}^{2}r^{d-1}}{\bar{m}}\right)}{(d-1)^{2}\sqrt{\bar{m}}\Gamma\left(\frac{d+2w_{\parallel}-5}{2(d-1)}\right)}.$
(A.8)

Notice that the function $h(r)$ is monotonically increasing, which is also
confirmed from $h^{\prime}=-\frac{1}{\psi}=({m(t,r)\over
r^{d-3}}+H_{0}^{2}r^{2})^{-{1\over 2}}$. Note that the gravitational
evaporation can happen but the gravitational collapse does not happen when the
space expands exponentially in this model as in the case with a perfect fluid
inside studied in the previous section.

Let us comment on a solution of (A.4) when the cosmological constant vanishes
and the dimensions are four. In this situation (A.4) can be solved in terms of
$h^{\prime}(r)$ by using some hypergeometric functions. The integration
constant can be determined by requesting the solution to reduce in the limit
$y\to 0$ to the expression of (A.8) with $H_{0}\to 0$. The result is

$\displaystyle
h^{\prime}(r)|_{H_{0}=0}=\frac{\sqrt{r}\,_{1}{F}_{1}\left(\frac{4}{3};\frac{2(w_{\parallel}+2)}{3};-\frac{2r^{3/2}y}{3\sqrt{\bar{m}}}\right)}{\sqrt{\bar{m}}(2w_{\parallel}+1)\,_{1}{F}_{1}\left(\frac{1}{3};\frac{2w_{\parallel}+1}{3};-\frac{2r^{3/2}y}{3\sqrt{\bar{m}}}\right)}.$
(A.9)

Note that we could not integrate this further.

###  A.2 Longitudinal component

We first consider the case that the pressure in the equation of state consists
only of the longitudinal one, $p_{\perp}=w_{\perp}\rho$, which is imposed at a
layer specified by $t=h(r)+\Delta(1-\sigma-\varepsilon)$. Employing (2.10) we
can rewrite this as ${\dot{m}(t,r)/\psi}=(w_{\perp}+1)m(t,r)^{\prime}.$ We can
compute this further by employing (3.24) as
$h^{\prime}=-1/(w_{\perp}+1){\bar{\psi}},$ where ${\bar{\psi}}$ is given below
(A.4). This can be solved as

$h(r)=c_{1}+\frac{2\mathrm{arcsinh}\left(\frac{H_{0}}{\sqrt{\bar{m}}}r^{\frac{d-1}{2}}\right)}{(d-1)H_{0}(w_{\perp}+1)}.$
(A.10)

We fix an integration constant $c_{1}$ by requesting this solution to reduce
to the one for the transverse pressure (A.8) in the pressureless limit
$w_{\perp},w_{\parallel}\to 0$. The result is $c_{1}=0$.

###  A.3 Mixed pressure

Finally we consider the general case (A.2). Using (2.10) and (3.24) we can
rewrite this as follows.

$\displaystyle-
yh^{\prime}(h^{\prime}+\frac{1}{\bar{\psi}})+(h^{\prime\prime}-\frac{\bar{\psi}^{\prime}}{\bar{\psi}^{2}})+\frac{v}{r\bar{\psi}}=-\frac{w}{r}h^{\prime},$
(A.11)

where $v=(d-2)v_{0}$, $w=(d-2)(w_{0}+v_{0})$, $y$ is defined below (A.4). As
in subsection A.1 we solve this in the case with $y=0$. In this situation the
above differential equation reduces to

$\displaystyle(h^{\prime\prime}-\frac{\bar{\psi}^{\prime}}{\bar{\psi}^{2}})+\frac{v}{r\bar{\psi}}=-\frac{w}{r}h^{\prime},$
(A.12)

We first solve this in terms of $h^{\prime}(r)$:

$\displaystyle h^{\prime}(r)|_{y=0}=$ $\displaystyle
c_{1}r^{-w}-\frac{r^{\frac{d-3}{2}}\left(2(w-v)\sqrt{\frac{H_{0}^{2}r^{d-1}}{\bar{m}}+1}\,_{2}F_{1}\left(\frac{1}{2},\frac{d+2w-3}{2(d-1)};\frac{-3d-2w+5}{2-2d};-\frac{H_{0}^{2}r^{d-1}}{\bar{m}}\right)-d-2w+3\right)}{(d+2w-3)\sqrt{H_{0}^{2}r^{d-1}+\bar{m}}},$
(A.13)

where $c_{1}$ is an integration constant. Requesting this solution to reduce
to the previous one in the limit $v_{0}\to 0$ we can fix it as $c_{1}=0$. We
can further integrate it as follows.

$\displaystyle h(r)|_{y=0}=$
$\displaystyle\frac{2\mathrm{arcsinh}\left(\frac{H_{0}}{\sqrt{\bar{m}}}r^{\frac{d-1}{2}}\right)}{(d-1)H_{0}}-\frac{2\sqrt{\pi}r^{\frac{d-1}{2}}(w-v)\Gamma\left(\frac{-3d-2w+5}{2-2d}\right)\,_{3}\tilde{F}_{2}\left(\frac{1}{2},\frac{d+2w-3}{2(d-1)},\frac{1}{2};\frac{-3d-2w+5}{2-2d},\frac{3}{2};-\frac{H_{0}^{2}r^{d-1}}{\bar{m}}\right)}{(d-1)\sqrt{\bar{m}}(d+2w-3)}.$
(A.14)

When the cosmological constant vanishes and the dimensions are four, (A.11)
can be solved in terms of $h^{\prime}(r)$ by using some hypergeometric
functions. We fix an integration constant by requesting the solution to reduce
in the limit $y\to 0$ to (A.14) with $H_{0}\to 0$. The result is

$\displaystyle
h^{\prime}(r)|_{H_{0}=0}=\frac{\sqrt{r}(2v+1)\,_{1}{F}_{1}\left(\frac{2(v+2)}{3};\frac{2(w+2)}{3};-\frac{2r^{3/2}y}{\sqrt{\bar{m}}}\right)}{3\sqrt{\bar{m}}(2w+1)\,_{1}{F}_{1}\left(\frac{2v+1}{3};\frac{2w+1}{3};-\frac{2r^{3/2}y}{3\sqrt{\bar{m}}}\right)}.$
(A.15)

## References

  * [1] R. Genzel, F. Eisenhauer, and S. Gillessen, “The Galactic Center massive black hole and nuclear star cluster,” Reviews of Modern Physics 82 (Oct., 2010) 3121–3195, 1006.0064.
  * [2] Event Horizon Telescope Collaboration, K. Akiyama et al., “First M87 Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole,” Astrophys. J. Lett. 875 (2019) L1, 1906.11238.
  * [3] J. R. Oppenheimer and H. Snyder, “On Continued Gravitational Contraction,” Physical Review 56 (Sept., 1939) 455–459.
  * [4] B. Datt, “Über eine Klasse von Lösungen der Gravitationsgleichungen der Relativitüt,” Zeitschrift fur Physik 108 (1938) 314–321.
  * [5] C. W. Misner and D. H. Sharp, “Relativistic equations for adiabatic, spherically symmetric gravitational collapse,” Phys. Rev. 136 (1964) B571–B576.
  * [6] R. Penrose, “Gravitational collapse and space-time singularities,” Phys. Rev. Lett. 14 (1965) 57–59.
  * [7] B. K. Harrison, K. S. Thorne, M. Wakano, and J. A. Wheeler, Gravitation Theory and Gravitational Collapse.  1965\.
  * [8] P. C. Vaidya, “An Analytical Solution for Gravitational Collapse with Radiation,” Astrophys. J. 144 (June, 1966) 943.
  * [9] G. C. McVittie, “An Example of Gravitational Collapse in General Relativity,” Astrophys. J. 143 (Mar., 1966) 682.
  * [10] H. Nariai, “A Simple Model for Gravitational Collapse with Pressure Gradient,” Progress of Theoretical Physics 38 (07, 1967) 92–106, https://academic.oup.com/ptp/article-pdf/38/1/92/5366896/38-1-92.pdf.
  * [11] P. S. Joshi and D. Malafarina, “Recent developments in gravitational collapse and spacetime singularities,” Int. J. Mod. Phys. D 20 (2011) 2641–2729, 1201.3660.
  * [12] Supernova Search Team Collaboration, A. G. Riess et al., “Observational evidence from supernovae for an accelerating universe and a cosmological constant,” Astron. J. 116 (1998) 1009–1038, astro-ph/9805201.
  * [13] SNLS Collaboration, P. Astier et al., “The Supernova Legacy Survey: Measurement of $\Omega_{M}$, $\Omega_{\Lambda}$ and ${\cal w}$ from the first year data set,” Astron. Astrophys. 447 (2006) 31–48, astro-ph/0510447.
  * [14] A. A. Penzias and R. W. Wilson, “A Measurement of excess antenna temperature at 4080-Mc/s,” Astrophys. J. 142 (1965) 419–421.
  * [15] Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, 1807.06209.
  * [16] S. Weinberg, “The Cosmological Constant Problem,” Rev. Mod. Phys. 61 (1989) 1–23. [,569(1988)].
  * [17] D. Garfinkle and C. Vuille, “Gravitational collapse with a cosmological constant,” General Relativity and Gravitation 23 (Apr., 1991) 471–475.
  * [18] M. Cissoko, J. C. Fabris, J. Gariel, G. Le Denmat, and N. O. Santos, “Gravitational dust collapse with cosmological constant,” gr-qc/9809057.
  * [19] Y. Yamanaka, K.-i. Nakao, and H. Sato, “Motion of a dust shell in the space time with a cosmological constant,” Prog. Theor. Phys. 88 (1992) 1097–1106.
  * [20] D. Markovic and S. Shapiro, “Gravitational collapse with a cosmological constant,” Phys. Rev. D 61 (2000) 084029, gr-qc/9912066.
  * [21] K. Lake, “Gravitational collapse of dust with a cosmological constant,” Phys. Rev. D 62 (2000) 027301, gr-qc/0002044.
  * [22] T. A. Madhav, R. Goswami, and P. S. Joshi, “Gravitational collapse in asymptotically anti-de Sitter/de Sitter backgrounds,” Phys. Rev. D 72 (2005) 084029, gr-qc/0502081.
  * [23] S. G. Ghosh and D. W. Deshkar, “Higher dimensional dust collapse with a cosmological constant,” Astrophys. Space Sci. 310 (2007) 111–117, gr-qc/0607142.
  * [24] M. Sharif and Z. Ahmad, “Gravitational Perfect Fluid Collapse With Cosmological Constant,” Mod. Phys. Lett. A 22 (2007) 1493–1502, gr-qc/0610033.
  * [25] G. Darmois, Les ’equations de la gravitation einsteinienne.  No. 25 in M’emorial des sciences math’ematiques. Gauthier-Villars, 1927\.
  * [26] S. O’Brien and J. L. Synge, “Jump conditions at discontinuities in general relativity.” Commun. Dublin Inst. Advanced Stud., Ser. A 9, 20 p. (1952)., 1952.
  * [27] A. Lichnerowicz, Theories relativistes de la gravitation et de l’electromagnetisme. Relativite generale et theories unitaires.  1955\.
  * [28] W. Israel, “Discontinuities in spherically symmetric gravitational fields and shells of radiation,” Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 248 (1958), no. 1254, 404–414.
  * [29] C. W. Misner, K. S. Thorne, and J. A. Wheeler, Gravitation.  W. H. Freeman, San Francisco, 1973.
  * [30] L. D. Landau and E. M. Lifschits, The Classical Theory of Fields, vol. Volume 2 of Course of Theoretical Physics.  Pergamon Press, Oxford, 1975.
  * [31] R. J. Adler, J. D. Bjorken, P. Chen, and J. S. Liu, “Simple analytic models of gravitational collapse,” Am. J. Phys. 73 (2005) 1148–1159, gr-qc/0502040.
  * [32] A. Gullstrand, Allgemeine Lösung des statischen Einkörperproblems in der Einsteinschen Gravitationstheorie.  Arkiv för matematik, astronomi och fysik. Almqvist & Wiksell, 1922\.
  * [33] P. Painlevé, “La mécanique classique et la théorie de la relativité,” Comptes Rendus Academie des Sciences (serie non specifiee) 173 (Jan., 1921) 677–680.
  * [34] A. J. S. Hamilton and J. P. Lisle, “The River model of black holes,” Am. J. Phys. 76 (2008) 519–532, gr-qc/0411060.
  * [35] Y. Kanai, M. Siino, and A. Hosoya, “Gravitational collapse in Painleve-Gullstrand coordinates,” Prog. Theor. Phys. 125 (2011) 1053–1065, 1008.0470.
  * [36] S. Yokoyama, “An analytic model for gravitational collapse of spherical matter under mixed pressure,” 2105.09676.
  * [37] F. Fayos, X. Jaén, E. Llanta, and J. M. M. Senovilla, “Matching of the Vaidya and Robertson-Walker metric,” Classical and Quantum Gravity 8 (1991) 2057–2068.
  * [38] F. Fayos, X. Jaén, E. Llanta, and J. M. M. Senovilla, “Interiors of vaidya’s radiating metric: Gravitational collapse,” Phys. Rev. D 45 (Apr, 1992) 2732–2738.
  * [39] D. M. Eardley and L. Smarr, “Time function in numerical relativity. Marginally bound dust collapse,” Phys. Rev. D 19 (1979) 2239–2259.
  * [40] D. Christodoulou, “Violation of cosmic censorship in the gravitational collapse of a dust cloud,” Commun. Math. Phys. 93 (1984) 171–195.
  * [41] K. Lake, “Naked singularities in gravitational collapse which is not self-similar,” Phys. Rev. D 43 (Feb, 1991) 1416–1417.
  * [42] S. L. Shapiro and S. A. Teukolsky, “Formation of naked singularities: The violation of cosmic censorship,” Phys. Rev. Lett. 66 (1991) 994–997.
  * [43] D. Christodoulou, “The Formation of Black Holes in General Relativity,” in 12th Marcel Grossmann Meeting on General Relativity, pp. 24–34.  5, 2008.  0805.3880.
  * [44] V. Fock, Theory of Space, Time, and Gravitation.  Pergamon Press, New York, USA, 1959.
  * [45] S. Aoki, T. Onogi, and S. Yokoyama, “Conserved charges in general relativity,” Int. J. Mod. Phys. A 36 (2021), no. 10, 2150098, 2005.13233.
  * [46] S. Aoki, T. Onogi, and S. Yokoyama, “Charge conservation, entropy current and gravitation,” Int. J. Mod. Phys. A 36 (2021), no. 29, 2150201, 2010.07660.
  * [47] F. Kottler, “Über die physikalischen Grundlagen der Einsteinschen Gravitationstheorie,” Annalen Phys 56 (1918) 401–462.
  * [48] H. Weyl, “Über die statischen, kugelsymmetrischen Lösungen von Einsteins ”kosmologischen” Gravitationsgleichungen,” Phys.Z. 20 (1919) 31–34.


# Fast Diffusion Probabilistic Model Sampling through the lens of Backward
Error Analysis

Yansong Gao1,∗, Zhihong Pan2, Xin Zhou2, Le Kang2, Pratik Chaudhari1
1University of Pennsylvania, ∗Work done during an internship, 2Baidu Research
USA

###### Abstract

Denoising diffusion probabilistic models (DDPMs) are a class of powerful
generative models. The past few years have witnessed the great success of
DDPMs in generating high-fidelity samples. A significant limitation of the
DDPMs is the slow sampling procedure. DDPMs generally need hundreds or
thousands of sequential function evaluations (steps) of neural networks to
generate a sample. This paper aims to develop a fast sampling method for DDPMs
requiring much fewer steps while retaining high sample quality. The inference
process of DDPMs approximates solving the corresponding diffusion ordinary
differential equations (diffusion ODEs) in the continuous limit. This work
analyzes how the backward error affects the diffusion ODEs and the sample
quality in DDPMs. We propose fast sampling through the Restricting Backward
Error schedule (RBE schedule) based on dynamically moderating the long-time
backward error. Our method accelerates DDPMs without any further training. Our
experiments show that sampling with an RBE schedule generates high-quality
samples within only 8 to 20 function evaluations on various benchmark
datasets. We achieved 12.01 FID in 8 function evaluations on the ImageNet
$128\times 128$, and a $20\times$ speedup compared with previous baseline
samplers.

##  1 Introduction

Denoising Diffusion Probabilistic Models (DDPMs) [29, 10] have emerged as a
powerful class of generative models. DDPMs model the data distribution through
an iterative denoising process. Recently, DDPMs have achieved impressive
performances on a variety of applications, including unconditional and
conditional image generation [31, 23, 5, 11], video generation [12], 3D point
cloud generation [22], text to speech [3, 15], and image super-resolution [25,
18]. Despite the remarkable success, inference of DDPMs ( e.g., sampling and
density estimation) takes much longer to produce high-quality results compared
with other generative models such as GANs [7]. To give a few examples, DDPMs
may take $4000$ steps [23] to generate one sample. Each step requires
evaluating the learned neural network once. Such inefficiency becomes a
critical bottleneck for the applications of DDPMs. For this reason, there is
an urgent request to design fast samples for DDPMs.

Several works aim at speeding up the inference process for DDPMs. One category
of methods includes knowledge distillation [26, 21], learning sampling
trajectory, and searching for better noise schedule [27, 17, 23, 34]. Such
methods require a training stage before the sampling stage. Hence their
applicability and effectivity might be limited when adapting their methods to
different models, datasets, and the number of sampling steps. The second
category of methods modifies/optimize the forward/backward noising process
[32] in a training-free manner, including dynamic programming [35] and
analytical generation process [1]. An important and effective instance is the
Denoising Diffusion Implicit Model (DDIM) [30], a non-Markovian noising
process. Another category of training-free methods modifies the numerical
solver for stochastic differential equations (SDEs) or ordinary differential
equations (ODEs) associated with the DDPMs [13, 33]. In particular, [32]
solves a marginal equivalent ODE known as the Probability Flow (PF) using a
black box ODE solver. [19, 20, 37] compute DDIM with high order methods and
achieve further acceleration. The deterministic DDIM can also be viewed as a
time discretization of the PF as it matches the latter in the continuous limit
[30].

This paper develops a simple and robust fast sampling method for DDPMs with
much fewer steps while retaining high sample quality. Our approach is inspired
by backward error analysis, a method for the numerical analysis of ordinary
differential equation (ODE) solvers [8]. Backward error analysis is a powerful
tool for studying the long-time behavior of numerical ODE solvers. The
critical insight of backward error analysis is that we can describe the bias
introduced when integrating an ODE with finite step sizes by introducing a
modified ancillary flow. This modified flow is derived to ensure that discrete
iterates of the original ODE lie on the path of the continuous solution to the
modified flow.

Our contributions are summarized as follows: 1) The inference process of DDPMs
approximates solving the corresponding diffusion ordinary differential
equations (diffusion ODEs) in the continuous limit [32]. This work analyzes
how the backward error affects the diffusion ODEs and the sample quality in
DDPMs; 2) we propose fast sampling through the Restricting Backward Error
schedule (RBE schedule) based on dynamically moderating the long-time backward
error, our method accelerates DDPMs without any further training; 3) the
experiments show that sampling with an RBE schedule generates high-quality
samples within only 8 to 20 function evaluations on various benchmark
datasets. We achieved 12.01 FID in 8 function evaluations on the ImageNet
$128\times 128$, and a $20\times$ speedup compared with previous baseline
samplers.

Figure 1:  Sampling through RBE schedule within $10$ forward evaluations. We
use the diffusion model [23] pre-trained on on LSUN bedroom 256×256.

##  2 Denoising Diffusion Probabilistic Models

Denoising Diffusion Probabilistic Models (DDPMs)[29, 10] are defined in terms
of a forward Markovian diffusion process $q$ and a learned reverse process
$p$. The forward diffusion process gradually adds Gaussian noise to a data
point $x_{0}$ through $N$ iterations,

$q(x_{1:N}\mid x_{0})=\prod_{i=1}^{N}q(x_{i}\mid x_{i-1}),$

and

$\displaystyle q(x_{i}\mid x_{i-1})$
$\displaystyle=\mathcal{N}(x_{i}\mid\sqrt{\alpha_{i}}x_{i-1},(1-\alpha_{i})\mathbf{I})$
$\displaystyle=\mathcal{N}(x_{i}\mid\sqrt{\frac{\gamma_{i}}{\gamma_{i-1}}}x_{i-1},\left(1-\frac{\gamma_{i}}{\gamma_{i-1}}\right)\mathbf{I}),$
(1)

where $\gamma_{0}=1,\ \gamma_{i}=\prod_{j=1}^{i}\alpha_{j}$. The scalar
parameters $\alpha_{1:N}$ determine the variance of the noise added at each
diffusion step, subject to $0<\alpha_{i}<1$. The learned reverse process aims
to model $q(x_{0})$ by inverting the forward process, gradually removing noise
from the signal starting from isotropic Gaussian noise $x_{N}$,

$\displaystyle x_{N}\sim\mathcal{N}(0,I),$ $\displaystyle
p_{\theta}(x_{0:N})=p(x_{N})\prod_{i=1}^{N}p_{\theta}(x_{i-1}\mid x_{i}),$
$\displaystyle p_{\theta}(x_{i-1}\mid
x_{i})=\mathcal{N}(x_{i-1}\mid\mu_{\theta}(i,x_{i}),\sigma_{i}^{2}I).$ (2)

The parameters $\theta$ of the reverse process can be optimized by maximizing
the following variational lower bound on the training set [10],

$\displaystyle L_{i-1}$
$\displaystyle=\mathbb{E}_{x_{0},\epsilon}\left[\left\|\sqrt{\frac{1}{\alpha_{i}}}\left(x_{i}(x_{0},\epsilon)-\frac{1-\alpha_{i}}{\sqrt{1-\gamma_{i}}}\right)\epsilon\right.\right.$
$\displaystyle\left.\left.\quad-\mu_{\theta}(i,x_{i}(x_{0},\epsilon))\right\|^{2}\right]$
$\displaystyle=\mathbb{E}_{x_{0},\epsilon}\left[\frac{(1-\alpha_{i})^{2}}{\alpha_{i}(1-\gamma_{i})}\left\|\epsilon-\epsilon_{\theta}(i,x_{i}(x_{0},\epsilon))\right\|^{2}\right],$
(3)

where $\gamma_{i}=\prod_{j=1}^{i}\alpha_{j}$,
$x_{i}(x_{0},\epsilon)=\sqrt{\gamma_{i}}x_{0}+\sqrt{1-\gamma_{i}}\epsilon$,
and $\epsilon\sim\mathcal{N}(0,I)$. $\epsilon_{\theta}$ is an estimate of the
noise $\epsilon$. The relationship between $\epsilon_{\theta}$ and
$\mu_{\theta}$ is
$\mu_{\theta}=\frac{1}{\sqrt{\alpha_{i}}}\left(x_{i}-\frac{1-\alpha_{i}}{\sqrt{1-\gamma_{i}}}\epsilon_{\theta}\right)$.

###  2.1 Diffusion SDEs

Before introducing our approach to efficiently sampling from DDPMs, it is
helpful to link DDPMs to continuous time affine diffusion processes, as it
reveals the intuition behind our approach. First we switch the domain of index
$i$ from $\\{0,1,2,...,N\\}$ to continuous interval $[0,T]$, in which
$t_{i}=\frac{iT}{N}$ and each $x_{i}$ is corresponding to the value at moment
$t_{i}$, e.g. $x_{i}=x_{t_{i}}$. As $N$ goes to infinity,
$\\{x_{t_{i}}\\}_{i=0}^{N}$ and $\\{\gamma_{t_{i}}\\}_{i=0}^{N}$ approach
continuous functions $x_{t}$ and $\gamma_{t}$ on interval $[0,T]$. Let
$x_{0}\sim q(x_{0})$ denote a data point drawn from the empirical distribution
of the train set. Let $q(x_{t}\mid x_{0})$ denote a stochastic process for
$t\in[0,T]$ defined through the following stochastic differential equation
(SDE),

$\displaystyle\mathrm{d}x_{t}=\frac{x}{2}\frac{\mathrm{d}\log\gamma_{t}}{\mathrm{d}t}\mathrm{d}t+\sqrt{-\frac{\mathrm{d}\log\gamma_{t}}{\mathrm{d}t}}\mathrm{d}w,$
(4)

where $\epsilon\sim\mathcal{N}(0,I)$, $\mathrm{d}w$ is Wiener process.
Equation 1 matches the differential Equation 4 in the continuous limit.
Meanwhile, the reverse time denoising process also coincides with a time-
continuous stochastic differential Equation [32] (diffusion SDE) starting from
time $T$ to $0$,

$\mathrm{d}x_{t}=\frac{\mathrm{d}\log\gamma_{t}}{\mathrm{d}t}\left(\frac{x}{2}-\frac{\epsilon_{\theta}(t,x_{t})}{\sqrt{1-\gamma_{t}}}\right)\mathrm{d}t+\sqrt{\frac{\mathrm{d}\log\gamma_{t}}{\mathrm{d}t}}\mathrm{d}w.$
(5)

Equation 5 is Variance Preserving stochastic differential equations (VP-SDEs)
[32]. We can generate samples by solving the diffusion SDE with numerical
solvers. This approach discretizes the diffusion SDE from $T$ to $0$. [32]
proved that the classical ancestral sampling method for DDPMs is a first-order
SDE solver for Equation 5. However, these first-order methods usually require
hundreds of or thousands of function evaluations to converge [32, 10],
resulting in extremely slow sampling speed.

###  2.2 Diffusion ODEs (Probability Flow)

By randomness of the Wiener process, the first-order SDE solvers always
require a limited step size [14]. Simulating SDEs with a large step size
(small number of steps) often results in non-convergence, especially in high-
dimensional spaces. For fast sampling, [32] shows that Equation 5 has a
marginal equivalent ODE known as the diffusion ODE (Probability Flow),

$\mathrm{d}x_{t}=\frac{\mathrm{d}\log\gamma_{t}}{\mathrm{d}t}\left(\frac{x}{2}-\frac{\epsilon_{\theta}(t,x_{t})}{2\sqrt{1-\gamma_{t}}}\right)\mathrm{d}t.$
(6)

We can draw samples by solving the diffusion ODE from $T$ to $0$. The
deterministic DDIM can also be viewed as a time discretization of the
diffusion ODEs as it matches the latter in the continuous limit [30]. In
contrast with SDEs, ODEs can be solved with larger step sizes as there is no
Wiener process. Diffusion ODEs show us a new possibility to accelerate the
reverse process. Furthermore, we can use efficient numerical ODE solvers to
accelerate the sampling.

[32] solves the diffusion ODEs using the RK45 ODE solver [6]. This approach
generates samples in 60 function evaluations on the CIFAR-10 dataset [16] and
achieves comparable quality with a 1000-step SDE solver for diffusion SDEs.
[20, 37] apply change of variable and solve diffusion ODEs using the numerical
methods called exponential integrators [24]. Their approaches generate
satisfactory samples in the few-step (10 steps) sampling regime. However, [20]
fails dramatically in the fewer-step (8 steps) sampling regime (Section 4). We
know there is still a lack of training-free samplers for DDPMs in the few-step
sampling regime. In addition, compared to the commonly-used GANs, diffusion
ODE solvers are still not fast enough for real-time applications.

##  3 A Backward Error Analysis of Diffusion ODEs

As we have discussed in Section 2.2, simulating SDEs with discrete updates is
generally difficult in high dimensions [14], and it is hard to converge within
a few steps. In contrast, solving ODEs is relatively easier. This fact yields
a potential for fast samplers and motivates us to design an efficient solver
for diffusion ODEs to enable fast and high-quality few-step sampling.

The critical insight of this work is inspired by Backward error analysis.
Backward error analysis has great potential to clarify the role of finite step
size and to help identify the implicit biases of different ODE solvers. We,
therefore, give a detailed introduction to the core methodology in Section 3.1
before deriving our diffusion ODE solvers in Section 3.2 and Section 3.3.

###  3.1  An Introduction to Backward Error Analysis

In numerical analysis, we often wish to integrate ODEs of the form
$\frac{\mathrm{d}x}{\mathrm{d}t}=f(t,x)$. This system usually cannot be solved
analytically, forcing us to approximate the continuous flow with discrete
updates, e.g., the Euler step

$x_{t+h}=x_{t}+h\cdot f(t,x_{t}).$ (7)

However, discrete updates introduce approximation error when the step size $h$
is finite. In order to study the bias introduced by this approximation error,
we assume the step size $h$ is relatively small, and introduce a modified flow
$\frac{\mathrm{d}x}{\mathrm{d}t}=\tilde{f}(t,x)$,

$\tilde{f}(t,x)=f(t,x)+hf_{1}(t,x)+h^{2}f_{2}(t,x)+...,$ (8)

where $f_{1},f_{2}...$ are the correction terms. The modified flow of
$\tilde{f}(t,x)$ converges to the original flow of $f(t,x)$ as $h$ goes to
$0$, but it differs from the original ODE if $h$ is finite. The goal of
backward error analysis is to choose the correction terms $f_{1}(t,x)$ such
that the iterates obtained from discrete updates of the original flow with
small finite step sizes lie on the path given by the continuous solution to
the modified flow with vanishing step sizes.

The standard derivation of backward error analysis follows by taking a Taylor
expansion in $h$ of the solution to the modified flow $x_{t+h}$ (see [8]). We
obtain the derivatives of $x_{t+h}$ recursively using the modified flow
equation $\frac{\mathrm{d}x}{\mathrm{d}t}=\tilde{f}(t,x)$, and we identify the
correction terms $f_{i}(t,x)$ by ensuring this Taylor expansion matches the
discrete update (e.g., $x_{t+h}=x_{t}+hf(t,x_{t})$ for all powers of $h$. This
section will build an intuition for the origin of the corrections terms. We
will clarify how we might apply this analysis to high-quality few-step
sampling. First, we will identify the path of the continuous modified flow by
solving it analytically in the continuous limit. Then we will compare this
continuous path to the single Euler step $x_{t+h}=x_{t}+h\cdot f(t,x_{t})$.
Imagine tracking the modified flow $\tilde{f}(t,x)$ with time step size $h$,

$\displaystyle x_{t+h}$
$\displaystyle=x_{t}+\int_{t}^{t+h}\tilde{f}(s,x_{s})\mathrm{d}s$
$\displaystyle=x_{t}+h\tilde{f}(t,x_{t})+\frac{h^{2}}{2}\frac{\mathrm{d}\tilde{f}(t,x_{t})}{\mathrm{d}t}+O(h^{3}).$
(9)

We arrived at Equation 9 by taking the Taylor expansion of $\tilde{f}$ then
integrating over time. Here we assume
$\frac{\mathrm{d}\tilde{f}}{\mathrm{d}t}$ does exist. Next, we replace
$\tilde{f}$ with its definition in Equation 8. As we will see below,

$\displaystyle x_{t+h}$
$\displaystyle=x_{t}+hf(t,x_{t})+h^{2}\left(f_{1}(t,x_{t})+\frac{1}{2}\frac{\mathrm{d}f(t,x_{t})}{\mathrm{d}t}\right)+O(h^{3}).$
(10)

Equation 10 is the key component of backward error analysis. It describes the
path taken while integrating the continuous modified flow $\tilde{f}$ over a
discrete time step of length $h$.

Notice that we have assumed that the Taylor expansion in Equation 10
converges, while the higher order terms $O(h^{3})$ will contain higher order
derivatives of the original flow $f$. Backward error analysis implicitly
assumes that $f$ is an analytic function in $t$ and $x$. We refer the reader
to [8] for a detailed introduction.

First order correction: we will now derive the first order correction $f_{1}$
of the modified flow $\tilde{f}$ in Equation 8. Recall that the discrete
updates follow by $x_{t+h}=x_{t}+hf(t,x_{t})$. In order to make sure that the
continuous modified flow coincides with the discrete updates, we need all
terms at $O(h^{2})$ and above in Equation 10 to vanish. As for the second
order terms, it implies that
$f_{1}(t,x_{t})+\frac{1}{2}\frac{\mathrm{d}f(t,x_{t})}{\mathrm{d}t}=O(h)$, and
yields the first order correction,

$\displaystyle f_{1}(t,x)=-\frac{1}{2}\left(\frac{\partial f(t,x)}{\partial
t}+\frac{1}{2}\nabla\|f(t,x)\|^{2}\right),$ (11)

where $\nabla$ is the gradient with respect to $x$.

We conclude that if the time step $h$ is sufficiently small such that we can
neglect higher order terms in Equation 10, then the discrete Euler updating
path Equation 7 lies close to the path of the modified ODE (see Fig. 2),

$\frac{\mathrm{d}x}{\mathrm{d}t}=f(t,x)-\frac{h_{t}}{2}\left(\frac{\partial
f(t,x)}{\partial t}+\frac{1}{2}\nabla\|f(t,x)\|^{2}\right),$ (12)

where $h_{t}$ is the finite step size at moment $t$.

Figure 2: The cartoon illustration for backward error analysis. Iterates
obtained from discrete Euler updates $(w_{k+1}=w_{k}+\epsilon f(w_{k}))$ of
the original flow $(\dot{w}=f(w))$ with small finite step sizes $\epsilon$ lie
close to the path given by the continuous solution to the modified flow
$(\dot{w}=\tilde{f}(w))$ with vanishing step sizes.

###  3.2  Dynamically Restricting the Backward Error

We now derive our first algorithm. Recall that the flow of the diffusion ODEs
satisfies,

$\frac{\mathrm{d}x}{\mathrm{d}t}=\frac{\mathrm{d}\gamma_{t}}{\mathrm{d}t}\left(\frac{x}{2\gamma_{t}}-\frac{\epsilon_{\theta}(t,x)}{2\gamma_{t}\sqrt{1-\gamma_{t}}}\right).$
(13)

By applying change-of-variable, we switch to consider the following ODE,

$\frac{\mathrm{d}x}{\mathrm{d}\gamma}=f_{\theta}(\gamma,x):=\frac{x}{2\gamma}-\frac{\epsilon_{\theta}(\gamma,x)}{2\gamma\sqrt{1-\gamma}},$
(14)

where $\gamma$ is between $0$ and $1$, and
$\epsilon_{\theta}(\gamma_{t},x)=\epsilon_{\theta}(t,x)$. The deterministic
DDIM can be viewed as a discrete version of Equation 14 as it matches the
latter in the continuous limit [30]. Followed by Equation 12, the discrete
flow lies close to the path of the modified ODE,

$\displaystyle\frac{\mathrm{d}x}{\mathrm{d}\gamma}$
$\displaystyle=\tilde{f_{\theta}}(\gamma,x)$
$\displaystyle=f_{\theta}(\gamma,x)-\frac{h_{\gamma}}{2}\left(\frac{\partial
f_{\theta}(\gamma,x)}{\partial\gamma}+\frac{1}{2}\nabla\|f_{\theta}(\gamma,x)\|^{2}\right),$
(15)

where $h_{\gamma}$ is the finite step size at the moment $\gamma$. If the
modified flow $\tilde{f_{\theta}}$ coincides with the original Diffusion ODE
$f_{\theta}$ at $\gamma$, as the moment goes from $\gamma$ to
$\gamma+h_{\gamma}$, the discrepancy between the modified flow and the
original flow approximately becomes

$\frac{h^{2}_{\gamma}}{2}\left(\frac{\partial
f_{\theta}(\gamma,x)}{\partial\gamma}+\frac{1}{2}\nabla\|f_{\theta}(\gamma,x)\|^{2}\right).$
(16)

Observing the discrepency 16 motivates us to design a solver for diffusion
ODEs to enable high-quality few-step sampling. Our first algorithm dynamically
crafts the step size $h_{\gamma}$ to maintain restricting the backward error.
Below we present our first method and name it as sampling through _Dynamically
Restricting the Backward Error (DRBE)_ schedule.

Algorithm 1 Improving the sampling quality by dynamically restricting the
backward error

threshold $r>0$; $f_{\theta}(\cdot,\cdot)$; initial state $x_{T}$.

$t\leftarrow T$

$x\sim\mathcal{N}(0,I)$

$\gamma\leftarrow\gamma_{t}$

$sche\leftarrow\\{t\\}$

while $t>0$ do

$h\leftarrow-\min\left\\{1-\gamma,\ \sqrt{\frac{r}{\left|\frac{\partial
f_{\theta}(\gamma,x)}{\partial\gamma}+\frac{1}{2}\nabla\|f_{\theta}(\gamma,x)\|^{2}\right|}}\right\\}$

$x\leftarrow x+hf_{\theta}(\gamma,x)$

$\gamma\leftarrow\gamma+h$

$t\leftarrow t^{-1}(\gamma)$

$sche\leftarrow sche\cup\\{t\\}$

end while

return $x,sche$

In Algorithm 1, $t^{-1}(\gamma)$ is the inverse function of $\gamma_{t}$.
$t^{-1}(\cdot)$ has an analytical formulation for both linear and cosine noise
schedules used in [23, 10]. The crafted threshold $r>0$ governs the sample
quality and sampling speed trade-offs. The average inference steps
increase/decrease if we accept a small/large threshold. We carefully fine-tune
the threshold to achieve high-quality sampling within 8, 10, 12, 15, and 20
steps.

###### Remark 1

Followed by Equation 16, the step size $h_{\gamma}$ is governed by the
magnitude of the correction term $\frac{\partial
f(\gamma,x)}{\partial\gamma}+\frac{1}{2}\nabla\|f(\gamma,x)\|^{2}$. The
correction term encodes the second-order derivatives information of the
diffusion ODEs. Intuitively, in numerical analysis, we can simulate the ODEs
with small/large step sizes while the curvature of the flow is
significant/minor.

###  3.3 Restricting the backward error learns effective inference schedule

As shown in the previous section, sampling through DRBE schedule requires one
function evaluation and one gradient evaluation per step. Despite the more
expensive steps, we observe a fun fact in practical experiments: even though
the initial states $x$ are randomly sampled from isotropic Gaussian
distribution, the output $sche$ in Algorithm 1 is roughly unchanged if we fix
the inference model $\theta$. We suspect that this phenomenon may stem from
the flow of diffusion ODEs being locally parallel. This observation motivates
us to design a cheaper but equally effective approach for fast sampling
compared with DRBE schedule. Suppose we have generated a small batch of $N$
samples using Algorithm 1 within $K$ inference steps. Our second method first
learns a newly crafted inference schedule by restricting the backward error
and taking the average of the collected $K-$steps schedules
$\\{sche^{i}\\}_{i=1}^{N}$. More precisely, for $k=1,2,...,K+1$, we compute
the following statistics

$t_{k}=\frac{1}{N}\sum_{i=1}^{N}sche^{i}[k],$ (17)

here $sche^{i}[k]$ is the $k$-th item in $sche^{i}$. We name
$\\{t_{1},t_{2},...,t_{K+1}\\}$ as _Restricting Backward Error schedule (RBE
schedule)_. Our second algorithm utilizes the learnt inference schedule
$\\{t_{1},t_{2},...,t_{K+1}\\}$ and only requires one function evaluation per
step, as shown in below,

Algorithm 2 Improving the sampling speed through the RBE schedule

$\gamma_{t}$; $f_{\theta}(\cdot,\cdot)$; RBE schedule
$\\{t_{1},t_{2},...,t_{K}\\}$.

$x\sim\mathcal{N}(0,I)$

for $k\leftarrow 1$ to $K$ do:

$h\leftarrow\gamma_{t_{k+1}}-\gamma_{t_{k}}$

$x\leftarrow x+hf_{\theta}(\gamma_{t_{k}},x)$

end for

return $x$

In Section 4, we show that sampling through the RBE schedule 2 achieves high-
quality sampling with in few function evaluations.

##  4 Experiments

In this section, we show that generating samples through DRBE schedule 1 and
RBE schedule 2 can significantly speed up the sampling of existing pre-trained
DDPMs, with both linear noise models [10] and cosine noise models [23]. We
vary the number of function evaluations (NFE), which is the number of calls to
the noise prediction model $\epsilon_{\theta}(t,x)$, and compare the sample
quality among RBE schedule, DRBE schedule, and other methods. For each
experiment, We draw $50K$ samples and use the widely adopted Fenchel Inception
Distance (FID [9]) to evaluate the sample quality, lower FID usually indicates
better sample quality. We refer to the Appendix for further implementation
details of the experiments.

###  4.1 Comparison with other training-free samplers

We generate samples through DRBE and RBE schedules using pre-trained DDPMs and
then compare our results with other DDPMs samplers, including DDPM [10],
DDIM[30], and Dpm-solvers [20]. We compare the sample quality by varying NFE
from $8$ to $1000$.

Specifically, we use the DDPMs trained by $L_{hybrid}$ in [23] on ImageNet
64x64 [4] with cosine noise schedule; the classifier guidance DDPMs in [5] on
ImageNet 128x128 [4] with linear noise schedule, whereas the classifier
guidance scale is $1.0$; the discrete-time model in [23] on LSUN 256x256 [36]
pre-trained with linear noise schedule. For the models trained on ImageNet, we
only use their “mean” model and omit the “variance” model.

[20] evidence that DPM-solvers present the superior efficiency in the few-step
regime around 10 to 15 NFE, where other samples such as GGDM [34], Analytic-
DDIM [1], exponential integrator [37] and Itô-Taylor [33] have rather large
discretization errors. However, our experiments show that Dpm-solvers [20]
fail dramatically in ”fewer” step regimes such as 8 NFE. Meanwhile, sampling
through our DEBE and RBE schedule shows the stability and superiority in a
broader range of NFEs. As shown in Fig 3, 4 and 5, on both ImagesNet and LSUN
datasets, sampling through DEBE and RBE schedule can obtain satisfactory
samples within 8 steps (FID 29.49 on ImageNet 64x64, FID 12.01 on ImageNet
128x128, and FID 27.18 on LSUN bedroom 256 x 256). Sampling with DRBE and RBE
schedule outperforms all the previous samplers in the regime of $8$ inference
steps, and achieves comparable performance with the previous fastest training-
free samplers such as DPM-solvers [20] in the regime of 12 to 20 inference
steps.

Figure 3: Sample quality measured by FID $\downarrow$ on ImageNet 64×64 with
the diffusion model [23], varying the number of function evaluations (NFE): 8,
10, 12, 15, 20, 50, 200, 1000. For DDIM and DDPM, we use uniform re-spacing
time steps for all the experiments. Sampling through our RBE schedule shows
uniform stability and superiority in a few inference step regimes. Our method
achieves FID 25.82 in 10 function evaluations and at least 20$\times$ speedup
compared with DDPM.  Figure 4: Sample quality measured by FID $\downarrow$ on
ImageNet 128×128 with the diffusion model [5], varying the number of function
evaluations (NFE): 8, 10, 12, 15, 20, 50, 100. For DDIM and DDPM, we use
uniform re-spacing time steps for all the experiments. Sampling through RBE
schedule achieves FID 12.01 within only 8 function evaluations.  Figure 5:
Sample quality measured by FID $\downarrow$ on LSUN bedroom 256×256 with the
diffusion model [23], varying the number of function evaluations (NFE): 8, 10,
12, 50, 200, 1000. For DDIM and DDPM, we use uniform re-spacing time steps for
all the experiments. Sampling through our RBE schedule achieves FID 17.9 in 10
function evaluations and at least 20$\times$ speedup compared with DDPM.

###  4.2  RBE schedule lies between linear and cosine noise schedule

As derived in section 3.3, RBE schedule is learned by restricting the backward
error through the sampling process. We have shown that sampling with RBE
schedule brings the efficiency of DDPMs samplers to a new level by producing
high-quality samples in the “few-step sampling” regime. In this section, we
empirically study the behavior of the RBE schedule. In particular, we
visualize the RBE schedule in Fig. 6 on both ImageNet 64×64 and ImageNet
128x128 experiments introduced in the previous section. We found that the RBE
schedule always lies between the popular choices of noise schedules, such as
linear and cosine schedules. RBE schedule is close to the linear noise
schedule in the 8 inference steps regime. On the other hand, RBE schedule
gradually approaches the cosine noise schedule as NFE increases. This
observation motivates us to future research on the design of the optimal noise
schedule.

Figure 6:  Comparing our RBE schedule with the linear and cosine schedule. The
$x-$axis represents the diffusion steps from $0$ to $T$, while the $y-$axis
represents the variance schedule $\gamma_{t}$. The left figure plots the RBE
schedules of ImageNet 64×64 experiments with the diffusion model [23], while
the right one presents the RBE schedules of ImageNet 128×128 experiments with
the diffusion model [5]. RBE schedule is closer to the linear noise schedule
in the 8-step regime. Meanwhile, RBE schedule gradually approaches the cosine
noise schedule as NFE increases.

##  5 Related work

Backward Error Analysis Backward error analysis is a method for the numerical
analysis of ordinary differential equation (ODE) solvers [8]. It is a powerful
tool for studying the long-time behavior of numerical ODE solvers. In the deep
learning community, [28, 2] analyzed the influence of finite learning rates on
the iterates of gradient descent (GD) and shed light on explaining why
stochastic gradient descent finds suitable solutions.

DDPMs samplers Several works aim at speeding up the inference process for
DDPMs. One category of methods includes knowledge distillation [26, 21],
learning sampling trajectory, and searching for better noise schedule [27, 17,
23, 34]. Such methods require a training stage before the sampling stage.
Hence their applicability and effectivity might be limited when adapting their
methods to different models, datasets, and the number of sampling steps. The
second category of methods modifies/optimize the forward/backward noising
process [32] in a training-free manner, including dynamic programming [35] and
analytical generation process [1]. An important and effective instance is the
Denoising Diffusion Implicit Model (DDIM) [30], a non-Markovian noising
process. Another category of training-free methods modifies the numerical
solver for stochastic differential equations (SDEs) or ordinary differential
equations (ODEs) associated with the DDPMs [13, 33]. In particular, [32]
solves a marginal equivalent ODE known as the Probability Flow (PF) using a
black box ODE solver. [19, 20, 37] compute DDIM with high order methods and
achieve further acceleration. The deterministic DDIM can also be viewed as a
time discretization of the PF as it matches the latter in the continuous limit
[30].

.

## References

  * [1] Fan Bao, Chongxuan Li, Jun Zhu, and Bo Zhang.  Analytic-dpm: an analytic estimate of the optimal reverse variance in diffusion probabilistic models.  arXiv preprint arXiv:2201.06503, 2022.
  * [2] David GT Barrett and Benoit Dherin.  Implicit gradient regularization.  arXiv preprint arXiv:2009.11162, 2020.
  * [3] Nanxin Chen, Yu Zhang, Heiga Zen, Ron J Weiss, Mohammad Norouzi, and William Chan.  Wavegrad: Estimating gradients for waveform generation.  arXiv preprint arXiv:2009.00713, 2020.
  * [4] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei.  Imagenet: A large-scale hierarchical image database.  In 2009 IEEE conference on computer vision and pattern recognition, pages 248–255. Ieee, 2009.
  * [5] Prafulla Dhariwal.  Diffusion models beat gans on image synthesis.  NeurIPS, 34:8780–8794, 2021.
  * [6] John R Dormand and Peter J Prince.  A family of embedded runge-kutta formulae.  Journal of computational and applied mathematics, 6(1):19–26, 1980\.
  * [7] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  Advances in neural information processing systems, 27, 2014.
  * [8] Ernst Hairer, Christian Lubich, and Gerhard Wanner.  Geometric numerical integration illustrated by the störmer–verlet method.  Acta numerica, 12:399–450, 2003.
  * [9] Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter.  Gans trained by a two time-scale update rule converge to a local nash equilibrium.  Advances in neural information processing systems, 30, 2017.
  * [10] Jonathan Ho.  Denoising diffusion probabilistic models.  NeurIPS, 33:6840–6851, 2020.
  * [11] Jonathan Ho and Tim Salimans.  Classifier-free diffusion guidance.  arXiv preprint arXiv:2207.12598, 2022.
  * [12] Jonathan Ho, Tim Salimans, Alexey Gritsenko, William Chan, Mohammad Norouzi, and David J Fleet.  Video diffusion models.  arXiv preprint arXiv:2204.03458, 2022.
  * [13] Alexia Jolicoeur-Martineau, Ke Li, Rémi Piché-Taillefer, Tal Kachman, and Ioannis Mitliagkas.  Gotta go fast when generating data with score-based models.  arXiv preprint arXiv:2105.14080, 2021.
  * [14] Peter E Kloeden and Eckhard Platen.  Stochastic differential equations.  In Numerical solution of stochastic differential equations, pages 103–160. Springer, 1992.
  * [15] Zhifeng Kong, Wei Ping, Jiaji Huang, Kexin Zhao, and Bryan Catanzaro.  Diffwave: A versatile diffusion model for audio synthesis.  arXiv preprint arXiv:2009.09761, 2020.
  * [16] Alex Krizhevsky, Geoffrey Hinton, et al.  Learning multiple layers of features from tiny images.  2009\.
  * [17] Max WY Lam, Jun Wang, Rongjie Huang, Dan Su, and Dong Yu.  Bilateral denoising diffusion models.  arXiv preprint arXiv:2108.11514, 2021.
  * [18] Haoying Li, Yifan Yang, Meng Chang, Shiqi Chen, Huajun Feng, Zhihai Xu, Qi Li, and Yueting Chen.  Srdiff: Single image super-resolution with diffusion probabilistic models.  Neurocomputing, 479:47–59, 2022.
  * [19] Luping Liu, Yi Ren, Zhijie Lin, and Zhou Zhao.  Pseudo numerical methods for diffusion models on manifolds.  arXiv preprint arXiv:2202.09778, 2022.
  * [20] Cheng Lu.  DPM-Solver: A fast ODE solver for diffusion probabilistic model sampling in around 10 steps.  2022.
  * [21] Eric Luhman and Troy Luhman.  Knowledge distillation in iterative generative models for improved sampling speed.  arXiv preprint arXiv:2101.02388, 2021.
  * [22] Zhaoyang Lyu, Zhifeng Kong, Xudong Xu, Liang Pan, and Dahua Lin.  A conditional point diffusion-refinement paradigm for 3d point cloud completion.  arXiv preprint arXiv:2112.03530, 2021.
  * [23] Alexander Nichol.  Improved denoising diffusion probabilistic models.  In ICML, pages 8162–8171. PMLR, 2021.
  * [24] Alexander Ostermann and M Hochbruck.  Exponential integrators.  Acta Numer, 19:209–286, 2010.
  * [25] Chitwan Saharia, Jonathan Ho, William Chan, Tim Salimans, David J Fleet, and Mohammad Norouzi.  Image super-resolution via iterative refinement.  arXiv preprint arXiv:2104.07636, 2021.
  * [26] Tim Salimans and Jonathan Ho.  Progressive distillation for fast sampling of diffusion models.  arXiv preprint arXiv:2202.00512, 2022.
  * [27] Robin San-Roman, Eliya Nachmani, and Lior Wolf.  Noise estimation for generative diffusion models.  arXiv preprint arXiv:2104.02600, 2021.
  * [28] Samuel L Smith, Benoit Dherin, David GT Barrett, and Soham De.  On the origin of implicit regularization in stochastic gradient descent.  arXiv preprint arXiv:2101.12176, 2021.
  * [29] Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli.  Deep unsupervised learning using nonequilibrium thermodynamics.  In International Conference on Machine Learning, pages 2256–2265. PMLR, 2015.
  * [30] Jiaming Song.  Denoising diffusion implicit models.  2020.
  * [31] Yang Song and Stefano Ermon.  Generative modeling by estimating gradients of the data distribution.  Advances in Neural Information Processing Systems, 32, 2019.
  * [32] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole.  Score-based generative modeling through stochastic differential equations.  arXiv preprint arXiv:2011.13456, 2020.
  * [33] Hideyuki Tachibana, Mocho Go, Muneyoshi Inahara, Yotaro Katayama, and Yotaro Watanabe.  It$\backslash$^$\\{$o$\\}$-taylor sampling scheme for denoising diffusion probabilistic models using ideal derivatives.  arXiv preprint arXiv:2112.13339, 2021.
  * [34] Daniel Watson.  Learning fast samplers for diffusion models by differentiating through sample quality.  In ICML, 2021.
  * [35] Daniel Watson, Jonathan Ho, Mohammad Norouzi, and William Chan.  Learning to efficiently sample from diffusion probabilistic models.  arXiv preprint arXiv:2106.03802, 2021.
  * [36] Fisher Yu, Ari Seff, Yinda Zhang, Shuran Song, Thomas Funkhouser, and Jianxiong Xiao.  Lsun: Construction of a large-scale image dataset using deep learning with humans in the loop.  arXiv preprint arXiv:1506.03365, 2015.
  * [37] Qinsheng Zhang and Yongxin Chen.  Fast sampling of diffusion models with exponential integrator.  arXiv preprint arXiv:2204.13902, 2022.