Thibaut Jacqmin

Bess Fang

Tarik Berrada

Tommaso Roscilde

Isabelle Bouchoule

Momentum distribution of 1D Bose gases at the quasi-condensation crossover: theoretical and experimental investigation

Keywords: numbers: 03, 75, Hh, 67, 10, Ba

We investigate the momentum distribution of weakly interacting 1D Bose gases at thermal equilibrium both experimentally and theoretically. Momentum distribution of single 1D Bose gases is measured using a focusing technique, whose resolution we improve via a guiding scheme. The momentum distribution compares very well with quantum Monte Carlo calculations for the Lieb-Liniger model at finite temperature, allowing for an accurate thermometry of the gas that agrees with (and improves upon) the thermometry based on in situ density fluctuation measurements. The quasicondensation crossover is investigated via two different experimental parameter sets, corresponding to the two different sides of the crossover. Classical field theory is expected to correctly describe the quasi-condensation crossover of weakly interacting gases. We derive the condition of validity of the classical field theory, and find that, in typical experiments, interactions are too strong for this theory to be accurate. This is confirmed by a comparison between the classical field predictions and the numerically exact quantum Monte Carlo calculations.

I. INTRODUCTION

Correlation functions are essential to describe manybody systems.

In particular, the first-order correlation function, or equivalently its Fourier transform, the momentum distribution, is an important observable, since it witnesses various phenomena such as the Bose-Einstein condensation transition, the Berezinskii-Kosterlitz-Thouless transition [1] or the Mott transition [2]. In the case of one-dimensional (1D) homogeneous gases, although one does not expect phase transitions, the correlation functions contain valuable information about the system. While the equation of state of 1D homogeneous Bose gases has been extensively studied in various regimes [3][4][5][6] and compared with exact theories [7], the momentum distribution has remained largely unexplored. It was measured for quasi-1D quasicondensates using Bragg spectroscopy [8][9][10], and more recently it was investigated using the focusing technique for quasi-1D Bose gases in the crossover from the ideal Bose gas to quasi-condensate regimes [4,11]. Density ripples that appear at the near field of a freely expanding quasi-1D quasi-condensate have also been studied [12,13]. Although related to the second-order correlation function, they also provide information about the first-order correlation function within the quasicondensate theory. The momentum distribution of an array of strongly correlated 1D Bose gases [14] constitutes the only measurement on a truly 1D system to our knowledge. From a theoretical point of view, the momentum distribution of 1D Bose gases with repulsive contact interactions is not known exactly over the entire phase diagram, but a few results have been established: the mean kinetic energy can be extracted from the exact Yang-Yang thermodynamics [11], and the short-range correlations are responsible for a momentum tail scaling as 1/p 4 [15][16][17]. The momentum distribution is known in the asymptotic regimes such as the ideal Bose gas and the quasi-condensate regimes, and a classical field approximation has been proposed to account for the crossover between them [11,[18][19][20].

In this work we investigate the momentum distribution of weakly interacting, purely or almost purely 1D Bose gases in various regimes around the quasi-condensation crossover. On the theory side, we find a criterion for the classical field theory [18,19] to be accurate at the quasi-condensation crossover. Interactions are however too strong in our experiment to fulfill this criterion, and we implement exact Quantum Monte Carlo (QMC) simulations for the equilibrium behavior of the Lieb-Liniger model. We find that the measured momentum distribution agrees very well with QMC calculations. Moreover, the temperature extracted from a fit of the momentum distribution to the QMC calculations is in agreement with that obtained from in situ density fluctuation measurements [3,6]. This paper is organised as follows. In Sec. II we recall the classical field approximation and investigate the conditions under which it is valid to describe the quasicondensation crossover. We show that, in most experimental situations, the gas is too strongly interacting for this description to be accurate. We then present our QMC calculations, and compare them with the classical field prediction in the asymptotic limit. In Sec. III, we describe the experimental setup and the focusing technique we use to measure momentum distribution. We discuss how we improve the resolution of this technique by means of a guiding scheme. In Sec. IV, we present our experimental results that span both the quasi-condensate regime (purely 1D) and the degenerate ideal Bose gas regime (almost 1D), and compare them to QMC results. We conclude and discuss the prospective work in Sec. V.

II. THEORETICAL PREDICTIONS

A. Model Hamiltonian and exact solution

We consider a 1D homogeneous Bose gas with repulsive contact interactions, described by the Lieb-Liniger Hamiltonian, whose grand-canonical expression is

Ĥ = dz - 2 2m Ψ † ∂ 2 ∂z 2 Ψ + g 2 Ψ † Ψ † Ψ Ψ -µ Ψ † Ψ , ( 1
)
where z is the position, Ψ ≡ Ψ(z) is the field operator, g is the coupling constant, m is the mass of one boson and µ is the chemical potential. The thermal equilibrium state of the system is completely characterised by the density ρ and the temperature T , or by the dimensionless interaction parameter, γ = mg 2 ρ , and the reduced temperature, t = 2 2 kB T mg 2 . The phase diagram [3] in the parameter space (γ, t) is given in Fig. 1, where lines should be understood as crossovers. In this article, we focus on the quasi-condensation crossover, which occurs, for γ ≪ 1, around the line γ co = t -2/3 (solid line). The Lieb-Liniger model has the remarkable property of being integrable, giving access to exact results valid over the entire phase diagram. For example, the equation of state is known exactly through the Yang-Yang theory [7]. Yet, the calculation of the first order correlation function g (1) (z) = Ψ † (z) Ψ(0) from the Bethe-Ansatz solution for finite γ and at finite temperature is still a subject of active research. Approximate analytic theories for the correlation functions exist in the limiting cases of ideal Bose gas regime (where the ideal Bose gas theory applies) and quasi-condensate regime (where the Bogoliubov theory applies [17]). However, the behavior of g (1) (z) at the crossover between those two regimes lacks a quantitative description.

B. Classical field theory

The full quantum many-body problem of Eq. ( 1) is notoriously complex. However, since thermal fluctuations are expected to dominate at the quasi-condensation crossover, it seems appropriate to simplify the problem using a classical field approximation, where the quantum field operators Ψ(z) and Ψ † (z) are replaced by the complex fields Ψ(z) and Ψ * (z) [18,19]. This strategy has been pursued in [11,20,[START_REF] Cockburn | [END_REF]. Such an approach is particularly attractive for 1D homogeneous gases: contrary to Phase diagram of the 1D Bose gas in the parameter space (γ, t). The ideal Bose gas to quasi-condensation crossover occurs for γ ≃ γco = t -2/3 (solid line). Within the ideal Bose gas regime (light gray area), the crossover from the classical to the degenerate sub-regime occurs for γ ≃ γ d = t -1/2 (dashed line). The region γ, t ≫ 1 (dark gray area) is the fermionised regime. The segments correspond to the regions explored by the two sets of experimental data (see Sec. IV). Note that the effect of quantum fluctuations on the g (1) function in the quasi-condensate regime are noticeable only for extremely low temperature [START_REF]Quantum fluctuations lead to an algebraic decay of g (1) as (z/ξ) -√ γ/2π for distances ξ < z < 2 /(mkBT ξ), where ξ = / √ mρg is the healing length. For the decay to be noticeable over this range, one requires kBT ∝ ρg exp(-2π/ √ γ)[END_REF], which is experimentally irrelevant.

higher dimensions, no ultraviolet divergence occurs when calculating g (1) (z), so that no energy cut-off is required for this particular quantity. This problem can then be mapped onto a time-independent 2D Schrödinger problem, describing a single particle evolving in imaginary time [18], thus leading to fast and easy-to-implement calculations.

The classical field model is parametrised by a single dimensionless variable ν = γ/γ co = γt 2/3 , provided that lengths are rescaled by the correlation length l c = 2 ρ/(mk B T ) [18,23]. Simple analytic formulas are found in the two limits of the ideal Bose gas regime (ν ≫ 1) and the quasi-condensate regime (ν ≪ 1): the momentum distribution is Lorentzian in both limits, with its full width at half maximum (FWHM) being ∆p = 2 /l c in the ideal gas limit and ∆p = /l c in the quasi-condensate limit. One recovers here the expressions derived from the highly degenerate ideal Bose gas model and the high-temperature Bogoliubov theory [17] respectively. Results of numerical calculations spanning the crossover between those two regimes are reported in Fig. 2, where one assumes a fixed value t = 1000 and parametrises the system with γ instead of ν. We define the crossover region as the domain for which ∆p (shown in Fig. 2(a) as the dashed line), differs by more than 10% from both asymptotic limits ( /l c and 2 /l c ), and we find that the crossover extends over about one order of magnitude in ν, or equivalently in γ at fixed t. Looking at the shape of the momentum distribution, we recover the Lorentzian distribution in both asymptotic regimes and, for any value of ν, we find momentum tails decreasing as 1/p 2 , as expected from the analog single-particle problem [24].

C. Validity of the classical field theory at the quasi-condensation crossover

The classical field theory is valid only if the population of the relevant modes is large. Within the ideal Bose gas regime, this requires ρ ≫ ρ d = √ mk B T / , i.e. a highly degenerate gas. This condition writes equivalently γ ≪ γ d where γ d = t -1/2 is the interaction parameter at degeneracy. Thus, the classical field theory only correctly describes the quasi-condensation crossover (which occurs at γ ≃ γ co ) provided that γ co ≪ γ d . The last condition translates into t 1/6 ≫ 1. Since crossovers span typically about one order of magnitude in γ, one requires that t 10 6 . This value is very difficult to achieve experimentally on cold-atoms experiments while maintaining a temperature sufficiently low to ensure the 1D condition, k B T ≪ ω ⊥ , unless extremely weak atomic interactions are reached, using for instance a Feshbach resonance [25].

In the experiment presented here, we have t < 1000. In this case, according to the above argument, the classical field approach is expected to be inaccurate at the quasicondensation crossover. This can be seen in Fig. 2(a), which shows the FWHM of the momentum distribution as a function of γ according to the ideal Bose gas theory and the classical field prediction for t = 1000. The ideal Bose gas theory is parametrised by χ = γ/γ d ≃ ρ d /ρ, which quantifies the level of degeneracy: the FWHM of the momentum distribution within this theory goes from the Maxwell-Boltzmann prediction for χ ≫ 1 to 2 /l c for χ ≪ 1, and this crossover spans more than one order of magnitude in χ, or equivalently in γ at fixed t. The classical field theory correctly describes the quasi-condensation crossover only if it shares a common plateau with the ideal Bose gas theory at ∆p = 2 /l c in the degenerate ideal Bose gas regime, where γ co ≪ γ ≪ γ d . As seen in Fig. 2(a), this is not the case at all for t = 1000: the highly degenerate ideal Bose gas regime is not very well identified for such a small value of t as far as the momentum distribution is concerned, so that the classical field theory does not correctly describe the quasi-condensation crossover.

D. Quantum Monte Carlo approach

For t < 10 6 , one needs therefore a more accurate treatment than the above-mentioned classical field approximation. We make use of QMC simulations to provide numerically exact results for the Lieb-Liniger model at finite temperature. Bosons in continuous space can be efficiently simulated making use of the density matrix renormalization group in one dimension [26] and Monte Carlo techniques in arbitrary dimensions. Variational and diffusion Monte Carlo results for the first-order correlation function at T = 0 are reported in Ref. [27]. At finite temperature, one can in principle make use of path-integral Monte Carlo simulations [28,29], although we are not aware of explicit calculations for the Lieb-Liniger model. This latter technique requires the discretization of the imaginary-time dimension (Trotter approximation), and its influence on the simulation results has to be carefully removed. Here, we adopt a rather complementary approach, discretizing the spatial dimension instead [26]. Indeed, bosonic lattice models lend themselves to efficient QMC simulations free of any Trotter approximation. By discretizing the field operator ψ( †) (z) ≈ b( †) j / √ a, where a is the lattice spacing, bj , b † j are bosonic operators at site j, and z ≈ ja, we obtain the Bose-Hubbard Hamiltonian

ĤBH = j -J b † j+1 bj + h.c. + U 2 nj (n j -1) -(μ -v j ) nj , (2)
with nj = b † j bj , and parameters J = 2 /(2ma 2 ), U = g/a, μ = µ -2J. We also include the presence of an external potential V (z), discretized to give v j . The Lieb-Liniger parameters γ and t read

γ = U 2Jn , and t = 4k B T J J U 2 , (3)
where n = nj = ρa is the lattice filling.

We study the Bose-Hubbard approximation to the Lieb-Liniger Hamiltonian by making use of Stochastic Series Expansion (SSE) QMC [30,[START_REF] Sandvik | AIP Conf.Proc[END_REF], extensively used to investigate lattice bosons. The Lieb-Liniger model is correctly recovered when the lattice spacing is much smaller than the correlation length in the system. In the classical ideal Bose gas regime, this condition amounts to require λ dB /a = 4πJ/(k B T ) ≫ 1. On the other hand, in the degenerate ideal Bose gas regime and quasi-condensate regime, we need l c /a ≫ 1, implying 2J/(U γt) ≫ 1 and 4J/(U γt) ≫ 1 respectively. The first-order correlation function and the momentum distribution are efficiently sampled with SSE-QMC during the directed-loop update [30,[START_REF] Sandvik | AIP Conf.Proc[END_REF]. The lattice simulation reproduces faithfully the momentum distribution of the Lieb-Liniger model for momentum p ≪ π/a. Fig. 2 shows QMC calculations for a homogeneous gas at t = 1000, for values of γ that span the quasicondensation crossover. In Fig. 2(a) the FWHM of the momentum distribution calculated with QMC is compared to that of the classical field and of the ideal Bose gas models. The QMC results follow the ideal Bose gas prediction almost until the FWHM of the ideal Bose gas prediction crosses that of the classical field prediction at γ = γ c , and then converges towards the classical field prediction. The disagreement between QMC and the above theories never exceeds 20%. These results mitigate the "failure" of the classical field approximation: as long as the FWHM is concerned, an approximate model where one uses the ideal Bose gas prediction for γ > γ c and the classical field prediction for γ < γ c , would give the correct predictions within 20% error. Investigation of the full momentum distribution is shown in Fig. 2 FIG. 2. Momentum distribution of a homogeneous gas at a t = 1000. (a) FWHM according to the classical field theory (dashed line), ideal Bose gas theory (solid line), Maxwell-Boltzmann prediction (dash-dotted line) and QMC calculations (stars). The arrows show the interaction parameter at the quasi-condensation crossover (γco = t -2/3 ) and at the degeneracy crossover (γ d = t -1/2 ). (b) momentum distribution in log scale for three different value of γ. Shown are QMC results (solid lines), classical field predictions (dashed lines) and ideal Bose gas prediction (dotted lines). Error bars, shown for a single point on the side of the profile, corresponds to one standart deviation of the QMC calculation.

three different values of γ. On the ideal Bose gas side of the quasi-condensation crossover, for γ = 5γ co , the QMC momentum distribution follows closely the ideal Bose gas behaviour. The classical field result is recovered for γ ≈ γ co /3, while the momentum width is still about 10% above the Bogoliubov prediction. For γ = γ co , we find that the tails of the momentum distribution agree well with the ideal Bose gas prediction, while the width of the distribution is narrowed by about 20%. Finally, note that, although the highly degenerate ideal Bose gas prediction ∆p = 2 /l c has not been reached before the gas undergoes the quasi-condensation crossover, the momentum distribution at γ ≃ γ co is much narrower than the Maxwell-Boltzmann prediction (dash-dotted line in Fig. 2(a)), so that the effect of degeneracy within the ideal Bose gas regime is already substantial.

E. Local density approximation

While the previous results hold for homogeneous infinite systems, in the experiment described below, the gas is trapped longitudinally in a harmonic confinement of frequency ω/(2π). As long as the correlation length l c is much smaller than the cloud size L, however, a local density approximation (LDA) is valid: at a given position z, the local value of an observable O is that of a homogeneous gas (O hom ) at the chemical potential µ loc (z) = µ-mω 2 z 2 /2, and the global value for a trapped system (O trap ) can be obtained by adding the contributions of each slice of the gas, i.e.,

O trap (ω, µ) = 2 mω 2 ∞ 0 dμ √ μ O hom (µ -μ) . (4)
Therefore, within this approximation and for a given value of µ (fixing the peak density), the normalized momentum distribution is independent of ω, while the total atom number scales as 1/ω. For a gas whose peak linear density lies in the crossover from the degenerate ideal Bose gas regime to the quasi-condensate regime, the LDA condition l c ≪ L writes ω ≪ mg 2 k 2 B T 2 / 5 1/3 [25]. At much higher densities, when almost the whole cloud lies within the quasi-condensate regime, the LDA condition gives ω ≪ k B T 2mg/ρ/ 2 [8]. Both conditions are fulfilled in the data shown below, so that the LDA is expected to be valid. Using QMC, we compute the momentum distribution (normalised to the total atom number) of a harmonically trapped gas for different chemical potentials (i.e. different peak linear densities) and temperatures. To shorten the computational time, calculations are performed for trapping frequencies five times larger than the experimental values, while still satisfying the LDA condition, as verified numerically. The computed momentum distribution are multiplied by five to be compared to experimental data. Linear interpolation between the calculated distributions permits the calculation of the momentum distribution for any chemical potential µ and temperature T , and enables us to perform fits to the measured distribution. We check a posteriori that the interpolation is correct up to a few percent.

III. EXPERIMENTAL TECHNIQUES

One-dimensional atomic clouds are realised in our experiment using an atom-chip setup [START_REF] Reichel | Atom Chips[END_REF]. More specifically, we load 8 × 10 4 Rb 87 atoms in the modulated guide described in [3]. The longitudinal harmonic potential has an oscillation frequency of about 8 Hz and is obtained by DC currents running through wires I 3 and I 4 (see Fig. 3). The transverse potential is provided by three wires carrying a 200 kHz sinusoidal current. The transverse oscillation frequency can be varied from a few hundred Hz to tens of kHz by changing the current amplitude. An important feature of this design is the independence between the longitudinal and transverse confinements. After 1200 ms of radio-frequency evaporation, we let the cloud reach thermal equilibrium for 400 ms, after which no remaining breathing can be observed. To characterise the atomic cloud, we perform in situ density fluctuations measurements, analysing the noise in an ensemble of absorption images as described in [5].

To measure the momentum distribution, one can monitor the free expansion of the gas after all confining potentials have been turned off. The rapid transverse expansion amounts to an effective instantaneous suppression of the interactions, so that the longitudinal expansion reflects faithfully the initial momentum distribution in the trap. However, reaching the far field regime where the longitudinal density profile is homothetic to the momentum distribution requires unrealistically large field of view and long expansion time for our experimental parameters. We thus use the so called focusing technique, already applied to 1D [11,[START_REF] Shvarchuck | [END_REF]34] and 2D systems [35]. Its principle is described below. A strong longitudinal harmonic potential is applied for a very short time during which atoms do not have time to move but acquire a longitudinal momentum shift δp = -Az, proportional to their distance z to the trap center. Then, the atomic interactions are suppressed by switching off the transverse confinement, and the longitudinal confinement is removed so that the cloud starts a free evolution. After a focusing time τ = m/A, the longitudinal density distribution f (z) is homothetic to the initial momentum distribution n(p) = (m/τ )f (pτ /m) [36]. More precisely, the time sequence we use is drawn in Fig. 3. The longitudinal focusing potential is realized by a four-wire configuration that cancels out the non harmonic terms up to the 6th order, thus minimizing aberrations. It has an oscillation frequency ω kick /2π = 40 Hz and it is applied for about t kick = 0.6 ms. After the focusing pulse, we ramp the transverse frequency down to 700 Hz in 1 ms before completely turning off the transverse trapping potentials and letting the cloud evolve freely. The 1 ms transverse ramp is used to limit the transverse velocity spread, so that the transverse final expansion of the cloud is reduced and the signal to noise ratio (SNR) is increased. We verified that this ramp is however quick enough to leave the longitudinal velocity distribution unchanged. After free evolution during the focusing time τ ≃ 1/ω 2 kick t kick ≈ 27 ms, we take an absorption image of the cloud and extract the longitudinal density distribution, from which we deduce the initial momentum distribution. The focusing time τ is adjusted by minimising the width of the density distribution of cold samples. Moreover, the width of very cold clouds saturates at a value that gives us the momentum resolution ∆p/ = 0.2 µm -1 . This resolution is significantly smaller than the width of the data in Fig. 5.

On the other hand, a better resolution may be necessary for a momentum distribution that is much narrower. This can be achieved by a longer focusing time τ . However, the depth of field of our optical system limits the maximum free fall under gravity to about 3 mm, corresponding to a free fall time of about 25 ms. Nevertheless, larger focusing times can be obtained if a transverse potential is maintained during the focusing time, strong enough to overcome gravity (which acts in the transverse plane in our experiment), but weak enough so that atomic interactions are negligible. More precisely the guiding scheme we use is the following. After having performed the longitudinal focusing pulse and switched off the longitudinal potential, we ramp the transverse frequency down in 1 ms, similar to the procedure described before. The final transverse frequency is f l = 700 Hz, which is the lowest value that allows the atoms to remain transversally trapped. We then hold the cloud for a guiding time t g , at the end of which we switch off the transverse confinement and let the cloud fall freely for t f . The guiding does not affect the longitudinal free evolution as long as the effect of interactions stays negligible. We have checked that this is indeed the case for the explored parameter range by measuring the longitudinal density profile at focus f (z) of identical clouds, either without guiding and a focusing time τ = t f = 24.7 ms, or with the above guiding scheme and τ = t g + t f = 43.3 ms, where t g = 20.6 ms and t f = 22.7 ms. The functions f (mz/τ ), shown on Fig. 4 for a cloud of 7000 atoms at T ≈ 95 nK, give indeed the same momentum distribution. The validation of this technique paves the way towards high resolution measurements of momentum distribution of 1D gases.

IV. EXPERIMENTAL RESULTS

Fig. 5 shows experimental results referring to two different parameter sets (particle number, temperature and trapping potential). In both cases, the momentum distribution is obtained without the guiding scheme, with a focusing time of 27 ms.

A. On the quasi-condensate side of the crossover Data in Fig. 5(A) correspond to a cloud initially trapped with a transverse oscillation frequency ω ⊥ /(2π) = 6.4 kHz and a longitudinal frequency of 8.3 Hz, measured using parametric heating and dipole oscillations respectively. Since the cloud temperature and chemical potential are much smaller than the transverse level spacing, the system is well described by a 1D gas with a coupling constant g = 2 ω ⊥ a s , where a s = 5.3 nm is the 3D scattering length [37]. In Fig. 5 (A,12), we show the measured momentum distribution (points) with a QMC fit (solid lines) that yields T = 72 nK, corresponding to t = 76. Note that the shape is reproduced up to a few per cent. We report in Fig. 1 the segment [(γ m , t), (γ M , t)] corresponding to the domain explored by the data: γ m corresponds to the peak density, and the segment length is such that 80% of the atoms lie in a region of density larger than that associated to γ M . We find that γ m ≃ γ co /2, where γ co = t -2/3 is the interaction parameter at the quasi-condensation crossover. Thus the cloud lies quite deeply in the quasi-condensate regime. However, the measured momentum distribution is substantially broader than that predicted using the quasicondensate Lorentzian n(p) ∝ 1/(p 2 + ( /l c ) 2 ) together with the LDA, shown as the dashed line in Fig. 5(A,1). On the other hand, for such a value of γ m , the momentum distribution shows slowly decreasing tails, compatible with a Lorentzian behavior within a few percent in a wide range of momenta. This is seen in Fig. 5(A,2) where the dashed-dotted line shows a 1/p 2 decrease. For our SNR of about 80, data are compatible with a Lorentzian ) is the prediction using the LDA and the Lorentzian momentum distribution in the quasi-condensate limit. In 2), the dasheddotted lines give the 1/p 2 behavior, while the momentum associated with the healing length ξ = / √ mρg indicates that sample (A) is more in the quasi-condensate, with ∆p ≪ /ξ, and sample (B) is more in the ideal Bose gas, with ∆p ≫ /ξ. In addition, the dotted lines indicate the contribution from the transverse excited states, which is insignificant for data (A) and hence neglected in our analysis. 3) in situ density fluctuations, with Yang-Yang prediction in solid lines, evaluated at the above-mentioned temperatures. The gray dashed and dotted lines shows Yang-Yang predictions with 30% temperature deviations. The dash-dotted lines indicate the poissonian fluctuation level. 4) in situ profile, with the QMC profile in dashed lines. The dotted lines again show the contribution from transverse excited states. We also plot in solid lines the density profile fitted with Yang-Yang calculation, which yields (A) T = 111 nK and (B) T = 76 nK.

