A Double Legendre Polynomial Order N Benchmark Solution for the 1D Monoenergetic Neutron Transport Equation in Plane Geometry

: As more and more numerical and analytical solutions to the linear neutron transport equation become available, veriﬁcation of the numerical results becomes increasingly important. This presentation concerns the development of another benchmark for the linear neutron transport equation in a benchmark series, each employing a diﬀerent method of solution. In 1D, there are numerous ways of analytically solving the monoenergetic transport equation, such as the Wiener–Hopf method, based on the analyticity of the solution, the method of singular eigenfunctions, inversion of the Laplace and Fourier transform solution, and analytical discrete ordinates in the limit, which is arguably one of the most straightforward, to name a few. Another potential method is the PN (Legendre polynomial order N ) method, where one expands the solution in terms of full-range orthogonal Legendre polynomials, and with orthogonality and series truncation, the moments form an open set of ﬁrst -order ODEs. Because of the half-range boundary conditions for incoming particles, however, full-range Legendre expansions are inaccurate near material discontinuities. For this reason, a double PN (DPN) expansion in half-range Legendre polynomials is more appropriate, where one separately expands incoming and exiting ﬂux distributions to preserve the discontinuity at material interfaces. Here, we propose and demonstrate a new method of solution for the DPN equations for an isotropically scattering medium. In comparison to a well-established fully analytical response matrix/discrete ordinate solution (RM/DOM) benchmark using an entirely diﬀerent method of solution for a non-absorbing 1 mfp thick slab with both isotropic and beam sources, the DPN algorithm achieves nearly 8-and 7-place precision, respectively


Introduction
Boltzmann's equation of particle transport indeed presents a significant challenge and noteworthy opportunities to solve because of its complexity and the wide range of physical phenomena it describes [1,2].Originally, the non-linear integro-differential equation, as prescribed by the kinetic theory of particle motion, seemed 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 has changed [3].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 the distributions admitted, specifically in one dimension.Alongside the development of analytical solutions were numerical solutions as well, such as the Monte Carlo and discrete ordinate methods, which, coupled with increasing computational performance, enabled the practical use of Boltzmann's equation in nuclear reactor design.With numerical solutions and applications to particle physics, there arose a need for numerical methods verification, which required 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 these were 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 hold; however, the opposite is true.Even the most advanced numerical method for solving the transport equation can contain undetected errors that a benchmark, even for a simple problem, can easily find.The following presentation is concerned with developing a high precision (DPN) benchmark in one dimension as another of the author's high precision benchmarks.Appendix A summarizes three previous benchmarks to provide a backdrop upon which to compare the DPN benchmark.
The DPN method is equivalent to the response matrix/discrete ordinate method (RM/DOM), with double Gauss quadrature [4,5] for the exiting intensities in directions [−1,0], [0,1] [see explanation of Equation (1)].The RM/DOM will also serve as the standard of comparison to verify the number of places of significance for the DPN benchmark.Rather than discretizing, we proceed analytically with an expansion of the angular flux in half-range Legendre polynomials in order to separate the forward and backward flow streams, which is entirely different from the benchmarks discussed in Appendix A. We introduce the expansions into the transport equation, which, when projected over the halfrange polynomials, yields two coupled first-order vector equations for the moments of expansion in the positive and negative directions.With re-arrangement, the equations take on a parity form, decoupling one equation into a second-order vector ODE, which we solve in the same way as for the RM/DOM but now for the moments not fluxes.The angular flux results in a truncated series expansion from the linear combination of the parities.The numerical implementation then follows the RM/DOM benchmark.The description of the DPN method follows.

