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

Abstract

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 to form an approximate transport operator. This approximate operator has eigenvalues that are mathematically 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 contribute the most to the change in the solution over a given time range, and the eigenvalue with the largest real part is not necessarily important to the system evolution at short and intermediate times.

Keywords:

• Time eigenvalues
• neutron transport
• time-dependent transport
• dynamic mode decomposition

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.⁴

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 $\alpha$ that makes the equivalent $k$-eigenvalue problem exactly critical.⁵ This is accomplished by subtracting $\alpha$ 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–8 or form the full discretization matrices to find eigenvalues.⁹ 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¹² 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

(1) $\frac{\mathrm{\partial }y}{\mathrm{\partial }t}=A(r)y(r,t),$(1)

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),\dots ,y(r,t_{N-1}),y(r,t_N)$, separated by a time $\mathrm{\Delta }t$. These solutions are formally determined using the exponential of the operator $A(r)$ via the relationship

$y(r,t_n)=e{}^{A\mathrm{\Delta }t}y(r,t_{n-1}),n=1,\dots ,N.$

We can write a single equation relating the solutions at each time level as

(2) $[y(r,t_N),y(r,t_{N-1}),\dots ,y(r,t_1)]=e{}^{A\mathrm{\Delta }t}[y(r,t_{N-1}),y(r,t_{N-2}),\dots ,y(r,t_0)].$(2)

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

(3) $\frac{\mathrm{\partial }\mathbf{y}}{\mathrm{\partial }t}=\mathbf{A}\mathbf{y}(t).$(3)

We will say that ${\mathbf{y}}_n$ is of length $M>N$ and $\mathbf{A}$ is an $M\times M$ matrix. In this case, the solutions are related through the matrix exponential:

(4) $[{\mathbf{y}}_N,{\mathbf{y}}_{N-1},\dots ,{\mathbf{y}}_1]=e^{\mathbf{A}\mathrm{\Delta }t}[{\mathbf{y}}_{N-1},{\mathbf{y}}_{N-2},\dots ,{\mathbf{y}}_0].$(4)

In shorthand, we can define the $NXM$ matrices

${\mathbf{Y}}_{+}=[{\mathbf{y}}_N,{\mathbf{y}}_{N-1},\dots ,{\mathbf{y}}_1]$

and

${\mathbf{Y}}_{-}=[{\mathbf{y}}_{N-1},{\mathbf{y}}_{N-2},\dots ,{\mathbf{y}}_0]$

as the matrices formed by appending the column vectors ${\mathbf{y}}_n$. This leads to the relation

(5) ${\mathbf{Y}}_{+}=e^{\mathbf{A}\mathrm{\Delta }t}{\mathbf{Y}}_{-}.$(5)

Equation (5) is exact; however, the matrix $\mathbf{A}$ may be too large to compute the exponential $e^{\mathbf{A}\mathrm{\Delta }t}$. Therefore, we desire to use just the solution to estimate the eigenvalues of $e^{\mathbf{A}\mathrm{\Delta }t}$.

To this end, we will use the solution vectors collected in ${\mathbf{Y}}_{+}$ and ${\mathbf{Y}}_{-}$ to produce an approximation to $\mathbf{A}$. We compute the thin SVD of the matrix ${\mathbf{Y}}_{-}$:

(6) ${\mathbf{Y}}_{-}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^{\ast },$(6)

where

$\mathbf{U}$ =	=	$M\times N$ unitary matrix
$\mathbf{V}$ =	=	$N\times N$ unitary matrix
$\mathbf{\Sigma }$ =	=	$N\times N$ diagonal matrix with nonnegative elements,

and the asterisk (*) denotes the conjugate transpose of a matrix.

Typically, some of the diagonal elements of $\mathbf{\Sigma }$ are effectively zero. Therefore, we make $\mathbf{\Sigma }$ the $r\times r$ matrix that contains all $r$ values greater than some small, positive $ϵ$.

Substituting Eq. (6) into Eq. (5), we get

${\mathbf{Y}}_{+}=e^{\mathbf{A}\mathrm{\Delta }t}\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^{\ast }.$

Rearranging Eq. (6) gives

(7) ${\mathbf{U}}^{\ast }{\mathbf{Y}}_{+}\mathbf{V}{\mathbf{\Sigma }}^{-1}={\mathbf{U}}^{\ast }{e}^{\mathbf{A}\mathrm{\Delta }t}\mathbf{U}\equiv \tilde{\mathbf{S}}.$(7)

