Saad Yalouz

Vincent Pouthier
email: vincent.pouthier@univ-fcomte.fr

Continuous time quantum walk on an extended star graph: disorder enhanced trapping process

Keywords: numbers: 64.60.aq, 63.22.-m, 71.35.-y, 71.38.-k

Using a tight binding model, we investigate the dynamics of an exciton on a disordered extended star graph whose central site acts as an energetic trap. Our results highlight the beneficial role of the disorder to improve the excitonic absorption at the core of the graph. In this context, we evidence the existence of an optimal disorder allowing to strongly minimize the absorption time. Numerical proofs are given to demonstrate that this enhancement originates from a disorder-induced restructuring process that positively affects the system eigenstates. In a complementary way, we also evidence the existence of an optimal value of absorption rate allowing to reduce even more the aborption time when the disorder is optimal. The resulting super optimized trapping process is interpreted as the positive interplay between the action of the disorder and the arising of the so-called superradiance transition.

I. INTRODUCTION

Exciton trapping is a key issue for explaining many phenomena in molecular lattices [1][2][3][4][5] . Examples among many are the time evolution of fluorescence processes 6 and the excitonic transport efficiency 7 , to cite but a few. Although most seminal works were initially devoted to translational invariant networks [8][9][10] , recent studies have been conducted to investigate excitonic trapping in more complex molecular networks. In this context, two kinds of systems have been considered simultaneously.

First, exciton trapping has been studied in realistic molecular networks including the Fenna-Matthews-Olson (FMO) complex 11 and polymer structures such as dendrimers 12 . The FMO complex, encountered in green sulfur bacteria, is a pigment-protein complex that involves seven bacteriochlorophyll-a molecules. In this protein, the capture of light generates excitons that converge toward a reaction center acting as an irreversible trap. On this specific site, the trapping of an exciton initiates the photochemistry that ultimately leads to the chemical storage of energy [13][14][15][16] . Similarly, a dendrimer whose terminal groups are functionalized by chromophores, can play the role of an artificial light-harvesting complex [17][18][19][20][21][22] . In this context, the absorption of photons by the peripheral chomophores generates excitons that can propagate through the branches of the molecule to reach the core site which contains either a fluorescent trap, a reaction center, or a chemical sensor 23,24 .

Second, from a more formal point of view, excitonic trapping process has also been studied in various complex networks including for example hyperbranched fractals 25 , Sierpinsky fractals 26 , cycle graphs with longrange interactions 27 , chains and rings [28][29][30] , and random networks 31 . These works, intimately connected to graph theory, exploit the formal resemblance between the excitonic delocalization occuring in nature and the continuous time quantum walk (CTQW) 32 . As the quantum analog of classical random walk, CTQW on complex networks has been extensively studied during the past few years due to its various applications ranging from the research of a perfect quantum transport 33 to the development of high-performance algorithms 34,35 .

In that context, to judge the efficiency of the excitonic trapping process on a network, the fundamental question arises whether the exciton can delocalize coherently between several sites, or conversely, if the latter localizes. Although this problem still remains open, recent works suggest that the localized or delocalized nature of the exciton depends on general features such as the excitontrap coupling strength, the network symmetry and the presence of disorder 32 .

First, the coupling with traps has a very pronounced impact on the efficiency of the excitonic propagation [36][37][38] . Resulting from the interaction with an external continuum, the trapping effect is usually taken into account using a non-hermitian exciton Hamiltonian with complex eigenvalues. The real parts of these eigenvalues define the excitonic energies, whereas the imaginary parts specify the energy widths, i.e. decay rates. In that case, a detailed study of these decay rates has revealed the occurrence of a general phenomenon called superradiance transition (ST) 39 . When the exciton is weakly coupled to traps on a network, all the excitonic eigenstates get similarly affected and exhibit quite similar decay rates. However, as the exciton-traps coupling increases, an eigenstate restructuring occurs. Only a few short-lived states, called superradiant states, exhibit cooperatively enhanced decay rates. These states are accompanied by subradiant eigenstates, that is, long-lived states almost decoupled from the traps. Superradiant states tend to localize on the traps, thus preventing the occurrence of an efficient trapping process.

In a different context, as shown by Mulken et al. 40 , localization processes may also result from the degeneracy of the excitonic spectrum that originates in the symmetry of the network. In this case, when the excitonic quantum state initially expands over few highly degenerate eigenstates, specific quantum self-interferences arise. The propagation of the exciton is thus stopped and the latter remains confined in the neighbourhood of the excited region of the network. Note that this feature has been observed in many types of networks including compact dendrimers 41 , star graphs 42 and Apollonian networks 43 .

Finally, as evidenced in solid-state physics, localization processes may also arise on networks that exhibit random defects. Nevertheless, in this case the influence of the disorder is quite subtle and different situations can occur. On the one hand, the disorder can acts as a negative ingredient that completely hinders the excitonic propagation. This feature has been observed for example in linear chains 44 , discrete rings 45 and binary trees 46 , to cite but a few. On the other hand, the disorder can also behave as a positive ingredient that enhances the exciton delocalization. For instance, in tree graphs and dendrimers, the presence of a weak disorder yields extended states through fluctuation-enabled resonances between states that initially may appear to be localized [49][50][51] . Similarly on a Watts-Strogatz network with small-world behavior 52 , the efficiency of the excitonic propagation is enhanced when the rewiring probability turns on 53 . Finally, it has been shown that the addition of random bonds to a star graph allows the exciton wave function to spread more 54 .

In light of all these features, the goal of the present paper is simple : we aim here at characterizing how the interplay between the excitonic spectral degeneracy, the exciton-trap coupling, and the presence of disorder can affect the evolution of the trapping process occuring on a complex network.

To do so, as a logical following of our precedent paper 57 , we consider in this work a network that has been extensively characterized in the absence of disorder : the extended stargraph with a central core occupied by a trap. Organized around a central core, this graph exhibits the local tree symmetry of irregular and complex networks. In this context, our previous study has shown that the exciton dynamics on this graph is governed by two kinds of eigenstates, namely many eigenstates associated to highly degenerate real eigenenergies (insensitive to the trap action) and three decaying eigenstates characterized by complex energies (sensitive to the trap action). As a result, when the exciton is initially located on a peripheral site of the graph, degeneracy induced localization favors the confinement of an important part of the exciton population in the neighborhood of the excited site. Only, a small part of the population is then absorbed at the core of the network. Nevertheless, it has been shown that when the size of the network and the strength of the trap are judiciously chosen, the efficiency of the transfer is optimized at ST, resulting in the minimization of the absorption time.

All these results being obtained on the extended star graph with no disorder, we address here the fundamen-tal question whether disorder-induced symmetry breaking modifies the evolution of the excitonic propagation and affects the trapping process.

Therefore, the present paper is organized as follows. In Sec. II, the disordered extended star graph with a central trap is introduced and the exciton hamiltonian defined. Then, the relevant ingredients required for characterizing the dynamics are described. The problem is then solved numerically in Sec. III, where a detailed analysis of the excitonic absorption is performed. Special attention is paid for characterizing the way the disorder affects the absorption process. The numerical results are finally discussed and interpreted in Sec. IV.

II. THEORETICAL BACKGROUND A. Model Hamiltonian

As displayed in Fig. 1, the network we consider is an extended star graph representing a two-generation tree-like structure built around a central core. Emanating from this central core, N 1 branches do the connexion to the central nodes of N 1 peripheral stars. Each peripheral star involves N 2 external branches so that the total number of sites is

N = 1 + N 1 + N 1 N 2 .
To give a precise description of all of the sites, we introduce two indexes ( , s) such as : ( = 0, s = 0) stands for the central core, ( , s = 0), with = 1, ..., N 1 , refers to the central node of the th peripheral star, and ( , s), with = 1, ..., N 1 and s = 1, ..., N 2 , characterizes the sth peripheral site of the th star.

On this network, the exciton dynamics is modelized as follows. First, we consider that each site ( , s) is occupied by a molecular subunit whose internal dynamics is described by a two-level system with a Bohr frequency ω ,s . Let | , s defines the state in which the ( , s)th two-level system occupies its first excited state, all the other twolevel systems remaining in their ground state. Then, we assume that a trap is located on the central core of the graph. This trap is responsible for an irreversible decay of the excitonic energy according to the decay rate Γ. Within this model, the exciton dynamics is governed by the non-hermitian effective Hamiltonian 55,56 defined as (with ≡ 1)

H = ω 0,0 -i Γ 2 |0, 0 0, 0| + N1 =1 N2 s=0 ω ,s | , s , s| + N1 =1 Φ (|0, 0 , 0| + | , 0 0, 0|) + N1 =1 N2 s=1 Φ(| , 0 , s| + | , s , 0|), ( 1
)
where Φ is the exciton hopping constant. According to the standard Anderson model 58,59 , the disorder is introduced by assuming that the site energies {ω ,s } are independent random variables uniformly distributed in the interval [ω 0 -W/2, ω 0 + W/2]. Here, ω 0 represents the mean value of each site energy. The amplitude W represents the strength of the disorder. It also defines the variance W 2 /12 of each site energy.

