Comparison of Direct and Adjoint k and α -Eigenfunctions

: Modal expansions based on k -eigenvalues and α -eigenvalues are commonly used in order to investigate the reactor behaviour, each with a distinct point of view: the former is related to ﬁssion generations, whereas the latter is related to time. Well-known Monte Carlo methods exist to compute the direct k or α fundamental eigenmodes, based on variants of the power iteration. The possibility of computing adjoint eigenfunctions in continuous-energy transport has been recently implemented and tested in the development version of TRIPOLI-4 ® , using a modiﬁed version of the Iterated Fission Probability (IFP) method for the adjoint α calculation. In this work we present a preliminary comparison of direct and adjoint k and α eigenmodes by Monte Carlo methods, for small deviations from criticality. When the reactor is exactly critical, i.e., for k 0 = 1 or equivalently α 0 = 0, the fundamental modes of both eigenfunction bases coincide, as expected on physical grounds. However, for non-critical systems the fundamental k and α eigenmodes show signiﬁcant discrepancies.


Introduction
Assessing the asymptotic behaviour of a nuclear system is intimately related to computing the dominant direct (Here and throughout the manuscript, the terms 'direct' and 'forward' are equivalently used.) eigenmode and the associated dominant eigenvalue of the configuration under analysis, which is central in several applications in reactor physics, including pulsed neutron reactivity measurements [1] and reactor period analysis [2][3][4]. Furthermore, for several key reactor coefficients, such as kinetics parameters, sensitivities to nuclear and material data, or perturbations [5][6][7][8][9][10], bilinear forms involving both the direct and the adjoint fundamental eigenfunctions are required. The most common bases for eigenfunction expansions are those related to the k-eigenvalues and to the α-eigenvalues, respectively [11]. Calculations of dominant k or α eigenvalues/eigenfunctions try to assess the asymptotic reactor behaviour, each with a distinct point of view: the former basically determines the shape of the neutron population after a large number of fission generations, whereas the latter after a sufficiently long time. When the reactor is exactly critical, i.e., for k 0 = 1 or equivalently α 0 = 0, the fundamental modes of both eigenfunction bases coincide, as expected on physical grounds [11]. However, for systems far from criticality the fundamental k and α eigenmodes show discrepancies and (with the possible exception of very simple cases involving single-speed transport) are not related to each other in any trivial manner [11,12]. Such discrepancies have been observed even for very small deviations from criticality [13]. Since both direct and adjoint eigenmodes are involved in the calculation of key reactor parameters, the investigation of the behaviour of the eigenfunctions might shed some light on the behaviour of such parameters at and close to the critical point.
In the context of Monte Carlo simulation, the analysis of the prompt direct fundamental k and α-eigenmodes (neglecting delayed neutron contributions) has been previously carried out for a set of sub-and super-critical systems based on homogeneous and heterogeneous cores with fast and thermal spectra [12]. In this work, we will revisit and extend these findings on some benchmark configurations, namely Godiva-like spheres inspired by the seminal work of D.E. Cullen [12] and the CROCUS core [14,15], with a twofold aim. First, we will explicitly include the effects of delayed neutrons, which have been originally neglected for the sake of simplicity. We will thus determine whether the presence of delayed contributions has an impact on the shape of the eigenfunctions, and in particular whether the discrepancies between k and α-eigenmodes increase or decrease. Then, we will examine the behaviour of the fundamental adjoint eigenmodes for both k and α-eigenvalue problems, whose comparison has not been addressed so far, to the best of our knowledge.

