Davide Spirito

Dominique Coquillat

Sergio L De Bonis

Antonio Lombardo

Matteo Bruna

Andrea C Ferrari

Vittorio Pellegrini

Alessandro Tredicucci

Wojciech Knap

Miriam S Vitiello
email: miriam.vitiello@sns.it

High performance bilayer-graphene Terahertz detectors

We report bilayer-graphene field effect transistors operating as THz broadband photodetectors based on plasma-waves excitation. By employing wide-gate geometries or buried gate configurations, we achieve a responsivity~1.2V/W (1.3 mA/W) and a noise equivalent power~2×10 -9 W/√Hz in the 0.29-0.38 THz range, in photovoltage and photocurrent mode. The potential of this technology for scalability to higher frequencies and the development of flexible devices makes our approach competitive for a future generation of THz detection systems.

Generation and detection of radiation across the far infrared or Terahertz (THz) region of the electromagnetic spectrum is promising for a large variety of strategic applications, 1 ranging from biomedical diagnostics 2 to process and quality control, 3 homeland security 1 and environmental monitoring. 4 Due to its non-ionizing nature, 1 THz radiation can penetrate many commonly used dielectrics, 1 otherwise opaque for visible and mid-infrared light, allowing detection of specific spectroscopic features 5 with a sub-millimeter diffractionlimited resolution. 1 Quantum cascade (QC) lasers 6 are nowadays capable to reach more than 120mW output powers in continuous wave, and to operate in a broad frequency range (1.2 -4.7 THz). 7 However, a low-cost roomtemperature (RT) THz detector technology with hundreds V/W responsivity, and fast (GHz modulation frequency) is still missing. Conventional THz and sub-THz detection systems based on incoherent (pyroelectric, Golay cell, Si bolometers) [START_REF] Sizov | THz detectors[END_REF] or coherent (heterodyne mixers) [START_REF] Sizov | THz detectors[END_REF] approaches are either very slow ( 100 Hz modulation frequency) [START_REF] Sizov | THz detectors[END_REF] or require deep cryogenic cooling, [START_REF] Sizov | THz detectors[END_REF] or, if exploiting fast (hundred MHz) nonlinear electronics, show a significant drop of performance above 1 THz. [START_REF] Sizov | THz detectors[END_REF] Recently, [9][10][11][12][13] detectors based on the gate-modulation of the channel conductivity by the incoming radiation have been realized, employing III-V materials in high-electron-mobility transistors (HEMT), 9,10 two dimensional electron gas field effect transistors (FET), 9 Si-metal-oxide-semiconductor FETs and complementary metal-oxide semiconductor architectures, 11 showing fast response times and high responsivities. This approach was also extended 12 to InAs nanowire one dimensional FETs, allowing RT detection of THz emission in the 0.3 -3 THz range. 12,13 The operating mechanism of a FET detector 14 can intuitively be interpreted as deriving from the nonlinearity arising from the simultaneous modulation of carrier density and drift velocity by the THz radiation. 14 When a THz beam is funneled onto the channel of a FET, it modulates the gate and source/drain potentials, generating a continuous source-drain voltage u (or current i) 14 which can be maximized by varying the gate bias V G . This modulation leads to excitation of plasma oscillations that can either propagate in the channel (plasma waves) 14 or be damped (overdamped plasma waves). 14 At THz frequencies plasma waves decay on a distance d w = s, [14] where s is the plasma wave velocity (typically assumed ~10 8 cm/s independently from the channel material) [14] and  the momentum relaxation time. If the FET gate length w G > d w , plasma waves excited at the source decay before reaching the drain, and broadband THz detection is achieved. 9 The same condition occurs in the low frequency regime 2 <1 ( being the frequency of the THz radiation), where the plasma oscillations decay with a frequency-dependent characteristic length 9,14 l ( ) .

The unique optoelectronic properties of graphene make it an ideal platform for a variety of photonic applications. [START_REF] Tredicucci | Device Concepts for Graphene-Based Terahertz Photonics[END_REF][START_REF] Bonaccorso | Graphene photonics and optoelectronics[END_REF] In particular, due to its high carrier mobility, gapless spectrum, and frequency-independent absorption, [START_REF] Bonaccorso | Graphene photonics and optoelectronics[END_REF][17][18] graphene is very promising material for the development of THz detectors, 19 still severely lacking in terms of solid-state devices.

Refs. 20, 21 theoretically suggested that THz detection in graphene could be mediated by electronhole pairs or interband photo-generation. However, we experimentally reported RT broadband THz detection as a result of overdamped plasma-waves excitation in the channel of a graphene-FET 19 operating in photovoltage mode, with an electrical output per optical input (responsivity) R v ~ 0.015 V/W and noise equivalent powers NEP = 30 nWHz 1/2 . Subsequently, Ref. 22 reported ultrafast (10 ps response time) THz detectors based on photo-induced bolometric effects in graphene, operating in photocurrent mode, with a responsivity ~ 5 nA/W. Ref. 23 then combined plasmonic/bolometric graphene FETs, reaching R v ~ 150 µV/W and operating in photovoltage mode at much higher frequencies (2.5 THz).

Here we report THz detectors operating at RT in either photovoltage or photocurrent mode with a significant enhancement of sensitivity and lowered NEPs compared to the state of the art. This is achieved by using either large (1 µm) gate lengths, or buried gate geometries on a bilayer graphene (BLG) FET. We use BLG instead of single layer graphene (SLG) since the modulation of carrier density was proved to be more effective in the former case, 19 allowing a higher responsivity at THz frequencies.

The devices are prepared as follows. Flakes are mechanically exfoliated from graphite 24 on an intrinsic Si substrate covered with 300nm SiO 2 . BLGs are selected and identified by a combination of optical microscopy 25 and Raman spectroscopy. 26,27 These are then used to fabricate two sets of FETs (samples A, B). In A, source (S) and drain (D) contacts are patterned by electron beam lithography (EBL) and metal evaporation (5nm Cr /80 nm Au). The S contact is connected to one lobe of a 50° bow-tie antenna, and D to a metal pad. After placement of 35nm HfO 2 by atomic layer deposition (ALD), the other lobe of the antenna is fabricated, and constitutes the FET gate (g). The channel length is L SD = 2.5 m, while the gate length w G =1 m (Fig. 1a). In B, one lobe of a log-periodic circular-toothed antenna is fabricated by EBL and metal evaporation to act as buried gate. After deposition of 35 nm HfO 2 by ALD, S and D electrodes are fabricated. Similar to sample A, S is the second lobe of the antenna, while the drain is a metal line connecting to a bonding pad. The channel length is 2.5 m, its width 7.5 m, while w G = 0.3m (Fig. 1b). A BLG flake is then placed onto the pre-fabricated electrodes by wet transfer. 28 Fig. 2 compares the Raman spectra of the flake prior and after placement onto the electrodes. The "2D" peak shows the characteristic multi-band shape of BLGs. 26,27 The absence of a significant Raman "D" peak prior and after transfer indicates the low-defect samples are not damaged by this process. Scanning electron microscope images of both devices are reported in Fig. 1c,d According to a diffusive model of transport, 14,30,31 for long gates (w G >> d w ) or in a low-frequency regime (2 <1), a second-order nonlinear response is expected in a FET when an oscillating THz is applied between G and S. In this case the photovoltage signal can be qualitatively extracted from the currentvoltage characteristics via the equation [14] , Du variation from negative to positive values around the CNP for A, consistent with ambipolar transport, 28 while only a small sign switch is recorder for B, in agreement with the evident asymmetry of the transfer characteristics. 18 On the right vertical axis of Figs 3c-d In all cases, the vertically polarized incoming radiation impinges on the FET mounted in a dual-in-line package, with an optical power P t = 1.8 mW (fixed frequency source) or P t = 0.34 mW (tunable source).

Figure 4a plots R ν extracted from the photovoltage measurements on A, while focusing the 0.292

THz radiation on it and keeping V SD =0. R ν is related to the measured photo-induced voltage √ via : R ν , where S t is the radiation beam spot area, S a is the detector active area and LIA the lock-in signal. 19 Loading effects due to the finite impedance of the preamplifier input are neglected. This procedure assumes that the entire power incident on the antenna is fully coupled to the FET, meaning that the present R ν is a lower limit. Taking into account our beam spot diameter d ≈ 4 mm, we have S t = πd 2 /4 = 12.6 × 10 -6 m 2 . Since the total detector area is smaller than the diffraction limited one S λ = λ 2 /4, the active area is taken equal to S λ . The dependence of on V G is in qualitative agreement with the trend shown in Fig. 3c, thereby proving that the detector still operates in the broadband overdamped regime. Note that the sign of the photovoltage changes abruptly in the vicinity of the CNP, following the switch of sign of the derivative in -. For a qualitative evaluation of the coupling efficiency of the THz radiation by the antenna, we can use the outcome of a fitting procedure to the responsivity of other FET detectors with similar antenna designs measured in the same experimental set-up, 31 giving = 1.25×10 -4 V 2 . By using the extrapolated K value, the predicted photovoltage (Fig. 3c) is in excellent agreement with Fig. 3a. The comparison with Δu * shows a good agreement with the measured R ν , with a small voltage offset possibly due to loading impedance and/or heating effects. We get Δu/Δu* ~ 0.25, as expected if one considers that ΔI×σ -1 is much less influenced by the loading than the photovoltage, since the former is a purely dc measurement, not affected by the device capacitive reactance.

Maximum responsivities of 1.2 V/W, i.e. more than one order of magnitude larger than in previous reports 10,19,23 are reached. This is interpreted as due to the larger gate length, as predicted by the photoinduced voltage equation: 14 * ( )+ where x is the gated region length and l the previously defined characteristic decay length of the ac voltage (for sample A is l ~ 1.2 µm). 19 The ratio w G /l ~ 0.8 is then more than one order of magnitude larger than in previous BLG-FET THz detectors (w G /l= 0.03), 19 meaning that any competing thermoelectric effects possibly arising from the ungated regions 19 are here less relevant. Also, no evidence of THz driven interband transitions 19 is seen.