Due to its non-hermiticity, the Hamiltonian H exhibits a set of N complex eigenvalues ωk = ω k -iγ k /2 (with k = 1 . . . N ) where the real parts ω k define effective energies whereas the imaginary parts γ k represent decay rates (i.e. widths). Each eigenvalue is associated to a corresponding couple of bi-orthogonal left/right eigenvectors called | Ψk and |Ψ k respectively. Note that, since H is symmetric, one has the property , s|Ψ k = Ψk | , s . Consequently, every corresponding left/right eigenvectors share a similar weight on each site of the graph i.e. | , s|Ψ k | 2 = | , s| Ψk | 2 . Also note that in the following, we will consider that each left/right eigenvector is normalized such as

,s | , s|Ψ k | 2 = ,s | , s| Ψk | 2 ≡ 1.
Considering the bi-orthogonality of the eigenvectors, a completeness relation can be introduced

k |Ψ k Ψk |/ Ψk |Ψ k = 1
, where 1 is the identity operator. Using this relation, the excitonic Hamiltonian H can be written in a diagonal form, as

H = N k=1 ωk |Ψ k Ψk | Ψk |Ψ k . (2)
B. Quantum dynamics Following our previous work 57 , the aim of the present study is to characterize the excitonic absorption process at the core of a disordered extended star graph assuming that an exciton is initially located on a peripheral site ( 0 = 1, s 0 = 1).

To study the effects of the disorder, we use a numerical approach which consists in estimating the mean behavior of the system properties for a fixed couple of parameters (Γ, W ). To better introduce this approach, let us consider that each disordered configuration (i.e. the set of random sites energies on the graph) can be labeled by a specific index (c). In this context, one numerically generates a serie of N C disordered configurations. For each configuration (c), the associated hamiltonian matrix H (c) is built and diagonalized. One thus obtains the complex eigenenergies ω(c)

k and eigenstates |Ψ (c) k (and | Ψ(c) k ).
Then, the dynamic is simulated by numerically building the evolution operator U (c) (t) ≡ exp(-iH (c) t) which allows to compute the excitonic state |Ψ (c) (t) for each time t as

|Ψ (c) (t) = U (c) (t)| 0 , s 0 . (3)
By following this method, it becomes possible to estimate any observable necessary for characterizing the absorption process. To proceed, one assumes that, for each disordered configuration (c), any generic observable O takes a specific value O (c) . In this context, two complementary ways will be followed. Most of the time, a statistical computation will be realized to estimate the average value of the observable O as a function of the model parameters. Nevertheless, in some situations a special attention will be paid to particular values O (c) obtained for specific configurations of the disorder in order to highlight interesting features.

III. NUMERICAL RESULTS

In this section, the previous formalism is applied for describing the excitonic absorption process on a disordered extended star graph with fixed parameters N 1 = N 2 = 4. Note that, in all this work the hopping constant Φ will be considered as the energy unit.

A. Absorption process

First, let us start by evidencing a fundamental property : on the extended star graph, the presence of disorder globally improves the efficiency of the excitonic absorption process. To observe this feature, let us introduce the absorbed population P (c) A (t) defined as

P (c) A (t) = 1 - ,s | , s|Ψ (c) (t) | 2 . ( 4
)
This quantity represents the probability for the exciton to be absorbed by the trap at time t (for a given disordered configuration (c)). Based on this definition, the time evolution of the absorbed population P (c) A (t) has been computed for particular disordered configurations. Fig. 2 displays the obtained results. In order to clearly observe the differences of behavior with/without disorder, the particular case of an ordered graph is also shown for (Γ = 1Φ, W = 0). When W = 0, the absorbed population increases over a short time ∼ 10Φ -1 to reach a maximum value equal to P max A = 6.25.10 -2 . In this case, the absorption process is not complete.

If now one considers the presence of disorder W > 0, a very contrasted behavior arises. In this case, whatever the strength of the disorder, the absorbed population always converges toward unity. Therefore, the presence of defects enhances the absorption process : the whole excitonic population is fully trapped at the core of the graph. Note that this new behavior is a general phenomenon that has been observed in all the simulations we have carried out.

Considering this new property, Fig. 2 clearly shows that the time needed to reach the complete regime of absorption strongly depends on the values of W and Γ. To characterize this feature, let us introduce the absorption time τ (c) . For each disordered configuration (c), τ (c) represents the time for which half of the total excitonic population is absorbed

P (c) A (τ (c) ) = 1 2 . ( 5
)
Based on this definition, let us focus on the different time evolutions shown in Fig. 2. For a weak disorder (Γ = 1Φ, W = 0.1Φ), the absorption process exhibits two phases. During the first phase, a population near from 6.25.10 -2 is fastly absorbed over a short time ∼ 10Φ -1 . Then, over a longer timescale, the absorbed population progressively increases to finally reach unity. The resulting absorption time is τ (c) ≈ 7.10 3 Φ -1 . In this situation, the absorption process is complete but not really efficient due to the long time needed to trap the exciton. Nevertheless, the absorption efficiency can be surprisingly improved if one considers a stronger disorder. Indeed, with (Γ = 1Φ, W = 1Φ), the absorbed population directly varies from zero to unity resulting in a shorter absorption time τ (c) ≈ 100Φ -1 . A more marked enhancement is present when one considers the last configuration with (Γ = 5Φ, W = 5Φ). In this case, the absorbed process is faster and the absorption time reduces to τ (c) 14Φ -1 . All these observations have motivated deeper numerical investigations to better understand the influence of W and Γ on the absorption time τ (c) . To proceed, the mean absorption time τ has been numerically estimated. This characteristic time is defined as the geometric mean over the values τ (c) as

τ = exp 1 N C N C c=1 ln(τ (c) ) . (6)
The reason why a geometric mean has been introduced to replace the standard arithmetic mean is simple to understand. When carrying out a statistic over the different disordered configurations, it turns out that the distribution of the absorption time exhibits a very long tail. In this context, the use of the standard arithmetic mean produces irrelevant results : the average value is not representative of the true central tendancy of the series. Such an inaccuracy is related to the great sensitivity of the arithmetic mean to the few particular values τ (c) that are very large. To overcome this problem, the solution consists in using another definition for the statistical mean. With long trail distributions, the geometric mean represents a very good option due to the lower sensitivity of this estimator regarding the large singular values of the series.

Based on this definition, the mean absorption time τ has been evaluated for different values of the parameters W and Γ. Fig. 3(a) displays the evolution of τ in the two-dimensional (Γ, W ) parameter space. The blue curves delimit the isochronal contours of τ whereas the black hatched area delimits the zone where the mean absorption time always exceeds the limit τ > 1000Φ -1 . To complete this representation, Fig. 3(b) shows three particular sectional views of τ in this space for different fixed values of Γ.

In figure 3(a), it is straightforward to see the presence of an optimal region of absorption where τ is strongly reduced. This optimal region is centered around the minimal value of the absorption time τ opt ≈ 50Φ -1 arising when the model parameters satisfy W opt ≈ Γ opt ≈ 6Φ. When one moves away from this optimal region, τ increases. However, such a behavior is not isotropic in the parameter space. Indeed, as illustrated by the isochronal curves, the contours correspond to cigar-shaped isolines more elongated along the Γ direction rather than along the W direction. For instance, along the W direction, the isochronal contour τ = 80Φ -1 extends over 10Φ and it varies from W = 3Φ to W = 13Φ. By contrast, along the Γ direction, the same isochronal contour extends over 20.5Φ and it varies from Γ = 1.5Φ to Γ = 22Φ. The mean absorption time is then typically twice as sensitive to the strength of the disorder as to the absorption rate.

Actually, the anisotropic behavior of τ hides a more general optimization phenomenon that is highlighted in Fig. 3(b). In this figure, the W dependance of τ is displayed for three fixed values of Γ. Comparing these three different curves, it clearly appears that, whatever the value of Γ, a minimum of τ is always reached when W ∼ W opt . In the case Γ = 1Φ, the minimum mean absorption time is τ ≈ 100Φ -1 . If now the absorption parameter is optimal Γ = Γ opt = 6Φ, the minimum mean absorption time is super-optimized and decreases to τ ≈ 50Φ -1 . However, for a larger absorption parameter Γ = 15Φ, the mean absorption time admits a larger minimum value that is τ ≈ 60Φ -1 .

B. Study of the active states

In order to better understand the evolution of the mean absorption time τ , the properties of the system eigenstates have been numerically studied. By doing so, an important feature has been observed : whatever the disordered configuration, only few eigenstates govern the excitonic dynamics. Due to their active partipation in the trapping process, these particular eigenstates have been renamed "active states".