Tools and Methods
For the numerical simulations presented in the following, we have used the development version of TRIPOLI-4 ® [16]. For direct k-eigenvalue calculations, the fundamental mode ϕ k0 satisfying is computed by using the standard power iteration method [16]. Here L denotes the net loss operator, F p the prompt fission operator, F d,j the delayed fission operator for precursor family j, and the other notation is standard. For α-eigenvalue problems, in the form in the development version of TRIPOLI-4 ® the fundamental mode ϕ α0 is determined by using the α-k power iteration [17]. The α-k method was originally proposed for prompt decay constants [18] and later extended to the general case with neutrons and precursors [2,17]. The basic idea is to iteratively search for the dominant α value that makes the α-eigenvalue equation exactly critical with respect to a fictitious k-eigenvalue applied to the production terms. For positive α, a "capture" cross section α/υ is taken into account while applying a modified power iteration [12,18]. For negative α, the contribution −α/υ was originally interpreted as a "production" term [12,18]. The standard implementation of this algorithm has been however shown to be numerically unstable, possibly leading to abnormal termination: an improved algorithm that overcomes the limitations for the case of negative α has been proposed [17], introducing a special "copy" operator (at the right hand side of Equation (2)) that preserves particle balance between the power iteration cycles, namely with a multiplicity υ α = 2. The operator F α is compensated by a "production" cross section -α/υ at the left hand side of Equation (2). The term υ α can be interpreted as the average number of (copy) neutrons produced by the α-production operator having a delta kernel. Details concerning the implementation and the possible convergence issues for the case of negative α are provided in [17]. Finally, the term λ j /(λ j + α), λ j being the decay constant for precursor family j, can be interpreted as a weight multiplier for the delayed neutron yield, observing that the dominant α satisfies α 0 > −min j λ j [2,17]. Concerning the adjoint k-eigenmodes ϕ † k , which satisfy in TRIPOLI-4 ® we have adopted the Iterated Fission Probability (IFP) method [19], by closely following the strategy proposed in [5,6]. The introduction of the IFP method has paved the way to obtaining the fundamental adjoint flux ϕ † k 0 for k-eigenvalue problems in continuous-energy Monte Carlo simulations: the adjoint flux ϕ † k 0 is equated to the neutron importance I k , which can be then estimated by running a direct calculation [5,6]. The neutron importance I k is obtained by recording the descendants after M latent generations for an ancestor injected into the system at coordinates r, Ω, E (neutrons are promoted to the next generation by fission events). By building upon these ideas, we have recently introduced a novel method to compute the fundamental adjoint flux ϕ † α 0 for α-eigenvalue problems, which solves by resorting to a generalized version of IFP (G-IFP) [13]. The fundamental adjoint flux ϕ † α 0 can be again equated to the neutron importance I α , i.e., can be estimated by recording the descendants after M latent generations for an ancestor injected into the system at coordinates r, Ω, E. The only difference with respect to the regular IFP method is that for α-eigenvalue problems additional events, other than fissions, may contribute to promoting the neutrons to the next generation. For positive α, the additional term α/υ acts as a sterile capture: neutrons can thus contribute to the importance only being promoted to the next generation by prompt and delayed fission. However, for negative α, neutrons can contribute to importance also via the α-production term, associated to the copy operator with cross section −α/υ. In both cases, the weight of the delayed neutrons is assigned a correction factor λ j /(λ j + α). For a thorough explanation of the G-IFP method, see [13].

Analysis of the Godiva-Like Spheres
In view of characterizing the behaviour of the direct and adjoint fundamental modes for k and α-eigenvalue problems, we have selected two simple benchmark configurations, both inspired from test-cases previously considered in the literature [12]. The first configuration consists in a bare sphere of uranium, whose specifications are taken from [12] (Problem I) and are very close to those of the standard Godiva benchmark [20]. In particular, the radius of the sphere is equal to 8.7407 cm and the uranium isotopic composition (normalized with respect to the uranium density) consists of 93.7695% atoms of U 235 , 5.2053% atoms of U 238 and 1.0252% atoms of U 234 . The system is spatially homogeneous, and the neutron spectrum is fast. The second configuration is also taken from [12] (Problem III) and corresponds to a sphere of uranium with equal radius and uranium isotopic composition, surrounded by a thick water reflector: the system is spatially heterogeneous with a total radius of 38.7407 cm and a strong thermal component. The water density is equal to 1 g/cm 3 with 2 atoms of H 1 and 1 atom of O 16 . For the sake of simplicity, we will call these configurations Problem I and Problem III.

Description of the Benchmark Configurations
For both cases, starting from the specifications given in [12], we have adjusted the uranium density in order to obtain slightly sub-critical and sightly super-critical configurations, with the aim of examining the effects of slight deviations from criticality on the shape of the direct and adjoint eigenmodes. The chosen values of uranium density for each configuration are shown in Table 1. Table 1. Uranium densities for benchmark configurations. Uranium isotopic mass fractions and water properties are equal to those described in the reference [12]. The simulation results displayed in the following have been obtained by resorting to TRIPOLI-4 ® . The forward simulations are performed via the power iteration method for the k-eigenvalue problem, and the α-k power iteration method for the α-eigenvalue problem. The corresponding numerical simulation parameters are presented in Table 2. The adjoint simulations are performed via the IFP method for the k-eigenvalue problem, and the G-IFP method for the α-eigenvalue problem. The corresponding numerical simulation parameters are presented in Table 3. All flux distributions have been scored into 281 energy meshes. Nuclear data for our calculations have been taken from the JEFF3.1.1 library, where all fissile isotopes have 8 families of precursors [21]. Table 2. Numerical simulation parameters for benchmark configurations: forward simulations.

