Double PN Benchmark Solution for the 1D Monoenergetic Neutron Transport Equation in Plane Geometry

1. Introduction

Boltzmann's equation of particle transport, indeed presents a significant challenge and noteworthy opportunity to solve because of its complexity and wide range of physical phenomena it describes. Originally, the non-linear integro-differential equation, as prescribed by kinetic theory of particle motion, was considered unsolvable. With time however, and advances in mathematics and physical applications, where, in some cases, non-linearity could be relaxed to give a linear equation, the situation changed. In the early to mid twentieth century, a flurry of analytical solutions were constructed for the linear and linearized Boltzmann equation primarily based on solving partial differential equations (PDEs) with distributions admitted, specifically for one-dimension. Alongside the development of analytical solutions were numerical solutions as well such as Monte Carlo, thus enabling the practical use of Boltzmann’s equation in nuclear physics experiments and nuclear reactor physics. With numerical solutions and applications to sensitive physics, there arose a need for numerical method’s verification, which led to the development of benchmarks and benchmarking. This, in turn, led to a host of numerical benchmark solutions to more relevant transport applications requiring sophisticated benchmarking techniques but still generally limited to model problems. Some readers may have the misconception that benchmarking comprehensive transport algorithms is a wasted exercise since only idealized cases can be treated; however, the opposite is true. Even the most advanced numerical method to solve the transport equation can contain unknown errors that a benchmark, even for a simple problem, can easily find. The following is about developing another high precision benchmark in one-dimension.

2. Proposed DPN Algorithm

The PN solution is a well-known solution to determine the angular flux as described by the transport equation. In the method, one expresses the angular flux to the transport equation in a 1D-plane slab as the full-range PN infinite series approximation. However, this approach is problematic for two fundamental reasons. The first concerns how the equations are closed since there is always one more unknown than equation. The second is how one best represents half-range boundary conditions by a continuous full-range expansion, which is not possible in an exact way. There have been several methodologies proposed on how to treat these outstanding issues [See Refs. 1-4]. The simplest approach to specifying boundary conditions is to avoid the difficulty altogether by resorting to the DPN moments approximation, where the expansions are split between forward and backward neutron directions. Now, the moments over the incoming fluxes on the slab near and far boundaries are exact in the limit. Since the DPN approximation also has the advantage of being a coupled set of ODEs for flux moments, one can solve for them following well-established methods. Furthermore, one way to close the system is through convergence acceleration as N approaches infinity.

2.1. DPN Moment and Parity Equations

We consider the simplest case of 1D plane geometry where scattering is isotropic, for which the transport equation is

$$\left[\mu \frac{\partial }{\partial x}+1\right]\varphi \left(x,\mu \right)=\frac{\omega }{2}{\displaystyle \underset{-1}{\overset{1}{\int }}d{\mu }^{\prime }}\varphi \left(x,{\mu }^{\prime }\right)+Q\left(x,\mu \right)$$

(1a)

$$\begin{array}{l}\varphi \left(0,\mu \right)=f\left(\mu \right)\\ \varphi \left(a,-\mu \right)=0.\end{array}$$

(1b)

$\varphi \left(x,\mu \right)$ is the neutron angular flux, x is its position, μ is the neutron direction cosine, ω is the probability of scattering and Q is an imposed external source. The DPN approximation assumes the flux representation

