Guillaume Truchet

Pierre Leconte

Alain Santamarina

Emeric Brun

Zohia Andrea

Andrea Zoia
email: andrea.zoia@cea.fr

Computing adjoint-weighted kinetics parameters in Tripoli-4 by the Iterated Fission Probability method

Keywords: Iterated Fission Probability, Kinetics parameters, Monte Carlo, Tripoli-4 R

Introduction

Knowledge of kinetics parameters is key to the study of nuclear reactor dynamics, in connection with transient behavior due to normal operation or accidents (Bell and Glasstone, 1970;Keepin, 1965). Usually, kinetics parameters stem from integrating the full time-dependent Boltzmann and precursors equations over the space, energy and angle variables: after various algebraic manipulations, they take the general form of a ratio of dot products, namely, ψ, Fϕ / ψ, Gϕ , where F and G are operators acting on the neutron flux ϕ, and ψ is a weighting function depending on specific choices associated to the system under analysis (Bell and Glasstone, 1970;Keepin, 1965;[START_REF] Ott | Nuclear Reactor Dynamics[END_REF]Henry, 1958;Cohen, 1958). Essential to reactor control are the mean generation time of a neutron within a reactor, which is expressed as

ψ, 1 v ϕ ψ, Pϕ , (1)
where P is the total fission operator (including prompt and delayed productions) and v is the neutron speed; and the delayed neutron fraction, which is expressed as

ψ, Bϕ ψ, Pϕ , ( 2
)
where B is the fission operator corresponding to delayed neutrons alone (Bell and Glasstone, 1970;Keepin, 1965). Most often, the weighting function ψ is either taken to be ψ = 1, in which case the kinetics parameters are said to be un-weighted, i.e.,

Λ 0 = 1, 1 v ϕ 1, Pϕ , β 0 = 1, Bϕ 1, Pϕ (3)
or ψ = ϕ † (ϕ † denoting the adjoint neutron flux), in which case the kinetics parameters are said to be adjoint-weighted, or effective, namely,

Λ eff = ϕ † , 1 v ϕ ϕ † , Pϕ , β eff = ϕ † , Bϕ ϕ † , Pϕ . (4)
In the context of reactor physics, Monte Carlo simulation is considered as the gold standard for the computation of physical quantities and is traditionally adopted so as to establish reference values for faster, but approximated, deterministic calculations [START_REF] Spanier | Monte Carlo Principles and Neutron Transport Problems[END_REF]Bell and Glasstone, 1970). Monte Carlo codes are intrinsically based on the simulation of forward random walks, with particles flowing from birth (fission) to death (spatial leakage or absorption): as such, they are ideally suited to provide accurate estimates of the direct flux ϕ and more generally of physical observables weighted by ϕ. Unweighted kinetics parameters can be straightforwardly obtained by running standard (forward) criticality simulations [START_REF] Spriggs | Tripoli-4 User guide[END_REF] and following their precise definition in Eq. (3). For instance, the mean generation time Λ 0 can be estimated by recording the neutron flight time from birth to fission, and the mean delayed neutron fraction β 0 can be estimated by recording the fraction of delayed neutrons emitted at fission [START_REF] Spriggs | Tripoli-4 User guide[END_REF]. The so-called removal time

T 0 = 1, 1 v ϕ 1, Lϕ = k eff Λ 0 , ( 5
)
where L is the net disappearance operator and k eff is the effective multiplication factor, can be similarly estimated by computing the neutron flight time from birth to disappearance by absorption or spatial leakage [START_REF] Spriggs | Tripoli-4 User guide[END_REF]Zoia et al., 2014a).