behavior at large p.

As shown in [3,6], an independent and reliable thermometry can be carried out using in situ density fluctuations, provided the pixel size ∆ is both much larger than the correlation length of density fluctuations, and much smaller than the typical length scale of the variation of the mean density profile : the atom number fluctuations are then given by the thermodynamic relation δN 2 = k B T ∆∂ρ/∂µ, where ρ(µ, T ) is the equation of state of homogeneous 1D gases, known exactly using the Yang-Yang calculation [6], and the derivative is computed for µ such that ρ(µ, T ) = N /∆ is the local density. The above criterion is well fulfilled in our experiment and the measured density fluctuations are shown in Fig. 5(A,3). For linear densities larger than 25 atoms per pixel, fluctuations are almost independent of the density, which confirms that the center part of the gas lies within the quasi-condensate regime. Indeed, in this regime µ ≃ gρ, so that the thermodynamic relation reduces to δN 2 = ∆k B T /g, which does not depend on N . The fluctuations computed using the thermodynamic relation and the Yang-Yang equation of state, for the temperature T = 72 nK (obtained by fitting the momentum distribution with the QMC results) are shown as solid line. The data are in agreement with this prediction. With our SNR, the density fluctuations based thermometry has a precision of about 20% (see gray curves in Fig. 5(A,3)), less precise than the thermometry obtained by fitting QMC calculations to the momentum distribution.