Configuration
Active Cycles Inactive Cycles Particles I sub-critical 5 × 10 3 10 3 10 5 I super-critical 5 × 10 3 10 3 10 5 III sub-critical 2 × 10 3 2 × 10 3 10 4 III super-critical 2 × 10 3 (*) 2 × 10 3 10 4 (*) 5 × 10 3 for the α-simulation including delayed neutron contributions. The fundamental eigenvalues k 0 and α 0 computed in the corresponding simulations for Problem I by including and neglecting the delayed neutron contributions are given in Tables 4 and 5, respectively. It is worth noting that the super-critical configuration in Problem I becomes sub-critical when delayed neutron contributions are neglected in the calculations. Moreover, the fundamental eigenvalue k 0 is reduced by approximately 650 pcm when the delayed contributions are not considered. Concerning the α-eigenvalue formulation, the absolute value of α 0 is smaller than 1 s −1 by including delayed neutrons, whereas the absolute value of this eigenvalue is larger than 10 5 s −1 when neglecting delayed neutrons. Neglecting the presence of the delayed neutrons in this fast spectrum system implies thus a decrease of the reactor period of approximately 5 orders of magnitude.  For illustration, the shapes of direct and adjoint eigenmodes ϕ α 0 (with and without delayed neutron contributions) for the sub-critical configuration and the super-critical configuration of Problem I are shown in Figures 1 and 2  For illustration, the shapes of direct and adjoint eigenmodes   For illustration, the shapes of direct and adjoint eigenmodes In order to quantitatively assess these differences, in Figures 3 and 4 we show the ratios between α-and keigenfunctions for the direct and adjoint problem, respectively. For the direct eigenfunctions, deviations are overall rather small (see Figure 3). In the sub-critical configuration, we have ϕ α 0 < ϕ k 0 in the fast region and vice-versa in the epithermal region, both with and without delayed neutrons (Figure 3a). For the super-critical configuration, the situation is different: in the fast region, we have ϕ α 0 < ϕ k 0 without delayed neutrons and ϕ α 0 > ϕ k 0 with delayed neutrons; in the epithermal region the behaviour is inverted. This is possibly due to the fact that in the super-critical configuration the sign of α 0 changes with or without delayed neutrons. In the thermal region, very few neutrons contribute to the direct eigenfunctions (as expected from a fast neutron spectrum system), although statistical uncertainty prevents from drawing solid conclusions.
neutrons contribute to the direct eigenfunctions (as expected from a fast neutron spectrum system), although statistical uncertainty prevents from drawing solid conclusions.
As for the adjoint eigenfunctions, deviations are somewhat stronger when delayed neutrons are not considered (see Figure 4), especially in the resonance region. On the contrary, when delayed neutron contributions are taken into account deviations of 0 † k ϕ from 0 † α ϕ become much smaller. Overall, neglecting the presence of delayed neutrons leads to 0 † α ϕ < 0 † k ϕ outside the resonance region and the difference between the two distributions is larger for the more sub-critical configuration.
Direct and adjoint eigenmodes computed by including delayed neutrons show similar distributions. In principle, the absolute value of α0 drops around 10 −1 s −1 , therefore the term (α0/υ) 0 α ϕ is significantly reduced. Moreover, the weight multipliers for both the kdelayed fission operator (1/k0) and the α-delayed fission operator around the unit value. For this reason, the k-and the α-eigenvalue problems presented in Equations (1) and (2) for the direct formulation and in Equations (4) and (5) for the adjoint formulation would be close to each other.   configuration, with (red squares) and without (blue circles) delayed contributions.

Problem III
The fundamental eigenvalues k0 and α0 computed in the corresponding simulations for Problem III with and without the delayed neutron contributions are given in Tables 6  and 7, respectively. It is worth noting that similarly to Problem I the system described in Problem III is sub-critical if delayed neutron contributions are neglected in the calcula- As for the adjoint eigenfunctions, deviations are somewhat stronger when delayed neutrons are not considered (see Figure 4), especially in the resonance region. On the contrary, when delayed neutron contributions are taken into account deviations of ϕ † k 0 from ϕ † α 0 become much smaller. Overall, neglecting the presence of delayed neutrons leads to outside the resonance region and the difference between the two distributions is larger for the more sub-critical configuration.
Direct and adjoint eigenmodes computed by including delayed neutrons show similar distributions. In principle, the absolute value of α 0 drops around 10 −1 s −1 , therefore the term (α 0 /υ)ϕ α 0 is significantly reduced. Moreover, the weight multipliers for both the k-delayed fission operator (1/k 0 ) and the α-delayed fission operator λ j /(λ j + α 0 ) are both around the unit value. For this reason, the kand the α-eigenvalue problems presented in Equations (1) and (2) for the direct formulation and in Equations (4) and (5) for the adjoint formulation would be close to each other.