Extending the capabilities of Monte Carlo codes to adjointweighted kinetics parameters demands the adjoint flux ϕ † , which in principle requires the simulation of particles flowing backward from death to birth: in continuous-energy transport codes, this turns out to be a daunting task that is still a subject of active research (Hoogenboom, 2003). To overcome this obstacle, two strategies have been proposed so far: either resorting to perturbative techniques (Verboomen et al., 2006;Nagaya and Mori, 2010), or finding an estimate of ϕ † that can be computed during forward simulations. For this latter, the early attempts based on the so-called Next Fission Probability (NFP) approximation for ϕ † (Nauchi and Kameyama, 2005;Meulekamp and van der Marck, 2006) have been shown to be possibly flawed for large and/or strongly heterogeneous reactors [START_REF] Nagaya | [END_REF]Nauchi and Kameyama, 2010). A major breakthrough has been nonetheless made possible by the rediscovery (Feghni et al., 2007(Feghni et al., , 2008;;Nauchi and Kameyama, 2010;Kiedrowski, 2011b) of the so-called Iterated Fission Probability (IFP) interpretation of the adjoint flux ϕ † , originally formulated at the beginning of the nuclear era [START_REF] Ussachoff | Proceedings of the 1st UN conference of the peaceful uses of atomic energy[END_REF]Weinberg, 1952;Hurwitz, 1964;[START_REF] Soodak | Pile kinetics[END_REF]. As detailed in the following, IFP enables the exact calculation of adjoint-weighted quantities, thus establishing Monte Carlo simulation as a reference tool for the analysis of kinetics parameters (Nauchi and Kameyama, 2010;Kiedrowski, 2011b;[START_REF] Shim | [END_REF]. A number of production codes have integrated or are planning to integrate IFP capabilities: a non-exhaustive list includes MCNP5 (Kiedrowski, 2011a), SCALE (Perfetti, 2012), andSERPENT2 (Leppanen, 2014).

In this paper, we detail the implementation of the IFP method in Tripoli-4 R , the reference Monte Carlo code developed at CEA (Brun et al., 2014;Tripoli-4 Project Team, 2012), for the calculation of Λ eff and β eff in view of a future release 1 . We will first illustrate the details of the adopted algorithm in Sec. 2, and then verify and validate the method on several reactor configurations in Sec. 3. We will finally discuss the impact of nuclear data libraries, cycle length and inter-cycle correlations on IFP in Sec. 4. Conclusions will be drawn in Sec. 5.

IFP and adjoint-weighted kinetics parameters

We will briefly recall the connection between the IFP method and the adjoint flux ϕ † , for the sake of completeness. Let us define I( r, E, Ω) as the average number of descendant neutrons asymptotically produced in a distant generation g by a single neutron injected into the system with phase space coordinates ( r, E, Ω) at the initial time (in other words, I represents the 1 The IFP algorithm has been previously tested in Tripoli-4 R so as to compute exact perturbations for small-sample experiments (Truchet, 2014a) or sodium void effects (Truchet, 2014b asymptotic importance of the injected neutron). By following the definition of I, it is possible to establish an exact balance equation (Nauchi and Kameyama, 2010;[START_REF] Ussachoff | Proceedings of the 1st UN conference of the peaceful uses of atomic energy[END_REF]. Consider a small neutron path d s = rr = Ωds: then, along d s we have

I( r, E, Ω) = p nc I( r + Ωds, E, Ω) + Q( r + Ωds, E, Ω), (6)
where p nc = 1 -Σ t ds is the non-collision probability, Σ t being the total cross section, and Q is the average number of descendants for neutrons having a collision in d s. Denote now by q( r + Ωds, E → E , Ω → Ω ) the average number of neutrons undergoing a collision during the flight ds and entering the phase space element Ω , E , namely,

q( r + Ωds, E → E , Ω → Ω ) = Σ s ( r + Ωds, E → E , Ω → Ω )ds + 1 k eff χ t (E → E ) 4π ν t Σ f ( r + Ωds, E)ds, (7)
where Σ s is the scattering kernel, Σ f is the fission cross section, ν t is the average number of fission neutrons, and χ t is the fission spectrum. The term Q can be then expressed as

Q( r, E, Ω) = ∞ 0 dE 4π 0 dΩ q( r, E → E , Ω → Ω )I( r, E , Ω ). ( 8
)
Dividing ( 6) by ds and taking ds → 0 we obtain

0 = dI( r, E, Ω) ds -Σ t I( r, E, Ω) + Q( r + Ωds, E, Ω) ds (9)
By developing the total derivative along d s, Eq. ( 9) can be rewritten as

0 = Ω • ∇I( r, E, Ω) -Σ t I( r, E, Ω) + ∞ 0 dE 4π 0 dΩ Σ s ( r, E → E , Ω → Ω )I( r, E , Ω ) + ν t Σ f ( r, E) 4πk eff ∞ 0 dE 4π 0 dΩ χ t (E → E )I( r, E , Ω ). ( 10
)
By inspection, the quantity I( r, E, Ω) satisfies the adjoint Boltzmann equation for ϕ † , which implies that I( r, E, Ω) ∝ ϕ † ( r, E, Ω) (Nauchi and Kameyama, 2010;[START_REF] Ussachoff | Proceedings of the 1st UN conference of the peaceful uses of atomic energy[END_REF]. This offers a practical means of estimating ϕ † ( r, E, Ω) by Monte Carlo methods: the adjoint neutron flux due to a neutron born at phase space coordinates r, E, Ω at a given generation g is proportional to the number of descendants of such neutron that will be found in the system at an asymptotic generation g + M, when M is sufficiently large (Nauchi and Kameyama, 2010). This provides the theoretical basis of the IFP method. As a particular case, choosing M = 0 leads to the the NFP approximation for ϕ † (Nauchi and Kameyama, 2010).

The IFP algorithm

The IFP principles recalled above can be practically integrated in production Monte Carlo codes so as to estimate general adjoint-weighted scores (Nauchi and Kameyama, 2010;Kiedrowski, 2011b;[START_REF] Shim | [END_REF]Kiedrowski, 2011a;Perfetti, 2012;Leppanen, 2014), including sensitivity profiles based on first-order perturbations (Kiedrowski, 2011b) and exact perturbations (Truchet, 2014a,b). Here we will exclusively focus on the calculation of kinetics parameters, namely, the effective mean generation time Λ eff , the effective delayed neutron fraction β eff and the so-called Rossi alpha

α Rossi = - β eff Λ eff = - ϕ † , Bϕ ϕ † , 1 v ϕ , (11)
which is a key quantity occurring in reactivity measurements [START_REF] Pfeiffer | GA-A13079, General Atomic[END_REF]Persson et al., 2008;Hansen, 1985;Cao and Lee, 2008;Bell and Glasstone, 1970;Keepin, 1965).

The dot products involved in the definitions of effective kinetics parameters (Eqs. ( 4) and ( 11)) are evaluated by computing 1) the score associated to each fission neutron during the regular forward random walk of the particle at each generation g, and 2) the associated importance at a later generation g + M by the IFP method. A unique identifier i is attributed to each fission neutron at generation g: this identifier will be transferred to all the descendants of the original particle (if any) along the following M so-called latent generations. The span from g to g + M takes the name of IFP cycle, M being the cycle length.

Memory is also kept of whether or not this neutron was issued from a prompt or a delayed fission. During generation g + M, fission neutrons originating from the original ancestor with identifier i (if any) are recorded, and allow thus estimating the asymptotic importance. In particular, the importance is denoted (ifp p ) i if the ancestor neutron was issued from a prompt fission, and (ifp d ) i if the ancestor neutron was issued from a delayed fission. The adjoint-weighted delayed fission score associated to an history can be expressed as

1 k eff ϕ † , Bϕ history = ϕ † d1 w d1 + ϕ † d2 w d2 + ... = i ϕ † di w di = i (ifp d ) i w di w di = i (ifp d ) i . (12)
In Eq. ( 12), the delayed neutron simulation weight w di of each collision i cancels when being multiplied by the related delayed neutron adjoint flux ϕ † di . The adjoint-weighted total fission score for the same history can be similarly evaluated as the sum of prompt and delayed contributions, namely,

1 k eff ϕ † , Pϕ history = i (ifp p ) i + i (ifp d ) i = i (ifp) i . ( 13
)
When averaging over all contributions coming from the ensemble of simulated neutron histories, Eq. ( 12) converges to ϕ † , Bϕ /k eff and Eq. ( 13) to ϕ † , Pϕ /k eff , respectively. Then, the ratio of ( 12) over ( 13) yields β eff . Finally, the adjointweighted total neutron lifetime can be estimated as

ϕ † , 1 v ϕ history = i [(ifp p ) i + (ifp d ) i ]t i , (14)
where t i is the neutron lifetime at collision i, counted from its birth in generation g.

The central issue of the IFP algorithm is the tracking of the 60 neutron identifier from generation g to generation g + M. In Tripoli-4 R we have decided to start a new IFP cycle at each generation g, which gives rise to an overlapping structure of superposed IFP cycles (each having the same length M), as illustrated in Fig. 1

(bottom). A similar approach has been recently

65 proposed for the Monte Carlo code SERPENT (Leppanen, 2014) and shown to yield no significant increase in computer memory requirements with respect to the sequential non-overlapping IFP scheme originally proposed in (Kiedrowski, 2011a) and displayed in Fig. 1 (top). Each fission neutron that begins a 70 new history is attributed a unique identifier that belongs to the current generation g. Just before the assignment, an array keeps record of the former neutron identifier that had been assigned in the previous generation. At each subsequent generation, a new array of this type is created and the previous is deleted. It is 75 then possible to track the identifiers from generation to generation, up to a maximum number specified by the user. Starting a new IFP cycle at each generation improves the statistics of kinetics parameters by a factor √ M (for fixed number of total generations in the criticality calculation), but it also increases 80 inter-cycle correlations. This issue will be briefly addressed in Sec. 4.

Verification and validation tests

In this Section we present extensive verification and validation tests for the IFP algorithm described above.

Infinite system with two energy groups

As a first verification test, we address an homogeneous system of infinite size, with two energy groups v 1 (rapid) and v 2 (thermal) and two delayed families a and b. Assuming, as in (Kiedrowski, 2010), that there is no up-scattering, that prompt neutrons are emitted in group g = 1 alone, and that fissions can be induced only from neutrons colliding in group 0 0.2 0.4 0.6 0.8 g = 2, yields the following system of equations for the scalar flux ϕ

1 1.2 1.4 1.6 1.8 2 -4 -2 0 2 4 k eff (X -X Analytical )/σ X X = Λ eff X = β eff X = α Rossi
Σ r,1 ϕ 1 = 1 k eff (1 -β tot + ξ 1 )ν f,2 Σ f,2 ϕ 2 Σ r,2 ϕ 2 = Σ s,12 ϕ 1 + 1 k eff ξ 2 ν f,2 Σ f,2 ϕ 2 , (15)
where we have set ϕ g = ϕ(v g ) and Σ x,g = Σ x (v g ). Here Σ s,g j = Σ s (v g → v j ) is the differential scattering kernel, Σ r,g = Σ t,g -Σ s,gg the removal cross-section of group g, Σ f,g the fission cross-section of group g, ν f,g the number of neutrons produced by a fission in group g, χ i,g is the delayed neutron spectrum from delayed family i to energy group g, β i the delayed neutron fraction of family i, β tot = β a + β b , and ξ g = χ a,g β a + χ b,g β b . For this simple configuration, the explicit solutions for the forward and adjoint fluxes can be obtained by algebraic manipulations. Hence, the kinetics parameters β eff , Λ eff and α Rossi can be computed analytically as a function of the physical parameters, as described in (Kiedrowski, 2010). More precisely, we obtain

β eff = Σ s,12 Σ r,1 ξ 1 + ξ 2 Σ s,12 Σ r,1 (1 -β tot + ξ 1 ) + ξ 2 , Λ eff = 1 v 1 Σ s,12 Σ r,2 + 1 v 2 Σ s,12 Σ r,2 -ξ 2 ν f,2 Σ f,2 Σ s,12 Σ r,1 (1 -β tot + ξ 1 ) + ξ 2 ν f,2 Σ f,2 Σ s,12 Σ r,2 -ξ 2 ν f,2 Σ f,2 , ( 16
)
whence also

α Rossi = - ( Σ s,12 Σ r,1 ξ 1 + ξ 2 ) ν f,2 Σ f,2 Σ s,12 Σ r,2 -ξ 2 ν f,2 Σ f,2 1 v 1 Σ s,12 Σ r,2 + 1 v 2 Σ s,12 Σ r,2 -ξ 2 ν f,2 Σ f,2 . ( 17
)
These expressions provide a reliable verification benchmark for the IFP algorithm. To this aim, the continuous-energy code Tripoli-4 R has been modified so as to handle multi-group 16) and ( 17), we display 100 the ratio (X -X Analytical )/σ X , where X is β eff , Λ eff or α Rossi , and σ X is the statistical uncertainty of the Monte Carlo calculation. Since the medium is homogeneous and the energy dependence is weak (only two speeds), the convergence of the IFP algorithm is rather fast, and several tests allow concluding that 105 M = 5 latent generations are sufficient to attain the asymptotic convergence for this simple example. A good agreement between the Monte Carlo results and the reference analytical values is found for both Λ eff and β eff for reactor configurations in the region k eff > 0.5 (see Fig. 2). For reactor configura-110 tions pertaining to the region k eff < 0.5 the convergence of the power iteration in the criticality simulation is rather poor, and numerical results are affected by a large statistical uncertainty. The convergence of α Rossi appears to be affected by a slight bias with respect to the exact solution, which progres-115 sively vanishes by increasing the number of simulated neutrons per generation, as shown in Fig. 3. This behavior is due to the

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 -4 -2 0 2 4 k eff (X -X Analytical )/σ X X = Λ eff X = β eff X = α Rossi
g Σ rg ν fg Σ fg χ a→g χ b→g Σ s,g→1 Σ s,g→2 Σ totg 1 1.5 0 0 3/4 1/2 1/2 1/2 2 2 2 24/5 1 1/4 1/2 0 1 3
Table 1: Reference physical parameters adopted in the verification models.

fact that the Monte Carlo code for practical reasons computes the kinetics parameters in Eqs. ( 4) and ( 11) as averages of ratios of fluctuating quantities, instead of ratios of averages of fluctuating quantities: this introduces a systematic bias, whose relevance is expected to decrease as the number of particles per generation increases. Not surprisingly, the statistical bias becomes particularly evident for α Rossi , which is the ratio of the two most fluctuating quantities, namely, the delayed produc-125 tions ϕ † , Bϕ and the lifetime ϕ † , 1 v ϕ . As a matter of fact, it would be possible to obtain α Rossi as a ratio of the statistical series, measured at each generation. The uncertainty would be then provided by the standard deviation of the ratio of a sum, in the form σ( i X i / i Y i ), which can be calculated by using 130 bootstrapping or jackknife methods.

The rod model with two energy groups

The rod model is possibly the simplest example of space-and direction-dependent transport problem [START_REF] Montagnini | [END_REF]): particles move along a line (the rod) and undergo collision events at a rate v g Σ t,g . Because of the geometric constraints, only two directions of flight are allowed, namely forward (Ω = +) and backward (Ω = -); here, we furthermore assume that scattering and fission are isotropic, i.e., directions taken by the particles after a collision are sampled with equal probability, a single fissile isotope is present, there is no upscattering, prompt neutrons are emitted in group g = 1 alone, and fissions can be induced only from neutrons colliding in group g = 2. If we define ϕ(x, +) the angular flux in the positive direction and ϕ(x, -) in the negative direction, respectively, where x is the spatial coordinate, the transport equations for two energy groups become

