that its perturbed expansion at $\mathcal{C}$ is positive,
$\delta\vartheta>0$. To find the apparent horizon at first order, we
infinitesimally perturb $\mathcal{C}$ inwards by $\psi l^{a}$ and denote the
first order change of the outgoing null-expansion $\vartheta$ by $\delta_{\psi
l}\vartheta$. If we now choose $\psi$ so that $\delta_{\psi
l}\vartheta=-\delta\vartheta$, then to first order, this uniquely locates an
apparent horizon lying along the flow of $\psi l^{a}$ from $\mathcal{C}$. To
first order, the area of the apparent horizon will not depend on the choice of
$l^{a}$.

The computation of the location of the apparent horizon is most conveniently
done using Gaussian null coordinates based upon the cross-section
$\mathcal{C}$ and the Killing parametrization, as previously introduced in
(145). We take $l^{a}=(\partial/\partial r)^{a}$, so $l^{a}$ is the (past
directed) ingoing null normal to $\mathcal{C}$, normalized so that
$l^{a}\xi_{a}=1$. The perturbed expansion of $\mathcal{C}$ along the horizon
direction with respect to $\xi^{a}$ is then $\delta\tilde{\vartheta}=\kappa
V\delta\vartheta$. The equation for $\psi$ that we must solve is

$\delta_{\psi l}\tilde{\vartheta}=-\delta\tilde{\vartheta}=-\kappa
V\delta\vartheta.$ (170)

The variation $\delta_{\psi l}\tilde{\vartheta}=:W\psi$ is given by a linear
operator $W$ intrinsic to $\mathcal{C}$ acting on $\psi$ which has been given
in [38] (see also eqs. (5.12) and (5.34) of [39]) as181818The differences
between formula (171) for our operator defined by $\delta_{\psi
l}\tilde{\vartheta}=:W\psi$ and the corresponding eq. (1) of [38] result from
the facts that (i) we are considering a null displacement rather than a
spacelike displacement and (ii) they set $\tilde{\vartheta}=0$ since they are
considering a marginally outer trapped surface. Equation (171) is valid for
the change in the expansion of the orthogonal null geodesics of any $n-2$
dimensional surface $\mathcal{C}$ under an infinitesimal displacement of the
surface by $\psi l^{a}$, where $l^{a}$ is the other null normal field to
$\mathcal{C}$, normalized to have unit inner product with the null normal
defining the expansion.

$W\psi=-D^{a}D_{a}\psi+\beta^{a}D_{a}\psi+\Big{(}G_{ab}\xi^{a}l^{b}+\frac{1}{2}(R[\gamma]-\frac{1}{2}\beta^{a}\beta_{a}+D_{a}\beta^{a})-\tilde{\vartheta}\overline{\tilde{\vartheta}}\Big{)}\psi$
(171)

where $R[\gamma]$ denotes the scalar curvature of the induced derivative
$D_{a}$ on $\mathcal{C}$ associated with the induced metric $\gamma_{ab}$,
$\tilde{\vartheta}=\frac{1}{2}\gamma^{cd}{\mathscr{L}}_{\xi}\gamma_{cd}$ is
the outgoing expansion of $\mathcal{C}$, and
$\overline{\tilde{\vartheta}}=\frac{1}{2}\gamma^{cd}{\mathscr{L}}_{l}\gamma_{cd}$
the ingoing expansion of $\mathcal{C}$.

We now use Einstein’s equation $G_{ab}=0$ and the fact that $\mathcal{H}$ has
vanishing outgoing shear and expansion in the stationary background. Using the
fact that $\tilde{\alpha}~{}\hat{=}~{}\kappa$ and using Einstein’s equation
again, the scalar curvature of $\mathcal{C}$ is found to be (see also eq. (82)
of [32])

$\displaystyle R[\gamma]={}$ $\displaystyle
D_{a}\beta^{a}+\frac{1}{2}\beta^{a}\beta_{a}+2\kappa\overline{\tilde{\vartheta}}.$
(172)

Thus, $\psi$ is determined by

$-D^{a}D_{a}\psi+\beta^{a}D_{a}\psi+(D_{a}\beta^{a})\psi+\kappa\overline{\tilde{\vartheta}}\psi=-\kappa
V\delta\vartheta.$ (173)

The second and third terms on the left side combine to a total divergence
$D_{a}(\psi\beta^{a})$ intrinsic to $\mathcal{C}$. Integrating over
$\mathcal{C}$, we see that $\psi$ satisfies the relation

$\kappa\int_{\mathcal{C}}\overline{\tilde{\vartheta}}\psi\,\bm{\epsilon}^{(n-2)}=-\kappa\int_{\mathcal{C}}V\delta\vartheta\,\bm{\epsilon}^{(n-2)}.$
(174)

On the other hand, the change in area resulting from displacing $\mathcal{C}$
infinitesimally along $\psi l^{a}$ is

$\delta
A=\int_{\mathcal{C}}\psi\overline{\tilde{\vartheta}}\,\bm{\epsilon}^{(n-2)}.$
(175)

Thus, we see that the area of the apparent horizon differs from the area of
$\mathcal{C}$ by

$\delta
A=\int_{\mathcal{C}}\psi\overline{\tilde{\vartheta}}\,\bm{\epsilon}^{(n-2)}=-\int_{\mathcal{C}}V\delta\,\vartheta\bm{\epsilon}^{(n-2)}.$
(176)

Comparing with (89), we see that up to first order, we have

$S[\mathcal{C}]=\frac{A_{\rm app}[\mathcal{C}]}{4}$ (177)

where $A_{\rm app}[\mathcal{C}]$ is the area of the apparent horizon at the
time corresponding to the cross-section $\mathcal{C}$ of the horizon.

##  Appendix B Generalization of Theorem 1 to the Case of a Tensor Dotted with
Two $\xi$’s

Consider the hypotheses of Theorem 1, but with $\alpha^{a}$ in (51) replaced
by

$\displaystyle\beta^{ab}={}$
$\displaystyle\sum_{i=0}^{k}T_{(i)}^{abc_{1}\cdots
c_{i}de}\nabla_{(c_{1}}\cdots\nabla_{c_{i})}\delta g_{de}.$ (178)

We wish to obtain a strengthening of the conclusion (52) of Theorem 1 that is
valid for the form of $\beta_{ab}\xi^{a}\xi^{b}$. Specifically, we will show
that $\beta_{ab}\xi^{a}\xi^{b}$ can be written in the form

$\displaystyle\beta_{ab}\xi^{a}\xi^{b}~{}\hat{=}~{}{}$
$\displaystyle\sum_{i=0}^{k-2}\tilde{T}_{(i)}^{c_{1}\cdots
c_{i}de}\nabla_{(c_{1}}\cdots\nabla_{c_{i})}(\mathscr{L}_{\xi}-\kappa)\mathscr{L}_{\xi}\delta
g_{de}$ (179)

where the tensors $\tilde{T}_{(i)}^{b_{1}\cdots b_{i}cd}$ are smooth on
$\mathcal{H}$ and are locally and covariantly constructed from the metric,
curvature, and its derivatives as well as from $\xi^{a}$, and $N^{a}$, with
$\xi^{a}$ and $N^{a}$ appearing only algebraically.

To show this, we start with the highest derivative term in
$\beta_{ab}\xi^{a}\xi^{b}$ and proceed as in the proof of Theorem 1. Each term
in the basis expansion of $\xi_{a}\xi_{b}T_{(k)}^{abc_{1}\cdots c_{i}de}$
corresponding to (53) will now have at least two “extra” $\xi$’s, i.e., there
will be at least 2 more $\xi$’s than $N$’s in each term. We can perform the
exactly the same manipulations with one of these “extra” $\xi$’s as was done
to derive (61) in the proof of Theorem 1. We thereby obtain

$\xi_{a}\xi_{b}T_{(k)}^{abc_{1}\cdots
c_{k}de}\nabla_{(c_{1}}\cdots\nabla_{c_{k})}\delta
g_{de}~{}\hat{=}~{}\tilde{T}_{(k-1)}^{c_{1}\cdots
c_{k-1}de}\nabla_{(c_{1}}\cdots\nabla_{c_{k-1})}{\mathscr{L}}_{\xi}\delta
g_{de}+\sum_{i=0}^{k-1}{T^{\prime}}_{(i)}^{c_{1}\cdots
c_{i}de}\nabla_{(c_{1}}\cdots\nabla_{c_{i})}\delta g_{de}$ (180)

with $\tilde{T}_{(k-1)}^{c_{1}\cdots c_{k-1}de}$ and
${T^{\prime}}_{(i)}^{c_{1}\cdots c_{i}de}$ satisfying the same properties as
enumerated in Theorem 1 but with $\tilde{T}_{(k-1)}^{c_{1}\cdots c_{k-1}de}$
now satisfying the additional property that there is at least one “extra”
$\xi^{a}$ appearing in each term in its basis expansion. We may now perform
the same manipulations with this additional $\xi^{a}$ to obtain

$\xi_{a}\xi_{b}T_{(k)}^{abc_{1}\cdots
c_{k}de}\nabla_{(c_{1}}\cdots\nabla_{c_{k})}\delta
g_{de}~{}\hat{=}~{}\tilde{T}_{(k-2)}^{c_{1}\cdots
c_{k-2}de}\nabla_{(c_{1}}\cdots\nabla_{c_{k-2})}{\mathscr{L}}^{2}_{\xi}\delta
g_{de}+\sum_{i=0}^{k-1}{T^{\prime\prime}}_{(i)}^{c_{1}\cdots
c_{i}de}\nabla_{(c_{1}}\cdots\nabla_{c_{i})}\delta g_{de}$ (181)

where ${T^{\prime\prime}}_{(i)}^{c_{1}\cdots c_{i}de}$ incorporates all of the
additional lower derivative terms arising from this second set of
manipulations. Now, the left side of (181) is $O(V^{2})$ as $V\to 0$, where
$V$ denotes an affine parameter on the horizon (with $V=0$ corresponding to
the bifurcation surface). However, the first term on the right side is only
$O(V)$ since

$\displaystyle\tilde{T}_{(k-2)}^{c_{1}\cdots
c_{k-2}de}\nabla_{(c_{1}}\cdots\nabla_{c_{k-2})}{\mathscr{L}}^{2}_{\xi}\delta
g_{de}={}$
$\displaystyle{\mathscr{L}}^{2}_{\xi}\left[\tilde{T}_{(k-2)}^{c_{1}\cdots
c_{k-2}de}\nabla_{(c_{1}}\cdots\nabla_{c_{k-2})}\delta g_{de}\right]$
$\displaystyle={}$ $\displaystyle\kappa V{\mathscr{L}}_{k}\left(\kappa
V{\mathscr{L}}_{k}\left[\tilde{T}_{(k-2)}^{c_{1}\cdots
c_{k-2}de}\nabla_{(c_{1}}\cdots\nabla_{c_{k-2})}\delta g_{de}\right]\right)$
$\displaystyle={}$
$\displaystyle\kappa^{2}(V^{2}{\mathscr{L}}^{2}_{k}+V{\mathscr{L}}_{k})\left[\tilde{T}_{(k-2)}^{c_{1}\cdots
c_{k-2}de}\nabla_{(c_{1}}\cdots\nabla_{c_{k-2})}\delta g_{de}\right]$ (182)

and the factor with $V{\mathscr{L}}_{k}$ is only $O(V)$. We can rectify this
by replacing ${\mathscr{L}}^{2}_{\xi}\delta g_{de}$ in (181) by
$[{\mathscr{L}}^{2}_{\xi}-\kappa{\mathscr{L}}_{\xi}]\delta g_{de}$ and
compensating for this by adding the quantity $\tilde{T}_{(k-2)}^{c_{1}\cdots
c_{k-2}de}\nabla_{(c_{1}}\cdots\nabla_{c_{k-2})}\kappa{\mathscr{L}}_{\xi}\delta
g_{de}$ to the second term on the right side of (181). Since this quantity
involves at most $(k-1)$ derivatives of $\delta g_{de}$, it does not change
the character of this term. Thus, have shown that

