\section{Introduction}
The simplified spherical harmonics (SP$_N$) approximation was initially developed in the 1960s to reduce the number of degrees of freedom required to solve transport problems in multiple dimensions using moment-based methods such as the spherical harmonics (P$_N$) expansion \cite{Gelbard1960:SPn,Gelbard1961:SPn,Gelbard1962:SPn}. Initially, the SP$_N$\  approximation was ``derived" by Gelbard in an ad-hoc way by manipulating the P$_N$ expansion from 1-D slab geometry to have a 3-D form.
The theoretical study of the SP$_N$\  equations underwent a great awakening in the 1990s. Larsen, Morel and McGhee \cite{Larsen1993:SPn} made an asymptotic derivation of the SP$_N$\  equations and showed that the accuracy of the approximation is dependent upon the scattering ratio of the system. Subsequent works extended the asymptotic analysis to anisotropic or time-dependent problems \cite{Larsen1996:SPn,Tomavsevic1996:SP2,Brantley2000:SP3,Frank2007:SPn,Olbrant2013:SPn}. However, the use of a Neumann series expansion to the unbounded streaming operator makes the asymptotic analysis a formal proof.

Another way to the SP$_N$\  equations follows the signposts of the solid harmonics. By generalizing the angular domain from $\mathbb{S}^2$ (i.e., the unit sphere) to $\mathbb{R}^3$ Ackroyd et al. used properties of solid harmonics to introduce additional harmonics, perform Galerkin projection, restrict the relaxed equations back to $\mathbb{S}^2$, and thereby derive the SP$_N$\  equations for steady-state, anisotropic scattering problems in homogeneous media \cite{Ackroyd1999:IsotropicSHPn,Ackroyd1999:AnisotropicSHPn}. Later, Hanu\v{s} used solid harmonics in tensor form and derived the SP$_3$ equations for steady-state isotropic-scattering problems \cite{Hanus2014:MCPn}; Chao was also able to use the solid harmonics approach to show that Pomraning's angular flux trial space, which can lead to the SP$_N$\  equations \cite{Pomraning1993:SPn}, is a particular solution to the steady-state isotropic-scattering transport equation \cite{Chao2016:SPn}.

In this contribution, we would extend the solid harmonics derivation to time-dependent anisotropic-scattering homogeneous-media problems and numerically demonstrate the equivalence between SP$_N$\  and P$_N$\  equations. As such, we extend, albeit incrementally, the domain of applicability of the SP$_N$\ \ equations to time-dependent, homogeneous problems with arbitrary scattering, a result that does not appear in the extant literature.



\section{Solid Harmonics for Time-Dependent Transport Problems}
As introduced in the previous section, the derivation below is based on solid harmonics which were also used by Reference\ \cite{Ackroyd1999:IsotropicSHPn} and \cite{Ackroyd1999:AnisotropicSHPn}. Some of the properties of the solid harmonics that we use can be found in the papers, but for clarity in exposition are repeated herein.



\subsection{Expansion of the Transport Equation}\label{SPnDerivationSection:TheGroupedPnExpansionOfTransportEquation}
For monoenergetic neutron transport in homogeneous media, we begin with the equation \cite{Bell_Glasstone,2007PhDT}
\begin{equation}
    \label{TransportEquation}
    \frac{1}{v}\frac{\partial}{\partial t}\psi(\vec{r},\hat{\Omega},t)+\hat{\Omega}\cdot\nabla\psi(\vec{r},\hat{\Omega},t)+\sigma_t\psi(\vec{r},\hat{\Omega},t)=\oint_{\mathbb{S}^2}\psi(\vec{r},\hat{\Omega}',t)\sigma_s(\hat{\Omega}\cdot\hat{\Omega}')d\Omega'+q(\vec{r},\hat{\Omega},t),
\end{equation}
where the integration over $\mathbb{S}^2$ is often written as an integral over $4\pi$ in the transport literature.

The P$_N$\  approximation expands the angular flux and the anisotropic external source in terms of surface spherical harmonics
\begin{align}
    \psi(\vec{r},\hat{\Omega},t)&=\sum_{l=0}^\infty\frac{2l+1}{4\pi}\psi_l(\vec{r},\hat{\Omega},t),\\
    q(\vec{r},\hat{\Omega},t)&=\sum_{l=0}^\infty\frac{2l+1}{4\pi}q_l(\vec{r},\hat{\Omega},t),
