Johannes-Geert Hagmann

Lasha Tkeshelashvili

Kurt Busch

Guido Schneider

Far-off-resonant wave interaction in one-dimensional photonic crystals with quadratic nonlinearity

Keywords: numbers: 42.70.Qs, 42.65.-k, 03.50.De

We extend a recently developed Hamiltonian formalism for nonlinear wave interaction processes in spatially periodic dielectric structures to the far-off-resonant regime, and investigate numerically the three-wave resonance conditions in a one-dimensional optical medium with χ (2) nonlinearity. In particular, we demonstrate that the cascading of nonresonant wave interaction processes generates an effective χ (3) nonlinear response in these systems. We obtain the corresponding coupling coefficients through appropriate normal form transformations that formally lead to the Zakharov equation for spatially periodic optical media.

I. INTRODUCTION

During the last two decades, the emergence of periodically structured optical materials -commonly called photonic crystals (PCs) -has lead to substantial progress in the science and technology of optics and photonics. For a recent review, we refer to Ref. [1]. It has soon been realized that PCs with nonlinear constituent materials facilitate the realization of novel effects and functionalities that are very hard to realize with other systems. However, a mathematically consistent description of wave-interaction processes in PCs is complicated due to a number of formal and as yet unsolved problems relating to the complex structure of the equations of motion. Recently, based on a Hamiltonian formulation of Maxwell's equations, an effective field theory for nonlinear PCs has been developed by Volkov and Sipe [2]. This approach allows a systematic study of nonlinear effects beyond the commonly used multiplescales analysis [START_REF] Dodd | Solitons and Nonlinear Wave Equations[END_REF]. In particular, in Ref. [2], the formalism has been applied to the problem of resonant interactions of wave packets in PCs with quadratic nonlinear response. The purpose of this article is to extend this formalism to the far-off-resonant, i.e. nonresonant, case in one-dimensional systems. In particular, we reduce the nonlinear equations of motions of corresponding nonlinear PCs to the simplest form possible. This reduction is facilitated through the systematic elimination of nonresonant interaction terms from the classical Hamiltonian via appropriate normal form transformations. This procedure results in the emergence of an effective third-order response of the system. We would like to note that similar techniques have been successfully employed in the analysis for magnetic systems [START_REF] L'vov | Wave Turbulence under Parametric Excitation: Applications to Magnets[END_REF] as well as for ideal fluid hydrodynamic waves [START_REF] Zakharov | [END_REF]. In the latter case, the result is the so-called Zakharov equation for surface water waves, which has recently been studied experimentally [6].

However, in contrast to the case of spin or water waves, where the underlying linear dispersion relations prohibit nontrivial three-wave resonances, such resonances do appear in PCs for certain sets of interacting wave numbers. This leads to nonremovable singular terms in the normalized Hamiltonian. Consequently, additional constraints on the wave packets are required in order to mathematically justify the validity of the simplified equations for finite times. In realistic physical systems, where losses due to fabricational tolerances and/or residual material absorption, or simply the finite sample size provide natural upper bounds on the time scales, these constraints may be less severe.

The article is organized as follows. In Sec. II, we extend the Hamiltonian formalism of Volkov and Sipe [2] for one-dimensional photonic crystals with χ (2) nonlinearity to the nonresonant case. In Sec. III, we analyze numerically the lowest-order resonance conditions and compare the results to the water wave problem for which the spectrum is known analytically. In Sec. IV, the normal form reduction of the Hamiltonian is performed and the Zakharov equation for nonlinear PCs is formally derived for the case when three-wave interaction processes are suppressed. In particular, the expressions for the effective third-order nonlinear response are derived and compared with the case of homogeneous dispersive χ (2) materials. We conclude with a discussion of the limitations of our approach, and provide an outlook on possible future work in Sec. V.

II. CLASSICAL HAMILTONIAN FORMULATION

In this section we derive a classical Hamiltonian formulation for nonlinear PCs. The derivation closely follows Ref. [2], however, we emphasize that, as compared to Ref. [2], we have chosen a different set of fundamental fields for the decomposition into a complete set of eigenfunctions.

We consider infinite nonlinear PCs with nonmagnetic, dispersion-and lossless constituent materials. Since the linear dispersion provided by photonic crystals' microstructure is typically much larger than the material dispersion of its constituent materials, this is a very good approximation. If material dispersion dominates -such as would be the case for quasiphase matched gratings where only the nonlinear properties vary periodically in space -one would have to resort to a slowly varying envelope approximation [7]. Assuming further that the material response is local in space and time, we can express the polarization P(x, t) that connects the dielectric displacement D(x, t) = ǫ 0 E(x, t)+P(x, t) with the electric field E(x, t) as

D i (x, t) = ǫ 0 E i (x, t) + χ (1) ij (x)E j (x, t) + χ (2) ijk (x)E j (x, t)E k (x, t) + . . . . (1)
Here, ǫ 0 denotes the free-space dielectric permittivity. The susceptibility tensors χ (n) of rank n + 1 are symmetric under permutations of their indices and -here as well as in the remainder of the article -summation over repeated indices is implied. Instead of the first-order susceptibility χ (1) , we employ the dielectric tensor ǫ(x). In addition, in the following, we will restrict ourselves to the case of isotropic linear response, so that we have χ (1) (x) = ǫ 0 (ǫ(x) -1). Throughout this paper, we identify the canonical momentum with the displacement field D rather than with the electric field E. Therefore, it is desirable to expand the electric field in terms of the displacement field according to

E i (x, t) = ǫ -1 0 ǫ(x) -1 D i (x, t) + Γ (2) ijk (x)D j (x, t)D k (x, t) + . . . . (2)
For the second-order susceptibility, the relations between χ (2) and Γ (2) reads

Γ (2) ijk (x) = - χ (2) ijk (x) (ǫ 0 ǫ(x)) 3 , (3)
and the relations for the higher-order susceptibilities can be found by recursion [8]. In PCs, the material properties are periodic functions with respect to the set of lattice vectors R, i.e., we have ǫ

(R + x) = ǫ(x), χ (n-1) i1i2...in (R + x) = χ (n-1) i1i2...in (x), Γ (n-1) i1i2...in (R + x) = Γ (n-1)
i1i2...in (x), etc. Within the Dzyaloshinski or generalized Coulomb gauge where the scalar potential φ(x, t) ≡ 0, we may express both, the electric field E and the magnetic field B, through the vector potential A according to

E(x, t) = -∂ t A(x, t), (4)
B(x, t) = ∇ × A(x, t).

(5)

For this system, the Hamiltonian that yields Maxwell's equations in the absence of external sources and currents reads

H(A, Π) = 1 2 d 3 x 1 ǫ 0 ǫ(r) D i (x, t)D i (x, t) + µ -1 0 B i (x, t)B i (x, t) + 1 3 d 3 xΓ (2) ijk (x)D i (x, t)D j (x, t)D k (x, t) + . . . = H 2 + H 3 + . . . . (6)
In the former expression, Π = -D denotes the canonical momentum associated with the vector potential A [8]. Consequently, the equal-time Poisson brackets are

{A i (x, t), Π j (x ′ , t)} A,Π = d 3 x ′′ δA i (x, t) δA k (x ′′ , t) δΠ j (x ′ , t) δΠ k (x ′′ , t) - δA i (x, t) δΠ k (x ′′ , t) δΠ j (x ′ , t) δA k (x ′′ , t) , = δ ij δ(x -x ′ ), (7)
where δ/δA i (x, t) and δ/δΠ i (x, t) denote functional derivatives (see Ref. [9] for a definition and properties). Theon first sight somewhat superfluous -indices A and Π on the left-hand side allow us to keep track of the fields, once a transformation is applied to the Hamiltonian (6). To decompose this Hamiltonian into symmetry-adapted basis functions (modes), it is useful to consider the solutions of the linearized equations of motion for periodic systems.

A. Bloch function decomposition

In the linear limit, D(x, t) = ǫ 0 ǫ(x)E(x, t), and the Maxwell equations reduce to two decoupled partial differential equations for the fields A(x, t) and D(x, t). Since we are interested in stationary solutions (modes), we employ the time-harmonic ansatz A(x, t) = A nk (x) exp(-iω nk t), D(x, t) = D nk (x) exp(-iω nk t), and obtain from the Maxwell equations

