# Numerical Analysis of a Hybrid Method for Radiation Transport††thanks: This
material is based upon work supported by the U.S. Department of Energy, Office
of Science, Office of Advanced Scientific Computing Research, as part of their
Applied Mathematics Research Program. The work was performed at the Oak Ridge
National Laboratory, which is managed by UT-Battelle, LLC under Contract No.
De-AC05-00OR22725. The United States Government retains and the publisher, by
accepting the article for publication, acknowledges that the United States
Government retains a non-exclusive, paid-up, irrevocable, world-wide license
to publish or reproduce the published form of this manuscript, or allow others
to do so, for the United States Government purposes. The Department of Energy
will provide public access to these results of federally sponsored research in
accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-
public-access-plan).

Cory, Victor  Andrés Galindo-Olarte Department of Mathematics, Michigan State
University, East Lansing, MI 48824, USA  Victor P. DeCaria Computer Science
and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831,
USA  Cory D. Hauck Computer Science and Mathematics Division, Oak Ridge
National Laboratory, Oak Ridge, TN 37831, USA and Mathematics Department,
University of Tennessee, Knoxville, TN 37996, USA

###### Abstract

In this work, we prove rigorous error estimates for a hybrid method introduced
in [15] for solving the time-dependent radiation transport equation (RTE). The
method relies on a splitting of the kinetic distribution function for the
radiation into uncollided and collided components. A high-resolution method
(in angle) is used to approximate the uncollided components and a low-
resolution method is used to approximate the the collided component. After
each time step, the kinetic distribution is reinitialized to be entirely
uncollided. For this analysis, we consider a mono-energetic problem on a
periodic domains, with constant material cross-sections of arbitrary size. To
focus the analysis, we assume the uncollided equation is solved exactly and
the collided part is approximated in angle via a spherical harmonic expansion
($\text{P}_{N}$ method). Using a non-standard set of semi-norms, we obtain
estimates of the form $C(\varepsilon,\sigma,\Delta t)N^{-s}$ where $s\geq 1$
denotes the regularity of the solution in angle, $\varepsilon$ and $\sigma$
are scattering parameters, $\Delta t$ is the time-step before
reinitialization, and $C$ is a complicated function of $\varepsilon$,
$\sigma$, and $\Delta t$. These estimates involve analysis of the multiscale
RTE that includes, but necessarily goes beyond, usual spectral analysis. We
also compute error estimates for the monolithic $\text{P}_{N}$ method with the
same resolution as the collided part in the hybrid. Our results highlight the
benefits of the hybrid approach over the monolithic discretization in both
highly scattering and streaming regimes.

##  1 Introduction

The radiation transport equation (RTE) [20, 19, 5] describes the movement of
particles through a material medium by means of a kinetic distribution
function that gives the density of particles with respect to the local phase
space measure. In a general setting, the phase space is six dimensional: three
dimensions for particle position and three for particle momentum, the latter
of which is typically decomposed into energy and direction (or angle) of
flight. Thus in the time-dependent setting, the RTE is defined over a seven-
dimensional domain.

The RTE describes two basic processes: particle advection and interactions
with the material medium. These interactions can be of various types and
include scattering and emission/absorption processes. The rate at which these
processes occur is determined by the properties of the material, expressed via
cross-sections. Material cross-sections may vary in space and depend on the
particle energy and, in situations that the material evolves, the cross-
sections may evolve as well. When cross-sections vary significantly, the RTE
may exhibit multiscale behavior. It is the combination of this multiscale
behavior with the high-dimensional phase space that makes simulating the RTE a
challenging task.

A well-known multi-scale feature of the RTE is the diffusion limit. In regions
where the scattering cross-section is large, the solution of the RTE can be
accurately approximated by its angular average [17, 4]. Moreover this average
is well-approximated by the solution of a diffusion equation. This solution to
diffusion equation has long been used as a cheap approximation the solution to
the RTE in scattering dominated regimes.

Another common limit is the absorption limit, which is characterized by a
complete lack of scattering. In this case, the RTE does not have a simple
asymptotic approximation. However, due the abscence of scattering, there is no
coupling between the angles and energies of the kinetic distribution. Thus,
with a proper discretiation, the RTE solution can be easily parallelized.

In problems for which the scattering cross-section varies dramatically, both
of the limits above can exist simultaneously, along with a range of transition
regimes in between. A consequence of this fact is that a monolithic numerical
treat of the RTE will require many degrees of freedom that are strongly
coupled. In practice, the time-dependent RTE is often updated in time with an
implicit scheme. In such cases, designing the linear solvers can be a
challenge.

A variety of approaches have been proposed for addressing the multiscale
challenges posed of the RTE. These include micro-macro decompositions [18],
high order-low order (HOLO) methods [6], diffusion-based acceleration [19, 1],
and preconditioned Krylov approaches [21]. In the current paper, we consider a
hybrid formulation [15] that is based on the notion of first-collision source
[2]. In this hybrid formulation, the RTE is split into two components: an
uncollided component that tracks the particles up to point of their first
material interaction and a collided component that track the particles that
remain. The resulting system is then approximated with two different angular
discretizations: a high-resolution discretization for the uncollided equation
and a low-resolution for the collided equation. The intuition that drives this
strategy is that scattering produces a smoother solution; hence the collided
equation should require less resolution to recover an accurate solution. The
uncollided equation, on the other, requires higher resolution; however it
takes the form of a purely absorbing RTE and can therefore be solved much more
efficiently the original RTE using the same number of degrees of freedom. The
efficiency of the hybrid approach for the RTE has been demonstrated in several
papers [8, 9, 10], including generalizations to hybrid energy discretizations
[22] and hybrid spatial discretizations [16].

A key component of the hybrid implementation for the time-dependent RTE is a
relabeling procedure that, after a given time step, maps the collided
numerical solution into the space of the uncollided numerical solution and
then uses the sum to re-initialize the uncollided equation. Meanwhile, the
collided equation is re-initialized to zero. This relabeling step is critical,
since otherwise the hybrid numerical solution would eventually convergence to
a low-resolution numerical solution of the collided equation.

Despite the intuitive motivation of the hybrid and the success of the hybrid
approach in numerical simulations, the method still lacks rigorous
justification. This is due in part to complications introduced by the
multiscale behavior of the RTE. For example, spectral approximations of the
RTE in angle are fairly straightforward to analyze [13], but a multiscale
analysis that takes into account the degree of scattering is significantly
more complicated [7]. The relabeling step of the hybrid formulation
complicates the situation even further.

In the current paper, we take a first step in analyzing the hybrid method for
the time-dependent mono-energetic version of the RTE with isotropic
scattering. We focus only on the angular discretization of the RTE, comparing
the standard spectral approximation ($\text{P}_{N}$) for the full system with
a discretization of the hybrid that features a spectral approximation of the
same resolution for the collided equation but assumes an exact solution for
the uncollided equation. Clearly, the hybrid formulated in this way is more
expensive than the monolithic approach. Thus the goal of the analysis is
determine what is gained from the extra work involves in a high-resolution
simulation of the uncollided equation, which in practice is computed with a
high-fidelity collocation method or with a Monte-Carlo method.

The remainder of this paper is organized as follows. In Section 2, we
introduce the RTE, reduce it to the purely scattering problem, recall the
$\text{P}_{N}$ method, and then describe the setup of the hybrid. Having
established the setting of of the problem, we then summarize the main results
of the paper. In Section 3, we derive error estimates for the $\text{P}_{N}$
equations. In Section 4, we analyze the hybrid problem. In Section 5, we
generalize results back to the original RTE with non-zero absorption. In
Section 6, we provide a short conclusion. The appendix contains some generic
results used for the estimates of the main paper.

##  2 Background

###  2.1 The radiation transport equation

We consider a time-dependent transport equation with periodic boundaries,
isotropic scattering, unit-speed particles, and diffusion scaling:

$\displaystyle\varepsilon\partial_{t}\Psi^{\varepsilon}+\Omega\cdot\nabla_{x}\Psi^{\varepsilon}+\frac{\sigma_{\rm{t}}}{\varepsilon}\Psi^{\varepsilon}=\left(\frac{\sigma_{\rm{t}}}{\varepsilon}-\varepsilon\sigma_{\rm{a}}\right)\overline{\Psi^{\varepsilon}}+\varepsilon
Q,\qquad\overline{\Psi^{\varepsilon}}=\frac{1}{4\pi}\int_{\mathbb{S}^{2}}\Psi^{\varepsilon}\,d\Omega,$
(1a) $\displaystyle\left.\Psi^{\varepsilon}\right|_{t=0}=g.$ (1b)

Here $\Psi^{\varepsilon}=\Psi^{\varepsilon}(x,\Omega,t)$ is a function of
position $x\in X=[0,2\pi)^{3}$, direction of flight $\Omega\in\mathbb{S}^{2}$,
and time $t>0$. It can be interpreted physically as the density of particles
at time $t$ with respect to the measure $d\Omega dx$. Particles interact with
a material background characterized by an absorption cross-section
$\sigma_{\rm{a}}\geq 0$, total cross-section
$\sigma_{\rm{t}}\geq\sigma_{\rm{a}}$ (which accounts for scattering and
absorption), and a known source $Q=Q(x,\Omega,t)$. The quantity
$\sigma_{\rm{t}}-\varepsilon^{2}\sigma_{\rm{a}}$ is the scaled scattering
cross-section, where the non-dimensional parameter $\varepsilon>0$
characterizes the strength of the scattering as well as the relevant time
scale. Indeed, it is well-known [17] that in the limit $\varepsilon\to 0$,
$\Psi^{\varepsilon}\to\Psi^{0}$ where $\Psi^{0}$ is independent of angle and
satisfies the diffusion equation of the form

$\displaystyle\partial_{t}\Psi^{0}-\nabla_{x}\cdot\left(\frac{1}{3\sigma_{\rm{t}}}\nabla_{x}\Psi^{0}\right)+\sigma_{\rm{a}}\Psi^{0}=0,\qquad\left.\Psi^{0}\right|_{t=0}=\frac{1}{4\pi}\int_{\mathbb{S}^{2}}g\,d\Omega.$
(2a)

For $g\in L^{2}(X\times\mathbb{S}^{2})$, $\sigma_{\rm{t}},\sigma_{\rm{a}}\in
L^{\infty}(X)$, $Q\in L^{2}(X\times\mathbb{S}^{2}\times[0,\infty))$ (1), is
known to have a semi-group solution $\Psi^{\varepsilon}\in
C([0,\infty);L^{2}(X\times\mathbb{S}^{2}))$ [12, Theorem XXI.2.3]. If in
addition, $\Omega\cdot\nabla_{x}g\in L^{2}(X\times\mathbb{S}^{2})$, then
$\Psi^{\varepsilon}\in C^{1}([0,\infty);L^{2}(X\times\mathbb{S}^{2}))$. We
assume this is the case for remainder of the paper.

In order to facilitate a clear stability and error analysis, we assume that
the cross-sections $\sigma_{\rm{a}}$ and $\sigma_{\rm{t}}$ are constant in
space. This assumption on $\sigma_{\rm{a}}$ allows us to convert (1) to a
purely scattering system for the function
$\psi=e^{\sigma_{\rm{a}}t}\Psi^{\varepsilon}$:111To reduce notation, we
suppress the dependence of $\psi$ on $\varepsilon$.

$\displaystyle\varepsilon\partial_{t}\psi+\Omega\cdot\nabla_{x}\psi+\frac{\sigma}{\varepsilon}\psi=\frac{\sigma}{\varepsilon}\overline{\psi}+\varepsilon
q,\qquad\overline{\psi}=\frac{1}{4\pi}\int_{\mathbb{S}^{2}}\psi\,d\Omega,$
(3a) $\displaystyle\left.\psi\right|_{t=0}=g,$ (3b)

where $q=e^{\sigma_{\rm{a}}t}Q$ and $\sigma:=\sigma_{\rm{t}}$. Henceforth, we
focus our analysis on (3). The results can then be translated back to the case
of non-zero absorption by undoing the transformation, which gives exponential
decay if $\sigma_{\rm{a}}>0$. This assumption is made for simplicity, but it
does introduce a measure of regularity into the solution that is not typical
in applications. Indeed, a more reasonable assumption is that the cross-
sections are piece-wise smooth and that the boundaries are equipped with
inflow data. Hence the analysis here can be viewed as a localized proxy for a
more realistic scenario. A more sophisticated analysis to include boundary and
interior layers is the subject of future work.

###  2.2 The $\text{P}_{N}$ approximation

Given $N\in\mathbb{N}_{\geq 0}$, the $\text{P}_{N}$ method is a spectral
discretization of the transport equation with respect to the angular variable
$\Omega$. Let $\\{m_{\ell,k}\\}_{\ell,k}$ be the real-valued, orthonormal
basis of spherical harmonics, where $\ell\geq 0$ denotes the degree and
$k\in\\{-\ell,\dots,\ell\\}$ denotes the order. For any $u\in
L^{2}(\mathbb{S}^{2})$, the angular moment $u_{\ell,k}$ is given by

$u_{\ell,k}=\int_{\mathbb{S}^{2}}m_{\ell,k}u\,d\Omega.$ (4)

For convenience, we collect the basis elements of degree $\ell$ into vectors
$\mathbf{m}_{{\ell}}=(m_{\ell,-\ell},\dots,m_{\ell,\ell})^{\mathsf{T}}$, and
we denote by
$\mathbf{u}_{\ell}=(u_{\ell,-\ell},\dots,u_{\ell,\ell})^{\mathsf{T}}$ the
vector of corresponding moments. Let $\mathbb{P}_{N}(\mathbb{S}^{2})\subset
L^{2}(\mathbb{S}^{2})$ to be the span of all spherical harmonics with degree
at most $N$. Then the orthogonal projections $\mathcal{P}_{N}\colon
L^{2}(\mathbb{S}^{2})\rightarrow\mathbb{P}_{N}(\mathbb{S}^{2})$ and
$\widetilde{\mathcal{P}}_{N}\colon
L^{2}(\mathbb{S}^{2})\rightarrow\mathbb{P}_{N}(\mathbb{S}^{2})$ are given by

$\mathcal{P}_{N}u=\sum_{\ell=0}^{N}\mathbf{m}_{{\ell}}^{\mathsf{T}}\mathbf{u}_{\ell}=\sum_{\ell=0}^{N}\sum_{k=-\ell}^{\ell}m_{\ell,k}u_{\ell,k}\quad\text{and}\quad\widetilde{\mathcal{P}}_{N}u=(\mathcal{I}-\mathcal{P}_{N})u=\sum_{\ell=N+1}^{\infty}\mathbf{m}_{{\ell}}^{\mathsf{T}}\mathbf{u}_{\ell}=\sum_{\ell=N+1}^{\infty}\sum_{k=-\ell}^{\ell}m_{\ell,k},$
(5)