$\xi_{a}\xi_{b}T_{(k)}^{abc_{1}\cdots
c_{k}de}\nabla_{(c_{1}}\cdots\nabla_{c_{k})}\delta
g_{de}~{}\hat{=}~{}\tilde{T}_{(k-2)}^{c_{1}\cdots
c_{k-2}de}\nabla_{(c_{1}}\cdots\nabla_{c_{k-2})}({\mathscr{L}}_{\xi}-\kappa){\mathscr{L}}_{\xi}\delta
g_{de}+\sum_{i=0}^{k-1}{T^{\prime\prime\prime}}_{(i)}^{c_{1}\cdots
c_{i}de}\nabla_{(c_{1}}\cdots\nabla_{c_{i})}\delta g_{de}$ (183)

where $\tilde{T}_{(k-2)}$ and ${T^{\prime\prime\prime}}_{(i)}$ are smooth and
are locally and covariantly constructed from the metric, curvature, $\xi^{a}$,
and $N^{a}$, with only algebraic dependence on $\xi^{a}$ and $N^{a}$. But now
the left side and the first term on the right side are both $O(V^{2})$ as
$V\to 0$. Therefore, the last term on the right side is also $O(V^{2})$ and
therefore ${T^{\prime\prime\prime}}_{(k-1)}^{c_{1}\cdots c_{i}de}$ must also
have two “extra” $\xi$’s in each term of its basis expansion. This allows us
to make an inductive argument to prove (179) in parallel with the proof of
Theorem 1.

######  Remark.

For the case of a tensor

$\displaystyle\gamma^{a_{1}\cdots a_{p}}={}$
$\displaystyle\sum_{i=0}^{k}T_{(i)}^{a_{1}\cdots a_{p}c_{1}\cdots
c_{i}de}\nabla_{(c_{1}}\cdots\nabla_{c_{i})}\delta g_{de},$ (184)

contracted into $p$ factors of $\xi^{a}$, the corresponding result is

$\displaystyle\gamma_{a_{1}\cdots
a_{p}}\xi^{a_{1}}\cdots\xi^{a_{p}}~{}\hat{=}~{}{}$
$\displaystyle\sum_{i=0}^{k-p}\tilde{T}_{(i)}^{c_{1}\cdots
c_{i}de}\nabla_{(c_{1}}\cdots\nabla_{c_{i})}(\mathscr{L}_{\xi}-\kappa(p-1))\cdots(\mathscr{L}_{\xi}-2\kappa)(\mathscr{L}_{\xi}-\kappa)\mathscr{L}_{\xi}\delta
g_{de}$ (185)

where the tensors $\tilde{T}_{(i)}^{c_{1}\cdots c_{i}de}$ are smooth on
$\mathcal{H}$ and are locally and covariantly constructed from the metric,
curvature, and its derivatives as well as from $\xi^{a}$, and $N^{a}$, with
$\xi^{a}$ and $N^{a}$ appearing only algebraically.

##  Appendix C Ambiguities of $S$ Beyond First Order

###  C.1 Non-dynamical Background Structures

As discussed in section IV.1, the definition of $S$ suffers from ambiguities
at second order and beyond. One source of ambiguities arises from the presence
of $\xi^{a}$ in our formula for $\textbf{B}_{\mathcal{H}}$ and the fact that,
going beyond linear order in perturbation theory, we can no longer treat
$\xi^{a}$ as the exact Killing field on the stationary background. The
approach taken in section IV.1 was to tie the metric $g_{ab}$ to $\xi^{a}$ by
postulating in effect that $g_{ab}$ is in a Gaussian null gauge (81) and is
related to $\xi^{a}$ by (82), in such a way that $\mathcal{B}$ is at
$\rho=0=V$.

In this section, we will present another possible approach wherein the class
of metrics $g_{ab}$ considered (i.e. effectively the gauge near
$\mathcal{H}$), as well as the form of $\xi^{a}$ off of $\mathcal{H}$ is
determined by a non-dynamical background structure given only _on
$\mathcal{H}$,_ as opposed to _off of $\mathcal{H}$_ as in (82). In
particular, off of $\mathcal{H}$, though not on $\mathcal{H}$, the $\xi^{a}$
considered in this appendix will depend on $g_{ab}$, unlike in section IV.1.

Our main reason for considering a different approach here is that with the
approach given below, it will be possible to make a precise comparison between
the ambiguities and structure of $S$ in our approach and that of recent
proposal by [17] (see section C.3). In fact, while [17] also works in a
Gaussian null gauge, just as we did to tie the gauge of $g_{ab}$ to $\xi^{a}$
in section IV.1, the gauge considered by [17] is tied to the cross section
$\mathcal{C}$ on which one wants to define $S[{\mathcal{C}}]$. By contrast, in
the present paper we have $\xi^{a}$ which is not tied from the outset to a
particular $\mathcal{C}$ on which we anticipate evaluating $S$.

In order to characterize the ambiguities in $S$ precisely, it is necessary to
specify the geometric structures underlying our construction. The rigidly
fixed, non-dynamical background structure considered in this appendix consists
of (a) an $(n-1)$-dimensional hypersurface $\mathcal{H}$ of $\mathcal{M}$
having the topological structure of a trivial fibre bundle with fibres
$\mathbb{R}$ and compact base; (b) a preferred cross section, $\mathcal{B}$,
on $\mathcal{H}$; (c) a vector field $\xi^{a}$ tangential to the fibres which
changes its direction at $\mathcal{B}$ (and vanishes only on $\mathcal{B}$);
(d) a positive real number $\kappa$; (e) a vector field $N^{a}$ transverse to
$\mathcal{H}$ that is defined for each point in
${\mathcal{H}}\setminus{\mathcal{B}}$, and is non-zero where defined.

######  Definition 1.

Given our rigidly fixed, non-dynamical background structure
$({\mathcal{H}},{\mathcal{B}},\xi^{a},N^{a},\kappa)$, we consider compatible
spacetime metrics $g_{ab}$ on $\mathcal{M}$ in the following sense. We require
that $\mathcal{H}$ is null with normal $\xi^{a}$ and, furthermore, is such
that $\xi^{a}\nabla_{a}\xi^{b}=\kappa\xi^{b},$ $N^{a}N_{a}=0$, and
$g_{ab}\xi^{b}N^{a}=1$ on $\mathcal{H}$. We call this space of metrics
$\mathcal{P}$. $\textbf{B}_{\mathcal{H}}$, as well as S are functionals on
$\mathcal{P}$.

The space $\mathcal{P}$ does not imply any physical restrictions onto the
metric. Indeed, let $\tilde{\mathcal{H}}$ be any null surface with compact
cross sections ruled by affinely parameterized null geodesics with tangent
$k^{a}$. Let $V$ be any affine parameter along the null geodesics vanishing on
some cut, $\mathcal{B}$ and set $\tilde{\xi}^{a}=\kappa Vk^{a}$. Then we have
$\tilde{\xi}^{a}\nabla_{a}\tilde{\xi}^{b}=\kappa\tilde{\xi}^{b}$ and
$\tilde{\xi}^{a}$ is obviously null and normal to $\tilde{\mathcal{H}}$.
Furthermore, let $\tilde{N}^{a}$ be a second null field defined on
$\tilde{\mathcal{H}}$ such that $\tilde{\xi}^{a}\tilde{N}_{a}=1$. By applying
a diffeomorphism $\phi$, we can clearly achieve that
$\phi^{*}\tilde{N}^{a}=N^{a},\phi^{*}\tilde{\xi}^{a}=\xi^{a}$ and
$\phi[\tilde{\mathcal{H}}]=\mathcal{H}$. Thus, $\phi^{*}g_{ab}$ is in
$\mathcal{P}$.

Due to the derivatives acting on $\xi^{a}$ in (69), we first define $N^{a}$
off of $\mathcal{H}$ by demanding that, besides relation (186) and
$N^{a}\xi_{a}=1,N^{a}N_{a}=0$ on $\mathcal{H}$, we have
$N^{a}\nabla_{a}N^{b}=0$. Since $\xi^{a}$ is hypersurface orthogonal, we have
$\nabla_{[a}\xi_{b]}~{}\hat{=}~{}2\kappa w_{[a}\xi_{b]}$ for some $w_{a}$. In
fact, we can and will adjust the extension such that $w_{a}=N_{a}$. To see
this, note that we are allowed to change $\xi^{a}\to\xi^{a}+fp^{a}$ where $f$
is any function on $\mathcal{M}$ vanishing on $\mathcal{H}$ and $p^{a}$ is any
vector field. Since $\xi^{a}$ is normal to $\mathcal{H}$, we must have
$\nabla_{a}f=c\xi_{a}$ on $\mathcal{H}$, where $c$ is a function that we can
choose as we please. This means that $w_{a}\to w_{a}-cp_{a}/\kappa$ which we
use to change $w_{a}$ to $N_{a}$.

Our requirement $\nabla_{[a}\xi_{b]}~{}\hat{=}~{}2\kappa N_{[a}\xi_{b]}$ fixes
the first derivative of $\xi^{a}$ on $\mathcal{H}$ only partially. Note that
we are not free to fix the symmetrized derivative of $\xi^{a}$ e.g. by (33)
since that relation used the Killing property of $\xi^{a}$, and we are no
longer requiring at this stage that $\xi^{a}$ is a Killing vector field with
respect to $g_{ab}.$ However, we may – and will – demand without loss of
generality that $\xi^{a}$ has been extended off of $\mathcal{H}$ such that

$\nabla_{[a}\xi_{b]}~{}\hat{=}~{}2\kappa
N_{[a}\xi_{b]},\quad(\nabla_{a}\xi_{b})N^{a}N^{b}~{}\hat{=}~{}0,\quad(\nabla_{(a}\xi_{b)})\xi^{b}~{}\hat{=}~{}0.$
(186)

We are allowed to change $\xi^{a}\to\xi^{a}+f\xi^{a}$, where $f$ is any
function on $\mathcal{M}$ vanishing on $\mathcal{H}$, without affecting the
first relation in (186). Since $\nabla_{a}f=c\xi_{a}$ on $\mathcal{H}$, we can
thereby change $\nabla_{(a}\xi_{b)}\to\nabla_{(a}\xi_{b)}+c\xi_{a}\xi_{b}$.
Since $c$ is a function that we can choose as we please and since
$N^{a}\xi_{a}=1$, we can clearly impose
$(\nabla_{a}\xi_{b})N^{a}N^{b}~{}\hat{=}~{}0$. The condition
$(\nabla_{(a}\xi_{b)})\xi^{b}~{}\hat{=}~{}0$ easily follows from
$\nabla_{[a}\xi_{b]}~{}\hat{=}~{}2\kappa N_{[a}\xi_{b]}$ and
$\xi^{a}\nabla_{a}\xi^{b}~{}\hat{=}~{}\kappa\xi^{b}$. Finally, using
$N^{a}\nabla_{a}N^{b}=0$ and
$(\nabla_{(a}\xi_{b)})\xi^{b}~{}\hat{=}~{}0~{}\hat{=}~{}(\nabla_{a}\xi_{b})N^{a}N^{b}$,
we have

$N^{a}\nabla_{a}(N^{b}\xi_{b})~{}\hat{=}~{}0,\quad\nabla_{a}(\xi^{b}\xi_{b})~{}\hat{=}~{}-2\kappa\xi_{a}.$
(187)