± ∂ϕ 1 ∂x (x, ±) + Σ t,1 ϕ 1 (x, ±) = Σ s,11 ϕ 1 (x) + 1 k eff (1 -β tot + ξ 1 )ν f,2 Σ f2 ϕ 2 (x), (18)
± ∂ϕ 2 ∂x (x, ±) + Σ t,2 ϕ 2 (x, ±) = Σ s,22 ϕ 2 (x) + Σ s,12 ϕ 1 (x) + 1 k eff ξ 2 ν f2 Σ f,2 ϕ 2 (x), (19)
where ϕ g (x) = [ϕ g (x, +) + ϕ g (x, -)]/2 is the scalar flux integrated over the allowed directions. Explicit solutions for the kinetics parameters of the rod model can be obtained in principle, but the calculations are quite cumbersome. Instead, we have computed β eff , Λ eff and α Rossi numerically by direct integration of the transport equations and the associated adjoint equations. Results for this example are illustrated in Fig. 4 for the physical parameters given in Tab. 1. The size L of the rod is varied in the interval 1 < L < 10, and we compare the reference values to the outcomes of the Monte Carlo simulations. Each reactor configuration is identified by the corresponding value of k eff . Similarly as in the previous example, the IFP algorithm converges to the reference values.

145

The Rossi alpha validation suite