Finally, we compare in Fig. 5(A,4) the measured in situ density profile to the QMC density profile (within LDA) for T = 72 nK (dashed line). We find a good agreement for most of the profile, although the measured data show higher wings. Since the central part of the cloud lies deep in the quasi-condensate regime, a large part of the cloud follows the Thomas-Fermi profile, and all the information on the temperature lies in the small wings. This renders Yang-Yang thermometry based on density profile less precise and extremely sensitive to the shape of the wings. Here, a Yang-Yang fit to the density profile gives T = 110 nK, a value 40% higher than that extracted from the QMC fit, incompatible with the measured momentum distribution or the density fluctuation measurements. This discrepancy, and thus the presence of the inflated wings, may come from the anharmonicity of the potential due to its residual roughness [38]. Alternatively, it may also indicate a lack of perfect thermal equilibrium.

B. On the ideal Bose gas side of the crossover While the previous results probe mainly the quasicondensate regime, we also probe the ideal Bose gas side of the quasi-condensation crossover, i.e. data with γ m > γ co = t -2/3 . The ideal Bose gas regime shows a very different behaviour from the trivial Maxwell-Boltzmann prediction only for a large t parameter, for which the quasi-condensation crossover occurs for an already highly degenerate gas. If one wishes to preserve the 1D condition k B T ≪ ω ⊥ , large t parameters can be accessed only by decreasing the transverse confinement [25]. We thus reduced the transverse confinement to 2.1 kHz. Data are shown in Fig. 5(B). No saturation of the density fluctuations is seen on the in situ fluctuation measurements, which indicates that the gas does not lie in the quasi-condensate regime. Fluctuations however rise well above the poissonian level (shown as dashed-dotted line on Fig. 5(B,3)) so that the gas is highly degenerate. Contrary to the data of Fig. 5(A), we now have k B T /( ω ⊥ ) = 0.8 which is of the order of unity, so that the transverse excited states contribute to the measured fluctuations and momentum distribution. We take into account the population of transverse excited states, assuming that they behave as independent ideal 1D Bose gases, while the transverse ground state is treated as a Lieb-Liniger gas. This modified Lieb-Liniger model has been used with success to describe density profiles [4] and density fluctuations [6], and has been applied to predict the rms width of momentum distribution [11].

Here we use our QMC calculations to describe the transverse ground state. A fit using the above strategy reproduces well the measured momentum distribution (see Fig. 5(B,1)), and yields the temperature T = 84 nK, corresponding to t = 840. The segment in the phase space (γ, t) explored by the data is shown in Fig. 1: the peak linear density is close to the quasi-condensation crossover density and most of the cloud lies in the degenerate ideal Bose gas regime. The contribution of the excited states to the momentum distribution is seen in Fig. 5(B,1): it is only 10% in the center but it rises to almost 50% in the wings around |p|/ ≃ 5 µm -1 . Fig. 5(B,2) shows that the theoretical momentum distribution decreases faster than 1/p 2 at large momenta. This is mainly due to the contribution of the excited states, which have approximately Gaussian momentum tails. However, with our SNR of about 50, no deviation from a 1/p 2 behavior can be identified. We report in Fig. 5(B,2) the fluctuations expected for the temperature T = 84 nK (obtained from the fit of the momentum distribution with QMC). They are in agreement with the measured fluctuations. Note however that, for these parameters, the uncertainty of thermometry based on fluctuations is about 30% (see gray lines in Fig. 5(B,3)). The density profile expected for T = 84 nK is not far from the measured one (see Fig. 5(B,4)). In contrast to case (A), Yang-Yang thermometry based on the profile is less sensitive to tiny modifications of the wings, since the profile of the high-density regions is also affected by the temperature. Thus, a thermometry based on the profile is expected to be quite precise and the potential roughness is expected to have a smaller effect. A Yang-Yang fit to the experimental profile yields a temperature of 76 nK, different only by 10% from the temperature deduced from the momentum distribution.

V. CONCLUSIONS AND PROSPECTIVES

We have shown that for t < 10 6 , a simple 1D classical field approximation where the momentum distribution is computed without energy cutoff, fails to describe the momentum distribution of weakly interacting 1D Bose gases in the crossover from ideal Bose gas to quasi-condensate. Thus, we performed QMC calculations to investigate the crossover. Experimentally, we measure the momentum distribution with the focusing technique, which is improved by a guiding scheme for enhanced resolution. We show that the temperature deduced from the momentum distribution is in agreement with an independent fluctuation-based thermometry. We find that the result of Bogoliubov theory is not appropriate for parameters that are close to the ones relevant to this article, even when the cloud lies quite deeply in the quasi-condensate regime. In [11], it has been proposed to use the mean kinetic energy deduced from measured momentum distribution in order to extract the temperature using the modified Yang-Yang calculation. For the data presented in this article, with a SNR of about 50, the measured momentum distribution shows tails compatible with a 1/p 2 behavior both in the quasi-condensate regime and in the degenerate ideal Bose gas regime. Therefore, extraction of the mean kinetic energy from the data is impossible. We argue that this technique was applicable in [11] because of the more 3D nature of the gas: for quasi 1D gases, as discussed in section IV B, the contribution of the transverse excited states leads to fast decaying tails, which enable the extraction of the mean kinetic energy per particle. We also remark that neither our measurements nor the QMC calculations display the asymptotic behaviour of 1/p 4 , which only appears at even larger momenta for our parameters. On the other hand, we believe that this behaviour could be more readily discernable in the strongly interacting regime, where the 1/p 4 tail should contain a larger proportion of the atoms. This work opens many perspectives for the study of 1D gases. First of all, we show that momentum distribution measurements provide a precise thermometry for 1D Bose gases. While in this article we use QMC calculations to fit the temperature of experimental data, a rough estimate of the cloud temperature could be performed with a lower numerical cost. Indeed, the comparison with the exact QMC momentum distribution shows that a model which combines the ideal Bose gas theory and the classical field approximation gives the correct width within 20% precision. Second, momentum distribution measurements are essential to characterize more complex systems. For instance, in the presence of a lattice, it enables the investigation of the correlation properties at the Mott and/or pinning transitions. Contrary to Bragg spectroscopy [10,14], where a single momentum component is probed at each shot, the focusing method as well as time of flight method gives access to the whole distribution at the same time, which allows for noise correlation measurements in momentum space [23,39]. Finally, recording the time evolution of the momentum distribution is essential to monitor the non-equilibrium dynamics and address the issue of thermalization in 1D closed quantum systems. For example, the dynamics resulting from a quench of the 1D coupling constant can be investigated. The measure of the momentum distribution of impurities would also permit the study of impurity dynamics, which is currently under intense investigation [40].

  FIG.1. Phase diagram of the 1D Bose gas in the parameter space (γ, t). The ideal Bose gas to quasi-condensation crossover occurs for γ ≃ γco = t -2/3 (solid line). Within the ideal Bose gas regime (light gray area), the crossover from the classical to the degenerate sub-regime occurs for γ ≃ γ d = t -1/2 (dashed line). The region γ, t ≫ 1 (dark gray area) is the fermionised regime. The segments correspond to the regions explored by the two sets of experimental data (see Sec. IV). Note that the effect of quantum fluctuations on the g(1) function in the quasi-condensate regime are noticeable only for extremely low temperature[START_REF]Quantum fluctuations lead to an algebraic decay of g (1) as (z/ξ) -√ γ/2π for distances ξ < z < 2 /(mkBT ξ), where ξ = / √ mρg is the healing length. For the decay to be noticeable over this range, one requires kBT ∝ ρg exp(-2π/ √ γ)[END_REF], which is experimentally irrelevant.

  (b) for