Figure 4b plots R ν extracted from the photovoltage measurements on B at both 0.292 and 0.358 THz.

The R ν vs V G trend is almost independent from the radiation frequency (confirming the broadband detection mechanism) with maximum R ν = 1.12 V/W. For comparison, Fig. 4c plots the responsivity values obtained while B is operating in photocurrent mode at 0.358 THz. The maximum R ν = 1.3 mA/W is in agreement with that recorded while the device is operating in the photovoltaic mode, once the 800Ω channel resistance is taken into account. The w G /l ratio is here ≈ 0.8, i.e. comparable to A, motivating the equivalent maximum responsivity values. Also in agreement with A, we have Δu* > Δu (see Fig. 3d), being Δu/Δu* ~ 0.5.

Figure 5a plots the photovoltage value measured at V G = 2.5V, while varying the frequency of the incoming THz beam. The frequency response has several peaks, as expected from the broadband nature of the antenna, with the broader and more intense band at 0.358 THz.

Figure 5b shows the Johnson-Nyquist noise extracted from the data of Fig. 3b, when B is operating

in photovoltage mode N V = √
, or photocurrent mode N I =√ . Assuming that the Johnson-Nyquist noise is the dominant contribution, 32 as typical for THz FET detectors, 12,13,19 minimum NEP ≈ 4×10 - 9 W/Hz 1/2 and NEP ≈ 2×10 -9 W/Hz 1/2 are reached in photovoltage and photocurrent mode, respectively, comparable with most commercially available RT THz detectors. [START_REF] Sizov | THz detectors[END_REF] In summary, by employing wide-gate geometries or buried gate configurations, we fabricated bilayer-graphene FETs with R v ≈ 1.2V/W (1.3 mA/W) and NEP ≈ 2×10 -9 W/√Hz. These value are competitive with the performance of commercially available detection systems, paving the way to a realistic exploitation for large-area, fast imaging of macroscopic samples, spectroscopic applications, process control and biomedical diagnostics. Raman Shift (cm -1 )

Figure captions

Intensity (arb-units) Δu (V) Δi (nA)

  .The devices are electrically characterized by means of two voltage generators to drive independently V G and the source-to-drain (V SD ) voltage. The drain contact is connected to a current amplifier converting the current into a voltage signal with an amplification factor of 10 4 V/A. The latter signal is then measured with an Agilent 34401A voltmeter reader. The current (I SD ) as a function of V G characteristics are measured while sweeping V G , keeping V SD = 0.001V. Two sets of experiments are performed in order to compare the FET transport behavior: one under 0.292 THz irradiation, the other without. Figs3a-bshow that the conductivity, goes through a minimum when the chemical potential below the gate crosses the charge neutrality point (CNP), as in Ref.19.Moreover, a shift of the CNP of ~ 0.1V is detected when the THz beam is focused on the sample. The channel resistance (R) varies from 2.8 to 4.2kΩ in A, and from 500 to 800Ω in B, corresponding to a mobility µ ≈ 3000 cm 2 /Vs and a carrier density n o ≈ 2×10 11 cm -2 , and µ ≈ 200 cm 2 /Vs, n o ≈ 70 ×10 11 cm -2 , respectively, as extracted from the transfer characteristics, following the procedure described in Ref.29. This implies a plasma wave decay distance d w ≈ 1.5µm ≥ w G (A) and l ≈ 350 nm ≥ w G (B), meaning that both are overdamped.

-

  where the constant K represents the coupling efficiency of the incoming radiation to the antenna together with a loading factor determined by the impedance of the transistor. The - trend shown on the left vertical axis of Figs. 2c,d indicates a clear

  photocurrent signal is measured with a transimpedance amplifier, having amplification factor G n = 10 7 V/A.

Figure 1 :

 1 Figure 1: (a-b) Schematic BLG-FETs with (a) top gate, (b) buried gate. (c-d) False-color SEM micrographs of the devices schematically shown in (c) panel a, and (d), panel b.

Figure 2 :

 2 Figure 2: Measured Raman spectra of the BLG flake prior (bottom curves) and after (upper curves) placement onto the electrodes.

Figure 3 :

 3 Figure 3: (a-b) Source-drain current (I SD ) measured as a function of V G in (a) sample A and (b) B, while keeping V SD = 1mV and with (red curve) or without (black curve) a 0.292 THz illumination. (c-d) Left vertical axis: Derivative of the conductivity multiplied by the resistance, as a function of V G , for (c) sample A, and (d) sample B. Right vertical axis: difference between the current measured while focusing the THz beam on the sample I SD,THz and without it I SD,no THz , multiplied by the channel resistance in case of (c) sample A, or (d) sample B. The dashed vertical lines indicate the charge neutrality point.

Figure 4 :

 4 Figure 4: (a) Room-temperature (RT) responsivity as a function of V G measured while sample A is operating in photovoltage mode and while keeping the polarization of the incoming 0.292 THz beam parallel to the antenna axis. The related measured photovoltage signal (blue curve) is shown on the right vertical axis. (b) RT responsivity as a function of V G measured while sample B is operating in photovoltage mode and while keeping the polarization of the incoming 0.292 THz (upper red curve) or 0.358 THz (pink curve) beams parallel to the antenna axis. The related measured photovoltage signals (blue and black curves) are shown on the right vertical axis. c) RT responsivity as a function of V G measured while sample B is operating in photocurrent mode (c) at 0.358 THz. The dashed line marks the zero value of output voltage.

Figure 5 :

 5 Figure 5: (a) Photocurrent (Δi) and photovoltage (Δu) signals measured in sample B as a function of frequency of the incoming THz beam while V G = 2.4 V and V SD = 1 mV. (b) Johnson-Nyquist noise extracted from the data of Fig. 2b for sample B while it is operating in photocurrent (N I ) or photovoltage mode (N V ).

  Figure 2

Figure 5

 5 Figure 5

Acknowledgments

We acknowledge funding from the Italian Ministry of Education, University, and Research (MIUR) through the program "FIRB -Futuro in Ricerca 2010" RBFR10LULP "Fundamental research on Terahertz photonic devices" and RBFR10M5BT PLASMOGRAPH by EU Graphene Flagship (contract no. CNECT-ICT-604391), ERC projects NANOPOTS, Hetero2D, EU projects GENIUS, CARERAMM, MEM4WIN, a Royal Society Wolfson Research Merit Award, EPSRC grants EP/K01711X/1, EP/K017144/1, the Italian Ministry of Economic Development through the ICE-CRUI project "Teragraph", the CNRS through the GDR-I project "Semiconductor sources and detectors of THz frequencies", the Region Languedoc-Roussillon Terahertz Platform.


Fernand G J Gounaris

Michel Renard

published or not. The documents may come     The Higgs boson structure functions as signposts of new physics


Jawher Jerray

Laurent Fribourg

Determination of limit cycles using stroboscopic set-valued maps

Keywords: parametric differential equations, periodicity, uncertainty

Given a dynamical system Σ p with a parameter p taking its values in a fixed interval Q, we present a simple criterion of set inclusion which guarantees that the Euler approximate solutions of Σ p0 for some value p 0 ∈ Q converge to a limit cycle E. Moreover, we characterize a compact set I containing E which is invariant for the exact solutions of Σ p whatever the value of p ∈ Q. We illustrate the application of our method on the example of a parametric Van der Pol system driven by a periodic input.

INTRODUCTION

Given a differential system Σ : dx/dt = f (x) of dimension d, an initial point x 0 ∈ R d , a real ε > 0, and a ball B 0 = B(x 0 , ε) 1 , we present here a simple method allowing to find a bounded invariant set of Σ containing the trajectories starting at B 0 . This invariant set has the form of a tube whose center at time t is the Euler approximate solution x(t) of the system starting at x 0 , and radius is a function δ ε (t) bounding the distance between x(t) and an exact solution x(t) starting at B 0 . The tube can thus be described as t≥0 B(t) where B(t) ≡ B(x(t), δ ε (t)).

To find a bounded invariant, we then look for a positive real T such that B((i + 1)T ) ⊆ B(iT ) for some i ∈ N. In case of success, the ball B(iT ) is guaranteed to contain the "stroboscopic" sequence {B(jT )} j=i,i+1,... of sets B(t) at time t = iT, (i + 1)T, . . . . It follows that the bounded portion t∈[iT,(i+1)T ] B(t) is equal to t∈[iT,∞) B(t), and thus constitutes the sought bounded invariant set.

We apply the finding of such a (forward) invariant set to the analysis of parameterized time-periodic differential systems. We illustrate our results on the example of a parametric Van der Pol system driven by a common periodic input.

1.0.0.1. Comparison with related work We explain here some similarities and differences of our method with several kinds of related work.

• There exists a trend of work on the generation of torus-shaped invariants using stroboscopic maps of quasi-periodic systems with possible parameters (see, e.g., [START_REF] Baresi | Quasi-periodic invariant tori of time-periodic dynamical systems: Applications 5 to small body exploration[END_REF]; [START_REF] Gómez | The dynamics around the collinear equilibrium points of the RTBP[END_REF]; [START_REF] Olikara | Numerical method for computing quasi-periodic orbits and their stability in the restricted three-body problem[END_REF]). A first difference is that these works consider stroboscopic