Problem III
The fundamental eigenvalues k 0 and α 0 computed in the corresponding simulations for Problem III with and without the delayed neutron contributions are given in Tables 6  and 7, respectively. It is worth noting that similarly to Problem I the system described in Problem III is sub-critical if delayed neutron contributions are neglected in the calculations. Moreover, the fundamental eigenvalue k 0 is reduced by approximately 700 pcm when the delayed contributions are not considered. Concerning the α-eigenvalue formulation, the absolute value of α 0 is again smaller than 1 s −1 by including delayed neutrons, whereas the absolute value of this eigenvalue is between 10 2 s −1 and 10 3 s −1 when neglecting delayed neutrons. Neglecting the presence of the delayed neutrons in this thermal system implies a decrease of the reactor period of approximately 3 orders of magnitude.
The shapes of the direct eigenmodes ϕ k 0 and ϕ α 0 (with and without delayed neutron contributions) for the sub-critical configuration of Problem III are shown in Figure 5a; the adjoint eigenmodes are also shown in Figure 5b. The direct and the adjoint distributions for the super-critical configuration of the same problem are shown in Figure 6. The behaviour of the adjoint eigenmodes shows noticeable differences with respect to what observed in spatially large systems, i.e., MOX and UOX assembly configurations [19], where the thermal component was higher and the fast component was lower as compared to Problem III. The strong impact of the fast component in our example, and the milder impact of the thermal component, can be justified by the fact that Problem III is strongly spatially heterogeneous, with a fast spectrum localized in the fissile lump and a thermal spectrum localized in the moderator. Due to normalization the amplitude of ϕ † shown in Figures 5 and 6 is smaller in the thermal region and larger in the fast region if compared to the results obtained from MOX and UOX assembly configurations. For comparison, all curves have been normalized.
Slight but significant deviations due to the kind of eigenfunction (either k-or α-) and to the presence of delayed contributions are again visible, especially for the adjoint fluxes. [19], where the thermal component was higher and the fast component was lower as compared to Problem III. The strong impact of the fast component in our example, and the milder impact of the thermal component, can be justified by the fact that Problem III is strongly spatially heterogeneous, with a fast spectrum localized in the fissile lump and a thermal spectrum localized in the moderator. Due to normalization the amplitude of † ϕ shown in Figures 5 and 6 is smaller in the thermal region and larger in the fast region if compared to the results obtained from MOX and UOX assembly configurations. For comparison, all curves have been normalized. Slight but significant deviations due to the kind of eigenfunction (either k-or α-) and to the presence of delayed contributions are again visible, especially for the adjoint fluxes.  pared to Problem III. The strong impact of the fast component in our example, and the milder impact of the thermal component, can be justified by the fact that Problem III is strongly spatially heterogeneous, with a fast spectrum localized in the fissile lump and a thermal spectrum localized in the moderator. Due to normalization the amplitude of † ϕ shown in Figures 5 and 6 is smaller in the thermal region and larger in the fast region if compared to the results obtained from MOX and UOX assembly configurations. For comparison, all curves have been normalized. Slight but significant deviations due to the kind of eigenfunction (either k-or α-) and to the presence of delayed contributions are again visible, especially for the adjoint fluxes.
(a) (b)  The corresponding ratios for the direct and adjoint eigenfunctions are displayed in Figures 7 and 8, respectively. For the direct eigenfunctions, when delayed neutron contributions are taken into account deviations are rather small (see Figure 7) for both the sub-and super-critical configurations. However, it is possible to notice the presence of a deviation between the kand α-eigenfunctions: for the sub-critical configuration, we have ϕ α 0 < ϕ k 0 for energy values larger than 0.6 MeV and ϕ α 0 > ϕ k 0 for energy values larger between 10 keV and 0.6 MeV; for the super-critical configuration an opposite behavior is noticeable in the same energy ranges. This inversion is justified by the transition from a sub-critical to a super-critical system. For a sub-critical configuration, the k-eigenvalue formulation hardens the energy spectrum by artificially increasing the amplitude of fission operator by a factor 1/k 0 , whereas the α-eigenvalue formulation promotes the thermal spectrum by the term α 0 /υ and at the same modifies the delayed fission operator by the factor γ = λ/(α 0 + λ). Observe that we have γ > 1 for negative α 0 and γ < 1 for positive α 0 . The presence of delayed neutrons shifts the behaviour where ϕ α 0 < ϕ k 0 towards 0.6 MeV which is around the average emission energy for delayed neutrons [22]. According to the α-eigenvalue formulation, the delayed fission operator is now rescaled by a factor λ/(λ + α 0 ), with λ = β/∑ j (β j /λ j ) = 0.0768 ± 0.0006 s −1 for U 235 [23]. This factor is still smaller than 1/k 0 , hence ϕ α 0 is still smaller than ϕ k 0 at fast energy range. The symmetric argument can be applied for the super-critical case. sion operator by a factor 1/k0, whereas the α-eigenvalue formulation promotes the thermal spectrum by the term α0/υ and at the same modifies the delayed fission operator by the factor γ = λ /(α0 + λ ). Observe that we have γ > 1 for negative α0 and γ < 1 for positive α0.
The presence of delayed neutrons shifts the behaviour where 0 α ϕ < 0 k ϕ towards 0.6 MeV which is around the average emission energy for delayed neutrons [22]. According to the α-eigenvalue formulation, the delayed fission operator is now rescaled by a factor λ /( λ + α0), with λ = β/∑j(βj/λj) = 0.0768 ± 0.0006 s −1 for U 235 [23]. This factor is still smaller than 1/k0, hence 0 α ϕ is still smaller than 0 k ϕ at fast energy range. The symmetric argument can be applied for the super-critical case.
The ratio shown for energies between 10 keV and 20 MeV for the super-critical configuration seems smaller compared to the one computed for the sub-critical configuration. This result is coherent with the former configuration being closer to the critical state (super critical configuration) with respect to the latter (sub-critical configuration). For energy ranges smaller than 10 keV, no significant differences are visible.  ϕ is shifted towards higher energies for a sub-critical system due to the 1/k0 factor which artificially increases the number of fissions. This behaviour is smoothed in the supercritical configuration (which is sub-critical, if delayed contributions are neglected) due to the system being closer to the critical state. As for the adjoint eigenfunc- The ratio shown for energies between 10 keV and 20 MeV for the super-critical configuration seems smaller compared to the one computed for the sub-critical configuration. This result is coherent with the former configuration being closer to the critical state (super critical configuration) with respect to the latter (sub-critical configuration). For energy ranges smaller than 10 keV, no significant differences are visible.
On the contrary, for the simulations excluding delayed neutrons deviations become larger: in the fast and epithermal region, we have ϕ α 0 < ϕ k 0 , whereas ϕ α 0 > ϕ k 0 in the thermal region. Again, if only prompt neutrons are considered, the fundamental eigenmode ϕ k 0 is shifted towards higher energies for a sub-critical system due to the 1/k 0 factor which artificially increases the number of fissions. This behaviour is smoothed in the supercritical configuration (which is sub-critical, if delayed contributions are neglected) due to the system being closer to the critical state. As for the adjoint eigenfunctions, strong deviations in the fast region (E > 0.1 MeV) are observed for the sub-critical configuration, when delayed neutrons are not considered (see Figure 8). In all the other configurations, no significant differences are visible.