FIG. 3 .

 3 FIG. 3. (a) Wire configuration: currents I1, I2, I3 and I4 realise the longitudinal trapping potential of frequency fz. The transverse confinement of frequency f ⊥ is ensured by the three wires carrying a current modulated at 200 kHz. (b,c) Focusing sequences. The focusing potential, harmonic to order 5 in z, is applied during t kick before the longitudinal potential is turned off. After a free evolution during the focusing time τ , an absorption image is taken. Without guiding, see (b), the transverse confinement is removed during the whole focusing time, though a 1 ms ramp down to f ⊥ = 700 Hz is carried out to decrease the transverse expansion. In the presence of guiding, see (c), f ⊥ is kept at 700 Hz for a time tg, during which the cloud does not fall. Numerical values in (b) and (c) correspond to data of Fig.5 and 4respectively.

5 FIG. 4 .

 54 FIG.4. Momentum distribution of a gas of about 7000 atoms at thermal equilibrium, in a 7.5 Hz longitudinal trap, and 4.6 kHz transverse trap. The solid distribution is obtained without guiding and with a free fall of 24.7 ms. The circles are the momentum distribution of the same sample, but with a guiding time of 20.6 ms followed by a free fall of 22.7 ms.

FIG. 5 .

 5 FIG.5. Data for atoms confined in a harmonic trap of transverse oscillation frequency of (A) 6.4 kHz and (B) 2.1 kHz, and a longitudinal oscillation frequency of (A) 8.3 Hz and (B) 7.6 Hz. The pixel size, in the object plane, for in situ data is ∆ = 2.7 µm. 1),2) momentum distribution in linear and log scales. QMC fits are shown in solid lines, yielding (A) T = 72 nK and (B) T = 84 nK. The dashed line in (A,1) is the prediction using the LDA and the Lorentzian momentum distribution in the quasi-condensate limit. In 2), the dasheddotted lines give the 1/p 2 behavior, while the momentum associated with the healing length ξ = / √ mρg indicates that

ACKNOWLEDGMENTS

The authors thank Karen Kheruntsyan for providing the data from the Yang-Yang calculations. This work was supported financially by the Triangle de la Physique, by the ANR grant ANR-09-NANO-039-04, the Austro-French FWR-ANR Project I607, the ARC Discovery Project Grant No. DP110101047, and the CoQuS Graduate school of the FWF.


Felipe Aguirre
email: felipe.aguirre-martinez@utc.fr

Mohamed Sallak
email: mohamed.sallak@utc.fr

Walter Schön
email: walter.schon@utc.fr

Fabien Belmonte
email: fabien.belmonte@transport.alstom.com

Application of Evidential Networks in quantitative analysis of railway accidents

Keywords: Railway accidents, Belief functions theory, Evidential Networks, epistemic uncertainty, Experts' opinion

Currently, a high percentage of accidents in railway systems are accounted to human factors. As a consequence, safety engineers try to take into account this factor in risk assessment. However, human reliability data are very difficult to quantify, thus, qualitative methods are often used in railway system's risk assessments. Modeling of human errors through probabilistic approaches has shown some limitation concerning the quantification of qualitative aspects of human factors. The proposed paper presents an original method to account for the human factor by using Evidential Networks and fault tree analysis.

I. INTRODUCTION

There is little doubt that human error is the most significant source of accidents or incidents in railway systems.

According to statistics of railway accidents in Korea, human errors have accounted for 61% of all train accidents from 1995 to 2004. The purpose of Human Reliability Analysis (HRA) is to include the likelihood of human actions that may cause hazardous events occurring during the qualitative or quantitative evaluation of risk. Several new techniques for HRA were invented in the last decades: THERP [START_REF] Swain | Handbook of human-reliability analysis with emphasis on nuclear power plant applications[END_REF], HEART [START_REF] Williams | HEART, a proposed method for assessing and reducing human error[END_REF], JHEDI [START_REF] Kirwan | A resources flexible approach to human reliability assessment for PRA[END_REF], etc. These techniques have been often attacked as being of dubious validity [START_REF] Reason | Human Error[END_REF], [START_REF] Hollnagel | Human reliability analysis: Context and control[END_REF]. One of the main critics is that human performance is not easy to quantify due to a large number of factors affecting it and the variability of one person over others. Nevertheless, studies have been made to prove the performance of current methods [START_REF] Kirwan | The validation of three human reliability quantification techniques THERP, HEART and JHEDI: Part I technique descriptions and validation issues[END_REF], [START_REF] Kirwan | The validation of three Human Reliability Quantification techniques THERP, HEART and JHEDI: Part II Results of validation exercise[END_REF], [START_REF] Kirwan | The validation of three human reliability quantification techniques THERP, HEART and JHEDI: part III Practical aspects of the usage of the techniques[END_REF], [START_REF] Bell | Review of human reliability assessment methods[END_REF], [BFB + 10]. One of their drawbacks is that performance factors are dependent of the working context and a validation of the experimental data must be done. This suggest that there is a need for new methods of HRA [START_REF] Hovden | Is there a need for new theories, models and approaches to occupational accident prevention?[END_REF].

Moreover, in the field of railway systems there has been few attempts to include human factors in safety analysis [START_REF] Belmonte | Interdisciplinary safety analysis of complex socio-technological systems based on the Functional Resonance Accident Model: An application to railway traffic supervision[END_REF]. Some of the factors which cause the occurrence of human error in the railway systems are classed in four categories [START_REF] Knaggs | Human reliability and its relation to accidents by railway and on sea[END_REF]: (1) automatic action-automatism, which is a kind of human habitude, (2) imperative idea, which is a kind of desire to do something without any reason, (3) the condition and capacity of the faculty of observation possessed by the observers, which contain mental alienation, toxic action of drugs or liquors, inefficient training, defect in the special sense like acuteness of hearing, color-blindness, acuteness of vision, etc, (4) the mental condition of the observers and their capacity at the time for the correct interpretation of what is to be observed, which contain overwork, anxiety, depressed conditions, impaired mental condition, etc. A method was developed by Vanderhaegen et al. to analyze the consequences of human unreliability on railway system safety [START_REF] Vanderhaegen | APRECIH: a human unreliability analysis method application to railway system[END_REF], [START_REF] Vanderhaegen | A non-probabilistic prospective and retrospective human reliability analysis method application to railway system[END_REF]. The human behavioral degradations were characterized by a behavioural model of human unreliability which included three behavioral factors: acquisition related factors, problem solving related factors and action related factors. Baysari et al. [START_REF] Baysari | Understanding the human factors contribution to railway accidents and incidents in Australia[END_REF] used the Human factors analysis and classification system (HFACS) method to identify errors in railway and human factors contribution in accidents and incidents in Australia. In this method, four errors types were proposed: unsafe acts, preconditions for unsafe acts, unsafe supervision and organizational influences.

Unsafe acts contain errors (skill based, decision, perception) and violations (routine, exceptional), preconditions for unsafe acts contain condition of operators (adverse mental states, adverse physiological states, physical/mental limitation), personnel factors (crew resource management, personal readiness) and environmental factors (physical environment, technological environment), unsafe supervision contains inadequate supervision, planned inappropriate operations, failed to correct a known problem and supervisory violations, organizational influences contain resource management (human resources, monetary/budget resources, equipment/facility resources), organizational climate (structure, policies, culture) and organizational process (operations, procedures, oversight). From the classification of error types, the causes of human errors can be deduced. In 2009, two human error identification techniques were used by Baysari et al. [START_REF] Baysari | Classification of errors contributing to rail incidents and accidents: A comparison of two human error identification techniques[END_REF] in railway systems: HFACS method and the Technique for the retrospective and predictive analysis of cognitive errors (TRACEr-Rail version). HFACS method identified errors of all operators in railway systems, but TRACE-Rail method identified only errors of drivers. In TRACE-Rail method, factors which cause driver errors were listed as follows: train/infrastructure/traffic, communications, procedures/documentation, information, training/knowledge/experience, workplace design/human-machine interaction/equipment, in-cab environment, personal factors and social and team factors. In [START_REF] Dong San | Development and evaluation of a computer-aided system for analyzing human error in railway operations[END_REF], kim et al. explained that human error analysis was a time-consuming task and a computer-aided system which helps to analysis human error in railway systems was developed. This system supported the hierarchy of error causes and the relations between these causes because of predefined links. In this system, thirteen categories of human error causes that influenced the occurrence of the human error were proposed: mental states of operators, physical states of operators, knowledge/experiences/ability of operators, task characteristics, tools/equipment, work environment, train/infrastructure, rules/procedures, human resource management, communication, team factors, supervision and organizational processes/policies/culture. More details can be added into each category (physical states of operators contain physical fatigue, physical illness, alcohol or drug use, temporary visual or hearing impairment, physical properties, general motor ability and age/gender, etc.). In [START_REF] Belmonte | Interdisciplinary safety analysis of complex socio-technological systems based on the Functional Resonance Accident Model: An application to railway traffic supervision[END_REF], Belmonte et al. presented an application of Functional Resonance Accident Models (FRAM) in safety analysis of railway traffic supervision using modern Automatic Train Supervision (ATS) systems. This study included not only technological, but also human and organizational components.

However, these works are mostly of a qualitative nature mainly because of the difficulty to quantify human behaviour. As stated by R.L. Boring et al. [BFB + 10] one of the lessons learned on benchmarking from the international human reliability analysis empirical study is the fact that taking into account both qualitative and quantitative elements in the prediction of HRA methods allowed a more complete understanding of the HRA methods' strengths and weaknesses in predicting equipments performance. Human behaviour is considered as surrounded by epistemic uncertainties, thus needing the use of proper theories to represent and propagate the uncertainty in risk analysis. As stated by M. Konstandinidoua et al. [START_REF] Konstandinidou | A fuzzy modeling application of CREAM methodology for human reliability analysis[END_REF] the limitations in the analysis of human actions in PRAs are always recognized as a constraint in the application of PRA results. The fundamental limitations are as follows:

• insufficient data.

• methodological limitations related to subjectivity of analysts and expert judgment.