In this context, it becomes clear that knowing the properties of these particular states could bring a key information to interpret the evolution of the absorption phenomenon. Consequently, we have proceeded to a deeper numerical characterization of the active states. To do so, we took advantage of the fact that these states always exhibited a consequent weight | 0 , s 0 |Ψ (c) k | 2 on the site ( 0 = 1, s 0 = 1) that is initially carrying the exciton. Based on this property, we have then introduced a simple but efficient method to detect them numerically : for each configuration (c), an eigenstate |Ψ (c) k was considered as active only if its weight on the initial peripheral site carrying the exciton respected the minimal condition | 0 , s 0 |Ψ (c) k | 2 ≥ 10%. Note that this choice of minimal weight of detection has been directly inspired by numerical observations. Indeed, in our study it turned out that each state owning a weight of only a few percents on the initial excited site never played a key role in the excitonic absorption and, consequently, behaved as a "non-active state". In this context, several tests revealed that using a threshold of 10% was a really efficient way to left apart all these states.

Proceeding with this numerical approach, a complete study of the active states properties has been realized. By doing so, an important thing has been noticed : the active states are always more sensitive to the strength of the disorder W than to the absorption rate Γ. This phenomenon clearly reminds the behavior of the mean absorption time τ , as seen previously. Consequently, we will here focus on the effect of the disorder for three situations Γ = 1Φ, Γ = Γ opt = 6Φ and Γ = 15Φ. In all the following figures, the curves with black circles correspond to results obtained for Γ = 1Φ, green triangles for Γ = Γ opt = 6Φ and blue squares for Γ = 15Φ (See color online).

To characterize the active states, the first property we have focused on is their mean number in the system. Their mean number N AS is defined as

N AS = 1 N C N C c=1 N (c) AS , (7) where N (c)
AS represents the total number of active states detected in the cth disordered configuration.

The W dependence of N AS is illustrated in Fig. 4. By comparing the three curves, we clearly see that, whatever Γ, a quite similar behavior is adopted as W increases. For small W values, the number of active states is the same, i.e. N AS ≈ 1.7. Then as W increases, N AS increases until it reaches a maximum value N AS max . For Γ = 1Φ and Γ = 15Φ, N AS max ≈ 2.6, whereas for Γ = Γ opt = 6Φ, N AS max ≈ 2.4. Note that this maximization phenomenon always arises in the neighbourhood of the optimal disorder W = W opt . Then, when W exceeds W opt the three curves decrease. They tend to a steady value N AS = 1 as one approaches the strong disorder limit W W opt . The next property we have focused on is the mean inverse participation ratio IP R AS of the active states. The inverse participation ratio gives a measure of the spatial extension of an eigenstate by counting the number of sites it occupies. For a particular eigenstate |Ψ

(c) k , it is defined as IP R(|Ψ (c) k ) = 1 ,s | , s|Ψ (c) k | 4 . ( 8
)
With this definition, an eigenstate |Ψ (c) k that is fully localized on one site is characterized by an IPR of 1. By contrast, the IPR of a state uniformely delocalized over the whole graph is equal to N . Considering this, the mean IPR of the active states IP R AS is written as

IP R AS = 1 N C N C c=1 k∈AS IP R(|Ψ (c) k ) N (c) AS , (9)
where k∈AS represents a sum over the active states. The W dependance of IP R AS is shown in Fig. 5. This figure reveals that the evolution of the mean IPR of the active states is quite similar to that of N AS . Indeed, for small disorder W , the three curves start with a same initial value IP R AS ≈ 2.2. Then, when W increases the three curves reach different maximum values IP R AS max that arise in the neighbourhood of the optimal value of disorder W = W opt . For Γ = 1Φ the maximum is IP R AS max ≈ 4.2, whereas for Γ = Γ = 6Φ it is IP R AS max ≈ 3.6. The last situation Γ = 15Φ presents a maximum of IP R AS max ≈ 3.2. After reaching these maxima, the three curves decrease. In the strong disorder limit W W opt , they converge to the same limit value IP R AS = 1.

To keep on characterizing the spatial extension of the active states, the next property we have estimated is their mean weight Π core AS over the central trap. This quantity is defined as

Π core AS = 1 N C N C c=1 k∈AS | 0, 0|Ψ (c) k | 2 N (c) AS . ( 10
)
The figure 6 shows the W dependance of Π core AS . Once more, the three curves present a similar progression as the disorder W increases. For small disorder, the mean weight of the active states over the core site is really low. Typically, for W = 0.1Φ the amplitude of the three Then as W increases, the curves progressively increase to reach different maxima all arising near from the optimal value of disorder W = W opt . For Γ = 1Φ the maximum is Π core AS max ≈ 7.2%, whereas for Γ = Γ opt = 6Φ the value is Π core AS max ≈ 2.4%, and finally for Γ = 15Φ, Π core AS max ≈ 0.5%. After reaching these maxima, the curves decrease to return, in the strong disorder limit W W opt , to lower values. To cite an example, when W = 100Φ the curves present an amplitude ranging like Π core AS ∈ [10 -2 %, 10 -1 %].

Finally, to complete our study we have investigated the W dependence of the mean decay rate γ AS of the active states. This quantity is defined as

γ AS = 1 N C N C c=1 k∈AS γ (c) k N (c) AS . ( 11
)
The figure 7 illustrates the behavior of γ AS . Here again, whatever Γ, a quite similar behavior is adopted as the disorder W increases. First, for a weak disorder, the mean decay rate is really low. To give an example, when W = 0.1Φ the amplitude of the three curves typically ranges in an interval γ AS ∈ [10 -5 Φ, 10 -4 Φ]. Nevertheless, as W increases, the curves progressively increase to reach different maxima all arising near from the optimal value of disorder W = W opt . For Γ = 1Φ the maximum reached is γ AS max ≈ 3.6 × 10 -2 Φ, whereas for Γ = Γ opt = 6Φ the value is γ AS max ≈ 7 × 10 -2 Φ. Finally for Γ = 15Φ one obtains γ AS max ≈ 4.1 × 10 -2 Φ. After reaching these maxima, the three curves decrease to return, in the strong disorder limit W W opt , to lower values. For instance, when W = 100Φ the three present a quite similar amplitude γ AS ∼ 10 -3 Φ.

IV. DISCUSSION

In an extended star graph, our numerical results reveal the strong influence the disorder has on the excitonic absorption process. Therefore, various new properties have been noted. Let us briefly summerize the latter.

First, when W > 0, it has been shown that the absorption process is always complete. Indeed, whatever the strength of the disorder, an exciton starting on a peripheral site of a disordered graph is always fully trapped at the core site. This is a feature that clearly contrasts with the incomplete trapping process observed in an ordered graph 57 .

Then, it has been evidenced that the efficiency of the absorption process strongly depends on the model parameters. As a matter of fact, the evolution of the mean absorption time τ in the parameter space (Γ, W ) reveals the existence of an optimal regime of absorption reached when W opt = Γ opt ≈ 6Φ. In this case, the trapping process is strongly enhanced and the absorption time τ reduces to a minimum approximately equal to 50Φ -1 . Around this optimal region, the study of the W dependence of τ has also revealed the arising of a more general optimization process. Indeed, whatever Γ, a minimization of the absorption time τ always occurs in the neighbourhood of the optimal disorder W ∼ W opt .

To understand these behaviors, a study of the active states properties has been realized. Four observables have been investigated : the mean number N AS of active states, their mean inverse participation ratio IP R AS , their mean weight over the central trap Π core AS , and their mean decay rate γ AS . In this context, it has been shown that all these observables reach a maximum value when W ∼ W opt .

To interpret these features, let us characterize the fundamental ingredients that govern the exciton dynamic, i.e. the system eigenstates. Indeed, to be absorbed at the core site (0, 0), an exciton initially situated on ( 0 , s 0 ) has to tunnel across the network to reach the trap.

According to the fundamental principles of quantum mechanics, this tunneling process, that is encoded in the evolution operator U (t), results from the sum over all the various paths that the exciton can follow to reach the core site at time t. A "path" defines a transition through an eigenstate |Ψ k involving specifically its projections on the initial site 0 , s 0 |Ψ k and on the core site 0, 0|Ψ k . Knowing this, it becomes clear that the time evolution of the absorption process strongly depends on the paths that actively support the excitonic delocalization to the trap, i.e. the "active states". Depending on whether the graph is ordered or disordered, as well as depending on the strength of the absorption rate Γ, the properties of the active states can change and thus generate various regimes of absorption.

Based on this idea, the scheme of our discussion is simple. As W increases, we will progressively detail the way the disorder modifies the system eigenstates. To proceed, let us begin with the case of an ordered graph.

