Calculating Time Eigenvalues of the Neutron Transport Equation with Dynamic Mode Decomposition

A novel method to compute time eigenvalues of neutron transport problems is presented based on solutions to the time dependent transport equation. Using these solutions we use the dynamic mode decomposition (DMD) to form an approximate transport operator. This approximate operator has eigenvalues that can be directly related to the time eigenvalues of the neutron transport equation. This approach works for systems of any level of criticality and does not require the user to have estimates for the eigenvalues. Numerical results are presented for homogeneous and heterogeneous media. The numerical results indicate that the method finds the eigenvalues that are most important to the solution evolution over a given time range, and the eigenvalue with the largest real part is not necessarily important to the system evolution.


I. INTRODUCTION
In scientific computing we are used to taking a known operator and making approximations to it. Usually, these approximations arise from the continuous operator and restrict it to some discrete representation. This is what is done in common methods for particle transport such as discrete ordinates where the continuous transport equation is replaced with equations for particular directions that are coupled through scattering via a quadrature rule.
Alternatively, it is possible to use the action of the operator to generate approximations rather than use the operator itself. This is what is done in, for example, Krylov subspace methods for solving linear systems where the action of a matrix is used to create subspaces of increasing size that are used to find approximations to the solution. The use of the known action of an operator, even if the operator is not known, is the basis for dynamic mode decomposition 1,2 (DMD).
The main idea behind DMD is that if we have a sequence of vectors generated by successively applying an operator, we can estimate properties of that operator. In fluid dynamics, DMD is used to find important modes in the evolution of a system, even when the system does not have an interesting steady state. 2,3 Additionally, because one does not need the operator, DMD can be applied to experimental measurements and quantitively compared to the DMD modes of a simulation. 4 In this paper, we use DMD to find time eigenvalues, also known as alpha eigenvalues, of the neutron transport equation using only the time-dependent solution for the angular flux. The calculation of alpha eigenvalues has traditionally been accomplished using iterative search procedures where an eigenvalue is determined by finding the value of α that makes the equivalent k-eigenvalue problem exactly critical. 5 This is accomplished by subtracting α divided by the neutron speed from the total interaction term. Unfortunately, if the alpha eigenvalue is negative, a negative total interaction term can result, leading to instabilities in most solution algorithms. Recently, there have been improvements to deterministic alpha-eigenvalue computation techniques that use specialized solvers to find positive and negative eigenvalues [6][7][8] or form the full discretization matrices to find eigenvalues. 9 Most of these methods either find only the eigenvalue with the largest real part (the rightmost eigenvalue in the complex plane) or require an accurate estimate to find other eigenvalues. Additionally, Monte Carlo can be used to find these eigenvalues with the transition rate matrix method. 10,11 The benefit of using the DMD method is that one can use standard transport solvers 12 to find any eigenvalues that are excited in a given calculation. The cost of the calculation, beyond the transport simulation, is the formation of a singular value decomposition (SVD) on the solution at several time steps. No development of transport solvers is required, and off-the-shelf linear algebra routines can be used to find the SVD. DMD will find the eigenvalues/eigenvectors that are the largest contributors to the dynamics of the system in a given time-dependent problem: This is a feature and not a bug. In many subcritical systems, the rightmost eigenvalue will be unimportant to the system behavior in a given experiment. For instance, if we consider a subcritical system struck by a pulse of neutrons, such as those in Ref. 13, there will be eigenmodes corresponding to the slowest neutrons traveling across the system 14,15 that will not impact the experiment. We will see an example of this later.
This paper is organized as follows. We begin with the presentation of the DMD in Sec. II and apply the method to time-eigenvalue problems in Sec. III. Numerical results are presented for a bare sphere in Sec. IV and for heterogeneous systems in Sec. V before presenting conclusions and future work in Sec. VI.