1 ǫ(x) ∇ × (∇ × A nk (x)) = ω 2 nk c 2 A nk (x), ∇ × ∇ × D nk (x) ǫ(x) = ω 2 nk c 2 D nk (x). (8)
In the above ansatz, we already made use of the Floquet-Bloch theorem that -in the case of lattice periodic functions ǫ(x) -applies to Eqs. (8) and, therefore, introduced the Bloch functions A nk (x) and D nk (x). These Bloch functions are labeled through the discrete band index n > 0 and the continuous wave vector k ∈ Ω b that lies in the first Brillouin zone Ω b . This labeling is commonly referred to as the reduced zone scheme [START_REF] Ashcroft | Solid State Physics[END_REF]. In general, the Bloch functions and their associated eigenvalues ω nk have to be determined numerically through one of the well-documented techniques [1]. For future reference, we would like to note that the Bloch functions have the forms of modulated plane waves

D nk (x) = dnk (x) exp(ikx), (9)
A nk (x) = ãnk (x) exp(ikx), (10)
where ãnk (x) and dnk (x) are lattice-periodic functions. For convenience, the Bloch functions D nk , A nk and the eigenvalues are chosen to be periodic in reciprocal lattice vectors G, i.e. D n(k+G) = D nk for all reciprocal lattice vectors G [START_REF] Lax | Symmetry principles in solid state and molecular physics[END_REF]. It is also useful to note that the time-reversal symmetry of our problem implies ω nk = ω n-k which, in turn, implies that together with the solutions (D nk , ∇ × A nk ), the corresponding pair (D ⋆ n-k , -∇ × A ⋆ n-k ) also fulfills the Maxwell equations for the same frequency. Therefore, we choose D nk = D ⋆ n-k which guarantees that the displacement field D(x, t) is real. The Bloch functions constitute a complete orthonormal set of basis functions. Since Eq. ( 4) implies

D nk (x) = iǫ 0 ǫ(x)ω nk A nk (x), (11)
we find the the orthogonality relations [10]

d 3 x ǫ 0 ǫ(x) ω nk A nk (x)A ⋆ n ′ k ′ (x) = 1 2 δ nn ′ δ(k -k ′ ), (12)
d 3 x (∇ × A nk (x)) • (∇ × A ⋆ n ′ k ′ (x)) = µ 0 ω nk 2 δ nn ′ δ(k -k ′ ), (13)
for k, k ′ ∈ Ω B if n = 0. For n = 0, k = 0, we have ω nk = 0 and the normalization has to be performed separately. This issue will be addressed below.

For the subsequent formal developments, we find it more transparent to represent the Bloch functions in the socalled extended zone scheme [START_REF] Ashcroft | Solid State Physics[END_REF], where the wave vector k varies over the entire reciprocal space. Formally, this is facilitated by defining the band index to be a function of the wave vector, i.e., to introduce

Āk (x) = A n(k)k (x), (14)
and an analogous definition for Dk and ωk . Whether such a mapping does or does not exist depends on the spectrum for the selected system, in particular on the occurrence of degeneracies. For one-dimensional PCs and wave propagation along the axis of periodicity, wave vector and Bloch functions reduce to scalars. As a result, the unfolding of the band structure is realized through n(k) := ⌊|k|d/π⌋ ∈ N, where ⌊•⌋ denotes the Gauss bracket defined via 0 < |k| -(π/d)n(k) ≤ (π/d) (d denotes the lattice constant of the PC). In two-or three-dimensional systems, this transformation between reduced and extended zone scheme is not straightforward but can, in principle, be carried out numerically. Therefore, we in the following restrict ourselves to the case of strictly one-dimensional PCs and assume that the system is periodic in the x direction. To simplify the notation, we assume in the following an isotropic nonlinear response of the constituent materials so that the relevant fields (A, D, and E) are all polarized along a certain transverse direction which we define as the y direction. As a result, we may suppress the vectorial character of all fields. In such a one-dimensional PCs, Eqs. ( 11) -( 13) read as

Dk (x) = iǫ 0 ǫ(x)ω k Āk (x), ( 15
)
dx ǫ 0 ǫ(x) ωk Āk (x) Ā⋆ k ′ (x) = 1 2 δ(k -k ′ ), ( 16
) dx ∂ x Āk (x) ∂ x Ā⋆ k ′ (x) = µ 0 ωk 2 δ(k -k ′ ), (17)
for k, k ′ ∈ R\{0}.

We now return to the above-mentioned normalization problem related to ωk=0 = 0 for which the eigenvalue equation reads as

∂ xx D0 (x) = 0. ( 18
)
Owing to the Bloch-Floquet condition and our choice of Dk = D⋆ -k , we obtain D0 (x) = a ∈ R. On the other hand, Dk must be a continuous function of k so that the scaling Dk ∝ √ ωk requires us to define

Dk=0 = 0, ( 19
) lim µ→0 Dµ ∝ lim µ→0 √ µ. (20)
Therefore, an absolute normalization in the sense of the orthogonality relations is impossible, yet our definition is justified from the physical point of view [START_REF] Ford | [END_REF]. In fact, physically the mode ωk=0 = 0 corresponds to a simple translation of our infinite PC and this clearly requires an infinite amount of energy. Consequently, this mode resists an absolute normalization. The above definition ensures that this is the case. Clearly, this discussion may be directly transferred to the case of two-and three-dimensional PCs. Finally, we expand the physical fields D(x, t) and A(x, t) of the nonlinear system into the orthogonal basis of the linearized problem

D(x, t) = R dk Dk (x)ā k (t) + c.c., (21)
A(x, t) = R dk Āk (x)ā k (t) + c.c., (22)
where the time-dependent amplitudes āk (t) and ā⋆ k (t) are the weights associated with each Bloch mode. In terms of the physical fields, the new set of variables (k = 0) are expressed as

āk (t) = dx 1 ωk D⋆ k (x)D(x, t) ǫ 0 ǫ(x) + ǫ 0 ǫ(x) ωk Ā⋆ k (x)A(x, t) . (23)
In the derivation of Eq. ( 23), the identity ( Dk , ∂ x Āk ) = ( D⋆ -k , -∂ x Ā⋆ -k ) has been employed and it can be verified that Eq. ( 23) represents a canonical transformation. Consequently, for arbitrary functions F, G in phase space,

{F, G} D,A = {F, G} ā,ā ⋆ . ( 24
)
holds. In terms of the new variables, the equal-time Poisson brackets now have the explicit form

{F, G} ā,ā ⋆ = i R dk δF δā k δG δā ⋆ k - δG δā k δF δā ⋆ k , (25)
so that

{ā k , ā⋆ k ′ } ā,ā ⋆ = iδ(k -k ′ ), ( 26
) {ā k , āk ′ } ā,ā ⋆ = 0, (27)
is immediate.

B. Hamiltonian in the new variables

The results of the previous section may now be utilized to rewrite the original Hamiltonian (6) in terms of the new field variables āk and ā⋆ k . The linear part, H 2 , is a quadratic form

H 2 = dk ωk āk ā⋆ k , (28)
while the lowest-order nonlinear term, H 3 , couples three Bloch modes and reads as

H 3 = dk 123 V (1) k1k2k3 āk1 āk2 āk3 G δ(G -k 1 -k 2 -k 3 ) + dk 123 V (2) k1k2k3 ā⋆ k1 āk2 āk3 G δ(G + k 1 -k 2 -k 3 ) + c.c. (29)
= H

(1) 3

+ H

(2) 3

+ c.c.

In the former expression, the integration is carried out over all three momenta, i.e., dk 123 = dk 1 dk 2 dk 3 and the sum over the entire set of reciprocal lattice vectors G ∈ {G ∈ R : G = m(2π/d) , m ∈ Z} reflects the fact that, in the extended zone description, in principle all bands are coupled. In the reduced zone scheme, this would include both the so-called Umklapp processes, i.e., when the sum of the coupled wave vectors leaves the Brillouin zone as well as the coupling between different bands without Umklapp processes. Indeed, this complication of classifying and keeping track of the various scattering processes in the reduced zone scheme has lead us to adopt the (perhaps less familiar) extended zone scheme. The corresponding coupling coefficients V