• uncertainty concerning the actual behaviour of people during accident conditions.

During the last decades, the reliability assessment community recognized that the distinction between different types of uncertainties plays an important role in reliability evaluation [START_REF] Apostolakis | The concept of probability in safety assessments of technological systems[END_REF], [START_REF] Aven | On the use of uncertainty importance measures in reliability and risk analysis[END_REF]. Uncertainty is considered of two types: aleatory uncertainty which arises from natural stochasticity and epistemic uncertainty which arises from incompleteness of knowledge or data [OHJ + 04]. The distinction is useful because epistemic uncertainty can be reduced by acquiring knowledge on the studied system.

The classical probabilistic approach was widely used to manage aleatory uncertainties in risk and reliability assessments [START_REF] Aven | Interpretations of alternative uncertainty representations in a reliability and risk analysis context[END_REF]. This approach was based on the definition given by Laplace of the probability of an event as the ratio of the number of cases favorable to it, to the number of all possible cases when all cases are equally possible [START_REF] Laplace | A Philosophical Essay on Probabilities[END_REF]. The frequentist probabilistic approach introduced by Venn [Ven66] which defined the event probability as the limit of its relative frequency in a large number of trials was also widely used to describe aleatory uncertainties.

To describe epistemic uncertainties, De Finetti [START_REF] Finetti | Theory of Probability[END_REF] introduced the subjective probabilities of an event to indicate the degree to which the expert believes it. Kaplan and Garrik [START_REF] Kaplan | On the quantitative definition of risk[END_REF] introduced the concept of probability of frequency to expand their definition of risk. Pate-Cornell [START_REF] Cornell | Uncertainties in risk analysis: six levels of treatment[END_REF] used the six level of uncertainty to obtain a family of risk curves in presence of both aleatory and epistemic uncertainties. The Bayesian approach proposed the use of subjective probabilities to represent expert judgement. The probability distributions representing the aleatory uncertainties are first proposed. The epistemic uncertainties about the parameter values of the distributions are then represented by prior subjective probability distributions [START_REF] Kaplan | On the quantitative definition of risk[END_REF]. The equation of Bayes is used to compute the new epistemic uncertainties in terms of the posterior distributions in case of new reliability data. Finally, the predictive distributions of the quantities of interest are derived by using the total probability law. The predictive distributions are subjective but they also take into account the aleatory uncertainties represented by the prior probability models [START_REF] Apostolakis | The concept of probability in safety assessments of technological systems[END_REF].

However, There are some critics about representing epistemic uncertainties using subjective probabilities. Particularly, in the case of components that fail only rarely such as railway systems or components that have not been operated long enough to generate a sufficient quantity of data. This is also the case of human errors. As stated by M. Konstandinidoua et al. [START_REF] Konstandinidoua | A fuzzy modeling application of cream methodology for human reliability analysis[END_REF] the limitations in the analysis of human actions in PRAs are always recognized as a constraint in the application of PRA results. The fundamental limitations are as follows:

• insufficient data.

• methodological limitations related to subjectivity of analysts and expert judgement.

• uncertainty concerning the actual behaviour of people during accident conditions.

In this work we are concerned with the problem of data insufficiency. For example, when there is few information about the value of a parameter α, the choice of probability distribution may not be appropriate. For example, there is a difference between saying that all that is known about the parameter α is that its value is located somewhere in an interval [x, y] and saying that a uniform distribution on [x,y] characterizes degrees of belief with respect to where the value of this parameter is located in the interval [x, y] [HJOS07], [START_REF] Aven | Interpretations of alternative uncertainty representations in a reliability and risk analysis context[END_REF]. Furthermore, in a situation of ignorance a Bayesian approach must equally allocate subjective probabilities over the frame of discernment. Thus there is no distinction between uncertainty and ignorance. A number of alternatives theories based on different notions of uncertainty were proposed to capture the imprecision in subjective probabilities.

Baudrit et al. [START_REF] Baudrit | Joint propagation and exploitation of probabilistic and possibilistic information in risk assessment[END_REF] explained that random variability can be represented by probability distribution functions, imprecision (or partial ignorance) is better accounted for by possibility distributions (or families of probability distributions) and thus propose a hybrid method which jointly propagates probabilistic and possibilistic uncertainty in risk assessment. Tucker et al. [START_REF] Tucker | Probability bounds analysis in environmental risk assessments[END_REF] propose probability bounds analysis which combines probability theory and interval arithmetic to produce probability boxes (p-boxes), structures that allow the comprehensive propagation of both variability and uncertainty through calculations in a rigorous way. The belief functions theory also known as the Dempster-Shafer or evidence theory is a generalization of the Bayesian theory of subjective probability. Whereas the Bayesian theory requires probabilities for each question of interest, belief functions allow us to base degrees of belief for one question on probabilities for a related question [START_REF] Smets | Belief functions: The disjunctive rule of combination and the generalized Bayesian theorem[END_REF]. To illustrate the idea of obtaining degrees of belief for one question from subjective probabilities for another, we propose an example in risk assessment inspired from the example of limb given by Shafer. Suppose we have subjective probabilities for the reliability of a risk expert A. The probability that A is reliable is 0.75, and the probability that A is not reliable is 0.25. The risk expert A reports us that a component i is failed. This information which must be true if A is reliable, is not necessarily false if A is not reliable. The risk expert testimony justifies a 0.75 degree of belief that the component i is failed, but only a 0 degree of belief (not a 0.25 degree of belief) that the component i is not failed. This value does not mean that we are sure that the component i is failed, as a zero probability would. It means that the risk expert's testimony gives us no reason to believe that the component i is failed. The 0.75 and the 0 constitute a belief function. Thus there is no requirement that belief not committed to a given proposition should be committed to its negation. The second point of evidence theory is that belief measures of uncertainty may be assigned to overlapping sets and subsets of hypotheses, events or propositions as well as to individual hypothesis. D-S theory which can be considered as an alternative approach to represent uncertainties has gained an increasing amount of attention both from the theoretical and the applied point of view [START_REF] Guth | A probabilistic foundation for vagueness and imprecision in fault tree analysis[END_REF], [START_REF] Inagaki | Independence between safety-control policy and multiple-sensors schemes via \textsc{D}empster-\textsc{S}hafer theory[END_REF], [START_REF] Simon | Evidential networks for reliability analysis and performance evaluation of systems with imprecise knowledge[END_REF], [START_REF] Sallak | Transferable belief model for reliability analysis of systems with data uncertainties and failure dependencies[END_REF]. In a finite discrete space, D-S theory is a generalization of probability theory where probabilities are assigned to sets instead of mutually exclusive singletons. This theory is still a young field compared to other theories and its main application is data fusion. This paper presents an original method to integrate, in a quantitative way, the human, organizational and technical factors to risk analysis in railway accidents using belief functions theory. We use a graphical model called Evidential Networks (ENs) on which we quantify beliefs given by experts and the relationships between the different factors through a valuation network. The advantage of such method is that it presents a natural and intuitive way for the experts to understand and use the model. To make an analogy, one can say that ENs are to belief functions as Bayesian networks (BNs) are to probability theory. They are both used to represent and propagate our beliefs but the deployment of Bayesian networks is sometimes inadequate for situations involving partial or total ignorance [START_REF] Simon | Bayesian networks inference algorithm to implement Dempster Shafer theory in reliability analysis[END_REF]. Take for example the study done by Contini et al. [START_REF] Contini | Description of the Project, Discussion of the Results and Conclusions[END_REF], they provided the same reference subject (an amonia storage facility) to 11 teams representing control authorities, research organizations, engineering companies and industries. After a complete probabilistic risk assessment was performed, the results differed over several orders of magnitude. This is a direct evidence that eliciting precise probabilities can conduct to abusively precise values as they can differ greatly from one expert to another. On the contrary, as ENs are based on belief functions theory, ignorance or imprecise knowledge can be easily modeled by using lower and upper estimates of the different factors. This way our models can be more conservative and more representative of our real state of knowledge.

More precisely, This paper proposes four contributions:

1) How to construct human reliability data in the belief functions framework (we added this section to the paper) from both observations and experts' opinion with some basic examples.

2) How to combine experts' opinion in the belief functions framework.

3) A general method to apply Evidential networks and FT in the evaluation of belief occurrence of accidents.

4) An application of the proposed approach in the evaluation of belief occurrence of a railway accident and a comparison with results obtained using classical approaches based on Monte-Carlo simulations.

The remainder of the paper is organized as follows. In Section 2, we provide a background of belief functions theory. In Section 3, some basic notions of Valuations Based Systems (VBSs) and ENs are presented. In Section 4, we describe the development of the combined EN and Fault Tree (FT) model. Section 4 demonstrates the proposed methodology on the example of a railway accident. Section 5 concludes the paper.

II. BELIEF FUNCTIONS THEORY

In the 1960's, Dempster [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF] gave birth to the DS theory also called evidence theory or belief functions theory with the study of upper and lower probabilities. Later on, it was extended by Shafer in 1976 [START_REF] Shafer | A Mathematical Theory of Evidence[END_REF].

Belief functions theory can be seen as an extension of the Bayesian probability theory. The first work using belief functions theory in reliability and risk assessment was presented by Dempster and Kong [START_REF] Dempster | Uncertain evidence and artificial analysis[END_REF]. In recent years, the belief functions theory was used by many researchers in order to quantify the uncertainty in reliability and risk

assessment studies [Alm95], [BGC04], [LS07], [UG08], [PS08], [SW09], [SSA10].

A. Representing information

The definition domain of the variable of interest x is called the frame of discernment Ω where all of the possible events are mutually exclusive elementary propositions. It is equivalent to the sample space in probability theory. A Basic Belief Assignment (BBA) on Ω, also called Basic Probability Assignment (BPA), is a function, m Ω : 2 Ω → [0, 1], which maps probability masses on events or subsets of events such that:

A∈2 Ω m Ω (A) = 1 (1)
An agent holding a piece of evidence allocates unitary amounts of beliefs to the different subsets of Ω. The number m Ω (A) represents the support of A given by the agent's belief [START_REF] Smets | The combination of evidence in the Transferable Belief Model[END_REF]. There is a distinction between probabilities and BBAs: probability distribution functions are defined on Ω and BBAs on the power set 2 Ω . Moreover the additivity hypothesis is not required in belief functions theory as it is in probability theory. The subsets A ⊂ Ω such that m Ω (A) > 0 are called focal sets of m Ω . Full knowledge is represented by a BBA having a singleton {x} (x ∈ Ω) as a unique focal set. A Bayesian BBA is a special case where all of the focal sets are singletons and is equivalent to probabilities. Complete ignorance is represented by a BBA having only one focal element equal to Ω and is called vacuous. Finally BBAs have some other properties, that distinguishes them from probability functions

[GF99]:

• It is not necessary that m(Ω) = 1. • It is not necessary that m(A) ≤ m(B) when A ⊆ B.
• m(A) + m( Ā) is not always equal to 1.

B. Belief and plausibility functions