When W = 0, the N 1 -fold symmetry of the network produces two kinds of eigenstates. On the one hand, N -3 degenerated eigenstates are insensitive to the effect of the trap and exhibit purely real eigenenergies. On the other hand, three non-degenerate decaying states are sensitive to the trap and exhibit complex eigenenergies. Due to their extension over the graph, the three decaying states are the only ones that produce an excitonic tunneling to the core and that participate in the trapping process. However, because of their weak weights on the initial site ( 0 , s 0 ), these eigenstates do not represent real "active states" : they only govern the dynamics of a weak part of the total excitonic population. Indeed, as demontrated in our previous work 57 , when an exciton starts from a peripheral site, its dynamic is mainly governed by the highly degenerate eigenstates which only interfer on the N 1 periphery stars of the network. The resulting interferences produce a localization of the exciton on its initial peripheral star which therefore hinders the propagation to the trap. As a consequence, the total absorbed population at the core is limited and reaches a maximum of

P max A = 1/N 1 N 2 . Note that with N 1 = N 2 = 4,
one obtains the value P max A = 6.25×10 -2 as observed in Fig. 2.

If now a weak disorder W = 0 + pertubates the graph, interesting features arise. In this context, the presence of local defects breaks the symmetry of the graph. As a result, all the degeneracies are raised and each eigenstate exhibits a non-zero decay rate. Indeed, due to the disorder, a spatial restructuring of the eigenstates occurs so that each state presents a non-vanishing weight on the core site. As a consequence, each state becomes sensitive to the trap and acquires a finite lifetime. The excitonic population initially located at the periphery of the network can thus entirely transfers to the trap, resulting in a complete absorption process, i.e. P max A = 1. Nevertheless, in the weak disorder limit, the absorption process remains inefficient. The reason is quite simple to understand. When W = 0 + , the eigenstates are weakly pertubed by the disorder so that the graph still keeps the memory of its properties in the absence of disorder. On the one hand, three short-lived eigenstates still remain delocalized over the graph and correspond to the relics of the three decaying states. On the other hand, N -3 long-lived eigenstates are still mainly localized on the peripheral stars of the graph and correspond to the relics of the degenerated states. Due to their very weak (but non-vanishing) weight on the central trap, these states exhibit very small decay rates compared with the three other short lived states.

In this context, among the N -3 long-lived states are the "active states", i.e. the eigenstates transmitting the main part of the excitonic population to the trap. Due to their long-lived nature, the resulting absorption process varies slowly and its time evolution exhibits two phases (as observed in Fig. 2). First, over a short timescale, the trapped population increases from zero to 1/N 1 N 2 owing to the action of the three "relic decaying states". Then, due to the influence of the long-lived active states, the absorption process enters in a second phase. Slowly, the active states transmit the remaining excitonic population to the trap. As a result, the absorbed population P A (t) progressively increases and it converges to the unity. Now, if W increases, the restructuring process induced by the disorder intensifies. As a result, the absorption efficiency is strongly affected and it exhibits a transition depending on the value of W . First, when W increases the absorption process is enhanced and the time needed to trap the exciton reduces until it reaches a minimum. This minimum occurs in the neighbourhood of the optimal value W = W opt . Second, when W exceeds this optimal value, the trapping process progressively loses in efficiency and the absorption time increases.

Understanding the precise origins of this transition is a very difficult task since the latter results from the interplay between different complex mechanisms.

However, deeper numerical investigations of the active states have allowed us to gather some key information to intepret this phenomenon.

To introduce these features, one will realize a simple thing : tracking the active states for a particular case of disordered configuration as a function of W . Of course, with this approach we do not pretend to be fondamentally exhaustive. However, based on a simple example, we want here to highlight some typical mechanisms that could explain the existence of the optimal disorder. In this context, Fig. 8 shows the results obtained for a given disordered configuration. On the left part of the figure, the evolution of the eigenenergies is shown (black curves). The red curves represent the eigenenergies of the active states (see color online). On the right part of the figure the spreading of the active states on the extended star graph is represented for a) a weak disorder, b) an optimal disorder and c) a strong disorder.

In the weak disorder regime, a single active state is present. The corresponding eigenenergy is located at the center of the excitonic spectrum in a place that originately contains a lot of degenerate states when W = 0. As a relic degenerate state, this active state remains mainly localized on the peripheral star that initialy carries the exciton (see Fig. 8.a).

When W increases, the restructuring process affecting the eigenstates strengthens which produces a shift of the eigenenergies. When W gets close to the optimal value W opt , the amplitude of the disorder becomes typically equal to the width of the excitonic spectral band known in absence of disorder 2Φ √ N 1 + N 2 ∼ 5.7Φ. In this case, the disorder effect is then strongly nonperturbative. The restructuring process affecting the eigenstates becomes so strong that all the eigenenergies drastically deviate in the spectrum. As a result, avoided crossings appear which indicates the emergence of a new phenomenon : the hybridization of the system eigenstates. Indeed, owing to the presence of the disorder some resonances may arise between states localized in different regions of the graph. As a consequence, a complex mixing between different relic eigenstates can occur, which generates new paths for the exciton. In this context, as shown in Fig. 8.b, the active states properties evolve : they become more numerous (three states detected when W = W opt ), and they better spread over the graph. In addition, their weight over the core site is globally enhanced as well as their widths, which reveals their better sensitivity to the action of the trap (not represented here). In this new dynamical context, the excitonic transport changes. The cooperative action of the new active states facilitates the excitonic delocalization from the initial peripheral site to the core of the graph. As a result, the absorption time strongly reduces and the trapping process gains in efficiency. At this point, let us note that all these properties clearly corroborate the maximization of all the observables N AS , IP R AS , Π core AS and γ AS obtained when W ∼ W opt . Now, if W exceeds W opt , the hybridization process progressively vanishes to give way to another phenomenon : the so-called Anderson localization 58 . Indeed, in the strong disorder limit, the local defects become so important that each eigenstate tends to localize on one particular site. As a consequence, as shown in Fig. 8.c, the number of active states decreases and only one active state subsists. Due to its strong localization on the site initialy excited, this state becomes weakly sensitive to the trap effect : its weight on the core site and its decay rate vanish when W increases (not represented here). In this context, the exciton remains stuck over a very long time on its initial site so that the propagation from the periphery to the core is hindered. As a result, the absorption time drastically increases and the trapping process loses its efficiency. These properties corroborate the limit behavior observed when W W opt for N AS , IP R AS , Π core AS and γ AS . Therefore, within this simple example we have depicted a typical scenario explaining how the eigenstates restructuring induced by the disorder can strongly enhance the absorption process when W ∼ W opt . Nevertheless, one last feature still has to be tackled : when W = W opt , the absorption process can be even more optimized if the absorption rate Γ tends to Γ opt . This process gives rise to the occurrence of the optimal region of absorption as presented in Fig. 3.

Actually, the emergence of the optimal region of absorption shall be seen as the result of the interplay between two phenomena. On the one hand, the positive retructuring of the eigenstates induced by the disorder. And on the other hand, the arising of the so-called superradiance transition [36][37][38][39] .

To better understand this feature, let us consider the optimal disorder W = W opt as a starting point. In this situation, we know that the restructuring of the eigenstates produces new active states exhibiting a maximal weight over the core of the graph (see Fig. 6). As a consequence, these states are in their best predisposition to delocalize the exciton to the trap. Therefore, if Γ starts to increase, the sensitivity of these states to the trap intensifies, which enhances their width (see in Fig. 7 the increasing of γ AS max when Γ varies from 1Φ to Γ opt ). However, due to the occurrence of a countering phenomenon, this enhancement exhibits a limit. Indeed, as Γ increases, a non-active state progressively localizes on the trap. As a consequence, all the active states tend to lose their extension over the core of the graph (see in Fig. 6 and 5 the decreasing of Π core AS max and IP R core AS max when Γ increases). When Γ < Γ opt , these effects are not significant so that the active states keep on improving their interaction with the trap, i.e. their widths increase with Γ. However, when Γ > Γ opt , the localization process of the non-active state is so strong that the trap starts to behave as a singular site almost disconnected from the rest of the graph. In this context, the beneficial interaction of the active states with the trap gets hindered and a segregation occurs in the decay rates of the eigenstates. On the one hand, the strongly localized non-active state becomes superadiant : its decay rate keeps on increasing with Γ. On the other hand, the active states become subradiant : their decay rate decreases with Γ (see in Fig. 7 γ AS max which diminishes when Γ > Γ opt ). This segregation process that is arising in the neighbourhood of Γ opt is the signature of a superradiance transition which gives birth the optimal region of absorption observed in Fig. 3.

V. CONCLUSION

In this paper, a tight binding model was introduced for studying the dynamics of an exciton moving on a disordered extended star graph whose central site is occupied by a trap. In a very surprising way, our work has evidenced how positively the disorder can influence the absorption process occuring at the core of the graph.