Analysis of the Effective Kinetics Parameters
Based on the observed discrepancies between the fundamental kand α-eigenmodes, it would be interesting to assess which is the practical impact on the reactor parameters that depend on these quantities. In this respect, a prominent example is represented by the kinetics parameters, which are indeed bilinear forms depending on both the forward and the adjoint eigenmodes. The expressions of the effective mean generation time Λ eff(α,k) and the effective delayed fraction β j eff(α,k) respectively read The kinetics parameters, in turn, influence the system reactivity, through the in-hour (Nordheim) formula [11,24]. The static reactivity ρ k and the dynamic reactivity ρ α follow from [11][12][13]24,25] and The static reactivity ρ k depends only be the fundamental eigenvalue k 0 , whereas the dynamic reactivity ρ α requires the computation of Λ e f f ,α and β j e f f ,α in addition to the fundamental eigenvalue α 0 . The IFP method and the G-IFP method allow the computation of the bilinear forms of the kind < ϕ † k , Aϕ k > and < ϕ † α , Aϕ α > respectively, given a generic operator A. Both k and α weighted effective kinetics parameters have been estimated by resorting the methods implemented in the development version of TRIPOLI-4 ® [7].
The effective kinetics parameters estimated for Problem I are shown in Tables 8 and 9 for the sub-critical configuration. Values for β eff weighted according to the k-formulation are statistically compatible to those weighted according to the α-formulation. Simulation parameters for the evaluation of these parameters are the same as those shown in Table  2, estimated by considering 20 latent generations. A slight difference is observed in the values of Λ eff , whereas a relatively larger difference is observed between dynamic and static reactivity. The latter deviation may be justified by the fact that the dynamic reactivity ρ α depends on the eigenfunction distributions integrated for the estimation of Λ eff,α and β eff,α (Equation (9)), whereas the static reactivity ρ k only depends on the fundamental eigenvalue k 0 . Significant differences are observed for both reactivity and effective mean generation time values when only prompt neutrons are considered. In this case, the difference on the adjoint eigenmodes from Figure 4a plays a significant role in weighting the kinetics parameters.  Table 9.
Effective kinetics parameters of Problem I, sub-critical configuration without delayed contributions.