where $\mathcal{I}$ is the identity operator.

The $\text{P}_{N}$ approximation of (3) seeks a function

$\psi^{N}(x,\Omega,t)=\sum_{\ell=0}^{N}\mathbf{m}_{\ell}^{\mathsf{T}}(\Omega){\bm{\uppsi}}^{N}_{\ell}(x,t)=\sum_{\ell=0}^{N}\sum_{k=-\ell}^{\ell}m_{\ell,k}(\Omega){\bm{\uppsi}}^{N}_{\ell,k}(x,t)$
(6)

such that

$\displaystyle\varepsilon\partial_{t}\psi^{N}+\mathcal{P}_{N}(\Omega\cdot\nabla_{x}\psi^{N})+\frac{\sigma}{\varepsilon}\psi^{N}=\frac{\sigma}{\varepsilon}\overline{\psi^{N}}+\varepsilon\mathcal{P}_{N}q,\qquad\overline{\psi^{N}}=\frac{1}{4\pi}\int_{\mathbb{S}^{2}}\psi^{N}\,d\Omega,$
(7a) $\displaystyle\psi^{N}|_{t=0}=\mathcal{P}_{N}g.$ (7b)

When expressed in terms of the moments ${\bm{\uppsi}}^{N}_{\ell}$, the
$\text{P}_{N}$ method yields the following linear, symmetric hyperbolic
system:

$\displaystyle\varepsilon\partial_{t}{\bm{\uppsi}}^{N}_{{0}}$
$\displaystyle+\sum_{i=1}^{3}a_{1}^{(i)}\partial_{x_{i}}{\bm{\uppsi}}^{N}_{{1}}$
$\displaystyle=\varepsilon\mathbf{q}_{0},$ $\displaystyle\text{for }\ell=0,$
(8a) $\displaystyle\varepsilon\partial_{t}{\bm{\uppsi}}^{N}_{{\ell}}$
$\displaystyle+\sum_{i=1}^{3}(a_{\ell}^{(i)})^{\mathsf{T}}\partial_{x_{i}}{\bm{\uppsi}}^{N}_{{\ell-1}}$
$\displaystyle+\sum_{i=1}^{3}a_{\ell+1}^{(i)}\partial_{x_{i}}{\bm{\uppsi}}^{N}_{{\ell+1}}$
$\displaystyle+\frac{\sigma}{\varepsilon}{\bm{\uppsi}}^{N}_{{\ell}}$
$\displaystyle=\varepsilon\mathbf{q}_{\ell},$ $\displaystyle\text{for
}1\leq\ell\leq N-1,$ (8b)
$\displaystyle\varepsilon\partial_{t}{\bm{\uppsi}}^{N}_{{N}}$
$\displaystyle+\sum_{i=1}^{3}(a_{N}^{(i)})^{\mathsf{T}}\partial_{x_{i}}{\bm{\uppsi}}^{N}_{{N-1}}$
$\displaystyle+\frac{\sigma}{\varepsilon}{\bm{\uppsi}}^{N}_{{N}}$
$\displaystyle=\varepsilon\mathbf{q}_{N},$ $\displaystyle\text{for }\ell=N.$
(8c)

Formulas for the elements in the matrices
$a_{\ell}^{(i)}\in\mathbb{R}^{(2\ell-1)\times(2\ell+1)}$ can be found in the
appendix of [13]. In the current work, we rely only on the fact the they are
bounded in the operator norm, specifically that $\|a_{\ell}^{(i)}\|_{2}\leq
4$.

The exact moments
${\bm{\uppsi}}_{{\ell}}=\int_{\mathbb{S}^{2}}\mathbf{m}_{\ell}\psi d\Omega$
satisfy an infinite system of equations with a structure similar to (8):

$\displaystyle\varepsilon\partial_{t}{\bm{\uppsi}}_{{0}}$
$\displaystyle+\sum_{i=1}^{3}a_{1}^{(i)}\partial_{x_{i}}{\bm{\uppsi}}_{{1}}$
$\displaystyle=\varepsilon\mathbf{q}_{0},$ $\displaystyle\text{for }\ell=0,$
(9) $\displaystyle\varepsilon\partial_{t}{\bm{\uppsi}}_{{\ell}}$
$\displaystyle+\sum_{i=1}^{3}(a_{\ell}^{(i)})^{\mathsf{T}}\partial_{x_{i}}{\bm{\uppsi}}_{{\ell-1}}$
$\displaystyle+\sum_{i=1}^{3}a_{\ell+1}^{(i)}\partial_{x_{i}}{\bm{\uppsi}}_{{\ell+1}}$
$\displaystyle+\frac{\sigma}{\varepsilon}{\bm{\uppsi}}_{{\ell}}$
$\displaystyle=\varepsilon\mathbf{q}_{\ell},$ $\displaystyle\text{for
}\ell\geq 1.$ (10)

In particular, the $\text{P}_{N}$ equations for ${\bm{\uppsi}}^{N}$ can be
obtained by truncating (9) at $\ell=N$ and then neglecting the moment
${\bm{\uppsi}}^{N}_{{N+1}}$ that would otherwise appear in (8c).

###  2.3 The hybrid method

The hybrid method is based on a separation of $\psi$ into a collided component
$\psi_{\mathrm{c}}$ and an uncollided component $\psi_{\mathrm{u}}$. These
components satisfy the coupled system

$\displaystyle\varepsilon\partial_{t}\psi_{\mathrm{u}}+\Omega\cdot\nabla_{x}\psi_{\mathrm{u}}+\frac{\sigma}{\varepsilon}\psi_{\mathrm{u}}$
$\displaystyle=\varepsilon q,$ (11a)
$\displaystyle\varepsilon\partial_{t}\psi_{\mathrm{c}}+\Omega\cdot\nabla_{x}\psi_{\mathrm{c}}+\frac{\sigma}{\varepsilon}\psi_{\mathrm{c}}$
$\displaystyle=\frac{\sigma}{\varepsilon}\overline{\psi_{\mathrm{c}}}+\frac{\sigma}{\varepsilon}\overline{\psi_{\mathrm{u}}},$
(11b)

where, as before, a bar denotes the angular average of $\mathbb{S}^{2}$. The
idea of the hybrid is to solve (11) using a high-resolution angular
discretization for $\psi_{\mathrm{u}}$ and a low-resolution angular
discretization for $\psi_{\mathrm{c}}$ over a time step $\Delta t=T/M$, where
$M\in\mathbb{N}_{>0}$ and then perform a reconstruction to reinitialize
$\psi_{\mathrm{u}}$ and $\psi_{\mathrm{c}}$ for the next time step. To
formalize this procedure, define a set of temporal grid points
$0=t_{0}<t_{1}<...<t_{M}=T$, and for $m\in\\{1,2,\dots,M\\}$, let
$f(t_{m}^{-})=\lim_{\delta\to 0^{+}}f(t_{m}-\delta)$ for any function $f$ of
$t$ that is continuous on $[t_{m-1},t_{m})$. Then for $m\in\\{1,2,\dots,M\\}$,
$\psi_{{\mathrm{u}},m}$ and $\psi_{{\mathrm{c}},m}$ satisfy the following
system of equations over the interval $[t_{m-1},t_{m})$

$\displaystyle\varepsilon\partial_{t}\psi_{{\mathrm{u}},m}+\Omega\cdot\nabla_{x}\psi_{{\mathrm{u}},m}+\frac{\sigma}{\varepsilon}\psi_{{\mathrm{u}},m}=\varepsilon
q,$ (12a)
$\displaystyle\varepsilon\partial_{t}\psi_{\mathrm{c},m}+\Omega\cdot\nabla_{x}\psi_{\mathrm{c},m}+\frac{\sigma}{\varepsilon}\psi_{\mathrm{c},m}=\frac{\sigma}{\varepsilon}(\overline{\psi_{\mathrm{u},m}}+\overline{\psi_{\mathrm{c},m}}),$
(12b)
$\displaystyle\left.\psi_{\mathrm{u},m}\right|_{t=t_{m-1}}=\begin{cases}g,&\quad
m=1,\\\
\psi_{\mathrm{u},m-1}(t_{m-1}^{-})+\psi_{\mathrm{c},m-1}(t_{m-1}^{-})&\quad
m>1.\end{cases}$ (12c)
$\displaystyle\left.\psi_{\mathrm{c},m}\right|_{t=t_{m-1}}=0.$ (12d)

The intuition behind this splitting is that (12a) can be discretized with a
high-resolution angular discretization but solved more efficiently than (3)
since angular unknowns are no longer coupled. Although (12b) features the same
type of angular coupling as (3), it can be solved with fewer degrees of
freedom because the source $\overline{\psi_{\mathrm{u},m}}$ is, in general,
more regular than $q$. However, because $\psi_{{\mathrm{u}},m}$ decays
exponentially whenever $\sigma>0$, the hybrid is only solved for a time step
$\Delta t$ before the relabeling in (12d)- (12c) is implemented. The enables
the hybrid to capture more high-resolution features than (12b) can do alone.

###  2.4 Angular discretization of the hybrid

We focus now on the angular discretization of (12). The strategy of the hybrid
is to discretize (12a) in angle with a high-resolution method and (12b) in
angle with a low-resolution method. In practice, there are a variety of
strategies and combinations available to do so. For the purposes of analysis,
we assume that (12a) is solved exactly and that (12b) is discretized with a
$\text{P}_{N}$ method. That is, we seek
$\psi^{N}=\psi_{\mathrm{u},m}^{N}+\psi_{\mathrm{c},m}^{N}$ where for each
$m\in\\{1,2,\dots M\\}$

$\left(\psi_{\mathrm{u},m}^{N},\psi_{\mathrm{c},m}^{N}\right)\in
C([t_{m-1},t_{m});X\times L^{2}(\mathbb{S}^{2}))\times
C([t_{m-1},t_{m});X\times\mathbb{P}_{N}(\mathbb{S}^{2}))$ (13)

satisfies

$\displaystyle\varepsilon\partial_{t}\psi_{{\mathrm{u}},m}^{N}+\Omega\cdot\nabla_{x}\psi_{{\mathrm{u}},m}^{N}+\frac{\sigma}{\varepsilon}\psi_{{\mathrm{u}},m}^{N}=\varepsilon
q,$ (14a)
$\displaystyle\varepsilon\partial_{t}\psi_{\mathrm{c},m}^{N}+\mathcal{P}_{N}\left(\Omega\cdot\nabla_{x}\psi_{\mathrm{c},m}^{N}\right)+\frac{\sigma}{\varepsilon}\psi_{\mathrm{c},m}^{N}=\frac{\sigma}{\varepsilon}(\overline{\psi_{\mathrm{u},m}^{N}}+\overline{\psi_{\mathrm{c},m}^{N}}),$
(14b)
$\displaystyle\left.\psi_{\mathrm{c},m}^{N}\right|_{t=t_{m-1}}=0,\quad\left.\psi_{\mathrm{u},m}^{N}\right|_{t=t_{m-1}}=\begin{cases}g,&\quad
m=1,\\\
\psi_{\mathrm{u},m-1}^{N}(t^{-}_{m-1})+\psi_{\mathrm{c},m-1}^{N}(t^{-}_{m-1})&\quad
m>1.\end{cases}$ (14d)

We compare the accuracy of the solution defined in (14) with the monolithic
$\text{P}_{N}$ method (7a), using the same value of $N$. To simplify the
numerical analysis of these two models, we keep the time and space variables
continuous.

Since (7a) and (14b) have the same computational complexity, the goal is to
asses the additional benefit of solving (14a). Clearly the additional cost of
solving (14a) involves both memory and run-time; both are fairly easy to
quantify. However, assessing the gains in accuracy is not as simple. Thus is
important to understand these gains in order to better quantify observed
improvements in run-time efficiency provided by the hybrid.

###  2.5 Preview of main results.

Let $s\geq 1$ be the number of angular $L^{2}$ derivatives in the solution
$\psi$ and let $N$ be an integer such that $N\geq s-1$. Let
$e^{N}=\psi-\psi^{N}$ be the error in the $\text{P}_{N}$ approximation (7a),
and let $e^{N}_{M}=\psi-(\psi_{\mathrm{u},M}^{N}+\psi_{\mathrm{c},M}^{N})$ be
the error in the hybrid at the $M$-th time step. The main results of the paper
are the $\text{P}_{N}$ error estimate in Theorem 13 and the hybrid error
estimate in Theorem 18. We compare these errors for two different regimes:
first, when $\sigma\asymp 1$222Recall that $a\asymp b$ if and only if $a=O(b)$
and $b=O(a)$ and $\varepsilon\to 0$ (the diffusion regime) and second, when
$\varepsilon\asymp(1)$ and $\sigma\to 0$ (the purely absorbing
regime).333Technically, this is the streaming regime for $\Psi$, since there
is no absorption.

When $\varepsilon\ll 1$ and $\sigma\asymp 1$, Theorem 13 implies that