First, in a marked contrast with what happens on an ordered graph, it has been shown that the absorption process is always complete on a disordered graph. The origin of this phenomenon is intimately linked to the symmetry breaking that affects the system. In this context, all the degenerated states that originately hindered the excitonic delocalization get perturbated so that the latter can participate to the excitonic transport. As a consequence, an exciton starting on a peripheral site can entirely transfer to the trap resulting in the arising of the full absorption process.

Second, the existence of optimal values of disorder strength and absorption rate have been evidenced. Indeed, when simultaneously both the parameters get typically equal to the width of the unperturbated excitonic spectrum, the time needed to fully trap the exciton strongly reduces. Our study has revealed that the origin of this enhancement results from the beneficial interaction of two complex phenomena. First, when the strength of the disorder becomes optimal, a mixing of several relic eigenstates occurs which produces new hybrid states. Thanks to their better spatial extension, these new hybrid states facilitate the delocalization of the exciton to the trap, which speeds up the trapping. Then, a complementary improvement can be done when the absorption rate also reaches its optimal value. In this context, the arising of a superradiance transition allows to maximally enhance the sensitivity of the new hybrid eigenstates to the trap action. As a result, the excitonic population reaching the core of the graph gets trapped even more fastly resulting in the arising of a superoptimal absorption process.

Naturally, the study introduced in this paper must be seen as a preliminary work that should open on some complementary investigations. For instance, an interesting question would be to determine how the size parameters of the graph can influence the optimization of the absorption process. In another context, it could also be wise to reconsider this study with a more realistic physical point of view. To do so, considering additional effects of the excitonic environment, such as exciton-phonon interactions and optical recombination, would represent the next logical step of this work.

FIG. 1 :

 1 FIG.1:The extended star graph with an absorbing core.

7 FIG. 2 :

 72 FIG. 2: Time evolution of the absorbed population P

A

  (t) obtained for particular disordered configurations generated with different values of W and Γ.

FIG. 3 :

 3 FIG. 3: (See color online) Evolution of the mean absorption time τ , (a) depending on the strength of the disorder W and the absorption rate Γ, and (b) for three fixed values of the absorption parameter Γ considering W as a free parameter. On (b), Γ = 1Φ is represented with black circles, Γ = Γopt = 6Φ with green triangles and Γ = 15Φ with blue squares. Each value τ is the result of a statistic realized over NC = 2000 configurations.

3 FIG. 4 :

 34 FIG. 4: Mean number NAS of the active states versus the disorder strength W for Γ = 1Φ (black circles), Γ = Γopt = 6Φ (green squares) and Γ = 15Φ (blue triangles). Each point of the graph is the result of a statistic realized over NC = 2000 configurations.

3 FIG. 5 :

 35 FIG. 5: Mean IPR IP RAS of the active states versus the strength W for Γ = 1Φ (black circles), Γ = Γopt = 6Φ (green squares) and Γ = 15Φ (blue triangles). Each point of the graph is the result of a statistic realized over NC = 2000 configurations.

3 FIG. 6 :

 36 FIG. 6: Mean weight on the core site Π core AS of the active states versus the disorder strength W for Γ = 1Φ (black circles), Γ = Γopt = 6Φ (green squares) and Γ = 15Φ (blue triangles). Each point of the graph is the result of a statistic realized over NC = 2000 configurations.

3 FIG. 7 :

 37 FIG. 7: Mean decaying rate γAS of the active states versus the disorder strength W for Γ = 1Φ (black circles), Γ = Γopt = 6Φ (green squares) and Γ = 15Φ (blue triangles). Each point of the graph is the result of a statistic realized over NC = 2000 configurations.

FIG. 8 :

 8 FIG. 8: Tracking of the active states for one particular configuration of disorder for Γ = 1Φ. Here, the local defects are preserved and scaled with respect to the amplitude W . On the left part of the figure, the graph represents in black the evolution of the real eigenenergies ω k as a function of W . The red curves indicate the energies of the active states. On the right part of the figure, an illustration of the spatial extension of the active states is given for three different cases : a) weak disorder W Wopt, b) optimal disorder W = Wopt and c) strong disorder W Wopt. The size of the balls stands for the weights of the active state on each node whereas the arrow indicates the initial excited site.


Amélie Rouchon
email: amelie.rouchon@cea.fr

Richard Sanchez
email: sanchezr@snu.ac.kr

Igor Zmijarevic
email: igor.zmijarevic@cea.fr

The new 3-D multigroup diffusion neutron noise solver of APOLLO3 ® and a theoretical discussion of fission-modes noise

Traditional neutron noise analysis addresses the description of time-dependent flux fluctuations induced by small global or local perturbations of the macroscopic cross-sections, which may occur in nuclear reactors due to stochastic density fluctuations of the coolant, to vibrations of fuel elements, control rods or any other structures in the core. Neutron noise equations are obtained by assuming small perturbations of macroscopic cross-sections around a steady-state neutron field and by subsequently taking the Fourier transform in the frequency domain. In this work, we present the new 3-D multigroup diffusion neutron noise solver implemented in APOLLO3 ® , the new multi-purpose deterministic nuclear code under development in CEA. We illustrate the capacities of this new 3-D diffusion neutron noise solver by performing two neutron noise simulations in a large pressurized water reactor with heavy baffle in three dimensions: a cross-sections oscillation and a traveling perturbation. Moreover, we give a separate analysis of the neutron noise anomaly recently observed in KWU PWRs while proving the existence of steady-state like global noise modes in the low frequency spectrum.

I. INTRODUCTION

Neutron noise appears as fluctuations of the neutron field induced by stochastic or deterministic changes in the macroscopic cross-sections [START_REF] Pázsit | Noise Techniques in Nuclear Systems[END_REF]. The latter may result from vibrations of fuel elements, control rods or any other mechanical structures in the core, as well as from global or local fluctuations in the flow, density or void fraction of the coolant [START_REF] Fry | Use of Neutron Noise of Diagnosis of In-vessel Anomalies in Light-Water Reactors[END_REF]. In power reactors, ex-core and in-core detectors can be used to monitor neutron noise with the aim of detecting possible anomalies and taking the necessary measures for continuous safe power production.

The general noise equations are obtained by assuming small perturbations around a steady-state neutron flux and by subsequently taking the Fourier transform in the frequency domain. The analysis is performed based on the neutron kinetic equations, including the coupling with neutron precursors. The outcome of the Fourier transform analysis is a fixed-source equation with complex operators for the perturbed neutron field, which can then be solved so as to predict noise measurements at detector locations. For each frequency, the neutron noise is a complex function having an amplitude and a phase.

Noise analysis relies upon the possibility of numerically simulating the behaviour of neutron noise and computing the changes in the neutron field produced by different representative sources of noise in reactor cores. These analysis are generally performed via analytical techniques [START_REF] Pázsit | Theoretical Investigation of the Neutron Noise Diagnostics of Two-Dimensional Control Rod Vibrations in a PWR[END_REF][START_REF] Jonsson | Analytical investigation of the properties of the neutron noise induced by vibrationg absorber and fuel rods[END_REF] or by resorting to diffusion theory with two energy groups [START_REF] Demazière | CORE SIM: A multi-purpose neutronic tool for research and education[END_REF][START_REF] Malmir | Development of a 2-D 2-group neutron noise simulator for hexagonal geometries[END_REF].

In this work, we present a new 3-D multigroup diffusion neutron noise solver implemented in APOLLO3 ® , the new multi-purpose deterministic nuclear code under development in CEA [START_REF] Schneider | APOLLO3: CEA/DEN Deterministic multi-purpose code for reactor physics analysis[END_REF], and an analysis of the neutron noise anomaly recently observed in KWU PWRs. This paper is organized as follows. In Sec. II, the general neutron noise theory will be briefly introduced. In Sec. III, the new neutron noise solver developed in APOLLO3 ® will be described and we will analyze a cross-sections oscillation and a traveling perturbation in a large pressurized water reactor with heavy baffle in three dimensions.

Finally, in Sec. IV, we give a separate analysis of the neutron noise anomaly recently observed in KWU PWRs. We show that in the low frequency spectrum all frequency components of the noise equations can be neglected and that to leading order the equation reduces to a steady-state equation for the noise flux with given source. This suggests that global modes close to the steady state ones can be excited as these low frequencies and prompts us to provide a method to identify a mode which has the experimentally observed symmetries of the noise anomaly. The appendix summarizes the properties of some simple analytical solutions of the one-group diffusion equation that are used as examples in our analysis.

Conclusions will be presented in Sec. V.

II. NEUTRON NOISE THEORY

Here we summarize the general theory of neutron noise in transport theory and we detail the final equations in diffusion theory. In diffusion theory, in addition to the Fick's law application, fluctuations of the diffusion coefficients are disregarded and the term 1 v ∂ t J is neglected (with J the scalar current). Note that the zero power reactor noise (fluctuations inherent to the branching process) is neglected in power reactor noise theory [START_REF] Pázsit | Noise Techniques in Nuclear Systems[END_REF].