The parameter α Rossi can be directly measured in low-power reactor configurations close to delayed criticality by resorting to neutron noise analysis techniques (Hansen, 1985;Keepin, 1965;Bell and Glasstone, 1970). When available, these mea-150 surements thus provide an invaluable means of validating the IFP algorithm. As a first test, we consider the Rossi alpha validation suite (Kiedrowski et al., 2011), which contains a collection of experimental measurements of α Rossi for a series of reactor mock-up configurations together with the corre-155 sponding benchmark specifications from (OECD NEA, 2010). The Monte Carlo calculations results obtained with the IFP method implemented in MCNP are also reported in the same document (Kiedrowski et al., 2011). For our simulations, we have chosen the ENDF/B-VII.0 nuclear data library in order to be 160 consistent with the choice of (Kiedrowski et al., 2011). The number of neutrons per generations is 40000, and the IFP cycle length is M = 20. The results obtained with Tripoli-4 R are resumed in Tab. 3, where they are compared to the experimental measurements and to those obtained with MCNP by using 165 the same nuclear data library (as reported in (Kiedrowski et al., 2011)). For all tested configurations, the values computed by Tripoli-4 R are in good agreement with the experimental measurements. Where slight discrepancies are observed, the α Rossi computed by Tripoli-4 R is in agreement with the value com-170 puted by MCNP: as pointed out in (Kiedrowski et al., 2011), the experimental uncertainty reported in the Rossi alpha suite does not include physical uncertainties, so that part of the discrepancies between Monte Carlo simulations and experimental data could be attributed to such factors as physical dimen-175 sions and material densities. Moreover, slight inaccuracies in a few isotopes of the nuclear data library (for instance 233 U and the unresolved resonance range of 235 U (Kiedrowski et al., 2011)) might play a role for some of the configurations. In all configurations, the statistical bias on the α Rossi estima-180 tor is small, thanks to the large number of neutrons per generation. For reference, in Fig. 5 we display the statistical bias -(α Rossi +β eff /Λ eff )/(β eff /Λ eff ) for the GODIVA configuration as a function of the number of particles per generation: for this ex-ample, the bias becomes negligible above 5000 simulated particles per generation.