k1k2k3 and V

(2

)
k1k2k3 are derived from the Bloch functions as given in Eqs. ( 9) and ( 10) (transferred to the extended zone notation) together with the identity

R exp (-ikR) = 2π d G δ(G + k), (30)
that allows one to reduce the integration over all space into an integration over the PC's unit cell

V cell V (1) k1k2k3 = 8π 3 3d V cell dx Γ (2) (x) Dk1 (x) Dk2 (x) Dk3 (x) , (31)
V (2) k1k2k3 = 8π 3 3d P c V cell dx Γ (2) (x) D⋆ k1 (x) Dk2 (x) Dk3 (x). (32)
In this expression, we have introduced the cyclic permutation operator P c that acts on any function F 123 with indices 123 according to P c F 123 = F 123 + F 213 + F 321 . Physically, the terms associated with V

k1k2k3 correspond to the simultaneous creation or annihilation of three waves, while the terms associated with V

(2) k1k2k3 describe wave mixing processes where a single wave is converted into two other waves. Clearly, a process associated to the matrix element V

(1) k1k2k3 does not conserve energy. In order to be efficient, a coupling mechanism to another physical system is required. Therefore, for conservative optical systems these processes are strongly suppressed.

The explicit form of the higher-order nonlinear contributions to the Hamiltonian, H j , j > 3, follow in strict analogy to the derivation of H 3 . For instance, the next nonlinear term H 4 contains both energy-suppressed and energetically allowed processes. Among all of those, the relevant part of the Hamiltonian corresponding to the third-order nonlinear susceptibility is

H (3) 4 = dk 1234 W (3) k1k2k3k4 āk1 āk2 ā⋆ k3 ā⋆ k4 G δ(G + k 1 + k 2 -k 3 -k 4 ), (33)
where processes in which two waves are scattered are described via the coupling coefficient

W (3) k1k2k3k4 = π 2d P 6 V cell dx Γ (3) (x) Dk1 (x, t) Dk2 (x, t) D⋆ k3 (x, t) D⋆ k4 (x, t). (34)
Here, the permutation operator P 6 acts on a function with indices F 1234 according to

P 6 F 1234 = F 1234 + F 1324 + F 3214 + F 1432 + F 4231 + F 4321 .
The equations of motion for the canonical variables are the Hamilton equations

ȧk (t) = {ā k (t), H} ā,ā ⋆ . ( 35
)
However, it is often preferable to work in a frame that "rotates" with the frequency of the free carrier wave, i.e., to employ the slow variables ᾱk (t) = āk (t) exp(-iω k t). In this case the equations of motions read

∂ t ᾱk = i dk 23 V (2) kk2k3 ᾱk2 ᾱk3 e i(ω k 2 +ω k 3 -ω k )t G δ(G + k -k 2 -k 3 ) (36) +{ ᾱk (t), H -H (2) 3 } ā,ā ⋆ (37)
The terms on the right-hand side of the equation of motion give an oscillatory contribution to the mode ᾱk unless certain resonance conditions are met. For the one-dimensional photonic crystals, these conditions are analyzed and compared to a related problem in ideal hydrodynamics in the following section.

III. THREE-WAVE INTERACTION

In many problems of interest in optics, the nonlinearities that are involved in the system dynamics are quite small. This suggests that in the equations of motion the nonlinear interaction terms can be considered as a perturbation to the linear part. The resonance conditions for nonlinear wave mixing processes are determined by the linear dispersion relation ω k . A profound example of such a resonant mixing process in optics is the second harmonic generation, where the frequency of an incident beam is doubled through a quadratic, i.e., χ (2) -nonlinear response of the system. The requirement for efficient conversion is the conservation of energy and the momentum in the process, ω + ω = 2ω and k + k = 2k, respectively. These equations are examples of the more general class of three-wave resonance conditions

s 1 ω(k 1 ) + s 2 ω(k 2 ) + s 3 ω(k 3 ) = 0 (38) s 1 k 1 + s 2 k 2 + s 3 k 3 , = 0 ( 39
)
where s i = ±1. These conditions are the requirements for an efficient energy transfer between three waves (modes) of wave numbers k 1 , k 2 , k 3 . For PCs, the second condition must be replaced by the conservation law for crystal momentum, i.e., reciprocal lattice vectors can be added to any of the momenta. For the moment, we leave this point aside but will return to it at a later stage of the discussion.

A detailed discussion of selection rules for nonlinear PCs can be found in Ref. [14]. We would like to note that in the case of non-negligible third-order response, i.e., χ (3) nonlinearities, similar conditions must hold for four-wave mixing processes. At first sight, it appears difficult to solve Eqs. (39) even if the linear spectrum ω k is known analytically such as in the case of surface waves in an ideal fluid. On the other hand, for fixed s i , Eqs. (39) contains six variables but only three of them can be chosen independently, and the two equations remove one degree of freedom. In this reduced parameters space and for any given ω k0 , one can easily carry out the analysis numerically and visualize the resonance conditions.

Upon combining the expressions in Eqs. (39), we obtain

ω(-s 1 s 2 k 2 -s 1 s 3 k 3 ) + s 1 s 2 ω(k 2 ) + s 1 s 3 ω(k 3 ) = 0 (40)
Equation ( 40) is invariant under the simultaneous inversion of the signs of all s i , so that only four out of eight cases have to be considered. Time reversal symmetry, i.e., ω(k) = ω(-k) ≥ 0, and the fact that in one-dimensional PCs the spectrum is nondegenerate, lead to further reduction. In fact, a semi-positive spectrum excludes the case of three identical signs and the time-reversal symmetry maps two other cases onto each other up to an exchange of indices. Consequently, in order to find nontrivial solutions to Eq. (40), we are left with only two distinct cases. These can readily be visualized by defining the two functions

a 1 (k 2 , k 3 ) = ω(k 2 -k 3 ) -ω(k 2 ) ω(k 3 ) + 1 (41) a 2 (k 2 , k 3 ) = ω(k 2 + k 3 ) -ω(k 2 ) -ω(k 3 ) + 1 (42) for k 2 ∈ R, k 3 ∈ R \ {k 0 = 0}, (k 3 = k 0 only yields a trivial resonance)
. With these definitions, the resonance conditions are fulfilled if and only if a i (k 2 , k 3 ) = 0. The first spectrum ω k for which we evaluate Eq. ( 42) is the dispersion relation for gravity-capillary waves on the surface of an ideal fluid with infinite depth

ω k = g|k| + σ|k| 3 , ( 43
)
where g is the gravitational acceleration and σ the capillarity of the fluid. Introducing k ′ = σg -1 k (see also the right panel of Fig. 1), transfers the previous expression into

ω k ′ σ 1 4 g -3 4 = |k ′ | + |k ′ | 3 (44)
for evaluation in a 1 (k 2 , k 3 ), a 2 (k 2 , k 3 ). Notice that both ω k ′ σ 1 4 g -

i (k ′ 2 , k ′ 3 ) = a i (-k ′ 2 , -k ′ 3 ). Besides the trivial resonances k ′ 2 = k ′ 3 for a 1 (k ′ 2 , k ′
3 ) (Fig. 1) and k ′ 2 = 0 for a 1 (k ′ 2 , k ′ 3 ) (Fig. 2), there are other sets of resonant and nearly resonant pairs of wave vectors that are located on ribbons. It is instructive to notice the influence of the capillary constant on the position of these ribbons by a scaling argument: Selecting fixed (k ′ 2 , k ′ 3 ) in the set of nearly resonant wave numbers, the location of the point moves in the direction of increasing (k 2 , k 3 ) with decreasing σ such that in the limit σ → 0 these resonances lie in the infinite for (k 2 , k 3 ). For pure gravity waves (σ = 0), there are only trivial resonances when one of the wave vectors involves the dc component (ω 0 = 0). The first nontrivial mixing processes occur at higher order (four-wave processes).