Transvecting $\nabla_{[a}\xi_{b]}~{}\hat{=}~{}2\kappa N_{[a}\xi_{b]}$ with
$N^{a}$, we get ${\mathscr{L}}_{N}\xi_{a}~{}\hat{=}~{}-2\kappa N_{a}$, which
clearly uniquely determines the first order change of $\xi^{a}$ off of
$\mathcal{H}$.

As we have argued, (186) defines $\xi^{a}$ to first order off of
$\mathcal{H}$, but we would like to define it to arbitrary orders off of
$\mathcal{H}$. In order to have a definition that gives us the Killing vector
field $\xi^{a}$ in the special case that $\mathcal{H}$ is a Killing horizon
for $g_{ab}$, we require $\xi^{a}$ to satisfy the geodesic deviation equation,

$N^{a}\nabla_{a}(N^{b}\nabla_{b}\xi_{c})=R_{cbad}N^{b}N^{a}\xi^{d},$ (188)

which is a viewed as a second order ordinary differential equation for
$\xi^{a}$ whose initial conditions are posed on $\mathcal{H}$. Since the
geodesic deviation equation implies that
$N^{a}\nabla_{a}(N^{b}\nabla_{b}(\xi^{c}N_{c}))=0$, the initial conditions
$\xi^{a}N_{a}~{}\hat{=}~{}1,N^{c}\nabla_{c}(N^{a}\xi_{a})~{}\hat{=}~{}0$,
imply that $N^{a}\xi_{a}=1$ in a neighborhood of $\mathcal{H}$.

Consider now a metric $\tilde{g}_{ab}$ having a stationary bifurcate Killing
horizon $\mathcal{H}$ with Killing field $\tilde{\xi}^{a}$ that vanishes on
$\mathcal{B}$. We have constructed in equations (33) to (38) a vector field
$\tilde{N}^{a}$ at points of $\mathcal{H}$ from the Killing field
$\tilde{\xi}^{a}$ and the metric $\tilde{g}_{ab}$. We can align these with
$\xi^{a},N^{a}$ as in definition 1 by a diffeomorphism $\phi$, so that
$g_{ab}=\phi^{*}\tilde{g}_{ab}\in\mathcal{P}$. Since any Killing field
satisfies (63), we automatically have
$N^{a}\nabla_{a}(N^{b}\nabla_{b}\xi_{c})=R_{cbad}N^{b}N^{a}\xi^{d}$, where the
derivatives and Riemann tensor refer to $g_{ab}$. Thus, our extension
procedure produces precisely the horizon Killing field if there is one.

We likewise require that any variation $\delta g_{ab}$ be compatible with our
the fixed (not varying) background structure
$({\mathcal{H}},{\mathcal{B}},\xi^{a},N^{a},\kappa)$ in this appendix, in the
sense that it is given as the derivative of a 1-parameter family of metrics in
the space $\mathcal{P}$ (definition 1), i.e. as an element in the tangent
space of $\mathcal{P}$. Then we have $\delta g_{ab}\xi^{a}=0$ on
$\mathcal{H}$. Furthermore, since $\delta\kappa=0$ as $\kappa$ is part of the
rigidly fixed background structure, the last equation in (187) implies that
$0=\delta[\nabla_{a}(\xi^{b}\xi_{b})]=\nabla_{a}(\delta
g_{bc}\xi^{c}\xi^{b})$. Thus, the previously considered gauge conditions (42),
(43) hold generally, not just off a stationary (bifurcate Killing-) horizon.

Our main first order variation formula for $S$ was (76) for any perturbation
$\delta g_{ab}$ satisfying gauge conditions (42), (43). In the construction
leading to (76), we previously considered $\xi^{a}$ as fixed not only on
$\mathcal{H}$, but also off of $\mathcal{H}$. By contrast, the algorithm we
have just given above for extending $\xi^{a},N^{a}$ off of $\mathcal{H}$
depended on the chosen metric, and therefore while
$\delta\xi^{a}~{}\hat{=}~{}0$, the derivatives of $\delta\xi^{a}$ off of
$\mathcal{H}$ will not in general vanish. Nevertheless, by a lengthy argument,
we have shown that (76) continues to hold. As a consequence we still have the
formulas expressing the second law obtained within the previous framework of
non-dynamical background structure in sections V.1 and V.2. Similarly, we
still have the second order variation formulas (110) and (111) in the present
setting with a varying $\xi^{a}$, with only minor modifications to the
definition of the boundary terms in the modified canonical energy (111).
However, we shall not present the details of these modifications here.

###  C.2 Structure and Ambiguities of $S$

As in the approach taken in section IV.1, there is substantial potential
ambiguity in $S$ beyond first order. The first ambiguity arises because it may
happen that two metrics $g_{ab}$ and $\phi^{*}g_{ab}$ may both be in
$\mathcal{P}$, where $\phi$ is a non-trivial diffeomorphism fixing
$\mathcal{H}$ and hence $\xi^{a}$ on $\mathcal{H}$. We therefore have
$(\phi^{*}g_{ab})\xi^{a}N^{b}=1$ by definition 1, so $g_{ab}$ satisfies the
definition 1 also for $\tilde{N}^{a}=\phi^{-1*}N^{a}$. Since $g_{ab}$ and
$\phi^{*}g_{ab}$ must be regarded as physically equivalent, this is telling us
that we have, in effect, changed $N^{a}\to\tilde{N}^{a}$ while keeping
$g_{ab}$ fixed, analogous to the freedom of choosing a different affine
parameterization in of $\mathcal{H}$ in section IV.1.

The second ambiguity arises from the fact that $\textbf{B}_{\mathcal{H}}$,
which enters $S$, was constructed to satisfy (70), which is a condition that
applies only to first order perturbations. We may add many terms to
$\textbf{B}_{\mathcal{H}}$ (and hence to
$\textbf{S}=\textbf{Q}-\xi\cdot\textbf{B}_{\mathcal{H}}$) that vanish
identically at first order but not second. The following lemma will enable us
to give a precise classification of such terms.

######  Lemma 1.

Let $\textbf{B}_{\mathcal{H}}$ be an $(n-1)$-form valued functional satisfying
1.–5. in section IV.1. Then it can be written as a sum of monomials of the
following factors i)–iii), contracted with $N^{a},\xi^{a},s^{a}_{i}$’s and
multiplied by powers of $\kappa$.

  * i)

$s_{i_{1}}^{a_{1}}\cdots
s_{i_{j}}^{a_{j}}\nabla_{(a_{1}}\cdots\nabla_{a_{j})}{\mathscr{L}}_{\xi}g_{ab}$;

  * ii)

$s_{i_{1}}^{a_{1}}\cdots
s_{i_{j}}^{a_{j}}\nabla_{(a_{1}}\cdots\nabla_{a_{j-1})}({\mathscr{L}}_{\xi})^{s}N_{a_{j}}$,
$s\geq 0$;

  * iii)

$\nabla_{(a_{1}}\cdots\nabla_{a_{j}}R^{b}{}_{a_{j+1}a_{j+2})}{}^{c},j>0$,
$R_{abcd}$.

Each contracted monomial in $\textbf{B}_{\mathcal{H}}$ is invariant under a
rescaling $\xi^{a}\to c\xi^{a}$, $N^{a}\to N^{a}/c$, $\kappa\to c\kappa$.
There must be at least one factor of type i), ii) with $s>0$ or a fully
contracted factor of type iii) with more $\xi^{a}$’s than $N^{a}$’s. The total
number of derivatives in each monomial must be less than the number of
derivatives in the Lagrangian, with each occurrence of $\kappa$ counting as
one derivative.

###### Proof.

By assumptions 1.–5. $\textbf{B}_{\mathcal{H}}$ is a sum of monomials
consisting of the following factors: (a)
$\nabla_{a_{1}}\dots\nabla_{a_{j}}N_{b}$, (b)
$\nabla_{a_{1}}\dots\nabla_{a_{j}}\xi_{b}$, (c)
$\nabla_{a_{1}}\dots\nabla_{a_{j}}R_{abcd}$, which are fully contracted into
the legs of the tetrad $N^{a},\xi^{a},s_{i}^{a}$. We will now show by
induction in $j\geq 0$ that each such term can be converted into one of the
terms i)–iii) listed in lemma 1. For term c) it follows from the results of
[40]. Thus, we only have to deal with the terms (a) and (b). For $j=0$, there
is nothing to show.

Assuming that the statement is true up to $j-1$ derivatives, let us prove it
for $j$ derivatives. We begin by considering all possible ways of dotting the
legs of our null tetrad into a term of the form (b),
$\nabla_{a_{1}}\dots\nabla_{a_{j}}\xi_{b}$. By commuting derivative operators
at the expense of Riemann tensors (giving terms of the form (b) with fewer
derivatives and terms of the form (c)), and by moving any $N^{a}$ factors into
the derivative operators at the expense of new terms of the form (a) with
fewer than $j$ derivatives, we may bring the expression into the form
$\nabla_{a_{1}}\dots\nabla_{a_{j-k}}(N^{c}\nabla_{c})^{k}\xi_{b}$, dotted into
$\xi^{a}$’s and $s_{i}^{a}$’s. For $k\geq 2$, we may use (188) to lower $k$ at
the expense of terms of the form (c), so only the case $k=0,1$ needs to be
considered. As an example of such a term for $k=1$, consider
$s_{i_{1}}^{a_{1}}\dots
s_{i_{s}}^{a_{j}}\xi^{a_{s+1}}\dots\xi^{a_{j-1}}\nabla_{a_{1}}\dots\nabla_{a_{j-1}}(N^{c}\nabla_{c}\xi_{b})$.
At the expense of terms of the form (a) when a derivative hits $N^{c}$, we can
move all derivatives tangent to $\xi^{d}$ onto $\xi_{b}$ and then use
$\xi^{d}\nabla_{d}\xi_{b}~{}\hat{=}~{}\kappa\xi_{b}$, since all derivatives
are at this stage tangent to $\mathcal{H}$. Thus, we can effectively assume
that the term we are dealing with is of the form

$\begin{split}\text{term (b)}&=s_{i_{1}}^{a_{1}}\dots
s_{i_{j-1}}^{a_{j-1}}\nabla_{a_{1}}\dots\nabla_{a_{j-1}}(N^{c}\nabla_{c}\xi_{b})\\\
&=s_{i_{1}}^{a_{1}}\dots
s_{i_{j-1}}^{a_{j-1}}\nabla_{a_{1}}\dots\nabla_{a_{j-1}}(N^{c}(\tfrac{1}{2}{\mathscr{L}}_{\xi}g_{cb}+2\kappa
N_{[c}\xi_{b]})).\end{split}$ (189)

The first term on the right side is of the form (i) in the lemma, whereas the
second has fewer than $j$ derivatives acting onto $\xi_{b}$ and at most $j$
derivatives acting onto $N^{c}$, i.e. it is a term of the form (a), after
symmetrizing derivatives at the expense of terms of the form (c) with Riemann
tensors. The case $k=0$ is treated similarly.

Next, we consider all possible ways of dotting the legs of our null tetrad
into a term of the form (a), $\nabla_{a_{1}}\dots\nabla_{a_{j}}N_{b}$. By
commuting derivative operators at the expense of Riemann tensors (giving terms
of the form (a) with fewer derivatives and terms of the form (c)), and by
moving any $N^{a}$ factors into the derivative operators at the expense of new
terms of the form (a) with fewer than $j$ derivatives, we may bring the
expression into the form
$\nabla_{a_{1}}\dots\nabla_{a_{j-k}}(N^{c}\nabla_{c})^{k}N_{b}$, dotted into
$\xi^{a}$’s and $s_{i}^{a}$’s. Such a term vanishes for $k>0$ because of
$N^{c}\nabla_{c}N_{b}=0$. We can move all derivatives tangent to $\xi^{d}$
onto $N_{b}$ at the expense of terms of the form (b) with fewer than $j$
derivatives. Thus, we can effectively assume that the term we are dealing with
is of the form