II. DYNAMIC MODE DECOMPOSITION
Consider an evolution equation over time that can be written in the generic form where yðr; tÞ is a function of a set of variables denoted by r, which could be space, angle, energy, etc., and time t. Consider the solution to the equation at a sequence of equally spaced times, yðr; t 0 Þ; yðr; t 1 Þ; . . . ; yðr; t NÀ1 Þ; yðr; t N Þ, separated by a time Δt. These solutions are formally determined using the exponential of the operator AðrÞ via the relationship yðr; t n Þ ¼ e AΔt yðr; t nÀ1 Þ; n ¼ 1; . . . ; N : We can write a single equation relating the solutions at each time level as ½yðr; t N Þ; yðr; t NÀ1 Þ; . . . ; yðr; t 1 Þ ¼ e AΔt ½yðr; t NÀ1 Þ; yðr; t NÀ2 Þ; . . . ; yðr; t 0 Þ : If we constrain ourselves to finite-dimensional problems, the solution is now a vector, and the operator is a matrix. In this case, the original equation has the form qy qt ¼ AyðtÞ : We will say that y n is of length M > N and A is an M Â M matrix. In this case, the solutions are related through the matrix exponential: ½y N ; y NÀ1 ; . . . ; y 1 ¼ e AΔt ½y NÀ1 ; y NÀ2 ; . . . ; y 0 : ð4Þ In shorthand, we can define the N X M matrices Y þ ¼ ½y N ; y NÀ1 ; . . . ; y 1 and Y À ¼ ½y NÀ1 ; y NÀ2 ; . . . ; y 0 as the matrices formed by appending the column vectors y n . This leads to the relation Equation (5) is exact; however, the matrix A may be too large to compute the exponential e AΔt . Therefore, we desire to use just the solution to estimate the eigenvalues of e AΔt .
To this end, we will use the solution vectors collected in Y þ and Y À to produce an approximation to A. We compute the thin SVD of the matrix Y À : Typically, some of the diagonal elements of Σ are effectively zero. Therefore, we make Σ the r Â r matrix that contains all r values greater than some small, positive .
Substituting Eq. (6) into Eq. (5), we get Rearranging Eq. (6) gives It has been shown 16 that an eigenvalue of e S is also an eigenvalue of e AΔt . To see this, we consider an eigenvalue λ and eigenvector v of e S. By definition, we have e Sv ¼ λv; which is equivalent to U Ã e AΔt Uv ¼ λv: Left multiplying this equation by U; we get e AΔt Uv ¼ λUv ; which shows that λ is an eigenvalue of e AΔt . Additionally, b v ¼ Uv is the associated eigenvector of e AΔt to eigenvalue λ.
The matrix e S is much smaller than that for e AΔt ; and we can form e S without forming A. To create e S; we need to know the result of e AΔt applied to an initial condition several times in succession. Then, we need to compute the SVD of the data matrix Y À . A direct computation requires OðM 2 NÞ operations, though iterative methods for computing the SVD exist. 17 As a comparison, the QR factorization of e AΔt requires OðM 3 Þ operations. Our formulation here requires a constant time step size, though this can be relaxed as shown by Tu et al. 16