Returning to three-wave mixing in PCs, the situation is very different from the previous case. In contrast to water waves, there is no explicit expression for ω k and the discussion is necessarily qualitative. In order to illustrate the situation, we assume that the PC consists of alternating layers with equal thickness but different refractive indices, n A and n B , respectively. Then, the only parameter that alters the spectrum is the refractive index contrast ∆n = n A -n B . To calculate the photonic band structure, ω k , we use a plane-wave expansion method. Of course, any other method for obtaining photonic band structures could be used [1]. Except for the case of a homogeneous dielectric, the spectrum takes the form of a band structure, as shown in right panel of Fig. 3. An example with an index contrast of ∆n = 2 is shown in Figs. 3 and4 the borders of successive Brillouin zones appear at half-integer values of the dimensionless variable. As in the case of water waves, the existence of ω 0 = 0 leads to trivial resonances in the diagonal of Fig. 3 and a horizontal line in Fig. 4. In addition, there are isolated resonant and nearly resonant subsets of wave vectors; their number and relative position depends on the refractive index contrast. However, it is hard to quantify this dependence. In the extreme case of a homogeneous material, ∆n = 0, the resonance conditions can trivially be fulfilled. For instance, for a 2 (k 2 , k 3 ) the two squares k 2 ≥ 0, k 3 > 0 and k 2 ≤ 0, k 3 < 0 will become identically zero. With increasing contrast, photonic band gaps form at the borders of the first Brillouin zone and the deviation from a linear spectrum is strongest in this region. This behavior is reflected by the pattern in the left panel of Fig. 4. As alluded to above, it is important to recall that in PCs only crystal momentum is conserved, while the a i (k 2 , k 3 ) in Figs. 3 and4 have been calculated on the basis of momentum conservation. Therefore, for a complete interpretation all points in the plane (k 2 , k 3 ) have to be shifted through the reciprocal lattice vector G = m2π/d, m ∈ Z, both vertically and horizontally. Consequently, there are more resonant zones than those shown in Figs. 3 and4 that have to be considered.

Although at this point, we have not yet formally constructed the normal-form Hamiltonian for the PC, we can already anticipate an important difficulty in its derivation. The resonance conditions, Eqs. (40), will appear in the denominators of the effective interaction coefficients. Whereas singularities resulting from trivial resonances are removable, the singular terms associated with nontrivial resonances are pertinent. Therefore, additional constraints on the distribution of wave vectors are required in order to ensure the applicability of the normal form transformation. Since a wave packet in a nonlinear medium generally acquires new wave vector components as time goes on, these constraints will translate into constraints on the time scales for which our theory will be applicable.

IV. REDUCTION TO NORMAL FORM

Normal form transformations represent a particularly useful tool for the study of nonresonant wave interaction processes in Hamiltonian systems. In essence, they facilitate the reduction of nonlinear dynamical systems to the simplest possible form by eliminating inessential terms or degrees of freedom from the original Hamiltonian. For physical systems with continuous degrees of freedom, a good illustration of this principle is the study of the aforementioned water wave problem which dates back to the early work of Zakharov [START_REF] Zakharov | [END_REF]. More recent highlights include the proof of the nonintegrability of the water wave problem with the help of normal form transformations by Craig and Worfolk [15,16].

In this section, we employ normal form transformations to the case of nonresonant wave interaction in χ (2) PCs. The goal is derive an effective χ (3) response of the system that originates from a cascading of χ (2) processes. This analysis shows how the engineering of the PC structure could, in principle, be used to tune the total χ (3) response of the system.

In our analysis, we have to consider boundary conditions which lead to uncountably many modes. In addition, in PCs the resonant set of arguments cannot be unambiguously separated from the nonresonant ones since the nontrivial resonances have to be determined numerically and the influence of the refractive index contrast on their positions cannot be quantified (see Sec. III). Similar to the cases of finite degrees of freedom with purely imaginary spectra of incommensurable eigenvalues (depending on whether or not un-removable resonances remain, the so-called Birkhoff of Birkhoff-Gustavson normal forms) we utilize -for our case of uncountably many degrees of freedom -normal forms as a formal tool to simplify the Hamiltonian and, hence, the equations of motion, but leave aside certain questions regarding the convergence of the transforms.

A. Normal form transformation -Lie series

In view of the discussion pertaining the efficiency of wave interaction processes in conservative systems (see Sec. II B), we assume in the following that the initial conditions, i.e. the initial distribution of modes a k (0), a ⋆ k (0) in the PC with quadratic nonlinearity, are such that the resonance condition of lowest order that can be satisfied is a four-wave mixing process of the form

ω k1 + ω k2 -ω k3 -ω k4 = 0. (45)
All other resonance conditions that may hold are assumed to be of order 5 or higher. This four-wave resonance condition, Eq. ( 45), is trivial in the sense that it can be fulfilled for any kind of dispersion relation and all possible initial conditions, even if the spectrum does not include the frequency ω = 0. Consequently, the system Hamiltonian H = H 2 + H 3 + H 4 + • • • comprises contributions in H 3 and H 4 that lead to oscillatory source terms in the equation of motions for the time-dependent amplitudes a k (t), a ⋆ k (t), as long as these modes are bound away from the resonant subsets. To simplify the structure of these equations, we search for a coordinate transformation that eliminates all terms in the Hamiltonian that correspond to resonance conditions which cannot be fulfilled. In constructing such a transformation, we follow the ideas outlined in [15,16,17].

In order to split the Hamiltonian into resonant and nonresonant parts, it is useful to introduce the linear mapping

ad F (G) = {F, G}, (46)
where {F, G} is the Poisson bracket of two functions F, G in phase space. In quantum optics, the operator ad F acting on any function G is known as the Liouville operator. As an illustration, let us consider the Hamiltonian presented in Sec. II. It can be easily verified that for H

(1) 3

the following relation holds:

ω k1 + ω k2 + ω k3 = 0 ∀(k 1 , k 2 , k 3 ) ∈ R 3 =⇒ ad H2 (H (1)
3 ) = 0. More generally, ad H2 splits the Hamiltonian into resonant H res m and nonresonant H nr m terms of order m, where H res m ∈ ker(ad H2 ) = {F m : ad H2 (F m ) = 0}. Now, our goal is to eliminate all nonresonant terms from the Hamiltonian and to replace the initial system containing resonant and nonresonant terms with a simpler system. In particular, we want to find a canonical (symplectic) transformation āk = F ( bk ), such that in the new variables bk of the transformed Hamiltonian H ′ (b) takes on the form

H ′ (b) = H 2 (b) + H ′ res 4 (b) + R ′′ 5 (b), (47)
where the residual R ′′ 5 (b) contains only processes of fifth or higher order in the amplitudes b k . Such reductions can be performed by means of the Lie series transformations, which in the nonlinear optics community received attention through the concept of the guiding center soliton in amplified fiber lines [START_REF] Hasegawa | Solitons in optical communications[END_REF]. We define the Lie series exp(ad K ) induced by the generating functional K as