$\text{term (a)}=s_{i_{1}}^{a_{1}}\dots
s_{i_{j-s}}^{a_{j-s}}\xi^{a_{j-s+1}}\dots\xi^{a_{j}}\nabla_{a_{1}}\dots\nabla_{a_{j}}N_{b}.$
(190)

We can now gradually express the number of derivatives into the
$\xi^{a}$-direction in terms of Lie derivatives ${\mathcal{L}}_{\xi}$ and
terms such that each factor has fewer than $j$ derivatives. So, we have
reduced the (a) terms to terms of the form

$\text{term (a)}=s_{i_{1}}^{a_{1}}\dots
s_{i_{j-s}}^{a_{j-s}}\nabla_{a_{1}}\dots\nabla_{a_{j-s}}({\mathcal{L}}_{\xi})^{s}N_{b},$
(191)

and we can symmtrize the $a_{i}$-indices at the expense of terms of the type
(c). If we dot $\xi^{b}$ into this expression, we get terms of type (b), since
$\xi^{b}N_{b}=1$. If we dot $N^{b}$ into this expression, we get terms of the
type (a) with fewer derivatives since $N^{b}N_{b}=0$ in a neighborhood of
$\mathcal{H}$. So we have effectively reduced attention to terms of the form
claimed in item ii). ∎

The main consequence of this lemma is as follows: If $\mathcal{H}$ is a
stationary bifurcate Killing horizon with Killing vector field $\xi^{a}$, then
any factor of type i), ii) with $s>0$ or any fully contracted factor of type
iii) with more $\xi^{a}$’s than $N^{a}$’s vanishes. By contrast, a factor of
type ii) with $s=0$, or a fully contracted factor of type iii) with no more
$\xi^{a}$’s than $N^{a}$’s factors will not, in general, vanish. It is thus
clear that the ambiguity in $\textbf{B}_{\mathcal{H}}$ consists in adding
monomials with at least _two_ factors of either type i), ii) with $s>0$, or a
fully contracted factor of type iii) with more $\xi^{a}$’s than $N^{a}$’s.

###  C.3 Comparison with HKR Proposal

In [17] Hollands, Kovacs, and Reall (HKR) constructed an entropy-current
$(n-2)$-form $\textbf{S}_{\rm HKR}$ in an “effective field theory” (EFT)
framework

$L=\frac{1}{8\pi}(2\ell^{-2}\Lambda+R+\ell^{2}L_{4}+\ell^{4}L_{6}+\dots),$
(192)

where each $L_{N}$ is a Lagrangian that is locally and covariantly constructed
out of the metric and $N$ derivatives. EFT is understood to mean that one
restricts attention to only those solutions in the theory truncated at some
given order $N$ such that, roughly speaking, the higher order terms191919The
cosmological constant term is also considered “higher order”, in the sense
that $\Lambda\ll 1$, cf. the cosmological constant problem. in the equations
of motion are locally small near $\mathcal{H}$. $\textbf{S}_{\rm HKR}$ is by
construction equal to the Dong–Wall entropy-current $(n-2)$-form [5, 6]
$\textbf{S}_{\rm DW}$ to first order off a solution with a stationary
bifurcate Killing horizon. [17, 21] showed that the ambiguities can also be
exploited, for each $N$, so that $\textbf{S}_{\rm HKR}$ satisfies the 2nd law
in the EFT sense, i.e. up to terms that are of the same order in $\ell$ as the
terms neglected when truncating the EFT order $N$.

A comparison between the ambiguities in the HKR approach and our approach may
now be given as follows: In [17, lem. 2.3] it was shown that the admissible
terms in $\textbf{S}_{\rm HKR}$ must be monomials in Gaussian null coordinate
components (see (81)) of the following factors of total boost weight zero (see
[17, def. 2.2]):

  * i’)

$D_{(a_{1}}\cdots D_{a_{j})}K_{ab}$;

  * ii’)

$D_{(a_{1}}\cdots D_{a_{j-1})}\beta_{a_{j}}$, $D_{(a_{1}}\cdots
D_{a_{j})}\bar{K}_{ab}$;

  * iii’)

$\nabla_{(a_{1}}\cdots\nabla_{a_{j}}R^{b}{}_{a_{j+1}a_{j+2})}{}^{c}$,
$R_{abcd}$;

where $n^{a}=(\partial_{\rho})^{a}$. To compare with lemma 1, we have

$\begin{split}N^{a}~{}\hat{=}~{}{}&\kappa^{-1}[\nabla^{a}\log
V+\tfrac{1}{2}\beta^{a}-\tfrac{1}{8}V\beta^{b}\beta_{b}k^{a}],\\\
\xi^{a}~{}\hat{=}~{}&\kappa Vk^{a},\\\
s_{i}^{a}~{}\hat{=}~{}{}&e_{i}^{a}-\tfrac{1}{2}Ve_{i}^{b}\beta_{b}k^{a},\\\
\tfrac{1}{2}{\mathscr{L}}_{\xi}g_{ab}~{}\hat{=}~{}{}&\kappa
VK_{ab}-\tfrac{1}{4}V^{3}\kappa
K_{cd}\beta^{c}\beta^{d}\,k_{a}k_{b},\end{split}$ (193)

where $\sum_{i}e_{i}^{a}e_{i}^{b}=\gamma^{ab}$. The scaling requirement in
lemma 1 is equivalent to the zero total boost weight requirement202020Here it
must be understood that the boost weights of $\kappa,V$ be $+1,-1$
respectively. Note that in the HKR scheme, no explicit factors of $\kappa$ nor
$V$ appear. in [17]. Furthermore, comparing i’)–iii’) in the light of (193)
with the terms i)–iii) in lemma 1, one can see some broad similarities between
iii) and iii’) and i) and i’), if we consider $x^{A}$ Gaussian null coordinate
components as analogous to $s^{a}_{i}$ tetrad components. Items ii) and ii’)
are also similar, in the following sense. If in ii), we have $s=0$ (no Lie-
derivative in the $\xi^{a}$-direction), then one can see from equations (193)
that we basically get a specific combination of the two types of terms in
ii’). When $s\geq 1$ in the term ii), we also get certain combinations of the
two types of terms in ii’). To see this, we first consider one Lie-derivative
in the $\xi^{a}$-direction on $N_{b}$ (193), to get

${\mathscr{L}}_{\xi}N_{a}~{}\hat{=}~{}\tfrac{1}{2}V\mathscr{L}_{k}\beta_{a}+\tfrac{1}{4}V^{2}(\beta^{c}\mathscr{L}_{k}\beta_{c}+K_{bc}\beta^{b}\beta^{c})k_{a}.$
(194)

This already is of the type of terms of the form i’), ii’) except for the term
involving $\mathscr{L}_{k}\beta_{a}$. However, by equations given in [17, sec.
A.1], we have

$\mathscr{L}_{k}\beta_{a}+K_{ab}\beta^{b}~{}\hat{=}~{}2R_{bcda}n^{b}k^{c}k^{d},$
(195)

so $\mathscr{L}_{\xi}N_{a}$ is in fact a combination of the terms in i’),
ii’), iii’) not involving $\bar{K}_{ab}$. This remains true by an inductive
argument for $(\mathscr{L}_{\xi})^{s}N_{a}$, and more generally for all our
terms of type ii) in lemma 1. Thus, our terms i), ii), iii) in lemma 1
correspond to specific combinations of the terms i’),ii’),iii’).

Thus, we see that the ambiguities present in our approach are more restrictive
than in the HKR approach. Correspondingly, our $S[\mathcal{C}]$ is more
rigidly fixed than $S_{\rm HKR}[\mathcal{C}]$. We do not believe that we would
have sufficient re-definition freedom in our approach to ensure that our
entropy flux (113) satisfies the 2nd law beyond general relativity. On the
other hand, it may still be possible to ensure that the second law holds for
specific EFTs, which may thus provide an interesting selection criterion to
discriminate unphysical EFTs.

It appears very unlikely that the more general freedom allowed by HKR will, in
general, yield an entropy $S_{\rm HKR}[\mathcal{C}]$ that is cross section
continuous [19]. To appreciate why, recall that, as shown in [17, sec. 2.1],
under a change of affine parameter $V=V^{\prime}\psi(x^{A})$, where $\psi>0$,
the quantities associated with a corresponding Gaussian null coordinate system
(see (81)) change as

$\begin{split}\beta^{\prime}_{a}~{}\hat{=}~{}&\,\beta_{a}+2D_{a}\log\psi-2VK_{a}{}^{b}D_{b}\log\psi,\\\
\bar{K}_{ab}^{\prime}~{}\hat{=}~{}&\,\psi^{-1}({\bar{K}}_{ab}-VD_{a}D_{b}\log\psi-V(D_{a}\log\psi)D_{b}\log\psi-V\beta_{(a}D_{b)}\log\psi),\quad
K_{ab}^{\prime}~{}\hat{=}~{}\,\psi K_{ab},\end{split}$ (196)

where $K_{ab}$ and $\bar{K}_{ab}$ are the extrinsic curvatures along $V$
respectively $\rho$. Let $\mathcal{C}$ be a cut of constant $V$ and
$\mathcal{C}^{\prime}$ a cut of constant $V^{\prime}$ and consider a $\psi$
that is extremely close to 1 but with extremely large angular derivatives (see
fig. 1). As an example, consider now the Ricci-squared theory, with Dong–Wall
entropy given by (151). It is easy to see that the first term in (151) is
cross section covariant because the definition of the Ricci tensor does not
involve any notion of foliation. Let us next focus on the $K\bar{K}$-term in
(151). Using (196), it is seen to change as

$K^{\prime}\bar{K}^{\prime}~{}\hat{=}~{}K\bar{K}-V^{\prime}D_{a}D^{a}\psi-V^{\prime}\beta^{a}D_{a}\psi$
(197)

under a change of affine parameter. Therefore, since
$\gamma_{ab}^{\prime}~{}\hat{=}~{}\gamma_{ab}$, we obtain

$\int_{{\mathcal{C}}^{\prime}}K^{\prime}\bar{K}^{\prime}\bm{\epsilon}^{\prime}{}^{(n-2)}=\int_{{\mathcal{C}}^{\prime}}\left(K\bar{K}-V^{\prime}D_{a}D^{a}\psi-V^{\prime}\beta^{a}D_{a}\psi\right)\bm{\epsilon}^{(n-2)}.$
(198)

The terms in this last expression involving derivatives of $\psi$ threaten to
blow up. However, these terms crucially are linear in $\psi$, and since
$V^{\prime}$ is constant on $\mathcal{C}^{\prime}$, we may pull it out of the
integral and then integrate the angular derivatives $D_{a}$ by parts. Hence,
it is clear that also the $K\bar{K}$ term in the Dong–Wall entropy given by
(151) is cross section continuous. On the other hand, it seems highly unlikely
that $\textbf{S}_{\rm HKR}$ would only produce terms that are linear in
$\psi$, since, in particular, it contains non-covariant terms starting at
sufficiently high EFT order [22]. Since our entropy is obtained from an
entropy $(n-2)$-form that is covariant in the metric and $\xi^{a}$, it
automatically yields a cross-section continuous entropy.

##  Appendix D Application to Dilaton Gravity Theories with Scalar Fields

Dong and Lewkowycz [30] have analyzed a dilaton gravity theory coupled to two
scalar fields and have obtained an entropy formula for this theory using
Dong’s approach. In this appendix we show that the linearization of their
entropy expression agrees with what would be calculated from (119). This lends
further support to the belief that Dong’s entropy agrees in general with the
entropy that would be obtained from Wall’s approach.

The Lagrangian of dilaton gravity in two dimensions coupled with two scalar
fields, $\sigma$ and $\omega$, considered in [30] is