III. ALPHA EIGENVALUES OF THE TRANSPORT OPERATOR
We will now demonstrate that we can estimate the alpha eigenvalues of a nuclear system by computing several time steps of a time-dependent transport equation and using the DMD theory presented above to form and compute the eigenvalues of e S. We begin by defining the alpha-eigenvalue transport problem without delayed neutrons.
Consider the time-dependent transport equation 18 where ψðx; Ω; E; tÞ is the angular flux at position x 2 R 3 , in direction Ω 2 S 2 , at energy E and time t. The transport operator A is given by with S and F the scattering and fission operators: where σ s ðΩ 0 ! Ω; E 0 ! EÞ = double-differential scattering cross section from direction Ω 0 and energy E 0 to direction Ω and energy E νσ f ðE 0 Þ = fission cross section times the expected number of fission neutrons at energy E 0 χðEÞ = probability of a fission neutron being emitted with energy E.
The scalar flux ϕðx; Ω 0 ; E 0 ; tÞ is defined as the integral of the angular flux over the unit sphere: Above, we used a continuous formulation of the transport problem. For our calculations later, we will use a discretized transport equation using the multigroup method 18 in energy, discrete ordinates in angle, and a spatial discretization. In this case the time-dependent transport equation can be written as a system of differential equations: where Ψ is a vector and A is a matrix that represents the discrete transport operator. To define alpha eigenvalues and eigenfunctions, consider a solution of the form b ψðx; Ω; EÞe αt , which, using Eq. (8), leads to the relation The values of α where this relation holds are called alpha eigenvalues and b ψ are the alpha eigenfunctions. In discrete form the alpha eigenvalue problem is where Ψ has the form b ψe αt . In general, the eigenvalues of the discrete problem are not the same as those for the continuous problem due to discretization. From here on, we consider the discrete problem.
In the alpha-eigenvalue problem, we are interested in the eigenvalues of A. We can use the DMD decomposition to form the operator e S and compute its eigenvalues and, as a result, the eigenvalues of e AΔt . To do this, we begin with an initial condition and compute the solution at N time steps. Then, we can form Y þ and Y À , compute the SVD, and get the eigenvalues of e AΔt .
We need a way to relate the eigenvalues of e AΔt to the alpha eigenvalues. The relationship is if ðα; vÞ is an eigenvalue/eigenvector pair of A; then e αΔt is an eigenvalue of e AΔt with eigenvector v. These facts can be seen through the definition of the matrix exponential. Consider an eigenvalue α with eigenvector v for the matrix A. Using the definition of an eigenvector, we can show that The definition of the matrix exponential gives where the last equality uses the Taylor series of the exponential function expðαΔtÞ around 0. Therefore, if λ is an eigenvalue of e S and, by construction, an eigenvalue of e AΔt , then is an alpha eigenvalue of the discrete transport operator. The discussion above suggests the following algorithm for estimating alpha eigenvalues of the discrete transport equation: 1. Compute N time-dependent steps starting from ψ 0 using a numerical method of choice and fixed Δt.
2. Compute the SVD of the resulting data matrix Y À , and form e S.
3. Compute the eigenvalues/eigenvectors of e S, and calculate the alpha eigenvalues from Eq. (14). This is an approximate method because the time steps typically will not be computed using the matrix exponential; rather, a time integration technique such as the backward Euler method will be used. The backward Euler algorithm estimates the matrix exponential as e AΔt % ðI À ΔtAÞ À1 : When we use the DMD method on a data matrix generated by the backward Euler method, we are computing eigenvalues of ðI À ΔtAÞ À1 : To relate these eigenvalues to the alpha eigenvalues, we use the relation This approximation will improve at first order as Δt ! 0.

III.A. Comparison with Existing Methods
Standard techniques for computing alpha eigenvalues require solving a series of k-eigenvalue problems. 5 The basis for these methods is that the alpha eigenvalues make the equivalent k-eigenvalue problem exactly critical when the total cross section is replaced with σ t ðEÞ þ αvðEÞ À1 . This approach will have problems when α is negative as it can cause negative absorption to arise in lower-energy groups.
To address this problem, other methods have been developed such as Rayleigh quotient methods, 6 the Arnoldi method, 7,19 and Newton-Krylov methods. 20 In these approaches, the equations that need to be solved are typically different from those required to solve time-dependent transport problems. The DMD method allows one to get both the time-dependent solution and eigenvalues as part of one calculation. Moreover, DMD

NEUTRON TRANSPORT EQUATION TIME EIGENVALUES · McCLARREN 857
provides an estimate for multiple eigenvalues based on the number of modes excited in the system and the number of steps used.

IV. RESULTS FOR PLUTONIUM SPHERE
Here, we present results for the prompt neutron solution for a sphere of 99 at. % 239 Pu and 1 at. % natural carbon using 12-group cross sections and a simple buckling model for leakage so that we can solve an infinite medium problem. The group structure is detailed in Table I. We will consider subcritical and supercritical systems by adjusting the radius of the sphere. Because we use a simple buckling model for this problem, we can directly form the matrix for the transport operator and compute exact eigenvalues for this model. For DMD, the time steps are computed using the backward Euler discretization for time integration. In this and subsequent sections, we consider only prompt neutrons.