1 B(x 0 , ε) is the set {z ∈ R d | z -x 0 ≤ ε} where • denotes the Euclidean distance.
mappings of points of R d , while our stroboscopic maps apply to sets of points. A second difference is that they often use a Fourier analysis in the frequency domain (using, e.g., the notion of "radii polynomials" Castelli and Lessard ( 2013)) while we remain in the time domain. • Our method makes use of a rigorous time-integration method in order to enclose the exact solutions with tubes, which is similar to what is done using high order of Taylor method in ODE integration as proposed by Lohner or Taylor models [START_REF] Lohner | Enclosing the solutions of ordinary initial and boundary value problems[END_REF]; [START_REF] Zgliczynski | C1 Lohner algorithm[END_REF]. Such methods are used in [START_REF] Castelli | Rigorous numerics in Floquet theory: computing stable and unstable bundles of periodic orbits[END_REF]; [START_REF] Kapela | Computer assisted proofs for nonsymmetric planar choreographies and for stability of the eight[END_REF] to rigorously compute the eigenvalues of a so-called "monodromy matrix", which allows to determine the linear stability of the equilibrium points of the system. Unlike these works, our stability analysis does not try to compute such eigenvalues of monodromy matrices.

• Our method shares also some common features with the works of [START_REF] Aylward | Stability and robustness analysis of nonlinear systems via contraction metrics and SOS programming[END_REF]; [START_REF] Aminzare | Contraction methods for nonlinear systems: A brief introduction and some open problems[END_REF]; van den Berg and Queirolo (2020), which aim at proving a contractivity property of the system (i.e., that any two solutions converge exponentially to each other). In [START_REF] Aminzare | Contraction methods for nonlinear systems: A brief introduction and some open problems[END_REF], contractivity amounts to the finding of a negative definitive quadratic form (which is equivalent to the existence of a Lyapunov function for the system). In [START_REF] Aylward | Stability and robustness analysis of nonlinear systems via contraction metrics and SOS programming[END_REF], "Squares-of-Sum programming is used to find ranges of uncertainty under which a system with uncertain perturbations is always contracting with the original contraction metric". In van den Berg and Queirolo (2020), they turn the stability problem into the contractivity of a fixed point operator that is checked with the assistance of a computer. The difference here is that we do not try to prove a contractivity property, but only the existence of two set-valued snapshots, one of which is included in the other one. This is a much weaker property and easier to prove.

Our method is simple, efficient, and allows us to construct a cyclic approximate solution E (using Euler's method) together with a (forward) invariant set I around E. Note however that, even if the approximate and exact solutions are very close in practice, our method does not guarantee the formal existence of cyclic exact solutions inside I, which requires more technical methods such as Poincaré maps.

1.0.0.2. Plan of the paper In Section 2, we present our method, then explain how to apply it to the analysis of parameterized systems in Section 3. We conclude in Section 4.

METHOD

Euler's method and error bounds

Let us consider the differential system:

dx(t) dt = f (x(t)),
with states x(t) ∈ R d . We will use x(t; x 0 ) (or sometimes just x(t)) to denote the exact continuous solution of the system at time t, for a given initial condition x 0 . We use x(t; y 0 ) (or just x(t)) to denote Euler's approximate value of x(t; y 0 ) (defined by x(t;

y 0 ) = y 0 + tf (y 0 ) for t ∈ [0, τ ],
where τ is the integration time-step).