A belief function Bel on Ω is a function Bel : 2 Ω → [0, 1] defined as follows:

Bel(A) = B⊆A m Ω (B), A ⊆ Ω (2)
As we impose the closed world hypothesis that states that the answer is strictly found inside our frame of discernment, we add two constraints to the belief function: Bel(Ω) = 1 and Bel(∅) = 0. The inverse formula called the Möbius transform of Bel is defined as follows:

m(A) = B⊆A (-1) |A|-|B| Bel(B) (3)
Where |A| represents the cardinality of A. A plausibility function P l on Ω is a function P l : 2 Ω → [0, 1] defined as follows:

P l(A) = B∩A =∅ m Ω (B), A ⊆ Ω (4)
Bel(A) is obtained by summing the BBAs of the propositions that totally agree with A whereas P l is obtained by summing BBAs of propositions that agree with A totally or partially (cf. Figure 1). [Bel(A), P l(A)] can be viewed as the interval that describes the uncertainty of A. Bel and P l are related by:

P l(A) + Bel( Ā) = 1 P l( Ā) + Bel(A) = 1 (5)
Where Ā is the negation of the event A.

C. Combination operations

The purpose of information aggregation related to event occurrence is to summarize data whether the data is coming from a single source or multiple sources. There are many aggregation techniques such as arithmetic averages, geometric averages, harmonic averages, maximum and minimum values. Combination rules are the types of aggregation techniques for data obtained from multiple sources which provide different assessments for the same event.

Consider two distinct and independent pieces of evidence m Ω i and m Ω j from two different and reliable sources i and j. In belief functions theory, the principal combination rules are the conjunctive and disjunctive combination rules [START_REF] Shafer | A Mathematical Theory of Evidence[END_REF]. The conjunctive combination is given by:

m Ω i∩j (H) = A∩B=H,∀A,B⊆Ω m Ω i (A)m Ω j (B), ∀H ⊆ Ω (6)
The disjunctive combination is given by:

m Ω i∪j (H) = A∪B=H,∀A,B⊆Ω m Ω i (A)m Ω j (B), ∀H ⊆ Ω (7)
The Dempster's rule of combination is the fundamental rule in belief functions theory for combining items of evidence. This rule is defined as the conjunctive combination of two BBAs followed by a normalization:

m Ω i⊕j (H) =    m Ω i∩j (H) 1-k if A = ∅ 0 otherwise (8)
With:

k = A∩B=∅,∀A,B⊆Ω m Ω i (A)m Ω j (B) (9)
k is called the conflict factor between the two pieces of evidence i and j.

D. Conditioning

The conditioned mass function is intended to represent the impact of additional information. The conditioning in belief functions theory is used both for prediction from observations and revision of uncertain information [START_REF] Eirik Albrechtsen | Conditioning in Dempster-Shafer theory: prediction vs. revision, Belief functions: theory and applications[END_REF].

When dealing with prediction, we have at our disposal a model of the world (e.g., belief occurrence degrees of principals causes of a human error) under the form of a BBAs distribution. Moreover we get some new observations B on the current state of the world (e.g., occurrence of one cause of human error). Then, one tries to predict some property A of the current world with its associated degree of belief (e.g. predict the human error occurrence). The conditional BBA m(A|B) (the belief of observation of A in context B) is used for estimating the degree of belief that the current world satisfies A. Conditioning in belief functions theory requires a sort of independence. Define the mass function m, the conditional of m given B for A ⊆ Ω is given by

m(A|B) = K C|C∩B=A m(C) ( 10
)
where K is the normalisation factor defined by

K -1 = C|C∩B =∅ m(C).
For example m(e i |θ i ) = 0.2 means that when an expert receives a new information θ i , his degree of belief that "an human error e i will occur" is 0.2.

III. BASICS OF VALUATION-BASED SYSTEMS (VBS) AND EVIDENTIAL NETWORKS (ENS

The VBSs were introduced by Shenoy in 1989 as general frameworks for managing uncertainty in expert systems [START_REF] Shenoy | A valuation-based language for expert systems[END_REF]. In rule-based languages, a domain knowledge is represented by variables and rules. The operator used to make inferences from the knowledge is modus ponens which can be be summarized as: p implies q; p is asserted to be true, so therefore q must be true. A categorical rule is a rule defined without uncertainty. In Valuation-Based Languages (VBL), a domain knowledge is represented by variables and valuations. The operators used to make inferences from the knowledge are called combination and marginalization. Thus the VBL have two parts: a static part that is concerned with representation of knowledge, and a dynamic part that is concerned with reasoning with knowledge. For more details, see [START_REF] Shenoy | A valuation-based language for expert systems[END_REF].

A. Valuation based language (VBL) Static part

Let X denote a variable. The set of its possible values will be denoted by Ω X and referred to as the frame of X.

If h is a non empty set of variables then Ω h denotes the Cartesian product of Ω X for X ∈ h:

Ω h = ×{Ω X |X ∈ h}.
We call Ω h the frame for h. The elements of Ω h are called configurations of h.

A valuation V h on h represents some knowledge about the variables in h. The set V represents the set of all valuations, i.e: V E = ∪{V h |h ⊆ E} where E is a finite set of variables.

Dynamic part

A crucial point in VBL is the use of local computations to optimize the propagation of beliefs. In order to compute marginals on the joint valuation, the local computations use two operations called combination and marginalisation.

The combination is a mapping ⊗ : V × V → V used to combine valuations G and H on g and h to produce

valuation G ⊗ H on g ∪ h.
The marginalization is a mapping ↓ h : ∪{V g |h ⊆ g} → V h used to produce valuations on subsets of the given set of variables such that if G is a valuation on g and h ⊆ g, then G ↓h is a valuation on h. G ↓h is called the marginal of G to h.

A VBL makes inference by first combining all valuations in the system in order to obtain the joint valuation. Then the marginalization operation is used to obtain the marginal of the variable of interest.

B. Evidential Systems (ES)

A VBS is a formal mathematical system for representation of and reasoning with knowledge. It consists of a 5-tuple {E, Ω E , V E , ⊗, ↓} where E denote a finite set of variables E, Ω E denote the set of frames, and V E the set of all valuations.

The graphical representation of a VBS is a graph called a valuation network where the nodes represent either valuations or variables. When the valuations are expressed using BBAs, the VBS is called an Evidential System (ES).

A categorical rule x 1 ∧x 2 ∧...∧x n → y is represented by a categorical BBA on the frame Ω X1 ×Ω X2 ×...×Ω Xn ×Ω Y .

Consider two variables X and Y and let Ω X = {x, x} and Ω Y = {y, ȳ} denote their frames. The categorical rule

x → y is then defined as follows [START_REF] Benavoli | Modelling Uncertain Implication Rules in Evidence Theory[END_REF]:

m ΩX ×ΩY ({(x, y), (x, y), (x, ȳ)}) = 1 (11)
On the other hand, if we hold some doubt about the rule x → y, our degree of belief can be quantified as follows:

m ΩX ×ΩY ({(x, y), (x, y), (x, ȳ)}) = p (12) m ΩX ×ΩY ({Ω X × Ω Y }) = 1 -p
Consider for example a study to evaluate the reliability of a risk expert E. We know that at least 80% of risk experts which are working in factory A are not reliable and at least 90% of risk experts which are working in We suppose that the variables A = a, B = b and R = r represent, respectively, the propositions "the risk expert E works in factory A", "the risk expert E works in factory B", "the risk expert E is reliable". The BBAs representing the fact that at least 80% of risk experts which working in factory A are not reliable are:

m ΩA×ΩR 1 ({(a, r), (ā, r), (ā, r)}) = 0.8 (13) m ΩA×ΩR 1 ({Ω A × Ω R }) = 0.2
The BBAs representing the fact that at least 90% of risk experts which working in factory B are reliable are:

m ΩB×ΩB 2 ({(b, r), ( b, r), ( b, r)}) = 0.9 (14) m ΩB×ΩB 2 ({Ω B × Ω R }) = 0.1
The obtained EN is shown in Figure 2. First we suppose that we have no prior BBAs about the fact that the expert E works in factory A or B, i.e:

m ΩA 0 ({a, ā}) = 1 (15) m ΩB 0 ({b, b}) = 1
The preliminary results obtained using Dempster rule of combination, without receiving any information, is that the belief interval that the expert E is reliable is [0,1] (cf. Table I). Then, we receive the information 1 that the expert E works in factory A. This information is represented by m ΩA 0 ({a}) = 1. The obtained belief interval that the expert is not reliable is [0.8,1]. Finally we received the second information that the expert E works also in factory B which is represented by m ΩB 0 ({b}) = 1. As we can see, this information weakens the fact that the expert E is not reliable which justify the fact that the obtained belief interval that the expert is not reliable is [0.2857,0.3571] (cf. Table I).

m In railway risk assessment, there are two different sources of information available about the human errors

Ω A 0 m Ω A,R 1 A m Ω B 0 m Ω B,R 2 B R Ω A = {a, a} Ω B = {b, b} Ω R = {r, r} Information 1 m Ω A ({a}) = 1 Information 2 m Ω B ({b}) = 1
• Expert opinion in the form of values or bounds on human errors in presence of some events.

• Data in form of number of human errors in presence of some events or time of observations.

A. Lower and Upper Expectations

We first need to introduce some basic concepts. Let B(R) be a Borel σ-algebra, and let m be a Basic Belief Assignment (BBA) and Bel and P l the associated belief and plausibility functions. Let Π be the set of probability measures compatibles with Bel and Pl. A probability measure P ∈ Π is said to be compatible with Bel and Pl if for each A ∈ B(R)

Bel(A) ≤ P (A) ≤ P l(A) (16)
The lower and upper expectations of a measurable, bounded and real-valued function f with respect to m are defined as follows [Was90]

E m (f ) = inf P ∈Π E p (f ) and E m (f ) = sup P ∈Π E p (f ) (17)
where E p (f ) = f (α)P (dα) denotes the expectation of function f with respect to P .

Lower expectations can be obtained from upper expectations through the expression

E m (f ) = -E m (-f ) (18)
It can be shown that in the continuous case [Was90]

E m (f ) = +∞ -∞ +∞ u f (u, v) inf u≤x≤v f (x)dvdu and E m (f ) = +∞ -∞ +∞ u f (u, v) sup u≤x≤v f (x)dvdu (19)
In the discrete case, we have [Dem67]

E m (f ) = n i=1 m i inf x∈Ii f (x) E m (f ) = n i=1 m i sup x∈Ii f (x) (20)

B. Construction of human errors data based on experts opinions

Let X be a indicator variable for the occurrence of a human error in the presence on an event. The occurrence of a human error is a Bernoulli process with a parameter p correspond to the occurrence of error and 1-p corresponds to the absence of error. Thus we have

Bel(X = 1) = E m (p) (21) = n i=1 m(I i ) min p∈Ii p = n i=1 m i u i P l(X = 1) = E m (p) (22) = n i=1 m(I i ) max p∈Ii p = n i=1 m i v i
Where

I i = [u i , v i ] and m(I i ) = m i .
The obtained BBAs are

m({1}) = n i=1 m i u i (23) m({0}) = n i=1 m i v i m({0, 1}) = 1 - n i=1 m i v i
Example 1: Suppose an expert provides his opinion that in the presence of noisy environment there is 0.3 a chance that the parameter p indicating the occurrence of human error lies between and 0.003 and 0.005, and a there is 0.5 chance that p lies between and 0.007 and 0.01. This suggests a mass function over Thus the occurrence of human error p in presence of noisy environment is given by: p ∈ [0.0044, 0.2065].

C. Construction of human errors data based based on observations

Let us now assume that, instead of eliciting expert opinion, we have made n independent observations X 1 , ..., X n of X. Given data of x human errors in n presence of noisy environment, what is the probability of human error in presence of noisy environment. The formulas for the upper and lower expectation of p are defined as follows:

E m [p] = x n + 1 (24) E m [p] = x + 1 n + 1 (25)
Example 2: Suppose in 15 observations, we have human errors in 3 observations. Then we have

E m [p 1 ] = 3 16 = 0.1875 (26) E m [p 1 ] = 4 16 = 0.2500 (27)
The occurrence of human error p 1 in presence of noisy environment is given by: p 1 ∈ [0.1875, 0.2500].

In the case we use the observations made in an other railway system (prior belief): in 35 observations, we have human errors in 5 observations

E m ′ [p 2 ] = 5 35 = 0.1429 (28)
E m ′ [p 2 ] = 6 35 = 0.1714 (29)
In this case, the occurrence of human error p 2 in presence of noisy environment is given by: p 2 ∈ [0.1429, 0.1714].

Finally, the aggregation of the two observations give: p 12 ∈ [0.1429, 0.2500].

D. Comparison with Laplace, LME (Maximum Likelihood Estimate) and Bayesian approaches

As shown in

A. Introduction

Within the context of this paper, we use the definition of human errors given by Swain and Guttman [START_REF] Swain | Handbook of human-reliability analysis with emphasis on nuclear power plant applications[END_REF]: An error is an out of tolerance action, where the limits of tolerable performance are defined by the system. The human error is then considered as a deviation from expected performance (signal passed at danger (SPAD) incidents, over speeding, etc.) and this deviation is defined by the consequence (immediate dangerous situation or accident).

Due to its capability in modeling the systems, ENs have been combined with FT analysis method to determine the belief occurrence of the railway accident. The proposed approach consists of five steps: 1) Constructing a Cause Effect diagram. At this stage we identify the principal causes and their sub-causes.

