\section{Introduction}

The stress-energy tensor of a quantized field is an extremely useful tool for studying quantum effects in curved space because it takes both particle
production and vacuum polarization into account.  It can be computed in a background spacetime to obtain the energy density, pressure, etc.
for a quantum field in that spacetime.  It can also be used in the context of semiclassical gravity to compute the backreaction of the quantum field on the
spacetime geometry.

For black holes in four-dimensional, 4D, spacetimes, the full stress-energy tensor must be computed numerically.  This is a difficult task that has to date only been done without other approximations for the cases of static spherically symmetric black holes~\cite{fawcett,howard-candelas,howard,jensen-ottewill,jmo,jensen-et-al,ahs1,ahs2,ahl,choag,abf,breen-ottewill,levi-ori,levi,Zilberman-Levi-Ori} and the stationary Kerr metric~\cite{duffy-ottewill,levi-et-al-kerr}.
However, because of the difficulty involved, to our knowledge, no one has numerically computed the full stress-energy tensor for a quantized field in a 4D spacetime in which a black hole forms from collapse.  This is important because there can be a significant difference between the stress-energy tensor for a quantum field in a 2D versus a 4D spacetime such as that found for a massless minimally coupled scalar field in an extreme Reissner-Nordstrom spacetime~\cite{trivedi,ahl}.

In this paper we present a  method to compute the renormalized stress-energy tensor, $\langle {\rm in}| T_{ab} | {\rm in} \rangle$, for a massless minimally coupled scalar field in the case that
a black hole forms from the collapse of a spherically symmetric null shell.  The method works in the region outside both the shell and the event horizon.

In the region outside the null shell, Birkhoff's theorem ensures that the metric is that for Schwarzschild spacetime~\eqref{metric-sch}.
In the region inside the shell the space is flat.  Thus in both regions the mode equation for the quantum field is separable and inside the shell its solutions are
known analytically.  This allows for a numerical computation of the stress-energy tensor for the field in which only ordinary differential equations need to be solved
numerically.

For the collapsing null shell model, the initial vacuum state of the quantum field is well defined and
the main complication that occurs is due to the propagation of the modes across the null shell surface.  The crux of our method involves the expansions of the {\it in}
modes in terms of a complete set of solutions to the mode equation in the region outside the shell.

The stress-energy tensor for the quantum field is obtained by expanding the quantum field in terms of a complete set of modes.
This expansion is substituted into the formula for the stress-energy tensor of the corresponding classical field and the expectation value is computed.
If the field is in the {\it in} vacuum state then the result is an expression which involves sums and integrals over the mode functions for the {\it in} state and their derivatives.  After the renormalization counterterms are subtracted off, the resulting stress-energy tensor is finite and can be computed. This
is straightforward inside the null shell since the mode functions are known analytically and for the {\it in} state, the result is that the stress-energy tensor is equal to zero.

Outside the null shell and the event horizon the {\it in} modes do not assume a simple form in 4D.  One way to extend them to this region would be to use their values on the null shell trajectory along with their values at past null infinity, $\mathscr{I}^{-}$, as initial data for direct numerical computations of these mode function using numerical methods for partial differential equations.  A second way, which is the one we adopt here, is expand each of the {\it in} modes in terms of a complete set of modes in Schwarzschild spacetime.  The complete set of modes can be obtained by solving ordinary differential equations since the mode equation is separable.  The radial parts of these modes and some of the matching parameters must be computed numerically.  The mode matching has been tested in the 2D case where there is no effective potential in the mode equation.  It has also been partially tested for spherically symmetric modes in 4D when the effective potential is modeled as a delta function and for the full effective potential.

When the expansions for the {\it in} modes are substituted into the formula for the unrenormalized stress-energy tensor one finds
a triple integral over the mode functions which form a complete set in Schwarzschild spacetime.  Renormalization can be accomplished by subtracting the corresponding
expression that occurs for the Unruh state, then adding that expression back and subtracting the renormalization counterterms.
  The result is the sum of two finite tensors, the difference between the stress-energy tensor for the {\it in} state and the Unruh state and the renormalized stress-energy tensor for the Unruh state.  The latter has been numerically computed for the masslesss minimally coupled scalar field in~\cite{levi-ori,levi}.
This type of renormalization scheme has been used to compute the stress-energy tensors in Schwarzschild spacetime in the Unruh state for the conformally coupled massless scalar field~\cite{elster,jmo} and for the massless spin $1$ field~\cite{jmo}.

We have tested our method by numerically computing the stress-energy tensor in the collapsing null shell spacetime in 2D where the answer can be compared with
previous analytic calculations for the Unruh state~\cite{Davies-Fulling-Unruh} and the {\it in} vacuum state for the collapsing null shell spacetime~\cite{hiscock, Fabbri:2005mw}.
Our results are in agreement with those calculations.

In Sec. II we introduce the collapsing null shell model and then discuss the modes for a massless minimally coupled scalar field in both the null shell spacetime and pure Schwarzschild spacetime.  A detailed description of the method, including the computation of the $in$ modes and the  renormalization of the stress-energy tensor is given in Sec. III.  General expressions for the matching coefficients in the 4D case are derived in Sec. IV followed by examples where the matching method is tested.
Formulas needed for the computation of the stress-energy tensor in the 4D case are derived in the first part of Sec. V.  In the second part the stress-energy tensor is computed using our method in the 2D case and compared with previous calculations.
Sec. VI contains a summary of our results.  The appendices contain some details of a proof and some derivations that are used in the 2D examples in Secs. IV and V.
Throughout the paper, we use the sign conventions of~\cite{Misner} and units are chosen such that $\hbar=c=G=1$.

\section{Collapsing Null Shell Model}

