<figure><img src="image_1.png" /></figure>

HAL
open science

# Analysis of Alpha Modes in Multigroup Diffusion

R. Sanchez, D. Tomatis, I. Zmijarevic, K. Dugan

► To cite this version:

R. Sanchez, D. Tomatis, I. Zmijarevic, K. Dugan. Analysis of Alpha Modes in Multigroup Diffusion. Nuclear Engineering and Technology, 2017, Special Issue on International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering 2017 (M&C 2017), 49 (6), pp.1259-1268. cea-02434003

HAL Id: cea-02434003

https://cea.hal.science/cea-02434003

Submitted on 9 Jan 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Analysis of Alpha Modes in Multigroup Diffusion

Richard Sanchez¹,ª Daniele Tomatis,ᵇ Igor Zmijarevicᵇ and Kevin Duganᵇ

ª Department of Nuclear Engineering, Seoul National University

1 Gwanak-ro, Gwanak-gu, Seoul, 151-744, Korea

ᵇ DEN, Service d'études des réacteurs et de mathématiques appliquées, CEA, Université Paris-Saclay

F-91191, Gif-sur-Yvette, France

sanchezr@snu.ac.kr, daniele.tomatis@cea.fr, igor.zmijarevic@cea.fr, kevin.dugan@cea.fr

**Abstract** - The alpha eigenvalue problem in multigroup neutron diffusion is studied with particular attention to the theoretical analysis of the model. Contrary to previous literature results, the existence of eigenvalue and eigenflux clustering is here investigated without the simplification of a unique fissile isotope or a single emission spectrum. A discussion about the negative decay constants of the neutron precursors concentrations as potential eigenvalues is provided. An in-hour equation is derived by a perturbation approach recurring to the steady state adjoint and direct eigenvalue problems of the effective multiplication factor and is used to suggest proper detection criteria of flux clustering. In spite of to prior work, the in-hour equation results for a necessary and sufficient condition for the existence of the eigenvalue-eigenvector pair. A simplified asymptotic analysis is used to predict bands of accumulation of eigenvalues close to the negative decay constants of the precursors concentrations. The resolution of the problem in one-dimensional heterogeneous problems shows numerical evidence of the predicted clustering occurrences and also confirms previous theoretical analysis and numerical results.

keywords: neutron kinetics, alpha modes, clustering.

# I. INTRODUCTION

The kinetic neutron equation with precursors is an efficient tool for the analysis of the neutron time evolution and, therefore, finds multiple applications in reactor control and the study of accident scenarios. The time-dependent flux depends on the initial conditions for neutrons and precursors, the boundary condition for neutrons, the presence of external sources and the time-dependent changes of the cross sections, such as those produced by rod motion and changes of Boron concentration or to small stochastic perturbation induced by the coolant flow and, in the long time, by radioactive decay and nuclide depletion and creation from fission.

However, if the cross sections and boundary conditions remain constant in time and the sources vanish, then the state of the system tends asymptotically in time to an exponential behavior, which is independent of the earlier changes of the system and, specially, of the initial conditions.[1]

The exponential behaviors that a given system can adopt are the solutions of an eigenvalue equation and can be used to describe the fast evolution of the system as well as to characterize the reactivity of the system. The solutions of this equation are known as time-dependent modes or, more simply, alpha modes, where “alpha” refers to the most frequently adopted symbol for these eigenvalues.

Alpha modes have been applied to formally derive different forms of the well-known in-hour equation, to obtain solutions of the kinetic equations by expansion techniques,[2] to develop numerical solution methods[3] and also as weighting fluxes to homogenize the kinetic equation.[4] These applications as well as the theoretical interest of alpha modes have

been the object of intense study. However, with the exception of a few scattered mathematical results,[5, 6, 7] a detailed description of these modes have not yet been given for the diffusion equation. The purpose of this paper is to give a detailed analysis of the alpha modes for a slab geometry using multigroup diffusion theory. We give ample numerical evidence for all the modes predicted from mathematical analysis as well as from physical arguments.

General equations for the alpha eigenvalue problem are discussed in Sec. II, including a perturbation expression in terms of reactivity. In Sec. III we summarize our results and observations for the multigroup one-dimensional slab. Some numerical results are illustrated in Sec. IV, while in Sec. V we compare with Asahi's analytical results for the one-group problem.[6] Conclusions follow in Sec. VI.

For simplicity we shall use a notation based on a continuous formulation and let to the reader the change to a fully discretized operator which is used in the numerical application.

# II. THE ALPHA EIGENVALUE EQUATIONS

Our starting point are the time-dependent kinetic diffusion equations coupled to the precursors equations in a heterogeneous domain $\mathcal{D}$:

$$ \left\{ \left( \frac{1}{v} \partial_t + \mathcal{B}_{pr} \right) \psi = \sum_p \chi_p \lambda_p C_p, \right. \quad (1a) $$

$$ (\partial_t + \lambda_p) C_p = \mathcal{F}_p \psi. \quad (1b) $$

Here $\psi(x, t)$ is the scalar flux, $x = (\mathbf{r}, E)$ stands for the phase space variables, $C_p(\mathbf{r}, t)$ is the concentration for precursor $p$ and

$$ \mathcal{B}_{pr} = \mathcal{L} - \mathcal{P}_{pr} $$

<footnote>*Corresponding author</footnote>
is the diffusion kinetic operator with prompt neutron production

$$ \mathcal{L} = -\nabla \cdot D\nabla + \Sigma - \mathcal{H}, \quad (2) $$

where $D(x, t)$ is the diffusion coefficient, $\Sigma(x, t)$ is the total macroscopic cross section, $\mathcal{H}$ stands for the scattering operator and $\mathcal{P}_{pr} = \sum_i \chi_i(1 - \beta_i)\mathcal{F}_i$ is the prompt fission operator with the sum in $i$ over the fissile isotopes.

Also in Eqs. (1b) and (2), $\mathcal{F}_p = \sum_i \beta_{i \to p} \mathcal{F}_i$ and $\mathcal{F}_i$ are the production operators for precursor $p$ and neutrons from isotope $i$, respectively. The latter is defined as $\mathcal{F}_i\psi = \int v\Sigma_{f,i}\psi dE$. We note $N_p$ and $N_i$ the number of precursors and fissile isotopes, respectively, and take the prompt and delayed fission spectra, $\chi_i$ and $\chi_p$, to be normalized so that $\int \chi_i dE = 1$.

By introducing the reduced scalar product over the energy $E$,

$$ \langle f, g \rangle = \int \bar{f} g dE, \quad (3) $$

where $\bar{f}$ denotes the complex conjugate of $f$, the two latter formulas can be written as $\mathcal{F}_i\psi = \langle v\Sigma_{f,i}, \psi \rangle$ and $\langle 1, \chi \rangle = 1$. Finally, $\lambda_p$ and $\beta_{i \to p}$ are the decay constant for precursor $p$ and its fission yield from isotope $i$, respectively, and $\beta_i = \sum_p \beta_{i \to p}$ is the total precursor yield for isotope $i$. We shall follow the usual convention of ordering the $\lambda_p$ in increasing values.

The kinetic equations in (1) are supplemented with initial conditions for both $\psi$ and $C_p$ and with boundary conditions for $\psi$. The cross sections, external source and boundary conditions may change in time. But at any time $t = t_*$ one might consider a new kinetic problem with the initial and boundary conditions at time $t_*$ and such that for $t > t_*$ the cross sections and boundary conditions remain constant and the external source vanishes. This defines a kinetic problem with constant cross sections and boundary conditions and with no external source which will evolve in time according solely to its initial conditions at $t_*$. The kinetic equations in (1) are supplemented with initial conditions for both $\psi$ and $C_p$ and with boundary conditions for $\psi$. The cross sections, external source and boundary conditions may change in time. But at any time $t = t_*$ one might consider a new kinetic problem with the initial and boundary conditions at time $t_*$ and such that for $t > t_*$ the cross sections and boundary conditions remain constant and the external source vanishes. This defines a kinetic problem with constant cross sections and boundary conditions and with no external source which will evolve in time according solely to its initial conditions at $t_*$. The alpha eigenvalue equations for this $t_*$-kinetic problem are obtained by introducing an exponential time behavior in (1),

$$ \psi(x, t) \sim e^{\alpha t} \psi_{\alpha}(x), \quad C_p(\mathbf{r}, t) \sim e^{\alpha t} C_{p,\alpha}(\mathbf{r}), $$

to obtain:

$$ \left\{ \begin{array}{l} \left(\frac{\alpha}{v} + \mathcal{B}_{pr}\right) \psi_{\alpha} = \sum_{p} \chi_{p} \lambda_{p} C_{p,\alpha}, \quad (4a) \\ (\alpha + \lambda_{p}) C_{p,\alpha} = \mathcal{F}_{p} \psi_{\alpha}. \quad (4b) \end{array} \right. $$

Here $\psi_{\alpha}(x)$ and $C_{p,\alpha}(\mathbf{r})$ are respectively the scalar flux and the concentration for precursor $p$ of an $\alpha$-mode, and all other quantities and operators are like those defined for the time-dependent kinetic equation but with the cross sections evaluated at time $t_*$.

More generally, by introducing the vector function $\vec{\Psi}_{\alpha} = \{\psi_{\alpha}, \vec{C}_{\alpha}\}$, where $\vec{C}_{\alpha}$ has components $\{C_{p,\alpha}, \alpha = 1, N_p\}$, the set of alpha equations can be written in the form of a classical linear eigenvalue problem:

$$ \mathcal{M}\vec{\Psi}_{\alpha} = \alpha\vec{\Psi}_{\alpha}, \quad (5) $$

where the structure of matrix operator $\mathcal{M}$, not explicitly here, can be inferred from Eqs.(4).

When appropriate we shall use Hilbert spaces for a precise formulation of results. We note $\mathcal{H}_{\chi}$ the Hilbert space of complex-valued functions defined over phase space $\chi$ and satisfying the boundary conditions of problem (4). This space is endowed with the scalar product

$$ (f, g) = \int_{\mathcal{X}} \bar{f}(x)g(x)dx, \quad (6) $$

where $dx = d\mathbf{r}dE$ is the volume element in $\chi$. We also note $\mathcal{H}_{\mathcal{D}}$ the Hilbert space of complex-valued functions defined over the geometrical domain $\mathcal{D}$ with scalar product

$$ (f, g)_{\mathcal{D}} = \int_{\mathcal{D}} \bar{f}(\mathbf{r})g(\mathbf{r})d\mathbf{r}. $$

Thus, we consider problem (5) in the Hilbert space of complex-valued functions $\mathcal{H}_{\chi} \times (\mathcal{H}_{\mathcal{D}})^{N_p}$ with scalar product

$$ ((\vec{\Psi}, \vec{\Psi}')) = (\psi, \psi') + \sum_{p} (C_{p}, C'_{p})_{\mathcal{D}}. $$

An useful form of the alpha-equations can be obtained by solving Eq. (4b) for the precursor concentrations:

$$ C_{p,\alpha} = \frac{\mathcal{F}_p \psi_{\alpha}}{\alpha + \lambda_p}, \quad \alpha \in \mathbb{C} \setminus \Lambda, \quad (7) $$

where $\Lambda = \{-\lambda_p, p = 1, N_p\}$. Next, replacing this result in Eq. (4a) yields an expression for the alpha-eigenvalue problem in terms of only the flux:

$$ \left(\frac{\alpha}{v} + \bar{\mathcal{B}}\right) \psi_{\alpha} = -\mathcal{K}_{\alpha} \psi_{\alpha}, \quad \alpha \in \mathbb{C} \setminus \Lambda, \quad (8) $$

$$ \text{where } \bar{\mathcal{B}} = \mathcal{L} - \bar{\mathcal{P}}, \qquad \bar{\mathcal{P}} = \sum_{i} \bar{\chi}_{i} \mathcal{F}_{i} \quad (9) $$

is the production operator with the steady-state fission spectrum

$$ \bar{\chi}_i = (1 - \beta_i)\chi_i + \sum_{p} \beta_{i \to p} \chi_p \quad (10) $$

and

$$ \mathcal{K}_{\alpha} = \sum_{p} \frac{\alpha}{\alpha + \lambda_{p}} \chi_{p} \mathcal{F}_{p} \quad (11) $$

comes from the delayed contribution to production. We consider this problem in $\mathcal{H}_{\chi}$.

We note that Eq. (8) is a non linear eigenvalue problem for $(\alpha, \psi_{\alpha})$. Moreover, as we discuss next, except for very pathological cases, $\Lambda$ is in the resolvent set of the linear operator in Eq. (5) and therefore Eqs. (4) and (8) are equivalent.

Let us analyze in which conditions $\alpha = -\lambda_p$ can be in the spectrum of Eqs. (4). If this is true, then one must have from Eq. (4b)

$$ \mathcal{F}_p \psi_{\alpha} = 0. \quad (12) $$

We note that one can write $\mathcal{F}_p = h\pi_p$ where $h(\mathbf{r}) = (\mathcal{F}_p 1)(\mathbf{r})$ and $\pi_p = h^{-1}\mathcal{F}_p$ is a projector $\pi_p : \mathcal{H}_{\chi} \to \mathcal{H}_{\mathcal{D}}$. Equivalently, constraint (12) can be written as

$$ \psi_{\alpha} \in \ker(\pi_{p}), $$
where $ker(\pi_p)$ is the null space of $\pi_p$. For the other precursors one obtains

$$C_{p',\alpha} = \frac{1}{\lambda_{p'} - \lambda_p} \mathcal{F}_{p'} \psi_{\alpha}, \quad p' \neq p, \qquad (13)$$

which yields a source equation for the flux:

$$\hat{\mathcal{B}}\psi_{\alpha} = \chi_{p}\lambda_{p}C_{p,\alpha}, \qquad (14)$$

where

$$\hat{\mathcal{B}} = \bar{\mathcal{B}} - \frac{\lambda_p}{v} - \sum_{p' \neq p} \frac{\lambda_p}{\lambda_{p'} - \lambda_p} \chi_{p'} \mathcal{F}_{p'}.$$

Hence, a solution of Eq. (14) with the constraint in (12) is a pair $(\psi_{\alpha}, C_{p,\alpha} = \hat{\mathcal{B}}\psi_{\alpha}/(\chi_{p}\lambda_{p}))^2$ iff<fnref n="3" /> (i) $0 \neq \psi_{\alpha} \in ker(\pi_p)$, (ii) $\hat{\mathcal{B}}\psi_{\alpha} \in \mathcal{E}_p$, where we have defined $\mathcal{E}_p = \{f(x) = \chi_p(E)h(\mathbf{r}); h(\mathbf{r}) \in \mathcal{H}_D\}$. In other words $-\lambda_p$ is an eigenvalue iff $\hat{\mathcal{B}}(ker(\pi_p)/\{0\}) \cap \mathcal{E}_p \neq \emptyset$. Although this iff condition does not apparently require operator $\hat{\mathcal{B}}$ to have an inverse, an equivalent formulation can be written as $\hat{\mathcal{B}}^{-1}\mathcal{E}_p \cap (ker(\pi_p)/\{0\}) \neq \emptyset$. Thus, a solution might exists if operator $\hat{\mathcal{B}}$ is degenerated or if it is invertible over $\mathcal{E}_p$; the former condition seems unrealistic but the second looks plausible because $v$ is large and $\Sigma$ contains scattering, capture and fission and therefore dominates $\mathcal{H}$. However, in both cases the solution of equation (14) must satisfy the orthogonality condition (12), which most likely will not be feasible. Therefore, with the exception of unrealistic cross sections, the existence of an alpha eigenvalue equal to minus a precursor decay constant can be excluded.

Clearly the precedent conclusion applies to the case of the one-group diffusion equation. But we give here a counter example for which $-\lambda_p$ is an alpha eigenvalue. Consider the case when the diffusion coefficient does not depend on the energy, $D = D(\mathbf{r})$, and take $\psi_{\alpha}(x) = f(E)h(\mathbf{r})$ with $\pi_p f = 0$ and $h(\mathbf{r})$ an eigenfunction of the Laplacian equation $-\nabla \cdot D\nabla h = \gamma h$ which satisfies the boundary conditions of the alpha problem. Then $\hat{\mathcal{B}}\psi_{\alpha} = h\hat{\mathcal{B}}\gamma f$, where $\hat{\mathcal{B}}\gamma$ equals $\hat{\mathcal{B}}$ with the replacement $-\nabla \cdot D\nabla \rightarrow \gamma$. Therefore, by setting $\chi_p = \hat{\mathcal{B}}\gamma f/c$, where $c(\mathbf{r}) = \int \hat{\mathcal{B}}\gamma f dE$, the pair $(\psi_{\alpha}, C_{p,\alpha} = (c/\lambda_p)h)$ is a solution of the alpha eigenvalue problem for $\alpha = -\lambda_p$. A problem with the definition of $\chi_p$ is that it depends on space, but this can be avoided by considering a homogeneous problem. However, even then, one would have to “adjust” the cross sections, perhaps with nonphysical values, so as to ensure that $\chi_p$ remains positive.

# 1. Perturbation formulation

Reactivity is a measure of the disequilibrium of a system and there are several ways to introduce a notion of reactivity in neutron kinetics[8, 9], in particular that derived with the help of the adjoint $\alpha$ equation.[10] In this work we define a generalized reactivity as $\rho_{\lambda} = 1 - 1/\lambda$, where $\lambda$ is any eigenvalue of the quasi-static eigenvalue equation obtained by setting $\alpha = 0$ in Eqs. (8):

$$\bar{\mathcal{B}}_{\lambda}\varphi_{\lambda} = 0, \qquad (15)$$

where

$$\bar{\mathcal{B}}_{\lambda} = \mathcal{L} - (1/\lambda)\bar{\mathcal{P}}. \qquad (16)$$

Here $\lambda \in \sigma(\mathcal{L}^{-1}\bar{\mathcal{P}})$, where $\sigma$ denotes the spectrum of the operator, is an eigenvalue of the quasi-static diffusion operator and $\bar{\mathcal{P}}$ is the steady-state production operator in Eq. (9). The boundary condition for Eq. (15) is that of the parent alpha eigenvalue equations. We note $E_{\lambda}$ and $N_{\lambda}$, respectively, the subspace of eigenfunctions of eigenvalue $\lambda$ and its dimension. We recall that the eigenvalue $\lambda = k$ with maximum absolute value in $\sigma(\mathcal{L}^{-1}\bar{\mathcal{P}})$ is real, non degenerate and has a positive eigenfunction $\varphi_k$ (a physical flux).

In order to analyze the solutions of the alpha-eigenvalue equation it is advantageous to write Eq. (8) in terms of the quasi-static operator $\bar{\mathcal{B}}_{\lambda}$ in (16) and cast the alpha-eigenvalue equation as

$$\bar{\mathcal{B}}_{\lambda}\psi_{\alpha} = \mathcal{K}_{\lambda,\alpha}\psi_{\alpha}, \qquad (17)$$

where $\lambda \in \sigma_p(\mathcal{L}^{-1}\bar{\mathcal{P}})$ and

$$\mathcal{K}_{\lambda,\alpha} = \rho_{\lambda} \bar{\mathcal{P}} - \left(\frac{\alpha}{v} + \mathcal{K}_{\alpha}\right). \qquad (18)$$

The appeal of formulation (17) is that (i) it offers a perturbation expression for the alpha eigenvalue problem in terms of the better known steady-state eigenvalue problem, and (ii) the operator $\bar{\mathcal{B}}_{\lambda}$ on the left of the equal sign is singular and therefore the “source” term $\mathcal{K}_{\lambda,\alpha}\psi_{\alpha}$ must satisfy a solvability condition.

According to the Riesz-Schauder theory,[11] the solvability condition requires the source term to be in the orthogonal complement $E^{\perp}$ of the eigenspace $E_{\lambda}^{\perp}$ associated to the eigenvalue $\bar{\lambda} \in \sigma((\bar{\mathcal{P}}\mathcal{L}^{-1})^{\dagger})$ of the adjoint operator $\bar{\mathcal{B}}_{\lambda}^{\dagger}$.<fnref n="4" /> More explicitly,

$$(\varphi^{\dagger}, \mathcal{K}_{\lambda,\alpha}\psi_{\alpha}) = 0 \quad \forall \varphi^{\dagger} \in E_{\lambda}^{\perp}. \qquad (19)$$

Thus, the solvability condition results in a system of $N_{\lambda}$ one-point-like “kinetic” alpha equations. For each function $\varphi^{\dagger} \in E_{\lambda}^{\perp}$ we have an inhour-type equation:

$$\rho_{\lambda} = \alpha \left( \Lambda + \sum_{p} \frac{\Gamma_{p}}{\alpha + \lambda_{p}} \right), \qquad (20)$$

where

$$\Lambda = \frac{(\varphi^{\dagger}, \frac{1}{v}\psi_{\alpha})}{(\varphi^{\dagger}, \bar{\mathcal{P}}\psi_{\alpha})}, \quad \Gamma_{p} = \frac{(\varphi^{\dagger}, \chi_{p}\mathcal{F}_{p}\psi_{\alpha})}{(\varphi^{\dagger}, \bar{\mathcal{P}}\psi_{\alpha})}. \qquad (21)$$

As an aside, we note that in Eqs. (21) one can write $\Gamma_p = \sum_i \beta_{i \to p} \gamma_{i,p}$ with the effectiveness factors[13] $\gamma_{i,p} = \int \langle \varphi_k^{\dagger}, \chi_p \rangle \langle \nu \Sigma_{f,i}, \psi_{\alpha} \rangle d\mathbf{r} / (\varphi_k^{\dagger}, \bar{\mathcal{P}}\psi_{\alpha})$ and $(\varphi_k^{\dagger}, \bar{\mathcal{P}}\psi_{\alpha}) = \sum_i \int \langle \varphi_k^{\dagger}, \bar{\chi}_i \rangle \langle \nu \Sigma_{f,i}, \psi_{\alpha} \rangle d\mathbf{r}$, recovering thus the formula derived for an infinite homogeneous medium in Ref. [14]. Also, because of the projection by $\varphi^{\dagger}$ the original eigenvalue $\alpha$ is dispersed among the $N_p + 1$ roots of Eq. (20).

<footnote n="2">²For simplicity, we assume that $\chi_p$ is strictly positive.</footnote>

<footnote n="3">³if and only if</footnote>

<footnote n="4">⁴It is understood that the adjoint operator is provided with appropriate “adjoint” boundary conditions [12].</footnote>
Under the solvability constraint (19), the general solution of (17) is an arbitrary multiple of

$$ \psi_{\alpha} = \varphi_{\lambda} + \psi^{\perp}, \qquad (22) $$

where $\varphi_{\lambda} = \sum_{n} a_{n} \varphi_{n}$ is written using an arbitrary basis $\{\varphi_{n}, n = 1, N_{\lambda}\}$ in $E_{\lambda}$ and $\psi^{\perp} \in E^{\perp}$. The latter condition entails

$$ a_{n} = (\varphi_{n}^{\dagger}, \psi_{\alpha}), \qquad (23) $$

where the $\{\varphi_{n}^{\dagger}, n = 1, N_{\lambda}\}$ is the dual basis in $E_{\lambda}^{\dagger}$, i.e., such that $(\varphi_{n}^{\dagger}, \varphi_{n'}) = \delta_{nn'}$.

We end this section with a comment regarding the effect of spatial symmetries on the structure of the in-hour equation. Consider a problem with a spatial involution, that is a symmetry $s$ such that $s^2 = 1$. Then $\psi_{\alpha}$ obeys the relation $s^2\psi_{\alpha} = \psi_{\alpha}$ and, since operator $s^2$ has eigenvalues 1 and -1, one can write $\psi_{\alpha}$ as the sum of an even (+) and an odd (-) components $\psi_{\alpha, \pm} = (s \pm 1)\psi_{\alpha}$ with the property that $\int d\mathbf{r}\psi_{\alpha,+}(\mathbf{r})\psi_{\alpha,-}(\mathbf{r})$ vanishes. Because both the steady-state transport equation and the alpha eigenvalue equation necessarily share the same geometry one concludes that the in-hour equation does not exist when $\varphi^{\dagger}$ and $\psi_{\alpha}$ have opposite symmetries.

### III. SUMMARY OF RESULTS

The spectrum of the alpha mode equations is conditioned by the nature of these equations and by the degree of complexity accepted in the delayed neutron term. In the continuous case the eigenmodes, $\tilde{\Psi}_{\alpha}(x)$, belong to a functional space of infinite dimensions, while in the discretized case they become vectors in a finite-dimensional space. On the other hand, the degree of complexity of the fission delayed contribution depends on the number of fissile isotopes $N_{i}$ and the number of precursors $N_{p}$. The operator in Eq. (11) can be written as a sum of projectors with $\mathbf{r}$-dependent coefficients:

$$ \mathcal{K}_{\alpha} = \sum_{p} \frac{\alpha}{\alpha + \lambda_{p}} a_{p} \pi_{p}, \qquad (24) $$

where $a_{p}(\mathbf{r}) = \mathcal{F}_{p}\chi_{p}$ and $\pi_{p} = (1/a_{p})\chi_{p}\mathcal{F}_{p}$ is a projector over the subspace generated by the $\{\chi_{p}, p = 1, N_{p}\}$. Therefore, for $\alpha \notin \{0 \cup \Lambda\}$, $\mathcal{K}_{\alpha}$ is a finite linear combination of projectors with finite range equal to the span of the functions $\{\chi_{p}, p = 1, N_{p}\}$. Moreover the null space of $\mathcal{K}_{\alpha}$ contains the set of all the functions that are orthogonal to the finite set $\{\nu\Sigma_{f,i}, i = 1, N_{i}\}$.

Simplifications of the fission delayed contributions are typically introduced to facilitate the theoretical analysis of the alpha spectrum. A simplifying assumption can be applied to the delayed fission spectra, to the fission production or to both. The first consists of replacing all the delayed neutron spectra with an averaged one, $\chi_{p} \sim \chi$, which leads to the simplified form $\mathcal{K}_{\alpha} = \chi\mathcal{F}_{\alpha}$, where $\mathcal{F}_{\alpha} = \sum_{p}[\alpha/(\alpha + \lambda_{p})]\mathcal{F}_{p}$. In the second one replaces the fission production for all fissile isotopes with a single one, $\nu\Sigma_{f,i} \sim w_{i}(\mathbf{r})\nu\Sigma_{f}$, where by convenience we assume $\sum_{i} w_{i}(\mathbf{r}) = 1$. This approximation yields the simplified form $\mathcal{K}_{\alpha} = \chi_{\alpha}\mathcal{F}$, where now $\chi_{\alpha} = \sum_{p}[\alpha/(\alpha + \lambda_{p})]\beta_{p}\chi_{p}$ with $\beta_{p}(\mathbf{r}) = \sum_{i} w_{i}(\mathbf{r})\beta_{i,p}$. Note that these approximations are automatically satisfied if one assumes one single precursor

($\chi = \chi_{1}$) or one single fissile isotope ($w_{1} = 1$ and $\beta_{p} = \beta_{1 \to p}$). Finally, when both approximations are used one has the simplified operator

$$ \mathcal{K}_{\alpha} = f_{\alpha}(\mathbf{r})\chi\mathcal{F}, \qquad (25) $$

where

$$ f_{\alpha}(\mathbf{r}) = \sum_{p} \frac{\alpha}{\alpha + \lambda_{p}} \beta_{p}(\mathbf{r}). \qquad (26) $$

Note that for a single fissile isotope $\beta_{p}(\mathbf{r}) = \beta_{p}$, while for a homogeneous medium $\beta_{p} = \sum_{i} w_{i}\beta_{i \to p}$. For both cases $f_{\alpha}(\mathbf{r})$ is independent of $\mathbf{r}$ and is an analytical function of $\alpha$ in $\mathbb{C} \setminus \Lambda$.

There is a sizable number of theoretical studies of the alpha spectrum in transport theory but, to our knowledge, there are only three communications pertaining to diffusion theory. Porsching[5] used algebraic techniques to analyze the spectrum of the alpha mode discretized diffusion equations in a homogeneous medium with a single fissile isotope. However his equations were based in a synthetic expression for the time-dependent flux and, in regard to the present standard alpha eigenvalue problem, Porsching's results apply only to the one-group diffusion equation. A generalization to the multigroup diffusion equation was carried out by Devooght[7] for the fully discretized diffusion equation in a heterogeneous medium with no upscattering and with a single fissile isotope. To simplify the delayed fission term, Devooght assumes that all the fissile neutron spectra are identical, i.e., not only that $\chi_{p} \sim \chi$ but also $\bar{\chi}_{1} \sim \chi$, which lead him to lump prompt and delayed fission in a single term and work with the function $f(\alpha) = 1 - f_{\alpha}$, which is also analytical in $\mathbb{C} \setminus \Lambda$. This function was used in the general study of the spectrum. The analysis of the fast spectrum, based on a perturbation result for the case of degenerated operators discussed in Kato,[15] was carried out for the case with no precursors, while a technique early introduced by Wings[1] was applied to the investigation of the degeneracy of the delayed spectrum. The result of the analysis is that the spectrum consists of $N_{g} + N_{p}$ clusters of $N_{r}$ eigenvalues each, where $N_{g}$ and $N_{r}$ are the number of groups and the number of spatial components, respectively. The delayed spectrum has $N_{p}$ clusters of real eigenvalues with the $p$-th cluster confined to the open interval$^5$

$$ I_{p} = (-\lambda_{p}, -\lambda_{p-1}) \qquad (27) $$

with the convention that $-\lambda_{0} = \infty$. The only theoretical result for the one-group continuous diffusion equation in one-dimensional heterogeneous slab geometry was given by Asahi[6] a few years earlier than Devooght's paper but seemingly unknown to the latter. Asahi demonstrated that for an arbitrary geometry the eigenvalues were real and bounded and that the $-\lambda_{p}$ values were not in the spectrum. He gave also a more detailed analysis for the case of the one-dimensional

<footnote>$^5$In reality, the result obtained by Devooght is that the $-\lambda_{p}$ also belongs to the $p$-th clusters and that, therefore, the precursor clusters contain $N_{r} + 1$ eigenvalues each. However, he recognized that this pathological result occurred because he had artificially replaced the precursor concentrations $C_{p}(\mathbf{r})$ with the energy dependent functions $C_{\uparrow}(\mathbf{r}, E) = \chi_{p}(E)C_{p}(\mathbf{r})$ resulting in the associated eigenvalue problem having $N_{r}N_{p}(N_{g} - 1)$ more eigenvalues than the original one</footnote>
slab, classifying the eigenfunctions by the number of nodes and showing that there are a countable infinite number of eigenvalues in each interval $I_p$ which accumulate at $-\lambda_p$ as well as in the open interval

$$I_{N_p}^L = (-\infty, -\lambda_{N_p}). \qquad (28)$$

The rightmost eigenvalue in each interval has an eigenfunction with positive flux (zero nodes) and the number of nodes of the eigenfunctions fluxes increases monotonously by one for the other eigenvalues in the interval. Finally, as shown in Appendix A, an analytical solution can be obtained for the one-group homogeneous case which confirms the structure of the prompt and delayed spectrum as well as the emergence of flux clustering. This result leads to a new derivation of Asahi's results for the one-group homogeneous finite slab (see Appendix B).

In the second part of this section we give general observations in the nature of the spectrum of the alpha eigenvalue problem derived from robust arguments and from numerical exploration of multigroup diffusion in a slab geometry. Our results match the known theoretical predictions in their restricted domain of validity.[5, 6, 7] Our discretization was done for the one-dimensional piecewise heterogeneous slab geometry by writing the linear alpha eigenvalue equations (5) in its multi-group form and by introducing a numerical discretization for the spatial dependence of both $\psi$ and $\bar{C}$. We observe that the number of degrees of freedom (DOF) of the discretized system of equations is

$$N_{DOF} = N_r(N_g + N_p). \qquad (29)$$

A natural ordering is build in the kinetic equations. This ordering stems from the observation that most of the neutrons resulting from fission are prompt and that, therefore, the contribution of delayed neutrons to the overall neutron production in the kinetic equations can be considered as small. In addition, the frequencies of emission of delayed neutrons are also much smaller than the neutron collision frequency. These two conditions can be expressed as

$$\beta_p \ll 1 \text{ and } \lambda_p \ll v\Sigma.$$

In the following we shall assume that all cross sections are $O(1)$ and we shall define $\epsilon_\beta = O(\beta)$ and $\epsilon_v = O(1/v)$. Physical data tell us that $\epsilon_\beta \sim 10^{-3} - 10^{-4}$ and $\epsilon_v \sim 10^{-5} - 10^{-9}$.

This ordering has direct implications in the distribution of the alpha spectrum, which can be sorted out into two subsets: delayed and prompt eigenvalues.

## 1. Delayed spectrum

These modes appear in the range $|\alpha| \sim O(\lambda_p)$, in $N_p$ clusters of $N_r$ eigenvalues each, for a total of $N_r N_p$ eigenvalues (see Fig. 1). The eigenvalues in the p-th cluster gather on the right of $-\lambda_p$:

$$-\lambda_p < \alpha_{p,N_r} < \dots < \alpha_{p,2} < \alpha_p < -\lambda_{p-1},$$

where $\alpha_{p,i}$ is the i-th eigenvalue of the cluster in the interval $I_p$. Also, we have written $\alpha_{p,1}$ simply as $\alpha_p$ to indicate the

special nature of the dominant mode in the cluster. While the other eigenvalues in the cluster are very close to $-\lambda_p$ and their values do not change much with the reactivity, the value of the dominant eigenvalue $\alpha_p$ increases with the reactivity towards $-\lambda_{p-1}$. The delayed eigenvalues are simple and real negative with the exception of $\alpha_1$ which becomes positive for positive reactivity. This latter eigenvalue admits a positive eigenvector (both fluxes and concentrations are positive). The remaining dominant eigenvalues ($\alpha_p, p > 1$) have positive fluxes but the concentrations change of sign.

<figure><img src="image_3.png" /><figcaption>Fig. 1. Clusterings of eigenvalues approaching the negative values of the precursor's decay constants.</figcaption></figure>

The p-th cluster can be analyzed using the expression

$$\alpha + \lambda_p = \epsilon_p \qquad (30)$$

and considering the eigenvalue problem for $|\epsilon_p| \to 0$. The analysis shows that there is a leading concentration for each of the $N_p$ clusters which because of the special variation degenerates into $N_r$ independent distributions. Hence, this part of the spectrum can also be named “precursor” spectrum. We have also made the conjecture that the eigenmodes concentrations in each cluster are linearly independent and this has been supported by all our calculations.

In order to analyze the clustering results we resort to an approach different from the ones adopted in the literature. We recognize that the behavior of the delayed spectrum is directly related to the changes undergone by the precursors concentrations and, therefore, replace Eq. (8) with an “equivalent” equation for the concentrations. We use Eq. (7) to write the delayed contribution in Eq. (8) in terms of the concentrations and by expliciting the $\psi_\alpha$ in the resulting equation we write

$$\psi_\alpha = -(\alpha/v + \bar{\mathcal{B}})^{-1}\alpha \sum_p \chi_p C_{p,\alpha}, \alpha \notin \sigma(-v\bar{\mathcal{B}}). \qquad (31)$$

Next operating on this equation with $(\alpha + \lambda_p)^{-1}\mathcal{F}_p$ yields a finite system of equations for the concentrations, which we choose to write:

$$(\alpha + \lambda_p)C_{p,\alpha}(\mathbf{r}) = -\alpha \sum_{p'} (\mathcal{T}_{pp',\alpha} C_{p',\alpha})(\mathbf{r}), \alpha \notin \sigma(-v\bar{\mathcal{B}}), \qquad (32)$$

where $\mathcal{T}_{pp',\alpha}$ is an integral operator over E and r with kernel

$$t_{pp',\alpha}(\mathbf{r}' \to \mathbf{r}) = \int dE \nu \Sigma_{f,p}(x) \int dE' \chi_{p'}(E') g_{\alpha}(x' \to x), \qquad (33)$$

where $\nu\Sigma_{f,p} = \sum_i \beta_{i \to p}$ is the fission production term for precursor $C_p$, $x = (E, \mathbf{r})$ and $g_\alpha(x' \to x)$ is the Green's function
of operator $(\alpha/v + \bar{B})^{-1}$. We note that Eq. (32) is not equivalent to Eq. (8) for the simple reason that the energy variable does not appear in the former equation. This can be stated more precisely for the discretized case where the number of DOF for Eq. (32) is $N_r N_p$ whereas Eq. (8) has $N_r N_g N_p$ DOF. This shows that, except for the one-group case, the two equations are not equivalent. However, if $\alpha$ is an eigenvalue of Eq. (32) with eigenvector $\vec{C}_\alpha(\mathbf{r})$ and $\alpha \notin \sigma(-v\mathcal{B})$, then $\alpha$ is also an eigenvalue of the full alpha problem in Eq. (8) with the eigenvector given by Eq. (31).

Thus, we can safely use Eq. (32) to analyze the delayed alpha spectrum. We introduce the assumption in Eq. (30) and use asymptotic analysis to study the behavior of the spectrum near $-\lambda_p$. A full analysis would consist of expanding the concentrations and the eigenvalue as a sum of increasing powers of $\epsilon_p$, but here we use a rough approximation and keep only the leading terms in the equation. However, the discussion is awkward when working with Eq. (32) and we shall, instead, take a step back and work with Eqs. (7) and (8). We notice first that, consistently with Eq. (7), we have $C_{p,\alpha} = \mathcal{F}_p \psi_\alpha / \epsilon_p \sim \epsilon_p^{-1}$ and, for $p' \neq p$, $C_{p',\alpha} = \mathcal{F}_{p'} \psi_\alpha / (\lambda_{p'} - \lambda_p + \epsilon_p) \sim O(1)$.

Hence, we can neglect the contribution from the $p' \neq p$ precursors in the sum in the right-hand-side of Eq. (8). Following the derivation of Eq. (32), this implies that the system of equations in (32) uncouples and, in particular, yields a single equation for $C_{p,\alpha}$. To simplify our notation we write this equation as follows

$$ \epsilon_p C_{p,\alpha}(\mathbf{r}) = \lambda_p (\mathcal{T}_\alpha C_{p,\alpha})(\mathbf{r}), \alpha \notin \sigma(-v\bar{\mathcal{B}}) \quad (34) $$

with $\epsilon_p$ playing now the role of eigenvalue. Still, this last equation is difficult to analyze and we should simplify the operator by keeping only the leading contributions. Going back to Eq. (8), we observe that because $\mathcal{F}_p \sim \epsilon_\beta$ the order of the right-hand term in this equation is $\epsilon_p^{-1} = \epsilon_\beta / \epsilon_p$ instead of $\epsilon_p^{-1}$. So, in order to have $\epsilon \ll 1$, we have to assume $\epsilon_p \ll \epsilon_\beta$. Having made sure that the right-hand-side of Eq. (8) is $\sim \epsilon_p^{-1}$, we observe that on the left-hand-side $(\alpha/v + \Sigma - \mathcal{H} - \bar{\mathcal{B}})\psi \sim O(1)$ so we have $\alpha/v + \bar{\mathcal{B}} \sim -\nabla \cdot D(E, \mathbf{r})\nabla$. The spectrum of this operator consists of a countable infinite number of real positive eigenvalues growing towards $+\infty$; these eigenvalues have finite multiplicity and the corresponding eigenfunctions form a complete basis in $L_2$.[16] It follows that the inverse of this operator is compact and definite positive and its Green's function is of the form $g_\alpha(E, \mathbf{r}' \to \mathbf{r})$, where $E$ appears only as a parameter. This implies that operator $\mathcal{T}$ is also compact, and therefore it has real eigenvalues, but it is not necessarily definite positive. However, it is easily verified that for a homogeneous medium as well as in the one-group case $\mathcal{T}$ is definite positive and therefore its eigenvalues are positive and accumulate at zero. This proves that the eigenvalues that the eigenvalues $\alpha$ which converge towards $-\lambda_p$ do this from the right $(\alpha > -\lambda_p)$.

## 2. Prompt spectrum

The remaining part of the spectrum lies on the complex plane $\mathcal{R}(\alpha) < -\lambda_{N_p}$ and consists of simple eigenvalues that are real or complex conjugated in the range $|\alpha| \sim O(v\Sigma)$; the

exception is the first dominant eigenvalue $\alpha_{N_p+1}$ which is real and, with increasing reactivity tends to $-\lambda_{N_p}$. The eigenvectors have negligible concentration values and they are dominated by energy modes.

To analyze these eigenvalues we use again an asymptotic approach by assuming that $\alpha \sim O(v)$. Hence, one can neglect the contribution of the delayed neutrons and write the alpha eigenvalue equation by keeping only the leading terms as

$$ \left(\frac{\alpha}{v} + \mathcal{B}_{pr}\right)\psi_{\alpha} \sim 0. \quad (35) $$

This shows that the prompt spectrum is closely related to the spectrum of the diffusion equation with prompt fission, $\alpha \in \sigma(-v\mathcal{B}_{pr})$.

## 3. Flux clustering

We use the expression flux clustering to refer to a set of distinct eigenvalues which have very close fluxes. Flux clustering for the alpha eigenvalue transport problem was discussed earlier by Gozani[17] and by Henry.[10] Gozani applied his analysis also to diffusion theory, but the most interesting result was Henri's analysis of the two-group P1 transport equation, which is closely related to diffusion; by using the analytical solutions of this problem, Henri was able to give a satisfactory argument, albeit non rigorous, in support of flux clustering.

The flux components of the $N_{p+1}$ dominant alpha modes can be shown to cluster about the dominant mode of the quasi-static equation $\varphi_k$. This includes the dominant modes in the delayed spectrum ($\alpha_p, p = 1, N_p$) plus the dominant mode of the prompt spectrum ($\alpha_{N_{p+1}}$). It can also be expected that non-dominant delayed modes can exhibit some form of flux clustering, like sharing the same number of spatial fluctuations in the slab case, and that this type of flux clustering might even apply to prompt modes. However, in our numerical calculations we have not observed yet any such clustering, except for the fact that for each energy group one might be able to identify $N_r$ modes that are associated, in the sense that they have a similar spatial behavior.

Our analysis of flux clustering is based on the perturbation approach discussed in Sec. 1, whose main results are the in-hour equation and the flux perturbation formula, Eqs. (20) and Eq. (22), respectively. For simplicity, hereafter we shall only consider the case when $\lambda$ is the dominant eigenvalue of the quasi-static diffusion operator, which we denote by $k$. This eigenvalue is simple, real positive and its eigenvector $\varphi_k$ is the so-called fundamental flux. Hence, the system of in-hour equations in (20) simplifies to a single equation which we write as

$$ \rho_k = \omega \left( \Lambda_\alpha + \sum_p \frac{\Gamma_{p,\alpha}}{\omega + \lambda_p} \right), \quad (36) $$

where, to stress the fact that the equation has $N_p + 1$ roots of which only one is the original alpha eigenvalue, we have replaced $\alpha$ with $\omega$ and we have used the lower index $\alpha$ to indicate that the coefficients $\Lambda$ and $\Gamma_p$ depend on the eigenvector $\psi_\alpha$.

The in-hour equation can have complex roots and we refrain here from considering this case. But, as shown in
Appendix C, if the coefficients $\Lambda_\alpha$ and the $\Gamma_{p,\alpha}$ in the in-hour equation have the same sign, which we take to be positive, then all the $N_p + 1$ roots are real. Therefore, since the associated steady-state adjoint flux $\varphi_k^\dagger$ is positive, we shall consider only the case when the alpha eigenvector $\psi$ is also positive. Note that this implies that the eigenvalue $\alpha$ is real. This is in fact the more clear-cut and best understood example of flux clustering.

A substantial argument in support of flux clustering can be made by looking for alpha modes with fluxes close to the fundamental static flux $\varphi_k$. We proceed by writing the flux perturbation formula, Eq. (22), in the form

$$ \psi_{\alpha} = (1 + \epsilon)\varphi_{k} + \psi^{\perp}, \quad (37) $$

where $|\epsilon| \ll 1$ and where we have adopted the normalizations $\|\varphi_k\| = \|\psi_\alpha\|$. Next, we use this expression to compute the coefficients in the in-hour equation Eq. (36) and observe that to leading term in $\epsilon$ the coefficients depend only on the static flux $\varphi_k$. The reason is that $\|\psi^\perp\| \sim \sqrt{\epsilon}$. As an aside, we point out that another reason for neglecting $\psi^\perp$ is that the coefficients depend on integrals of a product of a positive function times $\psi_\alpha$ and that the impact of $\psi^\perp$ in the value of the integrals is further reduced because, as opposed to $\varphi_k$, $\psi^\perp$ must necessarily change of sign.

Hence, to leading order in $\epsilon$, we might conclude that Eq. (36) is independent of $\psi_\alpha$ and also of small variations in $\epsilon$. Since this equation has $N_p + 1$ real roots we can safely assume that each one of these roots corresponds, to order $\epsilon$, to an alpha eigenvalue. However, there are two conundrums regarding the precedent conclusion. The first one is that if the in-hour equation were to depend only on $\varphi_k$, then this flux will be a solution of the equation, but the only case in which this might happen is for $k = 1$ which entails, of course, $\alpha = 0$. The second problem arises from the fact that the in-hour equation is only a necessary condition for $(\alpha, \psi_\alpha)$ to be an eigenpair of the alpha eigenvalue problem. Hence, the problem has to resolved with a simultaneous consideration of both the in-hour equation and the original alpha eigenvalue equation. For the present approach this involves showing that solutions of the form in Eq. (37) are truly solutions of the alpha eigenvalue problem. This is what, in a limited way, Henry did for the two-group P1 equation. For a more general approach may consist perhaps of using a direct perturbation analysis between the alpha eigenvalue equation and the equation satisfied by $\varphi_k$.

Finally, we observe that the leading eigenvalue $\alpha_1$ dominates the alpha spectrum and consequently it has a positive flux. It remains to show, however, that the fluxes of the other $N_p$ dominant eigenvalues $\alpha_p$, $p = 2, N_p$ are also positive, a fact that is confirmed by previous speculations and by all our calculations. In view of formula (37), this could perhaps be demonstrated by a dominance argument.

## IV. NUMERICAL RESULTS

We have used finite differences to discretize the multi-group diffusion alpha eigenvalue equation in 1D slab geometry as well as the eigenvalue equations for the quasi-static equation and its adjoint, and implemented three different criteria, all based in the $L_2$ norm $\|\psi_{ref} - \psi_\alpha\|_2 =$

$\sqrt{\sum_{g,r} (\psi_{ref} - \psi_\alpha)^2}_{g,r} < \epsilon_{clus}$, to explore flux clustering about alpha modes $(\psi_{ref} = \psi_{\alpha'}, \forall \alpha')$, about the modes of the quasi-static equation $(\psi_{ref} = \varphi_\lambda, \forall \lambda)$ and by using perturbation analysis $\|\psi_\perp\|_2, |1 - a| < \epsilon_{clus}$, where the latter formula has been obtained from Eqs. (22) and (23) for the case of a non degenerated eigenvalue ($N_\lambda = 1$). Calculations are performed in double precision by a Python module (v.2.7), using the numerical numpy library (v.1.11.0). The module solves the multi-group, ($S_N$), multi-region $\alpha$ and $\lambda$ eigenvalue problems with reflection, periodic and void boundary conditions; the full set of eigenpairs is computed.

We have analyzed a large number of problems for the infinite homogeneous case as well as for some finite cases with 1, 8 and 281 energy groups and for increasing values of $N_r$. Flux weighting was used to homogenize and collapse the cross sections for a PWR fuel assembly at different conditions (see Table I). Thus, the corresponding $\alpha$-sets are close. The calculations where done with 8 groups of precursors and 6 fissile isotopes ($N_p = 8$ and $N_i = 6$). We present here a sample of the results.

<table><caption>TABLE I. Physical conditions of the cross section preparation in the PWR fuel assembly.</caption><thead><tr><th>Case nb.</th><th>Tfuel (°C)</th><th>Tmod (°C)</th><th>Boron conc. (PPM)</th><th>ρ (pcm)</th></tr></thead><tbody><tr><td>1</td><td>286</td><td>286</td><td>0</td><td>+6880.6</td></tr><tr><td>2</td><td>286</td><td>286</td><td>400</td><td>+3405.6</td></tr><tr><td>3</td><td>286</td><td>286</td><td>600</td><td>+1690.3</td></tr><tr><td>4</td><td>286</td><td>286</td><td>800</td><td>-9.3</td></tr><tr><td>5</td><td>650</td><td>306</td><td>0</td><td>+4915.5</td></tr><tr><td>6</td><td>650</td><td>306</td><td>400</td><td>+1624.2</td></tr><tr><td>7</td><td>650</td><td>306</td><td>600</td><td>+0.1</td></tr><tr><td>8</td><td>650</td><td>306</td><td>800</td><td>-1611.3</td></tr><tr><td>9</td><td>1100</td><td>306</td><td>600</td><td>-1302.5</td></tr><tr><td>10</td><td>1800</td><td>306</td><td>600</td><td>-3037.1</td></tr></tbody></table>

For all cases in Table I, the neutron flux components of the $\psi_\alpha$ corresponding to the dominant eigenvalues $\alpha_p$ clustered within $\epsilon_{clus} = 10^{-2}$. The dominant eigenvalues for cases 3 and 4 are reported in Tables. II and III, respectively.

<table><caption>TABLE II. Comparisons of the fundamental and delayed reactor periods with 1, 8 and 281 energy groups for case 3 ($k_\infty = +1.017200$).</caption><thead><tr><th>α's</th><th colspan="2">rel. diff. (%)</th></tr><tr><th>281G</th><th>8G</th><th>1G</th></tr></thead><tbody><tr><td>+5.8029093620e+02</td><td>-25.514</td><td>+0.214</td></tr><tr><td>-1.2552000351e-02</td><td>-0.033</td><td>+0.002</td></tr><tr><td>-2.9640277111e-02</td><td>-0.235</td><td>-0.074</td></tr><tr><td>-4.3562900973e-02</td><td>-0.139</td><td>-0.005</td></tr><tr><td>-1.4048345489e-01</td><td>-0.282</td><td>-0.030</td></tr><tr><td>-3.2792305224e-01</td><td>-0.680</td><td>-0.118</td></tr><tr><td>-7.0546045788e-01</td><td>-0.379</td><td>-0.054</td></tr><tr><td>-1.7081304062e+00</td><td>-0.301</td><td>-0.041</td></tr><tr><td>-3.6271274472e+00</td><td>-0.139</td><td>-0.027</td></tr></tbody></table>
<table><caption>TABLE III. Comparisons of the fundamental and delayed reactor periods with 1, 8 and 281 energy groups for case 4 ($k_∞ = +0.999910$).</caption><thead><tr><th>α's</th><th colspan="2">rel. diff. (%)</th></tr><tr><th>281G</th><th>8G</th><th>1G</th></tr></thead><tbody><tr><td>-1.4397459518e-03</td><td>+4.003</td><td>+2.054</td></tr><tr><td>-1.4020587820e-02</td><td>+0.056</td><td>+0.171</td></tr><tr><td>-3.7787664954e-02</td><td>-0.011</td><td>-0.133</td></tr><tr><td>-8.2605497512e-02</td><td>+0.023</td><td>-0.063</td></tr><tr><td>-2.0561160002e-01</td><td>-0.016</td><td>+0.056</td></tr><tr><td>-5.9302199351e-01</td><td>-0.019</td><td>-0.016</td></tr><tr><td>-1.4945472270e+00</td><td>-0.034</td><td>-0.011</td></tr><tr><td>-3.4222886970e+00</td><td>-0.017</td><td>+0.014</td></tr><tr><td>-2.6351110679e+02</td><td>-35.019</td><td>-0.972</td></tr></tbody></table>

By increasing the number of regions the delayed eigenvalues cluster very close to the negative decay constants and are only clearly visible for the one-group calculation, as shown in Table IV. The table also shows that the dominant eigenvalue of the clustering, $α_p$, does not change with the number of regions.

As we show next, this follows from the fact that an infinite homogeneous medium sustains a uniform flux solution, regardless of the spatial discretization. The equation describing the energy spectrum of the uniform solution results by setting to zero spatial derivatives in Eq. (8):

$$ \left(\frac{\alpha}{v} + \Sigma - H - \bar{P}\right) \psi_{\alpha} = -K_{\alpha} \psi_{\alpha}, \quad \alpha \in \mathbb{C} \setminus \Lambda, $$

This equation is independent of the spatial discretization, be it continuous or not, and therefore its solutions are present for any value of $N_r$. Moreover, according with the analysis in Sec. III.3, these uniform solutions appear in a cluster which contains all the $N_p + 1$ leading modes.

Contrary to the few group cases, complex eigenvalues appear with 281 groups in the infinite homogeneous problem (i.e. one spatial cell and zero current at the boundary), and for different reactivities, see Fig. 2. The distribution of the prompt eigenvalues is in general weakly influenced by the reactivity.

## V. COMPARISON WITH ASAHI'S RESULTS

Asahi applied a continuous space description to the one-group, one precursor problem and computed results for a symmetric reflected critical slab reactor with a central core of half width $a = 10$ cm surrounded by two reflectors of width $b = 10$ cm.[6] He computed the first 15 eigenvalues in the delayed spectrum and the first 16 in the prompt spectrum by progressively increasing the number of zeros of the associated fluxes. Asahi presented his results with 6 and 3 digits, respectively, for the delayed and prompt eigenvalues.

We have used our numerical approach to recalculate Asahi's reactor by subdividing each cm into $M$ equal regions, for a total of $N_r = 40 \times M$ spatial DOF, and by gradually increasing $M$ until convergence. For $M = 32$ the corresponding delayed and prompt eigenvalues (15+16), out of a total of $N_{DOF} = 2560$ eigenvalues, converged to 6 digits and their

<table><caption>TABLE IV. Eigenvalue clustering in the delayed spectrum for increasing number of spatial cells, one energy group, case 4 ($α$ eigenvalues as $x.xx+y$ are to be read as $x.xxe+y$).</caption><thead><tr><th>N<sub>r</sub> = 1</th><th>N<sub>r</sub> = 2</th><th>N<sub>r</sub> = 4</th></tr></thead><tbody><tr><td>-1.3821081-3</td><td>-1.3821081-3<br>-1.2466691-2</td><td>-1.3821081-3<br>-1.2466692-2<br>-1.2466698-2<br>-1.2466699-2</td></tr><tr><td>-1.4012729-2</td><td>-1.4012729-2<br>-2.8291558-2</td><td>-1.4012729-2<br>-2.8291579-2<br>-2.8291665-2<br>-2.8291679-2</td></tr><tr><td>-3.7791806-2</td><td>-3.7791806-2<br>-4.2524306-2</td><td>-3.7791806-2<br>-4.2524320-2<br>-4.2524377-2<br>-4.2524386-2</td></tr><tr><td>-8.2586186-2</td><td>-8.2586186-2<br>-1.3304127-1</td><td>-8.2586186-2<br>-1.3304138-1<br>-1.3304182-1<br>-1.3304189-1</td></tr><tr><td>-2.0564488-1</td><td>-2.0564488-1<br>-2.9246435-1</td><td>-2.0564488-1<br>-2.9246477-1<br>-2.9246649-1<br>-2.9246678-1</td></tr><tr><td>-5.9313437-1</td><td>-5.9313437-1<br>-6.6648497-1</td><td>-5.9313437-1<br>-6.6648537-1<br>-6.6648702-1<br>-6.6648730-1</td></tr><tr><td>-1.4950512+0</td><td>-1.4950512+0<br>-1.6347760+0</td><td>-1.4950512+0<br>-1.6347767+0<br>-1.6347797+0<br>-1.6347803+0</td></tr><tr><td>-3.4228559+0</td><td>-3.4228559+0<br>-3.5545951+0</td><td>-3.4228559+0<br>-3.5545958+0<br>-3.5545988+0<br>-3.5545993+0</td></tr><tr><td>-3.5578919+2</td><td>-3.5578919+2<br>-1.1008588+7</td><td>-3.5578919+2<br>-1.2897302+7<br>-4.4033287+7<br>-7.5169271+7</td></tr></tbody></table>

associated fluxes exhibit the same number of zeros as those from Asahi's. Table V shows our 31 eigenvalues and the relative difference in % to the Asahi's values. Note that there are large differences in the prompt eigenvalues. We think that some of Asahi's values lacked precision. Our results prove that, as the number of DOF increases, our numerical technique approaches the more and more closely the full spectrum of the operator.

## VI. CONCLUSIONS

We have presented theoretical arguments about the alpha eigenvalue problem in multigroup neutron diffusion, in support
<figure><img src="image_4.png" /><figcaption>Fig. 2. Distribution on the complex plane of the prompt eigenvalues with 281 groups for the cases in Table I.</figcaption></figure>

<table><caption>TABLE V. Eigenvalues with the relative difference (%) to Asahi's results: delayed spectrum (left), prompt spectrum (right).</caption><thead><tr><th>α<sub>1,(·)</sub></th><th>r.d. (%)</th><th>α<sub>2,(·)</sub></th><th>r.d. (%)</th></tr></thead><tbody><tr><td>-0.000009</td><td>-</td><td>-6.593 × 10<sup>1</sup></td><td>+0.046</td></tr><tr><td>-0.095135</td><td>-3.307e-04</td><td>-1.340 × 10<sup>3</sup></td><td>+0.027</td></tr><tr><td>-0.098363</td><td>-2.959e-04</td><td>-3.957 × 10<sup>3</sup></td><td>-0.084</td></tr><tr><td>-0.099205</td><td>-4.805e-04</td><td>-8.107 × 10<sup>3</sup></td><td>-0.040</td></tr><tr><td>-0.099536</td><td>-2.496e-04</td><td>-1.382 × 10<sup>4</sup></td><td>+0.153</td></tr><tr><td>-0.099697</td><td>+3.656e-04</td><td>-2.098 × 10<sup>4</sup></td><td>-0.076</td></tr><tr><td>-0.099788</td><td>-3.637e-04</td><td>-2.790 × 10<sup>4</sup></td><td>+14.352</td></tr><tr><td>-0.099843</td><td>-2.955e-05</td><td>-2.961 × 10<sup>4</sup></td><td>+0.033</td></tr><tr><td>-0.099879</td><td>+2.620e-04</td><td>-3.205 × 10<sup>4</sup></td><td>-0.160</td></tr><tr><td>-0.099904</td><td>+3.194e-04</td><td>-4.114 × 10<sup>4</sup></td><td>+15.417</td></tr><tr><td>-0.099922</td><td>+3.322e-04</td><td>-5.168 × 10<sup>4</sup></td><td>-0.037</td></tr><tr><td>-0.099936</td><td>-2.921e-04</td><td>-5.736 × 10<sup>4</sup></td><td>-0.071</td></tr><tr><td>-0.099946</td><td>-9.027e-05</td><td>-6.038 × 10<sup>4</sup></td><td>+0.962</td></tr><tr><td>-0.099954</td><td>-1.342e-04</td><td>-6.950 × 10<sup>4</sup></td><td>-5.892</td></tr><tr><td>-0.099960</td><td>+1.894e-04</td><td>-8.366 × 10<sup>4</sup></td><td>-0.044</td></tr><tr><td></td><td></td><td>-9.781 × 10<sup>4</sup></td><td>+0.006</td></tr></tbody></table>

of the general behavior observed in our calculations including the distribution of the prompt and delayed spectrum as well as flux clustering. Our analysis and numerical results confirm the few theoretical and numerical results found in the literature regarding the diffusion problem.

We have used a perturbation formulation, based on the steady state eigenvalue problem of the effective multiplication factor, to derive a generalized form of the in-hour equation. As opposed to previous derivations of the in-hour equation, which have been derived by projection of the original alpha equation, in our formulation this equation arises as a necessary and sufficient condition for the existence of the alpha eigenvalue-eigenvector pair.

We have also conducted many numerical simulations for

one-dimensional homogeneous and heterogeneous configurations in multigroup diffusion theory. In all cases, our numerical results support the theoretical predictions.

The characterization of the alpha modes in multigroup neutron transport is deferred to a future study.

## APPENDIX A: ANALYSIS OF THE ONE-GROUP, HOMOGENEOUS CASE

The analysis for the one-group case is made especially simple because $\chi_i = \chi = 1$ and one can add all fissile isotopes into a single one with $\nu\Sigma_f = \sum_i \nu\Sigma_{f,i}$. Moreover, for a homogeneous medium we can also define $\beta_p = \sum_i \beta_{i \to p} \nu\Sigma_{f,i} / \nu\Sigma_f$ and cast Eq. (8) into the simplified form

$$f(\alpha)\psi = -\bar{B}\psi, \qquad (38)$$

where

$$f(\alpha) = \frac{\alpha}{\nu} + \sum_{p} \beta_{p} \frac{\alpha}{\alpha + \lambda_{p}} \nu\Sigma_{f}, \qquad (39)$$

and $\bar{B}$ is the diffusion operator with steady-state neutron fission spectrum (unity in the present case). Finally, using the boundary conditions one can construct the eigenfunctions of problem $\bar{B}\psi_{\zeta} = \zeta\psi_{\zeta}$ and then any $(\alpha, \psi_{\zeta})$, where $\alpha$ is a root of

$$f(\alpha) = \zeta \qquad (40)$$

is an alpha eigenpair.

Consider a homogeneous slab of thickness $a$. For the case with zero flux boundary conditions one has the eigenfunctions $\psi(z) \sim \sin(B_n z)$ and the eigenvalues

$$\zeta_n = \nu\Sigma_f - \Sigma_a - DB_n^2,$$

where $B_n = (n + 1)\pi/a$. Note also that

$$\zeta_n = -\rho_n \nu \Sigma_f,$$

where $\rho_n = 1 - 1/\lambda_n$ is the generalized reactivity associated to the classical reactor eigenvalues $\lambda_n$.

There are $N_p + 1$ roots of (40) and they split into $N_p + 1$ different groups; new roots are added progressively to each group as $n$ increases. We show next that the higher $N_p$ groups cluster on the right of each $-\lambda_p$, but for the group at the left of $-\lambda_{N_p}$, the eigenvalues decrease without bounds towards $-\infty$.

We examine first the behavior of the $\alpha$'s that are close to $-\lambda_p$ so that $\alpha + \lambda_p = \epsilon$ with $|\epsilon| \ll 1$ (this assumption will be proved by the following analysis). For these $\alpha$'s the critical condition (40) can be written to leading order in $\epsilon$ as $-\beta_p \lambda_p \nu \Sigma_f \epsilon^{-1} = \zeta_n$, where now $\epsilon$ plays the role of the eigenvalue $\alpha$. The solution is

$$\epsilon_n = -\frac{\beta_p \lambda_p \nu \Sigma_f}{\zeta_n} \rightarrow \alpha_n = -\lambda_p \left(1 - \frac{\beta_p \nu \Sigma_f}{\zeta_n}\right).$$

For increasing $n$, $\zeta_n$ becomes negative and increases as $(n+1)^2$ so that the factor between parentheses remains positive but diminishes very fast towards 1. Remember that we have assumed $\epsilon_n$ to be small, which implies $|\zeta_n| \gg \beta_p \lambda_p \nu \Sigma_f$. Clearly, this can only happen if $\zeta_n < 0$.
Consider, finally, the case with $\alpha \ll -\lambda_{N_p}$. Here we can neglect precursor contributions and write $\alpha = v\zeta$, which entails that $\zeta$ is negative. Note that in this case the first $\zeta$ of this group is the smallest root of (40) for $n = 0$, which we know that is $\ll -\lambda_{N_p}$. Here we have

$$\alpha_n - \alpha_{n-1} = vD \left(\frac{\pi}{a}\right)^2 (2n + 1).$$

This implies that the distance between two successive eigenvalues in this group increases with $n$ and that $\lim_{n \to \infty} \alpha_n = -\infty$. It follows that the eigenvalues in this group do not cluster.

## APPENDIX B: ASAHI PROBLEM REVISITED

The simple technique in the previous Appendix can be extended to a general analysis of a piecewise homogeneous slab geometry. The basic idea is to use the Laplacian eigenfunctions in each homogeneous piece and impose interface and boundary conditions in order to obtain a general solution for the problem which is then replaced in Eq. (38). Unfortunately, as we show next, this leads to a complicated system of trigonometric equations for the determination of the alpha spectrum.

As an example of this analytical approach we consider the Asahi problem discussed in Sec. V. This problem contains a central core of width $2a$ and a bilateral reflector of width $b$. The domain is symmetric with respect to the middle plane. We shall treat only the right half of the domain for $z \in [0, a+b]$, with vacuum boundary condition on the right at $z = a+b$ and either the current or the flux vanishing on the left at $z = 0$. The latter conditions are necessary because the initial problem is symmetric and therefore accepts even and odd solutions. Given the boundary conditions, the solutions are $\psi_c(z) \sim f_c(B_c z)$ for $z \in [0, a]$ and $\psi_r(z) \sim f_r(B_r(a+b-z))$ for $z \in [a, a+b]$. The core function $f_c$ can be either a cosine or a sine, depending on whether the alpha eigenfunction is symmetric or antisymmetric, while $f_r$ is a sine. We impose continuity of the flux and current at the core-reflector interface $z = a$:

$$A_c f_c(B_c a) = A_r f_r(B_r b), \\ A_c D_c B_c f'_c(B_c a) = -A_r D_r B_r f'_r(B_r b),$$

where the A's are two constants and $f'$ is the derivative of $f$. Dividing the second equation by the first we get the dispersion relation:

$$D_c B_c \hat{f}_c(B_c a) = B_r D_r \hat{f}_r(B_r b). \quad (41)$$

Here $\hat{f}_r$ is the cotangent function while $\hat{f}_c$ is the tangent function for the symmetric solution or minus the cotangent for the antisymmetric one.

We have now a solution that satisfies boundary and interface conditions with $B_c$ and $B_r$ constrained by Eq. (41). Next, we replace our solution in the alpha eigenvalue Eq. (38) to obtain

$$\left\{ \begin{array}{l} f(\alpha) = v\Sigma_f - \Sigma_{a,c} - D_c B_c^2, \quad (42a) \\ \frac{\alpha}{v} = -\Sigma_{a,r} - D_r B_r^2, \quad (42b) \end{array} \right.$$

where $f(\alpha)$ has been defined in (39).

Equation (41) and the two latter equations form a system of nonlinear equations for the triplet $(\alpha, B_c, B_r)$. We now proceed to eliminate $B_c$ and $B_r$ and obtain a final equation for $\alpha$. From the two equations in (42) we obtain:

$$\begin{aligned} B_r^2 &= -\frac{\alpha + v\Sigma_{a,r}}{vD_r}, \\ B_c^2 &= -\frac{f(\alpha) + v(\Sigma_{a,c} - v\Sigma_f)}{vD_c}, \end{aligned} \quad (43)$$

and use of these expressions in the dispersion relation results in a nonlinear equation $F(\alpha) = 0$ for $\alpha$. In practice one might proceed with a zero search for this nonlinear equation, where given a value of $\alpha \in \mathbb{R}$, one evaluates $B_r$ and $B_c$ via relations (43) and, then calculates $F(\alpha)$ as the difference between the left minus the right hand sides of Eq. (41).

Clearly, $B_r = 0$ is a singular point at which the solution vanishes in the reflector and in the core. Note that $B_r^2$ is negative for $\alpha > -v\Sigma_{a,r}$, in which case the cotangent $\hat{f}_r$ becomes a hyperbolic cotangent and the reflector flux, which can be taken as positive, decreases from the interface to vanish at the boundary. This case, for which the dispersion relation is easy to solve at least graphically, corresponds to the delayed spectrum calculated by Asahi, whose eigenfunctions have no zeros in the reflector. In the opposite case, with $B > 0$ or, equivalently, with $\alpha < -v\Sigma_{a,r}$, the flux in the reflector can change of sign. This is the case of Asahi's prompt spectrum. The numerical analysis fo Asahi's problem is further simplified because of the simple form of $f(\alpha)$ for the case of a single precursor, but we do not pursue this development further here.

## APPENDIX C: ROOTS OF THE GENERALIZED NORDHEIM EQUATION

Equation (20) is the condition for the existence of the eigenpair $(\alpha, \psi_\alpha)$. We write this equation in the form $\rho = \rho(\omega)$:

$$\rho = \omega \left( \Lambda + \sum_p \frac{\Gamma_p}{\omega + \lambda_p} \right),$$

where the $\Lambda(\psi_\alpha)$ and $\Gamma_p(\psi_\alpha)$ are defined in Eq. (23). The above equation can be written as a polynomial in $\omega$ of order $N_p + 1$ and therefore has as many roots but not all necessarily real. We assume that $\Lambda$ and the $\Gamma_p$'s are bounded and that $\Gamma_p \neq \Gamma_{p'}$ for $p \neq p'$. To simplify the analysis we shall assume that these coefficients are real and positive. Under these assumption we show that for any $\rho \in \mathbb{R}$ the equation has exactly $N_p + 1$ real distinct roots. As an aside, we recall that for $\rho = \rho_\lambda$ one of the roots must be $\omega = \alpha$.

Let $\mathbb{R}_0 = \mathbb{R} \setminus \Lambda$ be the real line excluding the $N_p$ points $\Lambda = \{-\lambda_p, p = 1, N_p\}$. We note that $\rho(\omega)$ is $C_\infty$ over $\mathbb{R}_0$. It is clear that every one of these points is a vertical asymptote. Let $\varepsilon > 0$, then for $\varepsilon \to 0$ we find that $\rho(-\lambda_p \pm \varepsilon) \to \mp\lambda_p\Gamma_p/\varepsilon$. Therefore, $\rho(\omega)$ goes to $-\infty$ at the right of the asymptote and to $\infty$ at its left. Moreover, the derivative

$$\partial_\omega \rho = \Lambda + \sum_p \frac{\lambda_p \Gamma_p}{(\omega + \lambda_p)^2}$$
shows that, $\rho(\omega)$ is a monotonously increasing function in $\mathbb{R}_0$. The roots of $\rho(\omega) = 0$ give the locations at which the function changes of sign. Note that there is a single root in each interval $(-\lambda_{p+1}, -\lambda_p)$ for $p = 1, N_p - 1$. Clearly $\omega = 0$ is another root and, since the function monotonously increases, there are no roots for strictly positive $\omega$. Next, we check the behavior as $\omega \to \pm\infty$ to find that, for large $\omega$, $\rho(\omega) \to \omega\Lambda$. Therefore there is a last root at some value smaller than $-\lambda_{N_p}$. The general aspect of the graph is depicted in Fig. 3.

<figure><img src="image_5.png" /><figcaption>Fig. 3. Sketch of generalized Nordheim in-hour equation for $N_p = 6$.</figcaption></figure>

Let us discuss what happens when not all the coefficients are positive. Assume that only one, say $\Lambda$, is negative. This changes the behavior for $\omega \to \pm\infty$ and now we find that $\rho(\omega) \to \mp\infty$ for $\omega \to \pm\infty$. This agrees with the fact that presently the derivative has two zeros. Hence, $\rho(\omega)$ has a positive minimum $\rho_{min} > 0$ for $\omega < -\lambda_{N_p}$ and a positive maximum $\rho_{max} < \rho_{min}$ for $\omega > -\lambda_1$. The number of roots is $N_p + 1$ except for $\rho \in (\rho_{max}, \rho_{min})$ where it is $N_p - 1$. Note that for $\rho = \rho_{max}$ or $\rho = \rho_{min}$ there is a double real root.

Finally, we consider the case when the coefficient $\Gamma_p$ is negative for one intermediary $p$ value, while all other coefficients remain positive. Here a positive minimum $\rho_{min} > 0$ appears in the interval $(-\lambda_{p+1}, -\lambda_p)$ while a positive maximum $\rho_{max} < \rho_{min}$ appears in the interval $(-\lambda_p, -\lambda_{p-1})$. The discussion is similar to the previous case. We left the analysis of more complicated cases, where more than one coefficient is negative, to the interested reader.

## REFERENCES

1. G. M. WING, *An introduction to transport theory*, Wiley (1962).

2. L. R. FOULKE and E. P. GYFTOPOULOS, "Application of the natural mode approximation to space-time reactor problems," *Nuclear Science and Engineering*, **30**, 3, 419–435 (1967).

3. R. BORSARI, T. TROMBETTI, and P. VESTRUCCI, "Time-varying modes in space-dependent nuclear reactor kinetics by a perturbation approach," *Transport Theory and Statistical Physics*, **24**, 1-3, 173–202 (1995).

4. K. DUGAN, I. ZMIJAREVIC, and R. SANCHEZ, "Cross-Section Homogenization for Reactivity-Induced Transient

Calculations," *Journal of Computational and Theoretical Transport*, **45**, 6, 425–441 (2016).

5. T. PORSCHING, "On the spectrum of a matrix arising from a problem in reactor kinetics," *SIAM Journal on Applied Mathematics*, **16**, 2, 301–317 (1968).

6. Y. ASAHI, "Theory of Omega-d Modes," *Journal of Nuclear Science and Technology*, **12**, 2, 92–106 (1975).

7. J. DEVOOGHT, "Spectrum of the multigroup-multipoint diffusion operator with delayed neutrons," *Nuclear Science and Engineering*, **67**, 2, 147–161 (1978).

8. T. GOZANI, "The concept of reactivity and its application to kinetic measurements," *Nukleonik (West Germany) Discontinued with vol. 12*, **5** (1963).

9. D. CACUCI, Y. RONEN, Z. SHAYER, J. WAGSCHAL, and Y. YEIVIN, "Eigenvalue-dependent neutron energy spectra: Definitions, analyses, and applications," *Nuclear Science and Engineering*, **81**, 3, 432–442 (1982).

10. A. HENRY, "The application of inhour modes to the description of non-separable reactor transients," *Nuclear Science and Engineering*, **20**, 3, 338–351 (1964).

11. K. YOSIDA, *Functional analysis*, Delhi., Narosa Publishing House (1965).

12. R. SANCHEZ, "Duality, Green's functions and all that," *Transport theory and statistical physics*, **27**, 5-7, 445–476 (1998).

13. G. R. KEEPIN, *Physics of nuclear kinetics*, Addison-Wesley Pub. Co. (1965).

14. A. DALL'OSSO, "Neutron spectrum kinetics in the infinite homogeneous reactor," *Annals of Nuclear Energy*, **85**, 662–669 (2015).

15. T. KATO, *Perturbation theory for linear operators*, Springer (1980).

16. J. PLANCHARD, *Méthodes mathématiques en neutronique*, Eyrolles (1995).

17. T. GOZANI, *Subcritical reactor kinetics and reactivity measurements*, Ph.D. thesis, Diss. Naturwiss. ETH Zürich, Nr. 3288, 0000. Ref.: Hälg, W.; Korref.: Schmid, P., 1962 (1962).