\end{align}
where
\begin{align}
    \psi_l(\vec{r},\hat{\Omega},t)&=\sum_{m=-l}^l\psi_l^m(\vec{r},t)Y_l^m(\hat{\Omega}),\\
    q_l(\vec{r},\hat{\Omega},t)&=\sum_{m=-l}^lq_l^m(\vec{r},t)Y_l^m(\hat{\Omega}),
\end{align}
and $Y_l^m$ are the surface spherical harmonics \cite{Bell_Glasstone}, commonly referred to without the ``surface'' adjective.

The scattering cross-section can be expanded into
\begin{equation}
    \sigma_s(\hat{\Omega}\cdot\hat{\Omega}')=\sum_{l=0}^\infty\frac{2l+1}{4\pi}\sigma_{sl}\sum_{m=-l}^lY_l^m(\hat{\Omega})Y_l^m(\hat{\Omega}'),
\end{equation}
so that, by defining
\begin{equation}
    \sigma_l:=\sigma_t-\sigma_{sl},
\end{equation}
Equation (\ref{TransportEquation}), by orthogonality of surface spherical harmonics, is written as the sum over $l$-moments
\begin{equation}
    \label{ExpandedTransportEquation}
    \sum_{l=0}^\infty(2l+1)\left[\left(\frac{1}{v}\frac{\partial}{\partial t}+\hat{\Omega}\cdot\nabla+\sigma_l\right)\psi_l(\vec{r},\hat{\Omega},t)-q_l(\vec{r},\hat{\Omega},t)\right]=0.
\end{equation}



\subsection{Galerkin Projection}\label{SPnDerivationSection:GalerkinProjection}
The next step is to split the streaming term $(2l+1)\hat{\Omega}\cdot\nabla\psi_l(\vec{r},\hat{\Omega},t)$ into harmonics of degree $l-1$ and $l+1$ respectively, and this is where solid harmonics come to play their role.

Harmonic polynomials, often referred to as solid harmonics, are polynomials that are solutions to the Laplace equation, and a surface spherical harmonic of degree $l$ is the restriction onto $\mathbb{S}^2$ of a homogeneous harmonic polynomial on $\mathbb{R}^3$ of degree $l$ whose collection is denoted as $\mathcal{H}_l(\mathbb{R}^3) $\cite{Axler2013:HarmonicFunctionTheory}. For $\vec{\Omega}\in\mathbb{R}^3$, solid harmonics can be obtained by rescaling surface harmonics by a factor of the length of $\vec{\Omega}$ to the $l$th power, $\Omega^l$, and the process relaxes the restriction of the domain from $\mathbb{S}^2$ to $\mathbb{R}^3$. In other words, if $Y_l^m(\hat{\Omega})$ is a surface harmonic, then $\Omega^l Y_l^m(\vec{\Omega})$ is a solid harmonic. For convenience, when $p_l$ is a surface harmonic, we denote the associated solid harmonic as $\tilde{p}_l$.

Now we make the following definitions for $\vec{\Omega}\in\mathbb{R}^3$
\begin{equation}
    \nabla_\Omega:=\left(\frac{\partial}{\partial\Omega_x},\frac{\partial}{\partial\Omega_y},\frac{\partial}{\partial \Omega_z}\right)^\text{T},
\end{equation}
\begin{equation}
    \mathcal{P}:=\nabla_\Omega\cdot\nabla,\qquad
    \mathcal{Q}:=\vec{\Omega}\cdot\nabla,
\end{equation}
and rewrite the streaming term $(2l+1)\hat{\Omega}\cdot\nabla\psi_l(\vec{r},\hat{\Omega},t)$ as
\begin{equation}
    \label{SplittingOfStreamingTerm}
    (2l+1)\hat{\Omega}\cdot\nabla\psi_l(\vec{r},\hat{\Omega},t)=\left[(2l+1)\mathcal{Q}\tilde{\psi}_l-\Omega^2\mathcal{P}\tilde{\psi}_l+\mathcal{P}\tilde{\psi}_l\right]_{\Omega=1}(\vec{r},\hat{\Omega},t),
\end{equation}
where the subscript $\Omega=1$ denotes that the term inside brackets is restricted to the unit sphere.

It is proved in Reference \cite{Ackroyd1999:IsotropicSHPn} and shown in \ref{Appendix:ProofOfTwoHarmonicPolynomials} that the term $\left[(2l+1)\mathcal{Q}\tilde{\psi}_l-\Omega^2\mathcal{P}\tilde{\psi}_l\right]$ is a solid harmonic of degree $l+1$ and $\mathcal{P}\tilde{\psi}_l$ is a solid harmonic of degree $l-1$. Ergo, when restricted back to $\mathbb{S}^2$, we get surface harmonics of degree $l-1$ and $l+1$.