We suppose that we know a bounded region S ⊂ R d containing the solutions of the system for a set of initial conditions B 0 and a certain amount of time. We now give an upper bound to the error between the exact solution of the ODE and its Euler approximation on S (see Le [START_REF] Le Coënt | Guaranteed optimal reachability control of reaction-diffusion equations using one-sided Lipschitz constants and model reduction[END_REF]; Le [START_REF] Coënt | Control synthesis of nonlinear sampled switched systems using Euler's method[END_REF]). Definition 1. Let ε be a given positive constant. Let us define, for t ∈ [0, τ ], δ ε (t) as follows: if λ < 0 :

δ ε (t) = ε 2 e λt + C 2 λ 2 t 2 + 2t λ + 2 λ 2 1 -e λt 1 2 if λ = 0 : δ ε (t) = ε 2 e t + C 2 (-t 2 -2t + 2(e t -1)) 1 2 if λ > 0 : δ ε (t) = ε 2 e 3λt + C 2 3λ 2 -t 2 - 2t 3λ + 2 9λ 2 e 3λt -1 1 2
where C and λ are real constants specific to function f , defined as follows:

C = sup y∈S L f (y) ,
where L denotes the Lipschitz constant for f , and λ is the "one-sided Lipschitz constant" (or "logarithmic Lipschitz constant" Aminzare and Sontag ( 2014)) associated to f , i. e., the minimal constant such that, for all y 1 , y 2 ∈ S:

f (y 1 ) -f (y 2 ), y 1 -y 2 ≤ λ y 1 -y 2 2 , ( H0
) where •, • denotes the scalar product of two vectors of S.

The constant λ can be computed using a nonlinear optimization solver (e. g., [START_REF] Cplex | V12. 1: User's manual for cplex[END_REF]) or using the Jacobian matrix of f (see, e. g., [START_REF] Aminzare | Contraction methods for nonlinear systems: A brief introduction and some open problems[END_REF]). Proposition 1. Le [START_REF] Coënt | Control synthesis of nonlinear sampled switched systems using Euler's method[END_REF] Consider the solution x(t; y 0 ) of dx dt = f (x) with initial condition y 0 and the approximate Euler solution x(t; x 0 ) with initial condition x 0 . For all y 0 ∈ B(x 0 , ε), we have:

x(t; y 0 ) -x(t; x 0 ) ≤ δ ε (t).
Proposition 1 underlies the principle of our set-based method where set of points are represented as balls centered around the Euler approximate values of the solutions. This illustrated in Fig. 1: for any initial condition x 0 belonging to the ball B(x 0 , δ(0)) with δ(0) = ε, the exact solution x 1 ≡ x(τ ; x 0 ) belongs to the ball B(x 1 , δ ε (τ )) where x1 denotes the Euler approximation x0 + τ f (x 0 ) at t = τ .

Figure 1. Illustration of Proposition 1

Systems with bounded uncertainty

Let us now show how the method extends to systems with "disturbance" or "bounded uncertainty". A differential system with bounded uncertainty is of the form dx(t) dt = f (x(t), w(t)),

with t ∈ R d ≥0 , states x(t) ∈ R d , and uncertainty w(t) ∈ W ⊂ R d ,
where W is a compact (i. e., closed and bounded) set. We assume that any possible disturbance trajectory is bounded at any point in time in the compact set W . We denote this by w(•) ∈ W, which is a shorthand for w(t) ∈ W, ∀t ≥ 0. The diameter of W (i.e., the maximal distance between two elements of W) is denoted by |W|. See Schürmann and Althoff (2017b,a) for details. We now suppose (see Le Coënt et al. (2017a)) that there exist constants λ ∈ R and γ ∈ R ≥0 such that, for all y 1 , y 2 ∈ S and w 1 , w 2 ∈ W:

f (y 1 , w 1 ) -f (y 2 , w 2 ), y 1 -y 2 ≤ λ y 1 -y 2 2 + γ y 1 -y 2 w 1 -w 2 (H1).
This formula can be seen as a generalization of (H0) (see Section 2.1). Recall that λ has to be computed in the absence of uncertainty (|W| = 0). The additional constant γ is used for taking into account the uncertainty w. Given λ, the constant γ can be computed itself using a nonlinear optimization solver (e. g., [START_REF] Cplex | V12. 1: User's manual for cplex[END_REF] Consider a point x 0 ∈ S and a point y 0 ∈ B(x 0 , ε).Let x(t; y 0 ) be the exact solution of the system dx(t) dt = f (x(t), w(t)) with bounded uncertainty W and initial condition y 0 , and x(t; x 0 ) the Euler approximate solution of the system dx(t) dt = f (x(t), 0) without uncertainty (|W| = 0) with initial condition x 0 . We have, for all w(•) ∈ W and t ∈ [0, τ ]:

x(t; y 0 ) -x(t; x 0 ) ≤ δ ε,W (t). with

• if λ < 0, δ ε,W (t) = C 2 -λ 4 -λ 2 t 2 -2λt + 2e λt -2 + 1 λ 2 Cγ|W| -λ -λt + e λt -1 + λ γ 2 (|W|/2) 2 -λ (e λt -1) + λε 2 e λt 1/2
(1)

• if λ > 0, δ ε,W (t) = 1 (3λ) 3/2 C 2 λ -9λ 2 t 2 -6λt + 2e 3λt -2 + 3λ Cγ|W| λ -3λt + e 3λt -1 + 3λ γ 2 (|W|/2) 2 λ (e 3λt -1) + 3λε 2 e 3λt 1/2
(2)

• if λ = 0, δ ε,W (t) = C 2 -t 2 -2t + 2e t -2 + Cγ|W| -t + e t -1 + γ 2 (|W|/2) 2 (e t -1) + ε 2 e t 1/2 (3)
We will sometimes write δ W (t) and δ(t) instead of δ ε,W (t) and δ ε (t) respectively.

Correctness

Consider a differential system Σ : dx/dt = f (x, w) with w ∈ W, an initial point x 0 ∈ R d , a real ε > 0 and a ball B 0 = B(x 0 , ε). Let B W (t) denote B(x(t), δ ε,W (t)) where x(t) is the Euler approximate solution of the system without uncertainty and initial condition x 02 . It follows from Proposition 2 that t≥0 B W (t) is an invariant set containing B 0 . We can make an a stroboscopic map of this invariant. by considering periodically the set B W (t) at the moments t = 0, T, 2T , etc., with T = kτ for some k (τ is the time-step used un Euler's method).

If moreover, we can find an integer i ≥ 0 such that B W ((i + 1)T ) ⊆ B W (iT ), then we have B W (iT ) = j=i,i+1,... B W (jT ) and t∈[iT,(i+1)T ] B W (t) = t∈[iT,∞) B W (t). The set t∈[iT,(i+1)T ] B W (t) is thus a bounded invariant set which contains all the solutions x(t) starting at B 0 , for t ∈ [iT, ∞). In the phase space, this invariant set has a "torus" shape. We have: Proposition 3. Suppose that there exist T = kτ with k ∈ N, and i ∈ N such that:

B W ((i + 1)T ) ⊆ B W (iT ).
Then I W ≡ t∈[0,T ] B W (iT +t) is a compact (i.e., bounded and closed) invariant set containing, for t ∈ [iT, ∞), all the solutions x(t) of Σ with initial condition in B 0 .

APPLICATION TO PARAMETRIC SYSTEMS

Let us now consider a family {Σ p } p∈P of differential systems Σ p of the form dx/dt = f p (x) involving a parameter p ∈ P (but no uncertainty). It is useful to find a subset Q of P and a system Σ : dx/dt = f (x, w) with uncertainty such that, for any p ∈ Q, Σ p is a particular form of Σ for an appropriate uncertainty function w(•) ∈ W. This is useful to infer certain common properties of the solutions of {Σ p } p∈Q from the analysis of the system Σ with uncertainty (cf. [START_REF] Aylward | Stability and robustness analysis of nonlinear systems via contraction metrics and SOS programming[END_REF], Section 5). Example 1. Consider the Van der Pol system Σ p with parameter p ∈ R driven by a periodic input of the form 0.015cos(t) (see, e.g, Capinski et al. (arXiv:1905.08116, 2020)), with initial condition in B 0 = B(x 0 , ε) for some given x 0 ∈ R 2 and ε > 0 (see van den Berg and Queirolo (2020)):

du 1 /dt = u 2 du 2 /dt = p(1 -u 2
1 )u 2 -u 1 + 0.015cos(t). Consider now the system Σ with uncertainty w(•) ∈ W = [-0.1, 0.1] (hence |W| = 0.2) and initial condition x 0 :

du 1 /dt = u 2 du 2 /dt = (p 0 + w(t))(1 -u 2
1 )u 2 -u 1 + 0.015cos(t), with p 0 = 0.4. It is easy to see that each solution of Σ p with p ∈ Q = [p 0 -0.1, p 0 + 0.1] = [0.3, 0.5] is a particular solution of Σ . Theorem 1. (cyclic approximate solution) Given a system Σ : dx/dt = f (x, w) with uncertainty w ∈ W and a set of initial conditions B 0 ≡ B(x 0 , ε), suppose that any solution x(t) of the parametric system Σ p : dx/dt = f p (x) with p ∈ Q = [p 0 -|W|/2, p 0 +|W|/2] and initial condition x 0 is a solution of Σ with initial set of conditions B 0 for some uncertainty function w(•) ∈ W. Let us suppose besides that, for Σ , there exist i ∈ N and T = kτ with k ∈ N such that:

B W ((i + 1)T ) ⊆ B W (iT ).
Then:

(1) When n → ∞, the sequence

{c n ≡ x((i + n)T ; x 0 )} converges to a limit c * ∈ R d satisfying x(T ; c * ) = c * . (2)
The curve E ≡ t≥0 x(t; c * ) is cyclic with period T , i.e.: x(t + T ; c * ) = x(t; c * ) for all t ≥ 0.

(3) When t → ∞, x(t; c 0 ) converges to E. More precisely, for all s ≥ 0, there exists c * s ∈ E : x(s + nT ; c 0 )c * s → 0 as n → ∞. This means that E is an attractive limit cycle for the approximate solution x(t; c 0 ) (whence x(t; x 0 )) as t → ∞. (4) The set I 0 = t≥0 B(x(t; c 0 ), δ(iT + t)) includes E and is (forward) invariant for Σ p0 , i.e.: all solution of Σ p0 starting from I 0 remains in I 0 . (5) For each p ∈ Q, the set

I W = t∈[0,T ] B W (iT + t)
is (forward) invariant for Σ p , i.e.: all solution of Σ p starting from I W remains in I W .

Proof. (sketch) The arguments of the 5 items of the theorem are as follows:

(1) As the balls B W (i + nT ) are included in each other, this implies that their radii ε n are decreasing and converge to a limit ε * . This also implies that the centers c n of these balls converge to a point c * ∈ R d , and that balls B W (i + nT ) converge to the ball B(c * , ε * ). This limit ball is a fixed point of the function mapping a ball B(c, ρ) to B(x(T ; c), δ ρ,W (T ), and satisfies:

B(x(T ; c * ), δ ε * ,W (T ) = B(c * , ε * ).
In particular, we have:

x(T ; c * ) = c * . (2) Follows from x(T ; c * ) = c * (see above) and x(t + T ; c * ) = x(t; x(T, c * )) = x(t; c * ). (3) For all s ∈ [0, T ), x(s + nT ; c 0 ) = x(s; x(nT, c 0 )) = x(s; c n ). It follows, that for all s ∈ [0, T ), x(s + nT ; c 0 ) -c * s = x(s; c n ) -c * s → 0 as n → ∞ with c *
s ≡ x(s; c * ) (thanks to the continuity of x(s; c) w.r.t. c).

(4) Follows from Proposition 1.

(5) Follows from Proposition 3.

Remark 1. Theorem 1 states the existence of an approximate cyclic solution E of Σ p0 ; in practice, when the integration time-step τ is sufficiently small, E is very close to the cyclic exact solution if such a cycle exists. But Theorem 1 does not guarantee by itself the existence of an exact cyclic solution. Although we know that certain exact solutions of Σ p0 lie in I 0 close to E, we do not know a priori if some of them are cyclic. To prove formally the existence of an exact cyclic solution, one could have to resort to specific notions of contractivity (see, e.g. [START_REF] Manchester | Transverse contraction criteria for existence, stability, and robustness of a limit cycle[END_REF]; [START_REF] Aminzare | Contraction methods for nonlinear systems: A brief introduction and some open problems[END_REF]; Capinski et al. (arXiv:1905.08116, 2020)). 3 The implementation has been done in Python and corresponds to a program of around 500 lines. The source code is available at lipn.univ-paris13.fr/ ~jerray/ parameter/. In the experiments below, the program runs on a 2.80 GHz Intel Core i7-4810MQ CPU with 8 GiB of memory. Given x 0 , τ, ε, |W|, the program searches heuristically for values of i ∈ N and T = kτ (for some k ∈ N) in order to verify B W ((i + 1)T ) ⊆ B W (iT ). This is illustrated in the following example. Example 2. Consider the system Σ p of Example 1 and the system Σ with uncertainty |W| = 0.2, and let p 0 = 0.4. Consider the set of initial conditions B(x 0 , ε) with x 0 = (1.46898897, -1.12538766) and ε = 0.1. For T = 6.368 = kτ with τ = 10 -3 , we find:

x(0) = (1.46898897, -1.12538766), δ W (0) = 0.1 3 In dimension d = 2, the existence of such a cycle is guaranteed by the Bendixon-Poincaré theorem.

x(T ) = (1.45389133, -1.10970492),

δ W (T ) = 2.437975771955677 x(2T ) = (1.45195751, -1.11004367), δ W (2T ) = 2.637252198149563 x(3T ) = (1.45230073, -1.11065593), δ W (5T ) = 2.652749677582014 (• • • ) x(13T ) = (1.49176459, -1.09165581), δ W (13T ) = 2.645700321656063 x(14T ) = (1.49470193, -1.09069935), δ W (14T ) = 2.644674179178997 x(15T ) = (1.49666192, -1.09054882), δ W (15T ) = 2.642640758083935 x(16T ) = (1.49747851, -1.09133825), δ W (16T ) = 2.6411743966251087 x(17T ) = (1.49699953, -1.09318874), δ W (17T ) = 2.631805337566251 x(18T ) = (1.49508659, -1.09620741), δ W (18T ) = 2.6266533499722065 x(19T ) = (1.49161476, -1.10048667), δ W (19T ) = 2.621622599187829 x(20T ) = (1.48647167, -1.10610402), δ W (20T ) = 2.614143707597656.
We have: B W ((i + 1)T ) ⊂ B W (iT ) for i = 14. The computation took 1120s of CPU time. This shows that the invariant I W has a torus shape for t ∈ [14T, 15T ], and contains the solutions of Σ p for each p ∈ [p 0 -|W|/2, p 0 + |W|/2]. Fig. 2 shows the cyclic approximate solution E of Σ p0 in the plane phase (top), and the radius δ(t) of the invariant I 0 of Σ p0 (bottom). As δ(t) becomes rapidly neglectible, the exact solutions of Σ p0 and E practically coincide within I 0 .

The invariant I W with |W| = 0.2 is depicted in green on Fig. 3, and it can be seen that I W encloses all the solutions of Σ p for p = p 0 , p 0 ± |W|/2 (represented in red, cyan and blue) for t ∈ [15T, 20T ], as stated by Proposition 2.

CONCLUSION

Given a parametric differential system parameter Σ p (with a parameter p taking its values in an interval of the form Q = [p 0 -|W|/2, p 0 + |W|/2]), we have introduced a simple condition of inclusion of sets (Euclidean balls) which guarantees that the approximate Euler solutions of Σ p0 are attracted by a limit cycle E, which is itself a cyclic approximate solution of Σ p0 . Moreover, we have introduced a toric set I 0 around E which is invariant for Σ p0 , and whose radius δ(t) becomes in practice quickly very small. This shows that the exact solutions of Σ p0 (bottom). This shows that the approximate and exact solutions of Σ p0 are rapidly almost identical.

exist in the close neighborhood of the cycle E and have themselves an almost cyclic behavior. Finally, we have constructed a compact set Σ W centered also on E, whose radius δ W (t) is now non negligible, which is invariant for every system Σ p (p ∈ Q). We hope that this method, illustrated on the example of a parametric Van der Pol system, opens an alternative practical way to the complex techniques based on contractivity, Lyapunov functions or Poincare maps that are used presently to show the existence of attractive limit cycles.

Figure 2 .

 2 Figure 2. Top: the approximate Euler solution of Σ p0 (p 0 = 0.4) in the phase plan near the limit cycle E. Bottom: the error bound δ(t) for t ∈ [0, 20T ] when |W| = 0. The invariant I 0 is a torus centered around E (top), and its radius δ(t) becomes rapidly neglectible(bottom). This shows that the approximate and exact solutions of Σ p0 are rapidly almost identical.

  ). Instead of computing them globally for S, it is advantageous to compute λ and γ locally depending on the subregion of S occupied by the system state during a considered interval of time. We now give a version of Proposition 1 with bounded uncertainty w(•) ∈ W, originally proved in LeCoënt et al. (2017a).

	Proposition 2. Le Coënt et al. (2017a) Consider a
	system Σ with bounded uncertainty of the form dx(t) dt f (x(t), w(t)) satisfying (H1).	=

Note that B W (0) = B 0 because x(0) = x 0 and δ ε,W (0) = ε.


Ci Liang
email: ci.liang@railenium.eu

Mohamed Ghazel

Olivier Cazier

El-Miloudi El-Koursi

Risk analysis on level crossings using a causal Bayesian network based approach

Keywords: Level crossings, causal statistic risk assessment, Bayesian networks, statistical analysis, railway safety;

Safety is a core issue in railway operation. In particular, Level Crossing (LX) safety is one of the most critical issues that railway stakeholders need to deal with. Accidents at European LXs account for about one-third of the entire railway accidents. They result in more than 300 deaths every year in Europe. However, due to non-deterministic causes, complex operation background and the lack of thorough statistical analysis based on detailed accident/incident data, the risk assessment of LXs remains a challenging task. In this paper, a general approach of Causal Statistic Risk Assessment based on hierarchical Causal Bayesian Networks (CSRA-CBN) is developed to analyze the various impacting factors which may cause accidents, and identify the factors which contribute most to the accidents at LXs, thus allowing for risk quantification. The detailed statistical analysis is carried out firstly according to the accident/incident database, then the CBN risk model is established based on the statistical results. In order to validate the effectiveness of this approach, we apply the CSRA-CBN approach on the basis of the accident data from SNCF, the French national railway operator. The CBN model allows for quantifying the causal relationships between the risk of accident occurring and the impacting factors considered. Moreover, the hierarchical structured modeling offers interesting benefits in terms of clarity, which makes it possible to highlight the complex factors influenced by a mass of parameters and identify the factors that contribute most to LX accidents. In addition, the main output of such a modeling system is conducive to improving safety at LXs, meanwhile, allowing for efficiently focusing on the effort/budget to make LXs safer.

Introduction

Literature review on level crossing study

Accidents at railway level crossings (LXs) often give rise to serious material and human damage and seriously hamper railway safety reputation, although the majority of accidents are caused by road user violations. LX safety is one of the most critical issues in the railway field which needs to be improved urgently [START_REF] Ghazel | Using stochastic Petri nets for level-crossing collision risk assessment[END_REF][START_REF] Mekki | Validation of a new functional design of automatic protection systems at level crossings with modelchecking techniques[END_REF][START_REF] Liu | Model-Based Diagnosis of Multi-Track Level Crossing Plants[END_REF]. In 2012, there were more than 118,000 LXs in the 28 countries of the European Union (E.U.) which corresponds to an average of 5 LXs per 10 line-km (ERA, 2014). Accidents at European LXs account for about one-third of the entire railway accidents. They result in more than 300 deaths every year in Europe [START_REF] Liu | Model-Based Diagnosis of Multi-Track Level Crossing Plants[END_REF]. In some European countries, accidents at LXs account for up to 50% of railway accidents (Ghazel and El-Koursi, 2014;Evans, 2011b). In the entire E.U. area, the overall number of fatalities per fatal accident in railways from 1990 to 2009 is 4.10, with no apparent long-term change over time (Evans, 2011a). In France, the railway network shows more than 15,000 LXs for 30,000 km of railway lines, which are crossed daily by 16 million vehicles on average, and around 13,000 LXs show heavy road and railway traffic (SNCF Réseau, 2011). In 2013, € 32 million was spent by the State, SNCF (the French national railway operator) and local authorities, to improve LX safety, and about € 36 million was invested in 2014. In recent years, a list of 300 prioritized LXs has been established and updated every year by SNCF Réseau, the French national railway infrastructure manager. Despite numerous measures taken to improve LX safety, SNCF Réseau counted 100 collisions at LXs leading to 25 deaths in 2014. This was half the total number of collisions at LXs a decade ago, but still too many (SNCF Réseau, 2015). If we want to significantly reduce the accidents and cost at LXs, it is necessary to develop further safety solutions to reach the most effective result.

In the literature pertaining to LX safety, a number of works can be found. In recent years, a systems analysis framework (Leveson, 2011;[START_REF] Read | Walking the line: Understanding pedestrian behaviour and risk at rail level crossings with cognitive work analysis[END_REF][START_REF] Wilson | Fundamentals of systems ergonomics/human factors[END_REF] and a psychological schema theory [START_REF] Salmon | The crash at Kerang: Investigating systemic and psychological factors leading to unintentional non-compliance at rail level crossings[END_REF][START_REF] Stanton | Exploring the psychological factors involved in the Ladbroke Grove rail accident[END_REF] have been used to analyze the contributory factors underlying the accidents occurring at LXs. The study presented by [START_REF] Salmon | The crash at Kerang: Investigating systemic and psychological factors leading to unintentional non-compliance at rail level crossings[END_REF] described a collision between a loaded semi-trailer truck and a train, occurring in North Victoria, caused by the truck driver crossing the LX while the barriers are closed. According to the investigation of OCI, the Office of the Chief Investigator, the truck driver in this study was not aware of the train and the activated state of the level crossing until it was too late to stop. A study conducted by [START_REF] Davey | The experiences and perceptions of heavy vehicle drivers and train drivers of dangers at railway level crossings[END_REF] discussed the intentional violation of vehicle drivers crossing LXs, particularly focusing on vehicle driver complacency due to the high level of familiarity. [START_REF] Tey | Measuring driver responses at railway level crossings[END_REF] conducted an experiment to measure vehicle driver response to LXs equipped with stop signs (passive), flashing lights and half barriers with flashing lights (active) respectively. In this study, the vehicle driver responses result from both the field survey and a driving simulator.

Based on the above investigations, the main available studies dealing with LX safety have tended to take a qualitative approach to understand the potential reasons causing accidents at LXs. They employed surveys, interviews, focus group methods or driving simulators, rather than collecting real data. However, the aforementioned approaches show some limits. For instance, the reaction of vehicle drivers in simulation scenarios could differ from that in reality, due to the different levels of feeling of danger. Although those qualitative approaches are beneficial to explore the potential factors causing accidents, they cannot quantify the contribution degree of these factors. Thereby, quantitative approaches are indispensable if we want to understand the impacting factors thoroughly and enable the identification of practical design and improvement recommendations to prevent accidents at LXs.

Literature review on risk modeling

Current risk analysis methodologies are required to deal with increasingly complex systems with a large amount of configuration parameters. Therefore, such methodologies shall be equipped with the following characteristics:

 Strong modeling ability  Easy to specify a risk scenario or a system  High computational efficiency In the domain of risk assessment, various methodologies are adopted for modeling and analyzing the systems or scenarios. Due to the combination of qualitative and quantitative analysis, the Fault Tree Analysis (FTA) developed by H.A. Watson at Bell Laboratories has been widely used for risk analysis since 1962 [START_REF] Ericson | Fault Tree Analysis[END_REF]. FTA is a deductive and top-down method aimed at analyzing the effects of initiating faults and events on a complex system and offering the designer an intuitive high-level abstraction of the system, which can also be efficiently applied through cooperating with techniques such as Binary Decision Diagrams [START_REF] Boudali | A discrete-time Bayesian network reliability modeling and analysis framework[END_REF]. Compared with Failure Mode and Effects Analysis (FMEA), which is an inductive and bottom-up analysis method aimed at analyzing the effects of single component or function failures on equipment or subsystems, FTA is good at showing how resistant a system is to single or multiple initiating faults. However, one obvious disadvantage of FTA is that it is not clear on failure mechanism, since the causal relationship between events is not a simple YES or NO (1 or 0). Therefore, FTA is prone to missing the possible initiating faults. In addition, traditional static fault trees cannot handle the sequential interaction and functional dependencies between components. Consequently, it is necessary to develop dynamic methodologies to overcome these weaknesses. Markov Chains (MCs) and their extensions have been mainly used for modeling complex dynamic system behavior and dependability analysis of dynamic systems. Two-state Markov switching multinomial logit models are introduced to explain unpredictable, unidentified or unobservable risk factors in road safety by [START_REF] Malyshkina | Markov switching multinomial logic model: an application to accident-injury severities[END_REF]. Although MCs can elaborate the statistical state transition of different variables, they cannot formalize causal relationships between various events.

Afterwards, risk analysis based on formal modeling flourished. In order to compare the effectiveness of two main Automatic Protection Systems (APSs) at LXs: two-half-barrier APS and four-half-barrier APS, Generalized Stochastic Petri Nets (GSPNs) were used to analyze the aleatory fluctuations of various parameters involved in the dynamics within the LX area by Ghazel and El-Koursi (2014). Over the last few years, Bayesian network (BN), a method of reasoning using probabilities, has been an increasingly popular method used for risk analysis of safety-critical systems or large and complicated dynamic systems [START_REF] Chemweno | Development of a risk assessment selection methodology for asset maintenance decision making: An analytic network process (ANP) approach[END_REF]. To obtain proper and effective risk control, risk planning should be performed based on risk causality which can provide more risk information for decision making. Thus, a model using BNs with causality constraints (BNCC) for risk analysis was propose by [START_REF] Hu | Software project risk analysis using Bayesian networks with causality constraints[END_REF]. The BN technique shows the flexibility of modeling, strong modeling power, high computational efficiency and, most importantly, the outstanding advantage involving causality analysis of both forward inference and reverse inference [START_REF] Weber | Overview on Bayesian networks applications for dependability, risk analysis and maintenance areas[END_REF].

Therefore, based on the above investigation of risk analysis methodologies, a general approach of Causal Statistic Risk Assessment based on hierarchical Causal Bayesian Network (CSRA-CBN) is developed in this article to solve the risk assessment issue of LXs. We make a thorough statistical analysis on the accident/incident data involving LXs from SNCF and use the statistical results as the input source of CBN, then quantify the risk level of various potential factors and identify the factors which contribute most to the accidents at LXs, as well as their combined impact on LX safety.

The rest of this paper is organized as follows: In Section 2, the theory of probability which is the basic principle of CBN is stated firstly, then the integrated framework of CSRA-CBN is elaborated. In Section 3, we apply the CSRA-CBN approach on the LXs in France based on the accident/incident data from SNCF to validate the effectiveness of this approach. Statistical analysis is performed, risk model is established and simulation results are discussed in this section. Here, the GeNIe is used as a CBN modeling tool to formalize the causal relationships between accidents and the considered parameters. Finally, some concluding remarks and further directions are given in Section 4.

Causal statistic risk assessment based on causal Bayesian networks

Causal Bayesian networks

In the concept of Socio-technical systems (Leveson, 2011), potential hazards which may lead to accidents, including equipment failures, human errors, and other uncertain factors, such as environment aspects. In fact, the causality between the accidents and these uncertain factors exists, see Fig. 1. Causality plays an important role in the process of inference. A functional intelligent identification model shall have the ability of making inference based on the causal knowledge. Causal Bayesian Network (CBN) employed to model causal relationships is a graphical model which consists of a structure and a set of parameters [START_REF] Yehezkel | Bayesian network structure learning by recursive autonomy identification[END_REF]. CBN = (P, G), where P represents the parameters of prior (conditional) probabilities that quantify the arcs and G represents the model structure. G = (V, A), which is a Directed Acyclic Graph (DAG) comprised of nodes Vs and arcs As. Nodes Vs represent random variables and directed arcs As between pairs of nodes represent dependencies between the Vs [START_REF] Langseth | Bayesian networks in reliability[END_REF].

CBN is developed on the theory of probability. Assume that there is a set of mutually exclusive events B1, B2, …,Bn and event A, in addition, P(A) can be expressed as follows:

1 ( ) ( ) ( | ) n i i i P A P B P A B    (1) According to Bayes formula: 1 ( ) ( | ) ( | ) ( ) ( | ) i i i n j j j P B P A B P B A P B P A B   
(2) Therefore, Eq. ( 2) can be converted into:

(

) ( | ) ( | ) () i i i P B P A B P B A PA 
(3) P (Bi) is the Prior Probability, P (Bi | A) is the posterior (conditional) probability. When evidence of the nodes is given in a CBN, the posterior probability can be calculated.

For any set of random variables in a CBN, the joint distribution can be calculated through conditional probabilities using the chain rule shown as follows:

1 1 n 1 1 n 1, ( , , = ) ( = | = , , = ) n n v v v v n v v V P X x X x P X x X x X x     (4)
Due to the conditional independence, Xv only relates to its parent nodes Pa (Xv) and is independent of other nodes. So the last equation can be written as:

1 1 n 1, ( , , = ) ( = | ( )) n n v v v v v V P X x X x P X x Pa X    (5)
For more details about CBN, the reader can refer to a tutorial book of Bayesian networks edited by [START_REF] Pourret | Bayesian networks: a practical guide to applications[END_REF].

CSRA-CBN Framework

A paradigm for modeling should view an influential network not merely as passive parsimonious codes for storing factual knowledge, also as a computational architecture for reasoning about the knowledge. The CSRA-CBN is developed based on this principle. It means that the links in the network should be treated as the only pathways and activation units that direct and propel the flow of data in the process of querying and updating causal knowledge.

The framework of the CSRA-CBN approach can be illustrated in Fig. 2. The process of applying the CSRA-CBN approach consists in following five stages:

1) Identifying the risk scenario that we wish to focus on; Before performing risk assessment, in order to set our research target, we must make a clear definition of the risk scenario boundary, and concentrate on this defined risk scenario to ensure that the follow-up study does not deviate from our original intention.

2) Collecting accident/incident data related to the defined risk scenario; For predicting accident probability, assessing risk level and making improvement measures, accident/incident data related to the defined risk scenario need to be collected through monitoring the risk scenario. These data should be recorded in a dedicated database and as the source of the statistical analysis in step 3).

3) Developing statistical analysis on the accident/incident data; Statistical analysis on the accident/incident data aims to identify impacting factors, further explore their statistic characteristics. The outcomes of statistical analysis should be as the importing evidence of CBN risk model. 4) Establishing the CBN risk model based on the statistical analysis;

In this step, we need to build the CBN model structure with respect to the causal relationships between identified impact factors and accidents, then employ the statistical results as the input evidence of the risk model.

5) Assessing risk level based on risk models; In this step, accident probability will be quantified and risk level will be assessed. Moreover, impacting factors which contribute most to the entire risk we concerned will be identified to provide improvement measures for reducing the risk concerned and the cost brought by accidents. In the fourth stage, the hierarchical network is adopted when establishing the risk model. The causal relationships between top nodes and intermediate nodes which cause the top nodes directly are established in the first layer of CBNs. In the second layer, the direct causal relationships between intermediate nodes and their root nodes are established. There is no indirect causality existing in any layer. Moreover, the network structure in each layer is single connected (SC), such as "Start" topology or "Tree" topology, which means one and only one path exists between any two nodes (without regard to the arc direction).