$\|{e^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(T)\lesssim\frac{T}{(N+1)^{s}}\left[e^{-\sigma
T/\varepsilon^{2}}\sum_{i=0}^{s-1}\frac{T^{i}}{\varepsilon^{i+1}}+\frac{\varepsilon^{s-1}}{\sigma^{s}}+O\left(\frac{\varepsilon}{\sigma}\right)\right].$
(15)

Thus with sufficient regularity, the $\text{P}_{N}$ approximation is
spectrally accurate and grows linearly in time for large $T$. The first and
third term in brackets depend on whether or not $q$ and $g$ are isotropic
(independent of angle). The first term is due to initial layers when $g$ is
non-isotropic, and the third terms arises when $q$ is non-isotropic. When $g$
and $q$ are isotropic, the $\text{P}_{N}$ error reduces to (see Corollary 2)

$\|{e^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(T)\lesssim\frac{\varepsilon^{s-1}T}{\sigma^{s}(N+1)^{s}}.$
(16)

For the hybrid method the initial condition and source are always isotropic.
Thus Theorem 18 gives the following compact bound

$\|{e_{M}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(T)\lesssim\frac{\varepsilon^{s-1}T}{\sigma^{s}(N+1)^{s}}.$
(17)

Thus the hybrid method comes equipped with a better error estimate than the
$\text{P}_{N}$ method, and the estimates agree when the data is isotropic.
This suggests that the hybrid method is at least as accurate as the
$\text{P}_{N}$ approximation when $\varepsilon\ll 1$ and $\sigma\asymp 1$. In
addition, the hybrid estimate is independent of the time step $\Delta t$ used
for re-initialization in this regime.

For both the hybrid method and $\text{P}_{N}$ approximation, the errors
converge to zero as $\varepsilon\to 0$ whenever $s>1$ (modulo the initial
layer in (15)). This fact is consistent with the fact that the $\text{P}_{N}$
method recovers the diffusion limit (2a) whenever $N\geq 1$ [14].

For the second regime of interest, $\sigma\ll 1$ and $\varepsilon\asymp 1$,
Theorem 13 gives the $\text{P}_{N}$ error estimate

$\|{e^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(T)\lesssim\frac{\mathfrak{p}_{s}(T/\varepsilon)}{(N+1)^{s}},$
(18)

where $\mathfrak{p}_{s}(\omega)=\sum_{i=0}^{s+1}c_{i}(T)\omega^{i}$ is a
polynomial of degree $s+1$ with non-negative coefficients
$c_{i}(T)=a_{i}T+b_{i}$, $a_{i},\,b_{i}\geq 0$. Here the hybrid estimate
provides a significant improvement:

$\|{e_{M}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(T)\lesssim\frac{\Delta
t^{s}T}{\varepsilon^{s+1}(N+1)^{s}}\min(1,\frac{\Delta
t\sigma}{\varepsilon^{2}})$ (19)

In particular, $\|{e_{M}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(T)=0$, when
$\sigma=0$. This result is expected since in that case the uncollided solution
and the transport solution agree. As expected, the error is monotonic in
$\Delta t$; however, small time steps require more evaluations of the
uncollided equation and more reintializions. In practice, this additional cost
must be taken into account.

##  3  $\text{P}_{N}$ Analysis

###  3.1 Spherical harmonics preliminaries

A natural space to analyze the transport equation and the $\text{P}_{N}$
approximation is the Sobolev space $H_{\circ}^{s}(\mathbb{S}^{2})$. To
describe this space, we recall some elementary facts about spherical harmonics
which can be found, for example, in [3, 11]. For $u,v\in
L^{2}(\mathbb{S}^{2})$, let

$(u,v)_{L^{2}(\mathbb{S}^{2})}=\sum_{\ell=0}^{\infty}\mathbf{u}_{\ell}^{\top}\mathbf{v}_{\ell}\qquad\text{and}\qquad\|u\|_{L^{2}(\mathbb{S}^{2})}^{2}=\sum_{\ell=0}^{\infty}\|\mathbf{u}_{\ell}\|^{2},$
(20)

where

$\mathbf{u}_{\ell}^{\top}\mathbf{v}_{\ell}=\sum_{k=-\ell}^{\ell}u_{\ell,k}v_{\ell,k}\qquad\text{and}\qquad\|\mathbf{u}_{\ell}\|^{2}=\mathbf{u}_{\ell}^{\top}\mathbf{u}_{\ell}=\sum_{k=-\ell}^{\ell}|u_{\ell,k}|^{2}.$
(21)

A standard way to define Sobolev spaces on the sphere is via the Laplace-
Beltrami operator $\Delta_{\circ}$, which is the spherical component of the
Laplacian and for which the spherical harmonics are eigenfunctions:
$-\Delta_{\circ}Y_{\ell,j}=\ell(\ell+1)Y_{\ell,j}.$ For even integers $s$, the
usual norm is

$\|u\|_{H^{s}_{\circ}(\mathbb{S}^{2})}:=\left\|\left(\frac{1}{4}+\Delta_{\circ}\right)^{s/2}u\right\|=\left(\sum_{\ell=0}^{\infty}b_{\ell}\|\mathbf{u}_{\ell}\|^{2}\right)^{1/2},\qquad
b_{\ell}=\left(\frac{1}{2}+\ell\right)^{2s}.$ (22)

The definition of this inner product extends naturally to all
$s\in\mathbb{R}$, and the space $H^{s}_{\circ}(\mathbb{S}^{2})$ is then the
completion of smooth functions under the $H^{s}_{\circ}(\mathbb{S}^{2})$ norm
[3].

Rather than working directly with the $H^{s}_{\circ}(\mathbb{S}^{2})$ norm in
(22), it is convenient in the analysis below to use an equivalent norm. For
$s\geq 0$, define the $H^{s}(\mathbb{S}^{2})$ semi-norm and norm by

$\displaystyle|u|_{H^{s}(\mathbb{S}^{2})}:=\left(\sum_{\ell=s}^{\infty}b_{\ell}\|\mathbf{u}_{\ell}\|^{2}\right)^{1/2}\quad\text{and}\quad\|u\|_{H^{s}(\mathbb{S}^{2})}=\left(s\|u\|_{L^{2}(\mathbb{S}^{2})}^{2}+|u|^{2}_{H^{s}(\mathbb{S}^{2})}\right)^{1/2},$
(23)

respectively, where the sum in the semi-norm definition in (23) begins at $s$
for technical arguments that are used in the proof of Lemma (5) below. When
$s=0$, the norms coincide:
$\|u\|_{H^{0}(\mathbb{S}^{2})}=\|u\|_{H^{0}_{\circ}(\mathbb{S}^{2})}=\|u\|_{L^{2}(\mathbb{S}^{2})}$.
More generally, the following equivalence holds.

######  Lemma 1 (Norm equivalence).

For any $s\geq 0$,

$c_{1}(s)\|u\|_{H^{s}(\mathbb{S}^{2})}\leq\|u\|_{H^{s}_{\circ}(\mathbb{S}^{2})}\leq
c_{2}(s)\|u\|_{H^{s}(\mathbb{S}^{2})}.$ (24)

where

$c_{1}(s)=\begin{cases}1&\text{if}\quad s=0,\\\
\frac{1}{\sqrt{3s}}&\text{if}\quad s\geq 1\end{cases}\quad\text{and}\quad
c_{2}(s)=\begin{cases}1&\text{if}\quad s=0,\\\
\sqrt{\frac{5}{s}}\left(s-\frac{1}{2}\right)^{s}&\text{if}\quad s\geq
1\end{cases}.$ (25)

###### Proof.

When $s=0$, the norms are equal, so (24) holds trivially. Thus assume that
$s\geq 1$. The first inequality in (24) follows from the fact that

$\begin{split}\|u\|_{H^{s}(\mathbb{S}^{2})}^{2}&=s\sum_{\ell=0}^{\infty}\|\mathbf{u}_{\ell}\|^{2}+\sum_{\ell=s}^{\infty}\left(\frac{1}{2}+\ell\right)^{2s}\|\mathbf{u}_{\ell}\|^{2}\leq
3s\sum_{\ell=0}^{\infty}\left(\frac{1}{2}+\ell\right)^{2s}\|\mathbf{u}_{\ell}\|^{2}=\frac{1}{[c_{1}(s)]^{2}}\|u\|_{H^{s}_{\circ}(\mathbb{S}^{2})}^{2}.\end{split}$
(26)

To prove the second inequality in (24), we use the elementary inequality

$\frac{s}{4}\leq\left(s-\frac{1}{2}\right)^{2s},\quad s\geq 1,$ (27)

to conclude that

$\begin{split}\|u\|_{H^{s}_{\circ}(\mathbb{S}^{2})}^{2}&=\sum_{\ell=0}^{s-1}\left(\frac{1}{2}+\ell\right)^{2s}\|\mathbf{u}_{\ell}\|^{2}+|u|_{H^{s}(\mathbb{S}^{2})}^{2}\leq\frac{1}{s}\left(s-\frac{1}{2}\right)^{2s}s\|{u}\|_{L^{2}(\mathbb{S}^{2})}^{2}+|u|_{H^{s}(\mathbb{S}^{2})}^{2}\\\
&\leq\frac{1}{s}\left(s-\frac{1}{2}\right)^{2s}s\|{u}\|_{L^{2}(\mathbb{S}^{2})}^{2}+\frac{4}{s}\left(s-\frac{1}{2}\right)^{2s}|u|_{H^{s}(\mathbb{S}^{2})}^{2}\\\
&\leq\frac{5}{s}\left(s-\frac{1}{2}\right)^{2s}\left(s\|{u}\|_{L^{2}(\mathbb{S}^{2})}^{2}+|u|_{H^{s}(\mathbb{S}^{2})}^{2}\right)=[c_{2}(s)]^{2}\|u\|_{H^{s}(\mathbb{S}^{2})}^{2}.\qquad\qquad\qquad\end{split}$
(28)

∎

######  Lemma 2 (Approximation property).

For $s\geq 0$ and $N\geq\max\\{0,s-1\\}$,

$\|(\mathcal{I}-\mathcal{P}_{N})u\|_{L^{2}(\mathbb{S}^{2})}\leq\frac{1}{(N+1)^{s}}|(\mathcal{I}-\mathcal{P}_{N})u|_{H^{s}(\mathbb{S}^{2})}\leq\frac{1}{(N+1)^{s}}|u|_{H^{s}(\mathbb{S}^{2})}.$
(29)

###### Proof.

From the definition of the projection $\mathcal{P}_{N}$ in (5),

$\displaystyle\underbrace{\sum_{\ell=N+1}^{\infty}\|\mathbf{u}_{\ell}\|^{2}}_{=\|(\mathcal{I}-\mathcal{P}_{N})u\|^{2}_{L^{2}(\mathbb{S}^{2})}}\leq\frac{1}{(N+1)^{2s}}\underbrace{\sum_{\ell=N+1}^{\infty}\left(\frac{1}{2}+\ell\right)^{2s}\|\mathbf{u}_{\ell}\|^{2}}_{=|(\mathcal{I}-\mathcal{P}_{N})u|^{2}_{H^{s}(\mathbb{S}^{2})}}\leq\frac{1}{(N+1)^{2s}}\underbrace{\sum_{\ell=s}^{\infty}\left(\frac{1}{2}+\ell\right)^{2s}\|\mathbf{u}_{\ell}\|^{2}}_{=|u|^{2}_{H^{s}(\mathbb{S}^{2})}}$
(30)

Taking square roots of each term above yields the desired result. ∎

For vector-valued functions of space, we define the usual $L^{2}(X)$ inner-
product and norm by

$(\mathbf{v},\mathbf{w})_{L^{2}(X)}=\int_{X}\mathbf{v}(x)^{\top}\mathbf{w}(x)\,dx\quad\text{and}\quad\|{\mathbf{v}}\|^{2}_{L^{2}(X)}=(\mathbf{v},\mathbf{v})_{L^{2}(X)}=\int_{X}\|\mathbf{v}(x)\|^{2}\,dx.$
(31)

The space $H^{0,s}(X\times\mathbb{S}^{2})=L^{2}(X;H^{s}(\mathbb{S}^{2}))$ is
the space of measurable functions $u\colon X\times\mathbb{S}^{2}\to\mathbb{R}$
with the semi-inner product

$\displaystyle(u,v)_{H^{0,s}(X\times\mathbb{S}^{2})}=\int_{X}\sum_{\ell=s}^{\infty}\left(\frac{1}{2}+\ell\right)^{2s}\mathbf{u}_{\ell}(x)^{\top}\mathbf{v}_{\ell}(x)\,dx=\sum_{\ell=s}^{\infty}\left(\frac{1}{2}+\ell\right)^{2s}(\mathbf{u}_{\ell},\mathbf{v}_{\ell})_{L^{2}(X)}$
(32)

such that the semi-norm

$\displaystyle|u|^{2}_{H^{0,s}(X\times\mathbb{S}^{2})}=\int_{X}\sum_{\ell=s}^{\infty}\left(\frac{1}{2}+\ell\right)^{2s}|\mathbf{u}_{\ell}(x)|^{2}dx=\sum_{\ell=s}^{\infty}\left(\frac{1}{2}+\ell\right)^{2s}\|\mathbf{u}_{\ell}\|_{L^{2}(X)}^{2}$
(33)

is bounded.

For $r\in\mathbb{N}_{\geq 0}$ and $s\geq 0$ we define
$H^{r,s}(X\times\mathbb{S}^{2})$ to be the space of functions $u\colon
X\times\mathbb{S}^{2}\to\mathbb{R}$ such that the semi-norms.

$\displaystyle|u|_{H^{\varsigma,s}(X\times\mathbb{S}^{2})}:=\sum_{i_{1},i_{2},\dots,i_{\varsigma}=1}^{3}|\partial_{x_{i_{1}}x_{i_{2}}\dots
x_{i_{\varsigma}}}u|_{H^{0,s}(X\times\mathbb{S}^{2})}.$ (34)

are bounded for all positive integers $\varsigma\leq r$. These semi-norms
above are equivalent to the standard semi-norms, but are more convenient in
the context of the RTE since

$|\Omega\cdot\nabla_{x}u|_{H^{r,s}(X\times\mathbb{S}^{2})}\leq\sum_{i=1}^{3}|\partial_{x_{i}}u|_{H^{r,s}(X\times\mathbb{S}^{2})}=|u|_{H^{r+1,s}(X\times\mathbb{S}^{2})},$
(35)

which will be used in the analysis below. Henceforth, the domain of
integration for $H^{s}$ and $H^{r,s}$ will be left off when there is no
ambiguity, i.e.,

$|u|_{H^{s}}:=|u|_{H^{s}(\mathbb{S}^{2})}\quad\text{and}\quad|u|_{H^{r,s}}:=|u|_{H^{r,s}(X\times\mathbb{S}^{2})}$
(36)

Finally, for $p\geq 1$ and measurable functions $u\colon
X\times\mathbb{S}^{2}\times[\alpha,\beta]\to\mathbb{R}$, we denote the space-
time semi-norms by

$|u|_{L^{p}([\alpha,\beta];H^{r,s})}=\left(\int_{\alpha}^{\beta}|u|_{H^{r,s}}^{p}\,d\tau\right)^{1/p}\quad\text{and}\quad|u|_{L^{\infty}([\alpha,\beta]);H^{r,s})}=\operatorname*{ess~{}sup}_{t\in[\alpha,\beta]}{|u|_{H^{r,s}}}.$
(37)

###  3.2 Stability of the $\text{P}_{N}$ system

In this section, we derive estimates on high-order semi-norms that arise in
the subsequent error analysis. The analysis requires iterated inequalities of
Grönwall type (see Lemma 21 in the Appendix), so we define the following
notation to simplify integrals that arise. Let $0\leq\alpha\leq
t\leq\beta<\infty$, and for any function $f\in L^{1}([\alpha,\beta])$, define
the bounded linear operator $\mathcal{A}_{\alpha}\colon
L^{1}([\alpha,\beta])\rightarrow C^{0}([\alpha,\beta])$ and its powers
$\mathcal{A}^{k}_{\alpha},k\in\mathbb{N}^{\geq 0}$, by

$\mathcal{A}_{\alpha}[f](t):=\frac{1}{\varepsilon}\int_{\alpha}^{t}e^{-\sigma(t-\tau)/\varepsilon^{2}}f(\tau)d\tau\quad\text{and}\quad\mathcal{A}^{k}_{\alpha}[f](t)=\begin{cases}\mathcal{I},&k=0,\\\
\mathcal{A}_{\alpha}[\mathcal{A}^{k-1}_{\alpha}[f]](t),&k\geq 1.\end{cases}$
(38)

In addition, let $\mathbbm{1}\in L^{1}([\alpha,\beta])$ be the function that
is identically one, and let

$F_{\alpha}(t)=e^{-\sigma(t-\alpha)/\varepsilon^{2}},\quad
t\in[\alpha,\beta].$ (39)

It is clear from the definition in (38) $\mathcal{A}_{\alpha}$ is monotonic;
that is, if $0\leq f(t)\leq g(t)$ for a.e. $t\in[\alpha,\beta]$, then
$\mathcal{A}_{\alpha}[f](t)\leq\mathcal{A}_{\alpha}[g](t)$ for all
$t\in[\alpha,\beta]$.

######  Lemma 3.

Let $F_{\alpha}$ be given as in (39). Then for all $t\geq\alpha$ and every
$k\in\mathbb{N}_{\geq 0}$,

$\mathcal{A}^{k}_{\alpha}[\mathbbm{1}](t)\leq\min\left(\frac{\varepsilon^{k}}{\sigma^{k}},\frac{1}{k!}\left(\frac{t-\alpha}{\varepsilon}\right)^{k}\right)\quad\text{and}\quad\mathcal{A}^{k}_{\alpha}[F_{\alpha}](t)=\frac{(t-\alpha)^{k}e^{-\sigma(t-\alpha)/\varepsilon^{2}}}{k!\varepsilon^{k}}$
(40)

###### Proof.

We first prove the bound in (40). Since
$e^{-\sigma(t-\tau)/\varepsilon^{2}}\leq 1$, a direct calculation gives

$\mathcal{A}^{k}_{\alpha}[\mathbbm{1}](t)\leq\frac{1}{\varepsilon^{k}}\int_{\alpha}^{t}\int_{\alpha}^{\tau_{k-1}}\cdots\int_{\alpha}^{\tau_{1}}\mathbbm{1}d\tau_{0}\cdots
d\tau_{k-1}=\frac{1}{k!}\left(\frac{t-\alpha}{\varepsilon}\right)^{k}.$ (41)

On the other hand, it follows directly from the definition of
$\mathcal{A}_{\alpha}$ that

$\mathcal{A}_{\alpha}[\mathbbm{1}](t)=\frac{\varepsilon}{\sigma}\left(1-F_{\alpha}(t)\right)\leq\frac{\varepsilon}{\sigma}\quad\implies\quad\mathcal{A}^{k}_{\alpha}[\mathbbm{1}](t)\leq\frac{\varepsilon^{k}}{\sigma^{k}}$
(42)

Together (41) and (42) yield (40).

We prove the second statement in (40) by induction on $k$. When $k=0$, the
statement follows trivially,

$\mathcal{A}_{\alpha}^{0}[F_{\alpha}](t)=F_{\alpha}(t)=\frac{(t-\alpha)^{0}e^{-\sigma(t-\alpha)/\varepsilon^{2}}}{0!e^{0}},$
(43)

Now let us assume the statement is true for arbitrary $k$, and

$\mathcal{A}_{\alpha}^{k+1}[F_{\alpha}](t)=\mathcal{A}_{\alpha}\left[\mathcal{A}_{\alpha}^{k}[F_{\alpha}]\right](t)=\frac{1}{\varepsilon}\int_{\alpha}^{t}e^{-\sigma(t-\tau)/\varepsilon^{2}}\frac{(\tau-\alpha)^{k}e^{-\sigma(\tau-\alpha)/\varepsilon^{2}}}{k!\varepsilon^{k}}\,d\tau=\frac{(t-\alpha)^{k+1}e^{-\sigma(t-\alpha)/\varepsilon^{2}}}{(k+1)!\varepsilon^{k+1}}$
(44)

∎

####  3.2.1 Estimates for the $\text{P}_{N}$ system

We derive evolution equations for the $H^{r,s}$ semi-norms, defined in (34).
For $r=0$, we test each equation in (8) by
$b_{\ell}{\bm{\uppsi}}^{N}_{{\ell}}$, integrate by parts, and then sum over
$\ell\in\\{s,s+1,...,N\\}$. This gives

$\begin{split}\frac{\varepsilon}{2}\partial_{t}|\psi^{N}|_{H^{0,s}}^{2}&+\frac{\sigma}{\varepsilon}|\psi^{N}|_{H^{0,s}}^{2}=\varepsilon(\mathcal{P}_{N}q,\psi^{N})_{H^{0,s}}\\\
&\quad-\sum_{i=1}^{3}\sum_{\ell=s}^{N}\left(\left(\ell+\frac{1}{2}\right)^{2s}-\gamma_{s,\ell}\left(\ell-\frac{1}{2}\right)^{2s}\right)\left({\bm{\uppsi}}^{N}_{{\ell}},\left(a_{\ell}^{(i)}\right)^{T}\partial_{x_{i}}{\bm{\uppsi}}^{N}_{{\ell-1}}\right)_{L^{2}(X)},\end{split}$
(45)

where $\gamma_{s,\ell}=(1-\delta_{s,\ell})$ is used to handle the first non-
zero term in the sum over $\ell$. For the special case $s=0$, (45) recovers
the usual $L^{2}$ energy equation:

$\frac{\varepsilon}{2}\frac{d}{dt}\|{\psi^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}^{2}+\frac{\sigma}{\varepsilon}\|{\psi^{N}-\overline{\psi^{N}}}\|_{L^{2}(X\times\mathbb{S}^{2})}^{2}=\varepsilon(\mathcal{P}_{N}q,\psi^{N})_{L^{2}(X\times\mathbb{S}^{2})}^{2}.$
(46)

To find a closed estimate with respect to the $H^{0,s}$ semi-norms, we focus
on the summation in (45).

######  Lemma 4.

Let $s\geq 1$ and $\ell\geq s$. Then

$\left(\ell+\frac{1}{2}\right)^{2s}-\gamma_{s,\ell}\left(\ell-\frac{1}{2}\right)^{2s}\leq
2es\left(\ell+\frac{1}{2}\right)^{s}\left(\ell-\frac{1}{2}\right)^{s-1}.$ (47)

###### Proof.

We first establish an elementary inequality. Since $\ell\geq s$,

$\ell+\frac{1}{2}=\left(\ell-\frac{1}{2}\right)\left(1+\frac{1}{\ell-\frac{1}{2}}\right)\leq\left(\ell-\frac{1}{2}\right)\left(1+\frac{1}{s-\frac{1}{2}}\right)\leq\left(\ell-\frac{1}{2}\right)e^{1/(s-1/2)}$
(48)

Therefore

$\left(\ell+\frac{1}{2}\right)^{s-1}\leq
e^{\frac{s-1}{s-1/2}}\left(\ell-\frac{1}{2}\right)^{s-1}\leq
e\left(\ell-\frac{1}{2}\right)^{s-1}.$ (49)

We use (49) to show (47). There are two cases:

Case 1 ($\ell=s$):

In this case, $\gamma_{s,\ell}=0$. Since $2s>s+1/2$ and by (49),

$\left(\ell+\frac{1}{2}\right)^{2s}=\left(s+\frac{1}{2}\right)\left(\ell+\frac{1}{2}\right)^{s}\left(\ell+\frac{1}{2}\right)^{s-1}\leq
2es\left(\ell+\frac{1}{2}\right)^{s}\left(\ell-\frac{1}{2}\right)^{s-1}$ (50)

Case 2 ($\ell>s$):

In this case, $\gamma_{s,\ell}=1$. Applying a binomial expansion and then (49)
gives

$\begin{split}&\left(\ell+\frac{1}{2}\right)^{2s}-\left(\ell-\frac{1}{2}\right)^{2s}=\sum_{k=0}^{2s}\binom{2s}{k}\left(\ell-\frac{1}{2}\right)^{k}-\left(\ell-\frac{1}{2}\right)^{2s}=\sum_{k=0}^{2s-1}\binom{2s}{k}\left(\ell-\frac{1}{2}\right)^{k}\\\
&\quad=\sum_{k=0}^{2s-1}\frac{2s}{2s-k}\binom{2s-1}{k}\left(\ell-\frac{1}{2}\right)^{k}\leq
2s\sum_{k=0}^{2s-1}\binom{2s-1}{k}\left(\ell-\frac{1}{2}\right)^{k}\\\
&\quad=2s\left(\ell+\frac{1}{2}\right)^{2s-1}=2s\left(\ell+\frac{1}{2}\right)^{s}\left(\ell+\frac{1}{2}\right)^{s-1}\leq
2es\left(\ell+\frac{1}{2}\right)^{s}\left(\ell-\frac{1}{2}\right)^{s-1}\end{split}$
(51)

∎

######  Lemma 5 (Semi-norm recurrence).

Let $s\geq 1$, $q\in L^{1}([0,T];H^{r,s})$ and $g\in H^{r,s}$. Then for all
$t\in[0,T]$,

$\displaystyle|\psi^{N}|_{H^{r,s}}(t)\leq
Cs\mathcal{A}_{0}[|\psi^{N}|_{H^{r+1,s-1}}](t)+|\mathcal{P}_{N}g|_{H^{r,s}}F_{0}(t)+\varepsilon\mathcal{A}_{0}[|\mathcal{P}_{N}q|_{H^{r,s}}](t),$
(52)

where $F_{0}$ is defined in (39) and $C$ is a constant independent of the
data.

###### Proof.

We assume first that $r=0$ and focus on the last term in (45). It follows from
(i) the induced norm bound $\|a_{\ell}^{(i)}\|_{2}\leq 4$ [13], (ii) the
bounds in Lemma 4, (iii) the Cauchy-Schwarz inequality, and (iv) the $H^{r,s}$
semi-norm definitions in (33) and (34) that

$\begin{split}&\sum_{\ell=s}^{N}\left(\left(\ell+\frac{1}{2}\right)^{2s}-\gamma_{s,\ell}\left(\ell-\frac{1}{2}\right)^{2s}\right)\left({\bm{\uppsi}}^{N}_{{\ell}},\left(a_{\ell}^{(i)}\right)^{T}\partial_{x_{i}}{\bm{\uppsi}}^{N}_{{\ell-1}}\right)_{L^{2}(X)}\\\
&\quad\leq
2es\sum_{\ell=s}^{N}\left(\ell+\frac{1}{2}\right)^{s}\left(\ell-\frac{1}{2}\right)^{s-1}\|a_{\ell}^{(i)}\|_{2}\enspace\|{\bm{\uppsi}}^{N}_{{\ell}}\|_{L^{2}(X)}\enspace\|\partial_{x_{i}}{\bm{\uppsi}}^{N}_{{\ell-1}}\|_{L^{2}(X)}\\\
&\quad\leq
Cs\left(\sum_{\ell=s}^{N}\left(\ell+\frac{1}{2}\right)^{2s}|{\bm{\uppsi}}^{N}_{{\ell}}|^{2}_{L^{2}(X)}\right)^{1/2}\left(\sum_{\ell=s}^{N}\left(\ell-\frac{1}{2}\right)^{2(s-1)}|\partial_{x_{i}}{\bm{\uppsi}}^{N}_{{\ell-1}}|^{2}_{L^{2}(X)}\right)^{1/2}\\\
&\quad\leq
Cs|\psi^{N}|_{H^{0,s}}|\partial_{x_{i}}\psi^{N}|_{H^{0,s-1}}.\end{split}$ (53)

Applying the bound above to the right-hand side of (45) and applying Lemma 21
gives (52). The case $r>1$ can be handled by differentiating the
$\text{P}_{N}$ equations in space and then repeating the arguments above. ∎

For a general function $\phi\in L^{p}([\alpha,\beta];H^{r,s})$, let
$|\phi|_{L^{p}([\alpha,\bullet];H^{0,r+s})}:[\alpha,\beta]\rightarrow\mathbb{R}$
be the map defined

$|\phi|_{L^{p}([\alpha,\bullet];H^{r,s})}(\tau)=|\phi|_{L^{p}([\alpha,\tau];H^{r,s})}.$

######  Lemma 6 (Stability of higher-order semi-norms).

Let $q\in L^{1}([0,T];H^{r,0})$ and $g\in H^{r,0}$. If $t\in[0,T]$, then

$\displaystyle|\psi^{N}|_{H^{r,0}}(t)\leq|\mathcal{P}_{N}g|_{H^{r,0}}+|\mathcal{P}_{N}q|_{L^{1}([0,t];H^{r,0})}.$
(54)

If, in addition, $s\geq 1$, $q\in L^{1}([0,T];H^{i,j})$ and $g\in H^{i,j}$ for
each $i,j$ such that $0\leq i\leq r$, $0\leq j\leq s$, and $i+j=r+s$,

$\begin{split}|\psi^{N}|_{H^{r,s}}(t)&\leq
C_{s}s!\mathcal{A}_{0}^{s}\left[|\mathcal{P}_{N}g|_{H^{r+s,0}}+|\mathcal{P}_{N}q|_{L^{1}([0,\bullet];H^{r+s,0})}\right](t)\\\
&+C_{s}\sum_{i=0}^{s-1}\frac{s!}{(s-i)!}\mathcal{A}_{0}^{i}[F_{0}|\mathcal{P}_{N}g|_{H^{r+i,s-i}}](t)+C_{s}\varepsilon\sum_{i=0}^{s-1}\frac{s!}{(s-i)!}\mathcal{A}_{0}^{i+1}[|\mathcal{P}_{N}q|_{H^{r+i,s-i}}](t),\end{split}$
(55)

where $F_{0}$ is given in (39) where $C_{s}$ is a constant depending only on
$s$.

###### Proof.

First we will prove (54) for $s=0$, in which case
$H^{r,0}={L^{2}(X\times\mathbb{S}^{2})}$. From (46) and the Cauchy-Schwarz
inequality,

$\frac{\varepsilon}{2}\frac{d}{dt}\|{\psi^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}^{2}\leq\varepsilon(\mathcal{P}_{N}q,\psi^{N})_{L^{2}(X\times\mathbb{S}^{2})}\leq\varepsilon\|{\mathcal{P}_{N}q}\|_{L^{2}(X\times\mathbb{S}^{2})}\|{\psi^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}$
(56)

an application of Lemma 21, gives

$\displaystyle\|{\psi^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t)\leq\|{\mathcal{P}_{N}g}\|_{L^{2}(X\times\mathbb{S}^{2})}+\|{\mathcal{P}_{N}q}\|_{L^{1}([0,t];L^{2}(X\times\mathbb{S}^{2}))}.$
(57)

For $r>0$, $\phi^{N}=\partial_{x_{i_{1}}x_{i_{2}}\dots x_{i_{r}}}\psi^{N}$
satisfies (7a), with initial condition
$\left.\phi^{N}\right|_{t=0}=\mathcal{P}_{N}(\partial_{x_{i_{1}}x_{i_{2}}\dots
x_{i_{r}}}g)=\partial_{x_{i_{1}}x_{i_{2}}\dots x_{i_{r}}}\mathcal{P}_{N}g$ and
source $\mathcal{P}_{N}(\partial_{x_{i_{1}}x_{i_{2}}\dots
x_{i_{r}}}q)=\partial_{x_{i_{1}}x_{i_{2}}\dots x_{i_{r}}}\mathcal{P}_{N}q$.
Repeating the argument above gives, in analogy with (57),

$\displaystyle\|{\partial_{x_{i_{1}}x_{i_{2}}\dots
x_{i_{r}}}\psi^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t)\leq\|{\partial_{x_{i_{1}}x_{i_{2}}\dots
x_{i_{r}}}\mathcal{P}_{N}g}\|_{L^{2}(X\times\mathbb{S}^{2})}+\|{\partial_{x_{i_{1}}x_{i_{2}}\dots
x_{i_{r}}}\mathcal{P}_{N}q}\|_{L^{1}([0,t];L^{2}(X\times\mathbb{S}^{2}))},$
(58)

Summing (58) over each permutation of the $r$ derivatives
$\partial_{x_{i_{1}}}\partial_{x_{i_{2}}}\cdots\partial_{x_{i_{r}}}$ and using
the definitions in (34) recovers (54).

We next prove (55). Let

$\displaystyle
b_{r,s}(t)=|\mathcal{P}_{N}g|_{H^{r,s}}F_{0}(t)+\varepsilon\mathcal{A}_{0}[|\mathcal{P}_{N}q|_{H^{r,s}}](t)\quad\text{and}\quad
c_{r,s}(t)=|\psi^{N}|_{H^{r,s}}(t),$ (59)

where $F_{0}$ is defined in (39). Then Lemma 5 gives the following recursion
relation:

$c_{r,s}(t)\leq Cs\mathcal{A}_{0}[c_{r+1,s-1}](t)+b_{r,s}(t),\qquad s\geq 1.$
(60)

Unrolling this recursion in $s$ gives

$c_{r,s}(t)\leq
s!C^{s}\mathcal{A}^{s}_{0}[c_{r+s,0}](t)+\sum_{i=0}^{s-1}\frac{s!}{(s-i)!}C^{i}\mathcal{A}_{0}^{i}[b_{r+i,s-i}](t),$
(61)

and translating back to the semi-norms with (59) gives

$\begin{split}|\psi^{N}|_{H^{r,s}}(t)&\leq
s!C^{s}\mathcal{A}_{0}^{s}[|\psi^{N}|_{H^{r+s,0}}](t)\\\
&+\sum_{i=0}^{s-1}\frac{s!}{(s-i)!}C^{i}|\mathcal{P}_{N}g|_{H^{r+i,s-i}}\mathcal{A}_{0}^{i}[F_{0}](t)+\varepsilon\sum_{i=0}^{s-1}\frac{s!}{(s-i)!}C^{i+1}\mathcal{A}_{0}^{i+1}[|\mathcal{P}_{N}q|_{H^{r+i,s-i}}](t).\end{split}$
(62)

We then apply (54) (with $r$ replaced with $r+s$) to the first term on the
right hand side of (62).

∎

####  3.2.2 Estimates for continuous system

In this section we extend the stability results for $|\psi^{N}|_{H^{r,s}}$ to
$|\psi|_{H^{r,s}}$, using the infinite moment hierarchy in (9). These
estimates will be useful in deriving consistency estimates. For the case
$r=0$, as in the previous section, we test each equation in (9) by
$b_{\ell}{\bm{\uppsi}}_{{\ell}}$, integrate by parts, and sum over $\ell\geq
s$. The result is analogous to (45), namely

$\begin{split}\frac{\varepsilon}{2}\partial_{t}|\psi|_{H^{0,s}}^{2}&+\frac{\sigma}{\varepsilon}|\psi|_{H^{0,s}}^{2}=\varepsilon(q,\psi)_{H^{0,s}}\\\
&\quad-\sum_{i=1}^{3}\sum_{\ell=s}^{\infty}\left(\left(\ell+\frac{1}{2}\right)^{2s}-\gamma_{s,\ell}\left(\ell-\frac{1}{2}\right)^{2s}\right)\left({\bm{\uppsi}}_{{\ell}},\left(a_{\ell}^{(i)}\right)^{T}\partial_{x_{i}}{\bm{\uppsi}}_{{\ell-1}}\right)_{L^{2}(X)},\end{split}$
(63)

As before, when $s=0$, (63) recovers the usual $L^{2}$ stability result:

$\frac{\varepsilon}{2}\frac{d}{dt}\|{\psi}\|_{L^{2}(X\times\mathbb{S}^{2})}^{2}+\frac{\sigma}{\varepsilon}\|{\psi-\overline{\psi}}\|_{L^{2}(X\times\mathbb{S}^{2})}^{2}=\varepsilon(q,\psi)_{L^{2}(X\times\mathbb{S}^{2})}^{2}$
(64)

The following results are continuous analogues to Lemmas 5 and 6. Their proofs
are nearly identical, so we only give a brief summary.

######  Lemma 7 (Semi-norm recurrence for the continuous system).

Let $s\geq 1$, $q\in L^{1}([0,t];H^{r,s})$ and $g\in H^{r,s}$. Then for all
$t\in[0,T]$,

$\displaystyle|\psi|_{H^{r,s}}(t)\leq
Cs\mathcal{A}_{0}[|\psi|_{H^{r+1,s-1}}](t)+|g|_{H^{r,s}}F_{0}(t)+\varepsilon\mathcal{A}_{0}[|q|_{H^{r,s}}](t),$
(65)

where $C$ is a constant independent of the data.

###### Summary of Proof.

The proof follows the same lines and with same constants as in Lemma 5 after
changing ${\bm{\uppsi}}^{N}_{{\ell}}$ by ${\bm{\uppsi}}_{{\ell}}$ and then
taking all the sums to infinity. To generalize the proof for $r>1$, we
differentiate the system (9) in space and repeat the process.∎

######  Lemma 8 (Stability of higher order semi-norms for the continuous
system).

If $q\in L^{1}([0,T];H^{r,0})$ and $g\in H^{r,0}$, then for all $t\in[0,T]$,

$\displaystyle|\psi|_{H^{r,0}}(t)\leq|g|_{H^{r,0}}+|q|_{L^{1}([0,t];H^{r,0})}.$
(66)

If, in addition, $q\in L^{1}([0,T];H^{i,j})$ and $g\in H^{i,j}$ for each $i,j$
such that $0\leq i\leq r$, $0\leq j\leq s$, and $i+j=r+s$,

$\begin{split}|\psi|_{H^{r,s}}(t)&\leq
C_{s}s!\mathcal{A}_{0}^{s}\left[|g|_{H^{r+s,0}}+|q|_{L^{1}([0,\bullet];H^{r+s,0})}\right](t)\\\
&+C_{s}\sum_{i=0}^{s-1}\frac{s!}{(s-i)!}\mathcal{A}_{0}^{i}[F_{0}|g|_{H^{r+i,s-i}}](t)+C_{s}\varepsilon\sum_{i=0}^{s-1}\frac{s!}{(s-i)!}\mathcal{A}_{0}^{i+1}[|q|_{H^{r+i,s-i}}](t),\end{split}$
(67)

where $C_{s}$ is a constant depending only on s.

######  Corollary 1 (Isotropic data and zero initial condition for the
continuous systems).

Let $s\geq 1$. In the special case that $g=0$ and $q$ is isotropic, then if
$t\in[0,T]$,

$\displaystyle|\psi|_{H^{r,s}}(t)\leq
C_{s}s!\mathcal{A}_{0}^{s}[|q|_{L^{1}([0,\bullet];H^{r+s,0})}](t).$ (68)

###### Proof of Lemma 8.

With the obvious changes, the proof is the same line by line as proof in Lemma
6, with some key steps replace by their continuous counterparts, namely, in
the initial step we use (64) instead of (46) and we invoke Lemma 7 instead of
Lemma 5. ∎

###  3.3 $\text{P}_{N}$ error analysis

In this section we will analyze the error produced by the solution of (3) when
the $\text{P}_{N}$ approximation is used.

######  Definition 1 ($\text{P}_{N}$ Error).

The $\text{\emph{P}}_{N}$ _error_ is

$e^{N}(t)=\psi(t)-\psi^{N}(t)=\eta^{N}(t)+\xi^{N}(t),$ (69)

where $\eta^{N}=\psi-\mathcal{P}_{N}\psi$ is the _consistency error_ and
$\xi^{N}=\mathcal{P}_{N}\psi-\psi^{N}$ is the _stability error._

######  Lemma 9.

For all $t\in[0,T]$, $\xi^{N}$ is controlled by $\eta^{N}$ via the following
estimate:

$\|{\xi^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t)\leq\frac{1}{\varepsilon}\int_{0}^{t}\|{\mathcal{P}_{N}(\Omega\cdot\nabla_{x}\eta^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(\tau)\,d\tau.$
(70)

###### Proof.

Applying the projection $\mathcal{P}_{N}$ to (3) and subtracting (7a) yields a
$\text{P}_{N}$ equation for $\xi^{N}$ with a source that depends on
$\eta^{N}$:

$\varepsilon\partial_{t}\xi^{N}+\mathcal{P}_{N}(\Omega\cdot\nabla_{x}\xi^{N})+\frac{\sigma}{\varepsilon}\xi^{N}=\frac{\sigma}{\varepsilon}\overline{\xi^{N}}-\mathcal{P}_{N}(\Omega\cdot\nabla_{x}\eta^{N}),\qquad\left.\xi^{N}\right|_{t=0}=0.$
(71)

Thus (71) follows immediately from the bound (57), replacing $g$ by zero and
$q$ by ${\varepsilon^{-1}\mathcal{P}_{N}(\Omega\cdot\nabla_{x}\eta^{N})}$. ∎

######  Lemma 10.

Let $t\in[\alpha,\beta)\subseteq[0,T]$, then

$\|{\eta^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t)\leq
e^{-\frac{\sigma(t-\alpha)}{\varepsilon^{2}}}\|{\eta^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(\alpha)+\varepsilon\mathcal{A}_{\alpha}[\|{\widetilde{\mathcal{P}}_{N}q}\|_{L^{2}(X\times\mathbb{S}^{2})}](t)+\mathcal{A}_{\alpha}[\|{\widetilde{\mathcal{P}}_{N}(\Omega\cdot\nabla_{x}\psi)}\|_{L^{2}(X\times\mathbb{S}^{2})}](t)$
(72)

###### Proof.

Applying the projection $\mathcal{P}_{N}$ to (3), and subtracting it from (3),
we see that $\eta^{N}$ satisfies

$\partial_{t}\eta^{N}-\frac{1}{\varepsilon}\mathcal{P}_{N}(\Omega\cdot\nabla_{x}\eta^{N})+\frac{\sigma}{\varepsilon^{2}}\eta^{N}=\widetilde{\mathcal{P}}_{N}q-\frac{1}{\varepsilon}\left[\Omega\cdot\nabla_{x}\psi-\mathcal{P}_{N}(\Omega\cdot\nabla_{x}\mathcal{P}_{N}\psi)\right]$
(73)

Since for any $\phi\in L^{2}(\mathbb{S}^{2})$,
$(\mathcal{P}_{N}\phi,\eta^{N})_{L^{2}(\mathbb{S}^{2})}=0$, testing the
equation above against $\eta^{N}$ gives

$\begin{split}\frac{1}{2}\frac{d}{dt}\|{\eta^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}^{2}+\frac{\sigma}{\varepsilon^{2}}\|{\eta^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}^{2}&=(\widetilde{\mathcal{P}}_{N}q,\eta^{N})_{L^{2}(X\times\mathbb{S}^{2})}-\frac{1}{\varepsilon}(\Omega\cdot\nabla_{x}\psi,\eta^{N})_{L^{2}(X\times\mathbb{S}^{2})}\\\
&=(\widetilde{\mathcal{P}}_{N}q,\eta^{N})_{L^{2}(X\times\mathbb{S}^{2})}-\frac{1}{\varepsilon}(\widetilde{\mathcal{P}}_{N}(\Omega\cdot\nabla_{x}\psi),\eta^{N})_{L^{2}(X\times\mathbb{S}^{2})}\\\
&\leq(\|{\widetilde{\mathcal{P}}_{N}q}\|_{L^{2}(X\times\mathbb{S}^{2})}+\frac{1}{\varepsilon}\|{\widetilde{\mathcal{P}}_{N}(\Omega\cdot\nabla_{x}\psi)}\|_{L^{2}(X\times\mathbb{S}^{2})})\|{\eta^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}\end{split}$
(74)

The conclusion then follows from Lemma 21. ∎

An immediate corollary of Lemma 10 is the following:

######  Lemma 11.

Let $s\geq 0$ and $N\geq\max\\{0,s-1\\}$. If $g\in H^{0,s}$ and $q\in
L^{\infty}([0,T];H^{0,s})$. Then we have

$\|{\eta^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(T)\leq\frac{e^{-\sigma
T/\varepsilon^{2}}}{(N+1)^{s}}|g|_{H^{0,s}}+\frac{\varepsilon}{(N+1)^{s}}|q|_{L^{\infty}([0,T];H^{0,s})}\mathcal{A}_{0}[\mathbbm{1}](T)+\frac{T}{\varepsilon(N+1)^{s}}\sup_{\tau\in[0,T]}|\psi|_{H^{1,s}}(\tau).$
(75)

###### Proof.

We apply Lemma 10 with $\alpha=0$. In this case
$\|{\eta^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(\alpha)=\|{\widetilde{\mathcal{P}}_{N}g}\|$.
Meanwhile, by Lemma 2,

$\|{\widetilde{\mathcal{P}}_{N}g}\|_{L^{2}(X\times\mathbb{S}^{2})}\leq\frac{1}{(N+1)^{s}}|g|_{H^{0,s}}\quad\text{and}\quad\|{\widetilde{\mathcal{P}}_{N}q}\|_{L^{2}(X\times\mathbb{S}^{2})}\leq\frac{1}{(N+1)^{s}}|q|_{H^{0,s}}.$
(76)

Thus since
$\mathcal{A}_{\alpha}[f](t)\leq\varepsilon^{-1}(t-\alpha)\sup_{\tau\in[\alpha,t]}f(\tau)$
(cf. (40) with $k=1$), another application of Lemma 2 gives

$\begin{split}&\mathcal{A}_{0}[\|{\widetilde{\mathcal{P}}_{N}(\Omega\cdot\nabla_{x}\psi)}\|_{L^{2}(X\times\mathbb{S}^{2})}](T)\leq\frac{T}{\varepsilon}\sup_{\tau\in[0,T]}\|{\widetilde{\mathcal{P}}_{N}(\Omega\cdot\nabla_{x}\psi)}\|_{L^{2}}(\tau)\\\
&\qquad\qquad\leq\frac{T}{\varepsilon(N+1)^{s}}\sup_{\tau\in[0,T]}\,|\Omega\cdot\nabla_{x}\psi|_{H^{0,s}}(\tau)\leq\frac{T}{\varepsilon(N+1)^{s}}\sup_{\tau\in[0,T]}|\psi|_{H^{1,s}}(\tau)\end{split}$
(77)

Plugging the preceding bound into Lemma 10 gives the result. ∎

######  Lemma 12 (a priori estimate).

Let $s\geq 1$ and $N\geq s-1$. Let $q\in L^{\infty}([0,T];H^{i,j})$ and $g\in
H^{i,j}$ for each $i,j$ such that $0\leq i\leq r$, $0\leq j\leq s$, and
$i+j=r+s$. Then for all $t\in[0,T]$,

$\begin{split}|\psi|_{H^{r,s}}(t)&\leq
C_{s}\Bigg{\\{}\left[|g|_{H^{r+s,}}+t|q|_{L^{\infty}([0,t);H^{r+s,0})}\right]\min\left(\frac{\varepsilon^{s}s!}{\sigma^{s}},\left(\frac{t}{\varepsilon}\right)^{s}\right)+e^{-\sigma
t/\varepsilon^{2}}\sum_{i=0}^{s-1}|g|_{H^{r+i,s-i}}\binom{s}{i}\frac{t^{i}}{\varepsilon^{i}}\\\
&\quad+\varepsilon\sum_{i=0}^{s-1}|q|_{L^{\infty}([0,t];H^{r+i,s-i})}\frac{s!}{(s-i)!}\min\left(\frac{\varepsilon^{i+1}}{\sigma^{i+1}},\frac{1}{(i+1)!}\left(\frac{t}{\varepsilon}\right)^{i+1}\right)\Bigg{\\}}.\end{split}$
(78)

###### Proof.

Recall the stability estimate from Lemma 8:

$\begin{split}|\psi|_{H^{r,s}}(t)&\leq
C_{s}s!\mathcal{A}_{0}^{s}\left[|g|_{H^{r+s,0}}+|q|_{L^{1}([0,\bullet];H^{r+s,0})}\right](t)\\\
&+C_{s}\sum_{i=0}^{s-1}\frac{s!}{(s-i)!}\mathcal{A}_{0}^{i}[F_{0}|g|_{H^{r+i,s-i}}](t)+C_{s}\varepsilon\sum_{i=0}^{s-1}\frac{s!}{(s-i)!}\mathcal{A}_{0}^{i+1}[|q|_{H^{r+i,s-i}}](t).\end{split}$
(79)

Substituting the bounds for $\mathcal{A}_{0}[\mathbbm{1}]$ and the formula for
$\mathcal{A}_{0}[F_{\alpha}]$ from Lemma 3 into the above estimate yields the
stated result. ∎

######  Theorem 13 ($\text{P}_{N}$ error).

Let $s\geq 1$ and $N\geq s-1$. Let $q\in L^{\infty}([0,T];H^{i,j})$ and $g\in
H^{i,j}$ for each $i,j$ such that $0\leq j\leq s$, $i+j\leq s+1$. Then

$\begin{split}\|{e^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(T)&\leq\frac{e^{-\sigma
T/\varepsilon^{2}}}{(N+1)^{s}}|g|_{H^{0,s}}+\frac{1}{(N+1)^{s}}|q|_{L^{\infty}([0,T];H^{0,s})}\min\left(\frac{\varepsilon^{2}}{\sigma},T\right)\\\
&+\frac{2C_{s}}{(N+1)^{s}}\Bigg{\\{}\Big{(}|g|_{H^{s+1,0}}+T|q|_{L^{\infty}([0,T];H^{s+1,0})}\Big{)}\min\left(\frac{\varepsilon^{s-1}s!T}{\sigma^{s}},\left(\frac{T}{\varepsilon}\right)^{s+1}\right)\\\
&\quad+e^{-\sigma
T/\varepsilon^{2}}\sum_{i=0}^{s-1}|g|_{H^{1+i,s-i}}\binom{s}{i}\frac{T^{i+1}}{\varepsilon^{i+1}}\\\
&\quad+\sum_{i=0}^{s-1}|q|_{L^{\infty}([0,T];H^{1+i,s-i})}\frac{s!}{(s-i)!}\min\left(\frac{\varepsilon^{i+1}T}{\sigma^{i+1}},\frac{1}{(i+1)!}\frac{T^{i+2}}{\varepsilon^{i+1}}\right)\Bigg{\\}}\end{split}$
(80)

###### Proof.

By the triangle inequality, Lemma 9,

$\begin{split}\|{e^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(T)&\leq\|\eta^{N}\|_{L^{2}(X\times\mathbb{S}^{2})}(T)+\|\xi^{N}\|_{L^{2}(X\times\mathbb{S}^{2})}(T)\\\
&\leq\|\eta^{N}\|_{L^{2}(X\times\mathbb{S}^{2})}(T)+\frac{1}{\varepsilon}\int_{0}^{T}\|{\mathcal{P}_{N}(\Omega\cdot\nabla_{x}\eta^{N})}\|_{L^{2}(X\times\mathbb{S}^{2})}(\tau)\,d\tau\\\
&\leq\|\eta^{N}\|_{L^{2}(X\times\mathbb{S}^{2})}(T)+\frac{1}{\varepsilon}\int_{0}^{T}\|{\nabla_{x}\eta^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(\tau)\,d\tau\\\
&\leq\|\eta^{N}\|_{L^{2}(X\times\mathbb{S}^{2})}(T)+\frac{1}{(N+1)^{s}}\frac{T}{\varepsilon}\sup_{\tau\in[0,T]}|\psi|_{H^{1,s}}(\tau)\\\
&\leq\frac{e^{-\sigma
T/\varepsilon^{2}}}{(N+1)^{s}}|g|_{H^{0,s}}+\frac{\varepsilon}{(N+1)^{s}}|q|_{L^{\infty}([0,T];H^{0,s})}\mathcal{A}_{0}[\mathbbm{1}](T)+\frac{2T}{\varepsilon(N+1)^{s}}\sup_{\tau\in[0,T]}|\psi|_{H^{1,s}}(\tau)\end{split}$
(81)

In the last two lines, we applied spectral estimate in Lemma 2. In the last
line we used Lemma 11. Applying Lemma 12 with $r=1$ yields the result. ∎

######  Corollary 2 ($\text{P}_{N}$ error for Isotropic data).

Let $s\geq 1$ and $N\geq s-1$. Let $q\in L^{\infty}([0,T];H^{s+1,0})$ and
$g\in H^{s+1,0}$. If $g$ and $q$ are isotropic, then

$\|{e^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(T)\leq\frac{2C_{s}}{(N+1)^{s}}\Big{(}|g|_{H^{s+1,0}}+T|q|_{L^{\infty}([0,T];H^{s+1,0})}\Big{)}\min\left(\frac{\varepsilon^{s-1}s!T}{\sigma^{s}},\left(\frac{T}{\varepsilon}\right)^{s+1}\right)$
(82)

##  4 Hybrid error analysis

###  4.1 A priori estimates of the uncollided component

Since our goal is to derive error estimates which only depend on data, and
$\psi_{\mathrm{u},m}$ appears as a source in the collided equation, we require
the following a priori estimates on $\psi_{\mathrm{u},m}$ to bound
$|\psi_{\mathrm{c},m}|_{H^{1,s}}$ in the proof of Theorem 18.

######  Lemma 14 (Stability of the uncollided component).

Let $1\leq m\leq M$, $q\in L^{\infty}([0,T];H^{r,0})$, and $g\in H^{r,0}$ for
some $r\geq 0$. Then for all $t\in[t_{m-1},t_{m})$,

$\displaystyle|\psi_{\mathrm{u},m}|_{H^{r,0}}(t)\leq
e^{-\sigma(t-t_{m-1})/\varepsilon^{2}}|g|_{H^{r,0}}+e^{-\sigma(t-t_{m-1})/\varepsilon^{2}}|q|_{L^{1}([0,t_{m-1}];H^{r,0})}+\varepsilon\mathcal{A}_{t_{m-1}}[|q|_{H^{r,0}}](t),$
(83a)
$\displaystyle|\psi_{\mathrm{u},m}|_{H^{r,0}}(t_{m})+\frac{\sigma}{\varepsilon^{2}}|\psi_{\mathrm{u},m}|_{L^{1}([t_{m-1},t_{m}];H^{r,0})}\leq|g|_{H^{r,0}}+|q|_{L^{1}([0,t_{m}];H^{r,0})},$
(83b)
$\displaystyle|\psi_{\mathrm{u},m}|_{L^{1}([t_{m-1},t];H^{r,0})}(t)\leq\varepsilon(|g|_{H^{r,0}}+t|q|_{L^{\infty}([0,t];H^{r,0})})\mathcal{A}_{t_{m-1}}[\mathbbm{1}](t).$
(83c)

###### Proof.

We prove the result only for $r=0$ since the other cases are obtained by
applying the same techniques to (12a) differentiated $r$ times in space.
Testing (12a) with $\psi_{\mathrm{u},m}$ and applying Cauchy-Schwarz
inequality, we obtain the following differential inequality

$\frac{1}{2}\frac{d}{dt}\|{\psi_{\mathrm{u},m}}\|^{2}+\frac{\sigma}{\varepsilon^{2}}\|{\psi_{\mathrm{u},m}}\|^{2}\leq\|{q}\|\|{\psi_{\mathrm{u},m}}\|.$
(84)

An application of (110) in Lemma 21, and the fact that
$\psi_{\mathrm{u},m}(t_{m-1})=\psi(t_{m-1})$ gives

$\displaystyle\|{\psi_{\mathrm{u},m}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t)\leq
e^{-\sigma(t-t_{m-1})/\varepsilon^{2}}\|{\psi}\|_{L^{2}(X\times\mathbb{S}^{2})}(t_{m-1})+\varepsilon\mathcal{A}_{t_{m-1}}[\|q\|_{L^{2}(X\times\mathbb{S}^{2})}](t).$
(85)

Using Lemma 8 on $\|{\psi}\|_{L^{2}(X\times\mathbb{S}^{2})}(t_{m-1})$ gives
the first result (83a). An application of (109) in Lemma 21 over
$[t_{m-1},t_{m})$ gives

$\|{\psi_{\mathrm{u},m}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t_{m})+\frac{\sigma}{\varepsilon^{2}}\|\psi_{\mathrm{u},m}\|_{L^{1}([t_{m-1},t_{m}],L^{2}(X\times\mathbb{S}^{2}))}\leq\|q\|_{L^{1}([t_{m-1},t_{m}],L^{2}(X\times\mathbb{S}^{2}))}+\|{\psi}\|_{L^{2}(X\times\mathbb{S}^{2})}(t_{m-1}),$
(86)

and then another application of Lemma 8 to bound
$\|{\psi}\|_{L^{2}(X\times\mathbb{S}^{2})}(t_{m-1})$ gives (83b). Estimate
(83c) is obtained from integrating the first estimate (83a) from $t_{m-1}$ to
$t$. ∎

###  4.2 Hybrid error estimates

In this section we will analyze the error in the hybrid method using the
formulation in (14).

######  Definition 2 (Hybrid errors).

Let $1\leq m\leq M$ and $t\in[t_{m-1},t_{m})$. The $m$-th _hybrid error_ is

$e_{m}^{N}(t)=e_{\mathrm{u},m}^{N}(t)+e_{\mathrm{c},m}^{N}(t),$ (87)

where
$e_{\mathrm{u},m}^{N}(t)=\psi_{\mathrm{u},m}(t)-\psi_{\mathrm{u},m}^{N}(t)$
and
$e_{\mathrm{c},m}^{N}(t)=\psi_{\mathrm{c},m}(t)-\psi_{\mathrm{c},m}^{N}(t)$
are the $m$-th errors in the uncollided and collided components. The collided
error can be further decomposed as

$\eta_{\mathrm{c},m}^{N}(t)=\psi_{\mathrm{c},m}(t)-\mathcal{P}_{N}\psi_{\mathrm{c},m}(t)\qquad\text{and}\qquad\xi_{\mathrm{c},m}^{N}=\mathcal{P}_{N}\psi_{\mathrm{c},m}(t)-\psi_{\mathrm{c},m}^{N}(t)$
(88)

so that
$e_{\mathrm{c},{m}}^{N{}}(t)=\eta_{\mathrm{c},m}^{N}(t)+\xi_{\mathrm{c},m}^{N}(t).$
Here $\eta^{N}_{\mathrm{c},m}$ is the _ $m$-th collided consistency error_ and
$\xi^{N}_{\mathrm{c},m}$ is the _$m$ -th collided stability error_. The error
$e^{N}_{M}$ is simply called the _hybrid error_.

This next lemma gives a one-step analysis of the growth of the error in the
uncollided and collided components from $t_{m-1}$ to $t_{m}$.

######  Lemma 15.

Let $1\leq m\leq M$, then if $t\in[t_{m-1},t_{m})$, the $m$-th uncollided and
collided errors satisfy, respectively,

$\displaystyle\|{e_{\mathrm{u},m}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t)$
$\displaystyle\leq
e^{-\sigma(t-t_{m-1)}/\varepsilon^{2}}\|{e_{{m-1}}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t^{-}_{m-1}),$
(89a)
$\displaystyle\|{e_{\mathrm{c},m}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t)$
$\displaystyle\leq\left(1-e^{-\sigma(t-t_{m-1)}/\varepsilon^{2}}\right)\|{e_{{m-1}}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t^{-}_{m-1})$
$\displaystyle\quad+\|{\eta_{\mathrm{c},m}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t)+\frac{1}{\varepsilon}\|{\Omega\cdot\nabla_{x}\eta_{\mathrm{c},m}^{N}}\|_{L^{1}([t_{m-1},t];L^{2}(X\times\mathbb{S}^{2}))}.$
(89b)

###### Proof.

Subtracting (14a) from (12a), yields the following evolution equation for
$e^{N}_{\mathrm{u},m}$,

$\displaystyle\varepsilon\partial_{t}e^{N}_{\mathrm{u},m}+\Omega\cdot\nabla_{x}e^{N}_{\mathrm{u},m}+\frac{\sigma}{\varepsilon}e^{N}_{\mathrm{u},m}=0,\qquad\left.e^{N}_{\mathrm{u},m}\right|_{t=t_{m-1}}=e^{N}_{m-1}(t_{m-1}^{-}),$
(90)

where $e^{N}_{0}(t_{0}^{-})=0$. Thus applying (83a) from Lemma 14, with $r=0$
and a zero source term yields (89a). To prove (89b), we subtract from (14b)
the projection applied to (12b). This gives the following $\text{P}_{N}$
equation for $\xi_{\mathrm{c},m}^{N}$

$\displaystyle\varepsilon\partial_{t}\xi_{\mathrm{c},m}^{N}+\mathcal{P}_{N}(\Omega\cdot\nabla_{x}\xi_{\mathrm{c},m}^{N})+\dfrac{\sigma}{\varepsilon}\xi_{\mathrm{c},m}^{N}=\dfrac{\sigma}{\varepsilon}(\overline{\xi_{\mathrm{c},m}^{N}}+\overline{e_{\mathrm{u},m}^{N}})-\mathcal{P}_{N}(\Omega\cdot\nabla_{x}\eta_{\mathrm{c},m}^{N}),$
(91a) $\displaystyle\left.\xi_{\mathrm{c},m}^{N}\right|_{t=t_{m-1}}=0.$ (91b)

We apply Lemma 6 to (91a) with zero initial data and source
$\varepsilon^{-2}\sigma\overline{e_{\mathrm{u},m}^{N}}-\varepsilon^{-1}\mathcal{P}_{N}(\Omega\cdot\nabla_{x}\eta_{\mathrm{c},m}^{N})$.
Combined with bound (89a), the estimate on $\xi_{\mathrm{c},m}^{N}$ becomes

$\begin{split}\|{\xi_{\mathrm{c},m}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t)&\leq\frac{\sigma}{\varepsilon^{2}}\int_{t_{m-1}}^{t}\|{e_{\mathrm{u},m}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(\tau)\,d\tau+\frac{1}{\varepsilon}\int_{t_{m-1}}^{t}\|{\Omega\cdot\nabla_{x}\eta_{\mathrm{c},m}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(\tau)\,d\tau\\\
&\leq\frac{\sigma}{\varepsilon^{2}}\|{e_{{m-1}}^{N}}\|(t_{m-1}^{-})\int_{t_{m-1}}^{t}e^{-\sigma(\tau-
t_{m-1})/\varepsilon^{2}}\,d\tau+\frac{1}{\varepsilon}\int_{t_{m-1}}^{t}\|{\Omega\cdot\nabla_{x}\eta_{\mathrm{c},m}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(\tau)\,d\tau\\\
&=\left(1-e^{-\sigma(t-t_{m-1)}/\varepsilon^{2}}\right)\|{e_{{m-1}}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t_{m-1}^{-})+\frac{1}{\varepsilon}\int_{t_{m-1}}^{t}\|{\Omega\cdot\nabla_{x}\eta_{\mathrm{c},m}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(\tau)\,d\tau.\end{split}$
(92)

Adding $\|{\eta_{\mathrm{c},m}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t)$ to
the both sides recovers (89b).∎

Now we can state an error for all time that only depends on the approximation
properties of the spherical harmonic discretization on the solution.

######  Theorem 16.

The hybrid error $e_{M}^{N}$ satisfies

$\displaystyle\|{e_{M}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t_{M}^{-})\leq\sum_{m=1}^{M}\left(\|{\eta_{\mathrm{c},m}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t_{m}^{-})+\frac{1}{\varepsilon}\|{\Omega\cdot\nabla_{x}\eta_{\mathrm{c},m}^{N}}\|_{L^{1}([t_{m-1},t_{m}];L^{2}(X\times\mathbb{S}^{2}))}\right).$
(93)

###### Proof.

Adding the inequalities in Lemma 15 and taking the limit $t\to t_{m+1}^{-}$
gives

$\begin{split}\|{e_{M}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t^{-}_{M})\leq\|{e_{{M-1}}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t^{-}_{M-1})&+\|{\eta_{\mathrm{c},{M}}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t^{-}_{M})\\\
&+\frac{1}{\varepsilon}\|{\Omega\cdot\nabla_{x}\eta_{\mathrm{c},{M}}^{N}}\|_{L^{1}([t_{M-1},t_{M}];L^{2}(X\times\mathbb{S}^{2}))}\end{split}$
(94)

Exhausting this recursion until $e_{0}^{N}(t^{-}_{0})=0$ yields the result. ∎

######  Lemma 17.

Let $s\geq 0$, $N\geq\max\\{0,s-1\\}$, and $1\leq m\leq M$. Then the $m$-th
projection error $\eta^{N}_{\mathrm{c},m}$ satisfies,

$\|{\eta_{\mathrm{c},m}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t^{-}_{m})\leq\frac{\Delta
t}{\varepsilon(N+1)^{s}}\sup_{\tau\in[t_{m-1},t_{m}]}\left|\psi_{\mathrm{c},m}\right|_{H^{1,s}}(\tau).$
(95)

###### Proof.

An application of Lemma 10 with $\alpha=t_{m-1}$, $\beta=t_{m}$,
$\eta_{\mathrm{c},m}^{N}(t_{m-1})=0$, an isotropic source
$q=\frac{\sigma}{\varepsilon^{2}}\overline{\psi_{\mathrm{u},m}}$, along with
the fact that
$\mathcal{A}_{\alpha}[f](t)\leq\varepsilon^{-1}(t-\alpha)\sup_{\tau\in[\alpha,t]}f(\tau)$
(cf. (40) with $k=1$)) and the spectral estimate in Lemma 2, gives

$\begin{split}&\|{\eta_{\mathrm{c},m}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t_{m}^{-})\leq\frac{\Delta
t}{\varepsilon}\sup_{\tau\in[t_{m-1},t_{m}]}\|{\widetilde{\mathcal{P}}_{N}(\Omega\cdot\nabla_{x}\psi_{\mathrm{c},m}))}\|_{L^{2}(X\times\mathbb{S}^{2})}(\tau)\,d\tau\\\
&\quad\leq\frac{\Delta
t}{\varepsilon(N+1)^{s}}\sup_{\tau\in[t_{m-1},t_{m}]}\left|\Omega\cdot\nabla_{x}\psi_{\mathrm{c},m}\right|_{H^{0,s}}(\tau)\,d\tau\leq\frac{\Delta
t}{\varepsilon(N+1)^{s}}\sup_{\tau\in[t_{m-1},t_{m}]}\left|\psi_{\mathrm{c},m}\right|_{H^{1,s}}(\tau).\end{split}$
(96)

∎

###  4.3 Estimating hybrid error in terms of the data

Finally, we will apply the approximation properties and stability estimates to
the estimate in Theorem 16 to obtain an estimate that depends only on the
regularity of the data.

######  Theorem 18.

Let $s\geq 1$ and $N\geq s-1$. If $q\in L^{1}([0,T];H^{s+1,0})$ and $g\in
H^{s+1,0}$, then

$\displaystyle\|{e_{M}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(T^{-})\leq\frac{2C_{s}}{(N+1)^{s}}\left(\left|g\right|_{H^{s+1,0}}+T\left|q\right|_{L^{\infty}([0,T];H^{s+1,0})}\right)\min\left(\frac{\varepsilon^{s-1}s!T}{\sigma^{s}},\frac{\Delta
t^{s}T}{\varepsilon^{s+1}}\min\left(1,\frac{\Delta
t\sigma}{\varepsilon^{2}}\right)\right).$ (97)

###### Proof.

Using
$\|\Omega\cdot\nabla_{x}\psi_{\mathrm{c},m}\|_{L^{2}(X\times\mathbb{S}^{2})}\leq\|\nabla_{x}\psi_{\mathrm{c},m}\|_{L^{2}(X\times\mathbb{S}^{2})}$,
and applying Lemma 2 and Lemma 17 to Theorem 16 yields, for $N\geq s-1$,

$\displaystyle\|{e_{M}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t^{-}_{M})$
$\displaystyle\leq\frac{1}{(N+1)^{s}}\sum_{m=1}^{M}\left(\frac{\Delta
t}{\varepsilon}\sup_{\tau\in[t_{m-1},t_{m}]}\left|\psi_{\mathrm{c},m}\right|_{H^{1,s}}(\tau)+\frac{1}{\varepsilon}|\psi_{\mathrm{c},m}|_{L^{1}([t_{m-1},t_{m}];H^{1,s})}\right)$
$\displaystyle\leq\frac{2\Delta
t}{\varepsilon(N+1)^{s}}\sum_{m=1}^{M}\sup_{\tau\in[t_{m-1},t_{m}]}|\psi_{\mathrm{c},m}|_{H^{1,s}}(\tau)$

Applying Corollary 1 with
$q=\frac{\sigma}{\varepsilon^{2}}\overline{\psi_{\mathrm{u},m}}$ and $r=1$,

$\|{e_{M}^{N}(t^{-}_{M})}\|_{L^{2}(X\times\mathbb{S}^{2})}\leq
2C_{s}s!\frac{\sigma}{\varepsilon^{3}}\frac{\Delta
t}{(N+1)^{s}}\sum_{m=1}^{M}\sup_{\tau\in[t_{m-1},t_{m}]}\mathcal{A}_{t_{m-1}}^{s}\left[|\psi_{\mathrm{u},m}|_{L^{1}([t_{m-1},\bullet];H^{s+1,0})}\right](\tau).$
(98)

For the summand above, it follows from (83b), the monotonicity of
$\mathcal{A}_{\alpha}$, and Lemma 3 that

$\begin{split}\sup_{\tau\in[t_{m-1},t_{m}]}\mathcal{A}_{t_{m-1}}^{s}\left[|\psi_{\mathrm{u},m}|_{L^{1}([t_{m-1},\bullet];H^{s+1,0})}\right](\tau)&\leq\sup_{\tau\in[t_{m-1},t_{m}]}|\psi_{\mathrm{u},m}|_{L^{1}([t_{m-1},t_{m}];H^{s+1,0})}\mathcal{A}_{t_{m-1}}^{s}\left[\mathbbm{1}\right](\tau)\\\
&\leq\frac{\varepsilon^{2}}{\sigma}\left(\left|g\right|_{H^{s+1,0}}+\left|q\right|_{L^{1}([0,T];H^{s+1,0})}\right)\mathcal{A}_{t_{m-1}}^{s}\left[\mathbbm{1}\right](t_{m})\\\
&\leq\frac{\varepsilon^{2}}{\sigma}\left(\left|g\right|_{H^{s+1,0}}+\left|q\right|_{L^{1}([0,T];H^{s+1,0})}\right)\min\left(\frac{\varepsilon^{s}}{\sigma^{s}},\frac{1}{s!}\left(\frac{\Delta
t}{\varepsilon}\right)^{s}\right).\end{split}$ (99)

Pluggin the above bound into (98) yields (since $T=M\Delta t$)

$\begin{split}\|{e_{M}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t^{-}_{M})&\leq
2C_{s}s!\frac{\Delta
tM}{\varepsilon(N+1)^{s}}\left(\left|g\right|_{H^{s+1,0}}+\left|q\right|_{L^{1}([0,T];H^{s+1,0})}\right)\min\left(\frac{\varepsilon^{s}}{\sigma^{s}},\frac{\Delta
t^{s}}{s!\varepsilon^{s}}\right)\\\
&=\frac{2C_{s}}{(N+1)^{s}}\left(\left|g\right|_{H^{s+1,0}}+\left|q\right|_{L^{1}([0,T];H^{s+1,0})}\right)\min\left(\frac{\varepsilon^{s-1}s!T}{\sigma^{s}},\frac{\Delta
t^{s}T}{\varepsilon^{s+1}}\right).\end{split}$ (100)

On the other hand, it follows from (83c) that

$\begin{split}&\sup_{\tau\in[t_{m-1},t_{m}]}\mathcal{A}_{t_{m-1}}^{s}\left[|\psi_{\mathrm{u},m}|_{L^{1}([t_{m-1},\bullet];H^{s+1,0})}\right](\tau)\\\
&\qquad\qquad\leq\varepsilon(\left|g\right|_{H^{s+1,0}}+T\left|q\right|_{L^{\infty}([0,T];H^{s+1,0})})\sup_{\tau\in[t_{m-1},t_{m}]}\mathcal{A}_{t_{m-1}}^{s+1}[\mathbbm{1}](\tau)\\\
&\qquad\qquad\leq\varepsilon(\left|g\right|_{H^{s+1,0}}+T\left|q\right|_{L^{\infty}([0,T];H^{s+1,0})})\min\left(\frac{\varepsilon^{s+1}}{\sigma^{s+1}},\frac{1}{(s+1)!}\left(\frac{\Delta
t}{\varepsilon}\right)^{s+1}\right).\end{split}$ (101)

Plugging this into (98) yields,

$\displaystyle\|{e_{M}^{N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(t^{-}_{M})\leq\frac{2C_{s}}{(N+1)^{s}}\left(\left|g\right|_{H^{s+1,0}}+T\left|q\right|_{L^{\infty}([0,T);H^{s+1,0})}\right)\min\left(\frac{\varepsilon^{s-1}s!T}{\sigma^{s}},\frac{\Delta
t^{s+1}\sigma T}{\varepsilon^{s+3}}\right).$ (102)

Taking a minimum of the right hand sides of (100) and (102) yields the result.
∎

##  5 Return to the original transport model

In this section we will show error estimates for the model (1). The analogous
discretizations for the models are the following. For the non-splitting
$\text{P}_{N}$ discretization we seek $\Psi^{\varepsilon,N}\in
C([0,T);X\times\mathbb{P}_{N}(\mathbb{S}^{2}))$, satisfying

$\displaystyle\varepsilon\partial_{t}\Psi^{\varepsilon,N}+\mathcal{P}_{N}(\Omega\cdot\nabla_{x}\Psi^{\varepsilon,N})+\frac{\sigma_{\rm{t}}}{\varepsilon}\Psi^{\varepsilon,N}=\left(\frac{\sigma_{\rm{t}}}{\varepsilon}-\varepsilon\sigma_{\rm{a}}\right)\overline{\Psi^{\varepsilon,N}}+\varepsilon\mathcal{P}_{N}Q,$
(103a) $\displaystyle\Psi^{\varepsilon,N}|_{t=0}=\mathcal{P}_{N}g$ (103b)

and for the hybrid, we seek
$\Psi_{m}^{\varepsilon,N}=\Psi_{\mathrm{u},m}^{\varepsilon,N}+\Psi_{\mathrm{c},m}^{\varepsilon,N}$
where for each $m\in\\{1,2,\dots,M\\}$

$\left(\Psi_{\mathrm{u},m}^{\varepsilon,N},\Psi_{\mathrm{c},m}^{\varepsilon,N}\right)\in
C([t_{m-1},t_{m});X\times L^{2}(\mathbb{S}^{2}))\times
C([t_{m-1},t_{m});X\times\mathbb{P}_{N}(\mathbb{S}^{2}))$ (104)

satifies

$\displaystyle\varepsilon\partial_{t}\Psi_{{\mathrm{u}},m}^{\varepsilon,N}+\Omega\cdot\nabla_{x}\Psi_{{\mathrm{u}},m}^{\varepsilon,N}+\frac{\sigma_{\rm{t}}}{\varepsilon}\Psi_{{\mathrm{u}},m}^{\varepsilon,N}=\varepsilon
Q,$ (105a)
$\displaystyle\varepsilon\partial_{t}\Psi_{\mathrm{c},m}^{\varepsilon,N}+\mathcal{P}_{N}\left(\Omega\cdot\nabla_{x}\Psi_{\mathrm{c},m}^{\varepsilon,N}\right)+\frac{\sigma_{\rm{t}}}{\varepsilon}\Psi_{\mathrm{c},m}^{\varepsilon,N}=\left(\frac{\sigma_{\rm{t}}}{\varepsilon}-\varepsilon\sigma_{\rm{a}}\right)\overline{\Psi_{\mathrm{u},m}^{\varepsilon,N}}+\overline{\Psi_{\mathrm{c},m}^{\varepsilon,N}},$
(105b)
$\displaystyle\left.\Psi_{\mathrm{c},m}^{\varepsilon,N}\right|_{t=t_{m-1}}=0,\quad\left.\Psi_{\mathrm{u},m}^{\varepsilon,N}\right|_{t=t_{m-1}}=\begin{cases}g,&\quad
m=1,\\\
\Psi_{\mathrm{u},m-1}^{\varepsilon,N}(t^{-}_{m-1})+\Psi_{\mathrm{c},m-1}^{\varepsilon,N}(t^{-}_{m-1})&\quad
m>1.\end{cases}$ (105d)

We define the correspondent $\text{P}_{N}$ error and $m$-th hybrid error
respectively as,

$e^{\varepsilon,N}=\Psi^{\varepsilon}-\Psi^{\varepsilon,N}\quad\text{and}\quad
e^{\varepsilon,N}_{m}=\Psi_{m}^{\varepsilon}-\Psi^{\varepsilon,N}_{m}$

Since $\psi^{N}=e^{\sigma_{\rm{a}}t}\Psi^{\varepsilon,N}$ and
$\psi_{m}^{N}=e^{\sigma_{\rm{a}}t}\Psi_{m}^{\varepsilon,N}$, applying Theorems
13 and 18 gives the following estimate for the $\text{P}_{N}$ error for the
original transport model (1).

######  Theorem 19.

Let $s\geq 1$ and $N\geq s-1$, $Q\in L^{\infty}([0,T];H^{i,j})$ and $g\in
H^{i,j}$ for each $i,j$ such that $0\leq j\leq s$, $i+j\leq s+1$. Then

$\begin{split}\|{e^{\varepsilon,N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(T)&\leq\frac{e^{-(\sigma_{\rm{t}}+\varepsilon^{2}\sigma_{\rm{a}})T/\varepsilon^{2}}}{(N+1)^{s}}|g|_{H^{0,s}}+\frac{1}{(N+1)^{s}}|Q|_{L^{\infty}([0,T];H^{0,s})}\min\left(\frac{\varepsilon^{2}}{\sigma_{\rm{t}}},T\right)\\\
+2C_{s}\frac{1}{(N+1)^{s}}&\Bigg{(}\left[e^{-\sigma_{\rm{a}}T}|g|_{H^{s+1,0}}+T|Q|_{L^{\infty}([0,T];H^{s+1,0})}\right]\min\left(\frac{\varepsilon^{s-1}s!T}{\sigma_{\rm{t}}^{s}},\left(\frac{T}{\varepsilon}\right)^{s+1}\right)\\\
&\quad+e^{-(\sigma_{\rm{t}}+\varepsilon^{2}\sigma_{\rm{a}})T}\sum_{i=0}^{s-1}|g|_{H^{1+i,s-i}}\binom{s}{i}\frac{T^{i+1}}{\varepsilon^{i+1}}\\\
&+\sum_{i=0}^{s-1}|Q|_{L^{\infty}([0,T];H^{1+i,s-i})}\binom{s}{i}\min\left(\frac{\varepsilon^{i+1}T}{\sigma_{\rm{t}}^{i+1}},\frac{1}{(i+1)!}\frac{T^{i+2}}{\varepsilon^{i+1}}\right)\Bigg{)}\end{split}$
(106)

Meanwhile for the hybrid approximation.

######  Theorem 20.

Let $s\geq 1$ and $N\geq s-1$, $Q\in L^{1}([0,T];H^{s+1,0})$ and $g\in
H^{s+1,0}$, we have

$\|{e_{M}^{\varepsilon,N}}\|_{L^{2}(X\times\mathbb{S}^{2})}(T^{-})\leq\frac{2C_{s}}{(N+1)^{s}}\left(e^{-\sigma_{\rm{a}}T}\left|g\right|_{H^{s+1,0}}+T\left|Q\right|_{L^{\infty}([0,T];H^{s+1,0})}\right)\\\
\times\min\left(\frac{\varepsilon^{s-1}s!T}{\sigma_{\rm{t}}^{s}},\frac{\Delta
t^{s}T}{\varepsilon^{s+1}}\min\left(1,\frac{\Delta
t\sigma_{\rm{t}}}{\varepsilon^{2}}\right)\right).$ (107)

##  6 Conclusion

In this paper, we have derived multiscale error estimates for the
$\text{P}_{N}$ approximation of the RTE and for a hybrid approximation for the
RTE that is built using the $\text{P}_{N}$ approximation. By construction, the
hybrid is is more expensive; we use these error estimates to understand the
benefits of the additional expense for different parameter regimes. At each
time step in the hybrid approximation, the collided equation is equipped with
isotropic initial conditions and zero initial condition. In scattering
dominating regimes, this property is key to improved estimates over the
monolithic $\text{P}_{N}$ approach. Meanwhile, in purely absorbing regimes,
the hybrid captures the RTE solution exactly.

In the future work, we intend to revisit the current analysis for more general
problems on non-periodic domains, with non-constant cross-sections and inflow
boundary conditions. In addition, we intend to explicitly examine the effects
of angular discretization errors in the treatment of the uncollided equation,
which for the purposes of the current paper was assumed to be solved exactly.

##  Appendix A Other estimates

######  Lemma 21.

Assume $\chi$ is a non-negative continuous function on $[\alpha,\beta]$.
Assume $\phi\in C^{1}([\alpha,\beta])$, $\phi\geq 0$, and satisfies the
following differential inequality

$\frac{1}{2}(\phi(t)^{2})^{\prime}+\kappa\phi(t)^{2}\leq\chi(t)\phi(t),\quad\phi(\alpha)=\phi_{\alpha}\geq
0.$ (108)

Then for all $t\in[\alpha,\beta]$,

$\phi(t)+\kappa\int_{\alpha}^{t}\phi(\tau)\,d\tau\leq\phi_{\alpha}+\int_{\alpha}^{t}\chi(\tau)\,d\tau.$
(109)

Furthermore

$\phi(t)\leq
e^{-\kappa(t-\alpha)}\phi_{\alpha}+\int_{\alpha}^{t}e^{-\kappa(t-\tau)}\chi(\tau)\,d\tau.$
(110)

###### Proof.

We prove first (109). Since $\phi$ and $\chi$ are non-negative functions, it
follows that for any arbitrary $\delta>0$, the following differential
inequality holds

$\phi(t)\phi^{\prime}(t)+\kappa\phi(t)^{2}\leq\chi(t)(\phi(t)+\delta),$ (111)

dividing both sides of the inequality by $\phi+\delta$, and integrating in
time, we arrive at

$\displaystyle\phi(t)+\kappa\int_{\alpha}^{t}\phi(\tau)\,d\tau$
$\displaystyle\leq\phi_{\alpha}+\int_{\alpha}^{t}\chi(\tau)\,d\tau+\delta\ln\left|\frac{\phi(t)+\delta}{\phi_{\alpha}+\delta}\right|+\kappa\delta\int_{\alpha}^{t}\frac{\phi(\tau)}{\phi(\tau)+\delta}\,d\tau,$
$\displaystyle\leq\phi_{\alpha}+\int_{\alpha}^{t}\chi(\tau)\,d\tau+\delta\ln\left|\frac{\phi(t)+\delta}{\phi_{\alpha}+\delta}\right|+\kappa\delta(t-\alpha),$

the conclusion follows taking $\delta\to 0^{+}$.

We next prove (110). When $\kappa=0$ the result follows immediately from
(109):

$\phi(t)\leq\phi_{\alpha}+\int_{\alpha}^{t}\chi(\tau)\,d\tau.$ (112)

For the general case we multiply (108) by $e^{2\kappa t}$, obtaining

$\frac{1}{2}\left[\Phi(t)^{2}\right]^{\prime}\leq e^{\kappa t}\chi(t)\Phi(t),$
(113)

where $\Phi(t)=e^{\kappa t}\phi(t)$. Applying (112) to $\Phi$ and undoing the
transformation yields (110). ∎

## References

  * [1] Marvin L Adams and Edward W Larsen.  Fast iterative methods for discrete-ordinates particle transport calculations.  Progress in nuclear energy, 40(1):3–159, 2002.
  * [2] Raymond E Alcouffe.  A first collision source method for coupling monte carlo and discrete ordinates for localized source problems.  In Monte-Carlo Methods and Applications in Neutronics, Photonics and Statistical Physics: Proceedings of the Joint Los Alamos National Laboratory-Commissariat à l’Energie Atomique Meeting Held at Cadarache Castle, Provence, France April 22–26, 1985, pages 352–366. Springer, 2006.
  * [3] Kendall Atkinson and Weimin Han.  Spherical harmonics and approximations on the unit sphere: an introduction, volume 2044.  Springer Science & Business Media, 2012.
  * [4] Alain Bensoussan, Jacques L Lions, and George C Papanicolaou.  Boundary layers and homogenization of transport processes.  Publications of the Research Institute for Mathematical Sciences, 15(1):53–157, 1979.
  * [5] K.M. Case and P.F. Zweifel.  Linear Transport Theory.  Addison-Wesley series in nuclear engineering. Addison-Wesley Publishing Company, 1967.
  * [6] Luis Chacon, Guangye Chen, Dana A Knoll, C Newman, H Park, William Taitano, Jeff A Willert, and Geoffrey Womeldorff.  Multiscale high-order/low-order (holo) algorithms and applications.  Journal of Computational Physics, 330:21–45, 2017.
  * [7] Zheng Chen and Cory Hauck.  Multiscale convergence properties for spectral approximations of a model kinetic equation.  Mathematics of Computation, 88(319):2257–22 93, 2019.
  * [8] Michael M Crockatt, Andrew J Christlieb, C Kristopher Garrett, and Cory D Hauck.  An arbitrary-order, fully implicit, hybrid kinetic solver for linear radiative transport using integral deferred correction.  Journal of Computational Physics, 346:212–241, 2017.
  * [9] Michael M Crockatt, Andrew J Christlieb, C Kristopher Garrett, and Cory D Hauck.  Hybrid methods for radiation transport using diagonally implicit runge–kutta and space–time discontinuous galerkin time integration.  Journal of Computational Physics, 376:455–477, 2019.
  * [10] Michael M Crockatt, Andrew J Christlieb, and Cory D Hauck.  Improvements to a class of hybrid methods for radiation transport: Nyström reconstruction and defect correction methods.  Journal of Computational Physics, 422:109765, 2020.
  * [11] Feng Dai and Yuan Xu.  Approximation Theory and Harmonic Analysis on Spheres and Balls.  Springer, 01 2013.
  * [12] Robert Dautray and Jacques-Louis Lions.  Mathematical analysis and numerical methods for science and technology: volume 6 evolution problems II, volume 6.  Springer Science & Business Media, 1999.
  * [13] Martin Frank, Cory Hauck, and Kerstin Küpper.  Convergence of filtered spherical harmonic equations for radiation transport.  Communications in Mathematical Sciences, 14(5):1443–1465, 2016\.
  * [14] Cory D Hauck and Robert B Lowrie.  Temporal regularization of the $\text{P}_{N}$ equations.  Multiscale Modeling & Simulation, 7(4):1497–1524, 2009.
  * [15] Cory D Hauck and Ryan G McClarren.  A collision-based hybrid method for time-dependent, linear, kinetic transport equations.  Multiscale Modeling & Simulation, 11(4):1197–1227, 2013.
  * [16] Vincent Heningburg and Cory D Hauck.  Hybrid solver for the radiative transport equation using finite volume and discontinuous galerkin.  arXiv preprint arXiv:2002.02517, 2020.
  * [17] Edward W Larsen and Joseph B Keller.  Asymptotic solution of neutron transport problems for small mean free paths.  Journal of Mathematical Physics, 15(1):75–81, 1974.
  * [18] Mohammed Lemou and Luc Mieussens.  A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit.  SIAM Journal on Scientific Computing, 31(1):334–368, 2008.
  * [19] Elmer Eugene Lewis and Warren F Miller.  Computational methods of neutron transport.  John Wiley and Sons, Inc., New York, NY, 1984.
  * [20] Gerald C Pomraning.  The equations of radiation hydrodynamics.  Courier Corporation, 2005.
  * [21] James S Warsa, Todd A Wareing, and Jim E Morel.  Krylov iterative methods and the degraded effectiveness of diffusion synthetic acceleration for multidimensional sn calculations in problems with material discontinuities.  Nuclear science and engineering, 147(3):218–248, 2004.
  * [22] Ben Whewell, Ryan G McClarren, Cory D Hauck, and Minwoo Shin.  Multigroup Neutron Transport Using a Collision-Based Hybrid Method.  Nuclear science and engineering, 2023.