IV.A. Subcritical Case
We consider a sphere of radius 4.77178 cm with an associated k eff in our model of 0.95000. The fundamental mode for this reactor is shown in Fig. 1a along with several alpha eigenmodes. The alpha eigenvalues for this system have a fast-decaying mode with a large number of neutrons in the fastest energy group, and the most slowly decaying mode closely follows the fundamental mode.
To test the DMD estimation of alpha eigenvalues, we run a time-dependent problem where at time zero the system has 1000 neutrons in the energy group corresponding to 14.1 MeV. This is a crude approximation to an experiment where a pulse of deuterium-tritium (D-T) fusion neutrons irradiates the sphere. The problem is run in time-dependent mode out to various final times with uniform time steps, and the time steps are used in the DMD procedure to estimate alpha eigenvalues. The alpha eigenvalues computed by DMD are shown in Table II and compared to the exact eigenvalues computed from the matrices generated by the buckling approximation. The number of neutrons in the system as a function of time is shown in Fig. 2, where one can see that subcritical multiplication is happening in the first 0.002 μs of the problem. As we argue next, DMD finds the eigenvalues that are important in the time-dependent solution over the timescales considered and that are resolved by the time-step size.
From Table II, we can see that during the phase where subcritical multiplication is occurring (before t ¼ 0:002 μs), DMD accurately computes to six digits the alpha eigenmode that corresponds to a large population of 14.1-MeV neutrons. This is the mode most excited by the initial condition. It also accurately computes the eigenvalues with magnitudes larger than 200 to several digits. However, we note that the dominant or most slowly decaying eigenmode is not detected by the DMD algorithm, indicating that its contribution at this early time is insignificant or cannot be distinguished from other slowly decaying modes. This indicates an important phenomenon in time-dependent transport: The most slowly decaying eigenvalue may not be important in a given problem.
As we look at simulations run to later time, more eigenvalues are identified using DMD. Running the simulation to intermediate times, 0.02 and 0.2 μs, we see that DMD finds all of the eigenvalues in the problem to several digits of accuracy. In both of these solutions, DMD does not find the eigenvalue near À28:85 μs À1 . This eigenmode has more neutrons in the energy ranges in the thermal and epithermal energy ranges relative to the other modes. Given that this problem has very little thermalization due to the small amount of carbon, this mode is not important at these intermediate times relative to other modes.
At a much later time, 2 μs, DMD identifies all of the slowly decaying modes but cannot find the rapidly decaying modes. This is because the larger time steps used make it so that the solution cannot resolve the timescale where these modes are important. As a result, DMD estimates a pair of complex eigenvalues with a real part that does not correspond to an actual eigenvalue. There are versions of DMD that allow variable time steps to be used, 16

IV.B. Supercritical Case
We consider a sphere of radius 5.029636 cm with an associated prompt k eff in our model of 1.000998; the eigenvectors for this problem are shown in Fig. 1b. We perform the same calculations as performed before on the subcritical sphere. Table III compares

NEUTRON TRANSPORT EQUATION TIME EIGENVALUES · McCLARREN 859
where modes associated with the fusion neutrons contributed to the growth of the neutron population is also present in this supercritical system. However, there are very few neutrons emitted in the fusion energy range from fission (χ 1 % 1:37 Â 10 À4 ), so these modes decay away. As the solution time increases, the DMD-estimated eigenvalues agree well with the true values. This is most evident in the solution computed up to 0:2 μs, where 11 of 12 eigenvalues are computed accurately to two digits. The exponentially growing mode is correctly estimated at later times; for the simulation run to the latest time, the eigenvalue is estimated accurately to six digits by DMD. At very late times, the rapidly decaying modes are not correctly estimated, and a complex eigenvalue is estimated, as we saw before in the subcritical case, but this is likely due to the large time step used.