$\displaystyle{L}_{ab}={}$
$\displaystyle-\frac{1}{2}{\epsilon}_{ab}\Big{[}\phi
R-(\nabla\sigma)^{2}-(\nabla\omega)^{2}+\lambda\omega\nabla_{c}\nabla_{d}\sigma\nabla^{c}\nabla^{d}\sigma\Big{]},$
(199)

where $\lambda$ is a constant. The equations of motion for $g_{ab}$ are

$\displaystyle\begin{split}({E}_{G})_{ab}={}&-\frac{1}{2}\Big{\\{}\nabla_{a}\sigma\nabla_{b}\sigma+\nabla_{a}\omega\nabla_{b}\omega+\nabla_{a}\nabla_{b}\phi-
g_{ab}(g^{ed}\nabla_{e}\nabla_{d}\phi+\frac{1}{2}[(\nabla\sigma)^{2}+(\nabla\omega)^{2}-\lambda\omega(\nabla\nabla\sigma)^{2}])\\\
{}&+\lambda[\nabla^{m}(\omega\nabla_{m}\nabla_{b}\sigma)\nabla_{a}\sigma+(a\leftrightarrow
b)]-\lambda\nabla^{c}(\omega\nabla_{a}\nabla_{b}\sigma\nabla_{c}\sigma)\Big{\\}}\end{split}$
(200)

The pullback of the constraints to the horizon can be put in the form

$\displaystyle\xi^{a}{{C}}_{ab_{2}}={}$ $\displaystyle
2{\epsilon}^{(1)}_{b_{2}}[\kappa\mathscr{L}_{\xi}P-\mathscr{L}_{\xi}^{2}P]$
(201)

where

$\displaystyle P(\psi,\xi;\delta\psi)={}$
$\displaystyle\frac{1}{2}\kappa(\delta\phi-\lambda\omega\mathscr{L}_{\xi}\delta\sigma
N^{a}\nabla_{a}\sigma)$ (202)

where $\psi$ collectively refers to the dynamical fields
$g,\phi,\sigma,\omega$.

The Noether charge 0-form is

$\displaystyle Q={}$
$\displaystyle-\frac{1}{2}{\epsilon}_{ae}A^{aebd}\nabla_{[b}\xi_{d]}+\Big{[}{\epsilon}_{ae}\nabla_{f}A^{fbae}+{\epsilon}_{af}{B}^{a(fb)}\Big{]}\xi_{b}$
(203)

where

$\displaystyle A^{abcd}={}$
$\displaystyle-\frac{1}{2}\phi(g^{ac}g^{bd}-g^{ad}g^{bc})$ (204)
$\displaystyle\tilde{B}^{abc}={}$
$\displaystyle\frac{1}{2}\nabla_{d}\phi(g^{ac}g^{bd}-g^{ad}g^{bc})+\frac{1}{2}\lambda\omega\nabla_{(m}\nabla_{n)}\sigma\nabla_{t}\sigma(2g^{ma}g^{nc}g^{bt}-g^{mb}g^{nc}g^{at}).$
(205)

The symplectic potential $\bm{\theta}$ is given by

$\displaystyle\begin{split}{\theta}_{b_{2}}={}&{\epsilon}_{ab_{2}}\Big{(}A^{abcd}\nabla_{d}\delta
g_{bc}+\tilde{B}^{abc}\delta
g_{bc}+D^{a}\delta\sigma+F^{af}\nabla_{f}\delta\sigma+M^{a}\delta\omega\Big{)}\end{split}$
(206)

where

$\displaystyle D^{a}={}$ $\displaystyle\nabla^{a}\sigma+\lambda
g^{me}g^{na}\nabla_{e}(\omega\nabla_{(m}\nabla_{n)}\sigma)$ (207)
$\displaystyle F^{af}={}$ $\displaystyle-\lambda\omega
g^{ma}g^{nf}\nabla_{(m}\nabla_{n)}\sigma$ (208) $\displaystyle M^{a}={}$
$\displaystyle\nabla^{a}\omega$ (209)

The pullback of $\bm{\theta}$ to the horizon takes the form

$\displaystyle\underline{\bm{\theta}}~{}\hat{=}~{}{}$
$\displaystyle\delta\textbf{B}_{\mathcal{H}}$ (210)

where

${B}_{\mathcal{H}b}=\frac{1}{2}{\epsilon}^{(1)}_{b}\Big{(}2\lambda\omega
N^{a}\nabla_{a}\sigma\mathscr{L}_{\xi}\sigma\Big{)}$ (211)

Our entropy 0-form ${S}$ (see (71)) is given by

$\displaystyle S={}$
$\displaystyle\frac{2\pi}{\kappa}(Q-\xi\cdot\textbf{B}_{\mathcal{H}})$ (212)
$\displaystyle={}$
$\displaystyle\frac{2\pi}{\kappa}[-\kappa\phi+\kappa\lambda\omega
N^{a}\nabla_{a}\sigma\mathscr{L}_{\xi}\sigma-\mathscr{L}_{\xi}(-\phi+\lambda\omega
N^{a}\nabla_{a}\sigma\mathscr{L}_{\xi}^{2}\sigma)]$ (213)

The Dong–Wall entropy (see (119)) is

$\displaystyle S_{\rm DW}={}$
$\displaystyle\frac{2\pi}{\kappa}(Q-\xi\cdot\textbf{B}_{\mathcal{H}})-\frac{4\pi}{\kappa}P(\psi,\xi;\mathscr{L}_{\xi}\psi)$
(214) $\displaystyle={}$ $\displaystyle 2\pi[-\phi+\lambda\omega
N^{a}\nabla_{a}\sigma\mathscr{L}_{\xi}\sigma],$ (215)

This agrees212121To compare, one must take into account—as one also must do in
comparing eq. (14) of [6] with eq. (1.3) of [5]—that Dong works in Euclidean
signature and uses a complex $z$ coordinate. with the linearization of (A.44)
of [30].

## References

  * Bardeen  _et al._ [1973] J. M. Bardeen, B. Carter, and S. W. Hawking, The four laws of black hole mechanics, Communications in Mathematical Physics 31, 161 (1973).
  * Hawking [1975] S. W. Hawking, Particle creation by black holes, Communications in Mathematical Physics 43, 199 (1975).
  * Bekenstein [1973] J. D. Bekenstein, Black holes and entropy, Phys. Rev. D 7, 2333 (1973).
  * Iyer and Wald [1994] V. Iyer and R. M. Wald, Some properties of the noether charge and a proposal for dynamical black hole entropy, Phys. Rev. D 50, 846 (1994).
  * Dong [2014] X. Dong, Holographic entanglement entropy for general higher derivative gravity, Journal of High Energy Physics 2014, 44 (2014).
  * Wall [2015] A. C. Wall, A second law for higher curvature gravity, International Journal of Modern Physics D 24, 1544014 (2015), https://doi.org/10.1142/S0218271815440149 .
  * Wald and Zoupas [2000] R. M. Wald and A. Zoupas, General definition of “conserved quantities” in general relativity and other theories of gravity, Phys. Rev. D 61, 084027 (2000).
  * Gao and Wald [2001] S. Gao and R. M. Wald, “physical process version” of the first law and the generalized second law for charged and rotating black holes, Phys. Rev. D 64, 084020 (2001).
  * Wald [1984] R. M. Wald,  _General Relativity_ (Chicago Univ. Pr., Chicago, USA, 1984).
  * Prabhu [2017] K. Prabhu, The first law of black hole mechanics for fields with internal gauge freedom, Classical and Quantum Gravity 34, 035011 (2017).
  * Hollands and Wald [2013] S. Hollands and R. M. Wald, Stability of black holes and black branes, Communications in Mathematical Physics 321, 629 (2013).
  * Iyer and Wald [1995] V. Iyer and R. M. Wald, Comparison of the noether charge and euclidean methods for computing the entropy of stationary black holes, Phys. Rev. D 52, 4430 (1995).
  * Seifert and Wald [2007] M. D. Seifert and R. M. Wald, General variational principle for spherically symmetric perturbations in diffeomorphism covariant theories, Phys. Rev. D 75, 084029 (2007).
  * Râcz and Wald [1992] I. Râcz and R. M. Wald, Extensions of spacetimes with killing horizons, Classical and Quantum Gravity 9, 2643 (1992).
  * Rácz and Wald [1996] I. Rácz and R. M. Wald, Global extensions of spacetimes describing asymptotic final states of black holes, Classical and Quantum Gravity 13, 539 (1996).
  * Rahman and Wald [2020] A. A. Rahman and R. M. Wald, Black hole memory, Phys. Rev. D 101, 124010 (2020).
  * Hollands  _et al._ [2022] S. Hollands, Á. D. Kovács, and H. S. Reall, The second law of black hole mechanics in effective field theory, Journal of High Energy Physics 2022, 258 (2022).
  * [18] A. C. Wall and Z. Yan, Linearised second law for higher curvature gravity and non-minimally coupled vector fields, to appear.
  * Chen  _et al._ [2022] P.-N. Chen, D. E. Paraizo, R. M. Wald, M.-T. Wang, Y.-K. Wang, and S.-T. Yau, Cross-section continuity of definitions of angular momentum, Classical and Quantum Gravity 40, 025007 (2022).
  * Wald [1993] R. M. Wald, Black hole entropy is the noether charge, Phys. Rev. D 48, R3427 (1993).
  * Davies and Reall [2023a] I. Davies and H. S. Reall, A non-perturbative second law of black hole mechanics in effective field theory, (2023a), arXiv:2312.07659 [hep-th] .
  * Davies and Reall [2023b] I. Davies and H. S. Reall, Dynamical Black Hole Entropy in Effective Field Theory, JHEP 05, 006, arXiv:2212.09777 [hep-th] .
  * Wald [1979] R. M. Wald, Entropy and black-hole thermodynamics, Phys. Rev. D 20, 1271 (1979).
  * Chandrasekaran  _et al._ [2018] V. Chandrasekaran, É. É. Flanagan, and K. Prabhu, Symmetries and charges of general relativity at null boundaries, Journal of High Energy Physics 2018, 125 (2018).
  * Rignon-Bret [2023] A. Rignon-Bret, Second law from the noether current on null hypersurfaces, Phys. Rev. D 108, 044069 (2023).
  * Ciambelli  _et al._ [2023] L. Ciambelli, L. Freidel, and R. G. Leigh, Null Raychaudhuri: Canonical Structure and the Dressing Time, (2023), arXiv:2309.03932 [hep-th] .
  * Azeyanagi  _et al._ [2009] T. Azeyanagi, G. Compère, N. Ogawa, Y. Tachikawa, and S. Terashima, Higher-Derivative Corrections to the Asymptotic Virasoro Symmetry of 4d Extremal Black Holes, Progress of Theoretical Physics 122, 355 (2009), https://academic.oup.com/ptp/article-pdf/122/2/355/9681053/122-2-355.pdf .
  * Wald [1990] R. M. Wald, On identically closed forms locally constructed from a field, Journal of Mathematical Physics 31, 2378 (1990).
  * Bueno  _et al._ [2017] P. Bueno, P. A. Cano, V. S. Min, and M. R. Visser, Aspects of general higher-order gravities, Phys. Rev. D 95, 044010 (2017).
  * Dong and Lewkowycz [2018] X. Dong and A. Lewkowycz, Entropy, extremality, euclidean variations, and the equations of motion, Journal of High Energy Physics 2018, 81 (2018).
  * Bhattacharyya  _et al._ [2021] S. Bhattacharyya, P. Dhivakar, A. Dinda, N. Kundu, M. Patra, and S. Roy, An entropy current and the second law in higher derivative theories of gravity, JHEP 09, 169, arXiv:2105.06455 [hep-th] .
  * Hollands  _et al._ [2007] S. Hollands, A. Ishibashi, and R. M. Wald, A Higher dimensional stationary rotating black hole must be axisymmetric, Commun. Math. Phys. 271, 699 (2007), arXiv:gr-qc/0605106 .
  * Unruh and Wald [1982] W. G. Unruh and R. M. Wald, Acceleration radiation and the generalized second law of thermodynamics, Phys. Rev. D 25, 942 (1982).
  * Wall [2012] A. C. Wall, Proof of the generalized second law for rapidly changing fields and arbitrary horizon slices, Phys. Rev. D 85, 104049 (2012).
  * Bousso  _et al._ [2016] R. Bousso, Z. Fisher, S. Leichenauer, and A. C. Wall, Quantum focusing conjecture, Phys. Rev. D 93, 064044 (2016).
  * Ceyhan and Faulkner [2020] F. Ceyhan and T. Faulkner, Recovering the qnec from the anec, Communications in Mathematical Physics 377, 999 (2020).
  * Casini and Huerta [2007] H. Casini and M. Huerta, A c-theorem for the entanglement entropy, J. Phys. A 40, 7031 (2007), arXiv:cond-mat/0610375 .
  * Andersson  _et al._ [2005] L. Andersson, M. Mars, and W. Simon, Local existence of dynamical and trapping horizons, Phys. Rev. Lett. 95, 111102 (2005).
  * Engelhardt and Wall [2019] N. Engelhardt and A. C. Wall, Coarse graining holographic black holes, Journal of High Energy Physics 2019, 160 (2019).
  * Müller  _et al._ [1997] U. Müller, C. Schubert, and A. van de Ven, A closed formula for the riemann normal coordinate expansion, arXiv preprint gr-qc/9712092 (1997).