The Proposed DPN Algorithm
The PN method, also called the spherical harmonics solution, is a well-known solution for determining angular flux, as described by the transport equation [6,7].Wick [8] and Marshak [9] first applied the method to the neutron transport equation.In addition, Mark [10] developed further details to solve 1D-plane slab problems, which is our focus.In the PN method, one expresses the angular flux in the monoenergetic transport equation for a 1D-plane slab as a full-range infinite series.However, this approach is problematic for two fundamental reasons.The first concern is how the equations close since there is always one more unknown than equation.The second is how one best represents discontinuous half-range boundary conditions by continuous full-range expansion, which is not possible.There have been several methodologies proposed on how to treat these outstanding issues (see Refs. [6][7][8][9][10][11][12]). The simplest approach to specifying boundary conditions is to avoid the difficulty altogether by resorting to the D(ouble)PN moment approximation, where the expansions are split between the forward and backward neutron streams.Now, the expansion over the moments of incoming and exiting fluxes on the slab's near (x = 0) and far (x = a) boundaries are separate exact infinite series.Since the DPN approximation also has the advantage of approximating with a coupled set of ODEs for the flux moments, one can solve for them through well-established methods, of which the method to follow is not.Furthermore, one way to close the system is through convergence acceleration by recalculating the solution as N approaches infinity and noting convergence.
For further clarity, the theoretical and numerical steps for the determination of the DPN benchmark are as follows: 1. Express the neutron angular flux of the transport equation as a series of orthogonal half-range Legendre polynomials and spatial moments for neutrons moving in positive and negative direction cosines (±µ) [Equation (2a,b)].

2.
Since the scattering integral on the RHS of Equation (1a) is expressed as a full-range moment, it is re-expressed as a sum of equivalent half-range moments [Equation ( 5)]. 3. The transport equation is then projected over half-range Legendre polynomials in the forward and backward directions to define two coupled first-order ODEs for the spatial moments in each direction [Equation (7e) ± ]. 4.After the application of the closure condition to the last moment (with the derivative set to zero), the ODEs are put in vector form, added, and subtracted to form the parity equations [Equation (12a,b)].5.By differentiation, the even parity equation becomes a single inhomogeneous secondorder ODE, and the odd parity equation remains first-order [Equation ( 18)].6.The even parity equation solution is expressed as the sum of the solution to the homogenous and particular parts.Note that the solution to the homogeneous part is constructed from matrix functions with assumed (known) boundary conditions [Equation (24c)].7.With the homogeneous solution known, the particular solution comes from the variation of parameters [Equation (25a)].8.The exiting spatial flux moments are recovered by deriving the response matrix connecting the input moments to the output moments across the slab [Equation ( 34)].9.The slab's interior spatial moments then come from adding the parity components (36a) ± .10.With known spatial moments, the last step is numerical evaluation of the half-range Legendre series using the Clenshaw algorithm and Wynn-epsilon convergence acceleration.

DPN Moments and Parity Equations
We consider the simplest case of 1D plane geometry, where scattering is isotropic, for which the neutron transport equation representing the neutron balance equation in the phase space dxdµ for a slab of thickness a, is as follows: ( ) ( ) ( ) , x φ µ is the neutron angular flux, x is its position in units of the mean free path (mfp), µ is the neutron direction cosine in directions [−1,1], ω is the probability of scattering and therefore not being absorbed, Q is an imposed external source, and f is the flux entering the surface at x = 0, with none entering at x = a.The first term represents neutron streaming between collisions, with the second term accounting for loss due to absorption and scattering into a different direction interval and the third term accounting for the gain from scattering into the µ-direction interval.
The DPN solution assumes the flux representations [13] ( ) l P µ ± is the half-range Legendre polynomial of degree l.
( ) are the DPN lth spatial moments by directional projection.When projected over half-range polynomials ( ) the coefficients in Equation (2) become the moments over the positive and negative directions, ( ) where (4) ± indicates + for the top equation and − for the bottom equation.
The integral on the RHS of Equation (1a), called the scalar flux, in terms of moments, is and the transport equation to solve becomes To find a recurrence relation for the moments, we first multiply Equation ( 6) by ( ) from the recurrence of the half-range Legendre polynomials [14], we substitute Finally, for the projection over intervals [ ] 0,1 and [ ] 1, 0 − , respectively, apply the orthogonality of Legendre polynomials [14] ( ) To arrive at the intermediate result, change j to l and multiply by ( ) and for closure at l = N-1, By defining the vectors of length N + 1, ( ) ( ) Expressing the RHS more conveniently in terms of the δ matrix gives where the matrix δ is The key to solving for the flux vector is the even/odd parity equations derived next