V. HETEROGENEOUS MEDIA
The plutonium sphere example required computing only the solution to infinite medium problems. We will now investigate how the DMD approach to estimating eigenvalues performs on heterogeneous problems in slab geometry. Our numerical solutions are computed using the discrete ordinates (S N ) method with diamond difference for the spatial discretization and backward Euler for time integration. 12

V.A. Heterogeneous, One-Speed Slab Problem
The first heterogeneous problem we solve is based on benchmark problems published by Kornreich and Parsons 21 as solved by the Green's function method (GFM). Their work defines a slab problem for single-speed neutrons (i.e., one group) consisting of an absorber surrounded by a moderator and fuel; see Fig. 3. They define configurations of this problem that are symmetric and asymmetric, as well as subcritical and supercritical versions. In the symmetric version of the problem, the total width of the slab is 9, whereas in the asymmetric version, the width is 9.1. The total cross section is one throughout the problem, and the scattering cross sections are σ s ðxÞ ¼ 0:8 x 2 fuel or moderator 0:1 x 2 absorber : The value of νσ f in the fuel is either 0:3 or 0:7 for the subcritical and supercritical cases, respectively. We solve this problem using DMD with 200 cells per mean free path and a 196-angle Gauss-Legendre quadrature set. We use a time step size of Δt ¼ 0:1 and run the problem for 500 time steps to a time of t ¼ 50. For initial conditions, we used two approaches: a symmetric initial condition where the solution is nonzero and inwardly directed in the outermost cells in the problem, and a random initial condition. In Table IV, results from the DMD calculations are compared with the GFM results. We use the nomenclature of "Fundamental" for the alpha eigenvalue that is rightmost in the complex plane to coincide with the published results; the "Second" eigenvalue in Table IV is the eigenvalue that is just left of the fundamental eigenvalue in the complex plane. The results in Table IV show that the DMD results were able to  3. Layout for the multiregion slab problem from Kornreich and Parsons. 21 The total width of the problem, X , can be either 9 or 9:1.

NEUTRON TRANSPORT EQUATION TIME EIGENVALUES · McCLARREN 861
reproduce the GFM eigenvalues within 10 À5 (1 pcm). Except for the second eigenvalue in the symmetric case, all the DMD eigenvalues agreed to greater than 1 pcm precision using both initial conditions. The DMD results in Table IV for the fundamental eigenvalue were the same for both initial conditions to six significant digits. We also have found that the eigenvalues found in the solution are insensitive to the number of time steps used in the DMD procedure, as long as any initial transients have died out (about 5 mean free times in this problem). Using 400 or 100 time steps in the eigenvalue estimate gave the same eigenvalue estimates to six significant digits. However, the second eigenvalue was not present in the solution for the symmetric initial condition on the symmetric problems. This is because the second eigenmode is asymmetric in space, and therefore, this mode is not excited by the symmetric initial condition. The DMD eigenvectors for the four configurations of this problem are shown in Figs. 4 and 5. The fundamental and second eigenvectors match the published plots for the νσ f ¼ 0:7 within the width of the lines. In the DMD results, we found a third, real-valued eigenvalue, α ¼ À1:02158875. This eigenvalue is part of the continuum spectrum for the transport operator for this problem. The fact that it is found by DMD is an artifact of the approximations made in the method.
We note that in the original paper by Kornreich and Parsons 21 they give results from the PARTISN discrete ordinates code 22 using 96 quadrature points (about half of what we used) and 2000 mesh cells per mean free path (ten times higher resolution than in our case). The PARTISN results agreed with the GFM results to within 0.1 pcm using this much finer spatial grid. Nevertheless, PARTISN was not able to estimate the second eigenvalue in the asymmetric cases, whereas the DMD results are as expected. Furthermore, the MCNP Monte Carlo transport code 23 was not able to estimate eigenvalues for any of the νσ f ¼ 0:3 cases. Recently, Betzler et al. 10 published Monte Carlo results for these cases using the Monte Carlo Markov Transition Rate Matrix Method.