exp(ad K )(H) = ∞ m=0 1 m! ad m K (H). (48)
Here, K is a minimally third-order functional in the dynamic variables. For the purposes of the present work, only a few of the remarkable properties of these series will be used. Among them are [17,[START_REF] Gröbner | Die Lie-Reihen und ihre Anwendungen[END_REF]:

exp(ad K ){F, G} = {exp(ad K )F, exp(ad K )G} (49) (exp(ad K )) -1 = exp(ad -K ), (50) exp(ad K )F (u) = F (exp(ad K )u), (51)
where F, G are functions in phase space and u are the phase space variables. The first equality establishes the property of being a symplectic transformation for arbitrary generating functional K and can be verified by using the Jacobi identity for the Poisson bracket. The second equality makes a connection to the inverse Lie series transformation. Finally, the last equation allows us to interchange the order of mapping.

The explicit form of the polynomial K can be determined by additional requirements on the normalization of the Hamiltonian. Writing K = K 3 + K 4 , where K 3 , K 4 , respectively, are -yet to be determined -functionals of order 3 and 4 in the dynamical variables, and expanding the exponential yields

exp(ad K )H = exp(ad K ) [H 2 + H 3 + H res 4 + H nr 4 + R 5 ] = H 2 + H 3 + H res 4 + H nr 4 + R 5 + {K, H 2 } +{K, H 3 } + 1 2 {K, {K, H 2 }} + R ′ 5 . (52)
This suggests that the choice

K 3 = ad -1 H2 (H 3 ), (53)
K 4 = ad -1 H2 [H 4 + 1 2 ad H3 (K 3 )] nr (54)
realizes the desired simplification, i.e., we obtain exp(ad

K )H = H 2 + H ′ res 4 + R ′′ 5 . (55)
In Eq. ( 54), [•] nr indicates that we only consider the nonresonant terms contained within the bracketed expression. Alternatively, we can apply a transformation from the original set of coordinates āk , ā⋆ k to a new set bk , b⋆ k according to

a k = exp(ad K(b) ) b k , ( 56
)
where K is the functional defined in Eq. ( 54) and we have -and will continue to do so in the remainder of this paper -dropped the bars over the coordinates. More explicitly, we choose

K 3 = dk 123 A (1) k1k2k3 b k1 b k2 b k3 G δ(G -k 1 -k 2 -k 3 ) (57) + dk 123 A (2) k1k2k3 b ⋆ k1 b k2 b k3 G δ(G + k 1 -k 2 -k 3 ) + c.c., (58)
as well as

K 4 = dk 1234 B (1) k1k2k3k4 b k1 b k2 b k3 b k4 G δ(G -k 1 -k 2 -k 3 -k 4 ) (59)
+ dk 1234 B

(2)

k1k2k3k4 b ⋆ k1 b k2 b k3 b k4 G δ(G + k 1 -k 2 -k 3 -k 4 ) + c.c..
One can note that since ({F, G} a,a ⋆ ) ⋆ = {F ⋆ , G ⋆ } a,a ⋆ , the requirement K = K ⋆ follows immediately from the definition of the coordinate transformation. The unknown coefficients A (i) in K 3 are determined from the nonresonance condition, Eq. ( 54), and are readily found to be

A (1) k1k2k3 = i V (1) k1k2k3 ω k1 + ω k2 + ω k3 , (60)
A (2) k1k2k3 = i V (2) k1k2k3 -ω k1 + ω k2 + ω k3 . (61)
Equation ( 54) also fixes the coefficients B (i) , which can be expressed in terms of V (i) and A (i) , i = 1, 2. However, their explicit form is not required if we are only interested in terms at fourth order to which K 4 does not contribute. Finally, we note that the expressions for A (i) become singular if the resonance condition is fulfilled.

In these new coordinates, the transformed Hamiltonian reads

H ′ (b) = H 2 (b) + H ′ res 4 (b) + R ′′ 5 (b) = dk ω k b k b ⋆ k + dk 1234 S k1k2k3k4 b k1 b k2 b ⋆ k3 b ⋆ k4 G δ(G + k 1 + k 2 -k 3 -k 4 ) +R ′′ 5 (b). ( 62
)
With the help of the corresponding Poisson brackets

{b k , b ⋆ k ′ } = iδ(k -k ′ ), (63) {b k , b k ′ } = 0, (64)
we finally obtain the transformed equations of motion

d dt b k = iω k b k +i dk 123 (S k1k2k3k + S k1k2kk3 )b k1 b k2 b ⋆ k3 G δ(G + k 1 + k 2 -k 3 -k) + R4 (b), (65)
where R4 (b) is a residual containing higher-order terms. Equation (65) represents the generalization of the Zakharov equation to the case of periodic media. This equation, which Zakharov reduced in his original article [START_REF] Zakharov | [END_REF] to the Nonlinear Schrödinger equation (NLSE), has been the subject of extensive studies for waves in homogeneous systems and has been tested experimentally for the case of surface water waves in channels [6].

B. Effective fourth-order coefficient

Using the coefficients A (i) , we now proceed to calculate {K(b), H 3 (b)} and its contribution to the transformed Hamiltonian, which contains the renormalized fourth-order matrix element. Indeed, the effective fourth-order element S k1k2k3k4 consists of the term originating from the third-order nonlinear susceptibility χ (3) of the system and three additional terms N (i) , i = 1, 2, 3 that arise from nonresonant χ (2) processes

S k1k2k3k4 = W (3) k1k2k3k4 + 1 2 N (1) k1k2k3k4 + N (2) k1k2k3k4 + N (3) k1k2k3k4 (66)
waves with the same frequency via the static field ω 0 or via waves at symmetry points ω G . In the field of nonlinear optics of dispersive homogeneous media, effective higher-order nonlinear contributions due to phase mismatch are generally referred to as cascaded processes. These effects are known to effectively modify the nonlinear third-order response of homogeneous materials with noncentrosymmetric nonlinear susceptibilities [START_REF] Ch | [END_REF]. The phase mismatch in these materials results from intrinsic material related dispersion, which -due to Kramers-Kronig relations -always appears together with absorption. Although some progress has recently been made on the Hamiltonian description of the optical properties of nonlinear dispersive materials [21], the Lie series transformation does not directly apply to this class, because the detailed knowledge of interaction coefficient relies on a modal description. The key aspect about PCs which enable the formal construction of a normal form is that nontrivial dispersion relations result from the formation of a photonic band structure. Consequently, if the validity of the equation could be justified for finite times with an appropriate error bound, the result suggest that it is possible to achieve cascading without using strongly dispersive and, therefore, lossy absorbing materials. In addition, the judicious engineering of PCs would allow to tailor the effective third-order nonlinear response of these systems through both the dispersion relation ω k and the matrix elements V

(1) k1k2k3 and V

k1k2k3 . Taken together, these two aspects enlarge the possibilities for the realization of optical materials with tailored nonlinear properties.

Although some results on the validity for finite times of envelope equations that have been derived under nonresonance assumptions in the presence of resonances are known [22], these turn out to be inapplicable to the present case, as outlined in the following section. for k, m → 0, the system truncated after the third order terms cannot be expected to be well posed. In detail, we have ω k ∼ k plus some correction terms coming from the periodic coefficients. It is well known [START_REF] Eastham | The Spectral theory of periodic differential equations[END_REF] that the corrections decay asymptotically with the regularity of the associated coefficients. If the coefficients are n-times differentiable, then for suitably chosen k and m, we obtain r(k, m) ≤ C |k| -n + |k -m| -n + |m| -n .

(71)

Since r(k, m) appears in the denominator in the normal form transformation, there is a substantial loss of regularity associated with the normal form transformation. Thus, the less regular the coefficients, the better are the conditions for well posed truncated third-order systems.

The same problems occur if the truncated third-order system is justified with some approximation result similar to that of Ref. [START_REF] Busch | [END_REF]. In this work, the NLSE has been justified in the above sense for semi linear wave equations with periodic coefficients. An additional difficulty occurs in our situation due to the quasilinearity of the problem that always leads to some loss of regularity in normal form transformations.

However, on physical grounds (see Sec. III) we expect that for suitably chosen initial conditions, the dynamics of the system is dominated by nonresonant processes for not too long times. Quantifying this conjecture by methods of asymptotic analysis is a challenging mathematical problem that has not been addressed yet.

V. CONCLUSIONS

In this article, we have extended a recently developed Hamiltonian formalism to the analysis of nonresonant wave interaction in one-dimensional PCs with quadratic nonlinearity. In particular, we have shown that an appropriate normal form transformation leads to a generalization of Zakharov's equation for periodic media. For the simplest case of monochromatic wave, the resulting effective coupling coefficient in the transformed equations of motion is the result of cascading second-order nonlinear processes. Moreover, we derived explicit expressions for the emerging effective nonlinear coupling terms. This may allow one to tailor the effective nonlinear response of the composite materials through judiciously engineering of its linear properties and could, in turn, suggest possibilities for the realization of optical materials with customized nonlinear properties. The approach presented here can also be applied to nonresonant higher-order interaction processes. Owing to the linear dispersion relation in the long-wavelength limit certain mathematical difficulties arise that limit the applicability of the method to finite times. However, we would like to emphasize that in experiments, various loss mechanisms or the finite sample size set maximal time scales so that this limitation may be less severe. In addition, full numerical simulations can provide further insight into the range of validity of our theory. A systematic study of the Zakharov equation's validity represent a challenging mathematical task and further research is indicated.

3 4

 3 FIG. 1: Left panel : |a1(k ′ 2 , k ′ 3 )| for gravity-capillary waves; Right panel: Corresponding dispersion relation ω k ′ . All variables are dimensionless, for notation and normalization see the text.

  FIG. 3: Left panel : |a1(k2, k3)| for a PC with refractive index contrast ∆n = 2; Right panel: Corresponding photonic band structure ω k in the extended zone scheme.

  D. Influence of resonancesAs a consequence of the strong asymptotic resonancer(k, m) = ω(k) -ω(k -m) -ω(m) → 0 (70)

Acknowledgements: JGH, LT, and KB acknowledge support from the DFG-Forschungszentrum Center for Functional Nanostructures (CFN) at the Universität Karlsruhe within project No. A1.2. The research of LT and KB is further supported through the DFG-Priority Program No. SPP 1113 Photonic Crystals within project Bu 1107/6-1. LT also acknowledges support from USA CRDF Grant No. GEP2-2848-TB-06. Finally, the authors would like to thank J. Niegemann, L. Shemer, and H. Uecker for fruitful discussions, and S.N. Volkov and J.E. Sipe for helpful correspondence.

where the explicit expressions for the N (i) are

In the above expressions, the denominators of

k1k2k3k4 contain the resonance conditions. In particular, frequencies which depend on the difference between two wave vectors appear. Note that S k1k2k3k4 describes coupling of the four waves b ki , and therefore, trivial resonances appear in the denominators for any possible distribution of wave vectors in the wave packets b ki . As a consequence, S k1k2k3k4 is ill defined, unless the singularity is removable. To analyze this issue in the next section, we further simplify the matrix elements by considering the simplest possible situation, i.e., the propagation of a monochromatic wave.

C. Simple case -monochromatic wave

In order to obtain further insight into the physics behind the coefficient, and to discuss the technical issue related to the appearance of resonances, we consider the simple situation when initially only a single monochromatic wave is excited in the transformed space

where B = const. ∈ C and k 0 = 0. In this case, the expression for the effective fourth-order element reduces to

At this point, the following comments are in order:

(i) The terms arising from the nonresonant interaction processes include a summation over all reciprocal lattice vectors G. Since the coefficients

k1k2k3 and frequencies ω ki are evaluated in the extended zone scheme, a summation over reciprocal lattice vectors G is equivalent to a summation over all band indices n.

(ii) For G = 0, the last term is potentially undefined since ω 0 = 0. Fortunately, this singularity is removable.

Indeed, as discussed in Sec. II, the Bloch functions are normalized as D k ∝ √ ω k . Therefore, we obtain

with some constant C(k, 0) for which 0

(iii) The denominator of the terms in the second line of Eq. (68) include the squared modulus of the initial threewave matrix element. From Eqs. (31) and (32) it is clear that the nonlinear susceptibility enters with a power of 2 into the numerators. We would like to note that the resonance conditions appearing in the denominators of N i correspond precisely to the three-wave interactions that have been eliminated from the Hamiltonian. As a result of this elimination, these processes give effective contributions to the interactions arising from nonlinear coupling via χ (3) . Considering the contribution of N

(2) k0k0k0k0 , the phase mismatch of second harmonic generation appears in the denominator, whereas for N (3) k0k0k0k0 , the sum of frequencies accounts for the interaction of two


Ranjan Kumar

Sevostian Bechta
email: bechta@kth.se

Pavel Kudinov

Florence Curnier
email: florence.curnier@cea.fr

Michel Marques
email: michel.marques@cea.fr

DYNAMIC HYBRID RELIABILITY STUDIES OF A DECAY HEAT REMOVAL SYSTEM

Keywords: Probabilistic safety assessment, dynamic reliability, decay heat removal system, stochastic hybrid automata, piecewise deterministic Markov process

Some critical safety systems exhibit the characteristics of hybrid stochastic class whose performance depends on the dynamic interactions of deterministic variables of physical phenomena and probabilistic variables of system failures. However, conventional probabilistic safety assessment (PSA) method involves static event and linked fault tree analysis and does not capture the dynamic interactions of such hybrid stochastic systems. Additionally, the existing dynamic PSA methods do not considers any repair possibility of some failed components during safety assessment. To address these issues, this paper presents a dynamic hybrid reliability assessment scheme for performance studies of repairable nuclear safety systems during a mission time. This scheme combines the features of reliability block diagram (RBD) for system compositions and partial differential equations for system physics using a customized stochastic hybrid automata tool implemented on Python platform. A case study of decay heat removal (DHR) systems has been performed using the introduced scheme. The impacts of failure rates and repair rates on sodium temperature evolution over a mission time have been analyzed. The results provide useful safety insights in mission safety tests of DHR systems. In sum, this work advances the dynamic safety assessment approach for complex system designs including nuclear power plants.

Introduction

The design of advanced nuclear safety systems, such as decay heat removal (DHR) systems are increasingly utilizing passive systems in addition to active ones to enhance the overall safety of nuclear power plants 1,[START_REF] Saez | Status of ASTRID nuclear island pre-conceptual design[END_REF] . As classified by IAEA-TECDOC-626, the category B safety systems such as reactor shutdown/emergency cooling systems, containment cooling systems, fluidic gates between process systems, etc. are well known passive systems [START_REF]Safety related terms for advanced nuclear plants[END_REF] . Although passive safety systems do not require any external input to operate, they often interact with physical processes such as thermal-hydraulic phenomena, flow of working fluids, etc. Therefore, such systems are prone to random failure due to mechanical and electrical parts and deterministic failure due to crossing the physical limits. Mathematically, the systems which involve continuous and discrete variables and subject to random process are termed hybrid stochastic systems [START_REF] Goebel | [END_REF][START_REF] Bouissou | Stochastic hybrid systems[END_REF] . The reliability evaluation of hybrid stochastic system (HSS) during a mission time is considered to be a complex task [START_REF] Labeau | [END_REF] of nuclear safety studies.

Conventionally, probabilistic safety assessment (PSA) methodology is applied for reliability studies of nuclear safety systems [START_REF]Nuclear Regulatory Commission, Reactor safety study, an assessment of accident risks in U.S. nuclear power plants[END_REF] . In practice, PSA level-1 is applied to assess the reliability of nuclear safety systems in terms of core damage frequencies (CDF) during both design and operation stages. The PSA level-1 method identifies initiating events (IE), generates accident sequences for each IE, and estimates consequence probabilities using the linked fault tree analysis of safety barriers. However, PSA level-1 studies are "static" in nature and do not take into account the dynamic states of the system. To overcome this issue, several works on "dynamic" methodologies for PSA of nuclear plants have been introduced in literature [START_REF] Aldemir | [END_REF] . These dynamic PSA methodologies can be divided into three main categories: continuous time methods, discrete time methods, graphical interface methods. The inputs of all dynamic methodologies are a time-dependent system model (such as RELAP5 9 , MELCOR 10 ), possible normal and abnormal system configurations, and transitions probabilities among these configurations.

Most of these dynamic PSA methods have been categorized under discrete dynamic event tree (DDET) as computational engines of dynamic reliability methods [6]. These are more focused on generating time-dependent event tree dynamically such as DYLAM [START_REF] Amendola | DYLAM-1, a software package for event sequence and consequence spectrum methodology[END_REF] , DETAM [START_REF] Acosta | [END_REF] , ADS 13 , MCDET 14 , ADAPT 15 , and the others. As opposed to typical dynamic reliability models, DDET methods do not take into account the real-time interaction between continuous physical processes and discrete events explicitly. Furthermore, none of the dynamic PSA methods described in Ref.8 take into account multiple initiating events, the possibility of repair of some failed parts during state transitions, functional dependencies and cascading behaviors, etc., in system modeling whereas needed to assess the safety of hybrid systems.

Nevertheless, the aforementioned methodological difficulties of hybrid systems safety assessment can be sorted out using dynamic reliability methods which provide a framework for explicitly capturing the influence of time and process dynamics on scenarios [START_REF] Labeau | [END_REF] . Dynamic reliability methods are capable of handling interactions among components and process variables explicitly. In principle, they constitute a more realistic modeling of systems for the purposes of safety analysis. Several methods such as Stochastic petri nets (SPN) [START_REF] Marsan | Modeling with generalized stochastic Petri Nets[END_REF] , Event sequence diagram (ESD) [START_REF] Swaminathan | [END_REF] , GO FLOW 18 , Dynamic Flow Method (DFM) 19 , Boolean logic Driven Markov Process(BDMP) [START_REF] Bouissou | Boolean Logic Driven Markov Processes: a powerful new formalism for specifying and solving very large Markov models[END_REF] , etc. have been proposed to assess dynamic reliability of complex systems. These methods are generally composed of input representation schemes and solution schemes and thus, do not allow direct coupling of physical behavior.

The recent works on hybrid stochastic methods demonstrate the capability to include physical equations directly in their modeling and simulations of dynamic reliability 21, 22, 23, 24, 25, and 26 . Most of these dynamic reliability approaches for hybrid stochastic systems are built upon the piecewise deterministic Markov processes (PDMP) [START_REF] Davis | [END_REF] . The PDMP allows semi-analytical resolution of the dynamic reliability. The use of hybrid stochastic automata (SHA) in PDMP provides compositional modeling capabilities to deal with large scale systems. The SHAs used in these methods are able to represent continuous and discrete variables and model the system structure and behavior. A Pythonic Object Oriented Hybrid Stochastic Automata (PyCATSHOO) [START_REF] Chraibi | Dynamic reliability modeling and assessment with PyCATSHOO: Application to a test case[END_REF] method provides the compositional modeling potentialities because it has two kinds of SHA, one discrete stochastic automaton for component behavior and the other hybrid stochastic automaton for physical phenomena. So, the compositional capabilities of PyCATSHOO make it a good candidate for modeling of stochastic hybrid behavior of complex nuclear safety systems.

Though compositional modeling of PyCATSHOO has Object Oriented, Multi-Agent Systems and functional programming paradigms, they lack inbuilt reliability block diagram (RBD) and Fault Tree (FT) features to estimate the top event probability of automaton states while selective interactions between physical and random variables.

Because in hybrid safety systems, some components interact with deterministic phenomena and are still connected to pure discrete system, i.e. support systems such as ventilation, I&C, and, electrical supply system, etc. A hybrid reliability model can solve the issue of modeling of inter-dependability of discrete systems and system dynamics 28 .

Therefore, this paper attempts to present a dynamic hybrid reliability scheme to combine traditional reliability tools such as RBD and FT with PyCATSHOO.

Rest of this paper is organized as follows. Section 2 introduces a dynamic hybrid reliability assessment scheme for large scale systems of nuclear safety systems. In this section, a detailed modeling procedure using a hybrid combination of reliability block diagram and PyCATSHOO is described. In Section 3, dynamic reliability studies of decay heat removal (DHR) systems of a sodium fast reactor type are performed. Section 4 presents results and discussion on the performed studies in Section 3. Finally, Section 5 concludes the works of this paper.

Dynamic hybrid reliability assessment approach

Dynamic reliability concept

The "dynamic reliability" studies the dynamic behavior of hybrid safety systems with random transitions and deterministic evolutions during a mission time. Mathematically, the dynamic reliability model expresses the system's dynamic behavior which is represented by a set of interacting continuous and discrete variables. The deterministic dynamics of continuous variables are formalized by a system of differential equations, and triggered behaviors of discrete variables are the transition states of the system. Probabilistic dynamics can allow interlacement between the dynamic evolution of a system and its random transitions 29 . Now, let us suppose that the states of a system with N components are defined by a vector, and the vector of the process variables at time is denoted by . If the system enters a state at time , the physical evolution of the system in that state is given by following vector equation:

,

Its formal solution is

, where is the initial value

To quantify deterministic and stochastic behavior jointly, the probability density function of the system to be in state at time with physical variable vector is introduced. Suppose, be the total transition rate out of the state and the be the specific transition rate between state to state , then the probability density obeys the Chapman-Kolmogorov equation under a Markovian hypothesis [29] and given by,

Thus, the dynamic reliability of a system can be expressed by the integral form of Eq.( 3)

as follow; 0 ( , )× [ 0 , ] × , , (5)
Where is a Dirac function representing possible trajectories of physical variables leading to at time. is the probability of a system to be in state at given the values of physical variables . is the reliability at time . The first term of Eq.( 5) corresponds to probability to be in discrete state over (0,t) without any transition and integrated over all possible deterministic trajectory given by system dynamics. The second term is made of all transitions from other states ( to state taking place before time and integrated over all possible deterministic trajectories and time instants. The transitions in second term occur according to the Markovian assumption i.e., without any memory of any prior trajectory.

Stochastic Hybrid Automata

For safety assessment of large scale systems like nuclear safety systems, analytical solution of dynamic reliability model, as explained in Eq. ( 5), is not workable [START_REF] Labeau | [END_REF] . The approximation methods such as PDMP and MC simulation with SHA are capable to solve dynamic reliability models of large systems 6, 8, 21, 22, 23, & 25 . Since Pythonic Object Oriented Hybrid Stochastic Automata (PyCATSHOO) method has two kinds of stochastic hybrid automata (SHA), this can be a competent method to deal with dynamic reliability studies of large scale system like nuclear safety systems [START_REF] Chraibi | Dynamic reliability modeling and assessment with PyCATSHOO: Application to a test case[END_REF] . A stochastic hybrid automaton can be defined as a 7-tuple: ,

Where :

 is a set of discrete states, locations or modes.

 is a continuous state space.

 is a function that gives, for each modem, , a subset of the continuous state space which represents the invariant condition of the mode .



is the set of random transitions.

 characterizes the continuous dynamics for every mode.

is an initial discrete state (initial mode) and , an initial continuous state.

Hybrid reliability assessment approach

The nuclear safety systems are usually complex and large in terms of components redundancies, interconnections, interdependency, and support systems. In such hybrid systems, some components can interact directly with the physical phenomena while others may not. Therefore, a hybrid approach for reliability assessment can be useful to study the dynamic behavior of the system. Fig. 1 exhibits the concept of hybrid reliability assessment approach in which the system consists of interacting subsystems

showing probabilistic and deterministic behavior and some probabilistic subsystems may not be repairable. Because PyCATSHOO lacks an in-built reliability analysis method, this paper introduces a scheme to create a hybrid system reliability models and solve it using a customized PyCATSHOO tool. The hybrid system reliability model formalizes the discrete variables using traditional reliability tools such as RBD and FT analysis and deterministic variables based on physical phenomena involved. Fig. 2 shows a flowchart of the iterative procedure to assess the dynamic reliability of large scale nuclear safety systems using the proposed approach. The following algorithm depicts the proposed approach:

Step-1: Hybrid system reliability is analyzed using RBD and FT analyses to find out different design options or minimal cut-sets to be studied under dynamic assessment over a mission time.

Step-2: Interacting probabilistic and deterministic variables are identified to feed into stochastic hybrid automata model.

Step-3: The stochastic hybrid automata modeling is done in two sub steps:

Step 3.1: A knowledge base is created by defining several parallel stochastic automata for implementing behavior of system components, a second automaton, termed PDMP controller, to monitor deterministic evolution, if any,

Hybrid system reliability (i) Design options (RBD &FT analysis) (ii) Governing equqations

Variables

Probabilistic: λ & μ; Deterministic: P, T, t

Start PyCATSHOO model and solution (i) Knowledge base creation (ii) System modeling

N histories

Results

Yes

No under this behavior, and message boxes for information exchange among automata and PDMP controller.

Step 3.2: A system model of large system is performed by instantiating the class created in knowledge base and connecting the automata and PDMP controller with message boxes according to architecture of system design under study.

Step-4: The developed model is run for N-histories and probabilistic and deterministic results are collected for further analysis.

A case study on decay heat removal systems

For testing the capability of introduced dynamic hybrid reliability approach to achieve the objectives, this paper applies this method to study the safety performance of DHR systems of a Sodium Fast Reactor (SFR). The SFR consists of four DHR systems namely BCS, RRA, RRB, and RRC, as shown in Fig. 3, to remove decay heat generated in the core and transferred through the hot pool [START_REF] Gauthé | Use of simplified PSA studies in support of the ASTRID design process[END_REF][START_REF] Curnier | Preliminary results from PSA studies to support the ASTRID design process[END_REF] . In normal conditions, the secondary systems (BCS) transfer the heat from the intermediate heat exchangers (IHX) to the power conversion system. Fig. 3: A schematic diagram of DHR System [31] After reactor shutdown, the residual heat can be removed through the secondary loops (BCS) bypassing turbine generator. The use of BCS for residual heat removal is limited to 3 days due to design and operational constraints. When BCS are not available, the decay heat removal relies on dedicated safety system known as RRA, RRB, and RRC (each composed of two redundant trains). Each train of RRA and RRB systems is composed of a Na/Na Decay Heat eXchanger (DHX) and an air/Na exchanger. The two trains of RRB operate in natural convection on both sodium and air sides while RRA trains operate in forced convection with electro-magnetic pumps. Another DHR system, RRC, is also investigated to complement the four DHX circuits.

System input data

This study includes three subsystems of DHR systems namely, RRA, RRB, and RRC for dynamic reliability assessment. The initial operating temperature is 500°C. DHR systems fail to function if the pool temperature exceeds 650°C temperature. Table 1 provides the DHR systems data used in this case studies. The top event failure probabilities of RRA, RRB, RRC have been calculated using RiskSpectrum® PSA tool.

Dynamic hybrid reliability analysis

Design options

As an example to test the PyCATSHOO applications, three design options were selected for this study as given in Table 2 below. The effective failure rates have been calculated using RiskSpectrum® PSA tool.

(MW, t in second)

3.2.3. Stochastic hybrid automata modeling

The knowledge base creation

There are three elements of knowledge base as follows and they are created according to the problem features:

DHR systems modeling

The system model is shown in Fig. 4. The system is modeled thanks to three instances of the class DHR and to an instance of the class POOL. Through its message box "sensor", the POOL sends current values of mass and temperature which are received by "Sensor" message box of each of the DHR system. In return, through its message box "POOL", a DHR sends the parameter of its current state and the limit of domain where such a state is valid. This message is received by "DHR" message box of the POOL.

Results

Quantitative results of dynamic reliability assessment of DHR systems using PyCATSHOO are shown in Table 3. Three combinations of those systems as worked out in Table 2 were analyzed for three sets of repair rates. Results include mean, standard deviation, and 95% confidence interval of the pool temperature for 100 histories during a mission time of 100 hours. The graphical results are shown in Fig. 5, Fig. 6, and Fig. 7 for three sets of repair rate for each design combinations. The x-axis and the y-axis show mission time and Natemperature, respectively. The lines marked in red, green, and blue show the 99%, mean, and 1% fractile of the Na-temperature for 100 histories during the mission time of 100 hrs. The 100 histories were taken because of computational limits available in this study.

This can be more with higher computational resources.

Discussion

The results obtained in Table 3 demonstrate that the proposed scheme of dynamic reliability evaluation of hybrid nuclear systems such as DHR systems provides an insight about its performance during a mission time. The hybrid combination helps to streamline the reliability composition using traditional RBD and FT analysis and dynamic transition using the customized PyCATSHOO. This approach can be applied during design and development stage to assess the performance of a design options with different mission time. This makes possible to correlate the value of probability with the deterministic behavior of the hybrid system. If the system components can be repaired within a grace period, this study enables the analysis of system behavior dynamically.

The graphical results, as shown in Fig. 4, Fig. 5, and Fig. 6, visualize the temperature evolution with respect to the failure and repair probabilities of each DHR systems combination. It can be observed that the repair possibilities, partial or complete, of the DHR systems enhance its performance to remove the decay heat successfully for a long mission time. The simulation provides the information on 95%, mean, and 1% fractile of 100 histories during the mission time which provides the uncertainty estimation for the DHR system behavior.

From the results, it can be found that the best combination among the three design options is 2 RRA, 2RRB, 1RRC in parallel for successful removal of decay heat even without any repair possibility. All the combinations without some repair possibility are not preferable for long mission time. Other combinations and dependency of support systems needs to be investigated for the safest possible systems design.

Though the studies on dynamic reliability evaluation provide an efficient tool to model and simulate large hybrid safety system, the probability density functions of component failure and repair have been assumed exponential for simplification. The customized PyCATSHOO tool allows including user's defined failure and repair law in modeling and simulation. Moreover, the future study on dynamic reliability of DHR systems should include the detailed modeling of DHR components and their dependencies on support systems. This may be done using dependency modeling of reliability theory.

Conclusion

To capture the dynamic interactions of probabilistic and deterministic variables during a mission time, this paper presents a hybrid scheme of dynamic reliability assessment scheme. This hybrid approach includes a combination of conventional reliability compositional analysis tools and the customized Pythonic objects oriented stochastic hybrid automata method. A case study on dynamic reliability studies of decay heat removal systems of a sodium fast reactor is performed using the introduced approach.

The results display that repair possibility of some components enhances the overall performance of decay heat removal systems during mission times. The dynamic reliability study results are more explicit than those from the static probabilistic safety assessment methods in mission safety test of decay heat removal systems. In future works, user's defined failure and repair laws with dependencies modeling can be included in similar studies of large scale safety systems.

Fig. 1 .

 1 Fig.1. Hybrid reliability concept for probabilistic and deterministic sub-systems

Fig. 2 .

 2 Fig.2. A scheme of dynamic hybrid reliability assessment of large scale systems

  (a) DHR class: Each DHR system is defined by automaton which includes declaration of the state with instruction and creation of the transitions. There are two states/modes of DHR class in this case. The time of occurrence of such failure is random variable which follows exponential distribution. (b) Message box: Message boxes provide the ability to the instances of such a class to interact with the other system objects. Message box declaration for DHR class got two message boxes. In this message box definition, a parameter is invoked when a new value for temperature is received and the results are reset for evaluation of DHR automata transition conditions. This box is used to send information to the POOL class when DHR state is changed. (c) POOL class: This class models the behavior of the POOL component. Because this class includes physical phenomena with continuous variables, Pool class is modeled by PDMP controller. This PDMP controller will manage a one dimensional state variable vector.

Fig. 4 :

 4 Fig.4: The PyCATSHOO model of the DHR systems

  (a) Repair rate=0.99 (b) Repair rate=0.

Fig. 5 :Fig. 6 :Fig. 7 :

 567 Fig.5: Pool temperature evolution with 2 RRA, 2RRB, 1RRC (in Series)

Table 1 :

 1 DHR systems data

		Top event failure probability	Operating time (hr)	Initial state
	RRA	3,63-02	100	ok
	RRB	3,12E-02	100	ok
	RRC	6,35E-03	100	ok

Table 2 :

 2 RBD analysis and effective failure rate

	S.	Combinations	RBD-Physical phenomena		Effective
	N.						failure rate
							(100 hrs)
	1	2 RRA, 2RRB,
		1RRC			RRA	RRB
		(in Series)			RRC	Na-pool	9,15E-5
					RRA	RRB
	2	2 RRA, 2RRB,
		1RRC			RRA
		(all	are	in
		parallel)		RRA
					RRB	Na-pool	2,00 E-10
					RRB
					RRA
	3	2 RRA, 2RRB,
		1RRC operating	are	RRA	Na-pool	RRB	RRA=1,84 E-5
		independently	RRA		RRB	RRB=7,683 E-6
							RRC=6,35
						RRC	E-5

Table 3 :

 3 Results with 100 histories

Acknowledgments

The authors thank Valentin Rychkov of EDF Clamart in Paris for his technical guidance on PyCATSHOO installation and usage for case studies conducted in this paper.