# PARALLELPROJ - An open-source framework for fast calculation of projections
in tomography

Georg Schramm  Radiological Sciences Lab, Stanford School of Medince, Stanford
University, Stanford, US  Department of Imaging and Pathology, Division of
Nuclear Medicine, KU Leuven, Belgium

> Abstract: In this article, we present a new open source framework, called
> parallelproj, for fast parallel calculation of projections in tomography
> using multiple CPUs or GPUs. This framework implements forward and back
> projection functions in sinogram and listmode using Joseph’s method, which
> is also extended for time-of-flight PET projections. In a series of tests
> related to PET image reconstruction using data from a state-of-the-art
> clinical PET/CT system, we benchmark the performance of the projectors in
> non-TOF and TOF, sinogram and listmode using a multi CPU, hybrid CPU/GPU and
> pure GPU mode. We find that the GPU mode offers acceleration factors between
> 20 and 60 compared to the multi CPU mode and that OSEM listmode
> reconstruction of real world PET data sets is possible within a couple of
> seconds using a single state-of-the-art GPU.

##  1 Introduction

For tomographic imaging techniques used in medicine such as X-ray computed
tomography (CT), positron emission tomography (PET) and single photon emission
tomography (SPECT), image reconstruction results are usually expected within
seconds or minutes after data acquisition, creating a severe computational
challenge when reconstruction data from state-of-the-art systems. With new
scanner generations, the problem size of this challenge is steadily growing,
since (i) the data size is increasing due to higher resolution detectors and
scanners with bigger field of view [1], and (ii) more advanced (iterative)
reconstruction algorithms are being used that try to exploit more information
from the acquired data, which usually necessitates the calculation of a huge
amount of projections. An example of the latter is the data-driven motion
correction in PET [2] where instead of reconstructing a single “static frame
”, many very short time frames are reconstructed and subsequently used for
motion estimation and correction. Another example for (ii) is the combination
of deep learning and tomographic image reconstruction [3, 4, 5], using, e.g.
unrolled networks, where during training also a tremendous number of
projections have to be calculated to evaluate the gradient of the data
fidelity term across a mini batch in every training epoch.

For most tomographic image reconstruction algorithms, the bottleneck in terms
of computation time is the evaluation of a linear forward model that describes
the physics of the data acquisition process. In CT, PET, and SPECT, the
forward model includes the computation of many (weighted) line or volume
integrals through an image volume, commonly called “projections” - which can
be slow when executed on a single processor. Fortunately, for most
reconstruction algorithms, the computation of projections can be executed in
parallel on multiple processors, e.g. using multiple CPUs or one or more
graphics processing units (GPUs). Note that the evaluation of the adjoint of
the forward model - commonly called “back projection“‘ - is computationally
more demanding, since race conditions usually occur. In recent decades, the
use of GPUs for faster calculation of projections in tomographic imaging has
been studied extensively; see, e.g. [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
17, 18, 19, 20, 21] or the reviews [22, 23] for the use of GPUs in PET
reconstruction. All of these articles conclude that the time needed to
calculate forward and back projections on state-of-the-art GPUs is usually
much shorter compared to using multiple CPUs.

Motivated by these findings and the recent availability of very powerful low-
and high-level GPU programming frameworks such as CUDA and cupy [24], we
developed a new open source research framework, called parallelproj, for fast
calculations of projections in tomography.

The objectives of the parallelproj framework are as follows:

  * •

To provide an open source framework for fast parallel calculation of
projections in tomographic imaging using multiple CPUs or GPUs.

  * •

To provide an accessible framework that can be easily installed without the
need for compilation of source code on all major operating systems (Linux,
Windows, and macOS).

  * •

To provide a framework that can be efficiently used in conjunction with
pytorch [25] GPU arrays to facilitate research on tomographic imaging methods,
including deep learning.

This article on the parallelproj framework is organized as follows. We first
review Joseph’s method for calculating projections, followed by a short
overview of the design choices and implementation of parallelproj.
Subsequently, we report the results of a few benchmarks related to image
reconstruction in PET with and without time-of-flight (TOF) information using
sinograms or listmode (LM) before ending the article with a detailed
discussion and conclusion. In this article, we focus on the performance of
parallelproj projectors for non-TOF and TOF PET reconstruction problems. Note,
however, that the non-TOF Joseph projectors could also be used in iterative CT
or SPECT reconstruction.

image plane $k$image plane $k+1$image plane
$k+2$$(i,j)$$(i+1,j)$$(i,j+1)$$(i+1,j+1)$TOF bin $t$TOF bin $t+1$TOF bin
$t+n$ray with or withoutsubdivision into TOF bins Figure 1: Illustration of
Joseph’s method for projecting rays through voxel images. In a ray-driven
approach the image volume is traversed plane by plane along a principal
direction. At every plane, the intersection point between the ray and image
plane is calculated (orange dot). The contribution of the four nearest voxels
(gray dots) to the line integral is modeled using bi-linear interpolation
weights. The method can also be easily extended to compute TOF-weighted
projections using a subdivision of the ray into TOF bins and by evaluation of
a TOF kernel. See text and [28] for more details. Figure not drawn to scale.

##  2 Materials and Methods

###  2.1 Joseph’s method for projecting rays through voxel images

Besides Siddon’s method [26] and the distance-drive method [27], Joseph’s
method [28] is a very efficient and popular way to calculate projections in
transmission and emission tomography. The basic idea of the original Joseph
method for calculating line integrals, that can be used to model non-TOF
projections, is shown in Figure 1. For a given ray, the algorithm first
determines the principal direction that is most parallel to the ray and then
steps through the image volume plane by plane along this principal direction.
In every plane, the intersection point between the ray and the plane is
calculated and the contribution of the image at that point to the line
integral is approximated using bi-linear of the four nearest neighbors around
the intersection point. In other words, only the four nearest neighboring
voxels are contributing to the line integral and their contributions are given
by the bi-linear interpolation weights. Finally, the contributions of all
planes are added and corrected for the incidence angle of the ray - see [28]
for more details.

An extension of the original Joseph method to calculate TOF-weighted
projections is straight forward. For every voxel contributing to the line
integral and every TOF bin along the ray, a TOF weight can be computed by
evaluating a TOF kernel that is a function of the Euclidean distance between
the voxel and the center of the TOF bin. The TOF kernel can be, e.g. modeled
as a Gaussian kernel convolved with the width of the TOF bins, resulting in
the evaluation of the difference of two error functions.

###  2.2 Design principles and implementation details

The application programming interface to the parallelproj framework was
designed such that:

  * •

The input to the low level projector functions are as generic as possible. In
practice, that means that these functions take a list of coordinates
representing the start and end point of the rays to be projected as input,
making the low-level functions agnostic to specific scanner geometries (or
symmetries). Thus, any scanner geometry can be modeled.

  * •

Projections can be performed in non-TOF or TOF mode.

  * •

In the TOF mode, optimized projections for sinogram and listmode are
available. In the former, the contributions to all available TOF bins along a
ray are computed while traversing the image volume plane by plane, whereas in
the latter only the contribution to one specific TOF bin (the TOF bin of a
given listmode event) is evaluated.

  * •

The back projections are the exact adjoint of the forward projections (matched
forward and back projections).

Parallelization across multiple processors was implemented in two different
ways. To enable parallelization across multiple CPUs, a first version of the
parallelproj library was implemented using C and OpenMP [29]
(libparalleproj_c). Furthermore, the exact same projector functions were
implemented in CUDA to enable parallelization on one or multiple GPUs
(libparalleproj_cuda). In the CUDA version, the input data is first
transferred from the host to all available GPU(s) followed by the parallel
execution of the projection kernels. After the kernel execution, the result is
transferred back to the host. To handle race conditions, all implementations
use atomic add operations in the back projections.

###  2.3 Availability of source code and precompiled libraries

parallelproj is an open source project and the source code is available at
https://github.com/gschramm/parallelproj under MIT license. Next to the
sources, we also offer precompiled libraries for all major operating systems
(Linux, Windows, and macOS) and various recent CUDA versions using the conda-
forge package manager.111The CUDA version of the parallelproj library is not
available for macOS. In addition to the precompiled libraries, the
parallelproj package also includes the source file of the CUDA projection
kernels such that they can be directly executed on GPU arrays using frameworks
that allow for just-in-time compilation of CUDA kernels such as, e.g. cupy
[24].

###  2.4 parallelproj computation modes

Using the two aforementioned projection libraries, as well as the CUDA
projections kernels, projections can be performed in the following three
different computation modes:

  1. 1.

CPU mode: Forward and back projections of image volumes (arrays) stored on the
host (CPU) can be performed using libparallelproj_c where parallelization
across all available CPUs is performed using OpenMP.

  2. 2.

hybrid CPU/GPU mode: Forward and back projections of image volumes (arrays)
stored on the host canbe performed using libparallelproj_cuda involving data
transfer from the host to all available GPUs, execution of projection kernels
on the GPUs, and transfer of the results back to the host.

  3. 3.

direct GPU mode: Forward and back projections of image volumes (arrays) stored
on a GPU can be performed by direct execution of the projection kernels using
a framework that supports just-in-time compilation of CUDA kernels, such as
cupy [24]. In contrast to the hybrid CPU/GPU mode, memory transfer between
host and GPU is avoided.

Figure 2: Results of the non-TOF sinogram benchmark tests. The non-TOF subset
sinogram contained 415 radial elements, 8 views and 1292 planes (1 out of 34
subsets). The image used in these tests contained (215,215,71) voxels with an
isotropic voxel size of 2.78 mm. The mean and the standard deviation estimated
from 10 runs are represented by the colored bars and the black error bars,
respectively. Note the different limits on the y axes. Figure 3: Results of
the TOF sinogram benchmark tests. The TOF subset sinogram contained 415 radial
elements, 8 views, 1292 planes and 29 TOF bins with a width of 169 ps. The
modeled TOF resolution was 375 ps. The image used in these tests contained
(215,215,71) voxels with an isotropic voxel size of 2.78 mm. The mean and the
standard deviation estimated from 10 runs are represented by the colored bars
and the black error bars, respectively. Note the different limits on the y
axes.