V.B. Multiregion, 70-Group System
As a final demonstration, we solve a problem consisting of two slabs of 239 Pu with high-density polyethylene (HDPE) between them and a reflector of HDPE on the outside. The initial condition has a pulse of D-T fusion neutrons striking the outer surface of the reflector, implemented as the angular flux for each angle directed toward the center being set to 1 in the outermost cell on each side for the initial condition. See Fig. 6 for  The behavior of the neutron population in time, as well as the three time intervals over which the eigenvalues were estimated, is shown in Fig. 8a. The time interval from 0.002 to 0.004 μs is during the subcritical multiplication phase of the simulation. It makes sense that during this phase the slowly decaying modes are not important in the solution. Later in time, these slowly decaying modes will dominate because the subcritical multiplication must end at some point given that the system is subcritical and does not have a fixed source.
In Fig. 8, we show the neutron spectrum at several points in space. The spectra shown are computed using time steps

NEUTRON TRANSPORT EQUATION TIME EIGENVALUES · McCLARREN 863
NUCLEAR SCIENCE AND ENGINEERING · VOLUME 193 · AUGUST 2019

864
McCLARREN · NEUTRON TRANSPORT EQUATION TIME EIGENVALUES from the indicated time ranges. From Fig. 8, we can see that early in the time, the solution is dominated by the presence of 14.1-MeV neutrons, though fission neutrons are present in the fuel and outer reflector. At late times, near 1 μs, the spectrum in the fuel and the reflector is close to the fundamental eigenmode of the k-eigenvalue problem. Nevertheless, the central moderator in the problem has not reached the fundamental k eigenmode, as there has not been enough time to fully thermalize the neutrons. Additionally, the eigenvalue for the most slowly decaying mode is associated with the travel time of the slowest neutrons crossing the moderator. This suggests that the problem would need to be run longer to relax to this mode. Moreover, it indicates that if this system were involved in an experiment, the neutrons produced in the first microsecond would give little information about the spectrum of the keigenvalue problem. Figure 9 shows the spatial distribution of neutrons. From Fig. 9, we see that at different times for the most slowly decaying mode, the DMD estimates correspond to the modes that are important to the dynamics during a time interval. Early in time, fast neutrons dominate; these fast neutrons then decay as more thermal neutrons are created from scattering. Nevertheless, near 1 μs, the the neutron density of epithermal neutrons is still larger than the density of thermal neutrons.

VI. DISCUSSION
Dynamic mode decomposition allows for the approximation of the eigenvalues present in a timedependent transport system from the solution at different times without a separate eigenvalue solve.

NEUTRON TRANSPORT EQUATION TIME EIGENVALUES · McCLARREN 865
The decomposition works for subcritical and critical systems and can give highly accurate (sub-pcm) estimates of eigenvalues. Our results from a variety of problem types indicate that the method is useful for general estimation of system eigenvalues, especially if one is interested in the modes driving the dynamics over a particular time interval. The problems we presented did not include delayed neutrons, but adding these to the DMD method is straightforward. Because DMD uses the solution from time-dependent transport to estimate eigenvalues, the time interval considered and the time-step size affect the eigenvalues found.
For instance, at early times of the simulation, there may be different modes present than at later times. DMD will not be able to accurately estimate modes that decay much more quickly than the time step size used to generate the time-dependent solution.
We note that DMD can be applied to nonlinear problems in the same fashion as we applied it to the linear problem of neutron transport. This could be useful for the situation where the neutron population dynamics are nonlinear. For instance, if we consider a system with negative feedback with regard to temperature, the dynamics of the neutron population would affect the temperature and the cross sections of the material. One could apply DMD to this problem, though the interpretation of the resulting eigenvalues would necessarily be different. Previous work, 1,24 has shown that the modes computed by DMD will be eigenfunctions of the Koopman operator, and the application of this type of analysis could be fruitful for understanding nuclear systems.