Analysis of the neutron noise induced by fuel assembly vibrations

The investigation of neutron noise is key to several applications in nuclear reactor physics, such as the detection of control rod or assembly vibrations and the diagnostic of coolant speed and void fraction. In this paper we will elucidate some aspects of the noise equations in the Fourier domain, for the case of periodic fuel rod vibrations with frequency ω0 in a small symmetrical system in which the perturbation is centrally located. We will in particular focus on the double frequency effect, i.e., the emergence of a noise component at 2ω0 (possibly stronger than the one at the fundamental frequency ω0). Our analysis will be carried out without truncating the noise source at the first order and in the context of a non-perturbative approach (i.e., without resorting to linearization). For this purpose, we will select a simple benchmark configuration that is amenable to accurate reference solutions obtained by solving the exact timedependent transport equations. The analysis carried out in this work suggests that the non-perturbative noise equations are mandatory in order to properly discriminate the possible emergence of double frequency effects in neutron noise, especially in view of comparing simulation results to experimental data.


Introduction
Power reactor noise is defined as the fluctuations of the neutron flux around the steady state, due to temperature or density changes and/or displacements of core components, which in turn induce fluctuations on the material cross sections (Pazsit and Demazière, 2010). Despite usually being an unwanted phe-5 nomenon, such fluctuations carry an information content that can be usefully extracted for the purpose of reactor diagnostics and measurements (Pazsit and Demazière, 2010;Williams, 1974;Thie, 1981). The analysis of the neutron flux variations at the local or global scale within the core can be used e.g. in order to detect anomalies (Behriguer et al., 1977;Fry et al., 1986), as in the case of control rod or assembly vibrations (Pazsit and Analytis, 1980), or infer the coolant speed and void fraction (Kosály, 1980). Traditionally, the equations governing the neutron noise are derived based on a small-perturbation formulation: a perturbation of the Boltzmann operator is assumed to induce a response of the same order for the neutron flux, and the 15 terms containing both the perturbed flux and the perturbed cross sections will be neglected by supposing that they will yield higher-order contributions (Pazsit and Demazière, 2010;Williams, 1974). This approach, dubbed the 'orthodox method' (Pazsit, 1984), allows obtaining simpler equations for the neutron noise, which are finally Fourier-transformed in order to establish the customary lin-20 earized noise equations in the frequency domain (Pazsit and Demazière, 2010;Williams, 1974). Such strategy has led to the development of ingenuous analytical tools, mostly relying on the use of Green's functions in diffusion theory, which have been successfully applied to the analysis of many practical diagnosis and monitoring problems (Pazsit and Analytis, 1980;Jonsson et al., 2012;25 Pazsit and Karlsson, 1996;Yamamoto and Sakamoto, 2019;Pazsit and Dykin, 2018); for a review, see, e.g., (Pazsit and Demazière, 2010;Williams, 1974). In the context of vibrating noise sources, a common approximation is to represent the spatial shape of the perturbation as an idealized delta function (the so-called Feinberg-Galanin-Williams model (Williams, 1970)), which might help in ob-30 taining closed-form solutions for simple one-dimensional configurations.
The orthodox method is known to yield to correct results for the case of absorbers of variable strength (Williams, 1974;Thie, 1981;Kosály, 1980), but generally fails for the case of absorber or fuel rod vibrations (Pazsit, 1977): these considerations were implicit already in the seminal work of (Weinberg and Chionis et al., 2020), as witnessed by the renewed interest stimulated by the CORTEX H2020 project (Demazière et al., 2018). For the case of frequency domain, newly developed transport (as opposed to diffusion) solvers (Rouchon, 2016;Yi et al., 2019) and Monte Carlo methods (Yamamoto, 2013;Rouchon et al., 2017) allow addressing realistic applications at the scale of fuel assem-80 blies or full reactor cores with unprecedented accuracy (Demazière, 2011;Rouchon, 2016;Yamamoto, 2018;Mylonakis et al., 2019;Rouchon et al., 2019). However, one must bear in mind that each numerical tool has specific potentialities and drawbacks, as well as implicit assumptions due to the choices made by the code developers, such as working with the linearized equations, introducing 85 the ε/d approximation, truncating the noise source, and so on. In this respect, the reliability of these codes inherently depends on a deeper understanding of the validity domain of the underlying models: the goal of the present work is to elucidate the behaviour of the noise equations, and in particular the double frequency effect, when the source is not truncated at the first order and in the 90 context of a non-perturbative approach (i.e., without resorting to linearization). For this purpose, we will use a simple benchmark configuration, namely, a single vibrating fuel rod centrally located in a host moderator, which avoids unnecessary complications and yet retains the key physical ingredients. This paper is organized as follows. In Sec. 2 we briefly recall the formal-95 ism of the neutron noise equations and the orthodox linearization. In Sec. 3 we introduce the noise source pertaining to the case of mechanical vibrations, and examine the contributions of the discrete harmonics of the source. In order to illustrate the behaviour of the noise field and the interplay of the different harmonics, in Sec. 4 we present a benchmark based on the rod model, a one-100 dimensional, single-speed neutron transport problem, and discuss our findings within the framework of the linearized noise theory. The solutions corresponding to the non-perturbative approach will be derived and contrasted to those of the linearized equations in Sec. 5. Conclusions will be finally drawn in Sec. 6.