Parameters <ϕ
45.6 ± 0.3 5.728 ± 0.002 The results obtained from the super-critical configuration of the same problem are shown in Tables 10 and 11. The system including delayed neutrons is super-critical and all kinetics parameters are statistically compatible. Conversely, a non negligible discrepancy is still noticeable between static and dynamic reactivities for the configuration without delayed contributions. Overall, minimal differences have been found between direct and adjoint eigenmodes according to the kand the α-eigenvalue formulations, which is mirrored in the small discrepancies observed in the effective kinetic parameters. For Problem III, the effective kinetics parameters computed for the sub-critical configuration are shown in Tables 12 and 13. The results from the simulation including delayed contributions shows minimal discrepancies for the values of Λ eff and β eff , whereas a clear difference is found between ρ k and ρ α . Conversely, the differences in the eigenmode distributions from Figure 7and   The parameters describing the super-critical configuration of Problem III are shown in Tables 14 and 15. The presence of delayed neutrons and the proximity to the critical state leads to statistically compatible values of the kinetics parameters. The results obtained from the simulation including only prompt neutrons show minimal differences for Λ eff values, a slightly larger discrepancy is found between static and dynamic reactivity.
12.93 ± 0.05 12.29 ± 0.05 In conclusion, we have investigated the effective kinetics parameters related to the Godiva-like benchmark problems. As expected, the differences between kand α-eigenmode distributions are mirrored in the discrepancies between the corresponding kinetics parameters. Overall, the reactivity ρ is more affected by the choice of the eigenvalue formulation than the other kinetics parameters Λ eff and β eff .
The kinetics parameters and the associated reactivities are crucial for the control and the safety of nuclear reactors. A comparison with existing measurements (which typically involve neutron noise detection combined with the application of a fitting procedure based on a formulation of the in-hour equation) might help in discriminating whether the k or a formulations have a prominent advantage over each other for the interpretation of these parameters.

Analysis of the Crocus Benchmark
In order to take into account a more realistic configuration, and ascertain whether the conclusions reached in the previous section hold true for larger reactor cores, in this section we will consider two configurations of the CROCUS critical facility, operated at the Ecole Polytechnique Fédérale de Lausanne (EPFL), Switzerland. Thanks to its detailed description and careful measurements, the CROCUS core has been selected as an international benchmark for reactivity, kinetics parameters and reactor period calculations [14,15].
CROCUS is an open-tank type zero-power reactor, characterized by two fuel regions and moderated by light water. A radial section of this reactor is shown in Figure 9: the outer fuel rods (green) are composed of metallic uranium at 0.947 wt% 235 U/U with 2.917 cm pitch, whereas the inner fuel rods (orange) are composed by UO 2 at 1.806 wt% 235 U/U with 1.837 cm pitch. Details on the number and positions of these fuel rods are given in reference [15]. The core can be modeled as a cylinder with a diameter of about 60 cm and a height of 100 cm. The critical state of the system is controlled by the level of light water filling the reactor from the bottom cadmium plate. The critical configuration is achieved at a water level of 91.66 cm.
k or a formulations have a prominent advantage over each other for the interpre these parameters.

Analysis of the Crocus Benchmark
In order to take into account a more realistic configuration, and ascertain wh conclusions reached in the previous section hold true for larger reactor cores, in tion we will consider two configurations of the CROCUS critical facility, operat Ecole Polytechnique Fédérale de Lausanne (EPFL), Switzerland. Thanks to its det scription and careful measurements, the CROCUS core has been selected as an tional benchmark for reactivity, kinetics parameters and reactor period cal [14,15]. CROCUS is an open-tank type zero-power reactor, characterized by two gions and moderated by light water. A radial section of this reactor is shown in the outer fuel rods (green) are composed of metallic uranium at 0.947 wt% 235 U 2.917 cm pitch, whereas the inner fuel rods (orange) are composed by UO2 at 1. 235 U/U with 1.837 cm pitch. Details on the number and positions of these fuel rods a in reference [15]. The core can be modeled as a cylinder with a diameter of abo and a height of 100 cm. The critical state of the system is controlled by the leve water filling the reactor from the bottom cadmium plate. The critical configu achieved at a water level of 91.66 cm.
A preliminary comparison for this facility between fundamental α-eigenp time dependent calculations has been carried out in [26]. Moreover, the kinetics ters of the CROCUS core for the k-eigenvalue formulation has been previously co in [4] and the associated reactivity was examined in [27]. In the following, we w tigate the behaviour of the fundamental forward and adjoint modes of the k and lations for this core, and we will then examine the impact of their respective shap kinetics parameters and on the reactivity. For this purpose, we will consider two ical configurations of the reactor, both obtained by lowering the water level. Co tion H1 is characterized by a water level equal to 90.8 cm, which induces a slig critical state. Configuration H2 is characterized by a water level equal to 80 cm corresponds to a more sub-critical condition. The fundamental eigenmodes will puted by using TRIPOLI-4 ® over 111 energy meshes and along 14 fuel pin positio the core center to the outer region, denoted by the red rectangular regions in Figu number of cycles and the particles per cycle simulated for these configurations a in Table 16 for the forward simulations and in Table 17 for the adjoint simulation  A preliminary comparison for this facility between fundamental α-eigenpairs and time dependent calculations has been carried out in [26]. Moreover, the kinetics parameters of the CROCUS core for the k-eigenvalue formulation has been previously computed in [4] and the associated reactivity was examined in [27]. In the following, we will investigate the behaviour of the fundamental forward and adjoint modes of the k and α formulations for this core, and we will then examine the impact of their respective shapes on the kinetics parameters and on the reactivity. For this purpose, we will consider two sub-critical configurations of the reactor, both obtained by lowering the water level. Configuration H1 is characterized by a water level equal to 90.8 cm, which induces a slightly sub-critical state. Configuration H2 is characterized by a water level equal to 80 cm, which corresponds to a more sub-critical condition. The fundamental eigenmodes will be computed by using TRIPOLI-4 ® over 111 energy meshes and along 14 fuel pin positions, from the core center to the outer region, denoted by the red rectangular regions in Figure 9. The number of cycles and the particles per cycle simulated for these configurations are listed in Table 16 for the forward simulations and in Table 17 for the adjoint simulations.