###  2.5 Benchmark tests

To evaluate the performance of the parallelproj projectors using the three
computation modes described above, we implemented a series of benchmarking
tests. All tests are related to a PET image reconstruction task and used the
geometry and properties of a state-of-the-art GE Discovery MI TOF PET/CT
scanner [30] with 20 cm axial FOV. This scanner consists of 36 detector
“rings”, where each “ring” has a radius of 380 mm and consists of 34 modules
containing 16 detectors each such that there are 16$\times$34$\times$36 =
19584 detectors in total. A non-TOF emission sinogram for this scanner in span
1 has 415 radial elements, 272 views, and 1292 planes, meaning that for a full
non-TOF sinogram projection, 415$\times$272$\times$1292 = 146e6 line integrals
have to be evaluated. For TOF data, each line of response (LOR) is subdivided
into 29 TOF bins using a TOF bin width of 169 ps. The reported TOF resolution
of the scanner is 375 ps FWHM [30]. In the TOF projectors of parallelproj, the
Gaussian TOF kernel is truncated beyond $\pm$3 standard deviations.

To evaluate the performance of parallelproj for projections in sinogram mode,
we measured the time needed for a forward and back projection of a span 1
subset sinogram containing 8 equally spaced views in non-TOF and TOF mode.
This is equivalent to the projection work required for an OSEM subset update
using 34 subsets in total, a setting that is used in many clinical
reconstructions. Since it is known that the in-memory data order severely
affects the computation time, especially on CUDA devices, we varied the order
of the spatial axis of the sinogram, as well as the orientation of the
symmetry axis of the scanner. E.g. in the sinogram order mode “PVR”, the
radial direction increased the fastest and the plane direction increased the
slowest in memory. For the sinogram order mode “VRP”, the plane direction
increased the fastest and the view direction the slowest. By varying the
orientation of the scanner symmetry axis, we could test the impact of
different image volume memory layouts. For example, the symmetry axis mode “2”
meant that the image volume memory increased the fastest in the axial
direction of the scanner, while “0” or “1” meant that one of the transaxial
directions increased the fastest.

In addition to the sinogram projection tests, we also evaluated the
performance of parallelproj for non-TOF and TOF projections in listmode as a
function of the number of acquired listmode events. Instead of randomly
generating the event coordinates, listmode events from an acquisition of a
NEMA image quality phantom were used, which guaranteed a more realistic event
distribution. In contrast to the projections in sinogram mode, where the ray
directions and memory access are somehow ordered, they are random for unsorted
listmode data. Similarly to the sinogram tests, the symmetry axis of the
scanner was also varied. For all sinogram and listmode projection benchmarks,
the coordinates of all LOR start and endpoints were precalculated such that
the overhead of calculating the LOR coordinates was not included in these
tests. All benchmarks were repeated 10 times and the mean and standard
deviation of the results were calculated and visualized.

Finally, we also measured the time needed for a complete listmode OSEM
iteration using 34 subsets as a function of the number of listmode events in
the NEMA acquisition. In contrast to the projection benchmarks, the LOR start
and end points were not precalculated. Instead, they were generated by
indexing a lookup table stored in CPU or GPU memory containing all detector
coordinates using the transaxial and axial detector numbers provided in the
original listmode file.

All tests used an image of size (215,215,71), an isotropic voxel size of 2.78
mm, and were performed on a workstation including an AMD Ryzen Threadripper
PRO 3955WX 16 core 32 thread CPU with 256 GB RAM, and an NVIDIA GeForce RTX
3090 GPU with 24 GB RAM on Ubuntu 20.04 LTS using CUDA v11.2 and parallelproj
v1.2.9. The scanner geometry, as well as the list mode OSEM algorithm, was
implemented in a minimal proof-of-concept Python package available at
https://github.com/gschramm/pyparallelproj. For the CPU and hybrid CPU/GPU
mode, Python’s ctypes module is used to project numpy arrays stored in CPU
(host) memory using a minimal interface to the low level projection functions
defined in libparallelproj_c and libparallelproj_cuda. In GPU mode, the CUDA
projection kernels were just in time compiled and directly executed on cupy
GPU arrays. Due to the interoperability between numpy and cupy the same high-
level listmode OSEM implementation could be used for both modes. Note that in
the latter, all operations needed for the OSEM update were executed directly
on the cupy GPU arrays, eliminating any memory transfer between the host and
the GPU during OSEM updates. In all listmode OSEM reconstructions, a shift-
invariant image-based resolution model was used, including a Gaussian kernel
of 4.5 mm FWHM.

##  3 Results

Figure 4: Results of the non-TOF listmode benchmark tests for different number
of listmode events. The image used in these tests contained (215,215,71)
voxels with an isotropic voxel size of 2.78 mm. The mean and the standard
deviation estimated from 10 runs are represented by the colored bars and the
black error bars, respectively. Note the different limits on the y axes and
that the x-axis scale is non-linear. Figure 5: Results of the TOF listmode
benchmark tests for different number of listmode events. The image used in
these tests contained (215,215,71) voxels with an isotropic voxel size of 2.78
mm. The mean and the standard deviation estimated from 10 runs are represented
by the colored bars and the black error bars, respectively. Note the different
limits on the y axes and that the x-axis scale is non-linear.

Figure 6: (top) Results for the timing of a complete LM OSEM iteration
including 34 subset updates for the NEMA image quality phantom acquisition.
The image used in these tests contained (215,215,71) voxels with an isotropic
voxel size of 2.78 mm. The mean and the standard deviation estimated from 6
iterations are represented by the colored bars and the black error bars,
respectively. Note the different limits on the y axes and that the x-axis
scale is non-linear. (bottom) Transaxial and coronal slice of a listmode OSEM
reconstruction of the NEMA image quality phantom with 40e6 events after 6
iterations with 34 subsets using a standard Gaussian post filter of 4 mm FWHM.
Note that for better visibility, the reconstructed image was cropped to the
center portion of the transaxial FOV.

Figures 2 and 3 show the results of the sinogram benchmarks in non-TOF and TOF
mode, respectively. In non-TOF mode, the best results in terms of the summed
time needed for the forward and back projection of one subset sinogram were
(compute mode, sinogram order mode, scanner symmetry axis): 1.50 s for (CPU,
PVR, 0), 0.053 s for (hybrid CPU / GPU, RVP, 2) and 0.025 s for (GPU, RVP, 2),
meaning that the pure GPU mode was approximately 60x faster than the CPU mode
and 2.1x faster than the hybrid mode. In TOF mode, the corresponding results
were: 8.99 s for (CPU, PVR, 0), 0.277 s for (hybrid CPU/GPU, RVP, 2), and
0.189 s for (GPU, RVP, 2), which means that the pure GPU mode was
approximately 48x faster than the CPU mode and 1.5x faster than the hybrid
mode. As expected, especially for the back projections where atomic operations
are used, the memory order in the sinogram as well as in the image has a
substantial impact on the results. The ratios between the best and worst
results for the combined projection times in terms of sinogram order and
symmetry axis (non-TOF, TOF mode) were: (2.5, 1.6) in the CPU mode, (4.2, 3.0)
in the hybrid mode, and (8.0, 4.0) in the GPU mode.

Figures 4 and 5 show the results of the listmode benchmarks in non-TOF and TOF
mode, respectively. For 40e6 events, the best results in terms of time needed
for the forward and back projection (non-TOF, TOF) were: (34.8, s, 10.6, s) in
CPU mode, (4.09 s, 0.75 s) in hybrid mode and (3.91 s, 0.58 s) in GPU mode.
For 1.25e6 events, the best results in terms of time needed for the forward
and back projection (non-TOF, TOF) were: (1.08 s, 0.34 s) in CPU mode, (0.11
s, 0.029 s) in hybrid mode, and (0.1 s, 0.018 s) in GPU mode. In non-TOF mode,
the pure GPU mode was approximately 10x faster than the CPU mode and 1.1x
faster than the hybrid mode. In TOF mode, the pure GPU mode was approximately
20x faster than the CPU mode and 1.6x faster than the hybrid mode. In contrast
to the sinogram benchmark results, the impact of the scanner symmetry axis
direction was small. In the CPU and GPU mode, the increase in projection time
as a function of the number of list-mode events was almost perfectly linear.
In hybrid mode at low number of events, the scaling was not linear due to the
overhead caused by the time needed for memory transfer.

Figure 6 shows the results for the timing of a complete TOF listmode OSEM
iteration, including 34 subset updates, as well as a reconstruction of the
NEMA image quality phantom data set using 40e6 total prompt events. The best
results for (40e6, 1.25e6) events were: (18.47 s, 3.77 s) in CPU mode, (8.48
s, 3.39 s) in hybrid mode, and (0.63 s, 0.065 s) in GPU mode which means that
for 40e6 events the pure GPU mode was approximately 29x faster than the CPU
mode and 13.4x faster than the hybrid mode.

##  4 Discussion

All results shown in our article demonstrate once more that parallel
computation of forward and back projections using a state-of-the-art GPU is
substantially faster compared to parallelization using OpenMP on a state-of-
the-art multicore CPU system. Certainly, the achievable GPU acceleration
factor strongly depends on the computational problem itself (e.g. sinogram or
listmode reconstruction) and the problem size. In our non-TOF and TOF sinogram
and listmode benchmark tests, we observed GPU acceleration factors between 10
and 60.

One important aspect that emerged from our sinogram benchmark tests - where
the projection data and memory access is ordered - is the fact that the
projection times varied substantially when using different memory layouts (up
to a factor of 8 in the GPU mode). This can be understood by taking into
account that the amount of race conditions that are created during the back
projection within a thread block heavily depends on the order and possible
intersections of rays to be projected within that block. Note that in pure GPU
mode, the time needed for sinogram forward projections also varied
substantially across the different memory layouts, which is probably due to
the way image memory is accessed and cached on CUDA GPUs.

Another interesting observation is the fact that in all compute modes the time
needed to calculate TOF sinogram projections was much longer than the times
needed to calculate non-TOF sinogram projections, whereas the situation was
reversed in listmode. For TOF sinogram projections, more floating point
operations have to be computed compared to non-TOF sinogram projections due to
the evaluations of the TOF kernels between the contributing voxels and a
number of TOF bins. In listmode, however, the computational work needed to
project a TOF event is much lower compared to projecting a non-TOF event. This
is the case because a TOF listmode event detected in a specific TOF bin is
only affected by a few voxels along the complete LOR in the image, where the
number of affected voxels is inversely proportional to the TOF resolution of
the scanner. That in turn means that with scanner TOF resolutions becoming
better and better, the gap between the TOF projection times in sinogram and
listmode will become bigger and bigger, strongly favoring listmode processing.
According to our experience, projection times in listmode are already much
faster for most standard clinical acquisitions (except for very long static
brain scans with high affinity tracers) on current PET systems with TOF
resolutions between 250-400 ps.222An alternative way to further accelerate
sinogram-based reconstructions is the use of dedicated sinogram rebinning
techniques. Extrapolating the timing results of one complete OSEM listmode
iteration of an acquisition with 40e6 counts in Figure 6, clinical listmode
OSEM reconstructions of a single bed position of a standard static FDG whole-
body acquisition using PET scanner with 20-25 cm axial FOV seem to be possible
in a couple of seconds and could even be faster than the acquisition
time.333This is obviously only true if all other necessary corrections, such
as scatter estimation, can be performed very quickly as well. For PET scanners
with a very long axial field of view (much higher sensitivity), the
reconstruction times could be substantially longer.