0 2 4 k eff (X -X Analytical )/σ Value X = Λ eff X = β eff X = α Rossi
For some of the reactor configurations of the Rossi alpha validation suite, experimental values of β eff are also available (Paxton, 1981). The comparison between the measured β eff and the results of Tripoli-4 R and MCNP by resorting to ENDF/B-VII.0 nuclear data library are shown in Tab. 4. For reference, in the same table we also display the values of β eff computed with Tripoli-4 R by using the NFP method and the average delayed neutron fraction β 0 . As expected on physical grounds, the results obtained by resorting to the IFP with a large M are in better agreement with experimental data as compared to un-weighted kinetics parameters. The NFP is in fairly good agreement with IFP in most cases, except for reflected configurations (such as FLATTOP), where large deviations are observed: in these cases, space and energy effects become important, and choosing M = 0 yields a poor approximation of the adjoint flux [START_REF] Nagaya | [END_REF]Nauchi and Kameyama, 2010). The un-weighted parameter β 0 typically shows large differences with respect to the IFP β eff .

Once the kinetics parameters α Rossi and β eff have been determined experimentally, the value of the mean neutron generation time Λ eff can be consistently computed as (Paxton, 1981). To complete our analysis, the values of Λ eff are displayed in Tab. 5. The results obtained by Tripoli-4 R by resorting to ENDF/B-VII.0 nuclear data library are in good agreement with the calculated Λ eff and with those obtained by MCNP. For comparison, we also display the values of Λ eff computed with Tripoli-4 R by using the NFP method and the un-weighted average neutron generation time Λ 0 . Similarly as in the case of β eff , the results obtained by resorting to the IFP with a large M are in better agreement with experimental data as compared to un-weighted kinetics parameters.