Application and validation of CSRA-CBN

In this section, the CSRA-CBN approach is applied on the LXs in France based on the accident/incident data collected by SNCF to assess the risk level and the probability of accidents occurring with regard to the various impacting factors taken into account, and further identify the important factors which contribute most to accidents occurring.

Risk scenario identification and data collection

There are four types of LXs in France [START_REF] Sncf | Research on the material of level crossing in[END_REF] according to the data of SNCF: a) SAL4: Automated LXs with four half barriers and flashing lights; b) SAL2: Automated LXs with two half barriers and flashing lights; c) SAL0: Automated LXs with flashing lights but without barriers; d) Crossbuck LXs, without automatic signaling. As shown in Table 1, SAL2 (more than 10,000) is the most widely used type of LX in France. Moreover, more than 4,000 accidents at SAL2 contributed most to the total number of accidents at LXs from 1974 to 2014. For the above reasons, "road users crossing French SAL2 LXs when trains are coming" is identified as the risk scenario. Detailed LX accident/incident data from 1974 to 2014 is collected by SNCF and input into a special database which we can access. Thereby, one analysis will consider the accident/incident data pertaining to the 40 years.

Statistical analysis

According to the statistics of SNCF, geographical region, transport mode and traffic moment are essential factors which should be considered when improving the LX safety. Nevertheless, these factors were not thoroughly analyzed by the previous works. In this subsection statistical analysis with respect to traffic moment, different kinds of transport mode and different geographical regions is performed to establish the potential relationships between these principal impacting factors and accidents, and explore their distribution characteristics.

