STUDY OF THE EIGENVALUE SPECTRA OF THE NEUTRON TRANSPORT PROBLEM IN PN APPROXIMATION

The study of the steady-state solutions of neutron transport equation requires the introduction of appropriate eigenvalues: this can be done in various different ways by changing each of the operators in the transport equation; such modifications can be physically viewed as a variation of the corresponding macroscopic cross sections only, so making the different (generalized) eigenvalue problems non-equivalent. In this paper the eigenvalue problem associated to the time-dependent problem (α eigenvalue), also in the presence of delayed emissions is evaluated. The properties of associated spectra can give different insight into the physics of the problem.


INTRODUCTION
In reactor physics, the study of the the neutronics of multiplying systems requires the solution of the transport equation in steady-state condition. In particular, the determination of the criticality condition is associated to the solution of an eigenvalue formulation of the neutron balance problem, thus allowing a non-trivial solution.
In general terms, such eigenvalue problem can be formulated in various alternative forms, each allowing different physical interpretations. The standard criticality problem is customarily approached by introducing a multiplication eigenvalue in correspondence of the fission term. In more recent times the α eigenvalue problem has gained larger interest [1,2], thanks to the associated significance in terms of asymptotic evolution of the system.
It must be pointed out as well that the solution of the eigenvalue problem for the neutron balance model can provide, in general, a set of eigenvalues associated to the corresponding eigenfunctions, which can constitute a useful basis for projection techniques of relevance in reactor physics analyses, such as higher-order perturbation theory [3,4]. It is therefore extremely interesting to approach the analysis of the eigenvalue formulation of the neutron transport model, focusing on both the possible different formulations which can be given to the problem and the methods adopted for its solution.
In this work, different formulations of the eigenvalue problems are compared in terms of characteristics of the resulting model and physical significance. The numerical solution is then achieved adopting a P N approximation for the angular variable, focusing on some specific features of the solution when comparing odd-and even-order expansions [5]. The solution approach adopted, applied to simplified configurations, allows to calculate a set of eigenvalues, thus accessing the information related to higher-order harmonics. The eigenvalue spectra for the different formulations of the problem will be compared and discussed.

Eigenvalue problem formulations
According to what is the physical phenomenon of interest for a certain system, the eigenvalue problem associated to the transport (or diffusion) equation, whereT is the time operator,L is the leakage,R is the removal,Ŝ is the scattering,F is the total fission and ϕ is the flux, can be formulated in different ways, each one providing a different interpretation of the problem [6].
The most popular approach to find the steady state neutron distribution is the one traditionally introduced by Fermi, i.e. the criticality eigenproblem, where the κ eigenvalue is introduced in the steady state balance equation as a "tuning" parameter for the fission operator, A similar approach for the stationary equation focuses on the collisions, dividing both the scattering and the fission operators by the γ eigenvalue, A more exotic and less explored alternative to study the steady state chain reaction is to look for the material density allowing the system to become critical: this objective is achieved through the introduction of a δ eigenvalue for all collision terms in the balance equation: As far as the time-dependent problem is concerned, assuming suitable initial conditions, it is possible to substitute the time operatorT with the so-called time eigenvalue, α, which has the physical meaning of the inverse of the reactor period: Such formulation, which can also be interpreted as a Laplace transform of the time-dependent balance model, is more correctly representing the asymptotic behavior of a nuclear system when delayed fission amissions are included, leading to: where C i is the i-th precursor concentration, λ i is the decay constant of the precursor family i, F d is the delayed fission operator andÊ is the delayed neutron emission density. This second formulation of the α eigenproblem has deserved recently a lot of attention, regarding both the algorithmic approach to its solution [1] and the characterization of the eigenvalue spectrum [2].
According to the kind of eigenproblem that is formulated and solved, one can retrieve different information on the system. The κ-modes are usually employed as a basis to reconstruct the flux of a perturbed system [7], while the α-modes are more useful for expanding the flux in timedependent problems. Other applications are related to the design and interpretation of experiments using the eigenvalues and their separation [8].
In addition to the practical applications, the study of the properties of the different eigenproblem spectra is significant also for theoretical interests: its properties -for instance the existence of complex eigenvalues [9], of a continuous part of the spectrum [10] or the eigenvalues spacing -are influenced by both the physical properties of the system under consideration and by the mathematical scheme employed to solve the problem.
In this work we present the main results about the eigenvalue formulations of the 1D plane geometry transport equation approximated numerically with the P N method [11]. In the following sections, an analysis of the convergence pattern of the various eigenvalue formulations with respect to the different types of boundary conditions and to the order employed to truncate the P N equations is carried out. Then, the properties of the spectrum of the different eigenproblem formulations are discussed, highlighting also the effects related to the medium heterogeneity and to the adoption of multi-group approximation.