Λ eff = -β eff /α Rossi
For reflected configurations, large deviations are again observed between NFP and IFP. The un-weighted parameter Λ 0 typically shows large differences with respect to the IFP Λ eff : to support our analysis, we finally display the values of Λ 0 computed by MCNP2 , which are in excellent agreement with those computed by Tripoli-4 R .

The CALIBAN reactor facility at CEA/Valduc

The CALIBAN reactor is a super prompt critical configuration operated by CEA/Valduc that is used for criticality analysis. The detailed reactor description and the corresponding experimental measurements of the α Rossi parameter are given in (Richard, 2013). We have performed Monte Carlo calculations for CALIBAN by using the JEFF3.1.1 nuclear data library with 40000 particles per generation and M = 20. The value of α Rossi estimated by Tripoli-4 R , α Rossi = -65.8 ± 0.8 s -1 , is in good agreement with the experimental measurements, α Rossi = -64 ± 1 s -1 , reported in (Richard, 2013).

The MISTRAL reactor facility at CEA/Cadarache

235

The MISTRAL facility is operated at CEA/Cadarache: MIS-TRAL1 has an enriched UO 2 (3.7 % of 235 U) core with about 750 PWR-type fuel pins; MISTRAL2 has 100 % MOX (7 % Pu) core with about 1600 fuel pins arranged in the same lattice pitch as MISTRAL1. The detailed description of the cores 240 and the corresponding experimental measurements of the β eff parameter are available in [START_REF] Santamarina | [END_REF]. We have performed Monte Carlo calculations for MISTRAL1 and MIS-TRAL2 by using the JEFF3.1.1 nuclear data library with 40000 particles per generation and M = 20. The values of 245 β eff estimated by Tripoli-4 R are coherent with the experimental values reported in [START_REF] Santamarina | [END_REF]) (see Tab. 6). For comparison, the results obtained by using the APOLLO-2 deterministic SHEM-MOC calculation scheme (Santamarina et al., 2013) are also collected in the same table. The ORPHEE facility is located at CEA/Saclay and jointly operated by CEA and CNRS3 . ORPHEE is a pool-type research reactor whose main goal is to produce neutron beams for neutron scattering experiments with a broad wavelength and 255 energy distribution. The core of ORPHEE is very small in size and highly enriched in 235 U, providing a high neutron density, with a heavy water reflector tank. A detailed description of the facility and the related α Rossi measurements performed during the initial rod worth calibration campaign are reported in the 260 nuclear safety report (Rapport de Sûreté, 1981). The measurements yield α Rossi -40 s -1 (the experimental uncertainty was unfortunately not reported in (Rapport de Sûreté, 1981)). Calculations performed with Tripoli-4 R by using JEFF3.1.1 nuclear data library, with 10 4 particles per generation and M = 20,

265

give α Rossi = -41.1 ± 0.3 s -1 (Zoia et al., 2014b), which is in good agreement with the experimental data. As for the delayed neutron fraction, the measurement quoted in the nuclear safety report is β eff = 744 pcm (Rapport de Sûreté, 1981), to be compared with the value β eff = 747.8 ± 3.9 pcm computed by Tripoli-4 R (Zoia et al., 2014b).

Discussion

Sensitivity to nuclear data library

The calculations of the IFP kinetics parameters in the previous Section have been performed by selecting a specific nuclear data library, either ENDF/B-VII.0 or JEFF3.1.1. In order to evaluate the impact of the nuclear data on the obtained results, we have re-run some of the configurations of the Rossi alpha validation suite by resorting to the JEFF3.1.1 library. The comparison with respect to ENDF/B-VII.0 for the α Rossi parameter is reported in Tab. 7. In most cases, the discrepancies on the α Rossi parameter are rather weak, and lie within a few percent. In order to single out the distinct contributions to the observed discrepancies, the cases of β eff and Λ eff are separately shown in Tabs. 8 and 9, respectively.

Choosing the IFP cycle length

The length M of the IFP cycle is an important parameter to be chosen before running the simulation. Typically, one assumes that M 10 up to M 20 latent generations are sufficient for the descendant neutrons to converge to the asymptotic adjoint flux. The value of M is of course problem-dependent: configurations where spectral and/or geometrical effects play a relevant role will generally speaking require a larger value of M to reach the asymptotic convergence of the adjoint flux. In Tripoli-4 R , we have implemented a diagnostic tool that allows the convergence of the IFP method to be displayed as a function of M. This tool benefits from the structure of superposed cycles so as to process the results of the IFP algorithm for different cycle lengths 1, 2, • • • , up to the maximum value M required by the user. An example of such analysis is illustrated in Fig. 6 (top) for the GODIVA configuration from the Rossi alpha validation suite. For this case, we have chosen the rather extreme value M = 100 and we have separately computed β eff and Λ eff as a function of the IFP cycle length M. It is immediately apparent that after approximately 10 latent generations the kinetics parameters converge to their respective asymptotic values. At the same time, the statistics of β eff and Λ eff deteriorates for growing M (for a fixed number of simulated particles), as expected on physical grounds: for very large M, the number of descendant neutrons that are still associated to the identifier of the first generation g might be considerably reduced, since at each generation there exists a finite probability that the current fission chain comes to an end. The same convergence analysis has been performed also on a full PWR reactor core of type N4, representative of the French 1450 MW PWR fleet [START_REF] Leconte | Proceedings of the PHYSOR 2010 405 conference[END_REF].