The noise equation(s) and the orthodox linearization
However, the source term does depend on the specific assumptions on the noise model for the system under consideration: the operator δB(ω) contains terms of the kind δΣ α (r, E, ω), expressing the Fourier-transformed cross section perturbation for a given reaction α.
In the following, we will assume that the neutron noise is induced by one or 145 several vibrating interfaces between homogeneous materials. Let us focus on the case where the vibration acts along a single axis, say x, and let us denote by x 0 the unperturbed position of the interface 2 . We define Σ L α (E) and Σ R α (E) the spatially homogeneous cross sections of the regions at the left (x < x 0 ) and the right (x > x 0 ), respectively, of the interface at x 0 . We assume that the vibration is periodic, 150 with frequency ω 0 and period T 0 = 2π/ω 0 , and with amplitude ε smaller than the linear size d of the region concerned by the interface perturbation. The position of the moving interface will be described by Under these assumptions, the general form of the interface perturbation in the time domain is where ∆Σ α (E) = Σ L α (E) − Σ R α (E) and H(z) is the Heaviside function, with H(z) = 1 for z ≥ 0, and H(z) = 0 otherwise. Observe that δΣ α (x, E, t) is non-zero only in the region traversed by the vibrating interface, namely, Before investigating the behaviour of the exact source term in Eq. (14), it is instructive to consider the ε/d approximation (Pazsit, 1977). When the amplitude ε of the perturbation is small as compared to d, Eq. (14) can be expanded in powers of ε, by using the Taylor series δ(z) and δ (z) being the delta distribution and its derivative, respectively. We 165 therefore have whose effect is spatially localized at the stationary position x 0 of the interface. Now, by using and we obtain the Fourier-transformed expansion where we have neglected the contributions associated to ω = 0 and to negative 170 frequencies, for reasons that will be discussed later. In the Fourier domain, the effect of an interface periodically vibrating at a frequency ω 0 is thus mirrored in an infinite sum of contributions corresponding to discrete frequencies, each Figure 1: The spatial shape of the noise source. Comparison between the ε/d approximation (plotted as a dashed line, using the nascent functions to represent the singular spatial behaviour located at the interface) and the exact representation (plotted as solid line). Left: the imaginary part of the component at ω 0 ; right: the real part of the component at 2ω 0 . multiple of the fundamental frequency ω 0 of the forcing function and with decreasing amplitudes.