The type of quantum field we consider is a massless minimally coupled scalar field.  In a static spherically symmetric spacetime it can be expanded in terms of a complete set of modes $f_{\omega \ell m}$ such that
\begin{equation} \phi = \sum_{\ell = 0}^\infty  \sum_{m = -\ell}^\ell \int_0^\infty [a_{\omega \ell m} f_{\omega \ell m} + a^\dagger_{\omega \ell m} f^{*}_{\omega \ell m}] \;, \label{Phi-expansion} \end{equation}
with $a_{\omega \ell m}$ an annihilation operator.  The modes are solutions to the equation
\begin{equation} \Box f_{\omega \ell m} = 0  \;, \label{Box-f} \end{equation}
and have the general form
\begin{equation} f_{\omega \ell m} = N \frac{Y_{\ell, m}(\theta, \phi)}{r } \psi(\tau,r) \;. \label{f-def-1} \end{equation}
The normalization constant $N$ is fixed by the condition
\begin{equation}  (f_{\omega \ell m}, f_{\omega' \ell' m'}) = \delta_{\ell,\ell'} \delta_{m,m'} \delta(\omega-\omega') \;. \label{f-norm} \end{equation}
The scalar product is defined by
\begin{eqnarray}
  (f_1, f_2) = -i\int_{\Sigma} d \Sigma \, n^\mu [f_1(x)\overset{\leftrightarrow}{\partial_\mu}f_2^{*}(x)]
 \;, \label{scalar-products}
\end{eqnarray}
where $n^{\mu}$ is a future-directed unit vector orthogonal to the spacelike (or null) hypersurface $\Sigma$ and $d\Sigma$ is the volume element in $\Sigma$. The hypersurface $\Sigma$ is taken to be a Cauchy surface and it is assumed that the spacetime is globally hyperbolic.

We consider a model in which a spherically symmetric black hole forms from the collapse of a null shell.  The Penrose diagram is shown in Fig.~\ref{BH_4D_penrose}.  The metric inside the shell is the flat space metric
\begin{equation}
ds^2 = -dt^2 + dr^2 + r^2 d \Omega^2    \;, \label{metric-flat}
\end{equation}
and the metric outside the shell is the Schwarzschild metric
\begin{equation}
ds^2 = -\left(1-\frac{2M}{r} \right) dt_s^2
+ \left( 1-\frac{2M}{r} \right)^{-1} dr^2 + r^2 d \Omega^2  \;. \label{metric-sch}
\end{equation}
It is useful inside the shell to define the radial null coordinates
\begin{subequations} \begin{eqnarray}
u &=& t - r \;, \label{u-flat}\\
 v &=& t + r  \;,\, \label{v-flat} \end{eqnarray}
\label{u-v-flat}
\end{subequations}
and outside the shell to define the radial null coordinates
\begin{subequations} \begin{eqnarray}
u_s &=& t_s - r_* \;, \label{u-sch-def}\\
  v &= & t_s + r_* \;, \label{v-sch-def}
 \end{eqnarray}   \label{u-v-sch}
\end{subequations}
where
\begin{eqnarray}
  r_* &= & r + 2M \log \left(\frac{r-2M}{2M} \right)  \;.
\label{rstar-def}
\end{eqnarray}
The angular coordinates, the $r$ coordinate and the $v$ coordinate are continuous across the shell.  The time coordinate is discontinuous as are $u$ and $u_s$ .
The relationship between the latter two is~\cite{m-p,Fabbri:2005mw}
\begin{equation}
u_s = u - 4M \log \left( \frac{v_H-u}{4 M} \right)  \;,
\label{us-u}
\end{equation}
with
\begin{equation}
v_H \equiv v_0 - 4M  \;.
\label{vH-def}
\end{equation}
Note that the value of $u$ on the event horizon is $v_H$.  Inverting, one finds that~\cite{mirror-bh}
\begin{equation}
u = v_H - 4 M \, W\left[\exp\left(\frac{v_H - u_s}{4 M}\right) \right]  \;,
\label{u-us}
\end{equation}
with $W$ the Lambert W function.

\begin{figure}[h]
\centering
\includegraphics [trim=0cm 20cm 0cm 0cm,clip=true,totalheight=0.3\textheight] {BH_4D_penrose.pdf}
\caption{Penrose diagram for a spacetime in which a null shell collapses to form a spherically symmetric black hole.  The trajectory of the shell (dashed blue curve) is $v = v_0$.  The horizon, $H^{+}$, is the dotted red curve.}
\label{BH_4D_penrose}
\end{figure}

For all of the mode functions that we consider the normalization constant is
\begin{equation} N = \frac{1}{\sqrt{4 \pi \omega}} \;. \label{Normalization-constant} \end{equation}
Inside the shell, using~\eqref{f-def-1} and the metric~\eqref{metric-flat} in~\eqref{Box-f} one finds mode functions of the form
\begin{subequations} \begin{eqnarray} f_{\omega \ell m} &=&  \frac{Y_{\ell m}(\theta,\phi)}{r \sqrt{4 \pi \omega}} \psi_{\omega \ell}(t,r)\;, \label{f-form-in} \\  \nonumber \\
    \psi_{\omega \ell} &=& e^{-i \omega t} \chi_{\omega \ell}(r) \;,  \\ \nonumber \\
     \frac{d^2 \chi_{\omega \ell}}{d r^2} &=& - \left[\omega^2 - \frac{\ell(\ell+1)}{r^2} \right] \chi_{\omega \ell} \;. \label{mode-eq-inside}
\end{eqnarray} \end{subequations}
Outside the shell, using~\eqref{f-def-1}, the metric~\eqref{metric-sch}, and~\eqref{rstar-def} in~\eqref{Box-f} one finds mode functions of the form
\begin{subequations} \begin{eqnarray} f_{\omega \ell m} &=&  \frac{Y_{\ell m}(\theta,\phi)}{r \sqrt{4 \pi \omega} } \psi_{\omega \ell}(t_s,r)\;,  \label{f-form-out}  \\ \nonumber \\
    \psi_{\omega \ell} &=& e^{-i \omega t_s} \chi_{\omega \ell}(r) \;, \label{psi-outside} \\ \nonumber \\
     \frac{d^2 \chi_{\omega \ell}}{d r_{*}^2} &=& -  \left[\omega^2 - \left(1 - \frac{2M}{r} \right) \left( \frac{2 M}{r^3}+  \frac{\ell(\ell+1)}{r^2} \right) \right] \chi_{\omega \ell} \;.\label{chi-eq-outside}  \end{eqnarray} \label{modes-outside} \end{subequations}
Note that in pure Schwarzschild spacetime outside of the past and future horizons, the form of the modes and the equation they satisfy are also given by~\eqref{modes-outside}.

The {\it in} vacuum state in the null shell spacetime is specified by
\begin{equation} \psi^{\rm in}_{\omega \ell} = e^{-i \omega v} \; \label{in-vacuum}\end{equation}
on past null infinity $\mathscr{I^{-}}$ and by requiring that $\psi_{\omega \ell}$ vanish at $r = 0$ in the region inside the shell so that $f_{\omega \ell m}$ is regular there.
The solution that has these properties is
\begin{equation} \psi^{\rm in}_{\omega \ell} =  C_{\ell} \, e^{-i \omega t}\, \omega r  j_\ell(\omega r) \;, \label{psi-in-a} \end{equation}
where $C_{\ell}$ is a normalization constant and $j_\ell $ is a spherical Bessel function.  The condition~\eqref{in-vacuum} fixes the value of $C_{ \ell}$.  For example, for $\ell = 0$, it is easy to show that $C_{0} = -2 i $ and
\begin{equation} \psi^{\rm in}_{\omega 0} = e^{-i \omega v} - e^{-i \omega u} \;. \label{psi-in-0} \end{equation}

\section{Method to compute the stress-energy tensor}
\label{sec:method-Tab}

The stress-energy tensor for the quantized massless minimally coupled scalar field, $\langle T_{ab} \rangle$, is to be computed for the {\it in} vacuum state in the region outside the null shell and outside the event horizon. The stress-energy tensor for the classical field is
\begin{equation} T_{ab} = \partial_a \Phi  \partial_b \Phi - \frac{1}{2} g_{ab} g^{cd} \partial_c \Phi  \partial_d \Phi  \;. \label{Tab-class} \end{equation}
To compute $\langle {\rm in}| T_{ab} |{\rm in} \rangle$, one can substitute~\eqref{Phi-expansion} into~\eqref{Tab-class}, use the complete set of modes for the {\it in} vacuum state $f^{\rm in}_{\omega \ell m}$, and compute the expectation value.  There are two things which make this difficult. One is computing the modes $f^{\rm in}_{\omega \ell m}$ in the region outside the shell and the other is renormalizing the stress-energy tensor.  Our method to compute the stress-energy tensor provides one way to overcome these difficulties.

First, we renormalize by subtracting from the unrenormalized expression for the stress-energy tensor for the {\it in} vacuum state, the unrenormalized stress-energy tensor for the Unruh state.  Since the renormalization counterterms are local and thus do not depend on the state of the quantum field, this quantity will be finite.  Then we add back the
 unrenormalized stress-energy tensor for the Unruh state and then subtract from it the renormalization counter terms.
 Schematically one can write
 \begin{eqnarray} \langle {\rm in}| T_{ab}| {\rm in} \rangle_{\rm ren} &=& \Delta \langle T_{ab} \rangle +   \langle U| T_{ab} | U \rangle_{\rm ren} \;,  \nonumber \\
               \Delta \langle T_{ab} \rangle  &=&  \langle {\rm in}| T_{ab} |{\rm in} \rangle_{\rm unren} - \langle U| T_{ab} | U \rangle_{\rm unren} \;. \label{Tab-Unruh-sub} \end{eqnarray}
The quantity $\langle U| T_{ab} | U \rangle_{\rm ren}$ has been numerically computed for a massless minimally coupled scalar field in Schwarzschild spacetime~\cite{levi-ori, levi}.
Thus what remains is to compute the difference between the unrenormalized expressions.  To do that it is necessary to discuss the computation of the mode functions for the quantum field that are relevant for the {\it in} and Unruh states.
It is worth pointing out that the computation of $\langle U| T_{ab} | U \rangle_{\rm unren}$ done in~\cite{levi-ori,levi} was done for pure Schwarzschild spacetime outside the event horizon.  However the computation we wish to do for
$\langle {\rm in}| T_{ab}| {\rm in} \rangle_{\rm ren}$ is for the null shell spacetime outside both the shell and the horizon.  The reason that there is no problem is that the renormalization counterterms are local and so are the same in this part of the null shell spacetime as they are in pure Schwarzschild spacetime.



Analytic expressions for the mode functions in the {\it in} vacuum state, $f^{\rm in}_{\omega \ell m}$
inside the shell are given in~\eqref{psi-in-a}.  However, it is not easy to continue these to the region outside the shell because the time coordinate $t$ and the right moving radial null coordinate $u$ are not continuous across the shell.  As a result, the modes $f^{\rm in}_{\omega \ell m}$, while still being solutions to~\eqref{Box-f}, do not have the form~\eqref{f-form-out} outside the shell. However, the known solutions inside the null shell along with their behavior on $\mathscr{I}^{-}$ can be used to fix the initial data on a Cauchy surface in the null shell spacetime.  The Cauchy surface we consider here, consists of the part of $\mathscr{I}^{-}$ with $v_0 \le v < \infty$ along with the trajectory of the null shell.  This initial data could be used for a numerical calculation of the partial differential equation satisfied by $f^{\rm in}_{\omega \ell m}$ outside the shell.  Alternatively, one can expand  $f^{\rm in}_{\omega \ell m}$ in terms of a complete set of modes in the region outside the shell and use the data on the Cauchy surface to determine the matching coefficients.

Here we take a variation of the latter approach by noting that the spacetime geometry outside the shell is the Schwarzschild geometry.  Because of this, it is possible to do the matching in the corresponding part of Schwarzschild spacetime.  The advantage of this is that the matching can be to a complete set of modes in the region outside the horizon in Schwarzschild spacetime.  These modes are well understood and straight-forward to work with numerically.  The disadvantage is that the relevant part of the Cauchy surface in the null shell spacetime discussed above does not form a Cauchy surface in the Schwarzschild spacetime.  This can be remedied by adding a segment along the future horizon with $-\infty < v \le v_0$.  The result is a Cauchy surface for the part of Schwarzschild spacetime that is outside of the past and future horizons.  It is illustrated in Fig.~\ref{fig:Cauchy}.  It is worth noting that the part of the Cauchy surface on the future horizon is not causally connected with the region outside the future horizon and outside the surface $v = v_0$.  The corresponding region in the null shell spacetime is the region where we want to compute the stress-energy tensor.  Thus any initial data
can be used for the mode function $f^{\rm in}_{\omega \ell m}$ on that surface so long as $f^{\rm in}_{\omega \ell m}$ is continuous at the point where the future horizon intersects the part of the Cauchy surface with $v = v_0$.

\begin{figure}[h]
\centering
\includegraphics   [trim=0cm 0cm 0cm 0cm,clip=true,totalheight=0.25\textheight]{Full-Schwarzschild-Penrose-diagram.pdf}
\caption{Penrose diagram for Schwarzschild spacetime showing the Cauchy surface used for matching the {\it in} modes in the null shell spacetime
to a complete set of modes in Schwarzschild spacetime in the region outside the past and future horizons.  The Cauchy surface is denoted by the dashed red curve. }
\label{fig:Cauchy}
\end{figure}

For the matching we find it most convenient to choose the complete set of modes that consists of modes that are positive frequency on the future horizon $H^{+}$ (labeled by $f^{H^{+}}_{\omega \ell m}$)  along with modes that are positive frequency on future null infinity $\mathscr{I}^{+}$ (labeled by $f^{\mathscr{I}^{+}}_{\omega \ell m}$).

The Unruh state is defined in Schwarzschild spacetime by considering modes that are positive frequency on past null infinity $\mathscr{I}^{-}$ and labeled by $f^{\mathscr{I}^{-}}_{\omega \ell m}$ along with modes labeled $f^{K}_{\omega \ell m}$ that are positive frequency with respect to the Kruskal null coordinate
\begin{equation} U = - \frac{e^{-\kappa u}}{\kappa}   \label{U-Kruskal}  \end{equation}
on the past horizon $H^{-}$.  Here $\kappa = (4 M)^{-1}$ is the surface gravity of the black hole and $M$ is its mass.  As shown below, the modes $f^{K}_{\omega \ell m}$  can be expanded in terms of modes that are positive frequency with respect to $u$ on $H^{-}$ and that we label as $f^{H^{-}}_{\omega \ell m}$.  Note that together the modes $f^{\mathscr{I}^{-}}_{\omega \ell m}$ and $f^{H^{-}}_{\omega \ell m}$ form a complete set of modes in the part of Schwarzschild spacetime that is outside the past and future horizons.


\subsection{Properties of mode functions in Schwarzschild spacetime}

The general form of the mode functions in Schwarzschild spacetime along with the equation satisfied by the radial mode functions is given in~\eqref{modes-outside}.  The two
complete sets of modes that we work with are defined for values of $\omega$ in the range $0 \le \omega < \infty$.
The radial functions for these two sets of modes can be defined in terms of two linearly independent solutions to the radial mode equation~\eqref{modes-outside} with the properties
\begin{subequations} \begin{eqnarray} \chi^\infty_R &\to& e^{i \omega r_*} \qquad r_* \to \infty \;,\label{chi-infinity-R-def} \\
\chi^\infty_L &\to& e^{ -i \omega r_*} \qquad r_* \to \infty \label{chi-infinity-L-def}\;. \end{eqnarray} \label{chiRL-def} \end{subequations}
Near the event horizon they have the behaviors~\cite{rigorous}
\begin{subequations} \begin{eqnarray} \chi^\infty_R &\to& E_R(\omega) e^{i \omega r_*} + F_R(\omega) e^{-i \omega r_*} \;, \label{chiR-hor} \\
      \chi^\infty_L &\to& E_L(\omega) e^{i \omega r_*} + F_L(\omega) e^{-i \omega r_*} \;, \label{chiL-hor} \end{eqnarray} \label{chiRL-hor} \end{subequations}
where $E_R$, $E_L$, $F_R$, and $F_L$ are scattering parameters that can be determined numerically.\footnote{The subscripts $r$ and $l$ in~\cite{rigorous} have been changed here to $R$ and $L$ respectively.}

As discussed above, one complete set of modes consists of modes that are positive frequency  on the past horizon and modes that are positive frequency on past null infinity.
On the past horizon $\psi_{\omega \, \ell}^{H^{-}} = e^{-i \omega u_s}$ and $\psi_{\omega \, \ell}^{\mathscr{I}^{-}} = 0$ while on past null infinity $\psi_{\omega \, \ell}^{\mathscr{I}^{-}} = e^{-i \omega v}$ and $\psi_{\omega \, \ell}^{H^{-}} = 0$.
 These modes  are scattered by the effective potential.  It is straight-forward to show that~\cite{rigorous}
\begin{subequations} \begin{eqnarray} \chi^{H^{-}}_{\omega \ell} &=& \frac{ \chi^\infty_R}{E_R} \;, \label{chi-H-minus-def} \\
         \chi_{\omega \ell}^{\mathscr{I}^{-}} &=&  \chi^\infty_L - \frac{E_L}{E_R} \chi^\infty_R  \;. \label{chi-scri-minus-def}
\end{eqnarray} \label{chi-H-scri-def}\end{subequations}

The other complete set of modes that we will use consists of modes that are positive frequency on the future horizon and those that are positive frequency on future null infinity.  On the future horizon $\psi_{\omega \, \ell}^{H^{+}} = e^{-i \omega v}$ and $\psi_{\omega \, \ell}^{\mathscr{I}^{+}} = 0$ while on future null infinity $\psi_{\omega \, \ell}^{H^{+}} = 0$ and $\psi_{\omega \, \ell}^{\mathscr{I}^{+}} = e^{-i \omega u_s}$.
Going backwards in time these modes scatter due to the effective potential.   It is straight-forward to show that
\begin{subequations} \begin{eqnarray} \chi^{H^{+}}_{\omega \ell} &=& \frac{1}{F_L} \chi^\infty_L \;, \label{chi-H-plus-def}  \\
\chi^{\mathscr{I}^{+}}_{\omega \ell} &=& \chi^\infty_R - \frac{F_R}{F_L} \chi^\infty_L  \;.  \label{chi-scri-plus-def} \end{eqnarray} \end{subequations}


Since the method used to compute the stress-energy tensor involves subtracting the unrenormalized stress-energy tensor for the Unruh state it is useful to write the modes associated with this state, $f^{H^{-}}$ and $f^{\mathscr{I}^{-}}$ in terms of $f^{H^{+}}$  and $f^{\mathscr{I}^{+}}$ .
\begin{subequations} \begin{eqnarray} \chi^{H^{-}}_{\omega \ell} &=& \frac{1}{E_R} ( F_R \,  \chi^{H^{+}}_{\omega \ell} +  \chi^{\mathscr{I}^{+}}_{\omega \ell} )  \;, \label{chi-H-minus-2} \\
\chi_{\omega \ell}^{\mathscr{I}^{-}} &=& \frac{1}{E_R} (\chi^{H^{+}}_{\omega \ell} - E_L \, \chi^{\mathscr{I}^{+}}_{\omega \ell})  \;. \label{chi-scri-minus-2}  \end{eqnarray} \label{chi-scri-H-minus} \end{subequations}
The modes $f^K$ which are positive frequency on $H^{-}$ with respect to the Kruskal time coordinate can be expanded in terms of the $f^{H^{-}}$ modes.  The result
is given in~\eqref{alphaK-betaK}.


\section{Matching Coefficients}

\subsection{General Formulas}

In this section general formulas are derived for the matching coefficients used in an expansion of the modes of a massless minimally coupled scalar field for the {\it in} vacuum state in the collapsing null shell spacetime in terms of a complete set of modes in Schwarzschild spacetime in the region outside the past and future horizons.  These can be used in the computation of the stress-energy tensor, $\left\langle in \middle| T_{ab} \middle| in \right\rangle$ for the scalar field in the part of the collapsing null shell spacetime that is outside of the shell and outside of the event horizon.

The expansion of the {\it in} mode functions has the form
\begin{eqnarray}
f^{in}_{\omega \ell m} = \sum_{\ell'=0}^{\infty}\sum_{m'=-\ell'}^{\ell'} & \int_{0}^{\infty} d \omega'
\Big[ A^{\mathscr{I^+}}_{\omega \ell m\omega' \ell' m'} f^{\mathscr{I}^+}_{\omega'\ell' m'} + B^{\mathscr{I^+}}_{\omega l m\omega' \ell' m'} (f^{\mathscr{I}^+}_{\omega' \ell' m'})^{*}
 \nonumber \\ & + A^{H^+}_{\omega \ell m \omega' \ell' m'} f^{H^+}_{\omega' \ell' m'} + B^{H^+}_{\omega \ell m\omega'\ell' m'} (f^{H^+}_{\omega' \ell' m'})^{*}\Big]\;. \label{General-in-modes}
\end{eqnarray}
The matching coefficients are found using the scalar product in~\eqref{scalar-products} and the orthonormality of the modes $f^{(\mathscr{I}^{+}, \, H^{+})}$ with respect to this scalar product.  The result is
\begin{subequations} \begin{eqnarray}
A^{(\mathscr{I^+}, H^+)}_{\omega \ell m\omega' \ell' m'} &=&( f^{\rm in}_{\omega \ell m} , f^{(\mathscr{I}^+, H^+)}_{\omega' \ell' m'}) \;, \\
B^{(\mathscr{I^+}, H^+)}_{\omega \ell m\omega' \ell' m'} &=& -( f^{\rm in}_{\omega \ell m} , (f^{(\mathscr{I}^+, H^+)}_{\omega' \ell' m'})^{*}) \;.
\end{eqnarray} \end{subequations}
For the Cauchy surface we consider, ~\eqref{scalar-products} reduces to integrals of the form
\begin{eqnarray}
    \int du\int d\Omega r^2\overset{\leftrightarrow}{\partial_u}\;, \quad  \quad \int dv\int d\Omega r^2\overset{\leftrightarrow}{\partial_v}\;.
\end{eqnarray}
On the hypersurfaces where these integrals are computed, the following properties for spherical harmonics can be used
\begin{subequations} \begin{eqnarray}
    \int d\Omega Y_{lm}(\theta, \phi)Y_{l'm'}^{*}(\theta, \phi)&=&\delta_{l,l'} \delta_{m,m'},  \\
    \int d\Omega Y_{lm}(\theta, \phi)Y_{l'm'}(\theta, \phi)&=&(-1)^{m}\delta_{l,l'} \delta_{m,-m'} \;.
\end{eqnarray} \end{subequations}
As a result we can write
\begin{subequations} \begin{eqnarray} A^{(\mathscr{I^+}, H^+)}_{\omega l m\omega' l' m'} & =& \delta_{l,l'} \delta_{m,m'} A^{(\mathscr{I^+}, H^+)}_{\omega \omega' \ell}  \;, \label{A-short-def} \\
          B^{(\mathscr{I^+}, H^+)}_{\omega l m\omega' l' m'} & =& (-1)^m \delta_{l,l'} \delta_{m,-m'} B^{(\mathscr{I^+}, H^+)}_{\omega \omega' \ell}  \;, \label{B-short-def}
\end{eqnarray} \label{A-B-short-def} \end{subequations}
and
\begin{equation}
f^{in}_{\omega \ell m} =  \frac{Y_{\ell m}}{r \sqrt{4 \pi}} \int_{0}^{\infty} \frac{d \omega'}{\sqrt{\omega'}}
\Big[ A^{\mathscr{I^+}}_{\omega \omega' \ell} \psi^{\mathscr{I}^+}_{\omega'\ell} + B^{\mathscr{I^+}}_{\omega \omega' \ell} (\psi^{\mathscr{I}^+}_{\omega' \ell})^{*}
 + A^{H^+}_{\omega \omega' \ell} \psi^{H^+}_{\omega' \ell} + B^{H^+}_{\omega \omega'\ell} (\psi^{H^+}_{\omega' \ell})^{*}\Big]\;. \label{General-in-modes-2}
\end{equation}
From this expression one can see that if, at small $\omega'$, the matching coefficients go like $\frac{1}{\sqrt{\omega'}}$ then there is an infrared divergence
in the integral and it is not obvious how to deal with it.  For this reason, we use integrations by parts in some of the computations of the matching coefficients below to avoid this difficulty.  For Schwarzschild spacetime in 4D, our results when substituted into~\eqref{General-in-modes-2}, do not give infrared divergences.