Results for the delayed neutron fraction β eff and the mean generation time Λ eff are displayed in Fig. 6 (bottom). The convergence pattern is quite similar to that obtained for GODIVA, in spite of the conspicuous difference in size of the two systems. In principle, it would be possible to establish some selec-320 tion criteria for an optimum value of M. For instance, one might identify the optimum cycle length as the smallest value of M whose associated statistics for the kinetics parameters is compatible within 1σ uncertainty with the values obtained for longer cycles. The results corresponding to the application of 325 these criteria are resumed in Tab. 2 for some of the reactor configurations described above. In most cases, it is apparent that the selection criteria are (not surprisingly) strongly problemdependent.

Inter-cycle correlations 330

The power iteration algorithm is intrinsically based on the idea of sampling the source of the fission neutrons for the next generation on the basis of the fission sites determined during the current generation. This in turn induces correlations between (local as well as global) scores computed by averaging the statistics over several generations (which are assumed to be stationary realizations of the same stochastic process). Fluctuations around average values appear because of Monte Carlo simulations being intrinsically based on the transport of a finite number of particles: by increasing the number of simulated histories, it is well known that the fluctuations affecting the computed scores are progressively reduced [START_REF] Lux | Monte Carlo particle transport methods: Neutron and photon calculations[END_REF]Brown, 2005). Another source of noise specific to Monte Carlo criticality calculations is related to fission events, because of two concurrent phenomena: i) fissions induce splitting of the trajectories, which increases the dispersion of the population number within a given volume; ii) the birth of a neutron occurring at a fission site induces correlations between particle positions [START_REF] Lux | Monte Carlo particle transport methods: Neutron and photon calculations[END_REF]Sjenitzer and Hoogenboom, 2011;[START_REF] Dumonteil | [END_REF]. As a consequence, the variance of the obtained Monte Carlo scores might be underestimated with respect to the ideal case of truly independent simulations. The kinetics parameters, which are computed by resorting to the power iteration, suffer from the same problem. Moreover, the fact of using superposed IFP cycles potentially increases the existing inter-generation correlations (nonetheless, the average value is not affected). In particular, we have

σ 2 = var(x i ) + 2 k cov(x i , x i+k ), (20)
where var(x i ) is the usual variance estimated from a sample of n terms, namely,

var(x i ) = n i (x i -x) 2 n (21)
In order to investigate the impact of inter-cycle correlations on the variance of kinetics parameters, N independent Monte Carlo simulations have been performed using different random generator seeds. For each run i ∈ N , n realizations of the kinetics parameters have been sampled, so as to compute a mean x i , a standard deviation σ i and a standard error σi = σ i √ n , where the values σ i are obtained by using Eq. ( 21). The effects of intercycle correlations are mirrored in σi , which can be underestimated (if the realizations are correlated) or overestimated (if the realizations are anti-correlated). By construction, the computed σi are strictly independent, so that they can be combined in order to obtain a global value called the apparent standard error σa , namely,

σa = i σ 2 i /N √ Nn = i σ 2 i /N √ Nn = i σ2 i N . (22)
Since Monte Carlo runs are independent, the standard error can be rigorously estimated as

σ = σ √ N = N i (x i -x) 2 N (23)
Then, the impact of inter-cycle correlations can be roughly assessed by computing the ratio f = σ/ σa as in [START_REF] Mervin | Proceedings of the M&C 2011 Conference[END_REF]. Several tests have been performed on a selection of the reactor configurations from the Rossi alpha validation suite presented above. For each configuration, N = 1000 independent 335 runs were performed, with n = 900 active IFP cycles. The number of neutrons per generation is relatively small (500), in order to enhance the effects of correlations. In Fig. 7, the quantity f is plotted as a function of the IFP cycle length. For both β eff and Λ eff , f shows only slight deviations from unit, which means 340 that the impact of inter-cycle correlations for the examples considered here is rather weak. The relevance of correlations increases for decreasing cycle length M: for small M, f rises up to f 1.1. The case of the N4 PWR has been also analyzed in Fig. 8, with N = 100 and n = 900. For this configuration, 345 the correlations on β eff are weak, whereas those on Λ eff rise up to f 1.6 for low cycle length M: this is almost surely due to the large size of the N4 core, which possibly enhances intercycle correlations due to the distribution of fission source sites. Rather surprisingly, Λ eff is only slightly affected by the number 350 of particles per generation.

Conclusions