We assume small perturbations of the macroscopic crosssections around the following critical steady state:

B 0 ψ 0 (r, Ω, E) = 0, (1)
where ψ 0 is the steady-state angular flux and B 0 = Ω • ∇ + Σ 0 -H 0 -P 0 the steady-state Boltzmann operator with Σ 0 the steady-state total cross-section, H 0 the steady-state scattering operator and P 0 the steady-state production operator. For the steady state, the effective multiplication factor is assumed to be k = 1. For a system in diffusion theory with one precursor group and one fissile isotope, the critical steady-state Boltzmann equation for each energy group g is (all notations are standard):

-∇ • D g (r)∇ + Σ g 0 (r) Φ 0,g (r) = g Σ g →g 0,s 0 (r)Φ 0,g (r) + χ g k g ν g Σ g 0, f (r)Φ 0,g (r), (2)
with Φ 0 the steady-state scalar flux. We impose a temporal perturbation of the macroscopic cross-sections, which yields the kinetic equation:

1 v ∂ t + B ψ(r, Ω, E, t) = 0, (3)
where ψ is the angular flux, v the neutron velocity and B = Ω • ∇ + Σ -H -P the kinetic Boltzmann operator with Σ the total cross-section, H the scattering operator and P the production operator containing prompt and delayed neutron contributions. We impose a periodic perturbation of the kinetic operator with a period T 0 , in the form:

B = B 0 + δB. ( 4
)
This perturbation is supposed to start at time t = -∞, so that we can reasonably assume that the asymptotic perturbation regime is attained. A similar decomposition is also used for the angular flux:

ψ(r, Ω, E, t) = ψ 0 (r, Ω, E) + δψ(r, Ω, E, t), (5)
where the perturbation term δψ is called "neutron noise". Finally, plugging expressions (4) and ( 5) into Eq. ( 3) leads to a kinetic source equation for the neutron noise:

1 v ∂ t + B δψ(r, Ω, E, t) = -δBψ 0 (r, Ω, E). (6)
The second order term δBδψ will be neglected, so that we obtain the traditional linearized kinetic equation:

1 v ∂ t + B 0 δψ(r, Ω, E, t) = -δBψ 0 (r, Ω, E). (7)
We want to determine the unique periodic solution of this equation. We apply the Fourier transform and we obtain the noise equation in the usual form:

B 0,ω δψ(r, Ω, E, ω) = -δB ω ψ 0 (r, Ω, E), (8)
where B 0,ω = i ω v +Ω•∇+Σ 0 -H 0 -P 0,ω is a modified (complex) Boltzmann operator, i the imaginary unit and ω = 2π f the angular frequency. The right hand side of Eq. ( 8) represents a (known) "noise source".

For a system in diffusion theory with one precursor group and one fissile isotope, the noise equation is:

-∇ • D g (r)∇ + Σ g 0 (r) δΦ g (r, ω) = -i ω v g δΦ g (r, ω) + g Σ g →g 0,s 0 (r)δΦ g (r, ω) + χ g p k g ν g p Σ g 0, f (r)δΦ g (r, ω) + χ g d k g ν g d,ω Σ g 0, f (r)δΦ g (r, ω) + S g (r, ω), (9)
with

ν g d,ω = λ 2 λ 2 + ω 2 -i λω λ 2 + ω 2 ν g
d and S g the noise source defined by: S g (r, ω) = -δΣ g (r, ω)Φ 0,g (r) + g δΣ g →g s 0 (r, ω)Φ 0,g (r)

+ 1 k g χ g p ν g p + χ g d ν g d,ω δΣ g f (r, ω)Φ 0,g (r), (10)
where δΣ g x (r, ω) is the Fourier transform of the perturbed term of the macroscopic cross-section Σ g x (r, t) = Σ g 0,x (r) + δΣ g x (r, t). Thus, because of the delayed neutrons, the production operator P 0,ω depends on the frequency. The real and imaginary components of the neutron noise are coupled by two terms: i ω v and the modified production operator P 0,ω .

III. APOLLO3 ® 3-D MULTIGROUP DIFFUSION NEUTRON NOISE SOLVER: DESCRIPTION AND EXAMPLES

APOLLO3 ® is the new multi-purpose deterministic nuclear code under development in the frame of the neutronics simulation project of the Nuclear Energy Division (DEN) of the CEA and with financial support from AREVA and EDF [START_REF] Schneider | APOLLO3: CEA/DEN Deterministic multi-purpose code for reactor physics analysis[END_REF]. In this section we will present the new 3-D multigroup diffusion neutron noise solver developed in this new generation deterministic code. Then, we will detail two neutron noise simulations in a large pressurized water reactor with heavy baffle in three dimensions: a cross-sections oscillation and a traveling perturbation.

Description of the 3-D multigroup diffusion neutron noise solver

One of the diffusion solvers implemented in APOLLO3 ® is a 3-D multigroup nodal diffusion solver based on the classical Nodal Expansion Method [START_REF]Piecewise Assembly Homogenization via Flux Discontinuity Factors[END_REF][START_REF] Gupta | Nodal methods for three-dimensional simulators[END_REF]. In order to solve Eq. ( 9) (which is an equation of a fixed-source problem with a coupling between the real and imaginary parts of δΦ), we choose to work with this diffusion solver: we apply the same iteration loops to the fission source (but the production operator is now complex) and to the scattering source as customary, and we add an iteration loop between the real and imaginary parts of the neutron noise equation to the solution of the one-group problem. Thus, we can use the same one-group diffusion solver as customary.

For the moment, the new 3-D multigroup diffusion neutron noise solver developed in APOLLO3 ® can be applied to homogeneous cartesian geometries.

This new solver was verified by comparaison with some analytical results in a homogeneous infinite medium and in an infinitely-long homogeneous cylinder exactly as in [START_REF] Yamamoto | Monte Carlo method with complexvalued weights for frequency domain analyses of neutron noise[END_REF].

Description of the 3-D core

In order to illustrate the capacities of the new 3-D multigroup diffusion neutron noise solver, we have performed neutron noise simulations in a large pressurized water reactor with heavy baffle in three dimensions. This 3-D core is composed of 17×17 homogeneous assemblies (x-and y-dimensions of each assembly are 21.5×21.5 cm) and its axial height is 420 cm. The 3-D core is decomposed into 21 axial planes of 20 cm: the first and the last planes are only composed by homogeneous heavy baffle assemblies (noted RE in Fig. 1) and all other planes are set by the same fuel loading map described in Fig. 1.

All simulations have been performed in diffusion theory with 2 energy groups (fast and thermal groups), 8 precursors groups, a quartic approximation of the volume flux and a corrected scattering anisotropy P * 0 (i.e., the diffusion coefficients are defined by 1/[3(Σ g 0,a + (Σ g 0,s 0 ) * )] where all macroscopic selfscattering cross-sections are corrected by (Σ g→g 0,s 0 ) * = Σ g→g 0,s 0g Σ g→g 0,s 1 with g the energy group). Figure 2 presents the fast and thermal steady-state fluxes obtained in diffusion theory with APOLLO3 ® . The eigenvalue of this system is close to 1 (k = 0.99674). Figure 3 details the steady-state flux the UOXGD assembly in position (5,7) (see Fig. 1, it is the future perturbed assembly) in function of the axial position z. As expected, the axial steady-state distribution follows the sine shape.

Analysis of a cross-sections oscillation

We present in this section the neutron noise induced by a cross-sections oscillation of the UOXGD assembly in position (5,7) (see Fig. 1). All macroscopic cross-sections of this homogeneous assembly are periodically perturbed on the whole height of the assembly (380 cm, the first and the last planes are excluded) around their steady-state value Σ g 0,x with an amplitude of 1% and a frequency of 3 Hz. Thus, the perturbed term δΣ g

x is defined by: with ω 0 the angular frequency of the perturbation (ω 0 = 6π rad/s). Note that the noise source is monochromatic so the neutron noise is also monochromatic. Figures 4 and5 detail the fast and thermal neutron noises moduli and phases at 3 Hz of the 5th axial plane. Observe that the neutron noise is localized around the perturbed assembly and the phases are almost constant in the whole core (the axial distribution of the phases is almost constant). Figure 6 presents the fast and thermal neutron noises moduli in the perturbed assembly in function of the axial position z. As expected, these axial distributions follow the sine shape as the steady-state distribution.

δΣ g x (r, t) = 0.01 × Σ g 0,x (r) cos(ω 0 t), (11)

Analysis of a traveling perturbation

We briefly analyze in this section a traveling perturbation in the 3-D core described previously. We perturb the same assembly as in the previous section (see Fig. 1) by a periodically traveling perturbation with an amplitude of 1% defined by:

δΣ g x (z, t) = δΣ g x (z = z i , t - z v e ) = 0.01 × Σ g 0,x (z = z i ) cos ω 0 (t - z v e ) , (12)
where z is the axial position, z i the beginning axial position of the perturbation (z i = 20 cm, i.e., the beginning of the fuel assembly), ω the angular frequency of the perturbation and v e = 190 cm/s the perturbation axial velocity (supposed constant whatever is the axial position z).

Here we only present how we can find the perturbation axial velocity v e thanks to the numerical results of the phase shift between two axial positions. Figure 7 details the phase shift versus frequency between axial positions z a = 60 cm and z b = 140 cm included in the perturbed assembly. According to the perturbation definition, the neutron noise in z b is phase shifted by -ωτ compared to the neutron noise in z a , with τ = (z bz a )/v e the propagating time of the perturbation. So, according to Fig. 7, τ ≈ 0.421 s and thus we find v e ≈ 189.9 cm/s.

IV. A COMMENT REGARDING ABNORMAL NOISE IN KWU PWRS

Recently, an anomaly of the Low-Frequency neutron noise Signals (LFNS) was observed in KWU PWRs [START_REF] Seidl | Review of the historic neutron noise behavior in German KWU built PWRs[END_REF]. The anomaly consisted of a period of increase and then decrease of the magnitude of the LFNS in the Von-Konvoi plants. A comparison of measurements done in Von-Konvoi and Konvoi plants showed that the auto-power spectral densities (APSD) at ex-core detectors were similar in shape and that in the LFSN domain ω ∈ [0.5, 10] Hz this shape showed a white noise behavior in 1/ω 2 .

More to the object of this discussion, measurements in inner and outer detectors revealed that along the axial direction the neutron noise at the top and bottom of the core was in phase, while it had opposite phase across opposite quadrants. This global noise response suggests that the perturbations excited one of the steady-state core modes which the same symmetries and, in particular, the simplest mode with the higher eigenvalue. Modes with opposite radial behavior have been calculated for a 2D PWR reactor [START_REF] Carney | Theory and applications of the fission matrix method for continuous-energy Monte Carlo[END_REF] (see, for instance, modes 1, 2, 6, 7, 13 and 14 in Fig. 1 of [START_REF] Carney | Theory and applications of the fission matrix method for continuous-energy Monte Carlo[END_REF]). A more compelling example is provided by analytical expressions for the one-group, homogeneous diffusion theory. For instance, the modes for a bare cylindrical core of radius R and height H are

Φ n r ,m,n z (ρ, ϕ, z) = J n r ( j n r ,m ρ R ) cos(n r ϕ + ϕ 0 ) sin(n z π z H ), ( 13
)
where j n r ,m is the m-th positive zero of the Bessel function of integer order J n r and ϕ 0 ∈ [0, 2π). The mode with the appropriate symmetries and highest eigenvalue is given by n r = m = 1 and n z = 1 and for a critical core its eigenvalue λ 1,1,1 = 1/(1 + Γ), where Γ = (D/νΣ f )( j2 1,1j 2 0,1 )/R 1, is very close to 1. Because the geometry is invariant under axial rotations this mode is degenerated with a continuum of states, as given by the angle ϕ 0 , so one must assume that only one of the modes is excited. For a bare core with the shape of a rectangular prism, there is also a double degeneracy of this mode and again one has to postulate that only one of the two modes are excited. For the reader convenience, a detailed discussion of these elementary diffusion results is given in the appendix.

In this section we give theoretical support to our assumption. We start with the transport (or diffusion) noise equation that we write as

( iω v + B 0 -P 0,ω )δψ ω = -δB ω ψ 0 , (14)
To obtain this equation we have Fourier transformed the classical kinetic equations with delayed neutrons. Then we have introduced a reference steady-state critical transport operator

B 0 = Ω • ∇ + Σ 0 -H 0 -P 0 (the diffusion case is obtained by the replacement Ω • ∇ → -∇ • D∇)
, where Σ 0 , H 0 and P 0 are the total cross section and the scattering and production operators, respectively, and

P 0,ω f = iω d χ d λ d + iω f β d, f (υΣ f f,0 , f ) (15)
is the a component of the contribution of delayed neutrons to the neutron production term. Also, in the last equation β d, f (10 -4 -10 -3 ) is the fission yield of fissile isotope f for precursors in group d and λ d and χ d denote the decay constant and the fission spectrum for delayed neutrons. To derive Eq. ( 14) we have followed the traditional technique and used the critical flux ψ 0 , associated to the reference steady-state operator B 0 , to introduce a perturbation formula so that B ω = B 0 + δB ω and ψ ω = ψ 0 + δψ ω , where δB ω and δψ ω are assumed small perturbations and, finally, we have neglected products of two small quantities. Also, in Eq. ( 14) δB ω = δΣ ω -δH ω -δP ω contains the Fourier transforms of the perturbations of the cross sections about their reference values in B 0 .

We note that for the LFNS the term in ω/v can be safely neglected in Eq. ( 14). Moreover, the production operator P 0,ω can also be neglected as compared with the steady-state production term P 0 . The resulting equation for the LFNS noise,

B 0 δψ ω -δB ω ψ 0 , (16)
has a degenerated operator on the left-hand-side and, according to Fredholm's alternative, a solution exist iff δB ω ψ 0 is orthogonal to ψ 0 : (ψ † 0 , δB ω ψ 0 ) = 0, where ψ † 0 is the adjoint flux and (•, •) denotes the scalar flux. Granted that this condition is satisfied, the solution of Eq. ( 16) can be written as

δψ ω = z 0 ψ 0 + ψ ⊥ , z 0 ∈ Z, (17)
where the latter flux component is orthogonal to ψ 0 : (ψ † 0 , ψ ⊥ ) = 0. Note that both z 0 and ψ ⊥ have to depend on the frequency ω.

Our approach to neutron noise equations for power reactors accounts for the fact that the averaged power is kept constant via the action of built in mechanisms or control rods. [START_REF] Sanchez R | Some comments in neutron noise theory[END_REF] This leads to a construction of operator B 0 so the iff condition (ψ † 0 , δB ω ψ 0 ) = 0 is automatically fulfilled with z 0 = 0 in Eq. ( 17). Therefore, we set z 0 = 0 and use an expansion of the orthogonal term over the eigenfunctions {λ n , ψ n } of the steady-state problem 2

δψ ω = n>0 z n ψ n , z n ∈ Z. (18)
We notice again that the z n ought to depend on ω. At this point we use the experimental evidence to assume that this dependence must be in 1/ω 2 . To simplify our notation hereafter we shall omit this common dependence on the frequency and concentrate on the spatial dependence. Hence, we shall consider Eq. ( 18) with the z n as simple complex constants. Finally, use of this expansion in Eq. ( 16) leads to a final expression for the noise source δS ω = -δB ω ψ 0 :

S ω = n>0 z n λ n ψ n , z n ∈ Z. (19)
Since the detailed spatial behavior of the neutron noise ψ ω is not known, our approach consists of assuming that all the modes in Eq. ( 18) are very small except the mode ψ * which is the simplest mode that reproduces all the experimentally observed symmetries with exclusion of other behavior. Among the modes which exhibit axial symmetry and radial antisymmetry this is the one with the highest eigenvalue. Hence, we assume that the source is

S ω = 1 ω 2 ψ * (r), (20)
where we have put back the frequency dependence. It remains, however, to characterize the original source of the noise, i.e., the δB. The explanation given in [START_REF] Seidl | Review of the historic neutron noise behavior in German KWU built PWRs[END_REF] points to inlet temperature fluctuations. Because the anomaly in the noise flux was coincidental with a change of fuel elements which grids that increased turbulence, we think that one important source of noise is due to the turbulence generated at the grids, which could generate a volume-distributed noise source as the one in Eq. ( 20) and lead to the excitation of the flux steady-state mode ψ * (r). To confirm or disapprove this thesis it remains, thus, to effectuate the pertinent numerical experimentation and this is not in the scope of this work. The observed increase of noise with burnup could provide a led by investigating first which are the components of δB ω which increases the most with burnup.

V. CONCLUSIONS

In this paper we have presented a new 3-D multigroup diffusion neutron noise solver implemented in the new generation deterministic code APOLLO3 ® developed at CEA. We illustrated the capacities of the new 3-D diffusion neutron noise solver by performing two neutron noise simulations in a large pressurized water reactor with heavy baffle in three dimensions. Future work will concern the development of a new 3-D transport neutron noise solver in APOLLO3 ® and the comparison with neutron noise simulations performed with the reference Monte Carlo code TRIPOLI-4 ® developed at CEA [START_REF] Brun | TRIPOLI-4 ® , CEA, EDF and AREVA reference Monte Carlo code[END_REF][START_REF] Rouchon | A new Monte Carlo method for neutron noise calculations in the frequency domain[END_REF].