Therefore, with Galerkin projection applied to Equation (\ref{ExpandedTransportEquation}), we get
\begin{equation}
    \label{PrimitiveProjectedGroupedTransportEquation}
    \left[(2n+1)\left(\frac{1}{v}\frac{\partial}{\partial t}+\sigma_n\right)\tilde{\psi}_n+\left((2n-1)\mathcal{Q}-\Omega^2\mathcal{P}\right)\tilde{\psi}_{n-1}+\mathcal{P}\tilde{\psi}_{n+1}-(2n+1)\tilde{q}_n\right]_{\Omega=1}=0.
\end{equation}
To proceed, we will define two more operators:
\begin{equation}
    \mathcal{F}\psi_n(\vec{r},\hat{\Omega},t):=\left[\left((2n+1)\mathcal{Q}-\Omega^2\mathcal{P}\right)\tilde{\psi}_n\right]_{\Omega=1}(\vec{r},\hat{\Omega},t),
\end{equation}
and
\begin{equation}\label{eq:Gn}
    \mathcal{G}\psi_n(\vec{r},\hat{\Omega},t):=\left[\mathcal{P}\tilde{\psi}_n\right]_{\Omega=1}(\vec{r},\hat{\Omega},t).
\end{equation}
With these definitions Eq.~(\ref{PrimitiveProjectedGroupedTransportEquation}) becomes
\begin{equation}
    \label{GalerkinProjectedGroupedTransportEquation}
    (2n+1)\left(\frac{1}{v}\frac{\partial}{\partial t}+\sigma_n\right)\psi_n(\vec{r},\hat{\Omega},t)+\mathcal{F}\psi_{n-1}(\vec{r},\hat{\Omega},t)+\mathcal{G}\psi_{n+1}(\vec{r},\hat{\Omega},t)=(2n+1)q_n(\vec{r},\hat{\Omega},t).
\end{equation}
It is this form of the moment equations that we will manipulate to derive time-dependent SP$_N$\  equations.


\subsection{Derivation of the SP$_N$\  Equations}\label{SPnDerivationSection:TheLastStep}
To further simplify Equation (\ref{GalerkinProjectedGroupedTransportEquation}), two useful identities from Reference \cite{Ackroyd1999:IsotropicSHPn} would be used
\begin{equation}
    \label{IntermediateIdentity(1)}
    \left(\mathcal{P}\mathcal{Q}-\mathcal{Q}\mathcal{P}-\nabla^2\right)\tilde{\psi}_n=0,
\end{equation}
\begin{equation}
    \label{IntermediateIdentity(2)}
    \mathcal{P}\Omega^2\tilde{\psi}_n=(\Omega^2\mathcal{P}+2\mathcal{Q})\tilde{\psi}_n.
\end{equation}
The two identities can be demonstrated using the definition of the underlying operators. Together, Eqs.~\eqref{IntermediateIdentity(1)} and \eqref{IntermediateIdentity(2)}  imply
\begin{equation}
    \label{DirectIntermediateIdentity}
    \mathcal{P}^{n+1}\left[(2n+1)\mathcal{Q}-\Omega^2\mathcal{P}\right]\tilde{\psi}_n=(n+1)^2\nabla^2\mathcal{P}^n\tilde{\psi}_n,
\end{equation}
as shown in \ref{Appendix:ProofOfDirectIntermediateIdentity}. From these we can use the definition in Eq.~\eqref{eq:Gn} to obtain
\begin{equation}
    \mathcal{G}^n\mathcal{F}\psi_{n-1}=n^2\nabla^2\mathcal{G}^{n-1}\psi_{n-1}.
\end{equation}
Therefore, applying $\mathcal{G}$ $n$ times to Equation (\ref{GalerkinProjectedGroupedTransportEquation}) yields
\begin{equation}
    \label{SPnIntermediateSubstitution(1)}
    (2n+1)\frac{1}{v}\frac{\partial}{\partial t}\mathcal{G}^n\psi_n+(2n+1)\sigma_n\mathcal{G}^n\psi_n+n^2\nabla^2\mathcal{G}^{n-1}\psi_{n-1}+\mathcal{G}^{n+1}\psi_{n+1}=(2n+1)\mathcal{G}^nq_n.