In this paper we have described the specific implementation of the IFP algorithm in the reference Monte Carlo code Tripoli-4 R developed at CEA, with the aim of computing the adjoint-355 weighted kinetics parameters β eff and Λ eff . We have adopted a scheme based on superposed IFP cycles, which has shown good performances in terms of memory requirements and computational cost. The proposed algorithm has been first verified on several simplified reactor configurations where exact values 360 of β eff and Λ eff were available: in all tested problems, a good agreement has been found between exact and computed kinetics parameters. Then, the IFP method implemented in Tripoli-4 R has been validated on a selection of reactor configurations, including the Rossi alpha validation suite benchmarks and some 365 experimental reactor mock-ups operated at CEA. The validation has concerned in particular β eff and α Rossi = -β eff /Λ eff , which can be measured experimentally in nuclear reactors operated at low power or during startup. Results obtained by comparing the Tripoli-4 R estimates to benchmark and/or experi-370 mental values have been satisfactory. Special attention has been paid to the impact of the nuclear data libraries on the Monte Carlo estimates. The choice of IFP cycle length and the relevance of inter-cycle correlations due to the use of a superposedcycles scheme have been also examined in detail.

  85

Figure 2 :

 2 Figure 2: Discrepancies between calculated and analytical values of Λ eff , β eff and α Rossi (normalized to respective standard deviations σ). Monte Carlo simulations have been performed with 2000 neutrons per generation and 2700 active IFP cycles. Each point identifies a reactor configuration where the parameters ν f,2 and Σ a,g have been modified around their reference values. The corresponding k eff value for each configuration is used in abscissas.

Figure 3 :

 3 Figure3: Discrepancies between calculated and analytical values of Λ eff , β eff and α Rossi (normalized to respective standard deviations σ). Monte Carlo simulations have been performed with 10000 neutrons per generation and 2700 active IFP cycles. Each point identifies a reactor configuration where the parameters ν f,2 and Σ a,g have been modified around their reference values. The corresponding k eff value for each configuration is used in abscissas.

Figure 4 :

 4 Figure 4: Discrepancies between calculated and analytical values of Λ eff , β eff and α Rossi (normalized to respective standard deviations σ). Monte Carlo simulations have been performed with 10000 neutrons per generation and 2700 active IFP cycles. Each point identifies a reactor configuration where the rod length L has been modified in 1 < L < 10. The corresponding k eff value for each configuration is used in abscissas.

Figure 5 :

 5 Figure 5: The GODIVA configuration of the Rossi alpha validation suite. The statistical bias on the α Rossi estimator as a function of the number of simulated neutrons per generation.

  The ORPHEE reactor facility at CEA/Saclay

Figure 6 :

 6 Figure 6: Top. The kinetics parameters β eff and Λ eff as a function of the IFP cycle length M for the GODIVA configuration of the Rossi alpha validation suite. Error bars represent statistical uncertainty at 1σ. Bottom. The kinetics parameters β eff and Λ eff as a function of the IFP cycle length M for the N4 PWR reactor. Error bars represent statistical uncertainty at 1σ.

Figure 7 :Figure 8 :

 78 Figure 7: Inter-cycle correlations on β eff and Λ eff for some of the reactor configurations in the Rossi alpha validation suite. The factor f = σ/ σa is plotted as a function of the IFP cycle length M. Large deviations from f 1 imply strong inter-cycle correlations.

  Top. IFP scheme with non-overlapping cycles. Bottom. IFP scheme with overlapping cycles. The arrows denote the link between the generation at which the scores are computed and the generation where the asymptotic neutron importance is estimated.

	1	2	3	4	5	6	7	8	9
		Cycle 1			Cycle 2			Cycle 3
			Cycle 3		Cycle 6
	1	2	3	4	5	6	7	8	9
	Cycle 1 Cycle 2	Cycle 4 Cycle 5		Cycle 7
	Figure 1:
	).

Table 2 :

 2 Optimal IFP cycle lengths M.

Preprint submitted to Annals of Nuclear EnergyDecember 17, 2014

Actually, the code MCNP computes the un-weighted neutron removal time T 0 (X-5[START_REF] Carlo | MCNP -A general Monte Carlo N-particle transport code, Version 5, Volume I: Overview and theory[END_REF], from which Λ 0 = T 0 /k eff can be easily deduced.

Laboratoire Léon Brillouin. See the website: http://www-llb.cea.fr/

Acknowledgments

The authors wish to thank AREVA and Électricité de France (EDF) for partial financial support.

The delayed neutron fraction β eff (pcm) for the reactor configurations MISTRAL1 and MISTRAL2. We compare the measured β eff to the values computed by Tripoli-4 R by resorting to the JEFF3.1.1 nuclear data library. For reference, the values computed by the APOLLO-2 deterministic code using the SHEM-MOC calculation scheme with the same nuclear data library are also displayed. Moreover, we show the β eff estimated by the NFP method and the un-weighted average delayed neutron fraction β 0 , as computed by Tripoli-4 R .