The Even/Odd DPN Parity Equations
We form the even and odd parity moment vectors with the flux moments recovered from ( ) ( ) ( ) Similarly, the source parity moments are as follows: The parity equations come from adding and subtracting Equation (10a) ± to give Combine Equation ( 12) by differentiating Equation (12a) as follows: Then, multiply by /2 A , ( ) ( ) ( ) and simplify, We find the two moments of the boundary conditions to apply to the solution of Equation (14b) from the incoming flux moments at the boundaries [see The DPN parity solution contains both known and unknown moments to accommodate the half-range condition; therefore, the BCs are unknown: This is where the response matrix enters the analysis.
It should be stated that the approach taken thus far is not unique.One can solve the first-order ODE, Equation (10a), in terms of the matrix exponential eigenfunctions [15][16][17].
Here, our approach is very different, where we convert the first-order equation into a set of second-order ODEs for ( ) , Equation (14c), and solve in terms of the normalized eigenfunctions.Instead, however, one could equally derive a set of second-order ODEs for ( ) , yielding the same solution in a different form.

Solution to the Parity Equations
Our task is to solve the following second-order inhomogeneous ODE, with the inhomogeneous term The solution naturally decomposes into additive solutions of the homogeneous and particular parts We begin the solution with the homogeneous parity equation, ( ) ( ) and apply matrix diagonalization to Γ [18], where T is a matrix of size N, whose columns are the independent eigenvectors Tk in Γ , corresponding to the eigenvalues k λ .With diagonalization, Equation (20a) becomes Equation ( 21a) is a scalar ODE along the diagonal with the convenient solution [18] ( ) y y y (22) chosen to directly incorporate the (unknown) boundary conditions.By re-inserting y from Equation (21b), the solution to Equation (20a) is where the matrix function,

( )
x H , forms the two independent solutions to the homogeneous part These are called customized eigenfunctions .
From Equation ( 19), we find for Equation (23a) Note that the solution contains the unknown boundary conditions in ( ) ( ) . With the two chosen complimentary solutions of Equation (23b,c) to the homogeneous part and from the variation of parameters, the particular solution is [18] ( ) with Wronskian ( ) Note that the particular solution has been constructed such that ( ) ( ) herefore, from Equation (24c), and from Equation ( 15), we have where The flux moments in the orthogonal expansion of Equation (2a,b) come from Equations (27), but since the parity vectors, ( ) ( ) contain the unknown exiting flux vectors ( ) ( ) , we must determine these first.

The Response Matrix
We define Â and B as and note the following: To recover the outgoing positive and negative moments at the boundaries, one introduces x = 0 and a into Equation (27b) to give or, from Equation (29a-d), Substituting Equation (11a,b) at x = 0 and a, respectively, yields equally by re-arranging, into matrix form where [ ].
x ± ≡ ± Ι γ When multiplied by the skew-symmetric block identity matrix, we find for Equation (32i) Finally, multiplying by the inverse of the leading matrix gives the exiting flux in terms of the response matrix, similarly found for the discrete ordinate method [4], where the response matrix is To complete the expression, we calculate On the left, in Equation (34a), 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.

The Final Moment Solution
One then recovers the N spatial moments from Equation (11c) ± as (where for (36a) ± + for top equation and − for bottom) and from Equation (27a,b) to give the general solution Or, with ( ) ( ) ( ) or introducing Equation (34a), x U is from Equation (27e).

The DPN Approximation
The DPN solution comes from the convergence of the series solution in Equation ( 2), ( ) ( ) where the DPN approximation is the N+1 term partial sums ( ) To incorporate the moment vector from Equation (38b), we let the negative equivalent of µ be −|µ| and trivially multiply by (+1) l in the partial sum for the positive value of µ to give for then Equation (40a) becomes the inner product We will evaluate the above in the next section.

Numerical Implementation for an Isotropic Source
Numerical implementation of the DPN algorithm is relatively straightforward, following the algorithm presented.In particular, matrix diagonalization is performed through an HQR procedure originating from the LAPACK routine [19] and is one of the most reliable numerical methods in use today.Once we know the eigenvalues, the matrix functions come from common matrix multiplication found in Equation (23b,c).LU decomposition [18] gives the matrix inversions.Our final consideration is convergence in DPN through a sequence of DPN approximations [Equation (40c)] and convergence acceleration.To develop a sequence of partial sums, we increment the solution for the number of moments by stride ∆N (usually ∆N = 5) and monitor the convergence of the sequence through sequential convergence and convergence acceleration until the relative error falls below a desired relative error or until convergence fails.If there is a failure, we restart with an increase in the number of sequence elements or change our desired relative error.In this way, convergence provides a measure of closure.There are numerous variants of acceleration to choose from.Here, we first apply common sequential convergence, similar to a sensitivity study, and then Wynn-epsilon (W-e) acceleration [20].Sequential convergence starts from double Legendre polynomial order N0 and compares the convergence standard, also called the engineering estimate, of the relative error between iterations m − 1 and m [i.e., orders ( ) to the desired relative error ε, ( ) If satisfied, the sequence has converged sequentially at (x,µ).W-e acceleration, on the other hand, is non-linear.Similarly, for sequential convergence, the same sequence of DPN approximations is the initial sequence list.W-e convergence essentially extrapolates a known sequence to give the next sequence element, including an estimate of the limit from the previous elements.The algorithm is recursive and for L + 1 iterations, in principle, improves the DPN partial sums result for k even The algorithm, conveniently, is the following tableau: (43) Every element of the even columns should give an improved estimate of the original DPN partial sum from column one.The last term in each even column at the end of the tableau along the diagonal, indicated by the arrow, is assumed to be the most precise and gives the relative error standard after L iterations which is considered converged when ( ) Then overall convergence is a competition between the two modes of partial sum convergence, where convergence occurs for the least relative error, ( . An isotropic source, f(µ) = 1, entering the near surface with no incoming flux entering the far surface gives the following surface flux moments vectors: From Equation (34a,b), the exiting flux moments vectors are, therefore, and the interior flux moment vectors are from Equation ( The flux vectors enter Equation (40a) re-written as the Clenshaw sum applies, where [21] ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ) We now turn to the issue of solution precision.To address precision, we compare the DPN solution to the well-established fully analytical response matrix/discrete ordinate method (RM/DOM) benchmark [4].We assume no absorption (ω = 1) for a slab of 1 mfp thickness.Table 1 shows, what are believed to be precise, angular fluxes to better than one unit in the 8th place from the RM/DOM.The DPN angular flux approximation for N = 100 agrees with all seven digits of the angular flux from the RM/DOM of quadrature order 250 and with all but three entries (the last digit is emboldened and underlined) in the 8th place when N is increased to 150 without the need for acceleration.Thus, the method presented does indeed successfully provide an extreme benchmark to nearly eight places.Figure 1 shows the effectiveness of W-e acceleration, which is the ratio of the relative errors (Re) without to with acceleration to converged over all directions at the seven spatial coordinates ( ) ( ) , ; / , ; .
One observes that the W-e relative errors at convergence are generally smaller than the errors of the original sequence for this benchmark solution.This is further confirmation of the high quality-to one unit in the seventh, and nearly one unit in the eighth place-of the proposed DPN algorithm for an isotropic source.Over half of the fluxes converged by W-e showing the significance of the Wynn-epsilon algorithm.

Numerical Implementation for a Beam Source
The final example is for a beam source entering the near surface and none at the far surface: where to be inserted into Equation (57b) and then into Equation (57a) and finally into (40a) for the angular flux.
The slab for a perpendicular beam source (µ0 = 1) is the same as for the isotropic source.Table 2 shows the results from the RM/DOM [4], which is a seven-place benchmark.Overlaid is the DPN result, where it is observed that DPN gives nearly all seven places, except for two entries in the last place.The beam source behaves differently from the isotropic source, however, as it is more sensitive to round-off errors.The DPN results shown in Table 2 required quadruple precision (QP) to achieve nearly seven places; whereas, the isotropic source required double precision.For the beam source, double precision only gave, at best, only four places.A future effort will attempt to determine where the loss of precision occurs; nevertheless, the DPN provides a solid six digits and nearly seven with QP.
All the results presented agree with all five digits of Ref. [22].

Conclusions
By expressing the solution to the 1D monoenergetic neutron transport equation in plane geometry as an infinite series of half-range Legendre polynomials, a first-order coupled set of ODEs for half-range moments in the positive and negative directions follows.Through adding and subtracting, the ODEs transform into a second-order form for the even and odd parity moments.At this point, one can choose to follow several solution schemes to solve for the parity moments.Here, we chose to establish an even parity solution to include the unknown boundary conditions.Through the manipulation of the matrix equations, one finds the relationship between the incoming and exiting fluxes through the response matrix.With the exiting fluxes known, the fluxes interior to the slab also become known.We form the solution in terms of the diagonalization of the matrix associated with the second-order ODE, or, in other words, it is expressible in eigenvectors and eigenvalues.Convergence of the infinite series occurs through the sequential convergence of the partial sums and through Wynn-epsilon convergence acceleration.Via two examples, we show that benchmarks of essentially seven and nearly eight digits are feasible.It should be noted that the precision quoted depends on the spatial positions and directions sampled but is thought to be generally representative.
to the parity equations, even in the case of zero absorption, which otherwise requires special attention.In addition, the solution features the response matrix (as the name suggests) over a homogeneous slab of any thickness and is dependent upon its material properties only.Slab responses combine to construct the response of contiguous heterogeneous media under the interaction principle, as described in [23].Additional features include convergence acceleration in quadrature order through Wynn-epsilon acceleration and faux quadrature, with the quadrature grid refinement also enabling the beam source to be in any direction.We achieved seven-place precision for the Cloud C1 300-term scattering kernel [24] with a quadrature order 2N of ~300.
2. Method of Doubling (MoD) [25] This solution uses the concepts developed by van de Hulst [26].The procedure is applied to the fully discretized 1D radiative transfer equation in both space and direction.The response of a thin slice of a discrete slab is continuously doubled until the entire discrete slab is covered.The full slab is then fully covered by adding via the interaction principle [24].Richardson and Wynn-epsilon convergence acceleration are applied to the initial thin slice and to the quadrature order, respectively.In addition, Richardson extrapolation is applied to the original slab discretization.The features found in the RM/DOM are also included.Heterogeneous media are a naturally consequence of adding.The Cloud C1 [24] case is found to be as precise as that of the RM/DOM.The advantage of doubling is that no matrix diagonalization with eigenvalues is necessary, and essentially, a host of spatial discretization schemes for evaluating the matrix exponential as a solution to a firstorder ODE are possible.Doubling demonstrates that the most basic of numerical finite difference schemes provides highly precise numerical solutions to the radiative transfer equation, thus demonstrating its simplicity.
3. The Matrix Riccati Equation Method (MREM) [27] The centerpiece of the Matrix Riccati Equation Method is a set of four non-linear matrix Riccati ODEs that represent the radiative transfer equation as a first-order ODE, similar to invariant embedding.Any two of these equations gives the interaction coefficients of reflectance and transmittance for a slab.These coefficients are found numerically from a Taylor Series (TS) representation to initiate doubling for a thin slab.Forming the interaction principle with the interaction coefficients and doubling then enables the slab response for any given thickness.With the response, the exiting and interior intensities are found.Thus, we have established a numerical method based on TS evaluation, doubling, and adding, which we wrap in the Wynn-epsilon convergence acceleration to give a highly precise solution.Results for the HAZE L and Cloud C1 [24] benchmarks are calculated to nearly seven-place precision in comparison to RM/DOM.

Figure 1 .
Figure 1.Ratio of unaccelerated to accelerated relative errors.
(scalar) component should not be confused with the zeroth vector a vector.The uncollided flux satisfies

Θ
is the Heaviside step function required to maintain the uncollided flux in the positive direction only.When we introduce Equation (52) into Equation (1a) with Equation (53c), 14c), for the collided component derived from Equation (53b),

Table 1 .
Angular flux for an entering isotropic source.