A somewhat unexpected result was the fact that the gap in the TOF projection
times between hybrid CPU/GPU and pure GPU mode was much bigger when timing the
execution of a complete listmode OSEM iteration compared to the pure
projection benchmark test when reconstruction 40e6 counts (approximately a
factor of 13 vs a factor of 1.6, respectively). After detailed profiling of a
listmode OSEM iteration in hybrid mode, it became obvious that the total time
spent for the 34 subset listmode forward and back projections (ca. 1.2 s) was
short compared to the time needed to calculate all other operations necessary
for the OSEM update. Profiling revealed that calculating the LOR start and end
coordinates using the lookup table approach took approximately 3.5 s and that
calculating all 68 Gaussian convolutions needed for image-based resolution
modeling took approximately 2.3 s. Keep in mind that those operations were
performed on the CPU in hybrid compute mode. While the overhead due to the
lookup of the LOR start and end coordinates can be completely eliminated by
precalculating those coordinates, reducing the overhead caused by the CPU-
based Gaussian convolution is more complicated, but probably possible as well,
which would bring the reconstruction time gap closer to the gap observed in
the pure projection benchmarks. An interesting lesson to be learned is that
once very fast GPU-based projectors are used, it should always be double-
checked whether other computational steps of any algorithm become new
bottlenecks.

An important limitation of our study is the fact that we only implemented and
benchmarked Joseph’s projection method. Compared to other methods such as the
distance-driven method, multiray models, or tube-of-response models, Joseph’s
method is inherently faster. Consequently, projection times are expected to be
somewhat slower for more advanced projectors, but a detailed investigation of
more advanced projectors is beyond the scope of this work and left for future
research. 444Since parallelproj is an open-source project, contributions of or
discussions on more advanced projectors from the reconstruction community are
more than welcome. Note, however, that according to our experience, combining
Joseph’s method with an image-based and / or sinogram-based resolution model
produces high-quality PET reconstructions.

Without a doubt, it is possible to further optimize the implementation of the
parallelproj projectors, especially the CUDA implementation. As an example, we
have decided not to use CUDA’s texture memory, which could substantially
accelerate the image interpolations needed in the Joseph forward projections,
or be also used to interpolate TOF kernel values based on a 1D lookup table.
The main reason for not using texture memory is the fact that it would only
accelerate the forward projections since writing into texture memory is not
possible and because reconstruction times are usually dominated by the back
projections. Another way to further improve the listmode projection times is
to pre-sort the listmode events to minimize race conditions during back
projection, as e.g. shown in [9].

##  5 Conclusion

parallelproj is an open-source, and easy accessible research framework for
efficient calculation of non-TOF and TOF projections in sinogram or listmode
on multiple CPUs or state-of-the-art CUDA GPUs. Conventional and advanced
research reconstructions (including deep learning) can be substantially
accelerated by using the hybrid and pure GPU compute modes of this framework.

## References

  * [1] Xuezhu Zhang et al.  “Quantitative image reconstruction for total-body PET imaging using the 2-meter long EXPLORER scanner”  In _Physics in Medicine and Biology_ 62.6 IOP Publishing, 2017, pp. 2465  DOI: 10.1088/1361-6560/aa5e46
  * [2] F Lamare et al.  “PET respiratory motion correction: quo vadis?”  In _Physics in Medicine and Biology_ 67.3 IOP Publishing, 2022, pp. 03TR02  DOI: 10.1088/1361-6560/ac43fc
  * [3] Ge Wang, Jong Chu Ye, Klaus Mueller and Jeffrey A Fessler  “Image Reconstruction Is a New Frontier of Machine Learning”  In _IEEE Transactions on Medical Imaging_ 37.2, 2018, pp. 1–14  DOI: 10.1109/TMI.2018.2833635
  * [4] Ge Wang, Jong Chul Ye and Bruno De Man  “Deep learning for tomographic image reconstruction” Publisher: Springer US  In _Nature Machine Intelligence_ 2.12, 2020, pp. 737–748  DOI: 10.1038/s42256-020-00273-z
  * [5] Andrew J. Reader et al.  “Deep Learning for PET Image Reconstruction”  In _IEEE Transactions on Radiation and Plasma Medical Sciences_ 5.1, 2021, pp. 1–25  DOI: 10.1109/TRPMS.2020.3014786
  * [6] Guillem Pratx et al.  “Fully 3-D List-Mode OSEM Accelerated by Graphics Processing Units” ISSN: 1082-3654  In _2006 IEEE Nuclear Science Symposium Conference Record_ 4, 2006, pp. 2196–2202  DOI: 10.1109/NSSMIC.2006.354350
  * [7] Guillem Pratx, Garry Chinn, Peter D. Olcott and Craig S. Levin  “Fast, Accurate and Shift-Varying Line Projections for Iterative Reconstruction Using the GPU” Conference Name: IEEE Transactions on Medical Imaging  In _IEEE Transactions on Medical Imaging_ 28.3, 2009, pp. 435–445  DOI: 10.1109/TMI.2008.2006518
  * [8] W. Barker, Shanthalaxmi Thada and William Dieckmann  “A GPU-accelerated implementation of the MOLAR PET reconstruction package” ISSN: 1082-3654  In _2009 IEEE Nuclear Science Symposium Conference Record (NSS/MIC)_ , 2009, pp. 4114–4119  DOI: 10.1109/NSSMIC.2009.5402353
  * [9] Jing-yu Cui, Guillem Pratx, Sven Prevrhal and Craig S. Levin  “Fully 3D list-mode time-of-flight PET image reconstruction on GPUs using CUDA” _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1118/1.3661998  In _Medical Physics_ 38.12, 2011, pp. 6775–6786  DOI: 10.1118/1.3661998
  * [10] J.. Herraiz et al.  “Fully 3D GPU PET reconstruction”  In _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 648, NIMA_4th International Conference on Imaging techniques in Subatomic Physics, Astrophysics, Medicine, Biology and Industry, 2011, pp. S169–S171  DOI: 10.1016/j.nima.2010.12.043
  * [11] J.. Herraiz et al.  “GPU-Based Fast Iterative Reconstruction of Fully 3-D PET Sinograms” Conference Name: IEEE Transactions on Nuclear Science  In _IEEE Transactions on Nuclear Science_ 58.5, 2011, pp. 2257–2263  DOI: 10.1109/TNS.2011.2158113
  * [12] Kyung Sang Kim and Jong Chul Ye  “Fully 3D iterative scatter-corrected OSEM for HRRT PET using a GPU”  In _Physics in Medicine and Biology_ 56.15, 2011, pp. 4991  DOI: 10.1088/0031-9155/56/15/021
  * [13] Jian Zhou and Jinyi Qi  “Fast and efficient fully 3D PET image reconstruction using sparse system matrix factorization with GPU acceleration”  In _Physics in Medicine and Biology_ 56.20, 2011, pp. 6739  DOI: 10.1088/0031-9155/56/20/015
  * [14] Cheng-Ying Chou et al.  “Accelerating Image Reconstruction in Dual-Head PET System by GPU and Symmetry Properties” Publisher: Public Library of Science  In _PLOS ONE_ 7.12, 2012, pp. e50540  DOI: 10.1371/journal.pone.0050540
  * [15] Shoko Kinouchi et al.  “GPU-Based PET Image Reconstruction Using an Accurate Geometrical System Model” Conference Name: IEEE Transactions on Nuclear Science  In _IEEE Transactions on Nuclear Science_ 59.5, 2012, pp. 1977–1983  DOI: 10.1109/TNS.2012.2201953
  * [16] Jingyu Cui, Guillem Pratx, Bowen Meng and Craig S. Levin  “Distributed MLEM: An Iterative Tomographic Image Reconstruction Algorithm for Distributed Memory Architectures” Conference Name: IEEE Transactions on Medical Imaging  In _IEEE Transactions on Medical Imaging_ 32.5, 2013, pp. 957–967  DOI: 10.1109/TMI.2013.2252913
  * [17] S. Ha, S. Matej, M. Ispiryan and K. Mueller  “GPU-Accelerated Forward and Back-Projections With Spatially Varying Kernels for 3D DIRECT TOF PET Reconstruction” Conference Name: IEEE Transactions on Nuclear Science  In _IEEE Transactions on Nuclear Science_ 60.1, 2013, pp. 166–173  DOI: 10.1109/TNS.2012.2233754
  * [18] P.. Markiewicz et al.  “High throughput CUDA implementation of accurate geometric modelling for iterative reconstruction of PET data”  In _2014 IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC)_ , 2014, pp. 1–4  DOI: 10.1109/NSSMIC.2014.7430963
  * [19] Jian Zhou and Jinyi Qi  “Efficient fully 3D list-mode TOF PET image reconstruction using a factorized system matrix with an image domain resolution model” Publisher: IOP Publishing  In _Physics in Medicine and Biology_ 59.3, 2014, pp. 541  DOI: 10.1088/0031-9155/59/3/541
  * [20] Moulay Ali Nassiri, Jean-François Carrier and Philippe Després  “Fast GPU-based computation of spatial multigrid multiframe LMEM for PET”  In _Medical and Biological Engineering and Computing_ 53.9, 2015, pp. 791–803  DOI: 10.1007/s11517-015-1284-9
  * [21] Tianyi Zeng et al.  “A GPU-accelerated fully 3D OSEM image reconstruction for a high-resolution small animal PET scanner using dual-ended readout detectors” Publisher: IOP Publishing  In _Physics in Medicine and Biology_ 65.24, 2020, pp. 245007  DOI: 10.1088/1361-6560/aba6f9
  * [22] Anders Eklund, Paul Dufort, Daniel Forsberg and Stephen M. LaConte  “Medical image processing on the GPU – Past, present and future”  In _Medical Image Analysis_ 17.8, 2013, pp. 1073–1094  DOI: 10.1016/j.media.2013.05.008
  * [23] Philippe Depres and Xun Jia  “A review of GPU-based medical image reconstruction | Elsevier Enhanced Reader”, 2017, pp. 76–92  DOI: 10.1016/j.ejmp.2017.07.024
  * [24] Ryosuke Okuta et al.  “CuPy: A NumPy-Compatible Library for NVIDIA GPU Calculations”  In _Proceedings of Workshop on Machine Learning Systems (LearningSys) in The Thirty-first Annual Conference on Neural Information Processing Systems (NIPS)_ , 2017  URL: http://learningsys.org/nips17/assets/papers/paper_16.pdf
  * [25] Adam Paszke et al.  “PyTorch: An Imperative Style, High-Performance Deep Learning Library”  In _Advances in Neural Information Processing Systems_ 32 Curran Associates, Inc., 2019  URL: https://proceedings.neurips.cc/paper/2019/file/bdbca288fee7f92f2bfa9f7012727740-Paper.pdf
  * [26] Robert L. Siddon  “Fast calculation of the exact radiological path for a three-dimensional CT array” _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1118/1.595715  In _Medical Physics_ 12.2, 1985, pp. 252–255  DOI: 10.1118/1.595715
  * [27] Bruno De Man and Samit Basu  “Distance-driven projection and backprojection in three dimensions.” ISBN: 0-7803-7636-6  In _Physics in medicine and biology_ 49.11, 2004, pp. 2463–2475  DOI: 10.1109/NSSMIC.2002.1239600
  * [28] Peter M. Joseph  “An Improved Algorithm for Reprojecting Rays Through Pixel Images”  In _IEEE Transactions on Medical Imaging_ 1.3, 1982, pp. 192–196  DOI: 10.1109/TMI.1982.4307572
  * [29] L. Dagum and R. Menon  “OpenMP: an industry standard API for shared-memory programming”  In _IEEE Computational Science and Engineering_ 5.1, 1998, pp. 46–55  DOI: 10.1109/99.660313
  * [30] David F.C. Hsu et al.  “Studies of a next-generation silicon-photomultiplier-based time-of-flight PET/CT system”  In _Journal of Nuclear Medicine_ 58.9, 2017, pp. 1511–1518  DOI: 10.2967/jnumed.117.189514