It has been shown¹⁶ that an eigenvalue of $\tilde{\mathbf{S}}$ is also an eigenvalue of $e^{\mathbf{A}\mathrm{\Delta }t}$. To see this, we consider an eigenvalue $\lambda$ and eigenvector $\mathbf{v}$ of $\tilde{\mathbf{S}}$. By definition, we have $\tilde{\mathbf{S}}\mathbf{v}=\lambda \mathbf{v},$ which is equivalent to ${\mathbf{U}}^{\ast }{e}^{\mathbf{A}\mathrm{\Delta }t}\mathbf{U}\mathbf{v}=\lambda \mathbf{v}.$ Left multiplying this equation by $\mathbf{U},$ we get

$e^{\mathbf{A}\mathrm{\Delta }t}\mathbf{U}\mathbf{v}=\lambda \mathbf{U}\mathbf{v},$

which shows that $\lambda$ is an eigenvalue of $e^{\mathbf{A}\mathrm{\Delta }t}$. Additionally, $\stackrel{\mathbf{ˆ}}{\mathbf{v}}=\mathbf{U}\mathbf{v}$ is the associated eigenvector of $e^{\mathbf{A}\mathrm{\Delta }t}$ to eigenvalue $\lambda$.

The matrix $\tilde{\mathbf{S}}$ is much smaller than that for $e^{\mathbf{A}\mathrm{\Delta }t},$ and we can form $\tilde{\mathbf{S}}$ without forming $\mathbf{A}$. To create $\tilde{\mathbf{S}},$ we need to know the result of $e^{\mathbf{A}\mathrm{\Delta }t}$ applied to an initial condition several times in succession. Then, we need to compute the SVD of the data matrix ${\mathbf{Y}}_{-}$. A direct computation requires $O(M^{2}N)$ operations, though iterative methods for computing the SVD exist.¹⁷ As a comparison, the QR factorization of $e^{\mathbf{A}\mathrm{\Delta }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.¹⁶

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 $\tilde{\mathbf{S}}$. We begin by defining the alpha-eigenvalue transport problem without delayed neutrons.

Consider the time-dependent transport equation¹⁸

(8) $\frac{\mathrm{\partial }\psi }{\mathrm{\partial }t}=A\psi ,$(8)

where $\psi (x,\mathrm{\Omega },E,t)$ is the angular flux at position $x\in {\mathbb{R}}^3$, in direction $\mathrm{\Omega }\in {\mathbb{S}}_2$, at energy $E$ and time $t$. The transport operator $A$ is given by

$A=v(E)(-\mathrm{\Omega }\cdot \mathrm{\nabla }+-{\sigma }_t+\mathcal{S}+\mathcal{F}),$

with $\mathcal{S}$ and $\mathcal{F}$ the scattering and fission operators:

(9) $\begin{array}{rl}& \mathcal{S}\psi ={\int }_{4\pi }d\mathrm{\Omega }\prime {\int }_0^{\mathrm{\infty }}dE\prime {\sigma }_s(\mathrm{\Omega }\prime \to \mathrm{\Omega },E\prime \to E)\\ & \times \psi (x,\mathrm{\Omega }\prime ,E\prime ,t)(9)\end{array}$(9)

and

(10) $\mathcal{F}\psi =\frac{\chi (E)}{4\pi }{\int }_0^{\mathrm{\infty }}dE\prime \nu {\sigma }_f(E\prime )\varphi (x,\mathrm{\Omega }\prime ,E\prime ,t),$(10)

where

${\sigma }_s(\mathrm{\Omega }\prime \to \mathrm{\Omega },E\prime \to E)$ =	=	double-differential scattering cross section from direction $\mathrm{\Omega }\prime$ and energy $E\prime$ to direction $\mathrm{\Omega }$ and energy $E$
$\nu {\sigma }_f(E\prime )$ =	=	fission cross section times the expected number of fission neutrons at energy $E\prime$
$\chi (E)$ =	=	probability of a fission neutron being emitted with energy $E$.

The scalar flux $\varphi (x,\mathrm{\Omega }\prime ,E\prime ,t)$ is defined as the integral of the angular flux over the unit sphere:

(11) $\varphi (x,E,t)={\int }_{4\pi }d\mathrm{\Omega }\psi (x,\mathrm{\Omega },E,t).$(11)

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¹⁸ 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:

(12) $\frac{\mathrm{\partial }\mathrm{\Psi }}{\mathrm{\partial }t}=\mathbf{A}\mathrm{\Psi },$(12)

where $\mathrm{\Psi }$ is a vector and $\mathbf{A}$ is a matrix that represents the discrete transport operator.

To define alpha eigenvalues and eigenfunctions, consider a solution of the form $\hat{\psi }(x,\mathrm{\Omega },E)e^{\alpha t}$, which, using Eq. (8), leads to the relation

$\alpha \hat{\psi }=A\hat{\psi }.$

The values of $\alpha$ where this relation holds are called alpha eigenvalues and $\hat{\psi }$ are the alpha eigenfunctions. In discrete form the alpha eigenvalue problem is

$\mathbf{A}\hat{\psi }=\alpha \hat{\psi },$

where $\mathrm{\Psi }$ has the form $\hat{\psi }{e}^{\alpha 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 $\mathbf{A}$. We can use the DMD decomposition to form the operator $\tilde{\mathbf{S}}$ and compute its eigenvalues and, as a result, the eigenvalues of $e^{\mathbf{A}\mathrm{\Delta }t}$. To do this, we begin with an initial condition and compute the solution at $N$ time steps. Then, we can form ${\mathbf{Y}}_{+}$ and ${\mathbf{Y}}_{-}$, compute the SVD, and get the eigenvalues of $e^{\mathbf{A}\mathrm{\Delta }t}$.

We need a way to relate the eigenvalues of $e^{\mathbf{A}\mathrm{\Delta }t}$ to the alpha eigenvalues. The relationship is if $(\alpha ,\mathbf{v})$ is an eigenvalue/eigenvector pair of $\mathbf{A},$ then $e^{\alpha \mathrm{\Delta }t}$ is an eigenvalue of $e^{\mathbf{A}\mathrm{\Delta }t}$ with eigenvector $\mathbf{v}$. These facts can be seen through the definition of the matrix exponential. Consider an eigenvalue $\alpha$ with eigenvector $\mathbf{v}$ for the matrix $\mathbf{A}$. Using the definition of an eigenvector, we can show that

${\mathbf{A}}^{\mathrm{\ell }}\mathbf{v}={\mathbf{A}}^{\mathrm{\ell }-1}(\alpha \mathbf{v})={\mathbf{A}}^{\mathrm{\ell }-2}({\alpha }^2\mathbf{v})=\dots ={\alpha }^{\mathrm{\ell }}\mathbf{v}.$

The definition of the matrix exponential gives

(13) $\begin{array}{rl}& e^{\mathbf{A}\mathrm{\Delta }t}\mathbf{v}=\left(\sum _{\mathrm{\ell }=0}^{\mathrm{\infty }}\frac{\mathrm{\Delta }{t}^{\mathrm{\ell }}}{\mathrm{\ell }!}{\mathbf{A}}^{\mathrm{\ell }}\right)\mathbf{v}\\ & =\left(\sum _{\mathrm{\ell }=0}^{\mathrm{\infty }}\frac{\mathrm{\Delta }{t}^{\mathrm{\ell }}}{\mathrm{\ell }!}{\alpha }^{\mathrm{\ell }}\right)\mathbf{v}\\ & =e^{\alpha \mathrm{\Delta }t}\mathbf{v},(13)\end{array}$(13)

where the last equality uses the Taylor series of the exponential function $exp(\alpha \mathrm{\Delta }t)$ around 0.

Therefore, if $\lambda$ is an eigenvalue of $\tilde{\mathbf{S}}$ and, by construction, an eigenvalue of $e^{\mathbf{A}\mathrm{\Delta }t}$, then

(14) $\alpha =\frac{log\lambda }{\mathrm{\Delta }t}$(14)

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 ${\psi }_0$ using a numerical method of choice and fixed $\mathrm{\Delta }t$.

2. Compute the SVD of the resulting data matrix ${\mathbf{Y}}_{-}$, and form $\tilde{\mathbf{S}}$.

3. Compute the eigenvalues/eigenvectors of $\tilde{\mathbf{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^{\mathbf{A}\mathrm{\Delta }t}\approx (\mathbf{I}-\mathrm{\Delta }t\mathbf{A}{)}^{-1}.$

When we use the DMD method on a data matrix generated by the backward Euler method, we are computing eigenvalues of $(\mathbf{I}-\mathrm{\Delta }t\mathbf{A}{)}^{-1}.$ To relate these eigenvalues to the alpha eigenvalues, we use the relation

$\alpha \approx \frac{1}{\mathrm{\Delta }t}\left(1-\frac{1}{\lambda }\right).$

This approximation will improve at first order as $\mathrm{\Delta }t\to 0$.

III.A. Comparison with Existing Methods

Standard techniques for computing alpha eigenvalues require solving a series of $k$-eigenvalue problems.⁵ 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 ${\sigma }_t(E)+\alpha v(E{)}^{-1}$. This approach will have problems when $\alpha$ 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,⁶ the Arnoldi method,^7,19 and Newton-Krylov methods.²⁰ 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 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.

TABLE I The Group Edges and Centers for the 12-Group Calculations in This Study

$g$	$E_g$ (MeV)	${\bar{E}}_g$ (MeV)
0	17
1	13.5	15.25
2	10	11.75
3	6.07	8.035
4	2.865	4.4675
5	1.353	2.109
6	0.5	0.9265
7	0.184	0.342
8	0.0676	0.1258
9	0.0248	0.0462
10	0.00912	0.01696
11	0.00335	0.006235
12	0.000454	0.001902

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.

Fig. 1. Fundmental $k$-eigenmode, and several alpha eigenmodes for the bare plutonium sphere problem with 12 groups in subcritical and supercritical configurations. The alpha eigenvalues are in units of inverse microsecond.

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.

TABLE II Alpha Eigenvalues for the Subcritical Sphere Computed Using DMD Using the Solution Obtained with Different Values of $\mathrm{\Delta }t$ and Final Times*

Exact	$t_{final}$ = 0.002 μs	$t_{final}$ = 0.02 μs	$t_{final}$ = 0.2 μs	$t_{final}$ = 2 μs
−17.7439	–	−17.5504	−17.7588	−17.7437
−28.8533	−24.5669	–	–	−28.8628
−34.4201	–	−35.7281	−34.1948	−34.3999
−48.4269	–	−46.6817	−48.0231	−48.4613
−75.0701	–	−75.7798	−75.2787	−74.9998
−132.352	–	−132.183	−132.197	−132.587
−261.942	−262.78	−261.974	−262.127	−260.218
−547.732	−531.575	−547.719	−547.11	−585.536
−893.385	−893.314	−893.399	−895.262	−763.974
−1368.92	−1335.16	−1368.90	−1362.45	–
−1732.99	−1721.75	−1733.01	−1725.84	−1708 ± $381i$
−1957.42	−1957.42	−1957.41	−1957.42	–
	$\mathrm{\Delta }t$ = 0.0002 μs	0.0002 μs	0.001 μs	0.01 μs

*The alpha eigenvalues are in units of inverse microseconds.

Fig. 2. The number of neutrons in the plutonium sphere in subcritical and supercritical configurations as a function of time. Because of subcritical multiplication, the peak number occurs about 0.002 μs into the simulation of the subcritical configuration.

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,¹⁶ and the use of adaptive time stepping should be investigated in future work in order to estimate the fast-decaying and slowly decaying modes.

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 the eigenvalues computed with DMD with the eigenvalues computed by solving the equivalent infinite medium problem. At an early time (0.002 μs), the DMD computation does not identify the exponentially increasing mode. Upon inspection of Fig. 2, we see that at this time the supercritical and subcritical systems have neutron populations that are very similar. The subcritical multiplication observed in the smaller sphere 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 (${\chi }_1\approx 1.37\times {10}^{-4}$), so these modes decay away.

TABLE III Alpha Eigenvalues for the Supercritical Sphere Computed Using DMD Using the Solution Obtained with Different Values of $\mathrm{\Delta }t$ and Final Times*

Exact	$t_{final}$ = 0.002 μs	$t_{final}$ = 0.02 μs	$t_{final}$ = 0.2 μs	$t_{final}$ = 2 μs
0.354439	−4.02079	0.332366	0.354291	0.354439
−28.2933	–	–	–	−28.2932
−33.1095	–	−32.8048	−33.1151	–
−46.0832	–	−45.3512	−45.817	−46.0703
−70.7945	–	−70.4805	−70.9448	−70.8261
−124.497	–	−124.568	−124.38	−124.381
−247.14	−247.914	−247.127	−247.281	−248.057
−521.689	−506.467	−521.693	−521.216	−507.684
−853.58	−853.733	−853.577	−855.008	–
−1309.4	−1279.91	−1309.4	−1305.12	−1050 $+23i$
−1659.02	−1649.68	−1659.02	−1655.16	–
−1872.99	−1872.99	−1872.99	−1872.98	−2059.43
	$\mathrm{\Delta }t$ = 0.0002 μs	0.0002 μs	0.001 μs	0.01 μs

*The alpha eigenvalues are in units of inverse microseconds.

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.¹²

V.A. Heterogeneous, One-Speed Slab Problem

The first heterogeneous problem we solve is based on benchmark problems published by Kornreich and Parsons²¹ 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

${\sigma }_s(x)=\left\{\begin{array}{cc}0.8& x\in fuelormoderator\\ 0.1& x\in absorber.\end{array}\right.$

Fig. 3. Layout for the multiregion slab problem from Kornreich and Parsons.²¹ The total width of the problem, $X$, can be either $9$ or $9.1$.

The value of $\nu {\sigma }_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 $\mathrm{\Delta }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 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.

TABLE IV Eigenvalues for the Benchmark as Computed via the GFM and the Difference Between the GFM and DMD Estimates in pcm (${10}^{-5}$)

Geometry	$\nu {\sigma }_f^{fuel}$	Fundamental $\alpha$ (GFM)	${\alpha }_{\mathrm{G}\mathrm{F}\mathrm{M}}-{\alpha }_{\mathrm{D}\mathrm{M}\mathrm{D}}$ (pcm)	Second $\alpha$ (GFM)	${\alpha }_{\mathrm{G}\mathrm{F}\mathrm{M}}-{\alpha }_{\mathrm{D}\mathrm{M}\mathrm{D}}$ (pcm)
Symmetric	0.3	−0.3196537	0.639	−0.3229855	0.694
	0.7	−0.006156369	0.7711	−0.006440766	0.7724
Asymmetric	0.3	−0.2932468	0.535	−0.3213939	0.666
	0.7	0.03759991	0.64	−0.006298843	0.7717

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 $\nu {\sigma }_f=0.7$ within the width of the lines. In the DMD results, we found a third, real-valued eigenvalue, $\alpha =-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.

Fig. 4. Fundamental and second eigenmodes for the one-group slab problem in the symmetric configurations.

Fig. 5. Fundamental and second eigenmodes for the one-group slab problem in the asymmetric configurations.

We note that in the original paper by Kornreich and Parsons²¹ they give results from the PARTISN discrete ordinates code²² 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²³ was not able to estimate eigenvalues for any of the $\nu {\sigma }_f=0.3$ cases. Recently, Betzler et al.¹⁰ 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 a schematic of the problem. The system is subcritical when the fuel regions are each 1.125 cm thick with a resulting $k_{eff}\approx 0.97$ and isotropic scattering is assumed. The fundamental mode has a large number of thermal neutrons in the middle of the problem as well as a fast peak in the fuel region.

Fig. 6. Problem layout for the 70-group test problem.

Running this problem out to a time of 1 μs with a time step size of ${10}^{-4}$ μs, $S_8$ quadrature, and 400 spatial zones, we use DMD to compute eigenvalues present in the solution over three different time windows: 0.002 to 0.004 μs, 0.09 to 0.1 μs, and 0.99 to 1 μs. These eigenvalues are shown in Fig. 7. The eigenvalues estimated by DMD at early time (0.002 to 0.004 μs) have a large imaginary component except for the rightmost value. As time progresses, the imaginary part of the eigenvalues decreases, and the real part moves rightward. This demonstrates a feature of the DMD method: Early in time there are many modes present in the solution, and the fast-decaying ones are governing the solution behavior early in time. As time goes on, only the slowly decaying modes are present, and DMD finds these later in time.

Fig. 7. Alpha eigenvalues for the 70-group test problem estimated by DMD over three different time intervals.

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.

Fig. 8. Neutron population and spectra in the outer reflector, fuel, and moderator averaged over the three time intervals. The time intervals are denoted by black lines in (a), and the fundamental $k$-eigenvalue spectra are shown in (b) through (d).

In Fig. 8, we show the neutron spectrum at several points in space. The spectra shown are computed using time steps 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 $k$-eigenvalue 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.

Fig. 9. Spatial distribution of neutrons for the fundamental mode of the $k$-eigenvalue problem, and the eigenvector for the rightmost alpha eigenvalue as estimated by DMD over different time intervals. Note that the alpha eigenvectors are not positive, so we plot the absolute value. In this figure thermal neutrons have energy below 5 eV, fast neutrons are above 0.5 MeV, and epithermal neutrons are in between.

VI. DISCUSSION

Dynamic mode decomposition allows for the approximation of the eigenvalues present in a time-dependent transport system from the solution at different times without a separate eigenvalue solve. 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.

Acknowledgments

The author would like to thank B. D. Lansrud-Lopez, C. D. Ahrens, and R. D. Baker for helpful discussions during the development of this work. Also, thanks are in order to S. R. Bolding for sharing some Python code for cross-section processing and infinite medium solutions.