In a separated analysis we have portrayed the lowfrequency neutron noise signal in some German PWRs as resulting from the excitation of a high-order steady-state like transport mode with the symmetries observed in the noise measurements. This leads to an inverse problem for the perturbations of the cross sections or to a trial and error search for the perturbations with guidance from the symmetries of the observed neutron noise. A result of this analysis is that large core modes can contribute to the LFNS noise with frequencies as high as 10 Hz.

VI. ACKNOWLEDGMENTS

APOLLO3 ® is a registered trademark of CEA. The authors gratefully acknowledge AREVA and EDF for their long term partnership and their support.

The final writing of Sec. IV was done at Seoul National University, where one of the authors (RS) has been invited under the Brain Pool Program. This work was supported by the Brain Pool Program through the Korean Federation of Science and Technology Societies (KOFST) funded by the Ministry of Science, ICT and Future Planning.

APPENDIX: ONE-GROUP DIFFUSION RESULTS FOR A BARE HOMOGENEOUS CORE

For the one-group homogeneous core with vacuum boundary conditions elementary diffusion theory gives the dispersion equation:

λ = k ∞ 1 + L 2 B 2 , ( 21
)
where k ∞ = νΣ f /Σ a is the infinite medium multiplication factor and L 2 = D/Σ a is the square of the diffusion length. The values of the so-called material bucklings B 2 are the eigenvalues of the Laplacian equation

-∆Φ = B 2 Φ, (22)
where ∆ is the Laplacian operator and Φ is the flux associated to eigenvalue λ.

It is well known that the latter problem has a real point spectrum and that the maximum eigenvalue is simple and its eigenfunction is nowhere negative. The λ eigenvalues increases with decreasing buckling. Hence, the dominant eigenvalue λ 0 corresponds to the minimum value B 2 0 of B 2 . We shall assume that the reactor is critical so that λ 0 = 1. Therefore, for any other eigenvalue λ one has

λ = λ 0 - Γ 1 + Γ = 1 1 + Γ , where Γ = L 2 k ∞ δB 2 with δB 2 = B 2 -B 2 0 .

Cylindrical core

We consider a homogeneous cylinder of height H and radius R. In cylindrical coordinates (ρ, ϕ, z) the Laplacian reads

∆ = ∂ 2 ρ + 1 ρ ∂ ρ + 1 ρ 2 ∂ 2 ϕ + ∂ 2 z . (23)
We want to solve eigenvalue problem (22) with vacuum boundary conditions:

Φ(R, ϕ, z) = 0, 0 < z < H, 0 ≤ ϕ < 2π, Φ(ρ, ϕ, 0) = Φ(ρ, ϕ, H) = 0, 0 ≤ ρ < R, 0 ≤ ϕ < 2π.
(24) Because the boundary conditions are independent of the azimuthal angle ϕ and apply separately to the radial ρ = R and the two horizontal z = 0 and z = H surfaces, we can write the general solution of problem (22) as the product of three functions of a single variable, each one or which depending on a different spatial coordinate:

Φ(ρ, ϕ, z) = R(ρ)F(ϕ)Z(z).
Substitution of this anszat in Eq. ( 22) gives the solution:

Φ ν,m,n r ,n z (ρ, ϕ, z) = J υ ( j ν,m ρ R ) cos(n r ϕ + ϕ 0 ) sin(n z π z H ).
where J ν is the Bessel function of order ν, j ν,m is the m-th root (excluding zero) of J ν , n r ≥ 0 and n z > 0 are integers and ϕ 0 ∈ [0, 2π). Under replacement in (22) and use of the boundary conditions (24) we find that the constraint

ν 2 ρ 2 - j 2 ν,m R 2 - n 2 r ρ 2 -(n z π H ) 2 = -B 2 ,
requires to select ν = n r ≥ 0 integer and results in the final relation

B 2 = j 2 ν,m R 2 + (n z π H ) 2 (25)
with the final eigenmodes of the form

Φ n r ,m,n z (ρ, ϕ, z) = J n r ( j n r ,m ρ R ) cos(n r ϕ + ϕ 0 ) sin(n z π z H ).
The critical mode λ 0 = 1 corresponds to n r = 0, m = 1 and n z = 1. This yields the criticality condition j 2 0,1

R 2 + ( π H ) 2 = k ∞ -1 L 2 , which implies k ∞ ≥ 1,
and a critical positive flux Φ 0,1,1 (ρ, z) = J 0 ( j 0,1 ρ R ) sin(π z H ).

All others λ eigenvalues corresponds to n z > 0 and either n r = 0 and m > 1 or n r > 0 and m > 0 and we have

δB 2 = j 2 n r ,m -j 2 0,1 R 2 + (n 2 z -1)( π H ) 2 .
We are interested in high λ modes which have opposite variations in opposed quadrants, Φ(ρ, ϕ + π, z) = -Φ(ρ, ϕ, z), and no opposition along z. Hence we must have n r odd and n z odd. The highest λ results from the choice n r = 1, m = 1 and n z = 1 and gives

Γ = 1 R 2 L 2 k ∞ ( j 2 1,1 -j 2 0,1 )
with j 2 1,1j 2 0,1 = 8.899. It follows that for a typical LWR this eigenvalue is very close to 1.

Rectangular prism core

For the simplest case of a bare rectangular prism-shaped core with square cross section the general eigenfunction is Φ(x, y, z) = sin(n x π x a ) sin(n y π y a ) sin(n z π z H ), with the eigenvalue

B 2 = (n 2 x + n 2 y )( π a ) 2 + (n z π H ) 2
and n x , n y , n z > 0 integers. The dominant λ eigenvalue corresponds to n x = n y = n z = 1 and gives the criticality condition

( 2 a 2 + 1 H 2 )π 2 = k ∞ -1 L 2 ,
where for higher modes ones have

δB 2 = (n 2 x + n 2 y -2)( π a ) 2 + (n 2 z -1)( π H ) 2 .
The highest modes with opposite variation at opposite quadrants, Φ(ax, ay, z) = -Φ(x, y, z), require either n x odd and n y even or n x even and n y odd. Therefore the highest mode is a degenerate eigenvalue with two possible states with n z = 1 and either n x = 1 and n y = 2 or n x = 2 and n y = 1. This gives

Γ = 3 L 2 k ∞ ( π a
) .

Validity of approximation (16) in the LFNS

Finally, in view of the results given by simplified diffusion theory, we analyze here the range of frequencies where one can neglect operator P 0,ω in the noise equation. We note that the fact that the fission precursor yields are small as compared to unity leads to the ordering O( P 0,ω ) ∼ ε β (ω/ λ 2 p + ω 2 )O(P 0 ), where typically ε β = O(β p ) ∼ 0.001. Let ε be the order required for the linear noise equation, i.e. O(δψ ω ) ∼ εO(δψ 0 ). This condition entails ω 2 ≤ r 2 λ 2 p /(1r 2 ), where r = ε/ε β . Hence, if ε ≥ ε β we might neglect P as compared to P 0 for all ω's or, otherwise, we can discard P only in the range ω ≤ rλ p / √ 1r 2 .

Fig. 1 .

 1 Fig. 1. Fuel loading map of the studied 3-D core. RE refers to homogeneous heavy baffle assembly, RL to homogeneous quasi-void assembly, UOX_N to homogeneous UOX fuel assembly and UOXGD_N to homogeneous UOX fuel assembly with Gadolinium control rods (with N the number of cycles passed in core).

  (a) Fast steady-state flux. (b) Thermal steady-state flux.

Fig. 2 .

 2 Fig. 2. Steady-state flux of the studied 3-D core (numerical results of the 5th axial plane starting from the bottom).

  (a) Fast steady-state flux. (b) Thermal steady-state flux.

Fig. 3 .

 3 Fig. 3. Steady-state flux of the studied 3-D core in the perturbed assembly in function of the axial position z.

Fig. 4 .

 4 Fig. 4. Moduli of the fast and thermal neutron noises at 3 Hz induced by a cross-sections oscillation in the UOXGD assembly in position (5,7) (numerical results of the 5th axial plane starting from the bottom).

  (a) Fast neutron noise phase. (b) Thermal neutron noise phase.

Fig. 5 .

 5 Fig. 5. Phases (in degree) of the fast and thermal neutron noises at 3 Hz induced by a cross-sections oscillation in the UOXGD assembly in position (5,7) (numerical results of the 5th axial plane starting from the bottom).

  (a) Fast neutron noise modulus. (b) Thermal neutron noise modulus.

Fig. 6 .

 6 Fig. 6. Moduli of the fast and thermal neutron noises at 3 Hz in the perturbed assembly in function of the axial position z.

Fig. 7 .

 7 Fig. 7. Phase shift versus frequency between axial positions z a = 60 cm and z b = 140 cm included in the perturbed assembly. Equation y = -151.66x + 10.37 refers to a linear regression with R 2 = 0.9998 the coefficient of determination.

(a) Fast neutron noise modulus. (b) Thermal neutron noise modulus.

this expansion should also include a contribution from the continuum spectrum which might be important at very high frequencies