2) Constructing ENs for only principal causes which have sub-causes. At this stage the relations between the principal causes and their sub-causes found at the previous step are translated into ENs.

3) Belief propagation (marginalization and combination operations). At this stage we propagate the prior BBAs of sub-causes and the conditional BBAs representing the relations between principal causes and their sub-causes through the network and compute the marginal BBAs of principal causes. 4) Constructing a FT which has top event (e.g. the railway accident) and basic events represented by principal causes. 5) Using generalized expressions of belief and plausibility measures to evaluate the top event belief occurrence from the marginal BBAs of principal causes computed using belief propagation.

B. Construction of Cause Effect diagram

In order to identify, sort, and display possible causes of the railway accident we use the Cause Effect diagram method which graphically illustrates the relationship between a given outcome and all the factors that influence the outcome. This type of diagram is sometimes called an "Ishikawa diagram" because it was invented by Kaoru For example consider a railway accident where the driver passes the signal at danger. The driver have mistakenly read the signal aspect for an adjacent line. The investigations concluded that the driver read an adjacent signal that was displaying a single yellow cautionary aspect due to three principal causes: driver route knowledge, infrastructure factors and driver alertness. Three sub-causes dr 1 , dr 2 and dr 3 were associated with the driver route knowledge, two sub-causes i 1 and i 2 were associated with the infrastructure factors, and three sub-causes da 1 , da 2 and da 3 were associated with the driver alertness. The Cause Effect diagram of this example is presented in Figure 3.

C. Construction of Evidential Network

The relations between each principal cause and its sub Figure 4). In this paper, due to insufficiency data, we use expert opinions to define the conditional and prior BBAs.

Then we use marginalization and Dempster combination to obtain BBAs of principal causes.

D. Combined FT with ENs

A FT represents how combinations of basic events lead to the occurrence of the undesired event (top event).

In this work, the principal causes are the basic events (cf. Figure 5) of the FT and we aim to compute the belief occurrence of top event using the BBAs of principal causes computed using the ENs.

Based on our previous works proposed in [START_REF] Sallak | Transferable belief model for reliability analysis of systems with data uncertainties and failure dependencies[END_REF], [START_REF] Aguirre | Generalized expressions of reliability of series-parallel and parallel-series systems using the Transferable Belief Model[END_REF], we can compute the interval belief of top event occurrence using the BBAs of basic events and the minimal cut sets of the FT. Let us consider that for each basic event e i , the state of belief on its occurrence is bounded by [Bel({e i }), P l({e i })] defined over Ω ei = {e i , e i } where e i and e i denote respectively the occurrence and the absence of event e i . The belief occurrence of the undesired event is obtained using Eqs. 30 and the following notation:

N C Number of minimal cuts in the FT.
C i Index set of the i th minimal cut set. Where size(C i ) denotes the number of basic events in the minimal cut set C i . Therefore, the belief occurrence of the top undesired event is bounded by the interval [Bel({e top }), P l({e top })]. In the previous example when the driver passes the signal at danger, we have three minimal cut sets:

Bel({e

top }) = NC i=1   1 - size(Ci) j=1 1 -m Ωe C i (j) (e Ci(j) )   P l({e top }) = NC i=1   1 - size(Ci) j=1 m Ωe C i (j) (e Ci(j) )   (30)
• C 1 : Driver route knowledge.

• C 2 : Infrastructure factors.

• C 3 : Driver alertness.

The belief and plausibility measures are obtained as follows:

Bel({e top }) = m Ωe C 1 ({e C1 })m Ωe C 2 ({e C2 })m Ωe C 3 ({e C3 }) P l({e top }) = (1 -m Ωe C 1 ({e C1 }))(1 -m Ωe C 2 ({e C2 )})(1 -m Ωe C 3 ({e C3 })
In some cases, risk engineers need to convert the interval beliefs of top event occurrences to a probability measures.

Such a transformation is called a probabilistic transformation. We define the probabilistic transformation as a mapping f : m → P where P denotes the probability distribution and m the BBA function [START_REF] Daniel | Probabilistic Transformations of Belief Functions[END_REF]. A probabilistic transformation f is:

• ulb-consistent (upper and lower bound consistent): if Bel(A) ≤ f (A) ≤ P l(A) for any A ⊆ Ω.

• p-consistent (or probabilistically consistent): if f (m) = m for any Bayesian BBA m.

• α-consistent: if f (αm 1 + (1 -α)m 2 ) = αf (m 1 ) + (1 -α)f (m 2 )
for any BBAs m 1 and m 2 .

The most known probabilistic transformation is the pignistic transformation BetP . It was introduced by Smets

and Kennes [START_REF] Smets | TheTransferable Belief Model[END_REF] and corresponds to the generalized insufficient reason principle: a BBA assigned to the union of n atomic sets is split equally among these n sets. The pignistic transformation BetP is the only one which has the three properties: α-consistency, p-consistency and ulb-consistency. It is defined for any set B ⊆ Ω and B = ∅ by the following:

BetP (B) = A =∅|A⊆Ω |A ∩ B| |A| m(A) 1 -m(∅) (31)
Where |A| denotes the cardinality of A ⊆ Ω. In the case of FT analysis and closed world hypothesis (i.e. m(∅) = 0), the pignistic top event occurrence BetP ({e top }) is given by:

BetP ({e top }) = m Ωe top ({e top }) + m Ωe top ({e top , ētop }) 2 (32)
It should be noted that other probabilistic transformations were also defined using different kinds of mappings either proportional to the plausibility, to the normalized plausibility, to all plausibilities, to the belief or a hybrid mapping [J. 06].

VI. CASE STUDY

A. Description of situation

This section presents a case study that illustrates the main principles of our methodology. The model is based on a situation on which there was a near head-to-head encounter between two trains. The event took place early in the morning after some maintenance works performed during the night-shift. The idea here is to study how several organizational, technical and human factors influenced the occurrence of this situation, i.e., we want to measure the risk of a head-to-head encounter given the surrounding conditions of the scenario. In a way, it can be seen as a sensitivity analysis of how the different factors (and its epistemic uncertainty) could have affected the occurrence (risk) of the head-to-head encounter. Although it is an hypothetical situation, it was inspired by the analysis of several reports of real case studies done by the BEATT ("Bureau d'Enquêtes sur les Accidents de Transport Terrestre" French land transport accident investigation bureau).

The scenario takes place in a railway section between three stations A, B and C (cf. Figure 6). Some maintenance works were scheduled on one side of the railway track heading west (Head 2). To allow the traffic to continue, all passing trains heading in this direction were redirected to the other track in the opposite direction (Head 1). In these kind of situations, protective measures are installed to avoid a potential accident. Moreover, the circulation agenda of the commercial trains is slightly modified. The works extended over 10 km, they took place over the night and were programmed to end early in the morning in order to avoid unnecessary delays in the commercial trains. The maintenance involved two trains (named engineering trains T T X1 and T T X2). The works were supervised by a foreman and the movement of the trains was controlled by a traffic agent. Fig. 6: Description of the scenario some moments before the near accident Early in the morning, when the engineering trains where performing the final maneuvers to clear the tracks, the train T T X2 encountered head to head with a commercial train T ER going in the opposite direction. The near accident was caused by a series of unexpected events. To start with, train T T X2 was initially planned to park in a zone near station C as soon as the maintenance works were finished. A delay on the works forced to change the initial parking position of train T T X2 so as to avoid unnecessary disturbances of the early morning commercial trains. When the proposed changes were going to be implemented, the foreman realized that the access to the newly reserved parking position (Zone 3) for train T T X2 was blocked by train T T X1, therefore, he was rushed into a second change of parking positions under complicated measures. Badly advised by the traffic agent, he chose to park the train in Zone 2 in a manner that forced it to pass by a protected zone running in an opposite direction (the place of the near accident). When the train T T X2 was going to engage the protected zone, he encountered a closed signal that forbade him from crossing. The driver of the train called the traffic agent who decided to grant him the permission to cross the closed signal without taking into account the early morning time table of the commercial trains. As a consequence of these series of events, the near accident happened (cf. Figure 6).

After a deep analysis of the circumstances, four basic events are identified as the precursors of the near accident, namely:

• Bad change of parking plans : The decision to change the parking plans of the engineering train T T X2 was taken by the foreman thanks to several factors. To start with, the planning was delayed at the beginning putting some extra pressure to finish the work on time. Moreover, the traffic agent validated the change of plans with a lack of knowledge of the different parking positions of the train station and their state of occupancy. In addition, the foreman was in a very noisy environment when he took the decision. Finally, there was over-familiarity between the foreman and the traffic agent.

• Crossing permission of closed signal : The traffic agent granted the permission to cross the closed signal to the train driver because he had the illusion of a "safe" situation. Indeed, as all of the signals in the working area were closed, he thought that there was no danger. He forgot that the signals were all closed on the moment that a train entered the protected area running in the opposite direction. He maybe was not aware about the reduced overlap between signals. As a consequence, he did not take into account the time table

of the commercial trains and did not realized that a train was heading in the direction of the accident. This situation is considered a consequence of the lack of experience of the traffic agent or due to a poor traffic agent interface design. He also based his decision on the fact that others had already granted the permission to cross the same signals. He thought, "Somebody already verified the situation, it should work this time".

• Blocked road : The train T T X1 was blocking the way of train T T X2, causing more delays on the parking maneuvers. The way the initial plan was stated, this situation should not have happened because the parking track of the trains headed in an opposite direction.

• Reduction of speed at time: Finally, the traffic agent could not ask the driver of the train to reduce speed at time due to a system communication failure.

B. The near accident model

The four basic events identified are considered as the principal causes of the near accident. Two of them are considered as technical, the blocked road and the Reduction of speed at time and the two others are considered as events influenced by human and organizational factors as indicated in the following descriptions:

1) Bad change of parking plans (BR)