There are 21 geographical administrative regions in mainland France. Accidents which are caused by main types of transport mode: 1) motorized vehicle (MV), 2) pedestrian or bicycle (PB), are considered respectively to allow for making statistical analysis in different regions. The number of train-MV collisions, train-PB collisions and the number of SAL2 in 21 regions are presented respectively in Fig. 3.

General risk analysis

The general frequency of accidents occurring at each SAL2 per year are used to represent the general risk level involving total accidents, MV accidents, and PB accidents in different regions during the last 40 years. We can calculate the general frequency through the Eq. ( 6) in Definition 1:

 Definition 1 _ , 1, 2, ..., 21; _ 2 _ i i G i Nb acc Fi Nb SAL Nb year   (6)
where 𝐹 G i represents the general frequency in i th region; Nb_acci represents the number of accidents occurring in i th region, Nb_SAL2i represents the number of SAL2 in i th region, and Nb_year represents the number of years of the considered period. Three kinds of general frequency related to total accidents, MV accidents, and PB accidents will be calculated. Correspondingly, the number of total accidents, MV accidents, and PB accidents are presented as Nb_acci respectively when calculating these three kinds of general frequency. Now that these three kinds of general frequency in each region are determined, maps of French regions with the frequency values presented are generated to show the frequency distribution in different regions. As shown in Fig. 4a, the general frequency value of total accidents in the red region (greater than 0.02) is the highest. The general frequency value of total accidents in the orange region (between 0.02 and 0.01) is the second highest, and the general frequency values of total accidents in the green region (less than 0.01) is the lowest. When we analyze the frequency figures in detail, we find that the risk is most serious in Île-de-France with a frequency of more than 0.02; Languedoc-Roussillon takes the second place with the frequency of about 0.017 followed by Provence-Alpes-Côte-d'Azur with the frequency of about 0.016. On the other hand, Limousin has the lowest risk with the frequency of about 0.005. Haute-Normandie and Basse-Normandie occupy the second and the third places of lowest risk successively. Turning to Fig. 3, the number of accidents occurring in Île-de-France is not the highest, but, when considering the corresponding number of LXs, this leads to the most serious risk. Conversely, although accidents in Rhône-Alpes are the most with the number of 521, the risk is not the highest due to a large number of LXs in this region. Limousin, Haute-Normandie and Basse-Normandie have significantly fewer accidents than any other region, less than 100 in 40 years. That is why they have the lowest risk.