\end{equation}
Notice that $\mathcal{G}^n\psi_n$ and $\mathcal{G}^nq_n$ are surface harmonics of degree 0 and thus have no angular dependence, and that the inverse Laplace operator is well-defined in an infinite-medium problem. Using these facts, we can make the following substitution
\begin{equation}
    F_n(\vec{r})=\frac{1}{n!}\nabla^{-2\ceil{\frac{n}{2}}}\mathcal{G}^n\psi_n,
\end{equation}
\begin{equation}
    G_n(\vec{r})=\frac{1}{n!}\nabla^{-2\ceil{\frac{n}{2}}}\mathcal{G}^nq_n.
\end{equation}
To these definitions we can relate the SP$_N$\  unknowns as, for $n$ a non-negative integer,
\begin{equation}
    \label{SPnMomentDefinition}
        \phi_{2n}=F_{2n},\qquad
        \vec{\phi}_{2n+1}=\nabla F_{2n+1},
\end{equation}
and write the source as
\begin{equation}
        Q_{2n}=G_{2n}, \qquad
        \vec{Q}_{2n+1}=\nabla G_{2n+1}.
\end{equation}
These definitions will
eventually lead us to the SP$_N$\  equations
\begin{equation}
    \label{SPnEquations:EvenEquation}
    \frac{1}{v}\frac{\partial}{\partial t}\phi_n+\sigma_n\phi_n+\frac{n}{2n+1}\nabla\cdot\vec{\phi}_{n-1}+\frac{n+1}{2n+1}\nabla\cdot\vec{\phi}_{n+1}=Q_n\ \ \ (n\text{  even}),
\end{equation}
\begin{equation}
    \label{SPnEquations:OddEquation}
    \frac{1}{v}\frac{\partial}{\partial t}\vec{\phi}_n+\sigma_n\vec{\phi}_n+\frac{n}{2n+1}\nabla\phi_{n-1}+\frac{n+1}{2n+1}\nabla\phi_{n+1}=\vec{Q}_n\ \ \ (n\text{  odd}).
\end{equation}

Therefore, we have started with the time-dependent transport equation in a homogeneous medium and demonstrated that we can derive the SP$_N$\  equations using properties of the solid harmonics. Moreover, due to the way  the derivation proceeded, a truncated P$_N$\  expansion should be equivalent to an SP$_N$\  solution truncated at the same level. Therefore, we expect that in a general problem, with homogeneous media, there to be an equivalence between the SP$_N$\  and P$_N$\  solutions.

\section{Numerical Examples}
In this section, we perform 2D simulations to demonstrate the equivalence between the SP$_N$\  and the P$_N$\  approximation for time-dependent, anisotropic scattering problems.

We solve the time-dependent SP$_N$\  and P$_N$\  equations under the implicit Euler scheme, so the problem is essentially a sequence of steady-state problems; the finite difference method is adopted so that all first-order derivatives are approximated by central differentiation; all moments are nodal values evaluated at grid points.
For all the numerical examples below, periodic boundary conditions are enforced, the external source is isotropic and scattering is anisotropic.



\subsection{An Anisotropic Box Problem}
The first test problem comes from Reference \cite{McClarren2010:SPnReview}, where McClarren used it to demonstrate the equivalence between steady-state P$_N$\  and SP$_N$\ , and Reference \cite{Seibold2014:StaRMAP} refers to it as the ``box" problem. The original problem has $\sigma_t=1\ cm^{-1}$ and isotropic scattering $\sigma_{s0}=0.1\ cm^{-1}$; here we add additional anisotropic scattering $\sigma_{s1}=0.1\ cm^{-1}=\sigma_{s2}=0.1\ cm^{-1}$ to show that such equivalence still holds for time-dependent anisotropic-scattering problems. The external source remains prescribed by
\begin{equation}
    Q(x,y)=
    \begin{cases}
    1   &   1.75\leq x\leq 2.25,1.75\leq y\leq 2.25\\
    1   &   2.75\leq x\leq 3.25,1.5\leq y\leq 2.5\\
    1   &   1.75\leq x\leq 2.25,2.75\leq y\leq 3.25\\
    1   &   3.5\leq x\leq 4.25,3.5\leq y\leq 3.75\\
    0   &   \text{otherwise}
    \end{cases}
\end{equation}