• Noisy environment: The foreman had to take the decision in a noisy environment which made the communication with the traffic agent complicated and created a harsh work environment.

• Traffic agent confirmation: The traffic agent confirmed that the new plan was feasible, but complicated.

In reality it was not, but the foreman trusted in his opinion.

• Delay on works: Several delays at the beginning and during the work forced the traffic agent to change the plans.

• Over-familiarity between the foreman and the traffic agent.

2) Crossing permission of closed signal (CR)

• Pre-existing confirmations: The traffic agent based his decision on the fact that others had already granted the permission to cross the same signals. He thought, "Somebody already verified the situation, it should • The traffic agent did not know of the reduced overlap between signals.

• Poor traffic agent interface design.

• Traffic agent on training: He was also in the final phase of a training period.

3) Reduction of speed at time. (RS) 4) Blocked road. (BR)

The Near Accident (N A) is illustrated with the cause effect diagram presented in Figure 7 and the fault tree presented in Figure 8. A list of the variables with their descriptions is given in Table III. As shown in the cause effect diagram, variables CR and BD were assossiated with four sub-causes each, i.e., factors that influenced the outcome of variables. On the other hand, variables RS and BR represent principal causes without sub-causes.

From the fault tree it can be seen that the four principal causes are related to the variable N A by an AND gate. where Ω BD = {BD, BD} and Ω N E = {N E, N E}.

C. Quantitative analysis of the near accident

First, Dempster combination rule is used to obtain the joint BBAs of the two ENs. Then marginalization operation is used to obtain the marginal BBAs of the variables BD and CR. Finally, using BBAs of basic events of the FT (BD, CR, RS and BR) the top event belief occurrence (the occurrence of NA) is computed using formula given in Eq. 30 and the pignistic values are computed using Eq. 32.

First, we consider that the experts give only conditional BBAs and they have no knowledge about the prior BBAs IV): Let us suppose the fact that we are certain that the traffic agent interface design was poor and there is a delay on work (Information 2). This will be represented by:

m Ωθ i ({θ i }) = 0 m Ωθ i ({θ i }) = 0 m Ωθ i ({θ i , θ i }) = 1, θ i ∈ {N E,
m ΩT C ({T C}) = 1; m ΩDW ({DW }) = 1
The risk of accident is then bounded by [0.00005, 0.06]. As we can see the belief over the occurrence of the near accident is increased because it depends significantly on the events DW and T C.

Finally suppose the fact that we are certain of the occurrence of all factors (Information 3). This will be represented by:

m ΩNE ({N E}) = 1 m ΩOF ({OF }) = 1 m ΩP C ({P C}) = 1 m ΩT R ({T R}) = 1 m ΩAT ({AT }) = 1 m ΩAT ({P T }) = 1
The risk of accident is then bounded by [0.00525, 0.06000]. The pignistic value of the risk of accident occurrence is 0.0326 (cf. Table VI). As we can expect, the belief and pignistic values of the near accident occurrence have increased.

In this probabilistic approach, we use Monte-Carlo sampling simulations to repeatedly sample component failure probabilities from the appropriate distributions, and to calculate and record the system's reliabilities. We choose uniform and normal distribution for the probabilities of events occurrence. The confidence interval of the probabilistic approach is 99%. Table VII presents the results computed by the belief functions and by probabilistic approaches. It

shows that the differences between results obtained using these two different approaches are very small. However, the width of the support defined by the belief functions approach is higher than the width of the support in the probabilistic approach because the belief approach is more conservative. By assuming a uniform and a normal distributions for probabilities of events occurrence, we introduce more uncertainty into the probabilistic approach and the results depends on the chosen law which is not the case in the belief function approach.

VII. CONCLUSIONS

This paper presents a first attempt to account for human factors in risk analysis using belief functions theory.

The ENs operations are used to elicit the masses of the basic events when they are influenced by several factors and finally, the belief over the basic events is propagated to the top undesired event.

The advantage of the presented method is that our state of belief about the conditional relationship between the different factors and the basic beliefs does not have to be perfect. Indeed, the method is well suited to account for ignorance and a priori knowledge about the basic events are not needed.

Another advantage is that the method is capable of taking into account human, organizational and technical factors in the risk analysis of railway systems.

Fig. 1 :

 1 Fig. 1: Belief and Plausibility functions

Fig. 2 :

 2 Fig. 2: EN for the risk expert example

  [0,1] expressing the expert's opinions about p with three focal elements [0.003,0.005], [0.007,0.01] and [0,1] which have mass values 0.3, 0.5, and 0.2 respectively. Then we have Bel(X = 1) = 0.003 * 0.3 + 0.007 * 0.5 + 0 * 0.2 = 0.0044 P l(X = 1) = 0.005 * 0.3 + 0.01 * 0.5 + 1 * 0.2 = 0.2065

Fig. 3 :

 3 Fig. 3: Cause Effect diagram for the example

  Fig. 4: EN for Driver route knowledge

Fig. 5 :

 5 Fig. 5: Fault Tree of railway accident

,Fig. 7 :

 7 Fig. 7: Cause Effect Diagram of the near accident

Fig. 9 :

 9 Fig. 8: Fault Tree of the near accident

ΩFig. 10 :

 10 Fig. 10: EN for Crossing Permission

  suppose that m ΩDW ({DW }) = 0.3, m ΩDW ({DW }) = 0.2 and m ΩDW ({DW, DW }) = 0.5 (Information 1). This is equivalent to the fact that degree of occurrence of the event DW is [0.3, 0.8]. The risk of accident is then bounded by [0, 0.042].

TABLE I :

 I Beliefs and plausibility measures for the example

		Initial situation	Information 1	Information 2
	R	Bel({R})	P l({R})	Bel({R})	P l({R})	Bel({R})	P l({R})
	{r}	0	1	0	0.2	0.6429	0.7143
	{r}	0	1	0.8	1	0.2857	0.3571
	{r, r}	1	1	1	1	1	1
	factory B are reliable.

TABLE II :

 II Expected p according to different methodsV. A COMBINED APPROACH BASED ON ENS AND FT ANALYSIS

	Methods	Expected p Expected p
		in Example 2
	Belief functions LME (Maximum Likelihood Estimate) x [ x n+1 , x+1 n+1 ] [0.1875,0.2500] 0.2 n Laplace x+1 0.2353 n+2 Jeffers x+1/2 n+1 0.2059

Table II, all three methods (Laplace, LME and Bayesian approaches)are contained within the upper and lower expectations given by belief functions. Furthermore, as n → +∞ the upper and lower bounds of the belief function converge, yielding a Bayesian estimate which is the same in all four methods.

TABLE III :

 III List of variables of the near accident model

	Variable Description
	NA	Near Accident
	BD	Bad change of parking plans
	NE	Noisy environment
	PT	Poor traffic agent interface design
	TC	Traffic agent confirmation
	DW	Delay on works
	OF	Over-familiarity between the foreman and traffic agent
	CR	Crossing permission of closed signal
	PC	Pre-existing confirmations
	TR	The traffic agent did not know the reduced overlap between signals
	AT	Agent on training
	RS	Reduction of speed at time
	BR	Blocked road
	work this time".

TABLE IV :

 IV Conditional and prior BBAs for the different factorsΩe i ({e i }|{θ i }) m Ω θ i ({θ i }) m Ω θ i ({θ i }) m Ω θ i ({θ i , θ i })According to Table V, in the case of no prior belief about BBAs of factors (Initial information), the risk of accident is bounded by [0, 0.06]. The pignistic value of the risk of accident occurrence is 0.03 (cf. TableVI). Now

			Conditional BBAs		Prior BBAs
	Human error e i m BD θ i NE 0.02	0	0	1
		TC	0.5	0	0	1
		DW 0.04	0	0	1
		OF	0.15	0	0	1
	CR	PC	0.1	0	0	1
		TR	0.05	0	0	1
		AT	0.2	0	0	1
		PT	0.05	0	0	1

TABLE V :

 V The obtained beliefs and plausibility measures after receiving informations

		Initial Situation	Information 1	Information 2	Information 3
	θ i	Bel({θ i })	P l({θ i })	Bel({θ i })	P l({θ i })	Bel({θ i })	P l({θ i })	Bel({θ i })	P l({θ i })
	BD	0	1	0.008	0.7	0.04	1	0.6	1
	CR	0	1	0	1	0.05	1	0.3502	1
	BR	0.5	0.6	0.5	0.6	0.5	0.6	0.5	0.6
	RS	0.05	0.1	0.05	0.1	0.05	0.1	0.05	0.1
	Near accident 0	0.06	0	0.042	0.00005	0.06	0.00525	0.06

TABLE VI :

 VI The Pignistic values after receiving informations

		Initial Situation Information 1 Information 2 Information 3
	θ i	BetP ({θ i })	BetP ({θ i })	BetP ({θ i })	BetP ({θ i })
	BD	0.5	0.35	0.52	0.8
	CR	0.5	0.5	0.55	0.6751
	BR	0.55	0.55	0.525	0.55
	RS	0.075	0.075	0.075	0.075
	Near accident 0.03	0.021	0.0	0.0326

TABLE VII :

 VII Comparison between belief functions and probabilistic approaches

		Initial Situation	Information 1	Information 2	Information 3
	Belief functions	[0, 0.06000]	[0, 0.04200]	[0.00005, 0.06000]	[0.0052, 0.06000]
	MC simulation (Uniform distribution) [0.00005, 0.03933]	[0.00005, 0.0390] [0.00041, 0.04007]	[0.00835, 0.04511]
	MC simulation (Normal distribution)	[0.00103, 0.02509]	[0.00119, 0.0259] [0.00194, 0.02576]	[0.01241, 0.03451]