In Fig. 4b and Fig. 4c, the general frequency of accidents caused by motorized vehicles, pedestrians and bicycles in different regions is shown respectively. Considering the accidents caused by motorized vehicles, the distribution of frequency in different regions in Fig. 4b is relatively consistent with the distribution shown in Fig. 4a. The only exception is Champagne-Ardenne. However, in Fig. 4c, as for the accidents caused by pedestrians and bicycles, the distributions of frequency in different regions are very different from the distribution shown in Fig. 4a. By analyzing the results shown in Fig. 4 and Fig. 5, we conjecture that the main mode causing accidents at SAL2 is the motorized vehicle. Moreover, as the frequency of accidents involving motorized vehicle collisions increases, the total accident frequency in the overwhelming majority of regions increases accordingly, as illustrated in Fig. 5. Further discussion and validation will be elaborated in section 3.3.

Normalized risk analysis based on the traffic moment

As mentioned in section 3.2.1, the motorized vehicle is the main transport mode causing accidents at SAL2. In this section, we will focus on the impact of traffic moment on the whole risk level related to MV accidents. We suspect that this parameter is one of the main parameters impacting LX risk level. Indeed, traffic moment gives the combined traffic (train / MV) at the LX and is defined as follows:

 Definition 2 = Traffic moment road traffic frequency railway traffic frequency 

In this section, we use "M" for short to denote the traffic moment. In order to analyze the impact of M, we use the square root of M at SAL2 to normalize the frequency of accidents caused by train-MV collisions, which can represent the normalized frequency in different regions. The formula we adopt to establish the normalized frequency is the following:

 Definition 3 _ 2 _ 1 __ , 1, 2, ..., 21; _ 2 _ _ i i Nb SAL reg j j j N i Nb acc MV M Fi Nb SAL reg Nb year     (8)
where F N 𝑖 represents the normalized frequency of accidents caused by train-MV collisions in i th region; Nb_acc_MVj represents the number of accidents caused by train-MV collisions at j th SAL2 in i th region; Mj represents the moment at the j th SAL2 in i th region; Nb_SAL2_regi represents the number of SAL2 in i th region, and Nb_year represents the number of years of the considered period.

As shown in Fig. 6, the normalized frequency caused by train-MV collisions occurring at SAL2 LXs in each region is computed. Unlike Fig. 4b, Fig. 6 illustrates that the risk level is the highest in Île-de-France with the normalized frequency of 7.31 ×10 -4 ; Alsace takes the second place followed by Languedoc-Roussillon and Provence-Alpes-Côted'Azur. Through detailed analysis, we find that the general frequency of MV accidents in these 4 regions (cf. Fig. 4a) is also the highest.

For the purpose of further analyzing the impact of M on the risk level in each region, SAL2 LXs are classified by the category of M. In order to make the number of SAL2 in each M group to be as far as possible similar to any other group in each region, the moment groups have been defined in such a way that the number of SAL2 LXs in every group belongs to the interval [102,155]. For example, in Table 2, there are three M categories in Auvergne, which are "0 ≤ M < 750", "750 ≤ M < 5,000", "5,000 ≤ M ≤ 401,850", with the corresponding numbers of SAL2: 143, 144, 147, respectively. Besides, 401,850 is the maximum M in this region. In this way, we can make risk analysis with regard to these SAL2 groups, thus making it possible to highlight the risk level related to different categories of M. Logically, if the same number of MV accidents or more MV accidents occur per year at the SAL2 with small M than at the SAL2 with large M, this indicates that the SAL2 with small M has higher risk level than the SAL2 with large M.

The normalized frequency of collisions regarding the different M categories in each region is shown in Table 2. The M categories marked in red color have the most serious risk in various regions. We find that SAL2 LXs in the category of smallest M have the highest risk in 86% of the regions, except in Auvergne, Champagne-Ardenne and Nord-Pas-de-Calais where the SAL2 LXs with the second smallest M have the highest risk. On the other hand, SAL2 LXs in the category of largest M have the lowest risk in 81% of the regions. These statistical results prove the logical inference mentioned above that the SAL2 with small M has higher risk level than the SAL2 with large M, when considering the ratio of the number of MV accidents occurring per year at the SAL2 to the traffic moment of the SAL2. Further discussion and validation will be elaborated in section 3.3.

M [0, 1,250) 0.000442 M [1,250, 7,200) 0.000134 M [7,200, 1428,000] 0.000152 M [a, b): a ≤ moment < b; M [c, d]: c ≤ moment ≤ d.

Risk model establishment and discussion

Risk model establishment

In this section, we establish the CBN risk model based on the analysis made so far, particularly the statistical results in section 3.2 to assess the risk level related to geographical region, transport mode and traffic moment. The "GeNIe" software tool is used as the CBN modeling mean.

Firstly, the structure of risk model is established with respect to the causal relationships between the impacting factors and accidents. As shown in Fig. 7, the initial hierarchical model contains two layers: Sub-layer 1 is the submodel in the"Risk_MV" node and Sub-layer 2 is the sub-model in the "Risk_PB" node. States of moment group are defined as "M_Region" (e.g., the "M_Alsace" node) of Sub-layer 1; states of MV accidents and PB accidents occurring in different regions are defined in region nodes of Sub-layer 1 and Sub-layer 2 (e.g., the "Alsace" node), respectively.

Assume that all nodes only have two states in CBNs. If one node has n parent nodes, the combinations of states in this node will be up to 2 (n+1) . In our model, some parent nodes contain more than two states. This means there are more than 2 22 combinations in "Risk_MV" and "Risk_PB" nodes. Moreover, more than half of the combinations are invalid (some cases do not exist in reality, e.g., riski only corresponds to the risk value in region i and it cannot exist when regarding more than two different regions). Thereby, in order to avoid the state space explosion (SSE), we optimize the structure of the model and eliminate the invalid combinations of states, then input evidence into the model based on the statistics in Section 3.2. As shown in Fig. 8, "Part 3" is used to replace "Part 1" in Fig. 7 and "Part 4" is used to replace "Part 2" in Fig. 7. Table 3 offers a summary of the states in each node of the optimized risk model. In order to avoid the SSE problem, we determine 3 groups of moment which contain all the SAL2 LXs and the number of SAL2 LXs in each group is almost equal.

Accident prediction and risk level assessment

We can estimate the probability of accidents occurring with regard to the impact of these factors interacting with each other on accidents. As shown in Fig. 9, the quantified probability is obtained through forward inference. The estimated probability of accidents occurring impacted by the interaction of geographical region, transport mode and traffic moment is lower than 0.0065, and this result is consistent with the statistics developed by SNCF.

Subsequently, "Risk_Accidents = True" is configured as the target evidence. In this way, we can assess the level of various impacting factors through reverse inference. Detailed results are offered in Fig. 10. We can find that when accidents have already occurred, the risk of accidents caused by motorized vehicles contributes 98% to the overall risk. However, the risk of accidents caused by pedestrians and bicycles only contributes 2% to the overall risk. Moreover, in terms of geographical region, the normalized risk level of MV accidents in Île-de-France is the highest, which induces 14% of the entire MV risk. The risk of MV accidents in Alsace takes the second place followed by Languedoc-Roussillon and Provence-Alpes-Côte-d'Azur. In addition, the risk level of PB accidents is also the highest in Île-de-France. These outcomes are consistent with the statistics in Section 3.2. Paris is indeed included in Île-de-France. In addition, Paris holds 6 main railway stations, among which the north railway station "Gare du Nord" is one of the most important railway stations in the world. This explains at least partially the fact that the risk of accidents occurring is the highest in Île-de-France.

In terms of moment group, the smallest moment group "M 0_1485" has the highest risk level due to inducing 61% of the risk caused by motorized vehicles and the largest moment group "M 10500_244000" has the lowest risk level for only inducing 19% of the risk caused by motorized vehicles. This result also strongly supports the logic inference mentioned in Section 3.2.2 that the SAL2 with small moment has higher risk level than the SAL2 with large moment when considering the ratio of the number of MV accidents occurring per year at the SAL2 to the moment value of the SAL2.

Part 3

Part 4

Fig. 8. The optimized risk model. Furthermore, the sensitivity analysis is performed to elaborate the sensitivity for "Risk_Accidents = True" related to various combination conditions. Through forward running the model, the probability of "Risk_Accidents = True" is obtained as 6.45749×10 -3 . Here, we use P (Risk_Accidents = True) and P (Risk_MV = False) to represent the probability of "Risk_Accidents = True" and "Risk_MV = False", respectively. In this case, the current value Vc of P (Risk_MV = False | Accidents_MV_region = Île-de-France 1, Moment = M 0_1485) is 0.95668.

The Sensitivity Tornado Diagram (STD) is offered by Fig. 11. The horizontal bars are viewed as two parts according to the central vertical axis P (Risk_Accidents = True) = 6.45749 ×10 -3 . P (Risk_Accidents = True) in the green part is lower than 6.45749 ×10 -3 and P (Risk_Accidents = True) in the red part is higher than 6.45749 ×10 -3 . As the value range of impacting parameter is extended to the interval [Vc -0.1Vc, Vc + 0.1Vc] (setting the parameter spread degree as 0.1), the values of P (Risk_Accidents = True) belong to an interval [Vmin, Vmax] The above analysis result is logical since the main transport mode causing accidents is MV and the probability of "no MV accidents occurring" is very high (after all, MV accidents occurring is a small probability event in reality). A small change in the probability of "Risk_MV = False" can cause a dramatic change in the probability of "Risk_Accidents = True". This analysis result further proves that, if we want to decrease the risk at SAL2 LXs significantly, MV accidents shall be reduced considerably. In particular, more attention shall be paid to MV accidents related to the smallest moment groups. In reverse, all the above outcomes have verified the validity of the risk model by comparing with the statistics in Section 3.2.