175
Although many authors have pointed out that the second harmonic might in some cases be relevant (Lucia et al., 1973;Antonopoulos-Domis, 1976;Pazsit, 1977Pazsit, , 1984, the Taylor expansion in Eq. (19) has been truncated at the first order in most of the previous works concerning neutron noise, thus including the contribution of the first harmonic alone. In view of the functional form of 180 Eq. (19), one might suppose that the contribution of the second harmonic at 2ω 0 is ε times smaller than that of the first harmonic at ω 0 . However, expression (19) is composed of a series of increasingly singular spatial functions localized at the interface position x 0 . A common approximation consists in smoothing these contributions over the spatial region of the vibration, which is expedient for a 185 numerical treatment (Vidal-Ferrándiz et al., 2020). For small ε, the delta distribution can be approximated by a suitable nascent delta function F ε (x) of the kind with F(z) being normalized 3 with respect to z. A widely adopted choice is to 3 By a change of variables, it follows that F ε (x) is also normalized with respect to x.
so that the ε term at the denominator cancels out with the prefactor of the firstorder contribution in the Taylor expansion. As for the derivative of the delta distribution in the expansion given in Eq. (19), a convenient choice is the hat function We therefore have so that ε 2 term at the denominator cancels out with the prefactor of the secondorder contribution in the Taylor expansion.
Within the limits of the use of nascent functions to describe the singular spatial behaviour of the ε/d approximation, the first harmonic ω 0 will be thus associated to a negative imaginary contribution of amplitude π/2, spatially flat 200 over the vibrating region, whereas the second harmonic 2ω 0 will be associated to a negative real contribution of amplitude π/4 (i.e., only half of the previous), spatially flat over the vibrating region on each side of the interface but changing sign across x 0 (see Fig. 1 for an illustration). These findings suggest that the effects of the second harmonic cannot be neglected a priori: under certain cir-205 cumstances, basically depending on the spatial shape of the fundamental flux ϕ c within the vibrating region 4 , there might be a subtle competition between the Fourier component at ω 0 and the one at 2ω 0 .

Analysis of the exact noise source
Bearing in mind the behaviour of the ε/d approximation, we can now examine the exact noise source for a periodically vibrating interface. Let us first consider the region x > x 0 , where In this case, the Heaviside function is equal to unity when and zero otherwise, i.e., between time and T 0 /2−τ(x, x 0 ) = π/ω 0 −τ(x, x 0 ). This condition yields a rectangular wave (Rouchon, 2016) whose Fourier transform yields Equation (28) defines an infinite sum of contributions corresponding to discrete frequencies ±nω 0 , n = 0, 1, · · · , multiple of the fundamental frequency ω 0 , in analogy with the findings for the ε/d approximation. The space-dependent Fourier coefficients appearing in Eq. (28) are given by and for n ≥ 1, whose amplitude decreases with increasing n. In particular, the coefficient corresponding to the first harmonic at ω 0 is negative and imaginary: for the second harmonic at 2ω 0 the coefficient is negative and real: More generally, all odd coefficients are imaginary, and all even coefficients are 220 real. The spatially-averaged contributions are given byc R 0 = 2, for n ≥ 2. This yields in particularc The amplitude of the source term corresponding to the second harmonic is about twice as small as the one of the first harmonic, which is again consistent with the 225 ε/d approximation. Observe thatc R 1 is the only non-vanishing spatially-averaged coefficient for odd n.
For the region x < x 0 , we have In this case, the Heaviside function is equal to unity when and zero otherwise, i.e., between time T 0 /2 + τ(x, x 0 ) and T 0 − τ(x, x 0 ). This condition defines a rectangular wave that is shifted by T 0 /2 with respect to the one for the region x > x 0 . By applying the Fourier transform, we have thus where c L n (x, x 0 ) = c R n (x, x 0 ) for n ≥ 1, and c L 0 (x, x 0 ) = c R 0 (x, x 0 ) − 2π. The spatially-230 13 averaged contributions are given bȳ The spatial behaviour of the exact noise source is illustrated in Fig. 1.

Multiple interfaces
So far we have examined the effects of a single vibrating interface. The 235 case of a region of size d vibrating into a host material can be dealt with by considering the linear superposition of the effects of two vibrating interfaces, each located at the boundaries of the region 5 . In this case, the terms ∆Σ α (E) corresponding to the respective boundaries will have opposite signs, so that one might possibly expect interference phenomena due to the simultaneous presence 240 of two vibrating interfaces. The intensity of the interference of the interfaces will depend on the interplay between the vibration amplitude ε and the linear separation d between the boundaries: these effects will be discussed in Sec. 4, with the help of a numerical example.

Impact of the source on the noise field 245
Equations (28) and (38) completely define the Fourier spectrum of the noise source corresponding to a vibrating boundary between two homogeneous regions, for a sinusoidal displacement at a single frequency ω 0 . The system will in principle respond to the vibration at several discrete frequencies nω 0 , despite the impulsion being monochromatic at ω 0 . The noise response for a given frequency ω = nω 0 will depend on the source intensity (i.e., on the Fourier coefficients c R,L n (x, x 0 ) and on the shape of the fundamental flux ϕ c ), and on the noise detector function h D (r, Ω, E), where the integral is defined over a region D of the phase space variables r, Ω, E. The noise source term can be generally written as with the appropriate noise source components δB n associated to the cross section perturbations δΣ α (ω), to be inferred from Eqs. (28) and (38). We thus expect the noise field to be of the form (Sanchez, 2015;Rouchon, 2016) where δϕ n (r, Ω, E) is the solution of the exact Fourier-transformed noise equation (10) corresponding to the noise source component at frequency nω 0 , namely Here we have denoted B k = B(kω 0 ) and δB k = δB(kω 0 ) the noise operator and the perturbation operator, respectively, evaluated at the discrete frequencies kω 0 of the source. This yields an infinite system of coupled equations for the noise components δϕ n , where the noise field at a given discrete frequency nω 0 depends 265 on the behaviour of the noise field for all the other frequencies mω 0 . For the orthodox approximation, Eq. (43) reduces to which is an infinite system of fully decoupled linear equations 6 for the noise components δϕ n . The noise source contributions due to negative frequencies n < 0 can be 270 dealt with by observing that δϕ −n = δϕ † n , where † denotes complex conjugation. Finally, the noise source contribution −B 0 ϕ c , which involves the coefficients c R,L 0 (x, x 0 ) and therefore corresponds to ω = 0 (i.e., the time-average of the perturbation), physically represents the offset due to the fact that the perturbation will in general introduce a non-vanishing reactivity effect in the system. For the 275 linearized equations (13) all the components are decoupled from each other, so that the equation for n = 0 will not influence the others and can be disregarded when solving for n ≥ 1 (Sanchez, 2015;Rouchon, 2016). The effect of −B 0 ϕ c on the non-perturbative Eqs. (43) is more involved and the discussion of this case is deferred to Section 5.

Analysis of the linearized noise equations: a benchmark model
In order to assess the linearized equations (13) combined with the exact model of the noise source (Eqs. (28) and (38)) or with the ε/d approximation (Eq. (19)), we consider the case of neutron noise due to vibrating interfaces in a benchmark configuration based on the rod model, which is possibly the sim-285 plest example of space-and direction-dependent transport problem (Wing, 1962;Montagnini and Pierpaoli, 1971): particles move at constant speed υ along a line (the rod) and undergo collision events with total cross section Σ t (x). Only two directions of flight are allowed, namely forward (Ω = +) and backward (Ω = −); we furthermore assume that the angular distributions for scattering and fission 290 are isotropic. For the sake of simplicity, we will consider a single precursor family.
If we define δϕ + (x) and δϕ − (x) to be the (angular) neutron noise at position x in the positive and negative direction, respectively, the linearized noise equation (13) takes the form of two coupled ordinary first-order differential equations where δϕ(x) = δϕ + (x) + δϕ − (x) is the scalar noise integrated over the directions. We will assume that the viable space is a segment [0, L], for some positive length L, with leakage boundary conditions δϕ + (0, ω) = 0 and δϕ − (L, ω) = 0. In principle, the Green's function of Eq. (45) could be determined exactly, by slightly adapting the analytical methods proposed in (Hoogenboom, 2008), and the solution would then follow by quadrature formulas with respect to the noise source term Q ± (x, ω). In practice, the resulting expressions are rather cumbersome, so that a numerical solution by discretization and matrix inversion has been chosen for this work. The noise source Q ± (x, ω) = −δB(ω)ϕ c is given by where ϕ ± c (x) is the fundamental mode corresponding to the stationary system, with associated eigenvalue k eff , and ϕ c (x) = ϕ + c (x)+ϕ − c (x). The quantity δΣ α (x, ω) is the Fourier-transformed, spatially dependent cross section perturbation for re- action α, whose functional form for a vibrating interface is given in Eqs. (28) and (38). The noise field δϕ ± (x, ω) will be decomposed as where the components δϕ ± n (x) satisfy which can be solved separately for each n ≥ 1 (Sanchez, 2015;Rouchon, 2016).

300
In practice, instead of working with complex variables we will further decompose the operators and the noise into real and imaginary parts, which is the strategy adopted by most numerical solvers for the Fourier-domain linearized noise equations, either deterministic or Monte Carlo (Rouchon, 2016).
To fix the ideas, we will consider a rod composed of three regions: a centrally 305 located fuel pin within the inner boundaries at x = x l and x = x r , and a host moderator in x ≤ x l and x ≥ x r , as illustrated in Fig. 2. This configuration is inspired by the Colibri experimental setup, developed in the framework of the CORTEX project, where one or several fuel rods of the Crocus zero-power pooltype reactor are oscillated at different frequencies for the purpose of neutron 310 noise analysis (Lamirand et al., 2020). The physical parameters for the three regions are the following: for the moderator, we take Σ s = 3 cm −1 and Σ c = 0.2 cm −1 ; for the fuel, we take Σ s = 1 cm −1 , Σ c = 1.2 cm −1 and Σ f = 0.98 cm −1 . The delayed neutron fraction is β = 7 × 10 −3 and the precursor decay constant is λ = 0.08 s −1 . The neutron speed is υ = 2.2 × 10 5 cm/s. For our simulations, 315 we will consider three symmetrical configurations with a decreasing size d for the fuel region, and the average number ν t = ν p + ν d of fission neutrons will be adjusted in order to make the system critical: for the first, we take d = 6 cm, with x l = 2 cm and x r = 8 cm, and ν t = 2.2711; for the second, d = 4 cm, with x l = 3 cm and x r = 7 cm, and ν t = 2.30635; for the third, d = 2 cm, with x l = 4 cm 320 and x r = 6 cm, and ν t = 2.4250. In all cases, the total segment length is L = 10 cm and the number of spatial meshes for numerical integration is 8000. The corresponding fundamental modes ϕ c , normalized to unit, are shown in Fig. 3. 17 As for the perturbation, we will assume that the forcing frequency is ω 0 /2π = 0.1 Hz, with a vibration amplitude ε = 0.5 cm, unless otherwise specified.

A single vibrating interface
We start our analysis by considering the somewhat artificial case of a single vibrating interface located at the right boundary of the fuel region (x = x r ). The main findings obtained from the application of the linearized equations coupled to the exact model of the noise source or to the ε/d approximation are illustrated 330 in Figs. 4-6 for the configuration with d = 4 cm. We have computed the first four components, from ω 0 to 4ω 0 . It turns out that the noise field δϕ(x, ω) is dominated by the components corresponding to the first and the second harmonic: this is due to the shape of the noise source, the first and the second harmonic having the smallest number of nodes, and thus the smallest compensation effects 335 between positive and negative contributions. Thus, the noise field for 3ω 0 and 4ω 0 will not be displayed.
Concerning the noise source, the real part of Q(x, ω) is dominated by the component due to 2ω 0 , whereas the imaginary part of Q(x, ω) is dominated by the component due to ω 0 . The real part of δϕ due to ω 0 has approximately the 340 same amplitude as the one due to 2ω 0 , but opposite sign. On the contrary, the imaginary part of δϕ due to ω 0 is considerably larger than the one due to 2ω 0 . This behaviour is mirrored in the shape of the absolute value of δϕ (Fig. 6): an hypothetical noise detector would show a peak at ω 0 whose intensity would be much stronger than at 2ω 0 at any spatial location within the rod.

345
An interesting observation is that the ε/d approximation yields an accurate estimate for the shape of the noise field, for both the first and the second harmonic, despite some apparent discrepancies concerning the shape of the noise source. Quite surprisingly, this holds true even for our choice of ε, which is somewhat large as compared to the linear size d of the fuel region.

Two vibrating interfaces
We discuss next the more physical case of two interfaces vibrating in phase at each end of the fuel region, each giving rise to the shapes discussed in Sec. 4.1, but with opposite signs because of ∆Σ α , which would mimic the effect of the fuel pin being oscillated into the moderator. The main simulation results obtained 355 from the application of the linearized equations coupled to the exact model of the noise source or to the ε/d approximation are illustrated in Figs. 7-9 for the configuration with d = 6 cm, in Figs. 10-12 for the configuration with d = 4 cm, and in Figs. 13-15 for the configuration with d = 2 cm. The computed noise field for 3ω 0 and 4ω 0 is again negligible with respect to that of the first two harmonics 360 and will not be displayed.
As expected, at each interface the qualitative behaviour of the noise source is similar to the case of the source corresponding to a single interface (see Figs. 7, 10 and 13). The real part of Q(x, ω) is dominated by the component due to 2ω 0 , whereas the imaginary part of Q(x, ω) is dominated by the component due to ω 0 . 365 However, the corresponding noise field δϕ(x, ω) has a distinct shape, because of the linear superposition of the noise sources, which results in non-trivial interference effects. The real part of δϕ(x, ω) is dominated by the component due to 2ω 0 (Figs. 8,11 and 14,left), which is symmetric with respect to the two interfaces and thus leads to a constructive interference; for the same cases, the contribution 370 of ω 0 to the real part is on the contrary suppressed by a destructive interference. Conversely, the imaginary part of δϕ(x, ω) is dominated by the component due to ω 0 , which is anti-symmetric with respect to the two interfaces, and thus leads to a destructive interference (Figs. 8,11 and 14,right). Observe in particular that the imaginary noise component due to ω 0 has now a node at x = L/2, corresponding 375 to the mid-point of the two interfaces. This behaviour is mirrored in the shape of the amplitude of δϕ(x, ω) (Figs. 9, 12 and 15, left): the component due to ω 0 is larger than the one due to 2ω 0 far from the vibrating region; conversely, the component due to 2ω 0 becomes larger than the one due to ω 0 within the vibrating region, where the first harmonic is suppressed. Again, the ε/d approximation 380 yields an accurate estimate for the shape of the noise field, for both the first and the second harmonic.
The linear separation d between x l and x r has a strong impact on the behaviour of the noise field. When reducing the size d of the fuel region by keeping fixed the amplitude ε of the vibration, a striking phenomenon occurs: the real 385 part of δϕ corresponding to the component due to 2ω 0 increases dramatically (Figs. 8, 11 and 14, left) and so does the corresponding amplitude of the noise field (Figs. 9, 12 and 15, left). Since the component due to ω 0 is comparatively much less sensitive to d, the amplitude of the noise field corresponding to the second harmonic may become everywhere larger (and eventually much larger) than 390 the one corresponding to the first harmonic, provided that the separation d of the two interfaces is sufficiently small (see Fig. 15, left, for the configuration with d = 2 cm). According to the linearized equations, an hypothetical noise detector would thus show a peak at 2ω 0 whose intensity would be much stronger than at ω 0 at any location within the rod. The same behaviour is observed when increas-395 ing the size of ε for a given d. These surprising features are consistent with the recent findings related to Monte Carlo simulations of the linearized noise equations mentioned before. In general, we would expect that a linear system under the effect of an external forcing function at a frequency ω 0 should primarily respond at the same frequency, except for special cases where the symmetries of 400 the forcing function and/or the detectors might suppress the response at ω 0 and thus artificially promote the response at 2ω 0 . One then wonders whether the predictions of the linear theory in the Fourier domain are physically sound for the configurations examined here.

Non-perturbative approach to the noise equations
In order to ascertain whether the linearized noise equations in the Fourier domain provide a faithful representation of the system behaviour in the presence of vibrating boundaries, reference solutions are required. For this purpose, we analyze the rod problem using a time-domain solver, the non-perturbative equations in the Fourier domain, and finally an average-corrected version of the 410 non-perturbative equations in which a compensation is included for the possible reactivity effect.

Reference solutions from a time-domain solver
We have developed a time-domain solver for the time-dependent rod model equations, namely, where ϕ ± (x, t) is the angular flux and ϕ( where ϕ(x, 0) = ϕ c (x) is the reactor fundamental mode. The functional behaviour of the time-dependent cross sections for reaction α is given by where ∆Σ α = Σ L α − Σ R α is the difference between the spatially homogeneous cross 420 sections at the left and the right of each interface located at x 0 .
Since the system of equations in (49) is stiff because of the large separation between the typical time scale of the neutrons and the one of the precursors, in order to determine the behaviour of the angular flux ϕ ± (x, t) the precursor equation has been discretized in time by using an explicit Euler scheme, whereas the a domain [0, T ], starting from the equilibrium conditions at t = 0. For our examples, we have taken T = 500 s. Once the scalar neutron field ϕ(x, t) has been determined at a collection of hypothetical detectors located at several spatial lo-430 cations 0 ≤ x j ≤ L within the domain, the neutron noise field is then obtained by taking Finally, the Fourier-transformed neutron noise field is computed by taking the FFT of the expression derived in Eq. (52), over the discrete grid of times 0 ≤ t i ≤ T for which the time-dependent noise field is available, and for each spa-435 tial location x j . We will thus have reference solutions for the noise field within the system, to be compared with the results stemming from the linearized equation (45)  regardless of the use of the exact model or the ε/d model for the noise source: the phase should have a step-wise variation within the vibrating region, and the amplitude should be much smaller. In particular, the abrupt increase of the second harmonic for small d seems to be an artefact induced by the orthodox linearization.

Non-perturbative equations in the Fourier domain
For the benchmark configurations discussed in this work, the reference solutions obtained from the time-domain analysis suggest that the linear approximation in the Fourier domain is appropriate for the first harmonic, but is not capable of reproducing the key properties of the amplitude and phase of the noise field 455 for the second harmonic. In this Section we will explore the impact of the nonperturbative noise equations (43) on the accuracy of the system response in the Fourier domain.
In practice, a numerical solution of Eq. (43) requires a truncation of the infinite sum over m, which can be obtained by applying a cut-off at m = ±M: this leads to a closed system, containing only terms for which an explicit equation is available. This system can be solved by iteration as follows (Sanchez, 2015;Rouchon and Sanchez, 2015;Rouchon, 2016): first, we solve the linearized problem in order to initialize the noise field δϕ (0) n for each frequency component nω 0 ; 465 then, for each iteration j ≥ 1 we compute the effective noise source finally, the noise field is determined by solving This procedure is iterated until appropriate convergence criteria on the norm ||δϕ || are met. At convergence, j → ∞, the noise field δϕ (∞) n formally satisfies a linearized equation where the operator B n is replaced by the modified 470 operator B n + δB 0 /2π and the source −δB n ϕ c is replaced by the effective source Q (∞) n (ω). The term δB 0 is kept on the left hand side of the equation: numerical investigations have shown that such operator shift can contribute to the stability of the iterations (Mancusi and Zoia, 2018). Since the equations are coupled, the components at each frequency nω 0 will also depend on δϕ 0 , which is related to 475 the reactivity offset and was dropped in the linearized equations: this term will need a distinct treatment, as shown in the following.

The rod model revisited
Based on the analysis of the linear system, it appears that the two major contributors to the behaviour of the noise field are ω 0 and 2ω 0 , the components 480 corresponding to n ≥ 3 being small for the configurations considered in this work. We can thus safely assume that the dominant portion of the corrections due to the convolution term in Eq. (43) will be carried by δϕ 1 and δϕ 2 , and choose the cut-off M = 4. For the time being, we also neglect the component at ω = 0 by artificially setting δϕ 0 = 0. 485 We revisit then the rod model with two vibrating interfaces. The convergence of the non-perturbative corrections is achieved in about 20 iterations for d = 6 cm, 30 iterations for d = 4 cm and 70 iterations for d = 2 cm. The numerical solutions for the non-perturbative noise equations are compared to those of the linearized equations in Figs. 16-18 for the configuration with d = 6, Figs. 19-21 for the cm. Concerning the effective noise source, the component at 2ω 0 is slightly but systematically affected by the additional terms given by the convolution product, whereas the component at ω 0 is almost unaffected (Figs. 16, 19 and 22).
Despite these modifications of the effective source being rather small, the 495 corresponding impact on the noise field is extremely strong, which shows that the symmetry-induced compensations occurring in linearized equations are broken by the presence of the coupling terms of the non-perturbative equations. In particular, the real part of the second harmonic is dramatically decreased with respect to the linearized noise equations (Figs. 17, 20 and 23, left), which in turn 500 is mirrored by a decrease of the amplitude of this harmonic (Figs. 18, 21 and24, left). Conversely, the component due to ω 0 is minimally modified by the nonperturbative terms. Eventually, within the framework of the non-perturbative noise equations the first harmonic dominates everywhere, except for the vibrating region, where the absolute value of δϕ 1 (x) is depressed because of the anti-505 symmetry of the imaginary part of the component due to ω 0 . The noise fields obtained by using the non-perturbative equations have been compared to the reference solutions in Figs. 18, 21 and 24. For the first harmonic, the corrections due to the non-perturbative terms are rather small, and the overall agreement with respect to the reference solutions are good, except for the case of the smallest fuel region (d = 2 cm), where a slight but systematic discrepancy is found far from the vibrating region. On the contrary, the impact of the non-perturbative terms on the second harmonic is very strong: for all the tested configurations, the non-linearized noise equations perform better than the linearized equations in reproducing the spatial behaviour of the second harmonic, 515 for both amplitude and phase (compare Figs. 18,21 and 24 against Figs. 9,12 and 15,respectively). However, for decreasing d the discrepancy with respect to the reference solution tends to increase. This is particularly apparent for the configuration with d = 2 cm: while for the other configurations the non-perturbative equations were able to correctly reproduce the phase jump of the second har-520 monic within the vibrating region (which the linearized equations failed to do), for this case the predicted phase is spatially flat, and the shape of the amplitude is correspondingly far from the reference curve. The origin of such discrepancies will be investigated in the next section.

525
In the customary derivation of the noise equations, the neutron noise is defined as the time-dependent deviation with respect to the steady-state flux ϕ c satisfying Eq. (1). This hypothesis underlies the entire formalism for both the linearized and the non-perturbative noise equations in the Fourier domain. Actually, 28 some recent works have suggested that the traditional decomposition of the timedependent neutron flux as in Eq. (4) might not be the wisest choice (Sanchez, 2015;Rouchon and Sanchez, 2015;Rouchon, 2016). The neutron flux ϕ c corresponds to the initial condition of the system, before the perturbation is introduced; once the perturbation is effective, the net reactivity induced by the perturbation will be compensated by control rod adjustments, physical feedbacks 535 such as Doppler effect, or any other external control mechanisms, which will prevent the system from drifting, and the reactor will asymptotically reach a new steady stateφ c ϕ c . In order to mimic these effects, we will simply assume that the average number of fission neutrons is rescaled byk eff (Sanchez, 2015), wherek eff k eff is the new fundamental eigenvalue associated to the steady state 540 including the perturbation (Sanchez, 2015;Rouchon and Sanchez, 2015;Rouchon, 2016). Correspondingly, we will define the neutron noise δφ aŝ which ensures that the noise will have a time-average equal to zero, i.e., where we have used the short-hand notation · for the time average. This proce-545 dure allows rigorously having δφ 0 = δφ(ω = 0) = 0, which minimizes the possible bias in the non-perturbative equations due to having artificially set δϕ 0 = 0 in the previous section. For this purpose, we decompose the perturbed Boltzmann operator as and we set and where A is a suitable time-independent operator, to be determined. We further requireφ c to satisfy the eigenvalue equationB cφc = 0, whereB c is the stationary operator associated toB(t). The resulting noise equation in the time domain reads 29 which can be then Fourier-transformed to yield for the noise components at discrete frequencies nω 0 . Then, a sufficient (albeit non-unique) condition to ensure that δφ 0 = 0 is Thus, we formally set the operator A to be (Sanchez, 2015) A Since A depends onφ c and δφ n , which are the solutions of the problem that 560 we are attempting to solve, an iterative strategy must be used in order to determine the noise field corresponding to the average-corrected non-perturbative noise equations (Rouchon and Sanchez, 2015;Rouchon, 2016). We first solve the non-perturbative noise equations without taking into account the average correction, as detailed in Sec. 5.2, and estimate δϕ n . Then, we determine A based on 565 Eq. (64) and compute the modified operatorsB n and δB n in the Fourier domain. This allows computing the new steady-state fluxφ c and successively updating the noise field δφ n . By iterating these steps, this procedure eventually converges towards the average-corrected δφ n satisfying δφ 0 = 0. We have numerically assessed the impact of using the average-corrected

Conclusions
In this work we have analyzed the behaviour of the neutron noise equations for a simple benchmark configuration consisting in a fuel rod periodically vibrating at the fundamental frequency ω 0 in a host moderator. In this framework, we have elucidated how the role of the noise source, and in particular the impact of the higher-order terms, is entangled with the structure of the noise equations. We have shown that the non-perturbative noise equations (as opposed to the widely used linearized noise equations) are in some case mandatory in or-600 der for the Fourier-domain analysis to be consistent with the reference solutions obtained by solving the time-dependent transport equations. A further improvement comes from the use of the average-corrected noise equations, which can take into account the reactivity offset due to the presence of the perturbation. The discrepancy between the predictions of these three families of noise equa-605 tions (linearized, non-perturbative and average-corrected) depends on the noise component δϕ n to be determined and on the system under investigation. For the benchmark selected in this work, the three equations provide consistent estimates for the spatial shape of the first harmonic δϕ 1 at ω 0 ; however, the differences between the three approaches are dramatic for the second harmonic δϕ 2 at 2ω 0 : the 610 linearized equation might hint that the second harmonic is everywhere larger than the first, which contradicts the reference solution; the non-perturbative equations allow obtaining closer, albeit still not entirely accurate, estimates; the averagecorrected equations finally yield spatial shapes in very good agreement with the reference solution.

615
The analysis carried out in this work seems thus to suggest that the nonperturbative, and even better the average-corrected, noise equations are mandatory in order to properly discriminate the double frequency effect of neutron noise possibly measured in detectors, especially in view of comparing the simulation results to experimental measurements. However, these conclusions hold 620 true for the configurations investigated here, namely a small symmetrical system in which the perturbation is centrally located, and their more broad validity for complex systems (three-dimensional geometries, realistic material compositions and cross sections, and multi-group or continuous-energy transport) should be carefully probed. In this context, an important issue concerns the role of 625 symmetries, as apparent from the following example. For the configuration examined here, the growth of the second harmonic with respect to the first for decreasing d appears to be an artefact of the linearized equations, and disappears when introducing the non-perturbative corrections. However, if we had taken a spatially-integrated detector symmetric with respect to the center of the rod, the dominant noise component would have been at 2ω 0 , independently of whether we had used the linearized or non-perturbative equations: since the noise field related to the first harmonic is anti-symmetric with respect to the center of the system, its integral over a symmetric detector is vanishing and the second harmonic would be the only surviving noise component. Similarly, the impact of the non-635 perturbative corrections on the first and second harmonic, which are intimately related to the behaviour of the fundamental neutron flux within the vibrating region, strongly depends on the symmetries of the system. Future work will be aimed at exploring more general cases, in order to further unveil the features the noise equations.