$$\varphi \left(x,\mu \right)=\left\{\begin{array}{c}{\displaystyle \sum _{l=0}^{\infty }\left(2l+1\right){\varphi }_l^{+}\left(x\right)}{P}_l\left(2\mu -1\right), \mu \ge 0\\ {\displaystyle \sum _{l=0}^{\infty }\left(2l+1\right){\varphi }_l^{-}\left(x\right)}{P}_l\left(2\mu +1\right), \mu \le 0,\end{array}\right.$$

(2a,b)

where $P_l\left(2\mu \pm 1\right)$ is the half-range Legendre polynomial of degree l. When projected over half-range polynomials $P_j\left(2\mu \pm 1\right)$ over intervals $\left[0,1\right]$ and $\left[-1,0\right]$ respectively and from orthogonality

$${\displaystyle \underset{0}{\overset{1}{\int }}d\mu }{P}_j\left(2\mu \pm 1\right)P_l\left(2\mu \pm 1\right)=\frac{1}{2l+1}{\delta }_{j,l},$$

(3)

the coefficients in Eqs(2) become the moments over the positive and negative directions,

$${\varphi }_j^{\pm }\left(x\right)={\displaystyle \underset{\left(0,-1\right)}{\overset{\left(1,0\right)}{\int }}d\mu }{P}_j\left(2\mu \mp 1\right)\varphi \left(x,\mu \right).$$

(4)

The integral on the RHS of Eq(1a), called the scalar flux, in terms of moments is

$$\varphi \left(x\right)={\displaystyle \underset{-1}{\overset{1}{\int }}d{\mu }^{\prime }}\varphi \left(x,{\mu }^{\prime }\right)={\displaystyle \underset{0}{\overset{1}{\int }}d{\mu }^{\prime }}\varphi \left(x,{\mu }^{\prime }\right)+{\displaystyle \underset{-1}{\overset{0}{\int }}d{\mu }^{\prime }}\varphi \left(x,{\mu }^{\prime }\right)={\varphi }_0^{+}\left(x\right)+{\varphi }_0^{-}\left(x\right)$$

(5)

and the transport equation to solve becomes

$$\left[\mu \frac{\partial }{\partial x}+1\right]\varphi \left(x,\mu \right)=\frac{\omega }{2}\left[{\varphi }_0^{+}\left(x\right)+{\varphi }_0^{-}\left(x\right)\right]+Q\left(x,\mu \right).$$

(6)

To find a recurrence relation for the moments, we first multiply Eq(6) over $P_j\left(2\mu \mp 1\right)$

$$\begin{array}{l}\left[\frac{\partial }{\partial x}\mu P_j\left(2\mu \mp 1\right)+P_j\left(2\mu \mp 1\right)\right]\varphi \left(x,\mu \right)=\\ =\frac{\omega }{2}\left[{\varphi }_0^{+}\left(x\right)+{\varphi }_0^{-}\left(x\right)\right]P_j\left(2\mu \mp 1\right)+Q\left(x,\mu \right)P_j\left(2\mu \mp 1\right);\end{array}$$

(7a)^±

and from the recurrence of half-range Legendre polynomials, substitute

$$\mu P_j\left(2\mu \mp 1\right)=\frac{1}{2}\left[\begin{array}{l}\frac{j+1}{2j+1}{P}_{j+1}\left(2\mu \mp 1\right)+\\ +\frac{j}{2j+1}{P}_{j-1}\left(2\mu \mp 1\right)\pm P_j\left(2\mu \mp 1\right)\end{array}\right],$$

(7b)^±

to find

$$\begin{array}{l}\frac{1}{2}\frac{\partial }{\partial x}\left[\frac{j+1}{2j+1}{P}_{j+1}\left(2\mu \mp 1\right)+\frac{j}{2j+1}{P}_{j-1}\left(2\mu \mp 1\right)\pm P_j\left(2\mu \mp 1\right)\right]\varphi \left(x,\mu \right)+\\ +P_j\left(2\mu \mp 1\right)\varphi \left(x,\mu \right)=\frac{\omega }{2}\left[{\varphi }_0^{+}\left(x\right)+{\varphi }_0^{-}\left(x\right)\right]P_j\left(2\mu \mp 1\right)+Q\left(x,\mu \right)P_j\left(2\mu \mp 1\right).\end{array}$$

(7c)^±

Finally, on projection over intervals $\left[0,1\right]$ and $\left[-1,0\right]$ respectively applying orthogonality

$$\begin{array}{l}\frac{1}{2}\left[\frac{j+1}{2j+1}\frac{d{\varphi }_{j+1}^{\pm }\left(x\right)}{dx}+\frac{j}{2j+1}\frac{d{\varphi }_{j-1}^{\pm }\left(x\right)}{dx}\pm \frac{d{\varphi }_j^{\pm }\left(x\right)}{dx}\right]+{\varphi }_j^{\pm }\left(x\right)=\\ =\frac{\omega }{2}\left[{\varphi }_0^{+}\left(x\right)+{\varphi }_0^{-}\left(x\right)\right]{\delta }_{j,0}+Q_j^{\pm }\left(x\right).\end{array}$$

(7d)^±

To arrive at the final result, change j to l and multiply by ${\left(\pm 1\right)}^l$ to give for l=0,1,…,N-1

$$\begin{array}{l}\frac{1}{2}\left[\frac{l+1}{2l+1}\frac{d{\tilde{\varphi }}_{l+1}^{\pm }\left(x\right)}{dx}+\frac{l}{2l+1}\frac{d{\tilde{\varphi }}_{l-1}^{\pm }\left(x\right)}{dx}+\frac{d{\tilde{\varphi }}_l^{\pm }\left(x\right)}{dx}\right]\pm {\varphi }_l^{\pm }\left(x\right)=\\ =\pm \frac{\omega }{2}\left[{\tilde{\varphi }}_0^{+}\left(x\right)+{\tilde{\varphi }}_0^{-}\left(x\right)\right]{\delta }_{l,0}\pm {\tilde{Q}}_l^{\pm }\left(x\right).\end{array}$$

(7e)^±

when

$${\tilde{\varphi }}_l^{\pm }\left(x\right)\equiv {\left(\pm 1\right)}^l{\varphi }_l^{\pm }\left(x\right),$$

(7f)^±

and for closure, we assume

$$\frac{d{\tilde{\varphi }}_{N+1}^{\pm }\left(x\right)}{dx}\equiv 0.$$

(7g)^±

By defining the vectors

$$\begin{array}{l}{\tilde{\varphi }}^{\pm }\left(x\right)\equiv {\left[\begin{array}{cccc}{\tilde{\varphi }}_0^{\pm }& {\tilde{\varphi }}_1^{\pm }& \mathrm{...}& {\tilde{\varphi }}_N^{\pm }\end{array}\right]}^T\\ {\tilde{Q}}^{\pm }\left(x\right)\equiv {\left[\begin{array}{cccc}{\tilde{Q}}_0^{\pm }& {\tilde{Q}}_1^{\pm }& \mathrm{...}& {\tilde{Q}}_N^{\pm }\end{array}\right]}^T,\end{array}$$

(8a,b)^±

Eq(7e) becomes

$$\begin{array}{l}\frac{1}{2}A\frac{d{\tilde{\varphi }}^{\pm }\left(x\right)}{dx}\pm {\tilde{\varphi }}^{\pm }\left(x\right)=\pm \frac{\omega }{2}\left[{\tilde{\varphi }}_0^{+}\left(x\right)+{\tilde{\varphi }}_0^{-}\left(x\right)\right]1_0\pm {\tilde{Q}}^{\pm }\left(x\right),\\ \end{array}$$

(9a)^±

where

$$\begin{array}{l}A\equiv \left\{{\delta }_{l-1,j}l/(2l+1), {\delta }_{l,j}, {\delta }_{l+1,j}\left(l+1\right)/(2l+1); l,j=0,\mathrm{..},N-1\right\}\\ 1_0\equiv {\left\{1 0 0....0\right\}}^T.\end{array}$$

(9b,c)

Expressing the RHS more conveniently in terms of $\psi \left(x\right)$ gives

$$\begin{array}{l}\frac{1}{2}A\frac{d{\tilde{\varphi }}^{\pm }\left(x\right)}{dx}\pm {\tilde{\varphi }}^{\pm }\left(x\right)=\pm \frac{\omega }{2}\delta \left[{\tilde{\varphi }}^{+}\left(x\right)+{\tilde{\varphi }}^{-}\left(x\right)\right]\pm {\tilde{Q}}^{\pm }\left(x\right),\\ \end{array}$$

(10a)^±

where the matrix $\delta$ is

$$\delta \equiv 1_0{1}_0{ }^T=\left\{{\delta }_{0,0}\right\}.$$

(10b)

The key to solving for the flux vector are the even-odd parity equations derived next.

2.2. The Even/Odd DPN Parity Equations

We form the even and odd parity moment vectors

$$\begin{array}{l}\psi \left(x\right)\equiv {\tilde{\varphi }}^{+}\left(x\right)+{\tilde{\varphi }}^{-}\left(x\right)\\ \chi \left(x\right)\equiv {\tilde{\varphi }}^{+}\left(x\right)-{\tilde{\varphi }}^{-}\left(x\right)\end{array}$$

(11a,b)

with the flux moments recovered from

$${\tilde{\varphi }}^{\pm }\left(x\right)=\frac{1}{2}\left[\psi \left(x\right)\pm \chi \left(x\right)\right].$$

(11c)^±

Similarly, the source parity moments are

$$q^{\pm }\left(x\right)\equiv Q^{+}\left(x\right)\pm Q^{-}\left(x\right).$$

(11d)^±

The parity equations come from adding and subtracting Eqs(10a)^± to give

$$\begin{array}{l}\frac{1}{2}A\frac{d\psi \left(x\right)}{dx}+\chi \left(x\right)=q^{-}\left(x\right)\\ \frac{1}{2}A\frac{d\chi \left(x\right)}{dx}+\psi \left(x\right)=\omega \delta \psi \left(x\right)+q^{+}\left(x\right)\end{array}$$

(12a,b)

On combining Eqs(12) by differentiating Eq(12a)

$$\frac{1}{2}A\frac{d^2\psi \left(x\right)}{d{x}^2}+\frac{d\chi \left(x\right)}{dx}=\frac{d{q}^{-}\left(x\right)}{dx}.$$

(13)

Then, multiplying by$A/2$

$$\frac{1}{4}{A}^2\frac{d^2\psi \left(x\right)}{d{x}^2}+\frac{1}{2}A\frac{d\chi \left(x\right)}{dx}=\frac{1}{2}A\frac{d{q}^{-}\left(x\right)}{dx}$$

(14a)

and introducing Eq(12b)

$$\frac{1}{4}{A}^2\frac{d^2\psi \left(x\right)}{d{x}^2}-\left[I-\omega \delta \right]\psi \left(x\right)=\frac{1}{2}A\frac{d{q}^{-}\left(x\right)}{dx}-q^{+}\left(x\right)$$

(14b)

and simplifying

$$\frac{d^2\psi \left(x\right)}{d{x}^2}-{\Gamma }^2\psi \left(x\right)=2A^{-1}\frac{d{q}^{-}\left(x\right)}{dx}-4A^{-2}{q}^{+}\left(x\right),$$

(14c)

where

$${\Gamma }^2\equiv 4A^{-2}\left[I-\omega \delta \right].$$

(14d)

Once $\psi \left(x\right)$ is found from Eq(14b), then from Eq(12a)

$$\chi \left(x\right)=-\frac{1}{2}A\frac{d\psi \left(x\right)}{dx}+q^{-}\left(x\right).$$

(15)

The two moments of boundary conditions for the solution of Eq(14b) are found from the incoming flux moments at the boundaries [see Eq(1b)]

$$\begin{array}{l}{\varphi }_l^{+}\left(0\right)={\displaystyle \underset{0}{\overset{1}{\int }}d\mu }{P}_l\left(2\mu -1\right)f\left(\mu \right)\\ {\varphi }_l^{-}\left(a\right)=0.\end{array}$$

(16a,b)

The exiting moments will be given by the DPN solution designed to accommodate the half-range condition; therefore, the BCs for $\psi \left(x\right)$,

$$\begin{array}{l}\psi \left(0\right)\equiv {\tilde{\varphi }}^{+}\left(0\right)+{\tilde{\varphi }}^{-}\left(0\right)\\ \psi \left(a\right)\equiv {\tilde{\varphi }}^{+}\left(a\right),\end{array}$$

(17a,b)

are not fully determined. This is where the response matrix enters the analysis.

It should be stated that the approach taken thus far is not unique. Instead of a second order ODE for $\psi \left(x\right)$, a second order ODE can equally be found for $\chi \left(x\right)$. However, there are additional technical issues with the second approach, and therefore will not be further considered.

3. Solution to the Parity Equations

Our task is to solve the following second order inhomogeneous ODE:

$$\frac{d^2\psi \left(x\right)}{d{x}^2}-{\Gamma }^2\psi \left(x\right)=Q\left(x\right),$$

(18a)

with inhomogeneous term

$$Q\left(x\right)\equiv 2A^{-1}\frac{d{q}^{-}\left(x\right)}{dx}-4A^{-2}{q}^{+}\left(x\right).$$

(18b)

The solution naturally decomposes into additive solutions to the homogeneous and the particular parts.

$$\psi \left(x\right)={\psi }_h\left(x\right)+{\psi }_p\left(x\right).$$

(19)

We begin with the homogeneous parity equation for ${\psi }_h\left(x\right)$,

$$\frac{d^2{\psi }_h\left(x\right)}{d{x}^2}-{\Gamma }^2{\psi }_h\left(x\right)=0,$$

(20a)

and apply diagonalization to $\Gamma$

$$\Gamma =Tdiag\left\{{\lambda }_k;k=1,\mathrm{...},N\right\}{T}^{-1},$$

(20b)

where T is a matrix of size N whose columns are eigenvectors T_k of $\Gamma$ corresponding to the eigenvalues ${\lambda }_k$. With diagonalization, Eq(20a) becomes

$$\frac{d^{2}y\left(x\right)}{d{x}^2}-diag\left\{{\lambda }_k^2\right\}y\left(x\right)=0,$$

(21a)

with

$$y\left(x\right)\equiv T^{-1}\psi \left(x\right).$$

(21b)

Eqs(21) are scalar ODEs along the diagonal with the convenient solution [5]

$$y\left(x\right)\equiv diag\left\{\frac{\sin \mathrm{h}\left({\lambda }_{k}x\right)}{\sin \mathrm{h}\left({\lambda }_{k}a\right)}\right\}y\left(a\right)+diag\left\{\frac{\sin \mathrm{h}\left({\lambda }_k\left(a-x\right)\right)}{\sin \mathrm{h}\left({\lambda }_{k}a\right)}\right\}y\left(0\right)$$

(22)

chosen to directly incorporate the (unknown) boundary conditions. By re-inserting y from Eq(21b), the solution to Eq(20a) is

$${\psi }_h\left(x\right)=H\left(x\right){\psi }_h\left(a\right)+H\left(a-x\right){\psi }_h\left(0\right),$$

(23a)

where the matrix function $H\left(x\right)$ forms the two solutions to the homogeneous part $H\left(x\right), H\left(a-x\right)$

$$H\left(x\right)\equiv Tdiag\left\{\frac{\sin \mathrm{h}\left({\lambda }_{k}x\right)}{\sin \mathrm{h}\left({\lambda }_{k}a\right)}\right\}{T}^{-1}$$

(23b)

$$H\left(a-x\right)\equiv Tdiag\left\{\frac{\sin \mathrm{h}\left({\lambda }_k\left(a-x\right)\right)}{\sin \mathrm{h}\left({\lambda }_{k}a\right)}\right\}{T}^{-1}.$$

(23c)

From Eq(19)

$$\begin{array}{l}{\psi }_h\left(0\right)=\psi \left(0\right)-{\psi }_p\left(0\right)\\ {\psi }_h\left(a\right)=\psi \left(a\right)-{\psi }_p\left(a\right),\end{array}$$

(24a)

we find for Eq(23a)

$$\psi \left(x\right)-{\psi }_p\left(x\right)=H\left(x\right)\left[\psi \left(a\right)-{\psi }_p\left(a\right)\right]+H\left(a-x\right)\left[\psi \left(0\right)-{\psi }_p\left(0\right)\right]$$

(24b)

or

$$\psi \left(x\right)=H\left(x\right)\psi \left(a\right)+H\left(a-x\right)\psi \left(0\right)+{\psi }_p\left(x\right)-{\psi }_p\left(a\right)-{\psi }_p\left(0\right).$$

(24c)

Note that the solution contains the unknown boundary conditions.

With the chosen two solutions of Eqs(23b,c) to the homogeneous part and from variation of parameters, the particular solution is [5]

$${\psi }_p\left(x\right)=W^{-1}\left[H\left(x\right){\displaystyle \underset{x}{\overset{a}{\int }}d{x}^{\prime }}H\left(a-x^{\prime }\right)Q\left(x^{\prime }\right)+H\left(a-x\right){\displaystyle \underset{0}{\overset{x}{\int }}d{x}^{\prime }}H\left(x^{\prime }\right)Q\left(x^{\prime }\right)\right]$$

(25a)

with Wronskian

$$W^{-1}\equiv -T\left[diag\left\{\frac{sinh\left({\lambda }_{k}a\right)}{{\lambda }_k}\right\}\right]T^{-1}.$$

(25b)

Note that the particular solution has been constructed such that

$${\psi }_p\left(0\right)={\psi }_p\left(a\right)=0.$$

(26)

From Eq(24c)

$$\psi \left(x\right)=H\left(x\right)\psi \left(a\right)+H\left(a-x\right)\psi \left(0\right)+{\psi }_p\left(x\right);$$

(27a)

and from Eq(15)

$$\chi \left(x\right)=-\frac{1}{2}A\left[H^{\prime }\left(x\right)\psi \left(a\right)+H^{\prime }\left(a-x\right)\psi \left(0\right)\right]+U\left(x\right),$$

(27b)

where

$$\begin{array}{l}{H}^{\prime }\left(x\right)\equiv Tdiag\left\{{\lambda }_k\frac{\cos \mathrm{h}\left({\lambda }_{k}x\right)}{\sin \mathrm{h}\left({\lambda }_{k}a\right)}\right\}{T}^{-1}\\ H^{\prime }\left(a-x\right)\equiv -Tdiag\left\{{\lambda }_k\frac{\cos \mathrm{h}\left({\lambda }_k\left(a-x\right)\right)}{\sin \mathrm{h}\left({\lambda }_{k}a\right)}\right\}{T}^{-1}\\ U\left(x\right)\equiv q^{-}\left(x\right)-\frac{1}{2}A{{\psi }^{\prime }}_p\left(x\right).\end{array}$$

(27c,d,e)

The flux moments in Eqs(2a,b) will be constructed from Eqs(27); but since the parity vectors $\psi \left(0\right),\psi \left(a\right)$

$$\begin{array}{l}\psi \left(0\right)={\tilde{\varphi }}^{+}\left(0\right)+{\tilde{\varphi }}^{-}\left(0\right)\\ \psi \left(a\right)={\tilde{\varphi }}^{+}\left(a\right)+{\tilde{\varphi }}^{-}\left(a\right)\end{array}$$

(28)

contain the unknown exiting flux vectors ${\tilde{\varphi }}^{-}\left(0\right),{\tilde{\varphi }}^{+}\left(a\right)$, these must be determined first.

4. Response Matrix

To begin, we define $\widehat{A},\widehat{B}$

$$\begin{array}{l}{\left.H^{\prime }\left(x\right)\right|}_0\equiv Tdiag\left\{\frac{{\lambda }_k}{\sin \mathrm{h}\left({\lambda }_{k}a\right)}\right\}{T}^{-1}\equiv \widehat{B}\\ {\left.H^{\prime }\left(x\right)\right|}_a\equiv Tdiag\left\{{\lambda }_k\frac{\cos \mathrm{h}\left({\lambda }_{k}a\right)}{\sin \mathrm{h}\left({\lambda }_{k}a\right)}\right\}{T}^{-1}\equiv -\widehat{A};\end{array}$$

(29a,b)

and note the following:

$$\begin{array}{l}{\left.H^{\prime }\left(a-x\right)\right|}_0\equiv -Tdiag\left\{{\lambda }_k\frac{\cos \mathrm{h}\left({\lambda }_k\left(a\right)\right)}{\sin \mathrm{h}\left({\lambda }_{k}a\right)}\right\}{T}^{-1}=\widehat{A}\\ {\left.H^{\prime }\left(a-x\right)\right|}_a\equiv -Tdiag\left\{\frac{{\lambda }_k}{\sin \mathrm{h}\left({\lambda }_{k}a\right)}\right\}{T}^{-1}=-\widehat{B}.\end{array}$$

(29c,d)

To recover the outgoing positive and negative moments at the boundaries, one introduces x = 0 and a into Eq(27b) to give

$$\begin{array}{l}\chi \left(0\right)=-\frac{1}{2}A\left[{\left.H^{\prime }\left(x\right)\right|}_0\psi \left(a\right)+{\left.H^{\prime }\left(a-x\right)\right|}_0\psi \left(0\right)\right]+U\left(0\right)\\ \chi \left(a\right)=-\frac{1}{2}A\left[{\left.H^{\prime }\left(x\right)\right|}_a\psi \left(a\right)+{\left.H^{\prime }\left(a-x\right)\right|}_a\psi \left(0\right)\right]+U\left(a\right)\end{array}$$

(30a,b)

or from Eqs(28a,b)

$$\begin{array}{l}\chi \left(0\right)=-\frac{1}{2}A\left[\widehat{B}\psi \left(a\right)+\widehat{A}\psi \left(0\right)\right]+U\left(0\right)\\ \chi \left(a\right)=\frac{1}{2}A\left[\widehat{A}\psi \left(a\right)+\widehat{B}\psi \left(0\right)\right]+U\left(a\right)\end{array}$$

(31a,b)

Substituting Eqs(11a,b) at x = 0 and a respectively yields

$$\begin{array}{l}{\tilde{\varphi }}^{+}\left(0\right)-{\tilde{\varphi }}^{-}\left(0\right)=-\frac{1}{2}A\left\{\widehat{B}\left[{\tilde{\varphi }}^{+}\left(a\right)+{\tilde{\varphi }}^{-}\left(a\right)\right]+\widehat{A}\left[{\tilde{\varphi }}^{+}\left(0\right)+{\tilde{\varphi }}^{-}\left(0\right)\right]+U\left(0\right)\right\}\\ {\tilde{\varphi }}^{+}\left(a\right)-{\tilde{\varphi }}^{-}\left(a\right)=\frac{1}{2}A\left\{\widehat{A}\left[{\tilde{\varphi }}^{+}\left(a\right)+{\tilde{\varphi }}^{-}\left(a\right)\right]+\widehat{B}\left[{\tilde{\varphi }}^{+}\left(0\right)+{\tilde{\varphi }}^{-}\left(0\right)\right]+U\left(a\right)\right\}\end{array}$$

(32a,b)

and if

$$\beta \equiv \frac{A\widehat{B}}{2}\gamma \equiv \frac{A\widehat{A}}{2},$$

(32c,d)

then Eqs(32a,b) become

$$\begin{array}{l}{\tilde{\varphi }}^{+}\left(0\right)-{\tilde{\varphi }}^{-}\left(0\right)=-\beta \left[{\tilde{\varphi }}^{+}\left(a\right)+{\tilde{\varphi }}^{-}\left(a\right)\right]-\gamma \left[{\tilde{\varphi }}^{+}\left(0\right)+{\tilde{\varphi }}^{-}\left(0\right)\right]+U\left(0\right)\\ {\tilde{\varphi }}^{+}\left(a\right)-{\tilde{\varphi }}^{-}\left(a\right)=\gamma \left[{\tilde{\varphi }}^{+}\left(a\right)+{\tilde{\varphi }}^{-}\left(a\right)\right]+\beta \left[{\tilde{\varphi }}^{+}\left(0\right)+{\tilde{\varphi }}^{-}\left(0\right)\right]+U\left(a\right);\end{array}$$

(32e,f)

and re-arranging

$$\begin{array}{l}\beta {\tilde{\varphi }}^{+}\left(a\right)-\left[{\rm I}-\gamma \right]{\tilde{\varphi }}^{-}\left(0\right)=-\beta {\tilde{\varphi }}^{-}\left(a\right)-\left[I+\gamma \right]{\tilde{\varphi }}^{+}\left(0\right)+U\left(0\right)\\ \left[{\rm I}-\gamma \right]{\tilde{\varphi }}^{+}\left(a\right)-\beta {\tilde{\varphi }}^{-}\left(0\right)=\left[I+\gamma \right]{\tilde{\varphi }}^{-}\left(a\right)+\beta {\tilde{\varphi }}^{+}\left(0\right)+U\left(a\right)\end{array}$$

(32g,h)

into matrix form

$$\left[\begin{array}{cc}\beta & -x^{-}\\ x^{-}& -\beta \end{array}\right]\left[\begin{array}{c}{\tilde{\varphi }}^{+}\left(a\right)\\ {\tilde{\varphi }}^{-}\left(0\right)\end{array}\right]=\left[\begin{array}{cc}-\beta & -x^{+}\\ x^{+}& \beta \end{array}\right]\left[\begin{array}{c}{\tilde{\varphi }}^{-}\left(a\right)\\ {\tilde{\varphi }}^{+}\left(0\right)\end{array}\right]+\left[\begin{array}{c}U\left(0\right)\\ U\left(a\right)\end{array}\right],$$

(32i,j)

where

$$x^{\pm }\equiv \left[{\rm I}\pm \gamma \right].$$

(32k)

When we multiply both sides by the skew symmetric block identity matrix

$$\left[\begin{array}{cc}-I& 0\\ 0& I\end{array}\right]\left[\begin{array}{cc}\beta & -x^{-}\\ x^{-}& -\beta \end{array}\right]\left[\begin{array}{c}{\tilde{\varphi }}^{+}\left(a\right)\\ {\tilde{\varphi }}^{-}\left(0\right)\end{array}\right]=\left[\begin{array}{cc}-I& 0\\ 0& I\end{array}\right]\left[\begin{array}{cc}-\beta & -x^{+}\\ x^{+}& \beta \end{array}\right]\left[\begin{array}{c}{\tilde{\varphi }}^{-}\left(a\right)\\ {\tilde{\varphi }}^{+}\left(0\right)\end{array}\right]+\left[\begin{array}{cc}-I& 0\\ 0& I\end{array}\right]\left[\begin{array}{c}U\left(0\right)\\ U\left(a\right)\end{array}\right],$$

(33a)

we find

$$\left[\begin{array}{cc}-\beta & x^{-}\\ x^{-}& -\beta \end{array}\right]\left[\begin{array}{c}{\tilde{\varphi }}^{+}\left(a\right)\\ {\tilde{\varphi }}^{-}\left(0\right)\end{array}\right]=\left[\begin{array}{cc}\beta & x^{+}\\ x^{+}& \beta \end{array}\right]\left[\begin{array}{c}{\tilde{\varphi }}^{-}\left(a\right)\\ {\tilde{\varphi }}^{+}\left(0\right)\end{array}\right]+\left[\begin{array}{c}-U\left(0\right)\\ U\left(a\right)\end{array}\right].$$

(33b)

Finally, multiplying by the inverse of the leading matrix gives the exiting flux

$$\left[\begin{array}{c}{\tilde{\varphi }}^{+}\left(a\right)\\ {\tilde{\varphi }}^{-}\left(0\right)\end{array}\right]=R\left[\begin{array}{c}{\tilde{\varphi }}^{-}\left(a\right)\\ {\tilde{\varphi }}^{+}\left(0\right)\end{array}\right]+{\left[\begin{array}{cc}-\beta & x^{-}\\ x^{-}& -\beta \end{array}\right]}^{-1}\left[\begin{array}{c}-U\left(0\right)\\ U\left(a\right)\end{array}\right],$$

(34a)

where the response matrix is

$$R\equiv {\left[\begin{array}{cc}-\beta & x^{-}\\ x^{-}& -\beta \end{array}\right]}^{-1}\left[\begin{array}{cc}\beta & x^{+}\\ x^{+}& \beta \end{array}\right].$$

(34b)

To complete the expression

$$\left[\begin{array}{c}-U\left(0\right)\\ U\left(a\right)\end{array}\right]=\left[\begin{array}{c}-q^{-}\left(0\right)\\ q^{-}\left(a\right)\end{array}\right]+\frac{1}{2}\left[\begin{array}{c}{{\psi }^{\prime }}_p\left(0\right)\\ -{{\psi }^{\prime }}_p\left(a\right)\end{array}\right],$$

(35a,b)

with

$$\left[\begin{array}{c}{{\psi }^{\prime }}_p\left(0\right)\\ -{{\psi }^{\prime }}_p\left(a\right)\end{array}\right]=W^{-1}\left[\begin{array}{c}\widehat{B}{\displaystyle \underset{0}{\overset{a}{\int }}d{x}^{\prime }}H\left(a-x^{\prime }\right)Q\left(x^{\prime }\right)\\ -\widehat{A}{\displaystyle \underset{0}{\overset{a}{\int }}d{x}^{\prime }}H\left(x^{\prime }\right)Q\left(x^{\prime }\right)\end{array}\right].$$

(35c,d)

On the left, we have the moments of the outgoing flux distribution and on the right the moments of the incoming flux distribution. Hence, for a known incoming distribution, the outgoing moments are now known.

5. Final Moments Solution

One then recovers the N spatial moments from Eqs(11c)^± as

$$\begin{array}{l}{\tilde{\varphi }}^{+}\left(x\right)=\frac{1}{2}\left[\psi \left(x\right)+\chi \left(x\right)\right]\\ {\tilde{\varphi }}^{-}\left(x\right)=\frac{1}{2}\left[\psi \left(x\right)-\chi \left(x\right)\right].\end{array}$$

(36a)^±

and from Eqs(27a,b) to give the general solution

$$\begin{array}{l}{\tilde{\varphi }}^{+}\left(x\right)=\frac{1}{2}\left[\begin{array}{l}\left[H\left(x\right)-\frac{1}{2}A{H}^{\prime }\left(x\right)\right]\psi \left(a\right)+\left[H\left(a-x\right)-\frac{1}{2}A{H}^{\prime }\left(a-x\right)\right]\psi \left(0\right)+\\ +{\psi }_p\left(x\right)+U\left(x\right)\end{array}\right]\\ {\tilde{\varphi }}^{-}\left(x\right)=\frac{1}{2}\left[\begin{array}{l}\left[H\left(x\right)+\frac{1}{2}A{H}^{\prime }\left(x\right)\right]\psi \left(a\right)+\left[H\left(a-x\right)+\frac{1}{2}A{H}^{\prime }\left(a-x\right)\right]\psi \left(0\right)-\\ +{\psi }_p\left(x\right)-U\left(x\right)\end{array}\right].\end{array}$$

(36b)^±

Or with

$$H^{\pm }\left(x\right)\equiv H\left(x\right)\pm \frac{1}{2}A{H}^{\prime }\left(x\right)$$

(37a)^±

and

$$P^{\pm }\left(x\right)\equiv {\psi }_p\left(x\right)\pm U\left(x\right),$$

(37b)^±

we have

$$\begin{array}{l}{\tilde{\varphi }}^{+}\left(x\right)=\frac{1}{2}\left[H^{-}\left(x\right)\psi \left(a\right)+H^{-}\left(a-x\right)\psi \left(0\right)+P^{+}\left(x\right)\right]\\ {\tilde{\varphi }}^{-}\left(x\right)=\frac{1}{2}\left[H^{+}\left(x\right)\psi \left(a\right)-H^{+}\left(a-x\right)\psi \left(0\right)\right]+P^{-}\left(x\right).\end{array}$$

(37c)^±

With additional manipulation

$$\left[\begin{array}{c}{\tilde{\varphi }}^{+}\left(x\right)\\ {\tilde{\varphi }}^{-}\left(x\right)\end{array}\right]=\frac{1}{2}\left[\begin{array}{cc}{H}^{-}\left(x\right)& H^{-}\left(a-x\right)\\ H^{+}\left(x\right)& H^{+}\left(a-x\right)\end{array}\right]\left\{\left[\begin{array}{c}{\tilde{\varphi }}^{+}\left(a\right)\\ {\tilde{\varphi }}^{-}\left(0\right)\end{array}\right]+\left[\begin{array}{c}{\tilde{\varphi }}^{-}\left(a\right)\\ {\tilde{\varphi }}^{+}\left(0\right)\end{array}\right]\right\}+\frac{1}{2}\left\{{\psi }_p\left(x\right)\left[\begin{array}{c}1\\ 1\end{array}\right]+U\left(x\right)\left[\begin{array}{c}1\\ -1\end{array}\right]\right\},$$

(38a)^±

or introducing Eq(34a)

$$\begin{array}{l}\left[\begin{array}{c}{\tilde{\varphi }}^{+}\left(x\right)\\ {\tilde{\varphi }}^{-}\left(x\right)\end{array}\right]=\frac{1}{2}\left[\begin{array}{cc}{H}^{-}\left(x\right)& H^{-}\left(a-x\right)\\ H^{+}\left(x\right)& H^{+}\left(a-x\right)\end{array}\right]\left\{\left[R+I\right]\left[\begin{array}{c}{\tilde{\varphi }}^{-}\left(a\right)\\ {\tilde{\varphi }}^{+}\left(0\right)\end{array}\right]+{\left[\begin{array}{cc}-\beta & x^{-}\\ x^{-}& -\beta \end{array}\right]}^{-1}\left[\begin{array}{c}-U\left(0\right)\\ U\left(a\right)\end{array}\right]\right\}+\\ +\frac{1}{2}\left\{{\psi }_p\left(x\right)\left[\begin{array}{c}1\\ 1\end{array}\right]+U\left(x\right)\left[\begin{array}{c}1\\ -1\end{array}\right]\right\}\end{array}$$

(38b)^±

with

$$\left[\begin{array}{c}{\tilde{\varphi }}^{-}\left(a\right)\\ {\tilde{\varphi }}^{+}\left(0\right)\end{array}\right]=\left[\begin{array}{c}0\\ {\displaystyle \underset{0}{\overset{1}{\int }}d\mu P_N\left(2\mu -1\right)f\left(\mu \right)}\end{array}\right]$$

(38c)^±

and

$$P_N\left(2\mu -1\right)\equiv {\left[\begin{array}{cccc}{P}_0\left(2\mu -1\right)& P_1\left(2\mu -1\right)& \mathrm{...}& P_N\left(2\mu -1\right)\end{array}\right]}^T$$

(38d,e)^±

$$1\equiv {\left\{1 1 ....1\right\}}^T.$$

and $U\left(x\right)$ is from Eq(27e).

6. The DPN Approximation

The DPN solution comes from the convergence of the series solution of Eqs(2)

$$\varphi \left(x,\mu \right)=\underset{N\to \infty }{\lim }\varphi \left(x,\mu ;N\right),$$

(39a)

where the DPN approximation is the partial sum

$$\varphi \left(x,\mu ;N\right)=\left\{\begin{array}{c}{\displaystyle \sum _{l=0}^N\left(2l+1\right){\varphi }_l^{+}\left(x\right)}{P}_l\left(2\mu -1\right), \mu \ge 0\\ {\displaystyle \sum _{l=0}^N\left(2l+1\right){\varphi }_l^{-}\left(x\right)}{P}_l\left(2\mu +1\right), \mu \le 0.\end{array}\right.$$

(39b,c)

To incorporate the moments vector of Eqs(38b), we let μ negative be −|μ| and trivially multiply by (+1)^l in the partial sum for μ positive to give for $\mu \ge 0$

$$\varphi \left(x,\pm \mu ;N\right)={\displaystyle \sum _{l=0}^N\left(2l+1\right){\tilde{\varphi }}_l^{\pm }\left(x\right)}{P}_l\left(2\left|\mu \right|-1\right),$$

(40a)

and if

$$L_N\equiv diag\left\{2l+1, l=1,\mathrm{...},N\right\},$$

(40b)

Eq(40a) becomes the N-term inner product

$$\varphi \left(x,\pm \mu ;N\right)=P_N^T\left(2\left|\mu \right|-1\right)L_N{\tilde{\varphi }}^{\pm }\left(x\right)$$

(40c)

to be evaluated in the next section.