Conclusions

The contributions of the present study are as follows: an integrated risk assessment approach called CSRA-CBN is developed. In the framework of CSRA-CBN, the risk scenario is firstly identified. Then, thorough statistical analysis is performed based on the accident/incident data of the risk scenario. Finally, the risk model is established based on the statistical analysis to quantify the probability of accidents occurring and the contribution of various impacting factors, thus allowing for determining the important contributors to the overall risk.

The application of CSRA-CBN to investigate French LX safety allows us to identify the main LX risk factors and quantify their respective contributions to the overall LX risk. Thanks to the CSRA-CBN approach, various significant factors pertaining LX safety, i.e., the involved road transport mode, geographical regions and the traffic moment, have been analyzed meticulously. More specifically, the statistical analysis on general frequency of the entire LX accidents, train-MV LX accidents and train-PB LX accidents in different French regions is performed. Then, further investigation is carried out to analyze the train-MV LX accident frequency normalized by the square root of "traffic moment" according to various "traffic moment" categories in different regions. Afterwards, the risk model is established to assess the risk level influenced by the geographical region, the transport mode and the traffic moment, based on the causal relationships between these impacting factors and accidents occurring and the aforementioned statistical results. The sensitivity analysis based on the risk model is performed to highlight the greatest contributor to the overall risk at LXs. The greatest contributor is identified among the combinations of various impacting factors in different situations. The effectiveness of CSRA-CBN approach is validated by the obtained statistical results, the risk model outputs and the sensitivity analysis outcomes of the whole risk model.

It is worth noticing that the CSRA-CBN is a general approach which can be applied on various risk assessment studies. In future works, we will consider further impacting factors and bring into play Bayesian risk models to quantify the causal relationships between the factors and accidents occurring, and assess their impact on the risk level. The first version of risk model will be enriched while integrating further factors according to real data and causality analysis that we will perform on the basis of fault events.

Fig. 1 .

 1 Fig. 1. Reasoning between hazards and accidents.

Fig. 2 .

 2 Fig. 2. Frame diagram of the CSRA-CBN approach.

Fig. 3 .

 3 Fig. 3. Accidents caused by different transport modes at SAL2 in 21 French regions from 1974 to 2014.

  Fig. 4. (a) Frequency of total accidents distributing in different regions; (b) Frequency of MV accidents distributing in different regions; (c) Frequency of PB accidents distributing in different regions.

Fig. 5 .

 5 Fig. 5. Comprehensive frequency of accident in different regions.

Fig. 6 .

 6 Fig. 6. Normalized frequency of MV accidents distributing in different regions.

Fig. 7 .

 7 Fig. 7. The initial hierarchical risk model.

Fig. 9 .

 9 Fig. 9. Estimated probability of accidents occurring.

  and distribute with respect to the values of impacting parameter in the whole range [Vc -0.1Vc, Vc + 0.1Vc], accordingly. The STD highlights the greatest contributor "Risk_MV = False | Accidents_MV_region = Île_de_France 1, Moment = M 0_1485" (the combination of various impacting factors) to the overall risk, that is no MV accidents occurring under the condition of "Accidents_MV_region = Île-de-France 1, Moment = M 0_1485". Here, the value range of P (Risk_MV = False | Accidents_MV_region = Île_de_France 1, Moment = M 0_1485) is in the interval [0.86101, 1], the value of P (Risk_Accidents = True) belongs to the interval [0.00586, 0.00778] distributing with respect to the values of P (Risk_MV = False | Accidents_MV_region = Île_de_France 1, Moment = M 0_1485) in the whole range [0.86101, 1]. The second greatest contributor is "P (Risk_MV = False | Accidents_MV_region = Alsace 4, Moment = M 0_1485)" (i.e., no MV accident occurring in the moment group M 0_1485 of Alsace); P (Risk_MV = False | Accidents_MV_region = Provence_Alpes 3, Moment = M 0_1485), that is, no MV accident occurring in the moment group M 0_1485 of Provence-Alpes-Côte-d'Azur, and P (Risk_MV = False | Accidents_MV_region = Languedoc_Roussillo 2, Moment = M 0_1485), that is, no MV accident occurring in the moment group "M 0_1485" of Languedoc-Roussillon, almost have the same degree of influence on the overall risk and take the third and fourth places. Besides, the overall risk is not very sensitive to other impacting parameters.

Fig. 10 .

 10 Fig. 10. Risk level assessment of various impacting factors.

Fig. 11 .

 11 Fig. 11. Sensitivity Tornado Diagram of "Risk Accidents = True".

Table 1 .

 1 Accidents at different types of LXs in France from 1974 to 2014.

	Type of LX	Number	# Accident
	SAL4	>600	>400
	SAL2	>10,000	>4,000
	SAL0	>60	>50
	Crossbuck LX	>3,000	>700

Table 2 .

 2 Normalized frequency of MV accidents at SAL2 in each region according to the category of "moment".

	Normalized

frequency of MV accidents at SAL2 in each region according to the category of "moment"

				0.000072		0.000095	0.000065
	Basse-	M [0, 2,100)	M [2,100,
	Normandie	0.000147	555,600]
			0.000067
	Poitou-	M [0, 500)	M [500, 2,200)	M [2,200,	M [10,500,
	Charentes	0.000302	0.000064	10,500)	1169,782]
				0.000098	0.000078
	Bourgogne	M [0, 1,000)	M [1,000,	M [8,000,
		0.000224	8,000)	720,000]
			0.000130	0.000061
	Limousin	M [0, 700)	M [700, 5,000)	M [5,000,
		0.000269	0.000058	337,524]
				0.000055
	Pays de la	M [0, 500)	M [500, 1,400)	M [1,400,	M [3,900, 15,000)	M [15,000,
	Loire	0.000175	0.000121	3,900)	0.000122	815,400]
				0.000111		0.000062
	Champagne-	M [0, 420)	M [420, 1,800)	M [1,800,	M [8,000,
	Ardenne	0.000129	0.000285	8,000)	256,410]
				0.000099	0.000185
	Centre-Val	M [0, 350)	M [350, 1,400)	M [1,400,	M [4,300, 18,000)	M [18,000,
	de Loire	0.000201	0.000121	4,300)	0.000087	688,320]
				0.000092		0.000059
	Aquitaine	M [0, 600)	M [600, 1,600)	M [1,600,	M [3,800, 11,000)	M [11,000,	M [40,000,
		0.000162	0.000135	3,800)	0.000093	40,000)	585,644]
				0.000095		0.000091	0.000055
	Bretagne	M [0, 800)	M [800, 4,200)	M [4,200,	M [18,000,
		0.000324	0.000117	18,000)	2244,000]
				0.000068	0.000044
	Rhône-Alpes M [0, 320)	M [320, 1,400)	M [1,400,	M [3,600, 11,000)	M [11,000,	M [35,000,	M
		0.000297	0.000122	3,600)	0.000101	35,000)	100,000)	[100,000,
				0.000106		0.000080	0.000054	1121,088]
								0.000047
	Franche-	M [0, 2,500)	M [2,500,
	Comté	0.000143	271,936]
			0.000101
	Picardie	M [0, 600)	M [600, 2,550)	M [2,550,	M [14,000,
		0.000314	0.000135	14,000)	817,000]
				0.000127	0.000101
	Nord-Pas-	M [0, 450)	M [450, 3,000)	M [3,000,	M [10,000,	M [45,000,
	de-Calais	0.000079	0.000203	10,000)	45,000)	1477,229]
				0.000082	0.000063	0.000059
	Lorraine	M [0, 1,350)	M [1,350,	M [8,000,
		0.000250	8,000)	402,500]
			0.000122	0.000099
	Provence-	M [0, 5,700)	M [5,700,
	Alpes-Côte-	0.000212	921,250]
	d'Azur		0.000094
	Île-de-	M [0, 11,500)	M [11,500,
	France	0.000285	752,854]
			0.000099
	Alsace	M [0, 12,000)	M [12,000,
		0.000513	767,088]
			0.000113
	Languedoc-
	Roussillon
	Auvergne	M [0, 750)	M [750, 5,000)	M [5,000,
		0.000049	0.000100	401,850]
				0.000078
	Haute-	M [0, 500)	M [500, 4,500)	M [4,500,
	Normandie	0.000190	0.000089	382,788]
				0.000072
	Midi-	M [0, 460)	M [460, 1,150)	M [1,150,	M [3,000, 9,000)	M [9,000,	M [35,000,
	Pyrénées	0.000202	0.000068	3,000)	0.000066	35,000)	1169,782]

Table 3 .

 3 States of each node in the optimized risk model.

	Nodes	States
	Accidents_MV_region	Île_de_France 1, Languedoc_Roussillon 2, Provence_Alpes 3, Alsace
	Accidents_PB_region	4, Rhone_Alpes 5, Lorraine 6, Nord_Pas_de_Calais 7, Picardie 8,
		Franche_Comte 9, Aquitaine 10, Champagne_Ardenne 11, Bretagne
		12, Midi_Pyrenees 13, Pays_de_la_Loire 14, Bourgogne 15,
		Centre_Val_de_Loire 16, Auvergne 17, Poitou_Charentes 18,
		Basse_Normandie 19, Haute_Normandie 20, Limousin 21;
	Moment	M 0_1485 (moment between 0 and 1485), M 1485_10500 (moment
		between 1485 and 10500), M 10500_244000 (moment between 10500
		and 244000);
	Risk_MV	False, True;
	Risk_PB
	Risk_Accidents

Acknowledgements

This work has been in the framework of "MORIPAN project: Modèle de RIsque pour les PAssages à Niveau" within the Railenium technological research institute , in cooperation with the National Society of French Railways Networks (SNCF Réseau) and the French Institute of Science and Technology for Transport, Development and Networks (IFSTTAR).