NUMERICAL SETUP OF THE PROBLEM AND VERIFICATION
Thanks to the simplified geometry, the different operators appearing in equation (1) can be assembled directly and manipulated in order to retrieve the various eigenvalue problem formulations. To perform these tasks, we exploited the great flexibility offered by the object-oriented programming paradigm implementing a set of classes in Matlab that discretises the transport equation, builds the operators, combines them and looks for a user-defined number of eigenpairs. The code has been verified against other numerical results found in literature for the different eigenvalue formulations in some cases like two-group homogeneous and heterogeneous media. The eigenvalues are in good agreement (<50 pcm for the fundamental one) with the reference ones [12][13] [14].
The implementation of the equation numerical solver arises two interesting facts about the P N approximation. The first one is related to the fact that the spatial discretisation, performed using the centred finite difference scheme, requires that two staggered grids are defined for the even and odd equations in order to ensure consistency with the continuous equations. between P N and S N +1 . On the other hand, the choice of µ i in case of even N leaves a certain freedom, since one can choose the set of µ i as roots of L N (i.e. the same of the previous odd order) or L N +1 , removing µ = 0 as it would not yield a meaningful condition since the particles moving along this directions would not cross the boundaries (i.e., the problem is no longer differential).
The observation on the Mark BCs in case of even N is related to the fact that, moving from the odd to the subsequent even approximation, the order of the system of equations does not change, as already pointed out in the past by Davison [15]. This feature allows to recast the even P N equations into a set of P N −1 equations with modified coefficients. The influence of both the approximation orders and the BCs choice is briefly discussed in the next section.

ANGULAR ORDER CONVERGENCE
The convergence behaviour of the eigenvalues κ, γ and α can be quite different depending on the boundary conditions imposed. An inspection of figures 1-3 would allow to draw some preliminary conclusions on the P N convergence properties. At first, it is interesting to notice that, moving Figure 1: κ, γ, and α eigenvalues convergence as a function of the angular approximation order with Mark BCs using the roots of L N .
from P 1 to P 2 , the absolute error of κ and γ with respect to the reference solution (evaluated using a P 63 for all the three BCs cases) grows for the case of Mark BCs evaluated in the roots of L N , decreases for the case of Mark BCs computed in the roots of L N +1 and is about the same with Marshak BCs. As regards α, the convergence trend from P 1 to P 2 shows that, on the contrary of what is common believed, the even order approximations can provide a significant improvement with respect to the previous odd order. This fact supports the idea that, even if the even P N has the same rank of the odd system, the modified coefficients are more representative of the physics described by the problem.

EIGENVALUE SPECTRUM
In this section, we use the code to investigate some of the features of the numerical spectrum of different eigenvalue problems presented in section 2, Results for δ are not reported, as their behaviour requires a deeper analytical investigation. This fact is in principle not so surprising if one recall the physical meaning of the various eigenvalue. For instance, the κ eigenvalue meaning  is the inverse of the fractional variation of the fissile density which is required to make the system critical. Then, a real κ value smaller than 1 is expected for a subcritical system and viceversa for a supercritical one. As far as δ eigenvalue is concerned -as it is related to the number of collisions -we simply remark its physical meaning is the inverse of the fractional variation of the global density of the system which is needed to make it critical: it is not obvious that a real positive value exists as a solution for this rather unusual problem.
In Figures 4-7 the spectra evaluated are reported, with a red dot indicating the fundamental eigenvalue. Figure 4 shows the behaviour of the κ spectrum for increasing values of the angular approximation order in case of a symmetrically reflected system, with a reflector thickness of 100 cm of H 2 O and 20 cm of fissile material (for the two-group data, please refer to the URRb-H2Oa(1)-2-0-SL benchmark in [14]). Increasing the P N order does not seem to change the behaviour of the spectrum, which attempts to reproduce the continuous spectrum for reflected systems The α delayed spectrum is evaluated for a homogeneous slab with Mark boundary conditions (for the two-group data employed, please refer to the PU-2-0-SL benchmark in [14]).  Figure 5: γ eigenvalue distribution for a symmetrically reflected system.
The comparison of the results highlights that the γ spectrum has real eigenvalues, as for κ, while the time-eigenvalue spectrum shows, as expected, a distribution in the complex plane with clustering of eigenvalues as observed in other situation [2].

CONCLUDING REMARKS
In this work, the eigenvalue problem associated to the solution of the neutron balance equation is analyzed, focusing on the different formulations that can be given to the problem, with their physical significance. The role of the method adopted for its solution is discussed, focusing on the P N method and pointing out at some peculiar effects regarding the angular expansion order which deserves further investigations to assess the suitability of even-order PN for the solution of reactor physics-relevant problems.

ACNOWLEDGMENTS
The Authors are grateful to the reviewer for their very insightful comments and useful suggestions.