We perform the computation with $\Delta x=0.03\ \text{cm}$ and $\Delta t=0.03\ \text{s}$. The SP$_3$ scalar flux at $T=0.75\ \text{s}$ is compared with that at $T=1.5\ \text{s}$ in Figure\ \ref{AnisotropicBoxProblem:SP3AtT=0.75AndT=1.5}; it can be seen that at $T=0.75\ \text{s}$, the problem is still in a transient state, thus we choose it as the final time of simulation.
\begin{figure}[!htb]
\center
\begin{minipage}{0.49\linewidth}
\center
\includegraphics[width=5cm]{AnisotropicBoxProblemSP3AtT=0d75.png}
\subcaption{$T=0.75\ \text{s}$}
\end{minipage}
\begin{minipage}{0.49\linewidth}
\center
\includegraphics[width=5cm]{AnisotropicBoxProblemSP3AtT=1d5.png}
\subcaption{$T=1.5\ \text{s}$}
\end{minipage}
\caption{Comparison of SP$_3$ scalar flux at $T=0.75\ \text{s}$ and at $T=1.5\ \text{s}$}
\label{AnisotropicBoxProblem:SP3AtT=0.75AndT=1.5}
\end{figure}

The scalar flux solutions for $N=1$ and $N=5$ at $T=0.75\ \text{s}$ are shown in Figure\ \ref{BoxProblem:2dPlotN=1And5}; we also plot the scalar flux along $y=x$ for $N=3$ and $N=7$ in Figure\ \ref{BoxProblem:DiagonalPlotN=3And7}.

\begin{figure}[!htb]
\center
\begin{minipage}{0.49\linewidth}
\center
\includegraphics[width=5cm]{AnisotropicBoxProblemSP1.png}
\subcaption{SP$_1$}
\end{minipage}
\begin{minipage}{0.49\linewidth}
\center
\includegraphics[width=5cm]{AnisotropicBoxProblemP1.png}
\subcaption{P$_1$}
\end{minipage}
\begin{minipage}{0.49\linewidth}
\center
\includegraphics[width=5cm]{AnisotropicBoxProblemSP5.png}
\subcaption{SP$_5$}
\end{minipage}
\begin{minipage}{0.49\linewidth}
\center
\includegraphics[width=5cm]{AnisotropicBoxProblemP5.png}
\subcaption{P$_5$}
\end{minipage}
\caption{Scalar flux computed at $T=0.75\ \text{s}$ for $N=1$ and $N=5$}
\label{BoxProblem:2dPlotN=1And5}
\end{figure}

\begin{figure}[!htb]
\center
\begin{minipage}{0.49\linewidth}
\center
\includegraphics[width=5cm]{AnisotropicBoxProblemDiagonalN=3.pdf}
\subcaption{$N=3$}
\end{minipage}
\begin{minipage}{0.49\linewidth}
\center
\includegraphics[width=5cm]{AnisotropicBoxProblemDiagonalN=7.pdf}
\subcaption{$N=7$}
\end{minipage}
\caption{Scalar flux computed at $T=0.75\ \text{s}$ along $y=x$ for $N=3$ and $N=7$}
\label{BoxProblem:DiagonalPlotN=3And7}
\end{figure}

As can be seen from Figure\ \ref{BoxProblem:2dPlotN=1And5}-\ref{BoxProblem:DiagonalPlotN=3And7}, despite the existence of oscillations, the SP$_N$\  solutions agree with the P$_N$\  solutions. In fact, with $N$ ranging from 1 to 7, the $L^\infty$ norm of the deviation from SP$_N$\  solution to P$_N$\  solution does not exceed $1.6653\times10^{-16}$, so we can claim that the two methods generate identical results.



\subsection{An Anisotropic Line Source Problem}
The line source (pulse) problem is a typical example to demonstrate the Gibbs phenomenon encountered by spectral methods, and it is also important as it is the Green's function for problems with isotropic scattering and source. Here we modify the problem configuration to make the scattering anisotropic such that $\sigma_t=\sigma_{s0}=1\ cm^{-1}$, $\sigma_{s1}=0.5\ cm^{-1}$, $\sigma_{s2}=0.25\ cm^{-1}$ and $\sigma_{s3}=0.125\ cm^{-1}$, but the initial angular flux remains isotropic. To avoid the the singularity in the initial condition, the initial scalar flux $\phi(\vec{r},0)$, as suggested by Reference \cite{Seibold2014:StaRMAP}, is replaced by the distribution:
\begin{equation}
    \label{TestProblem1:InitialScalarFlux}
    \phi(\vec{r},0)=\frac{1}{4\pi s}\exp\left(-\frac{r^2}{4s}\right)