Analysis of the Fundamental Eigenmodes
The fundamental eigenvalues k 0 and α 0 computed in the corresponding simulations for H1 and H2 configurations of the CROCUS reactor (with and without precursor contributions) are given in Tables 18 and 19, respectively. As expected, the computed values show a slight sub-critical level for configuration H1 and a larger sub-critical state for configuration H2. According to the results obtained from the k-eigenvalue formulation, the criticality level of both configurations decreases by approximately 780 pcm when delayed neutrons are neglected. Concerning the α-eigenvalue formulation, the absolute value of α 0 is smaller than 1.2 × 10 −2 s −1 by including delayed neutrons, whereas the absolute value of this eigenvalue is larger than 10 2 s −1 for the simulation without precursor contributions.  The shapes of direct eigenmodes ϕ k 0 and ϕ α 0 (with and without delayed neutron contributions) are shown as a function of energy ( Figure 10) and of the fuel pin position ( Figure 11) for H1 (a) and H2 (b) configurations. The adjoint eigenmodes are shown as a function of the fuel pin position in Figure 12. All curves have been normalized. No major differences can be easily spotted in these figures, so that we have computed the ratios ϕ α 0 /ϕ k 0 and ϕ † α 0 /ϕ †    Figure 13 shows the ratios of the direct fundamental eigenmodes in the energy domain. The results obtained for configuration H1 (Figure 13a) are similar to those previously discussed for the Problem III of the benchmark configurations. When delayed neutrons are considered, no differences are visible between 0 k ϕ (E) and 0 α ϕ (E). This is mainly due to the system reactivity being close to critical (about −50 pcm, as shown in Table 18): the k and α formulation are supposed to be very close to each other in this regime. If precursors are disregarded in the simulation, for this configuration the fundamental keigenmode is slightly different from the fundamental α-eigenmode: a shift towards high energy values is observed. The same behaviour is found for the H2 configuration ( Figure  13b), characterized by a larger effect due to a larger sub-critical level with respect to the previous case (about −1600 pcm, as shown in Table 18). The effect of precursors is clearly visible for this configuration: the k-eigenmode displays significant discrepancies with respect to the α-eigenmode towards high energies for this sub-critical system, but   Figure 13 shows the ratios of the direct fundamental eigenmodes in the energy domain. The results obtained for configuration H1 (Figure 13a) are similar to those previously discussed for the Problem III of the benchmark configurations. When delayed neutrons are considered, no differences are visible between ϕ k 0 (E) and ϕ α 0 (E). This is mainly due to the system reactivity being close to critical (about −50 pcm, as shown in Table 18): the k and α formulation are supposed to be very close to each other in this regime. If precursors are disregarded in the simulation, for this configuration the fundamental k-eigenmode is slightly different from the fundamental α-eigenmode: a shift towards high energy values is observed. The same behaviour is found for the H2 configuration (Figure 13b), characterized by a larger effect due to a larger sub-critical level with respect to the previous case (about −1600 pcm, as shown in Table 18). The effect of precursors is clearly visible for this configuration: the k-eigenmode displays significant discrepancies with respect to the α-eigenmode towards high energies for this sub-critical system, but ϕ α 0 (E) > ϕ k 0 (E) only for E > 0.6 MeV, which is again similar to the findings of the Problem III configuration. Precursor contributions minimize the discrepancies between the two eigenvalue formulations and the delayed spectra move the threshold for ϕ α 0 (E) > ϕ k 0 (E) at higher energy.
For the sake of completeness, we show the ratios of these eigenfunctions as a function of the fuel pin positions in Figure 14 for the direct formulation and in Figure 15 for the adjoint formulation. Overall, the behaviour of these two eigenmodes is similar: within uncertainty limits, no major differences can be detected with respect to the spatial coordinate.