The contribution to the matching coefficients from the three segments of the Cauchy surface in Fig.~\ref{fig:Cauchy} are
\begin{subequations} \begin{eqnarray}
A^{(\mathscr{I^+}, H^+)}_{\omega \omega' l} & =& \left(A^{(\mathscr{I^+}, H^+)}_{\omega \omega' l} \right)_{H^{+}} + \left(A^{(\mathscr{I^+}, H^{+})}_{\omega \omega' l} \right)_{v_0}
   +  \left(A^{(\mathscr{I^+}, H^+)}_{\omega \omega' l} \right)_{\mathscr{I}^{-}} \;, \label{A-parts} \\
\left(A^{(\mathscr{I^+}, H^+)}_{\omega \omega' l} \right)_{H^{+}} &=&
-\frac{i}{4\pi\sqrt{\omega\omega'}}\int_{-\infty}^{v_{0}}dv\;\psi^{in}_{\omega l}(u = v_H,v) \overset{\leftrightarrow}{\partial_v}[\psi_{\omega' l}^{(\mathscr{I^+}, H^+)}(u_s = \infty,v)]^{*} \;, \\
\left(A^{(\mathscr{I^+}, H^{+})}_{\omega \omega' l} \right)_{v_0} &=& -\frac{i}{4\pi\sqrt{\omega\omega'}} \int_{-\infty}^{v_{H}}du\;\psi^{in}_{\omega l}(u,v_0) \overset{\leftrightarrow}{\partial_u}[\psi_{\omega' l}^{(\mathscr{I^+}, H^+)}(u_s(u),v_0)]^{*} \;, \label{A-I+-v0-1} \\
 \left(A^{(\mathscr{I^+}, H^+)}_{\omega \omega' l} \right)_{\mathscr{I}^{-}} &=& -\frac{i}{4\pi\sqrt{\omega\omega'}}\int_{v_{0}}^{\infty}dv\;\psi^{in}_{\omega l}(u =-\infty,v) \overset{\leftrightarrow}{\partial_v}[\psi_{\omega' l}^{(\mathscr{I^+}, H^+)}(u_s=-\infty,v) ]^{*}\;, \end{eqnarray} \end{subequations}
and
\begin{subequations} \begin{eqnarray}
B^{(\mathscr{I^+}, H^+)}_{\omega \omega' l} & =& \left(B^{(\mathscr{I^+}, H^+)}_{\omega \omega' l} \right)_{H^{+}} + \left(B^{(\mathscr{I^+}, H^{+})}_{\omega \omega' l} \right)_{v_0}
   +  \left(B^{(\mathscr{I^+}, H^+)}_{\omega \omega' l} \right)_{\mathscr{I}^{-}} \;, \label{B-parts} \\
\left(B^{(\mathscr{I^+}, H^+)}_{\omega \omega' l} \right)_{H^{+}} &=&
\frac{i}{4\pi\sqrt{\omega\omega'}}\int_{-\infty}^{v_{0}}dv\;\psi^{in}_{\omega l}(u = v_H,v) \overset{\leftrightarrow}{\partial_v} \psi_{\omega' l}^{(\mathscr{I^+}, H^+)}(u_s = \infty,v) \;, \\
\left(B^{(\mathscr{I^+}, H^{+})}_{\omega \omega' l} \right)_{v_0} &=& \frac{i}{4\pi\sqrt{\omega\omega'}} \int_{-\infty}^{v_{H}}du\;\psi^{in}_{\omega l}(u,v_0) \overset{\leftrightarrow}{\partial_u}\psi_{\omega' l}^{(\mathscr{I^+}, H^+)}(u_s(u),v_0) \label{B-I+-v0-1} \;, \\
 \left(B^{(\mathscr{I^+}, H^+)}_{\omega \omega' l} \right)_{\mathscr{I}^{-}} &=& \frac{i}{4\pi\sqrt{\omega\omega'}}\int_{v_{0}}^{\infty}dv\;\psi^{in}_{\omega l}(u=-\infty,v) \overset{\leftrightarrow}{\partial_v}\psi_{\omega' l}^{(\mathscr{I^+}, H^+)}(u_s=-\infty,v) \;. \end{eqnarray} \end{subequations}
Note that the modes $\psi_{\omega' l}^{(\mathscr{I^+})} $ vanish on $H^{+}$ so
\begin{equation} \left(A^{\mathscr{I^+}}_{\omega \omega' l} \right)_{H^{+}} = \left(B^{\mathscr{I^+}}_{\omega \omega' l} \right)_{H^{+}} = 0 \;. \end{equation}

As discussed in Sec.~\ref{sec:method-Tab}, it is necessary to specify $\psi^{\rm in}_{\omega \ell}$ on the part of the future horizon with $v < v_0$.  This part is
causally disconnected from the region with $v > v_0$ outside the horizon so the only necessary constraint is that the mode functions should be continuous at the
point on the future horizon where $v = v_0$.  The simplest way to accomplish this is to take
\begin{equation} \psi^{\rm in}_{\omega \ell}(u = v_H, v) = \psi^{\rm in}_{\omega \ell}(u = v_H, v_0) \;. \end{equation}
With this choice it turns out to be useful to write the contribution to the matching coefficients from $H^+$ in the form
\begin{subequations} \begin{eqnarray} \left(A^{(\mathscr{I^+}, H^+)}_{\omega \omega' l} \right)_{H^{+}} &=&  \frac{i}{4 \pi \sqrt{\omega \omega'}} \psi^{\rm in}_{\omega \ell}(v_H,v_0) e^{i \omega' v_0}
  - \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{e^{i \omega' v_0}}{\omega' - i \epsilon} \psi^{\rm in}_{\omega \ell}(v_H,v_0)
   \;, \\
  \left(B^{(\mathscr{I^+}, H^+)}_{\omega \omega' l} \right)_{H^{+}} &=&  -\frac{i}{4 \pi \sqrt{\omega \omega'}} \psi^{\rm in}_{\omega \ell}(v_H,v_0) e^{-i \omega' v_0}
 + \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{e^{i \omega' v_0}}{\omega' + i \epsilon} \psi^{\rm in}_{\omega \ell}(v_H,v_0)  \;, \end{eqnarray} \end{subequations}
where for each integral an integration by parts has been done and an integrating factor $0 < \epsilon \ll 1$ has been included to make the integrals converge.

To avoid infrared divergences in~\eqref{General-in-modes-2} it is useful to subtract and then add back the quantity $e^{-i \omega v_0}$ from $\psi^{\rm in}_{\omega \ell}$
in~\eqref{A-I+-v0-1} and~\eqref{B-I+-v0-1}.  Then after integrations by parts the contributions from the surface $v = v_0$ can be written as
\begin{subequations} \begin{eqnarray}  \left(A^{ H^{+}}_{\omega \omega' l} \right)_{v_0} &=& - \frac{i}{4 \pi \sqrt{\omega \omega'} }  \psi^{\rm in}_{\omega \ell}(u_H,v_0) e^{i \omega' v_0} +  \frac{i}{4 \pi \sqrt{\omega \omega'}\, F_L^{*}(\omega',\ell)}  e^{-i (\omega-\omega') v_0}  \nonumber \\
   & & + \frac{i}{2 \pi \sqrt{\omega \omega'} } \int_{-\infty}^{v_H} du \, \left[\partial_u \psi^{\rm in}_{\omega \ell}(u,v_0) \right] \psi^{H^{+} *}_{\omega' \ell}(u,v_0) \;, \\
   \left(B^{ H^{+}}_{\omega \omega' l} \right)_{v_0} &=&  \frac{i}{4 \pi \sqrt{\omega \omega'} }  \psi^{\rm in}_{\omega \ell}(v_H,v_0) e^{-i \omega' v_0} -  \frac{i}{4 \pi \sqrt{\omega \omega'}\, F_L(\omega',\ell)}  e^{-i (\omega+\omega') v_0}  \nonumber \\
   & & - \frac{i}{2 \pi \sqrt{\omega \omega'} } \int_{-\infty}^{v_H} du \, \left[\partial_u \psi^{\rm in}_{\omega \ell}(u,v_0) \right] \psi^{H^{+}}_{\omega' \ell} (u,v_0)\;, \\
  \left(A^{ \mathscr{I}^{+}}_{\omega \omega' l} \right)_{v_0} &=& - \frac{i}{4 \pi \sqrt{\omega \omega'} } \frac{F^{*}_R(\omega',\ell)}{F^{*}_L(\omega', \ell)} e^{-i (\omega - \omega') v_0} \nonumber \\
 && + \frac{i}{2 \pi \sqrt{\omega \omega'}} \int_{-\infty}^{v_H} du \,  [\psi^{\rm in}_{\omega \ell}(u, v_0) - e^{-i \omega v_0}] \partial_u \psi^{\mathscr{I}^+ *}_{\omega \ell} \;, \\
  \left(B^{ \mathscr{I}^{+}}_{\omega \omega' l} \right)_{v_0} &=&  \frac{i}{4 \pi \sqrt{\omega \omega'} } \frac{F_R(\omega',\ell)}{F_L(\omega', \ell)} e^{-i (\omega + \omega') v_0} \nonumber \\
 && - \frac{i}{2 \pi \sqrt{\omega \omega'}} \int_{-\infty}^{v_H} du \, \left[ \psi^{\rm in}_{\omega \ell}(u, v_0) - e^{-i \omega v_0} \right] \partial_u \psi^{\mathscr{I}^+}_{\omega \ell} \;.
      \end{eqnarray} \end{subequations}

Since the spatial dependence of the modes is very simple on $\mathscr{I}^{-}$ it is possible to evaluate the integrals for the contributions to the matching parameters from there. After integrating by parts, we find
\begin{subequations} \begin{eqnarray}  \left(A^{ H^{+}}_{\omega \omega' l} \right)_{\mathscr{I}^{-}} &=& - \frac{i}{4 \pi \sqrt{\omega \omega'} \, F^{*}_L(\omega',\ell)} e^{-i (\omega-\omega') v_0}
 + \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{1}{F_L^{*}(\omega',\ell)}\frac{e^{i(\omega'-\omega)v_0}}{\omega'-\omega+i \epsilon}  \;, \\
\left(B^{ H^{+}}_{\omega \omega' l} \right)_{\mathscr{I}^{-}} &=& \frac{i}{4 \pi \sqrt{\omega \omega'} \, F_L(\omega',\ell)} e^{-i (\omega+\omega') v_0}
 - \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{1}{F_L(\omega',\ell)}\frac{e^{-i(\omega+\omega')v_0}}{\omega'+\omega-i \epsilon}  \;, \\
\left(A^{ \mathscr{I}^{+}}_{\omega \omega' l} \right)_{\mathscr{I}^{-}} &=& \frac{i}{4 \pi \sqrt{\omega \omega'}}\frac{F_R^{*}(\omega',\ell)}{ F_L^{*}(\omega',\ell)} e^{-i (\omega-\omega') v_0}
 - \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{F_R^{*}(\omega',\ell)}{F_L^{*}(\omega',\ell)}\frac{e^{-i(\omega-\omega')v_0}}{\omega'-\omega+i \epsilon}  \;, \\
\left(B^{ \mathscr{I}^{+}}_{\omega \omega' l} \right)_{\mathscr{I}^{-}} &=& -\frac{i}{4 \pi \sqrt{\omega \omega'}}\frac{F_R(\omega',\ell)}{ F_L(\omega',\ell)} e^{-i (\omega+\omega') v_0}
 + \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{F_R(\omega',\ell)}{F_L(\omega',\ell)}\frac{e^{-i(\omega+\omega')v_0}}{\omega'+\omega-i \epsilon}  \;.
\end{eqnarray} \end{subequations}

Combining these results together, the general formulas for the matching coefficients are
\begin{subequations} \begin{eqnarray} A^{ H^{+}}_{\omega \omega' l} &=& - \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{e^{i \omega' v_0}}{\omega' - i \epsilon} \psi^{\rm in}_{\omega \ell}(v_H,v_0)
   + \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{1}{F_L^{*}(\omega',\ell)}\frac{e^{i(\omega'-\omega)v_0}}{\omega'-\omega+i \epsilon}  \nonumber \\
   & & + \frac{i}{2 \pi \sqrt{\omega \omega'} } \int_{-\infty}^{v_H} du \, \left[\partial_u \psi^{\rm in}_{\omega \ell}(u,v_0) \right] \psi^{H^{+} *}_{\omega' \ell}(u,v_0) \;, \label{A-H-gen-mat}\\
      B^{ H^{+}}_{\omega \omega' l} &=& \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{e^{-i \omega' v_0}}{\omega' + i \epsilon} \psi^{\rm in}_{\omega \ell}(v_H,v_0)
      - \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{1}{F_L(\omega',\ell)}\frac{e^{-i(\omega+\omega')v_0}}{\omega'+\omega-i \epsilon} \nonumber \\
      & & - \frac{i}{2 \pi \sqrt{\omega \omega'} } \int_{-\infty}^{v_H} du \, \left[\partial_u \psi^{\rm in}_{\omega \ell}(u,v_0) \right] \psi^{H^{+}}_{\omega' \ell} (u,v_0)\;, \label{B-H-gen_mat}\\
     A^{ \mathscr{I}^{+}}_{\omega \omega' l} &=&  - \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{F_R^{*}(\omega',\ell)}{F_L^{*}(\omega',\ell)}\frac{e^{-i(\omega-\omega')v_0}}{\omega'-\omega+i \epsilon} \nonumber \\  & &
       -\frac{i}{2 \pi \sqrt{\omega \omega'}} \int_{-\infty}^{v_H} du \, \left[  \psi^{\rm in}_{\omega \ell}(u, v_0) - e^{-i \omega v_0} \right] \partial_u \psi^{\mathscr{I}^+  *}_{\omega' \ell} \;,  \label{A-I-gen-mat} \\
     B^{ \mathscr{I}^{+}}_{\omega \omega' l} &=&  \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{F_R(\omega',\ell)}{F_L(\omega',\ell)}\frac{e^{-i(\omega+\omega')v_0}}{\omega'+\omega-i \epsilon} \nonumber \\
  & &     + \frac{i}{2 \pi \sqrt{\omega \omega'}} \int_{-\infty}^{v_H} du \, \left[ \psi^{\rm in}_{\omega \ell}(u, v_0) - e^{-i \omega v_0} \right] \partial_u \psi^{\mathscr{I}^+}_{\omega' \ell} \;.
     \label{B-I-gen-mat} \end{eqnarray} \label{gen-mat}\end{subequations}

\subsection{Expansion of the Kruskal modes}

As discussed in Sec.~\ref{sec:method-Tab}, our method for renormalizing the stress-energy tensor involves subtracting the unrenormalized stress-energy tensor for the Unruh modes.  For this purpose the best way to do this is to express the modes $f^K_{\omega \ell m}$ that are positive frequency on the past horizon with respect to the Kruskal time coordinate in terms of the modes $f^{H^{-}}_{\omega \ell m}$ that are positive frequency with respect to the usual time coordinate $t_s$ on the past horizon. Then the relation~\eqref{chi-H-minus-2} can be used to express $f^K_{\omega \ell m}$ in terms of the modes $f^{(\mathscr{I}^{+}\,H^{+})}_{\omega \ell m}$.
The initial expansion can be written as
\begin{equation} f^K_{\omega \ell m} =  \int_0^\infty d \omega' \, \left[  \alpha^K_{\omega \omega' \ell} f^{H^{-}}_{\omega' \ell m} + \beta^K_{\omega \omega' \ell}   f^{H^{-}*}_{\omega' \ell m} \right]  \;. \label{fK} \end{equation}
 The Bogolubov coefficients can be obtained using the scalar product~\eqref{scalar-products}
   with a Cauchy surface consisting of the union of past null infinity and the past horizon in Schwarzschild spacetime.  Integrating over the angular coordinates one finds that
   the Bogolubov coefficients can be written in the form~\eqref{A-B-short-def} with $\alpha$ replacing $A$ and $\beta$ replacing $B$.
   Integrating the remaining integrals over $u_s$ by parts one finds that\footnote{This calculation was originally done in~\cite{BEC-2013} but note that there is a mistake in the results.  The expressions in that paper are missing a factor of $(4M)^{\pm i \omega'}$. }
\begin{subequations} \begin{eqnarray} \alpha^K_{\omega_K \omega' \ell} &=& - \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega_K}} (4 M)^{1+i 4 M \omega'}  \int_{-\infty}^0 d U_K e^{-i \omega_K U_K} (-U_K)^{-1 - i 4 M \omega'} \nonumber \\
&=& \frac{1}{2 \pi} \sqrt{\frac{\omega'}{\omega_K}} \, (4M)^{1 +i 4 M \omega'} \frac{\Gamma(\delta - i 4 M \omega)}{ (-i \omega_K + \epsilon)^{-i 4 M \omega'}} \;, \label{alphaK} \\
\beta^K_{\omega_K \omega' \ell} &=&  \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega_K}} (4 M)^{1-i 4 M \omega'}  \int_{-\infty}^0 d U_K e^{-i \omega_K U_K} (-U_K)^{-1 + i 4 M \omega'} \nonumber \\
&=& \frac{1}{2 \pi} \sqrt{\frac{\omega'}{\omega_K}} \, (4M)^{1 -i 4 M \omega'} \frac{\Gamma(\delta + i 4 M \omega)}{ (-i \omega_K + \epsilon)^{i 4 M \omega'}} \;. \label{betaK}
\end{eqnarray} \label{alphaK-betaK} \end{subequations}
Here $\delta$ and $\epsilon$ are integrating factors with $0 < \delta \ll 1$ and $0 < \epsilon \ll 1$.  Note that the Bogolubov coefficients are independent of the value of $\ell$.  This is because the effective potential vanishes on $H^{-}$ which is the surface where the integrals are being computed.

\subsection{2D example}

In this section we will illustrate the matching for the case of a 2D spacetime which has a perfectly reflecting mirror at $r = 0$.  The metric
inside the shell is the flat space metric
\begin{equation}
ds^2 = -dt^2 + dr^2    \;, \label{metric-flat-2D}
\end{equation}
and the metric outside the shell is the Schwarzschild metric
\begin{equation}
ds^2 = -\left(1-\frac{2M}{r} \right) dt_s^2
+ \left( 1-\frac{2M}{r} \right)^{-1} dr^2   \;. \label{metric-sch-2D}
\end{equation}
The Penrose diagram is the same as in the 4D case as is the definition of the radial null coordinates $u$, $u_s$, and $v$ and the relation between $u$ and $u_s$.

The general form of the mode functions is
\begin{equation} f_\omega = \frac{\psi_\omega}{\sqrt{4 \pi \omega}}  \;.  \label{f-form-2D} \end{equation}
There is no scattering for the massless minimally coupled scalar field modes in 2D so
\begin{equation} E_R = F_L = 1\;, \qquad E_L = F_R = 0 \;. \label{scattering-coeff-2D} \end{equation}

Inside the shell the {\it in} modes are
\begin{equation} \psi^{\rm in}_\omega =   e^{-i \omega v} - e^{-i \omega u} \;.  \label{psi-in-2D} \end{equation}
In the region outside the shell the spacetime is the 2D version of Schwarzschild spacetime and the modes are
\begin{subequations} \begin{eqnarray}   \psi^{\mathscr{I}^{+}}_\omega &=& \psi^{H^{-}}_\omega = e^{- i \omega u_s}  \;, \label{right-moving-2D} \\
           \psi^{H^{+}}_\omega &=& \psi^{\mathscr{I}^{-}}_\omega  = e^{-i \omega v} \;. \label{left-moving-2D} \end{eqnarray} \label{sch-modes-2D} \end{subequations}

The expansion for the {\it in} modes is similar to the 4D case except there are no parameters $\ell$ and $m$ related to the spherical harmonics.  Thus
\begin{eqnarray}
f^{in}_{\omega} & =  & \int_{0}^{\infty} d \omega'
\Big[ A^{H^+}_{\omega \omega' } f^{H^+}_{\omega' } + B^{H^+}_{\omega \omega'} (f^{H^+}_{\omega' })^{*}
   + A^{\mathscr{I^+}}_{\omega \omega'} f^{\mathscr{I}^+}_{\omega'} + B^{\mathscr{I^+}}_{\omega \omega' } (f^{\mathscr{I}^+}_{\omega' })^{*}\Big]\;.\label{2D-in-mode-expansion}
\end{eqnarray}
The matching coefficients are given by substituting~\eqref{scattering-coeff-2D},~\eqref{psi-in-2D}, and~\eqref{sch-modes-2D} into~\eqref{gen-mat}.  It is then easy to show that
\begin{eqnarray} \left[ B^{H^{+}}_{\omega,  \omega'} f^{H^+\, *}_{\omega' } \right]_{\omega' \to - \omega'} &=& - A^{H^{+}}_{\omega,  \omega'} f^{H^+}_{\omega' } \;, \nonumber \\
     \left[ B^{\mathscr{I}^{+}}_{\omega,  \omega'} f^{\mathscr{I}^{+} \, *}_{\omega' } \right]_{\omega' \to - \omega'} &=& - A^{ \mathscr{I}^{+}}_{\omega,  \omega'} f^{\mathscr{I}^{+}}_{\omega' } \;,
   \label{B-w-minus-w} \end{eqnarray}
where the quantities on the right hand side are to be evaluated at $\omega' < 0$.  As a result
\begin{eqnarray}
f^{in}_{\omega} & =  & \int_{-\infty}^{\infty} d \omega'
\Big[   A^{H^+}_{\omega \omega' } f^{H^+}_{\omega' } + A^{\mathscr{I^+}}_{\omega \omega'} f^{\mathscr{I}^+}_{\omega'} \Big]\;. \label{2D-in-mode-expansion-2}
\end{eqnarray}
Because $\psi^{H^{+}}_{\omega'}$ does not depend on $u$, the integral in~\eqref{A-H-gen-mat} is trivial to evaluate and one finds that
\begin{equation} A^{H^+}_{\omega \omega'} = - \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{e^{-i (\omega - \omega') v_0}}{\omega' - i \epsilon}
   + \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{e^{i(\omega'-\omega)v_0}}{\omega'-\omega+i \epsilon}
  \;. \label{A-H-2D} \end{equation}

To see what the contribution to $f^{\rm in}_\omega$ is from the $f^{H^{+}}_\omega$ modes, first substitute ~\eqref{A-H-2D} into~\eqref{2D-in-mode-expansion} along with~\eqref{f-form-2D}
and~\eqref{psi-in-2D} with the result
\begin{eqnarray} \left(f^{\rm in}_\omega \right)_{H^{+}}  &=&  \frac{i e^{-i \omega v_0}}{2 \pi \sqrt{4 \pi \omega}} \int_{-\infty}^\infty d w' \; \left[ e^{i \omega' (v_0 - v)} \left(-\frac{1}{\omega' - i \epsilon}  +   \frac{1}{\omega' - \omega + i \epsilon} \right) \right] \nonumber \\
 & = &  \frac{ e^{-i \omega v_0}}{\sqrt{4 \pi \omega}}  \theta(v_0-v) + \frac{ e^{-i \omega v}}{\sqrt{4 \pi \omega}}  \theta(v-v_0) \;. \label{f-in-H+-2D}\end{eqnarray}

We next consider the contribution of the $f^{\mathscr{I}^{+}}$ modes.  The matching coefficient in~\eqref{A-I-gen-mat} is
\begin{eqnarray} A^{\mathscr{I}^{+}}_{\omega \omega'} &=& -\frac{1}{2 \pi }\sqrt{\frac{\omega'}{\omega}} \int_{-\infty}^{v_H} du \, e^{-i \omega u}  e^{i \omega' u_s(u)} \frac{d u_s}{du} \nonumber \\
&=&
-\frac{1}{2 \pi} \sqrt{\frac{\omega'}{\omega}}  \int_{-\infty}^{v_H} du \, e^{-i (\omega - \omega') u}  \left(\frac{v_H-u}{4 M} \right)^{-i 4 M \omega'} \left[ 1 + \frac{4M}{v_h-u} \right]  \;. \end{eqnarray}
Changing variables to $x = v_H -u$ and performing an integration by parts gives
\begin{eqnarray} A^{\mathscr{I}^{+}}_{\omega \omega'} &=& \frac{i}{2 \pi} \sqrt{\omega \omega'} e^{-i(\omega-\omega') v_H} \frac{(4 M)^{1+i 4 M \omega'}}{i (\omega'-\omega)+\epsilon} \int_0^\infty dx \; e^{i(\omega-\omega')x-\epsilon x}  \;
  x^{-i 4 M \omega' -1 + \delta} \nonumber \\
  &=& \frac{i}{2 \pi} \sqrt{\omega \omega'} e^{-i(\omega-\omega') v_H} (4 M)^{1+i 4 M \omega'} \frac{\Gamma(\delta - i 4 M \omega')}{[i (\omega'-\omega)+ \epsilon]^{1-i4M\omega'} }  \;. \label{Ascrp2D} \end{eqnarray}
Note that two integrating factors have been used with $0 < \epsilon \ll 1$ and $0 < \delta \ll 1$.

To find the contribution to $f^{\rm in}_\omega$ from the $f^{\mathscr{I}^{+}}_\omega$ modes, first substitute~\eqref{Ascrp2D}  into~\eqref{2D-in-mode-expansion} with the result
\begin{eqnarray} \left(f^{\rm in}_\omega\right)_{\mathscr{I}^{+}} &=& \frac{i 4 M \sqrt{\omega}}{2 \pi \sqrt{4 \pi}} e^{-i \omega v_H}
\int_{-\infty}^\infty d \omega'   e^{i \omega' (v_H-u_s)} (4 M)^{i 4 M \omega'} \frac{\Gamma(\delta - i 4 M \omega')}{[i (\omega'-\omega) + \epsilon]^{1-i 4 M \omega'}} \;. \label{I1-a} \end{eqnarray}
Note that the denominator has an essential singularity in the upper half $\omega'$ plane while the Gamma function has simple poles in the lower half plane at
 \begin{equation}  \omega' = - \frac{i \delta}{4 M} \;, \end{equation}
  and
  \begin{equation} \omega' = - \frac{i n}{4 M} \;, \qquad n = 1, 2, \ldots \end{equation}
In the complex plane at large $|\omega'|$ Sterling's approximation gives
\begin{equation}  \Gamma( -i 4 M \omega') \approx \sqrt{2 \pi} e^{i 4 M \omega'} e^{(-i 4 M \omega' - 1/2) \log(-i 4 M \omega')}   \;. \end{equation}
Using the usual change of variables $\omega' = R e^{i \theta}$, with $R > 0$, it is straight-forward to show that the dominant contribution to the integrand of $I_1$ in the large $R$ limit comes from the factor $e^{4 M R \sin \theta \, \log R}$ and therefore one must close in the lower half plane.  This means there is no contribution from the essential singularity but there is a contribution from each pole of the Gamma function.  At these poles it is straight-forward to show that
\begin{equation} \Gamma(\delta -i 4 M \omega) \to \frac{(-1)^n}{n! (n - i 4 M \omega)}  \;, \qquad n = 0, 1, 2, \ldots  \end{equation}
Then
\begin{equation}  \left(f^{\rm in}_\omega\right)_{\mathscr{I}^{+}}  = \frac{4Mi\sqrt{\omega}}{\sqrt{4 \pi}} e^{- i \omega v_H} \sum_{n=0}^\infty \frac{(-1)^n}{n!} (n-i 4 M \omega)^{n-1}
    \left[\exp\left(\frac{(v_H-u_s)}{4 M} \right)\right]^n  \;. \label{finu} \end{equation}
Because the general solutions to the 2D mode equation in Schwarzschild spacetime are of the form $\psi = g(u_s) + h(v)$ with $g$ and $h$ arbitrary functions, the exact solution for the {\it in} modes is
\begin{equation} \left(f^{\rm in}_\omega\right)_{\mathscr{I}^{+}} = - \frac{e^{-i \omega u(u_s)}}{\sqrt{4\pi\omega}}  = -\frac {e^{-i \omega v_H}}{\sqrt{4\pi\omega}} \, \exp \left\{ i 4 M \omega\;  W\left[\exp\left(\frac{(v_H-u_s)}{4M}\right)\right] \right\} \;, \label{f(u)Lambert}\end{equation}
where~\eqref{u-us} has been used and $W(z)$ is the Lambert W function.
  To make a comparison between ~\eqref{finu} and ~\eqref{f(u)Lambert}, one needs to write the latter in terms of a series. This has been done in~\cite{Corless}.
  An alternative derivation is given in Appendix~\ref{appendix-A}.  The result is
\begin{equation} e^{- c W(z)}=\sum_{n=0}^\infty \frac {c (n+c)^{n-1}}{(n)!} (-z)^n\;. \label{final-Lambert}\end{equation}
A detailed derivation of this expression is presented in Appendix A.
Taking $c=-4iM\omega$ and $z=\exp\left(\frac{v_H-u_s}{4M}\right)$ in ~\eqref{final-Lambert}, one can see that~\eqref{f(u)Lambert} and ~\eqref{finu} are equivalent.

\subsection{Delta function potential}
\label{sec:delta-function}

In this section, we apply our matching method to the case where the potential term in ~\eqref{chi-eq-outside} is replaced by
\begin{equation} V = \lambda \delta(r_*)  \;, \label{V-delta} \end{equation}
with $\lambda$ a positive real constant.  This can serve as a model for the original potential which has a single peak and vanishes at the horizon and infinity.
The resulting mode equation can be solved analytically and the solutions are simple enough that the matching coefficients can be computed analytically. Some of these matching coefficients will be used to partially reconstruct the mode functions $f^{\rm in}_{\omega \ell}$ in the case that $\ell = 0$.



For $\ell = 0$ in 4D the {\it in} modes inside the null shell take on the particularly simple form~\eqref{psi-in-0}.
In the region outside the shell the mode functions in the complete set with $\ell = 0$  have the general form
\begin{equation} f^{(H^{+}, \mathscr{I}^{+})}_{\omega' 0 0} = \frac{Y_{00}}{r \sqrt{4 \pi \omega'}} \psi^{(H^{+}, \mathscr{I}^{+})}_{\omega'0} \;, \qquad \psi^{(H^{+}, \mathscr{I}^{+})}_{\omega'0} = e^{-i \omega' t_s} \chi^{(H^{+}, \mathscr{I}^{+})}_{\omega'0} \;.  \label{f-H-I-delta} \end{equation}
The radial parts of the modes satisfy the following equation
\begin{equation}
    \frac{d^{2}\chi}{dr_{*}^2}+(\omega^2-\lambda\delta(r_{*}))\chi=0 \;.
\label{DiracDeltaModeEquation} \end{equation}
In the region where $r_*>0$, two linearly independent solutions are
\begin{subequations} \begin{eqnarray}
\chi_R &=& e^{i\omega r_{*}}\;, \label{RightMoving1}\\
\chi_L &=& e^{-i\omega r_{*}}\;. \label{LeftMoving1}
\end{eqnarray} \label{modes-delta-rs-gt-0} \end{subequations}
 For $r_*<0$, $\chi_R$ and $\chi_L$ can be expressed in the following way
\begin{subequations} \begin{eqnarray}
  \chi_{R}&=& E_R e^{i\omega r_{*}}+ F_R e^{-i\omega r_{*}}\;,  \label{RightMoving2}\\
  \chi_{L}&=& E_L e^{i\omega r_{*}}+ F_L e^{-i\omega r_{*}}\;,  \label{LeftMoving2}
\end{eqnarray} \label{modes-delta-rs-lt-0}\end{subequations}
where $E_R$ and $F_L$ are reflection and $F_R$ and $E_L$ are transmission coefficients.
Imposing the continuity of the mode function and discontinuity of its first derivative in the usual way at the spacelike curve $r_{*}=0$, the following analytic expressions are found for the scattering coefficients
\begin{subequations} \begin{eqnarray}
E_R &=& 1+\frac{i\lambda}{2\omega}\; , \label{Dirac-ER}\\
F_R &=& -\frac{i\lambda}{2\omega}\; ,\label{Dirac-FR}\\
E_L &=& F_R ^{*}=\frac{i\lambda}{2\omega}\; ,\label{Dirac-EL}\\
F_L &=& E_R ^{*}=1-\frac{i\lambda}{2\omega}\; .\label{Dirac-FL}
\end{eqnarray} \label{Scattering-Coefficients} \end{subequations}

Then the mode functions that we are using for the matching are
\begin{subequations} \begin{eqnarray} \psi^{H^{+}}_{\omega'0}  &=& \theta(-r_*) \left[ e^{-i \omega' v} + \frac{\frac{i \lambda}{2}}{\left(\omega' - \frac{i \lambda}{2} \right)} e^{-i \omega' u_s} \right]
  + \theta(r_*) \frac{\omega'}{\left(\omega' - \frac{i \lambda}{2} \right)} e^{-i \omega' v} \;, \label{psi-H-delta} \\
  \psi^{\mathscr{I}^{+}}_{\omega'0} &=& \theta(-r_*) \frac{\omega'}{\left(\omega' - \frac{i \lambda}{2} \right)} e^{-i \omega' u_s}  + \theta(r_*) \left[  \frac{\frac{i \lambda}{2}}{\left(\omega' - \frac{i \lambda}{2} \right)} e^{-i \omega' v} + e^{-i \omega' u_s} \right] \;. \label{psi-I-delta} \end{eqnarray} \label{psi-delta} \end{subequations}

To verify that the matching coefficients can be used to reconstruct the original mode functions for the case $\ell = 0$ it is useful to break them up into
contributions that come from the term proportional to  $e^{-i \omega v}$ in~\eqref{psi-in-0} and the term proportional to $-e^{-i \omega u}$.  In what follows we compute
the matching coefficients for both terms but then focus only those that come from the term proportional to $e^{-i \omega v}$.
Substituting~\eqref{psi-delta}, and~\eqref{Scattering-Coefficients} into~\eqref{gen-mat} one finds the matching coefficients
\begin{subequations} \begin{eqnarray}
A^{H^{+}}_{\omega \omega' 0} &=& (A^{H^{+}}_{\omega \omega' 0})_v  + (A^{H^{+}}_{\omega \omega' 0})_u  \;, \nonumber \\
(A^{H^{+}}_{\omega \omega' 0})_v &=& - \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{e^{i \omega' v_0}}{\omega' - i \epsilon} e^{-i \omega v_0}
   + \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{1}{F_L^{*}(\omega',\ell)}\frac{e^{i(\omega'-\omega)v_0}}{\omega'-\omega+i \epsilon} \;, \label{AHv-delta} \\
 (A^{H^{+}}_{\omega \omega' 0})_u &=&   \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{e^{i \omega' v_0}}{\omega' - i \epsilon} e^{-i \omega v_H}
    - \frac{1}{2 \pi} \sqrt{\frac{\omega}{ \omega'} } \int_{-\infty}^{v_H} du \,  e^{-i w u}  \left[ \theta(r_{*}) \frac{\omega'}{\omega' + \frac{i \lambda}{2}} e^{i \omega' v_0} \right.  \nonumber \\
 & &  \left. \qquad    + \, \theta(-r_{*}) \left( e^{i \omega' v_0}
     - \frac{\frac{i \lambda}{2}}{\omega' + \frac{i \lambda}{2}} e^{i \omega' u_s(u)} \right) \right] \;, \label{AHu-delta}
\end{eqnarray} \end{subequations}
\begin{subequations} \begin{eqnarray}
A^{\mathscr{I}^{+}}_{\omega \omega' 0} &=& (A^{\mathscr{I}^{+}}_{\omega \omega' 0})_v  + (A^{\mathscr{I}^{+}}_{\omega \omega' 0})_u  \;, \nonumber \\
(A^{\mathscr{I}^{+}}_{\omega \omega' 0})_v &=&  - \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{\frac{i \omega' \lambda}{2}}{\omega' + \frac{i \lambda}{2}}\frac{e^{-i(\omega-\omega')v_0}}{\omega'-\omega+i \epsilon}\;,  \label{AIv-delta}   \\
  (A^{\mathscr{I}^{+}}_{\omega \omega' 0})_u &=&  -\frac{1}{2 \pi} \sqrt{\frac{\omega'}{ \omega}} \int_{-\infty}^{v_H} du  \frac{du_s(u)}{du} \, e^{-i \omega u} e^{i \omega' u_s(u)} \left[ \theta(r_{*})  +
  \theta(-r_{*}) \frac{\omega'}{\omega'-\frac{i \lambda}{2}}  \right]
   \;.
  \label{AIu-delta} \end{eqnarray} \end{subequations}
Note that the relations~\eqref{B-w-minus-w}  are satisfied by these matching coefficients so the relation~\eqref{2D-in-mode-expansion-2} also holds. Thus
 \begin{eqnarray}
 (f^{in}_{\omega 0 0})_v &=&  \int_{-\infty}^\infty d \omega'   \left[ (A^{H^{+}}_{\omega \omega' 0})_v  f^{H^+}_{\omega'00} +
    (A^{\mathscr{I}^{+}}_{\omega \omega' 0})_v  f^{\mathscr{I}^+}_{\omega'00} \right] \;.
    \label{in-modes-v-part}
    \end{eqnarray}
Substituting~\eqref{AHv-delta},~\eqref{AIv-delta}, and~\eqref{psi-delta} into~\eqref{in-modes-v-part} gives after some algebra
\begin{subequations} \begin{eqnarray}  (f^{in}_{\w00})_v &=& \frac{Y_{00}}{r \sqrt{4 \pi \omega}} \left[\theta(-r_*) I_1 + \theta(r_*) I_2 \right] \;, \nonumber \\
   I_1 &=& - \frac{i}{2 \pi } e^{-i \omega v_0} \int_{-\infty}^\infty d \omega' \left[ \frac{\frac{i \lambda}{2} e^{i \omega'(v_0-u_s)}}{\left(\omega' - \frac{i \lambda}{2} \right)}
     + \frac{e^{i \omega'(v_0-v)}}{(\omega'-i \epsilon)} - \frac{\omega' e^{i \omega' (v_0-v)}}{\left(\omega' + \frac{i \lambda}{2} \right) \left( \omega' - \omega + i \epsilon \right)} \right] \nonumber \\
     &=&  \theta(v_0-u_s) e^{-i \omega v_0} \left[e^{-\frac{\lambda}{2}(v_0-u_s)} - 1 \right] + \theta(v_0-v) e^{-i \omega v_0} \nonumber \\
   & &     + \frac{ \theta(v-v_0)}{\left( \omega + \frac{i \lambda}{2} \right)} \left[ \frac{i \lambda}{2} e^{-i \omega v_0} e^{- \frac{\lambda}{2}(v-v_0)} + \omega e^{- i \omega v} \right]\;, \\
   I_2 &=& - \frac{i}{2 \pi } e^{-i \omega v_0} \int_{-\infty}^\infty d \omega' \left[ \frac{\omega' e^{i \omega'(v_0-v)}}{(\omega'-i \epsilon)\left(\omega' - \frac{i \lambda}{2} \right)}
   -  \frac{ e^{i \omega'(v_0-v)}}{\omega'-\omega + i \epsilon} + \frac{\frac{i \lambda}{2} e^{i \omega'(v_0-u_s)}}{(\omega'-i\epsilon)\left(\omega' - \frac{i \lambda}{2} \right)} \right] \nonumber \\
   &=&   \theta(v_0-v) e^{-i \omega v_0} e^{-\frac{\lambda}{2}(v_0-v)} + \theta(v-v_0) e^{-i \omega v}  \nonumber \\
    & & - \theta(u_s-v_0) \frac{ \frac{i \lambda}{2}}{\left( \omega + \frac{i \lambda}{2} \right)} \left[
   - e^{-i \omega v_0} e^{-\frac{\lambda}{2}(u_s-v_0)} + e^{-i \omega u_s} \right]\;.
 \end{eqnarray} \label{f-in-v}\end{subequations}
It is easy to verify that~\eqref{f-in-v} gives the correct values for $(f^{in}_{\w00})_v$ on the future horizon for $v \le v_0$, on the null shell surface $v = v_0$, and  on past null infinity for $v \ge v_0$.


\subsection{Partial analytic results for the matching coefficients in 4D for $\ell = 0$}

Because of the simple form of the {\it in} modes for $\ell = 0$ inside the null shell~\eqref{psi-in-0}, it is possible to compute the matching coefficients
for the $e^{-i \omega v}$ part analytically.  To do so we begin by
substituting~\eqref{psi-in-0} into~\eqref{gen-mat} with the result
\begin{subequations} \begin{eqnarray} A^{ H^{+}}_{\omega \omega' 0} &=& - \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{e^{i \omega' v_0}}{\omega' - i \epsilon} (e^{-i \omega v_0} - e^{-i \omega v_H})
   + \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{1}{F_L^{*}(\omega',0)}\frac{e^{i(\omega'-\omega)v_0}}{\omega'-\omega+i \epsilon}  \nonumber \\
   & & - \frac{1}{2 \pi} \sqrt{\frac{\omega }{\omega'}}  \int_{-\infty}^{v_H} du \,  \omega e^{-i \omega u}  \psi^{H^{+} *}_{\omega' \ell}(u,v_0) \;, \label{A-H-ell-0-mat}\\
      B^{ H^{+}}_{\omega \omega' 0} &=& \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{e^{-i \omega' v_0}}{\omega' + i \epsilon} (e^{-i \omega v_0} - e^{-i \omega v_H})
      - \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{1}{F_L(\omega',0)}\frac{e^{-i(\omega+\omega')v_0}}{\omega'+\omega-i \epsilon} \nonumber \\
      & & + \frac{1}{2 \pi} \sqrt{\frac{\omega}{ \omega'}}  \int_{-\infty}^{v_H} du \,e^{-i \omega u} \psi^{H^{+}}_{\omega' \ell} (u,v_0)\;, \label{B-H-ell-0-mat}\\
     A^{ \mathscr{I}^{+}}_{\omega \omega' 0} &=&  - \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{F_R^{*}(\omega',0)}{F_L^{*}(\omega',0)}\frac{e^{-i(\omega-\omega')v_0}}{\omega'-\omega+i \epsilon}
       +\frac{i}{2 \pi} \sqrt{\frac{\omega }{\omega'}} \int_{-\infty}^{v_H} du \, e^{-i \omega u} \partial_u \psi^{\mathscr{I}^+  *}_{\omega \ell} \;,  \label{A-I-ell-0-mat} \\
     B^{ \mathscr{I}^{+}}_{\omega \omega' 0} &=&  \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{F_R(\omega',0)}{F_L(\omega',0)}\frac{e^{-i(\omega+\omega')v_0}}{\omega'+\omega-i \epsilon}
     - \frac{i}{2 \pi} \sqrt{\frac{\omega}{ \omega'}} \int_{-\infty}^{v_H} du \,e^{-i \omega u} \partial_u \psi^{\mathscr{I}^+}_{\omega \ell} \;.
     \label{B-I-ell-0-mat} \end{eqnarray} \label{ell-0-mat}\end{subequations}

Note that the integrals have to be computed numerically because the mode functions in Schwarzschild spacetime must be computed numerically.
However, because of the simple form that $\psi^{\rm in}_{\omega 0}$ takes it is possible to separate the matching coefficients into separate matching
coefficients for the part that goes like $e^{- i \omega v}$ inside the null shell and the part that goes like $e^{-i \omega u}$ there.  The matching coefficients
for $e^{-i \omega v}$ don't depend on the integrals.  In what follows we focus on these matching coefficients.  Examination of~\eqref{ell-0-mat} gives for these coefficients
\begin{subequations} \begin{eqnarray} \left(A^{ H^{+}}_{\omega \omega' 0}\right)_v &=& - \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{e^{i (\omega'-\omega) v_0}}{\omega' - i \epsilon}
   + \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{1}{F_L^{*}(\omega',0)}\frac{e^{i(\omega'-\omega)v_0}}{\omega'-\omega+i \epsilon}  \nonumber \\
   & & \;, \label{A-H-ell-0-v-mat}\\
      \left(B^{ H^{+}}_{\omega \omega' l}\right)_v &=& \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{e^{-i (\omega'+\omega) v_0}}{\omega' + i \epsilon}
   - \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{1}{F_L(\omega',0)}\frac{e^{-i(\omega+\omega')v_0}}{\omega'+\omega-i \epsilon} \nonumber \\
      & & \;, \label{B-H-ell-0-v-mat}\\
     \left(A^{ \mathscr{I}^{+}}_{\omega \omega' l}\right)_v &=&  - \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{F_R^{*}(\omega',0)}{F_L^{*}(\omega',0)}\frac{e^{-i(\omega-\omega')v_0}}{\omega'-\omega+i \epsilon}
      \;,  \label{A-I-ell-0-v-mat} \\
     \left(B^{ \mathscr{I}^{+}}_{\omega \omega' l}\right)_v &=&  \frac{i}{2 \pi} \sqrt{\frac{\omega'}{\omega}} \frac{F_R(\omega',0)}{F_L(\omega',0)}\frac{e^{-i(\omega+\omega')v_0}}{\omega'+\omega-i \epsilon}
     \;.
     \label{B-I-ell-0-v-mat} \end{eqnarray} \label{ell-0-v-mat}\end{subequations}

These matching coefficients can be used to reconstruct the part of the mode function which goes like $e^{-i \omega v}$ inside the shell by substituting the expressions into~\eqref{General-in-modes}.  To check them we shall compute the resulting integral on $H^{+}$.  Recall that we are working in
the exact Schwarzschild spacetime rather than the null shell spacetime when we do the matching.  The same applies to the reconstruction.  Thus the results for the reconstruction for which $v \ge v_0$ also apply to the null shell spacetime,  but the results for $v < v_0$ do not apply to the null shell spacetime.

Recall that the modes $f^{\mathscr{I}^{+}}$ vanish on $H^{+}$.
\begin{subequations} \begin{eqnarray} \left(f^{\rm in}_{\omega 0 0} \right)_v &=& \frac{Y_0^0}{r \sqrt{4 \pi \omega}} (I_1 + I_2)\;,  \nonumber \\
   I_1 &=&  - \frac{i}{2 \pi} e^{-i \omega v_0} \int_0^\infty d \omega' \,\left[ \frac{e^{i(v_0-v)\omega'}}{\omega'-i \epsilon} -  \frac{e^{-i(v_0-v)\omega'}}{\omega'+i \epsilon} \right]\;, \label{I1} \\
   I_2 &=& \frac{i}{2 \pi} e^{-i \omega v_0} \int_0^\infty d \omega' \, \left[\frac{1}{F_L^{*}(\omega',0)} \frac{e^{-i \omega'(v-v_0)}}{\omega'-\omega + i \epsilon}
   - \frac{1}{F_L(\omega',0)} \frac{e^{i \omega'(v-v_0)}}{\omega'+\omega - i \epsilon} \right] \;. \label{I2}
 \end{eqnarray} \label{fin-ell-0-v}\end{subequations}
If in the second term of $I_1$ a change of variables is made so that $\omega' \to - \omega'$ then one finds that
\begin{equation} I_1 =   - \frac{i}{2 \pi} e^{-i \omega v_0} \int_{-\infty}^\infty d \omega' \frac{e^{i(v_0-v)\omega'}}{\omega'-i \epsilon} = e^{-i \omega v_0}\, \theta(v_0-v) \;,\end{equation}
with $\theta$ the step function.  It is thus clear that the initial data on $H^{+}$ for $-\infty < v < v_0$ does not affect the mode functions on the
part of the future horizon for which $v_0 < v < \infty$.

It can be shown from the properties of the scattering coefficients given in~\cite{rigorous}, that $F_L(\omega')=F^{*}_L(-\omega')$. Using this identity and changing the variable of integration in the second integral in the same way as was done for $I_1$, one obtains
\begin{equation}
(f^{in}_{\omega 0 0})_v =\frac{e^{-i\omega v_0}}{\sqrt{4\pi\omega}}\theta(v_{0}-v)+\frac{i}{2\pi \sqrt{4\pi\omega}}\int_{-\infty}^{\infty} d\omega' \frac{e^{-i\omega'v}}{
F^*_{L}(\omega')}\frac{e^{iv_{0}(\omega'-\omega)}}{\omega'-\omega+i\epsilon} \;.\end{equation}
To compute this integral using complex integration techniques one must know the singularity structure of $\frac{1}{F_L^{*}}$ which is difficult since this scattering coefficient
must be computed numerically.  However, one can at least test whether it has one or more singularities in the complex plane by assuming it does not and computing the integral.
We'll call the result $f^{\rm test}$ because there is no guarantee that this method will give the correct answer.
The result of such an integration is
\[ f^{\rm test} =\frac{e^{-i\omega v_0}}{\sqrt{4\pi\omega}}\theta(v_{0}-v)+\frac{e^{-i\omega v}}{ F_{L}^{*}(\omega)\sqrt{4\pi\omega}}\theta(v-v_{0})\;.
\]
Here complex integration has been performed using a contour in the lower half of the complex plane.
It is obvious that at $v=v_{0}$, the continuity condition for $(f^{in}_{\omega 0 0})_v$ is not satisfied so $f^{\rm test} \ne (f^{in}_{\omega 0 0})_v$
which implies that $\frac{1}{F_{L}^{*}(\omega')}$ has one or more singularities in the complex plane.

Alternatively one work with $I_2$ in the form~\eqref{I2}, use the relation $(\omega'\mp\omega \pm i\epsilon)^{-1} = \mp i \pi \delta(\omega'\mp\omega) + (\omega' \mp w)^{-1}$, and compute the principle value parts of the integral numerically for $v > v_0$.    This has been done and the result is shown in Fig.~\ref{fig:I2}.  It is clear from the plots in this figure that on the future horizon $(f^{\rm in}_{\omega 0 0})_v$ is continuous at $v = v_0$.
\begin{figure}[h]
\centering
\includegraphics [totalheight=0.27\textheight]{RealinMode.pdf}
\includegraphics [totalheight=0.27\textheight]{ImaginaryinMode.pdf}
\caption{The real (left) and the imaginary (right) parts of  $\sqrt{\frac{4 \pi}{M}} (f^{\rm in}_{\omega 00})_v $ on the future horizon have been plotted. In both plots, $M\omega = 0.02$ and $\frac{v_0}{M}=3$. The plots clearly show that $(f^{\rm in}_{\omega 00})_v$ is continuous at $v = v_0$.}
\label{fig:I2}
\end{figure}



\section{Stress-Energy Tensor}

\subsection{Method in 4D}
\label{sec-4D-method}

For the massless minimally coupled scalar field the classical stress-energy tensor in a general curved spacetime is given in~\eqref{Tab-class} and a renormalization expression for $\langle T_{ab} \rangle$ is given in~\eqref{Tab-Unruh-sub}.
To compute $\langle T_{ab} \rangle$ using~\eqref{Tab-Unruh-sub} it is useful to begin with the points split and to write the stress-energy tensor in terms of derivatives of the
Hadamard Green's function
\begin{equation}  G^{(1)}(x,x') = \langle \{ \phi(x), \phi(x') \}\rangle  \;, \label{G1-def} \end{equation}
Then adopting the notation
\begin{equation} \Delta G^{(1)}(x,x') = \langle {\rm in}| \{ \phi(x), \phi(x') \} | {\rm in}\rangle - \langle {U}| \{ \phi(x), \phi(x') \} | {U}\rangle \;, \label{Delta-G-def} \end{equation}
with $|{\rm in} \rangle$ representing the {\it in} vacuum state and $|U \rangle$ the Unruh state.  The corresponding difference in the stress-energy tensors is then
\begin{equation}  \Delta \langle T_{ab} \rangle = \frac{1}{4} \lim_{x' \to x} \left[ \left( g_a^{c'} \Delta G^{(1)}_{;c'; b} + g_b^{c'}  \Delta G^{(1)}_{;a; c'} \right)  - \, g_{a b}\, g^{c d'} \Delta G^{(1)}_{;c ;d'} \right]  \;. \label{Delta-Tab-1} \end{equation}
Here the quantity $g_a^{b'}$ parallel transports a vector from $x'$ to $x$ and is called the bivector of parallel transport~\cite{christensen-76}.
To leading order when the point separation is small
\begin{equation} g_a^{b'} = g_a^b = \delta_a^b  \;. \label{gabp-leading}\end{equation}
The subleading orders all vanish in the limit $x' \to x$.
Since there are no ultraviolet divergences in the quantity $\Delta \langle T_{ab} \rangle$ one can use~\eqref{gabp-leading} in~\eqref{Delta-Tab-1} with the result
\begin{equation}  \Delta \langle T_{ab} \rangle =   \frac{1}{4}  \left[\lim_{x' \to x} \left( \Delta G^{(1)}_{;a' ;b} +   \Delta G^{(1)}_{;a ;b'} \right)  - g_{a b}\, g^{c d}  \lim_{x' \to x} \Delta G^{(1)}_{;c ;d'} \right]  \;, \label{Delta-Tab-2} \end{equation}
where a slight abuse of notation has been used for the implied sum over $d$ and $d'$ in the last term.  It is important to note that this expression is valid in both two and four dimensions.

Expanding the field in terms of modes as in~\eqref{Phi-expansion} one finds for the {\it in} modes
that
\begin{equation}   \langle 0 \, {\rm in} |  \{ \phi(x), \phi(x') \} |0 \, {\rm in} \rangle = \sum_{\ell = 0}^\infty \sum_{m = - \ell}^\ell \int_0^\infty d \omega \; [f^{\rm in}_{\omega \ell m}(x) (f^{\rm in}_{\omega \ell m}(x'))^{*} + f^{\rm in}_{\omega \ell m}(x') (f^{{\rm in}}_{\omega \ell m}(x))^{*}]  \;.  \label{G1-in} \end{equation}
The Unruh state in Schwarzschild spacetime consists of modes that are positive frequency with respect to the usual time coordinate on $\mathscr{I}^{-}$ along with modes that are positive frequency with respect to the Kruskal time coordinate on $H^{-}$ so that
\begin{eqnarray}   \langle U |  \{ \phi(x), \phi(x') \} |U \rangle &=& \sum_{\ell = 0}^\infty \sum_{m = - \ell}^\ell \left\{ \int_0^\infty d \omega_K \, [f^{K}_{\omega_K \ell m}(x) (f^{K}_{\omega_K \ell m}(x'))^{*} + f^{K}_{\omega_K \ell m}(x') (f^{{K}}_{\omega_K \ell m}(x))^{*}] \right.  \nonumber \\
  && \; \left.  + \,  \int_0^\infty d \omega \; [f^{\mathscr{I}^{-}}_{\omega \ell m}(x) (f^{\mathscr{I}^{-}}_{\omega \ell m}(x'))^{*} + f^{\mathscr{I}^{-}}_{\omega \ell m}(x') (f^{{\mathscr{I}^{-}}}_{\omega \ell m}(x))^{*}]  \right\}  \;.  \label{G1-U}  \end{eqnarray}

The next step is to find expansions for these two-point functions in terms of the complete set of modes $f^{(\mathscr{I}^{+}, H^{+})}$ that we are using.
For $ \langle 0 \, {\rm in} |  \{ \phi(x), \phi(x') \} |0 \, {\rm in} \rangle$ one can substitute~\eqref{General-in-modes} into~\eqref{G1-in}.  This results in integrals of the form
\begin{eqnarray}  &&\sum_{\ell = 0}^\infty \sum_{m = - \ell}^\ell \int_0^\infty d \omega \; \int_0^\infty d \omega_1  \; \int_0^\infty d\omega_2 \left\{ \left[  A^{(H^{+}, \mathscr{I}^{+})} f^{(H^{+}, \mathscr{I}^{+})} +  B^{(H^{+}, \mathscr{I}^{+})} (f^{(H^{+}, \mathscr{I}^{+})})^{*} \right] \right. \nonumber \\
&& \qquad \qquad \left.  \times \left[  (A^{(H^{+}, \mathscr{I}^{+})})^{*} (f^{(H^{+}, \mathscr{I}^{+})})^{*} +  (B^{(H^{+}, \mathscr{I}^{+})})^{*} f^{(H^{+}, \mathscr{I}^{+})} \right] \right\}
 \;, \label{triple-integral-form} \end{eqnarray}
where the subscripts on the matching coefficients and mode functions have been suppressed.
For $\langle U |  \{ \phi(x), \phi(x') \} |U \rangle$ one can first substitute~\eqref{fK} and~\eqref{alphaK-betaK} into~\eqref{G1-U} to obtain an expression in terms of
$f^{(H^{-},\mathscr{I}^{-})}$.  The relationship between these modes and $f^{(\mathscr{I}^{+}, H^{+})}$ can be deduced from~\eqref{chi-H-scri-def} and used to obtain an expression for $\langle U |  \{ \phi(x), \phi(x') \} |U \rangle$ that depends only on $f^{(\mathscr{I}^{+}, H^{+})}$.

\subsection{2D Example}
\label{sec:Tab-2D}

In this section the method discussed above is tested by using it to computing the stress-energy tensor for the scalar field in the corresponding 2D spacetime where the answer is known.  The computation will be done in the region $v > v_0$ outside the null shell and outside the horizon.
From~\eqref{f-in-H+-2D} it is clear that for $v > v_0$ the contribution from the $f^{H^{+}}_{\omega'}$ modes to $f^{\rm in}_\omega$ is
\begin{equation} \left(f^{\rm in}_\omega \right)_{H^{+}} = \frac{e^{-i \omega v}}{\sqrt{4 \pi \omega}}  = f^{H^{+}}_{\omega}  \;. \label{f-in-H+-2} \end{equation}
Thus
\begin{eqnarray} f^{\rm in}_\omega = f^{H^+}_\omega +  \int_0^\infty d\omega \; [A^{\mathscr{I}^{+}}_{\omega \omega'} f^{\mathscr{I}^{+}}_{\omega'}  + B^{\mathscr{I}^{+}}_{\omega \omega'} (f^{\mathscr{I}^{+}}_{\omega'})^{*}]\;, \label{fin-A-B}
\end{eqnarray}
with $A^{\mathscr{I}^{+}}_{\omega \omega'}$ given in~\eqref{Ascrp2D}. Using the relation $\Gamma(x) = \frac{\Gamma(1+x)}{x}$ one obtains the form used for the numerical computations
\begin{eqnarray}
A^{\mathscr{I}^{+}}_{\omega \omega'} &=& -\frac{1}{2 \pi} \sqrt{\frac{\omega}{ \omega'}} (4M)^{i 4M \omega'} e^{-i(\omega-\omega') v_H} \frac{\Gamma(1-i4M\omega')}{[-i(\omega-\omega')+\epsilon]^{1-i4M\omega'}}\;.  \label{AI2}
\end{eqnarray}
Then, using the relations~\eqref{B-w-minus-w} one finds
\begin{equation} B^{\mathscr{I}^{+}}_{\omega \omega'} =  \frac{1}{2 \pi} \sqrt{\frac{\omega}{ \omega'}} (4M)^{-i 4M \omega'} e^{-i(\omega+\omega') v_H} \frac{\Gamma(1+i4M\omega')}{[-i(\omega+\omega')+\epsilon]^{1+i4M\omega'}}  \;. \label{BI2}  \end{equation}
In what follows the superscript $\mathscr{I}^+$ on the matching coefficients A and B will be suppressed.

Next, with the aim of finding the components of the stress-energy tensor using~\eqref{Delta-Tab-2}, we construct the Hadamard form of Green's function which in 2D is
\begin{eqnarray}
 G^{(1)}(x,x') = \int_0^\infty d \omega \; [f^{in}_\omega(x) f^{in \;*}_{w}(x') + f^{in \;*}_\omega(x) f^{in}_{w}(x') ]\;.
  \label{Hadamard}
 \end{eqnarray}
Substituting ~\eqref{fin-A-B} into ~\eqref{Hadamard} gives
 \begin{eqnarray}
 G^{(1)}(x,x') &=& \int_0^\infty d \omega \; \Big\{ \Big[f^{H^{+}}_\omega(x) +  \int_0^\infty d\omega_1 \; [A_{\omega \omega_1} f^{\mathscr{I}^{+}}_{\omega_1}(x)  + B_{\omega \omega_1} f^{\mathscr{I}^{+}\;*}_{\omega_1}(x)]\Big].  \nonumber \\
 & & \left.   \Big[f^{H^{+} \;*}_\omega(x') +  \int_0^\infty d\omega_2 \; [A^{*}_{\omega \omega_2} f^{\mathscr{I}^{+}\;*}_{\omega_2}(x')  + B^{*}_{\omega \omega_2} f^{\mathscr{I}^{+}}_{\omega_2}(x')]\Big] \right.  \nonumber \\
 & & \left. + \Big[f^{H^{+}}_\omega(x') +  \int_0^\infty d\omega_1 \; [A_{\omega \omega_1} f^{\mathscr{I}^{+}}_{\omega_1}(x')  + B_{\omega \omega_1} f^{\mathscr{I}^{+}\;*}_{\omega_1})(x')]\Big] \right.  \nonumber \\
 & & \Big[f^{H^{+} \;*}_\omega(x) +  \int_0^\infty d\omega_2 \; [A^{*}_{\omega \omega_2} f^{\mathscr{I}^{+}\;*}_{\omega_2}(x)  + B^{*}_{\omega \omega_2} f^{\mathscr{I}^{+}}_{\omega_2}(x)]\Big] \Big\}\;. \label{G1-1} \end{eqnarray}
  Expanding the integrand of the integral over $\omega$ results in three types of expressions: an integral consisting of products of the modes $f^{H^+}_\omega$, which we  call $G_A$, another integral which includes cross products between the modes $f^{H^+}_\omega$ and  $f^{\mathscr{I}^+}_\omega$, which we call $G_B$, and finally an integral consisting of products of the modes  $f^{\mathscr{I}^+}_\omega$, which we called $G_C$\;.

 To renormalize we follow a procedure equivalent to that outlined in Sec.~\ref{sec-4D-method}.
 We begin by subtracting off the integrals with the integrand evaluated in the large $\omega$ limit.  When we add them back, we get contributions that are identical to those obtained for the Unruh state.  We are not quite subtracting off the Unruh modes because the large $\omega$ solutions have a dependence on $v_H$.  However, when the subtracted terms are added back and the integral over $\omega$ is computed, then factors of $\delta(\omega_1 - \omega_2)$ and $\delta(\omega_1 + \omega_2)$ are obtained.  Note that terms proportional to $\delta(\omega_1 + \omega_2)$ vanish.  For the ones that do not vanish, once one integrates over say $\omega_2$,  the dependence on $v_H$ vanishes.

 In Appendix~\ref{appendix-B} it is shown that when this method is applied to $ G^{(1)}(x,x')$, the $\Delta G_A$ term vanishes.  It is also shown that, while the $\Delta G_B$ term does not vanish, it does not contribute to the stress-energy tensor.  As a result, the only term that contributes to $\Delta \langle T_{ab} \rangle$ is $\Delta G_C(x,x')$ which has the form
\begin{eqnarray}
\Delta G_C(x,x')&=& \int_0^{\infty} d\omega_1 \int_0^{\infty} d\omega_2 \Big\{ [f_{\omega_1}^{\mathscr{I^+}}(x)f_{\omega_2}^{\mathscr{I^+}\; *}(x')+f_{\omega_1}^{\mathscr{I^+}}(x')f_{\omega_2}^{\mathscr{I^+}\;*}(x)]\Delta I_1 \nonumber \\ & &
+\; [f_{\omega_1}^{\mathscr{I^+}}(x)f_{\omega_2}^{\mathscr{I^+}}(x')+f_{\omega_1}^{\mathscr{I^+}}(x')f_{\omega_2}^{\mathscr{I^+}}(x)]\Delta I_2 \nonumber \\ & &
+\;[f_{\omega_1}^{\mathscr{I^+}\;*}(x)f_{\omega_2}^{\mathscr{I^+}\;*}(x')+f_{\omega_1}^{\mathscr{I^+}\;*}(x')f_{\omega_2}^{\mathscr{I^+}\;*}(x)]\Delta I_3 \nonumber \\ & &
+\;[f_{\omega_1}^{\mathscr{I^+}\;*}(x)f_{\omega_2}^{\mathscr{I^+}}(x')+f_{\omega_1}^{\mathscr{I^+}\; *}(x')f_{\omega_2}^{\mathscr{I^+}}(x)]\Delta I_4 \Big\} \;, \label{Del-GC}
\end{eqnarray}
with
\begin{subequations} \begin{eqnarray}
\Delta I_1 &=& \int _0^{\infty} d\omega \left\{A_{\omega \omega _1}A_{\omega \omega _2}^{*}-\mathcal{O}(A_{\omega \omega _1}A_{\omega \omega _2}^{*})\right\} \label{Del-I1}\;,\\
\Delta I_2 &=& \int _0^{\infty} d\omega \left\{ A_{\omega \omega _1}B_{\omega \omega _2}^{*}-\mathcal{O}(A_{\omega \omega _1}B_{\omega \omega _2}^{*})\right\} \label{Del-I2}\;,\\
\Delta I_3 &=& \int _0^{\infty} d\omega \left\{B_{\omega \omega _1}A_{\omega \omega _2}^{*}-\mathcal{O}(B_{\omega \omega _1}A_{\omega \omega _2}^{*})\right\}\label{Del-I3}\;,\\
\Delta I_4 &=& \int _0^{\infty} d\omega \left\{B_{\omega \omega _1}B_{\omega \omega _2}^{*}-\mathcal{O}(B_{\omega \omega _1}B_{\omega \omega _2}^{*})\right\}\; \label{Del-I4}.
\end{eqnarray}  \label{In}  \end{subequations}
Here $\mathcal{O}$ indicates the asymptotic behavior of the matching coefficients for $\omega \gg \omega_{1,2}$.

The integrals in~\eqref{In} can be computed analytically. Substituting the explicit expression for A from~\eqref{AI2} into~\eqref{Del-I1} gives
\begin{subequations} \begin{eqnarray}  \Delta I_{1} &=&  \frac{1}{4 \pi^2  \sqrt{ \omega_1 \omega_2}} (4 M)^{i 4 M (\omega_1-\omega_2)} e^{i v_H(\omega_1 - \omega_2)} \Gamma(1-i 4M\omega_1) \Gamma(1+i4M\omega_2)
  \Delta K_1 \;, \label{Del-I1-1} \\
  \Delta K_{1} &=&   \lim_{\Lambda \to \infty} (-i)^{i 4 M \omega_1} (i)^{-i 4 M \omega_2}\bigg\{ \left[ \int_{0}^\Lambda d\omega \frac{\omega}{ ( \omega-\omega_1 + i \epsilon_1)^{1-i4M\omega_1} ( \omega-\omega_2 - i \epsilon_2)^{1+i4M\omega_2}}  \right. \nonumber \\   && \left. \qquad
- \int_{1}^\Lambda d\omega \,\omega^{-1+i4M(\omega_1-\omega_2)} \right]- \int_0^1  d\omega \, \omega^{-1+i4M(\omega_1-\omega_2)} \bigg\} \;. \label{DelK1-1}
\end{eqnarray} \end{subequations}
First we compute the indefinite integrals and evaluate them at the limits. Since $\epsilon_1$ and $\epsilon_2$  go to $O^+$ at the end of the calculation, it is acceptable to add terms containing them to the exponents. The first indefinite integral is
 \begin{eqnarray} \Delta K_{1a} &=& (-i)^{i 4 M \omega_1} (i)^{-i 4 M \omega_2} \int_{0}^\Lambda d\omega \frac{\omega}{ ( \omega-\omega_1 + i \epsilon_1)^{1-i4M(\omega_1- i \epsilon_1)} ( \omega-\omega_2 - i \epsilon_2)^{1+i4M(\omega_2+i \epsilon_2)}}\nonumber \\
                        &=&  (-i)^{i 4 M \omega_1} (i)^{-i 4 M \omega_2} \left[-i \frac{(\omega-\omega_1+i \epsilon_1)^{-i4M(\omega_1- i \epsilon_1)} (\omega-\omega_2-i \epsilon_2)^{-i 4M(\omega_2+ i \epsilon_2)}}{4 M (\omega_1-\omega_2 - i \epsilon_1 - i \epsilon_2)} \right]_0^\Lambda  \nonumber \\
                        &=& -i (-i)^{i 4 M \omega_1} (i)^{-i 4 M \omega_2} \left[ \frac{(\Lambda-\omega_1)^{i4M(\omega_1)} (\Lambda-\omega_2)^{-i 4M(\omega_2)}}{4 M (\omega_1-\omega_2) - i \epsilon_1 - i \epsilon_2} \right.
                         \nonumber \\
                         & & \left. \; +  i \frac{(-\omega_1)^{i4M(\omega_1)} (-\omega_2)^{-i 4M(\omega_2)}}{4 M (\omega_1-\omega_2 - i \epsilon_1 - i \epsilon_2)} \right] \;.  \end{eqnarray}
Note that after evaluating the integral at the limits, $\epsilon_1$ and $\epsilon_2$ are set to zero in the exponents because they have no effect there.
Also, each term is a combination of a principle value and a term proportional to $\delta(\omega_1-\omega_2)$, thus
\begin{subequations} \begin{eqnarray} \Delta K_{1a} &=&   e^{2 \pi M (\omega_1+\omega_2)} \left[-i \frac{(\Lambda-\omega_1)^{i4M\omega_1} (\Lambda-\omega_2)^{-i 4M\omega_2}}{4 M (\omega_1-\omega_2)} + \frac{\pi}{4M} \delta(\omega_1-\omega_2) \right]
   \nonumber \\
 & &  +  e^{-2 \pi M (\omega_1+\omega_2)} \left[ i \frac{\omega_1^{i4M\omega_1} \omega_2^{-i 4M\omega_2}}{4 M (\omega_1-\omega_2)}  - \frac{\pi}{4M} \delta(\omega_1-\omega_2) \right] \;.  \label{K1a}
    \end{eqnarray}
Here we adopt the notation that the principle value of a term such as $\frac{1}{a \pm i \epsilon}$ is written as $\frac{1}{a}$.
   The second and third integrals in~\eqref{DelK1-1} are
  \begin{eqnarray} \Delta K_{1b} &=& -(-i)^{i 4 M \omega_1} (i)^{-i 4 M \omega_2} \int_{1}^\Lambda d\omega \omega^{-1+i4M(\omega_1-\omega_2)} \nonumber \\
     &=&      i \frac{e^{2 \pi M(\omega_1+\omega_2)}}{4 M(\omega_1-\omega_2)} \left[ \Lambda^{i 4M(\omega_1-\omega_2)} - 1\right]  \;, \label{Del-K1b} \\
           \Delta K_{1c} &=&  -(-i)^{i 4 M \omega_1} (i)^{-i 4 M \omega_2} \int_0^1  d\omega \, \omega^{-1+i4M(\omega_1-\omega_2)} \nonumber \\
              &=&   - e^{2 \pi M(\omega_1+\omega_2)} \int_{-\infty}^0 dz e^{[i 4M (\omega_1-\omega_2)+ \epsilon] z}  =  -  \frac{e^{2 \pi M(\omega_1+\omega_2)}}{i 4M (\omega_1-\omega_2)+ \epsilon}  =  \frac{i e^{2 \pi M(\omega_1+\omega_2)}}{4M (\omega_1-\omega_2) -i \epsilon} \nonumber \\  & & \qquad =   \frac{ie^{2 \pi M(\omega_1+\omega_2)}}{4M (\omega_1-\omega_2)} -e^{2 \pi M(\omega_1+\omega_2)} \frac{\pi}{4M} \delta(\omega_1-\omega_2) \;, \label{Del-K1c}
\end{eqnarray} \end{subequations}
where in the integral for $\Delta K_{1c}$ the change of variable $z = \log \omega$ has been made and an integrating factor $\epsilon$ has been inserted.
Combining these results, one finds
\begin{subequations} \begin{eqnarray} \Delta K_1 &=&    e^{-2 \pi M (\omega_1+\omega_2)} \left[ i \frac{\omega_1^{i4M\omega_1} \omega_2^{-i 4M\omega_2}}{4 M (\omega_1-\omega_2)}  - \frac{\pi}{4M} \delta(\omega_1-\omega_2) \right] \nonumber \\
   & =&  i  e^{-2 \pi M (\omega_1+\omega_2)} \frac{\omega_1^{i4M\omega_1} \omega_2^{-i 4M\omega_2}}{4 M (\omega_1-\omega_2 - i \epsilon)} \;. \label{Del-K1-1} \end{eqnarray}
Substituting ~\eqref{Del-K1-1} into~\eqref{Del-I1-1} gives
\begin{eqnarray} \Delta I_1 &=& \frac{i}{4 \pi^2  \sqrt{ \omega_1 \omega_2}} (4 M)^{i 4 M (\omega_1-\omega_2)} e^{i v_H(\omega_1 - \omega_2)}   \Gamma(1-i 4M\omega_1) \Gamma(1+i4M\omega_2) \nonumber \\
  & & \qquad \times \, e^{-2 \pi M( \omega_1+\omega_2)} \frac{w_1^{i4M\omega_1} \omega_2^{-i 4M\omega_2}}{4 M (\omega_1-\omega_2 - i \epsilon)} \;. \label{Del-I1-2}
       \end{eqnarray} \end{subequations}
Note that this is a finite contribution to $\Delta G_c$ because of the factor of $e^{-2 \pi M( \omega_1+\omega_2)}$.

Next consider  $\Delta I_4$ which is the other term with non-vanishing delta functions.
\begin{subequations} \begin{eqnarray}  \Delta I_{4} &=&  \frac{1}{4 \pi^2  \sqrt{ \omega_1 \omega_2}} (4 M)^{-i 4 M (\omega_1-\omega_2)} e^{-i v_H(\omega_1 - \omega_2)} \Gamma(1+i 4M\omega_1) \Gamma(1-i4M\omega_2)
  \Delta K_4 \;, \label{Del-I4-1} \\
  \Delta K_{4} &=&   \lim_{\Lambda \to \infty} (-i)^{-i 4 M \omega_1} (i)^{i 4 M \omega_2}\bigg\{ \left[ \int_{0}^\Lambda d\omega \frac{\omega}{ ( \omega+\omega_1 + i \epsilon_1)^{1+i4M\omega_1} ( \omega+\omega_2 - i \epsilon_2)^{1-i4M\omega_2}}  \right. \nonumber \\   && \left. \qquad
- \int_{1}^\Lambda d\omega \,\omega^{-1-i4M(\omega_1-\omega_2)} \right]- \int_0^1  d\omega \, \omega^{-1-i4M(\omega_1-\omega_2)} \bigg\} \;. \label{DelK4-1} \;
\end{eqnarray} \end{subequations}
The integrals in $ \Delta K_{4}$ can be computed analytically
\begin{eqnarray}  \Delta K_{4a} &=& e^{-2\pi M(\omega_1+\omega_2)} \int_{0}^{\Lambda} d\omega  \frac{\omega}{( \omega+\omega_1+i \epsilon_1)^{1+i4M(\omega_1+i \epsilon_1)} (\omega+\omega_2-i \epsilon_2)^{1-i4M(\omega_2-i \epsilon_2)}} \nonumber \\
              &=& e^{-2\pi M(\omega_1+\omega_2)} \frac{i}{4 M} \left[ \frac{(\Lambda+\omega_1)^{-i4M\omega_1} (\Lambda+\omega_2)^{i4M\omega_2}}{\omega_1-\omega_2+ i (\epsilon_1+\epsilon_2)}
               -  \frac{\omega_1^{-i4M\omega_1} \omega_2^{i4M\omega_2}}{\omega_1-\omega_2+ i (\epsilon_1+\epsilon_2)}   \right] \nonumber \\
              \Delta K_{4b}& =&  e^{-2\pi M(\omega_1+\omega_2)} \left[ -\frac{\Lambda^{-i 4 M  (\omega_1-\omega_2)}}{-i 4M(\omega_1-\omega_2)} + \frac{1}{-i 4M(\omega_1-\omega_2)} \right] \,, \nonumber \\
              \Delta K_{4c} &= &  - e^{-2\pi M(\omega_1+\omega_2)} \int_{-\infty}^0 dz \, e^{[-i 4M(\omega_1-\omega_2)+\epsilon]z} = -\frac{ e^{-2\pi M(\omega_1+\omega_2)}}{-i 4M(\omega_1-\omega_2)+\epsilon}
             \nonumber \\
              &=& -i \frac{e^{-2\pi M(\omega_1+\omega_2)}}{4 M (\omega_1-\omega_2) + i \epsilon}  = -i  \frac{e^{-2\pi M(\omega_1+\omega_2)}}{4 M (\omega_1-\omega_2)} - e^{-2\pi M(\omega_1+\omega_2)} \frac{\pi}{4M}   \delta(\omega_1-\omega_2)  \;.
                \end{eqnarray}
Both terms in $\Delta K_{4a}$ can be written in terms of their principle values added to a term proportional to $\delta(\omega_1-\omega_2)$. Combining these terms, the following expression for $\Delta K_4$ is obtained
\begin{subequations} \begin{eqnarray}  \Delta K_4 &=&   e^{-2\pi M(\omega_1+\omega_2)}\left[-i \frac{\omega_1^{-i4M\omega_1} \omega_2^{i4M\omega_2}}{4M(\omega_1-\omega_2)} - \frac{ \pi }{4M}   \delta(\omega_1-\omega_2)\right]
   \nonumber \\
     &=& - i e^{-2\pi M(\omega_1+\omega_2)} \; \frac{\omega_1^{-i4M\omega_1} \omega_2^{i4M\omega_2}}{4M(\omega_1-\omega_2 + i \epsilon)}  \;.
     \end{eqnarray}
Finally
\begin{eqnarray} \Delta I_4 &=& -\frac{i}{4 \pi^2  \sqrt{ \omega_1 \omega_2}} (4 M)^{-i 4 M (\omega_1-\omega_2)} e^{-i v_H(\omega_1 - \omega_2)}   \Gamma(1+i 4M\omega_1) \Gamma(1-i4M\omega_2) \nonumber \\
    & & \qquad \times\, e^{-2\pi M(\omega_1+\omega_2)} \; \frac{\omega_1^{-i4M\omega_1} \omega_2^{i4M\omega_2}}{4M(\omega_1-\omega_2 + i \epsilon)}  \;. \label{Del-I4-2} \end{eqnarray} \end{subequations}
Note that if we let $\omega_1 \leftrightarrow \omega_2$ in the expression~\eqref{Del-I1-2} for $\Delta I_1$, then we get $\Delta I_4$ in~\eqref{Del-I4-2}.
It is also true that if this switch is made in the entire contribution to the two-point function from $\Delta I_1$ then that is equal to the contribution from $\Delta I_4$.
Next consider $\Delta I_2$
\begin{subequations} \begin{eqnarray}  \Delta I_{2} &=&  -\frac{1}{4 \pi^2  \sqrt{ \omega_1 \omega_2}} (4 M)^{i 4 M (\omega_1+\omega_2)} e^{i v_H(\omega_1 + \omega_2)} \Gamma(1-i 4M\omega_1) \Gamma(1-i4M\omega_2)
  \Delta K_2 \;, \label{Del-I2-1} \\
  \Delta K_{2} &=&   \lim_{\Lambda \to \infty} (-i)^{i 4 M \omega_1} (i)^{i 4 M \omega_2}\bigg\{ \left[ \int_{0}^\Lambda d\omega \frac{\omega}{ ( \omega-\omega_1 + i \epsilon_1)^{1-i4M\omega_1} ( \omega+\omega_2 - i \epsilon_2)^{1-i4M\omega_2}}  \right. \nonumber \\   && \left. \qquad
- \int_{1}^\Lambda d\omega \,\omega^{-1+i4M(\omega_1+\omega_2)} \right]- \int_0^1  d\omega \, \omega^{-1+i4M(\omega_1+\omega_2)} \bigg\} \;, \label{DelK2-1}
\end{eqnarray} \end{subequations}
where the integrals in $ \Delta K_{2}$ can be computed analytically
\begin{subequations} \begin{eqnarray}  \Delta K_{2a} &=&  (-i)^{i 4 M \omega_1} (i)^{i 4 M \omega_2} \int_{0}^\Lambda d\omega \frac{\omega}{( \omega-\omega_1 + i \epsilon_1)^{1-i4M(\omega_1-i \epsilon_1)} (\omega+\omega_2-i \epsilon_2)^{1-i4M(\omega_2-i \epsilon_2)}}  \nonumber \\
   &=& -\frac{i}{4M} (-i)^{i 4 M \omega_1} (i)^{i 4 M \omega_2} \left\{ \frac{(\Lambda-\omega_1)^{i4M \omega_1} (\Lambda+\omega_2)^{i4M\omega_2}}{[\omega_1+\omega_2 - i(\epsilon_1 + \epsilon_2)]} \right. \nonumber \\  & & \left. -  \frac{(-\omega_1)^{i4M\omega_1} \omega_2^{i4M\omega_2-}}{[\omega_1+\omega_2 - i(\epsilon_1 + \epsilon_2)]} \right\} \nonumber \\
   &=& -\frac{i}{4M} e^{2 \pi M(\omega_1 -\omega_2)} \frac{(\Lambda-\omega_1)^{i4M \omega_1} (\Lambda+\omega_2)^{i4M\omega_2}}{[\omega_1+\omega_2 - i(\epsilon_1 + \epsilon_2)]}  \nonumber \\
     & &   + \frac{i}{4M} e^{-2 \pi M(\omega_1+\omega_2)} \frac{\omega_1^{i4M\omega_1} \omega_2^{i4M\omega_2}}{[\omega_1+\omega_2 - i(\epsilon_1 + \epsilon_2)]}\;,  \label{Del-K2a}  \\
  \Delta K_{2b} &=&  \frac{i}{4 M} e^{2 \pi M(\omega_1 -\omega_2)} \frac{\Lambda^{i 4M (\omega_1+\omega_2)}}{\omega_1+\omega_2}  - \frac{i}{4 M} e^{2 \pi M(\omega_1 -\omega_2)} \frac{1}{\omega_1+\omega_2} \;, \\
  \Delta K_{2c} &=&  - e^{2\pi M(\omega_1-\omega_2)} \int_{-\infty}^0 dz \, e^{[i 4M(\omega_1+\omega_2)+\epsilon]z} = -\frac{ e^{2\pi M(\omega_1-\omega_2)}}{i 4M(\omega_1+\omega_2)+\epsilon}
             \nonumber \\
              &=& i \frac{e^{2\pi M(\omega_1-\omega_2)}}{4 M (\omega_1+\omega_2) - i \epsilon}  = i  \frac{e^{2\pi M(\omega_1-\omega_2)}}{4 M (\omega_1+\omega_2)} - \frac{\pi}{4M}   \delta(\omega_1+\omega_2)  \;.
\end{eqnarray} \end{subequations}
Given that $\delta(\omega_1+\omega_2) = 0$ since the frequencies are all non-negative, one can set $\epsilon_1=\epsilon_2=0$. Then
\begin{subequations} \begin{equation} \Delta K_2 = i e^{-2 \pi M(\omega_1+\omega_2)} \frac{\omega_1^{i4M\omega_1} \omega_2^{i4M\omega_2}}{4M(\omega_1+\omega_2)} \;, \end{equation}
and
\begin{eqnarray} \Delta I_2 &=& -\frac{i}{4 \pi^2  \sqrt{ \omega_1 \omega_2}} (4 M)^{i 4 M (\omega_1+\omega_2)} e^{i v_H(\omega_1 + \omega_2)}  \Gamma(1-i 4M\omega_1) \Gamma(1-i4M\omega_2) \nonumber \\
   & & \times e^{-2 \pi M(\omega_1+\omega_2)} \frac{\omega_1^{i4M\omega_1} \omega_2^{i4M\omega_2}}{4M(\omega_1+\omega_2)}  \;. \label{Del-I2-2} \end{eqnarray}  \end{subequations}

Comparing $\Delta I_2$ in~\eqref{Del-I2} and $\Delta I_3$ in~\eqref{Del-I3}, one can immediately see that their contributions to the two-point function,~\eqref{Del-GC}, are the complex conjugate of each other if one also takes $\omega_1 \leftrightarrow \omega_2$ in the contribution from $\Delta I_2$
 \begin{eqnarray} \Delta I_3 &=&  (\Delta I_2)^{*}\nonumber\\
 &=&\frac{i}{4 \pi^2  \sqrt{ \omega_1 \omega_2}} (4 M)^{-i 4 M (\omega_1+\omega_2)} e^{-i v_H(\omega_1 + \omega_2)}  \Gamma(1+i 4M\omega_1) \Gamma(1+i4M\omega_2) \nonumber \\
   & & \times e^{-2 \pi M(\omega_1+\omega_2)} \frac{\omega_1^{-i4M\omega_1} \omega_2^{-i4M\omega_2}}{4M(\omega_1+\omega_2)} \;.  \label{Del-I3-2} \end{eqnarray}

Substituting~\eqref{Del-I1-2},~\eqref{Del-I4-2},~\eqref{Del-I2-2}, and~\eqref{Del-I3-2} into~\eqref{Del-GC} one finds
\begin{eqnarray}
\Delta G_C(x,x')&=&\mathfrak{R}\Bigg\{ \frac{i}{8\pi^3}\int_0^{\infty}\frac{d\omega_1}{\omega_1}\int_0^{\infty}\frac{d\omega_2}{\omega_2}e^{-2\pi M(\omega_1+\omega_2)} \nonumber \\ & &
\times \bigg\{ \Big[e^{-i\omega_1 u_s+i\omega_2 u'_s}+e^{-i\omega_1 u'_s+i\omega_2 u_s}\Big] \frac{(4M\omega_1e^{\frac{v_H}{4M}})^{4iM\omega_1}}{(4M\omega_2e^{\frac{v_H}{4M}})^{4iM\omega_2}}\frac{\Gamma(1-4iM\omega_1)\Gamma(1+4iM\omega_2)}{4M(\omega_1-\omega_2-i\epsilon)} \nonumber \\ & &
-\Big[e^{-i\omega_1 u_s-i\omega_2 u'_s}+e^{-i\omega_1 u'_s-i\omega_2 u_s}\Big](4M\omega_1e^{\frac{v_H}{4M}})^{4iM\omega_1}(4M\omega_2e^{\frac{v_H}{4M}})^{4iM\omega_2} \nonumber \\ & &
 \times \frac{\Gamma(1-4iM\omega_1)\Gamma(1-4iM\omega_2)}{4M(\omega_1+\omega_2)}\bigg\} \Bigg\}\;. \label{Delta-GC-fin}
 \end{eqnarray}
There are infrared divergences in this expression.  However, it is easy to see that the derivatives in the general formula for the stress-energy tensor~\eqref{Tab-class} bring down factors of $\omega_1$ and $\omega_2$ which remove these infrared divergences.  Recalling that $\Delta G_C$ is the only contribution to $\langle T_{ab} \rangle$ from  $\Delta G^{(1)}$, it is straight-forward to show using~\eqref{Delta-Tab-2},~\eqref{T_trB},~\eqref{rtpttp},~\eqref{DTrt-DTtt}, and~\eqref{Delta-GC-fin} that
\begin{eqnarray}
\Delta \langle T_{tt} \rangle&=&  -(1-\frac{2M}{r}) \lim_{x' \to x} \frac{1}{4}(\Delta G_{C\; ;t';r}+\Delta G_{C \;;t;r'})\nonumber \\
&=&\mathfrak{R}\Bigg\{ \frac{i}{8\pi^3}\int_0^{\infty}d\omega_1 \; \int_0^{\infty}d\omega_2 \; e^{-2\pi M(\omega_1+\omega_2)} \nonumber \\ & &
\times \bigg\{ e^{i(\omega_2-\omega_1) u_s} \frac{(4M\omega_1e^{\frac{v_H}{4M}})^{4iM\omega_1}}{(4M\omega_2e^{\frac{v_H}{4M}})^{4iM\omega_2}}\frac{\Gamma(1-4iM\omega_1)\Gamma(1+4iM\omega_2)}{4M(\omega_1-\omega_2-i\epsilon)} \nonumber \\ & &
+ e^{-i(\omega_2 +\omega_1)u_s}(4M\omega_1e^{\frac{v_H}{4M}})^{4iM\omega_1}(4M\omega_2e^{\frac{v_H}{4M}})^{4iM\omega_2} \nonumber \\ & &
 \times \frac{\Gamma(1-4iM\omega_1)\Gamma(1-4iM\omega_2)}{4M(\omega_1+\omega_2)}\bigg\} \Bigg\}\;.
\label{2D-Ttt} \end{eqnarray}

 The integral over $\omega_2$ of the first term inside the curly bracket can be written in the form
  \begin{eqnarray}
 \Delta \langle T_{tt} \rangle _1 &=& \int_0^{\infty} d\omega_2 \; \frac{f(\omega_2)}{\omega_1-\omega_2-i\epsilon}=  \int_0^\infty d\omega_2 \;\left[ \frac{f(\omega_2)}{\omega_1-\omega_2}+i\pi\delta(\omega_1-\omega_2) \right] \nonumber \\
 &=& \lim_{\epsilon \to 0^+} \left[ \int_0^{\omega_1-\epsilon} d \omega_2  \frac{f(\omega_2)}{\omega_1-\omega_2} + \int_{\omega_1+\epsilon}^\infty d \omega_2  \frac{f(\omega_2)}{\omega_1-\omega_2} \right] + i \pi f(\omega_1)
\;, \label{CPV}
  \end{eqnarray}
where the definition of the Cauchy principal value integral has been explicitly used.

Thus, extracting the explicit form of the $f(\omega_2)$  from ~\eqref{2D-Ttt} and substituting it into ~\eqref{CPV} yields
\begin{eqnarray}
\Delta \langle T_{tt} \rangle _1
 &=&\mathfrak{R}\Bigg\{ \frac{i}{8\pi^3}\int_0^{\infty} d\omega_1\; \Big[\int_0^{\infty} d\omega_2 \;  e^{-2\pi M(\omega_1+\omega_2)} e^{i(\omega_2- \omega_1) u_s} \nonumber \\ & & \times
 \frac{(4M\omega_1e^{\frac{v_H}{4M}})^{4iM\omega_1}}{(4M\omega_2e^{\frac{v_H}{4M}})^{4iM\omega_2}}\frac{\Gamma(1-4iM\omega_1)\Gamma(1+4iM\omega_2)}{4M(\omega_1-\omega_2)}\Big]\Bigg\}+
 \frac{1}{4\pi^2}\int_0^{\infty} d\omega_1 \; \frac{e^{-4\pi M\omega_1}}{4M} \nonumber \\ & &
\times \Gamma(1-4iM\omega_1)\Gamma(1+4iM\omega_1)\;.  \label{Del-Ttt}
 \end{eqnarray}

The stress-energy tensor for a massless minimally coupled scalar field in the 2D collapsing null shell spacetime has been previously computed analytically using a different method~\cite{hiscock}\cite{Fabbri:2005mw} and the stress-energy tensor for the Unruh state has also been computed analytically~\cite{Davies-Fulling-Unruh,Fabbri:2005mw}.  For the difference one finds
 \begin{eqnarray}
\Delta \langle T_{uu} \rangle &=& -\frac{1}{24\pi}\left[\frac{8M}{(u-v_0)^3}+\frac{24M^2}{(u-v_0)^4}\right] - \frac{1}{768 \pi M^2}  \;, \nonumber  \\
\Delta \langle T_{uv} \rangle &=& \Delta \langle T_{vv} \rangle = 0 \;, \nonumber \\
\Delta \langle T_{tt} \rangle &=& \Delta \langle T_{uu} \rangle + 2 \Delta \langle T_{uv} \rangle + \Delta \langle T_{vv} \rangle \nonumber \\
 & = & -\frac{1}{24\pi}\left[\frac{8M}{(u-v_0)^3}+\frac{24M^2}{(u-v_0)^4}\right]  - \frac{1}{768 \pi M^2}  \;.
\label{Hiscock-Result} \end{eqnarray}


Both terms in~\eqref{Del-Ttt} have been computed numerically. In the first integral, the numerical computation has been performed  by the symmetric removal of the neighborhood with radius $\epsilon$ about the singular points of the integrand, $\omega_1=\omega_2$. The integral of the second term in ~\eqref{2D-Ttt}  has also been found by a more straightforward numerical method.
 Our results for $\Delta \langle  T_{tt} \rangle$ in~\eqref{2D-Ttt} are shown in Figure~\ref{fig:Ttt}. Although it is not possible to detect this from the plot, our numerical results agree with the analytical results in~\cite{hiscock} \cite{Fabbri:2005mw} to more than ten digits.
 \begin{figure}[h]
\centering
\includegraphics [totalheight=0.45\textheight]{Stress_Tensor_2D_plot.pdf}
\caption{The quantity $10^{4}M^2\Delta \langle T_{tt} \rangle$ is plotted for the massless minimally coupled scalar field in the region exterior to the null shell
and to the event horizon. The dots correspond to the results of the numerical computations.  The solid curve represents the analytic results in~\eqref{Hiscock-Result}.   }
\label{fig:Ttt}
\end{figure}


It is worth mentioning that in 2D, once  $\Delta \langle T_{tt} \rangle$ is numerically computed,  $\Delta \langle T_{rr} \rangle$ and  $\Delta \langle T_{tr} \rangle$ can be easily determined from the relations~\eqref{T_rrB} and~\eqref{DTrt-DTtt}.

\section{Summary}

We have presented a method of numerically computing the stress-energy tensor for a massless minimally coupled scalar field in the case when a black hole that is formed from the
collapse of a spherically symmetric null shell in four dimensions.  There are two primary parts to the method.  The first is to expand the mode functions in the natural {\it in} vacuum state in terms of a complete set of mode functions in the part of Schwarzschild spacetime that is outside of the event horizon of the black hole. Expressions have been found for the matching coefficients that involve integrals of these mode functions over the trajectory of the null shell.

The second part of the method involves subtracting the unrenormalized expression for the stress-energy tensor in the Unruh state from the expression for the unrenormalized stress-energy tensor in the {\it in} vacuum state.  Since the ultraviolet divergences in the stress-energy tensor are independent of the state, this difference is finite.  Then one can add to this the renormalized expression for the stress-energy tensor in the Unruh state that has already been computed~\cite{levi-et-al-kerr} and the result is the full renormalized stress-energy tensor for the {\it in} vacuum state.

We have tested the first part of the method by  analytically  computing the matching coefficients in the 2-D case and reconstructing the mode functions for the {\it in} vacuum state.  We have also analytically computed the matching coefficients in 4D for the spherically symmetric mode functions (those with $\ell = 0$) in the {\it in} vacuum state for a simple model in which the effective potential in the mode equation is proportional to a Dirac delta function.  In this case it was possible to analytically compute the part of the mode function in the {\it in} vacuum state that is proportional to $e^{-i \omega v}$ inside the null shell and to verify that it gives the known result on the matching surface.  Finally, for the actual case of a collapsing null shell in 4D, we have analytically computed parts of the matching coefficients and used those parts to numerically compute part of one of the {\it in} modes on the future horizon and shown that it has the correct value at the point where the future horizon intersects the null shell trajectory.
 We also used the method to numerically compute the stress-energy tensor for the 2-D case and found that the result was in excellent agreement with the known analytic result.

These tests provide substantial evidence that the method will work and that it will be possible to numerically compute the exact renormalized stress-energy tensor for a massless minimally coupled scalar field in a 4D spacetime in which a black hole forms from the collapse of a spherically symmetric null shell.  Work on that computation is in progress.

\acknowledgments

P. R. A. would like to thank Eric Carlson, Charles Evans, Adam Levi, and Amos Ori for helpful conversations and Adam Levi for sharing some of his numerical data.
A.F. acknowledges partial financial support from the
Spanish Ministerio de Ciencia e Innovaci\'on grant
FIS2017-84440-C2-1-P and from the Generalitat Valenciana grant
PROMETEO/2020/079.  This work was supported in part by the National Science Foundation under Grants No.  PHY-1308325, PHY-1505875, and PHY-1912584 to Wake Forest University.
Some of the numerical work was done using the WFU DEAC cluster; we thank the WFU Provost's Office and Information
Systems Department for their generous support.