\end{equation}
where $s=3.2\times10^{-4}$.

We perform numerical computation for $T=0.5\ \text{s}$ with $\Delta x=0.01\ \text{cm}$ and $\Delta t=0.005\ \text{s}$. Figure\ \ref{AnisotropicLineSourceProblem:CompareSP3WithIsotropicProblemCutAlongY=0AtT=0.5} shows that the solution of the anisotropic line source problem is different from that of the original line source problem, where the scattering is isotropic. The $L^\infty$ norm of the difference is 1.1810.
\begin{figure}[!htb]
\center
\includegraphics[width=5cm]{IsotropicAndAnisotropicLineSourceProblemSP3Y=0PlotAtT=0d5.pdf}
\caption{Comparison of SP$_3$ scalar flux of the original and the isotropic line source problem on a cut along the positive $x$-axis at $T=0.5\ \text{s}$}
\label{AnisotropicLineSourceProblem:CompareSP3WithIsotropicProblemCutAlongY=0AtT=0.5}
\end{figure}

The scalar flux solutions for $N=3$ and $N=7$ are compared in Figure\ \ref{AnisotropicLineSourceProblem:Order3And7AtT=0.5}, where one can easily tell that the SP$_N$\  and P$_N$\  solutions are almost identical.

\begin{figure}[!htb]
\center
\begin{minipage}{0.49\linewidth}
\center
\includegraphics[width=5cm]{AnisotropicLineSourceProblemSP3AtT=0d5.png}
\subcaption{SP$_3$}
\end{minipage}
\begin{minipage}{0.49\linewidth}
\center
\includegraphics[width=5cm]{AnisotropicLineSourceProblemP3AtT=0d5.png}
\subcaption{P$_3$}
\end{minipage}
\begin{minipage}{0.49\linewidth}
\center
\includegraphics[width=5cm]{AnisotropicLineSourceProblemSP7AtT=0d5.png}
\subcaption{SP$_7$}
\end{minipage}
\begin{minipage}{0.49\linewidth}
\center
\includegraphics[width=5cm]{AnisotropicLineSourceProblemP7AtT=0d5.png}
\subcaption{P$_7$}
\end{minipage}
\caption{Scalar flux for $N=3$ and 7 at $T=0.5\ \text{s}$}
\label{AnisotropicLineSourceProblem:Order3And7AtT=0.5}
\end{figure}

We also plot the scalar flux along the positive $x$-axis in Figure\ \ref{AnisotropicLineSourceProblem:CutAlongY=0AtT=0.5} for odd $N$ from 1 to 7. Though not shown in these figures, the SP$_N$\  solutions also agree with the P$_N$\  solutions for $N$=2, 4 and 6. For $N$ ranging from 1 to 7, the maximum of the $L^\infty$ norms of the SP$_N$\  and P$_N$\  scalar flux difference is $1.9540\times10^{-14}$.

\begin{figure}[!htb]
\center
\includegraphics[width=6cm]{AnisotropicLineSourceProblemY=0PlotAtT=0d5.pdf}
\caption{Scalar flux on a cut along the positive $x$-axis at $T=0.5\ \text{s}$}
\label{AnisotropicLineSourceProblem:CutAlongY=0AtT=0.5}
\end{figure}



\FloatBarrier



\section{Discussions and Conclusions}
With the use of harmonic polynomials, we derived the SP$_N$\  equations starting from the monoenergetic transport equation in homogeneous media. The derivation implies that we can expect the same scalar flux solution from either the SP$_N$\  or P$_N$\  approximation.

The theory is then followed by two numerical examples which are modified from the box problem and the line source problem to have anisotropic scattering. In both examples, the SP$_N$\  and P$_N$\  solutions match very well as the theoretical prediction.

We hope that the work can extend the research and application of the SP$_N$\  method. In fact, the study of the SP$_N$\  equations is far from over because the strongest theoretical results exist for homogenous, infinite media. The boundary and interface conditions the give SP$_N$-P$_N$ equivalence in heterogeneous problems are a major, open problem. A less ambitious, but useful result, would be the development an algorithm that efficiently generates SP$_N$\  moments for anisotropic external source.



\section*{Acknowledgement}
This project is funded, in part, by Department of Energy NEUP research grant from
Battelle Energy Alliance, LLC- Idaho National Laboratory, Contract No: C12-00281.