Analysis of the Effective Kinetics Parameters
The effective kinetics parameters have been computed for the two CROCUS configurations by using the IFP and G-IFP of TRIPOLI-4 ® . Table 20 shows the parameters computed for the H1 configuration with precursor contributions: the results obtained for both eigenvalue formulations are statistically compatible. As expected, the proximity to the critical level of this configuration implies very close values of the k-weighted and α-weighted effective kinetics parameters. Similar results are found in Table 21 by neglecting precursor contributions. For the sake of completeness, we show the ratios of these eigenfunctions as a function of the fuel pin positions in Figure 14 for the direct formulation and in Figure 15 for the adjoint formulation. Overall, the behaviour of these two eigenmodes is similar: within uncertainty limits, no major differences can be detected with respect to the spatial coordinate.  For the sake of completeness, we show the ratios of these eigenfunctions as a function of the fuel pin positions in Figure 14 for the direct formulation and in Figure 15 for the adjoint formulation. Overall, the behaviour of these two eigenmodes is similar: within uncertainty limits, no major differences can be detected with respect to the spatial coordinate.

Analysis of the Effective Kinetics Parameters
The effective kinetics parameters have been computed for the two CROCUS configurations by using the IFP and G-IFP of TRIPOLI-4 ® . Table 20 shows the parameters computed for the H1 configuration with precursor contributions: the results obtained for both eigenvalue formulations are statistically compatible. As expected, the proximity to the critical level of this configuration implies very close values of the k-weighted and αweighted effective kinetics parameters. Similar results are found in Table 21 by neglecting precursor contributions.
The effective kinetics parameters of the H2 configuration with precursor contributions are shown in Table 22. A discrepancy is observed between the static and the dynamic reactivity, whereas the average values of all the other kinetics parameters are within one standard deviation for the two eigenvalue formulations. Table 23 shows the parameters obtained when neglecting precursor contributions: no significant discrepancies are observed in the computed values.   The effective kinetics parameters of the H2 configuration with precursor contributions are shown in Table 22. A discrepancy is observed between the static and the dynamic reactivity, whereas the average values of all the other kinetics parameters are within one standard deviation for the two eigenvalue formulations. Table 23 shows the parameters obtained when neglecting precursor contributions: no significant discrepancies are observed in the computed values. −1621 ± 2 −1638 ± 10 Λ eff [µs] 48.59 ± 0.03 48.39 ± 0.03

Conclusions
Inspired by the analysis originally proposed by D. E. Cullen [12], we have applied Monte Carlo methods for the estimation of k and α fundamental eigenmodes. We have extended the findings discussed in [12] in two directions, by addressing the evaluation of the fundamental adjoint eigenmodes and the influence of precursor contributions. Additional information regarding the discrepancies between the two eigenvalue formulations was found by assessing the effective kinetics parameters weighted by the kor by the αeigenmodes. We have focused our attention on the analysis of two Godiva-like benchmark configurations and the CROCUS reactor. In this way, we have explored thermal and fast spectra, homogeneous and heterogeneous media, simplified and realistic systems.
Significant, albeit globally small, differences have been detected, as expected on physical grounds based on previous investigations. In particular, we have recovered the same behaviour previously analyzed [12] in the energy domain for the distribution of the direct fundamental eigenfunctions ϕ k 0 and ϕ α 0 , and we have found similar discrepancies in the corresponding adjoint fundamental distributions ϕ † k 0 and ϕ † α 0 . Moreover, the presence of precursors has a non-trivial influence on the eigenfunctions, and this impact has been carefully examined for each configuration. Overall, the presence of precursors reduces the discrepancies between the two eigenmodes.
As a general remark, based on the configurations investigated here, it seems that the discrepancies between the kand α-eigenfunctions are enhanced by the presence of strong spatial heterogeneity, such as for a core surrounded by a thick moderator/reflector. In this case, the system will be characterized by multiple time scales (as shown in [12]), related to the different times required by the neutrons to explore the multiplying and the diffusing region. Then, it appears that the α-eigenvalue formulation is more sensitive to these different time scales than the k-eigenvalue formulation, which is coherent with α being related to the time behaviour of the system and k being related to the fission generation behavior. These discrepancies on the eigenmodes are mirrored in the kinetics parameters and on the reactivity. In this respect, it is interesting to remark that the CROCUS reactor can be basically considered as a homogeneous system, with minimal discrepancies between the α and k formulations.
An application related to the effective kinetics parameters is the estimation of the reactor period and the sub-criticality level of the system. Both quantities can also be measured during experiments by the pulsed neutron source method [1] and reactor noise analysis methods [2,4]. The choice of the optimal adjoint weighting function (ϕ † α or ϕ † k ) in order to compare the results computed from numerical simulations to those obtained from measurements depends on the procedure and the techniques adopted during the experiment. Moreover, mixing α and k weighted kinetics parameters can be considered for the estimation of α 0 and ρ [28].

Funding:
The authors wish to thank Electricité de France (EDF) for partial financial support.