Hassan Oukraou

Laura Vittadello

Virginie Coda

Charles Ciret

Massimo Alonzo

Andon A Rangelov

Nikolay V Vitanov

Germano Montemezzani

Control of adiabatic light transfer in coupled waveguides with longitudinally varying detuning

We study adiabatic light transfer in systems of two coupled waveguides with spatially varying detuning of the propagation constants, providing an analogy to the quantum phenomena of Rapid Adiabatic Passage (RAP) and two-state Stimulated Raman Adiabatic Passage (two-state STIRAP). Experimental demonstration using a photo-induction technique confirms the robust and broadband character of the structures that act as broadband directional couplers and broadband beam-splitters, respectively.

I. INTRODUCTION

The formal analogy between some quantum systems and classical photonic systems has motivated intensive research in recent years, both from a fundamental and application-oriented point of view [1,2]. In this context a strong attention was given to the coupling and passage of light between waveguides, which behaves similarly to the population transfer in quantum systems driven by an appropriate external field. Several studies in waveguides were inspired by the Stimulated Raman Adiabatic Passage process (STIRAP) [3][4][5], which represents a very e↵ective and robust way to coherently control the population of quantum states by means of a proper adiabatic temporal evolution of the coupling hamiltonian. In this way the system population can be e↵ectively driven from an initial state to a final target state. The first proposals [6,7] and the first experimental demonstration [8] of a STIRAP-like process in waveguides have stimulated several studies on adiabatic light passage in waveguides by slight modifications of this concept. These include theoretical and experimental studies related to the fractional STIRAP process [9], multi-state STIRAP [10][11][12], beamsplitting [9,[13][14][15][16], adiabatic mode conversion [11,17] the role of nonlinear e↵ects [18], or the use of such waveguide structures for photonic quantum gate operations [16]. Generally these approaches profit from the high robustness of the adiabatic process, leading for instance to a broadband behavior of the light spatial adiabatic passage process. However, since the necessary adiabatic condition is better fulfilled at longer than at shorter wavelengths [12], the use of such systems for spectral low-pass or high-pass filtering was also proposed [19,20]. ⇤ hassan.oukraou@centralesupelec.fr Previous studies principally rely on an adiabatic spatial evolution of the coupling constant C between neighboring waveguides. As follows from the coupled-mode theory [21], the latter is mainly connected to the distance between the evanescently coupled waveguides. Its counterpart in atomic physics is a temporal evolving Rabi frequency ⌦ coupling the quantum states. In the case of atomic physics the atom-laser detuning ! provides a useful and versatile parameter for the adiabatic control of the system evolution. This detuning is given by a temporal dependent mismatch of the frequency of the driving field from the transition frequency. The corresponding parameter in waveguides is represented by spatial evolution of a detuning of the propagation constants in the individual waveguides. Therefore, it is expected that also for waveguides the detuning parameter shall increase the possibilities for controlling the adiabatic evolution of the propagating light fields. Even though the interest of using tapered waveguides was recognized early for the optimization of various components [22][23][24], little attention was given to the combination with a spatially evolving coupling coe cient.

In the present work we address this issue in relation with existing quantum phenomena. We consider two coupled waveguides with spatially varying coupling coe cients and we study theoretically and experimentally the e↵ect of the detuning (z) (associated to a longitudinal variation of the index contrast) on the adiabatic light transfer between them. We choose two concrete examples. The first mimics the Rapid Adiabatic Passage (RAP) process [4] in quantum physics for which the detuning crosses zero while the coupling pulse is applied, which leads to a robust and rapid way to adiabatically invert a two-state system. Such a population transfer by adiabatic passage via a level crossing was initially implemented in nuclear magnetic resonance [START_REF] Abragam | The principles of nuclear magnetism[END_REF]. Laser-driven adiabatic passage in atoms and molecules was proposed [published in Phys. Rev A 95, 023811 (2017).]

by Treacy [START_REF] Treacy | [END_REF] and demonstrated first in the infrared by Stark-shifting the transition frequency [27], or by sweeping the laser frequency through resonance [28]. In the 80s, adiabatic passage was observed also in the nearinfrared [29] and with visible light [30]. Further details can be found in Refs. [4,[START_REF] Allen | Optical resonance and two-level atoms[END_REF][START_REF] Leidenbaum | [END_REF]. The second example mimics the so-called two-state STIRAP process [33,34], which leads finally to an equal coherent superposition of the two states. Two-state STIRAP was experimentally demonstrated by Yamazaki et al. [35] with a trapped 40 Ca + ion. In the case of waveguides the RAP process implements an achromatic directional coupler, while the two-state STIRAP process implements a broadband oneto-two beam splitter. The experiments are performed with the help of reconfigurable photoinduced waveguides recorded through the photorefractive e↵ect using the lateral illumination technique [12,36,37]. Two probing wavelengths more than 200 nm apart are used. Section II gives a general framework for the studies, section III addresses the case of the RAP-like process, while section IV treats the case of two-state STIRAP. The experimental results are in good agreement with the theoretical expectations and confirm the crucial role of the waveguide detuning.

II. GENERAL FRAMEWORK

We consider two neighboring waveguides with di↵erent longitudinal propagation constants 1 (z) and 2 (z) for their fundamental mode. In paraxial approximation, the propagation of a monochromatic light wave in such a structure can be analyzed in the framework of the coupled mode theory (CMT) [21] which treats the problem in a discrete way by involving the evanescent coupling between the waveguides. The corresponding evolution of the waves amplitudes can then be described by two coupled di↵erential equations, which, similarly as in [6,37,38], read as follows in the present case

i d dz  a 1 (z) a 2 (z) =  (z)/2 C(z) C(z) (z)/2  a 1 (z) a 2 (z)
.

Here the components of A(z) = [a 1 (z), a 2 (z)] T are, besides for an unimportant phase, proportional to the amplitudes of the fundamental modes in the two waveguides.

We call H(z) the coupling matrix in the above equation, which takes the role of the Hamiltonian and describes the interaction between the waveguide modes. It depends on the detuning parameter (z) = 2 (z) 1 (z) as well as on the space-dependent coupling constant C(z) between the two waveguides. Note that the latter corresponds to the geometrical average of the values for the coupling from waveguide 1 to 2 and from waveguide 2 to 1 [37],

C(z) = p C 12 (z)C 21 (z)
, which can be di↵erent due to the local asymmetry. Note also that the driving quantity for C(z) is the lateral distance d(z) between the waveguides. Even though an exact analytic expression cannot be given in our case, C(z) decreases nearly exponentially with d(z).

By substituting the space coordinate z for the time coordinate t in the Schrödinger-type equation ( 1), the two-waveguides system is fully equivalent to a coupled quantum two-level system, which population amplitude dynamics is described by the same equation within the rotating wave approximation [4,[START_REF] Shore | The theory of coherent atomic excitation[END_REF]. The role of the coupling constant C(z) is then taken by the Rabi frequency ⌦(t) and the waveguide detuning (z) is equivalent to the atom-laser detuning !(t).

Equation ( 1) is expressed in the natural basis formed by the modes of the two waveguides, called diabatic basis. However, we can rewrite it in another basis, called adiabatic basis, where the basis vectors are the local (instantaneous) eigenvectors of the space-varying Hamiltonian. This leads to

i d dz  b 1 (z) b 2 (z) =  ✏(z) i d✓/dz i d✓/dz ✏(z)  b 1 (z) b 2 (z) , (2)
where the quantities ✏ and ✓ are given by

✏(z) = C 2 (z) + ( /2) 2 1/2 , (3)
and

✓(z) = 1 2 arctan (2C(z)/ (z)) . (4)
The vector

B(z) = [b 1 (z), b 2 (z)]
T contains the modes amplitudes in the adiabatic basis. The amplitudes in the adiabatic and diabatic basis are connected by a space-dependent rotation operator such that A(z) = R(✓(z))B(z), with

R(✓(z)) =  cos ✓(z) sin✓(z) sin ✓(z) cos ✓(z) . (5)
For an adiabatic evolution, the amplitudes C(z) and (z) should vary slowly so that ✓(z) evolves smoothly in order to ensure that the adiabatic state vector B(z) remains fixed in the adiabatic basis. To do so, the following condition should be fulfilled [4] 1 2

@C @z C @ @z ⌧ C 2 + ✓ 2 ◆ 2 ! 3/2 . ( 6
)
When the adiabatic condition holds, the state vector B(z) remains fixed in the space-varying basis formed by the adiabatic states. In particular, if the state vector B(z) coincides with a single adiabatic state at some distance z, then it will remain in that adiabatic state as long as the evolution is adiabatic, the state vector B(z) will adiabatically follow this state.

The two eingenstates (the adiabatic states) of the rotation matrix (5) are where A 1 and A 2 are the two diabatic states of the system. In general, the state vector B(z) is aligned with an initial state at the beginning of the propagation and aligned with a specified target state at the end of the propagation. Under appropriate conditions, this alignement can be provided by a single adiabatic state or by a superposition of the adiabatic states. In the first situation, this state is called an adiabatic transfer state.

B (z) = A 1 cos ✓(z) A 2 sin ✓(z) , (7a) B + (z) = A 1 sin ✓(z) + A 2 cos ✓(z) , (7b)
In the next sections we will apply and verify the above general behavior to the specific cases of the RAP-like and two-state STIRAP-like light transfer. The waveguide structures being studied are shown schematically on Fig. 1. They are composed of two planar waveguides where the light is confined only in the plane of the drawing. One of the waveguide is straight (WG1) while the other is weakly curved (WG2), so that the distance separating them and the coupling constant C evolves with z. The propagation constant 1 of WG1 changes with z, which is achieved in our case by means of a variation of its refractive index contrast. The propagation constant 2 is left constant.

III. RAP-LIKE LIGHT TRANSFER

The waveguide structure used to simulate the RAP process is shown schematically in Fig. 1(a). The corresponding evolution of the parameters C(z) and (z) are shown in Fig. 1(b). These parameters mimic the experimental situation discussed below. For the RAP process to occur, it is essential that the detuning should cross zero during the evolution, which is the case for the function (z) in Fig. 1(b), that follows a tanh function with its zero corresponding to the maximum of C(z) (nearest distance between the waveguides). At early and FIG. 2. Schematics of the technique to create the reconfigurable waveguide structures. A control beam at the wavelength = 532 nm carries the structure of the waveguides. The latter is shaped by a spatial light modulator (SLM) and imaged to the surface of the SBN crystal, to which an electric field is applied along its polar axis. The probe beam (alternatively at = 633 nm or = 850 nm) propagates through the resulting waveguide structures to test their functionalities.

late distances the magnitude | (z)| is much larger than C(z), while the contrary is true near the zero-crossing point. Since the mixing angle ✓(z) depends strongly on the ratio C/ (see Eq. ( 4)), ✓(z) will evolve during the propagation. Initially (z = 0), this ratio tends to 0 and 2✓ is nearly equal to ⇡, so that ✓ ⇡ ⇡/2. With inscreasing z, C(z) and (z) increase ( (z) approaches zero from the negative side), which leads to a decrease of ✓(z). At half distance, (z) vanishes and C(z) reaches its maximum, so that 2✓ = ⇡/2 and ✓ = ⇡/4. In the following (z) continues to increase and becomes again much larger than C(z), which is decreasing, so that ✓ tends to zero at large distances. Consequently, starting initially from state A 1 (all light in WG1), the system follows adiabatically the adiabatic state B + (z) of Eq. (7b) and eventually ends up in state A 2 (all light in WG2). This implies that a complete light transfer should take place. This adiabatic process is intrinsically robust and is expected to take place over a wide range of the design parameters and of the propagating wavelengths for a given design, which means that the behavior should be highly achromatic.

We have verified the above expectations by using photoinduced dynamic waveguides recorded by a properly structured lateral control illumination on a Sr

x Ba 1 x Nb 2 O 6 (SBN) crystal with x = 0.61. A highly simplified conceptual scheme of this method is shown in Fig. 2. Essentially, the combination of a local illumination by the control beam and of the electric field applied to the photorefractive SBN crystal leads to a local modification of the refractive index landscape of the crystals. This modification is such that maxima of the refractive index are found at the positions of maximum illumina-tion, leading therefore to a guiding waveguide structure. Variations of the index contrast (and of the propagation constant) can be achieved by varying locally the intensity (grey level) of the control light being directed from the spatial light modulator (SLM) to the crystal. Since the recording and erasure of the structures takes only few tens of seconds in SBN, the technique provides the advantage of being able to produce reconfigurable structures in a unique crystal. Moreover, the dynamics of the waveguide recording process can be exploited to characterize the index contrast interferometrically. Since the technique was used earlier [12,36,[START_REF] Dittrich | [END_REF], we refer to the corresponding works for additional details. Note that, even though for the sake of simplicity we use here waveguides confined in one transverse dimension, the confinement in both transverse directions can be achieved too [36].

The design of the waveguide configuration used to demonstrate the RAP-like e↵ect is such that both C(z) and (z) have a smooth evolution to satisfy at best the adiabaticity criterion within the limits permitted by the physics of our experimental approach. The nominal z-dependences corresponding to Fig.

1

(b) are C(z) = C max exp[ ((z z 0 )/ 0 ) 2 ] and (z) = max tanh[2⇡(z z 0 )/L],
where C max is the maximum of the coupling constant reached at half the propagation distance z 0 = 11.5 mm, 0 = 5 mm is the 1/e half-width of the gaussian coupling-constant pulse, max is the maximum of detuning, and L = 23 mm is the full propagation distance in the crystal. The parameters C max and max depend not only on the design of the waveguide but also on the probe wavelength. Their values were estimated by carrying out preliminary experiments in non-modulated waveguides. While the observation of the coupling of waves in parallel waveguides leads to an estimation of the coupling constant, the estimation of requires to perform interferometric measurement during the formation process of the waveguides. For our experimental conditions we estimate the parameters within a precision of roughly 15%, C max = 0.35 ± 0.05 mm 1 and max = 0.52 ± 0.08 mm 1 for the probe wavelength of 633 nm, as well as C max = 0.49 ± 0.07 mm 1 and max = 0.32 ± 0.05 mm 1 for the probing at 850 nm. With the above parameters the expected evolution of the light intensity in the two-waveguide structure of Fig. 1(a)+(b) can be obtained from Eq. ( 1) and is shown in Figs. 3(a) and 3(b), for the wavelengths 633 and 850 nm, respectively. In both cases at the output of the structure one expects essentially a switch from WG1 to WG2. The corresponding experimental results are shown on Fig. 3(c) and 3(d), which depict the intensity profile at the output of the crystal. The dashed blue lines give the output of the reference waveguide, that is WG1 in absence of WG2, while the solid red lines give the output when the RAP-like structure including the WG2 is produced. The experimental profiles show some slight asymmetry and irregularities associated mainly to experimental e↵ects such as possible crystal imperfections and partial back-scattering of the recording light discussed in Ref. [41]. Nevertheless it can be recognized that for the RAP-like structure most of the light output essentially switches to WG2 for both wavelengths. This corresponds to the expectations and confirms the robustness of the approach.

It is worth verifying that the detuning is essential for the above RAP-like light transfer to occur. Figure 4 shows the expected and measured behavior for the same kind of two-waveguide structure leading to the results in Fig. 3, but in the absence of detuning ( (z) = 0). As seen in Fig. 4 in this case the RAP process does not take place and the results are not robust. An essentially complete return of the light to the straight waveguide is expected for =633 nm, while a strong transfer to the curved waveguide is expected for =850 nm (see Fig. 4(a) and(b)). This is confirmed by the experimental results of Fig. 4(c) and 4(d), respectively. Therefore, such a structure essentially leads to a damped Rabi-like oscillation between the two waveguides, with very di↵erent output distributions for the two probing wavelengths.

Even though the experiments on the RAP-like light transfer of Fig. 3 give satisfactory results, the corresponding experimental parameters are still at the edge of the zone for which the adiabatic condition is fulfilled at best in the (C max , max )-space. In order to illustrate this, we study the robustness with respect to these two parameters by plotting the expected intensity transferred to WG2 as a function of C max and max . Figure 5(a) gives a contour plot of the theoretically expected transferred intensity as a function of these two parameters. Clearly the experimental parameters (blue circle for 633 nm and green square for 850 nm) are at the border of the red region of maximum adiabaticity. The counter example in absence of detuning corresponds to points on the abscissa of the diagram of Fig. 5(a). It is also worth noting that upon an increase of the wavelength, one gets a higher coupling constant but a lower detuning. This means that the increase of wavelength corresponds to moving transversally with respect to the direction leading to a maximum fulfillment of the adiabatic condition.

The decrease in the detuning is due on one hand to the proportionality of on 1/ and, on the other hand, to a slight decrease of the electro-optic response and of the index contrast in the SBN material for longer wavelengths.

IV. TWO-STATE STIRAP-LIKE BEAM SPLITTING

The concept of two-state STIRAP [33] relies on a mathematical analogy between the equations governing the population dynamics of a two-level system with a temporally varying detuned coupling field and a resonantly coupled chainwise-linked three-level system leading to the standard STIRAP process. As discussed above, the two-level system is described by a time-dependent Schrödinger equation of two complex variables. The latter can be recast as three coupled equations of three real variables [START_REF] Shore | The theory of coherent atomic excitation[END_REF]42] in the form of the optical Bloch equation. This, on its turn, can be brought in a form equivalent to the Schrödinger equation applying to the standard STIRAP process in a three-level system. In this way, The red area correspond to the region where the adiabatic condition is fulfilled the best and the light transfer is complete. The blue circle and the green square correspond to the experimental conditions of Fig. 3 for probing at 633 nm and 850 nm, respectively. (b) Corresponding landscape for the two-state STIRAP process (case A) showing the relative intensity transferred from WG1 to WG2. The green area gives the region where the adiabatic condition is best fulfilled and the waveguide structure gives essentially a 50:50 beam splitter. The blue circle and the green square correspond to the experimental conditions of Fig. 7 for probing at 633 nm and 850 nm, respectively.

knowledge about the STIRAP process can be directly inherited for the chirped two-state excitation [5,33]. Ultimately the detuning of the two-level excitation can be identified with the Stokes pulse of a three-level STIRAP system (the Rabi-frequency pulse that couples the target state and the intermediate state), while the Rabifrequency of the two-level system is identified with the corresponding pump pulse of STIRAP, which couples the initial and the intermediate state.

For the analogy between the two-state STIRAP process and the coupling within a pair of detuned waveguides we can consider two distinct cases. The first one, Case A, is the one where the waveguide detuning (z) is spatially preceding the coupling constant pulse C(z). This case is equivalent to a counter-intuitive sequence in the standard three-state STIRAP process [3,5]. It corre-sponds to the situation shown in Fig. 1(c) and Fig. 1(d). The second case, Case B, has the coupling preceding the detuning and would correspond to an intuitive order of the pulses for an atomic three-level system.

We start by discussing Case A with the help of the general relationships given in Section II. Initially only (z) is present (see Fig. 1(d)) and its amplitude increases continuously, while C(z) is negligibly small. Therefore the ratio (z)/C(z)! + 1 and the mixing angle ✓(z) in ( 4) is initially equal to zero, as shown in Fig. 6(a). With (7a) this implies that the state B is initially aligned with A 1 . Subsequently, (z) increases towards its maximum value, as long as C(z) is still small, B deviated only little from the state A 1 . When (z) is behind its maximum and C(z) reaches comparable values, the angle ✓ rises and the adiabatic transfer state B is no longer aligned with A 1 but becomes a superposition of A 1 and A 2 . At the final stage, (L) tends to zero and (L)/C(L)! 0, so that ✓(L) tends to ⇡/4. Consequently, the system shall tend to B (L) = (A 1 A 2 )/ p 2, which leads to an equal intensity of the light in the two waveguides, as shown in Fig. 6(c). The Case B corresponds to the situation where the light is injected in WG1 from the right side of the struc-ture in Fig. 1(c). In this case, the ratio (z)/C(z) goes from 0 to +1 and ✓(z) (see figure 6(b)) from ⇡/4 to 0. Importantly, in this situation the initial state A 1 (injection in WG1) does not correspond to one of the adiabatic states B or B + , but to a superposition B ⇤ of them given by

B ⇤ (z) = 1 p 2 (B (z) + B + (z)) . (8)
Therefore, unlike for case A where only one adiabatic state is involved, here the system is expected to remain in the superposition state B ⇤ as long as the system evolves adiabatically. Since at the end of the structure ✓(L) = 0, the system shall converge towards B ⇤ (L) = (A 1 + A 2 )/ p 2, which is again an equal intensity distribution of the light in the two waveguides, as shown in Fig. 6(d). However, unlike for case A, the evolution towards this 50:50 splitting presents here some transient oscillations.

The curves shown in Fig. 6(a)-(d) are modeled according to the parameters used for the experiments given below and the following functions,

(z) = max exp[ ((z z 1 )/ 0 ) 2 ] and C(z) = C max exp[ ((z z 2 )/ 0 ) 2 ] ,
where z 1 = 9 mm and z 2 = 13.8 mm for case A, and z 1 = 13.8 mm and z 2 = 9 mm for case B. A common 1/e 2 half-width of 0 = 5 mm is used. The values of max and C max are those for the wavelength of 633 nm, which are shown as a blue point on the landscape applying to the two-state STIRAP case in Fig. 5(b). These are max = 1.2 ± 0.18 mm 1 , and C max = 0.35 ± 0.05 mm 1 . The corresponding values for the same structure but for the wavelength of 850 nm (green square in Fig. 5(b)) are max = 0.73 ± 0.11 mm 1 and C max = 0.52 ± 0.08 mm 1 , where the errors correspond to the experimental conditions. Figure 6 also contains numerical calculations of the beam propagation in waveguide structures corresponding to case A and case B using the Beam Propagation Method (BPM) [43,[START_REF] Kawano | Introduction to optical waveguide analysis solving Maxwell's equations and the Schrödinger equation[END_REF]. This second method has been widely used to design optical waveguides and predict the light evolution in such structures. The design structure only takes as an input the form and the index profile of the waveguides and does not assume anything with respect to coupled wave theory. As can be seen in Figs. 6(e) and 6(f), the BPM calculations also correctly predict the smooth evolution for case A, and the initial oscillations for case B, in agreement with the curves in Figs. 6(c) and 6(d) obtained from coupled wave theory.

Therefore, the theory predicts that at the output of the waveguide structures that mimic two-state STIRAP one should obtain a similar splitting irrespective of case A or B, and irrespective of the wavelength (see Fig. 5(b)), provided that the system fulfills su ciently well the adiabatic condition. Figure 7 gives the experimental output distributions for propagation in photoinduced waveguide structures corresponding to cases A and B at the probe wavelengths of 633 and 850 nm. It is seen that the beam splitting is satisfactorily obtained for both cases at both wavelengths, thus confirming the robustness of the pro- cess. In absence of detuning the behavior would be expected to be highly dependent on wavelength, similar to the case of Fig. 4 for RAP. Therefore, the observed broadband behavior is a direct consequence of the waveguide detuning, which permits to place the experiments outside the region near the abscissa in Fig. 5(b), for which Rabi-like oscillations strongly dependent on wavelength would be expected.

V. CONCLUSION

We have discussed the role of waveguide detuning for the adiabatic light transfer among a pair of coupled waveguides in analogy with adiabatic quantum population dynamics processes. By the examples of the Rapid Adiabatic Passage and the two-state STIRAP processes the theoretical expectations have been successfully verified experimentally using photo-induced reconfigurable waveguides recorded in a photorefractive crystal. The corresponding structures act as broadband directional couplers and broadband beam splitters, respectively. The robustness and achromaticity brought about by the adiabatic evolution mediated by the detuning and the coupling strength have been confirmed. Generally, the combination of the e↵ects of longitudinally varying coupling and longitudinally varying detuning permits to access a much wider range of possibilities for the adiabatic evolution of the light fields as compared to systems where only the coupling constant is modulated. As proposed recently for the case of three-waveguide couplers [38], this kind of concept can be extended to a larger number of waveguides, leading potentially to several new rich functionalities. The adiabatic approaches discussed here can be advantageously applied also to the case of nonclassical propagating light such as obtained from single photon sources. Also, by imposing an appropriate modulation of the mode propagation velocities, instead of standard dielectric waveguides also other types of guiding structures (e.g. photonic crystals or plasmonic waveguides) can be used in principle.

FIG. 1 .

 1 FIG. 1. (a) Schematics of the waveguide structure providing an optical analogy to the RAP process and, (b), corresponding longitudinal evolution of the coupling constant C(z) and of the detuning (z). (c) Schematics of the waveguide structure providing an optical analogy to the two-state STIRAP process and, (d), corresponding evolution of C(z) and (z) (case A). Note that on the schemes (a) and (c) the longitudinal and transverse dimensions are not in scale.

2 FIG. 3 .

 23 FIG. 3. Theoretical expectation and experimental verification of the RAP-like light transfer for the parameter given in the text. (a) Numerically calculated spatial evolution of the intensity in WG1 (blue dashed curve) and WG2 (solid red curve) for =633 nm. (b) Same for =850 nm. (c) Measured intensity output distribution with WG1 but in the absence of WG2 (reference, blue dashed line) and with the presence of both waveguides (solid red line). (d) Same for =850 nm and the same structure.

FIG. 4 .

 4 FIG. 4. Counter example to the RAP process in absence of the detuning ( (z) = 0). (a) Expected evolution of the intensity in WG1 (blue dashed curve) and WG2 (solid red curve) according to Eq. (1) for =633 nm. (b) Same for =850 nm. (c) Measured intensity output distribution with WG1 but in the absence of WG2 (reference, blue dashed line) and with the presence of both waveguides (solid red line). (d) Same for =850 nm and the same structure.

FIG. 5 .

 5 FIG. 5. (a) Landscape of the expected output intensity of WG2 for the RAP-process as calculated using coupled mode theory under variation of the parameters Cmax and max.

FIG. 6 .

 6 FIG. 6. Theoretical expectations for a two-state STIRAP-like waveguide structure modeled for propagation of light at the 633 nm wavelength. Panels (a), (c) and (e) correspond to case A, panels (b), (d) and (f) to case B. (a) and (b) show the evolution of the mixing angle ✓(z) with the propagation distance z. (c) and (d) show the intensity evolution in WG1 (dashed blue line) and WG2 (solid red line) obtained from coupled mode theory. Panels (e) and (f) show the corresponding evolution of the light distribution obtained with the alternative beam propagation method (BPM). For both cases the system leads to an equal splitting of the input intensity among the two waveguides.

FIG. 7 .

 7 FIG. 7. Experimental output intensity distribution for the case of two-state STIRAP as probed for cases A and B at two di↵erent wavelengths. (a) and (b): probing at = 633 nm, (c) and (d): probing at = 850 nm. Panels (a) and (c) correspond to case A, while panels (b) and (d) correspond to case B (see text).

ACKNOWLEDGMENTS

This work has been supported by the Bulgaria-France bilateral program RILA and by the Region Lorraine.


P Lefebvre

K Tse

Ve Koon

E Brusseau

S Nicolle

J.-F Pallierne

Simon A Lambert

D Grenier

Comparison of viscoelastic property characterization of plastisol phantoms with magnetic resonance elastography and high-frequency rheometry

I. INTRODUCTION

Changes in the mechanical properties of biological tissues can reveal the development of a pathology. This consideration has led to the development of elastographic techniques, such as Magnetic Resonance Elastography (MRE) or ultrasound elastography, that are dedicated to the investigation of these mechanical properties in vivo and noninvasively [START_REF] Muthupillai | Magnetic resonance elastography by direct visualization of propagating acoustic strain waves[END_REF]- [START_REF] Bercoff | Supersonic Shear Imaging: A new technique for soft tissue elasticity mapping[END_REF]. In MRE, shear waves are induced in the tissue of interest using an external driver. Tissue motion is encoded in phase images thanks to a motion-sensitive MR imaging sequence and these images are then processed to map and quantify the viscoelastic properties of the examined medium. Magnetic resonance elastography remains the object of many studies and is investigated for different organs, such as liver, brain or muscle [START_REF] Yin | Assessment of hepatic fibrosis with magnetic resonance elastography[END_REF]- [START_REF] Murphy | Decreased brain stiffness in Alzheimer's disease determined by magnetic resonance elastography[END_REF]. P. M. Lefebvre, K. Tse Ve Koon, S. A. Lambert, D. Grenier and E. Brusseau are with the University of Lyon, CREATIS, CNRS UMR5220, Inserm U1206, INSA-Lyon, University Claude Bernard Lyon 1, 69621 Villeurbanne, France. (corresponding author: 334-26-23-44-99; e-mail: pauline.lefebvre@creatis.insa-lyon.fr).

S. Nicolle is with the Univ Lyon, Université Claude Bernard Lyon 1, IFSTTAR, UMR_T9406, LBMC, F69622, Lyon, France.

J.-F. Palierne is with the University de Lyon, Laboratoire de Physique, CNRS UMR5672, ENS de Lyon, Lyon, France.

Generally speaking, assessing new imaging methods involves the use of phantoms mimicking the properties of biological tissues. For elastography in particular, various types of phantoms (agarose [START_REF] Muthupillai | Magnetic resonance elastography by direct visualization of propagating acoustic strain waves[END_REF], polyvinyl alcohol cryogel [START_REF] Sinkus | High-resolution tensor MR elastography for breast tumour detection[END_REF], plastisol [START_REF] Chakouch | Development of a phantom mimicking the functional and structural behaviors of the tight muscles characterized with magnetic resonance elastography technique[END_REF]- [START_REF] Leclerc | Characterization of a hyper-viscoelastic phantom mimicking biological soft tissue using an abdominal pneumatic driver with magnetic resonance elastography (MRE)[END_REF], copolymer-in-oil [START_REF] Oudry | Copolymer-in-oil phantom materials for elastography[END_REF]…) can be found in the literature, yielding a good range of elastic modulus. Deeper assessment can be performed by comparing the MRE results with either those of other elastographic techniques or with some reference mechanical tests [START_REF] Vappou | Magnetic resonance elastography compared with rotational rheometry for in vitro brain tissue viscoelasticity measurement[END_REF]- [START_REF] Okamoto | Viscoelastic properties of soft gels: a comparison of magnetic resonance elastography and dynamic shear testing in the shear wave regime[END_REF]. However, it is usually quite difficult to compare the techniques in a common frequency range. For example, in [START_REF] Vappou | Magnetic resonance elastography compared with rotational rheometry for in vitro brain tissue viscoelasticity measurement[END_REF] [START_REF] Vappou | Magnetic resonance elastography compared with rotational rheometry for in vitro brain tissue viscoelasticity measurement[END_REF], rheometry measurements were performed up to 10 Hz, while MRE data were acquired for frequencies between 80 Hz and 140 Hz. Similarly, [START_REF] Ringleb | shear wave magnetic resonance elastography: comparison to a dynamic shear material test[END_REF] [12] presented MRE results measured at frequencies between 100 Hz and 200 Hz, complemented with dynamic mechanical analyses from 10 Hz to 50 Hz. The study realized by [START_REF] Okamoto | Viscoelastic properties of soft gels: a comparison of magnetic resonance elastography and dynamic shear testing in the shear wave regime[END_REF] presents yet an overlapping frequency range: from 100 to 400 Hz for MRE and from 20 Hz to 200 Hz for dynamic shear test [START_REF] Okamoto | Viscoelastic properties of soft gels: a comparison of magnetic resonance elastography and dynamic shear testing in the shear wave regime[END_REF].

Working in a common frequency range allows a direct comparison of results: the purpose of this study is to present MRE and High-Frequency Rheometry (HFR) measurements in a common frequency range. Results are provided for three plastisol phantoms with different softener concentrations. These phantoms have indeed the main interest of presenting a viscoelastic behavior and are therefore good candidates to mimic biological tissues.

II. MATERIALS AND METHODS

A. Preparation of phantoms

Three homogeneous phantoms were prepared with different concentrations of standard plastisol (which is a suspension of polyvinyl chloride particles in a liquid plasticizer) and softener (which is a Di-2-ethyl-hexyl adipate) (Plastileurre Standard and Assouplissant, Bricoleurre, France).

The combinations of liquid preparations were: phantom (A) 100% standard plastisol, (B) 75% standard plastisol -25% softener, and (C) 50% standard plastisol -50% softener. Each solution was heated to 160°C, poured into cubic molds (35 x 35 x 35 mm) for MRE and on aluminum foil so as to get thin samples as required for rheological tests. They were then cooled until the phantom solidified (about one hour). The phantoms were kept at room temperature (20°C). MRE and HFR experiments were performed one week after phantom preparation.

Comparison of viscoelastic property characterization of plastisol phantoms with magnetic resonance elastography and highfrequency rheometry

P. M. Lefebvre, K. Tse Ve Koon, E. Brusseau, S. Nicolle, J.-F. Palierne, S. A. Lambert, D. Grenier

B. Magnetic resonance elastography data acquisition, image processing and data analysis

MRE experiments were run on a Bruker 4.7 T small animal MRI system. Phantoms were placed inside a volume coil for both radiofrequency emission and reception. The MRE sequence developed was a 2D gradient-echo sequence with the following parameters: repetition time = 250 ms, echo time = 20 ms, scan resolution = 128 x 128 pixels, field of view (FOV) = 4 cm. Four vibration phase offsets were recorded for each excitation frequency. Shear waves were induced in the phantoms using a noninvasive custom-made device maintaining the medium between two plates, the upper one being fixed while the lower one being activated by a piezoelectric actuator (CEDRAT Technologies), as shown in Fig. 1. This piezoelectric device can cover a large frequency range up to 15 kHz. Acquisitions were conducted at different excitation frequencies, ranging from 400 to 1,200 Hz: below 400 Hz, reconstruction is difficult due to long wavelengths (above 1 cm), and at 1,200 Hz, waves are strongly attenuated. Finally, during MRE experiments, temperature was roughly 22°C.

For each excitation frequency, 2D phase images were obtained twice at opposite motion-encoding gradient polarity and then subtracted. The resulting image was Fourier transformed and the complex wave image U(x,y,ω) was taken at the excitation frequency ω. Next a band-pass Butterworth filter was applied in order to remove noise and exclude wavelength values that are not in the range of possible values with the considered type of medium. Finally, under common assumptions (linearity, isotropy and local homogeneity), storage (G') and loss (G'') moduli (corresponding to the real and imaginary parts of the complex shear modulus G*) were determined by applying a 2D Helmholtz inversion algorithm, according to the following equation:

G = -ρ ω² U / ΔU, ()
ρ being the density of the materials (close to 1000 kg/m 3 ) and Δ the Laplacian operator, approximated using a simple central difference scheme.

For each phantom, a region of interest (ROI), excluding the edges of the phantoms, was defined on the obtained elastograms, and G' and G'' mean values as well as standard deviations were extracted from this ROI. The storage modulus G' represents the elasticity of the material whereas the loss modulus G'' is associated to the viscous properties of this material.

C. High-frequency rheometry measurements

A high-frequency rheometer [START_REF] Nicolle | Shear linear behavior of brain tissue over a large frequency range[END_REF]- [START_REF] Constantin | Highfrequency rheological behaviour of a multiconnected lyotropic phase[END_REF] was used for mechanical tests. The sample is contained between two glasses plates mounted on piezoelectric. The upper plate is vertically driven by the piezoelectric element in a sinusoidal motion. A squeezing flow is consequently induced in the sample and the stress transmitted to the second plate is measured by the other piezoelectric ceramic. The complex modulus G* is directly deduced from the pressure exerted by the sample on the captor and the amplitude of the motion [START_REF] Hébraud | Role of permeation in the linear viscoelastic response of concentrated emulsions[END_REF] and is thus rheological model independent. With this rheometer, the sample has to have a diameter of about one cm and has to be very thin (ideally less than one mm in thickness) in order to ensure an optimal measurement. In comparison with standard rheometers whose frequencies are usually up to 10 Hz for very soft materials, this piezorheometer is designed to perform measurements in a very large frequency range, from 0.1 Hz to 10 kHz. These measurements were realized at 25°C and in the same week than MRE experiments. For each phantom, four samples were tested and the mean and standard deviation were calculated. Results presented hereafter are those obtained at frequencies ranging from 160 to 630 Hz, providing thus an overlapping frequency range with MRE.

III. RESULTS

A. Wave propagation

Fig. 2 shows the obtained phase images at different excitation frequencies for phantom A. It illustrates clearly the wave propagation, and displays the dependence of wavelength and attenuation on frequency. Due to the viscous properties of the phantom, the amplitude of the wave decreases gradually and this is amplified at higher frequencies. As can be expected, the wavelength decreases when the frequency increases. B. Reproducibility of the magnetic resonance elastography measurements MRE tests were repeated three times on phantom (A). After an initial acquisition, which served as reference, a second acquisition was performed immediately after, the phantom remaining in the same position. The third acquisition (acquisition 3) was recorded ten days later. This last acquisition therefore takes also into account the variability induced by repositioning the phantom. As can be seen from Fig. 3, storage and loss moduli are similar for the three acquisitions.

C. Magnetic resonance elastography and high-frequency rheometry results

MRE and HFR results are presented in Fig. 4 for the three different phantoms (A), (B) and (C). Fig. 4a) displays the storage modulus G' whereas Fig. 4b) represents the loss modulus G''. Results with the two modalities are in quite good agreement for both G' and G''. For example, for phantom (A) at 400 Hz, G' = 33.5 ± 1.9 kPa with MRE technique and G' = 32.7 ± 0.7 kPa with HFR data. The G'' estimation presents significant standard deviations compared to its mean value.

Additionally, different ranges of storage and loss moduli were found for the three studied phantoms: with the MRE results, at 400 Hz, the storage modulus G' for phantom (A), (B) and (C) are respectively 33.5 ± 1.9 kPa, 11.5 ± 1.3 kPa and 3.1 ± 0.5 kPa. Concerning the loss modulus G'', it is respectively equal to 6.7 ± 2.2 kPa, 3.9 ± 1.1 kPa and 1.24 ± 0.5 kPa for each phantom.

IV. DISCUSSION

Results are presented for excitation frequencies varying between 400 Hz and 1,200 Hz for MRE, and 160 Hz and 630 Hz for HFR. In MRE, frequencies needed to be sufficiently high to be suited to the small-sized media investigated, and sufficiently low to avoid too strong shear waves attenuation. In HFR, although the device used can cover a large range of frequencies (up to 10 kHz), results are provided up to 630 Hz due to inertial effects with the plastisol phantoms at higher frequencies. It should also be added that sample thicknesses achieved were larger (between 1 and 2.5 mm) than the recommended size which can be responsible for less accurate results. Despite these difficulties, storage and loss modulus measurements over an overlapping frequency range were obtained, showing good agreement between MRE and HFR.

Acquisitions with both techniques were performed with as close as possible conditions: same phantom preparations, similar temperatures during acquisitions and experiments taking place the same week. However, the size of phantoms remains a main difference between the two techniques: cubic phantoms for MRE experiments and very thin samples for HFR. HFR phantoms need to be regular in thickness whereas plastisol phantoms are difficult to be cut properly, this explains why phantom preparations were poured in cubic molds for MRE, and on an aluminum foil for HFR and then allowed to freely spread out in order to obtain thin layers. These two processes may induce variations in the resulting viscoelastic properties of phantoms, as, for example, the cooling rate is different. A first test was led by considering two phantoms with different sizes (35 x 35 x 35 mm and 35 x 35 x 15 mm) cooled at room temperature, and a second one with three other cubic plastisol phantoms cooled respectively at room temperature (20°C), in a fridge (2°C) and in an oven (at 60°C). Storage and loss moduli were here determined with the MRE technique only. Similar results were obtained for all these five phantoms.

Viscosity property quantification in MRE is quite difficult [START_REF] Bohte | MR Elastography of the Liver: Defining Thresholds for Detecting Viscoelastic Changes[END_REF] in particular because G'' estimation is highly affected by noise. Thus the G'' estimations presented here have significant standard deviation compare to their mean values. However, these values are consistent with the ones obtained with HFR, suggesting that MRE is a technique that could give an estimation of the loss modulus G'' and thus of the medium viscosity.

Finally, different ranges of storage and loss moduli were found for the three studied phantoms, depending on the softener concentration. These values are in the storage and loss moduli range for different biological tissues, such as liver, brain or muscle [START_REF] Chakouch | Development of a phantom mimicking the functional and structural behaviors of the tight muscles characterized with magnetic resonance elastography technique[END_REF]. Porcine white matter storage modulus (measured in vitro) varies from 5 to 10 kPa between 100 and 1,000 Hz and loss modulus from 2 to 6 kPa [START_REF] Nicolle | Shear linear behavior of brain tissue over a large frequency range[END_REF]. The storage modulus in thin rat liver slices was measured with variable fibrosis extent: results ranged between 2.5 and 4.5 kPa, at 600 Hz, depending on the METAVIR score [START_REF] Ronot | Viscoelastic parameters for quantifying liver fibrosis: threedimensional multifrequency MR elastography study on thin liver rat slices[END_REF]. A fibrotic human liver was found to present a storage modulus varying between 2 to 30 kPa between 100 and 1,000 Hz and a loss modulus varying from 0.2 to 10 kPa [START_REF] Riek | Wide-range dynamic magnetic resonance elastography[END_REF]. Varying the concentrations in phantoms could be a way to easily mimic elasticity and viscosity of different healthy and pathological tissues.

V. CONCLUSION

This paper presents a comparison of magnetic resonance elastography and high-frequency rheometry measurements, within an overlapping frequency range and with similar experimental conditions. The results obtained using these two techniques satisfactorily match, confirming that the complex modulus of a viscoelastic phantom can be reliably assessed using MRE.

Moreover, different ranges of storage and loss moduli were obtained by varying softener concentrations during phantom preparation. These values are in the range of storage and loss moduli for different tissues. Plastisol phantoms could also be helpful to assess new elastographic methods and tools adapted to different soft tissues, as well as to evaluate MRE parameters before in vivo experiments.

Figure 1 .

 1 Figure 1. MRE driver used to generate shear waves in the phantoms with a piezoelectric actuator.

Figure 2 .

 2 Figure 2. Wave images for phantom (A) at different excitation frequencies: 500, 700, 900 and 1200 Hz.

Figure 3 .

 3 Figure 3. MRE results obtained for phantom A and three different acquisitions: Acquisition 1 is the reference acquisition; acquisition 2 was performed right after 1. Acquisition 3 was recorded ten days later.

Figure 4 .

 4 Figure 4. Storage (a) and loss (b) moduli for the three different phantoms and the two techniques, magnetic resonance elastography (dashed lines) and high-frequency rheometry (dotted lines). Results for phantoms (A), (B) and (C) appear in blue, green and orange, respectively.

ACKNOWLEDGMENT

This work was carried out within the framework of the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR).


Sophie Mercier
email: sophie.mercier@univ-pau.fr

Meier-Hirmer Mercier

Roussignol

Modelling track geometry by a bivariate Gamma wear process, with application to maintenance

This paper discusses the maintenance optimization of a railway track, based on the observation of two dependent randomly increasing deterioration indicators. These two indicators are modelled through a bivariate Gamma process constructed by trivariate reduction. Empirical and maximum likelihood estimators are given for the process parameters and tested on simulated data. The EM algorithm is used to compute the maximum likelihood estimators. A bivariate Gamma process is then fitted to real data of railway track deterioration. Preventive maintenance scheduling is studied, ensuring that the railway track keeps a good quality with a high probability. The results are compared to those based on both indicators taken separately, and also on one single indicator (usually taken for current track maintenance). The results based on the joined information are proved to be safer than the other ones, which shows the interest of the bivariate model.

INTRODUCTION

This paper is concerned with the maintenance optimization of a railway track, based on the observation of two dependent randomly increasing deterioration indicators. The railway track is considered as deteriorated when any of these two indicators is beyond a given threshold. The point of the paper is the study of preventive maintenance scheduling, which must ensure that, given some observations provided by inspection, the railway track will remain serviceable until the next maintenance action with a high probability.

Track maintenance is a very expensive task to accomplish. Consequently, it is essential to carry out maintenance actions in an optimal way, while taking into account many parameters: safety and comfort levels to be guaranteed, available logistic means, . . . The earlier the deterioration is detected, Longitudinal levelling (NL)

Transversal levelling (NT) the easier it is to schedule maintenance actions. The objective is therefore to develop a good prediction model. Deterioration of track geometry is characterized by the development of different representative parameters like, for example, the levelling of the track. Figure 1 shows the defects that are measured by two of these parameters: the longitudinal (NL) and transversal (NT) levelling indicators.

At the SNCF (French National Railways), track inspections are programmed annually on a national level. The interval between two inspections on high speed tracks is currently about two weeks, the inspections are carried out by a modified high-speed train. The collected time series are transformed into indicators that sum up the state of the track over each km. These new indicators are referred to as synthesized Mauzin data. Numeric Mauzin data are available since the opening of the French high-speed lines.

Usually, the synthesized Mauzin indicator of the longitudinal levelling (NL indicator) is used for maintenance issues: thresholds are fixed for this indicator in order to obtain a classification of the track condition and to fix dates for maintenance operations. For example, an intervention should be scheduled before the NL indicator exceeds 0.9.

Based on expert judgements, a Gamma process has been used in [START_REF] Meier-Hirmer | Modèles et techniques probabilistes pour l'optimisation des stratégies de maintenance[END_REF] both to model the evolution of the NL indicator and to plan maintenance actions. As noted by J.M. van Noortwijk in his recent survey [START_REF] Van Noortwijk | A survey of the application of gamma processes in maintenance[END_REF], this kind of process is widely used in reliability studies (see also [START_REF] Abdel-Hameed | A gamma wear process[END_REF], [START_REF] Grall | Continuous-time predictive-maintenance scheduling for a deteriorating system[END_REF] and [START_REF] Zuckerman | Optimal replacement policiy for the case where the damage process is a one-sided Lévy process[END_REF]). Various domains of applications exist, such as civil engineering ( [START_REF] Buijs | Timedependent reliability analysis of flood defences using gamma processes[END_REF], [START_REF] Ciampoli | Stochastic process for extrapolating concrete creep[END_REF]), highway engineering [START_REF] Nicolai | A comparison of models for measurable deterioration: an application to coatings on steel structures[END_REF] or railway engineering [START_REF] Meier-Hirmer | Maintenance optimisation for system with a gamma deterioration process and intervention delay: application to track maintenance[END_REF]. Gamma processes are also used in other domains, such as finance [START_REF] Joshi | Intensity Gamma: a new approach to pricing portfolio credit derivatives[END_REF] or risk analysis [START_REF] Dufresne | Risk theory with the gamma process[END_REF]. All these papers use univariate Gamma processes.

In the present case, as the two indicators NL and NT are dependent, the use of a bivariate model is required. For this purpose, different processes might be used, such as Bessel [START_REF] Newby | A bivariate process model for maintenance and inspection planning[END_REF] or Lévy processes [START_REF] Kallsen | Characterization of dependence of multidimensional Lévy processes using Lévy copulas[END_REF]. In this paper, the approach of F.A. Buijs, J.W. Hall, J.M. van Noortwijk and P.B. Sayers in [START_REF] Buijs | Timedependent reliability analysis of flood defences using gamma processes[END_REF] is used: a specific Lévy process called bivariate Gamma process is considered. This process is constructed from three independent univariate Gamma processes by trivariate reduction, and has univariate Gamma processes as marginal processes.

It is the first time that both NL and NT indicators are used conjointly to predict the optimal dates of maintenance actions. The objective is to analyse "book" -2009/9/22 -17:52 -page 3 -#3 the correlation between the two processes and to determine in what circumstances this bivariate process allows a better prediction of the maintenance times than the current univariate one, based only on the NL indicator.

The paper is organized in the following way: bivariate Gamma processes are introduced in Section 2. Empirical and maximum likelihood estimators for their parameters are provided in Section 3. An EM algorithm is proposed to carry out the maximum likelihood estimation. Both methods are tested on simulated data. Section 4 is devoted to the study of preventive maintenance planning and to the comparison of the results based on the bivariate and on the univariate models. Finally, a bivariate Gamma process is fitted to real data of railway track deterioration in Section 5 and it is shown that the preventive maintenance scheduling based on the two available deterioration indicators are clearly safer than those based on a single one, or on both taken separately.

THE BIVARIATE GAMMA PROCESS

Recall that an univariate (homogeneous) Gamma process (Y t ) t≥0 with parameters (α, b) ∈ R * 2

+ is a process with independent increments such that Y t is Gamma distributed Γ (αt, b) with probability density function (p.d.f.)

f αt,b (x) = b αt Γ (αt) x αt-1 e -bx 1 R+ (x) , E (Y t ) = αt b , Var(Y t ) = αt b 2
for all t > 0, and Y 0 ≡ 0 (see [START_REF] Van Noortwijk | A survey of the application of gamma processes in maintenance[END_REF] for more details). Following [START_REF] Buijs | Timedependent reliability analysis of flood defences using gamma processes[END_REF], a bivariate Gamma process X t t≥0 = X (1) t , X

(2) t t≥0 is constructed by trivariate reduction: starting from three independent univariate Gamma processes Y (i) t t≥0 with parameters (α i , 1) for i ∈ {1, 2, 3} and from b 1 > 0, b 2 > 0, one defines:

X (1) t = Y (1) t + Y (3) t /b 1 , and X (2) t = Y (2) t + Y (3) t /b 2 for all t ≥ 0. The process (X t ) t≥0 = X (1)
t , X

(2) t t≥0 is then a homogeneous process in time with independent increments and it is a Lévy process. The marginal processes of (X t ) t≥0 are univariate Gamma processes with respective parameters (a i , b i ), where a i = α i + α 3 for i = 1, 2.

For any bivariate Lévy process, the correlation coefficient ρ Xt of X

(1) t and X

(2) t is known to be independent of t. For a bivariate Gamma process, one obtains:

ρ = ρ Xt = α 3 √ a 1 a 2 and α 1 = a 1 -ρ √ a 1 a 2 , α 2 = a 2 -ρ √ a 1 a 2 , α 3 = ρ √ a 1 a 2 .
This entails

0 ≤ ρ ≤ ρ max = min (a 1 , a 2 ) √ a 1 a 2 . ( 1
)
"book" -2009/9/22 -17:52 -page 4 -#4

See [START_REF] Devroye | Non-Uniform Random Variate Generation[END_REF] section XI.3 for results on bivariate Gamma distributions. This leads to two equivalent parameterizations of a bivariate Gamma process:

(α 1 , α 2 , α 3 , b 1 , b 2 ) and (a 1 , a 2 , b 1 , b 2 , ρ).
With the parameterization (α 1 , α 2 , α 3 , b 1 , b 2 ), the joint p.d.f. of X t is:

g t (x 1 , x 2 ) = b 1 b 2 min(b1x1,b2x2) 0 f α1t,1 (b 1 x 1 -x 3 ) f α2t,1 (b 2 x 2 -x 3 ) f α3t,1 (x 3 ) dx 3 , = b 1 b 2 e -b1x1-b2x2 Γ (α 1 t) Γ (α 2 t) Γ (α 3 t) × • • • × min(b1x1,b2x2) 0 (b 1 x 1 -x 3 ) α1t-1 (b 2 x 2 -x 3 ) α2t-1 x α3t-1 3 e -x3 dx 3 . ( 2
)

PARAMETER ESTIMATION

The data used for the parameter estimation are values of the process increments for non overlapping time intervals on a single trajectory, and also on different independent trajectories. The data can then be represented as ∆t j , ∆X

j (ω) , ∆X (1)
j (ω) 1≤j≤n where ∆t j = t j -s j stands for a time increment and ∆X

(i) j = X (i) tj -X (i) sj for the associated i-th marginal increment (i = 1, 2). For different j, the random vectors ∆X (1) j , ∆X (2)
j are independent, but not identically distributed. The random variable ∆X (i)

j (i = 1, 2)
is Gamma distributed with parameters (a i ∆t j , b i ). The joint p.d.f. of the random vector ∆X

(1) j , ∆X

(2) j is equal to g ∆tj (., .), with ∆t j substituted for t in (2). In the same way as for parameter estimation of a (univariate) Gamma process, both empirical and maximum likelihood methods are possible in the bivariate case.

Empirical estimators

Using E ∆X (i) j = ai bi ∆t j and Var ∆X

(i) j = ai b 2 i ∆t j for i = 1, 2 and for all j, empirical estimators (â 1 , b1 , â2 , b2 ) of (a 1 , b 1 , a 2 , b 2
) are given in [START_REF] Ciampoli | Stochastic process for extrapolating concrete creep[END_REF] and [START_REF] Van Noortwijk | A stochastic deterioration process for time-dependent reliability analysis[END_REF], with:

âi bi = n j=1 ∆X (i) j t n and âi b2 i = n j=1 ∆X (i) j -âi bi ∆t j 2 t n -1 tn n j=1 (∆t j ) 2 , ( 3
)
where we set

t n = n j=1 ∆t j . Using Cov ∆X (1) j , ∆X (2) j = ρ √ a 1 a 2 b 1 b 2 ∆t j ,
"book" -2009/9/22 -17:52 -page 5 -#5

Modelling track geometry by a bivariate Gamma wear process a similar estimator ρ may be given for ρ, with:

ρ √ â1 â2 b1 b2 = n j=1 ∆X (1) j -â1 b1 ∆t j ∆X (2) j -â2 b2 ∆t j t n -1 tn n j=1 (∆t j ) 2 . ( 4
)
These estimators satisfy:

E âi bi = a i b i , E âi b2 i = a i b 2 i , E ρ √ â1 â2 b1 b2 = ρ √ a 1 a 2 b 1 b 2 .
If the time increments ∆t j are equal, these estimators coïncides with the usual empirical estimators in the case of i.i.d. random variables.

Maximum likelihood estimators

The parameter estimation of an univariate Gamma process is usually done by maximizing the likelihood function (see e.g. [START_REF] Meier-Hirmer | Modèles et techniques probabilistes pour l'optimisation des stratégies de maintenance[END_REF]). With this method, estimators āi and bi (i = 1, 2) of the marginal parameters are computed by solving the equations:

āi bi = n j=1 ∆X (i) j n j=1 ∆t j and n j=1 ∆t j × ln āi n j=1 ∆t j n j=1 ∆X (i) j + n j=1 ∆t j ln ∆X (i) j -ψ (ā i ∆t j ) = 0, where ψ (x) = Γ (x) Γ (x) , Γ (x) = ∞ 0 e -u u x-1 du
for all x > 0 (ψ is the Digamma function).

In order to estimate all the parameters of the bivariate process (α 1 , α 2 , α 3 , b 1 , b 2 ) (which are here prefered to (a 1 , b 1 , a 2 , b 2 , ρ)), the likelihood function associated with the data ∆t j , ∆X (1)

j , ∆X (2) j 1≤j≤n can be written as L(α 1 , α 2 , α 3 , b 1 , b 2 ) = n j=1 g ∆tj (∆X (1) j , ∆X (2)
j ). However, because of the expression of the function g t (., .), it seems complicated to optimize this likelihood function directly. An EM algorithm (see [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm (with discussion)[END_REF]) is then used, considering ∆Y

(3) j = Y (3) tj -Y (3) sj
1≤j≤n as hidden data. This procedure is still too complicated to estimate all the five parameters and does not work numerically. So, the procedure is restricted to the three parameters (α 1 , α 2 , α 3 ). For the parameters b 1 , b 2 , the values b1 , b2 computed using the maximum likelihood method for each univariate marginal process are taken.

In order to simplify the expressions, the values of the data ∆t j , ∆X

j , ∆X (2) j , ∆Y (3) j 1≤j≤n are denoted by t j , x (1) j , x (2) j , y (1)
(3) j 1≤j≤n , the associated n-dimensional random vectors by X (1) , X (2) , Y (3) and the associated n-dimensional data vectors by x (1) , x (2) , y (3) .

"book" -2009/9/22 -17:52 -page 6 -#6

Mercier, Meier-Hirmer & Roussignol

The joint p.d.f. of the random vector X

(1)

t , X (2) t , Y (3) t is equal to: b 1 b 2 f α1t,1 (b 1 x 1 -y 3 ) f α2t,1 (b 2 x 2 -y 3 ) f α3t,1 (y 3 ) = b 1 b 2 e -(b1x1+b2x2) Γ (α 1 t) Γ (α 2 t) Γ (α 3 t) (b 1 x 1 -y 3 ) α1t-1 (b 2 x 2 -y 3 ) α2t-1 y α3t-1 3 e y3 , with 0 ≤ y 3 ≤ min (b 1 x 1 , b 2 x 2 ), x 1 > 0 and x 2 > 0. Then, the log-likelihood function Q x(1) , x(2) , ȳ (3)
associated with the complete data x (1) , x (2) , y (3) is derived:

Q x(1) , x(2) , ȳ(3) = n (ln (b 1 ) + ln (b 2 )) -• • • n j=1 (ln Γ (α 1 t j ) + ln Γ (α 2 t j ) + ln Γ (α 3 t j )) -b 1 n j=1 x (1) j -• • • b 2 n j=1 x (2) j + n j=1 (α 1 t j -1) ln b 1 x (1) j -y (3) j + • • • (α 2 t j -1) ln b 2 x (2) j -y (3) j + (α 3 t j -1) ln y (3) j + y (3) j .
For the EM algorithm, the conditional log-likelihood of the complete data given the observed data is needed:

E Q X(1) , X(2) , Ȳ (3) | X(1) = x(1) , X(2) = x(2) = n (ln (b 1 ) + ln (b 2 )) -b 1 n j=1 x (1) j -b 2 n j=1 x (2) j + • • • n j=1 ((α 1 t j -1) E ln b 1 x (1) j -∆Y (3) j |∆X (1) j = x (1) j , ∆X (2)
j = x (2) j + (α 2 t j -1) E ln b 2 x (2) j -∆Y (3) j |∆X (1) j = x (1) j , ∆X (2) j = x
(2) j

+ (α 3 t j -1) E ln ∆Y (3) j |∆X (1) j = x (1) j , ∆X (2) j
= x (2) j +E Y (3) j |∆X (1) j = x (1) j , ∆X (2) j
= x (2) j - n j=1
(ln Γ (α 1 t j ) + ln Γ (α 2 t j ) + ln Γ (α 3 t j )) .

(

Finally, the conditional density function of

Y (3) t given X (1) t = x 1 , X (2) t = x 2 is equal to: f α1t,1 (b 1 x 1 -y 3 ) f α2t,1 (b 2 x 2 -y 3 ) f α3t,1 (y 3 ) min(b1x1,b2x2) 0 f α1t,1 (b 1 x 1 -x 3 ) f α2t,1 (b 2 x 2 -x 3 ) f α3t,1 (x 3 ) dx 3 = (b 1 x 1 -y 3 ) α1t-1 (b 2 x 2 -y 3 ) α2t-1 y α3t-1 3 e y3
min(b1x1,b2x2) 0

(b 1 x 1 -x 3 ) α1t-1 (b 2 x 2 -x 3 ) α2t-1 x α3t-1 3 e x3 dx 3
, "book" -2009/9/22 -17:52 -page 7 -#7

where 0 ≤ y 3 ≤ min (b 1 x 1 , b 2 x 2 ), x 1 > 0 and x 2 > 0.

Step k of the EM algorithm consists of computing new parameter values (α

(k+1) 1 , α (k+1) 2 , α (k+1) 3
) given the current values (α

(k) 1 , α (k) 2 , α (k)
3 ) in two stages:

• stage 1: compute the conditional expectations in (5) using the current set (α

(k) 1 , α (k) 2 , α (k)
3 ) of parameters, with:

f 1 j, α k) 1 , α (k) 2 , α (k) 3 = E ln b1 x(1) j - Ȳ (3) j | X(1) = x(1) j , X(2) = x(2) j , f 2 j, α k) 1 , α (k) 2 , α (k) 3 = E ln b2 x(2) j - Ȳ (3) j | X(1) = x(1) j , X(2) = x(2) j , f 3 j, α k) 1 , α (k) 2 , α (k) 3 = E ln Ȳ (3) j | X(1) = x(1) j , X(2) = x(2) j , h α k) 1 , α (k) 2 , α (k) 3 = n j=1 E Ȳ (3) j | X(1) = x(1) j , X(2) = x(2) j .
• stage 2: take for (α

(k+1) 1 , α (k+1) 2 , α (k+1) 3
) the values of (α 1 , α 2 , α 3 ) that maximize [START_REF] Zuckerman | Optimal replacement policiy for the case where the damage process is a one-sided Lévy process[END_REF], which here becomes:

g α 1 , α 2 , α 3 , α (k) 1 , α (k) 2 , α (k) 3 = n ln b1 + ln b2 -b1 n j=1 x (1) j -b2 n j=1 x (2) j + n j=1 (α 1 t j -1) f 1 j, α (k) 1 , α (k) 2 , α (k) 3 + (α 2 t j -1) f 2 j, α (k) 1 , α (k) 2 , α (k) 3 + (α 3 t j -1) f 3 j, α (k) 1 , α (k) 2 , α (k) 3 - n j=1 (ln Γ (α 1 t j ) + ln Γ (α 2 t j ) + ln Γ (α 3 t j )) + h α (k) 1 , α (k) 2 , α (k) 3 .
The maximization in stage 2 is done by solving the following equation with respect to α i :

∂g α 1 , α 2 , α 3 , α (k) 1 , α (k) 2 , α (k) 3 ∂α i = n j=1 t j f i j, α (k) 1 , α (k) 2 , α (k) 3 - n j=1 t j ψ (α i t j ) = 0 (6) for i = 1, 2, 3.
In the same way, it is possible to take the values ā1 , ā2 , b1 , b2 obtained by maximum likelihood estimation on the univariate marginal processes for (a 1 , a 2 , b 1 , b 2 ) and to estimate only the last parameter α 3 by the EM algorithm. In that case, α

(k+1) 3
is the solution of the equation:

n j=1 t j f 3 j, α (k) 1 , α (k) 2 , α (k) 3 -f 1 j, α (k) 1 , α (k) 2 , α (k) 3 -• • • f 2 j, α (k) 1 , α (k) 2 , α (k) 3 - n j=1 t j ψ (α 3 t j ) -ψ ((ā 1 -α 3 )t j ) -• • • ψ ((ā 2 -α 3 )t j ) = 0.

Tests on simulated data

500 time increments (t j ) 1≤j≤500 are randomly chosen similar to the data of track deterioration (the proposed methods will be used on these data in Section 5). Then, 500 values of a bivariate Gamma process are simulated corresponding to these time increments and with parameters a 1 = 0.33, a 2 = 0.035, b 1 = 13.5, b 2 = 20 and ρ = 0.5296. These parameter values have the same order of magnitude than those observed for track deterioration studied in Section 5. Three series of 500 data points are simulated independently. Results of parameters estimation are given in Tables 1, 2 and 3, each corresponding to a series of data. In these tables, one can find: the true values in column 2, the empirical estimators in column 3, the univariate maximum likelihood estimators of a 1 , b 1 , a 2 , b 2 in column 4, the EM estimator of the three parameters a 1 , a 2 , ρ in column 5, using the parameters b1 , b2 previously estimated by the univariate maximum likelihood method (from column 4), and the second EM estimator of the parameter ρ in column 6, using the estimated parameters ā1 , b1 , ā2 , b2 from column 4.

The initial values for the EM algorithm are different for the three tables. For Table 1, the EM algorithm has been initiated with α 2 and3, α and ρ (k) along the different steps of the EM algorithm, one may note that the parameters a (k) i stabilize more quickly than the parameter ρ (k) (about 5 iterations for a (k) i and between 20 and 30 iterations for ρ (k) ).

The conclusion of this section is that estimation of parameters (a i , b i ) by empirical and maximum likelihood methods give satisfactory results, with a slight preference to maximum likelihood. Empirical estimators of ρ have a good order of magnitude, but are sometimes not precise enough. Estimators of ρ obtained by EM are always reasonable. The estimation of the three parameters (α 1 , α 2 , α 3 ) (column EM1) seems to give slightly better results "book" -2009/9/22 -17:52 -page 9 -#9 than those of the estimation of the parameter α 3 alone (column EM2). The results obtained by the EM algorithm for parameters a i (column EM1) are good, with a quality quite similar to those obtained by univariate maximum likelihood estimation. Finally, the EM algorithm does not seem sensitive to initial values, at least if the initial value of α 3 is not too small.

Modelling track geometry by a bivariate Gamma wear process

PREVENTIVE MAINTENANCE PLANNING

A bivariate Gamma process

X t = X (1) t , X (2) t
is now used to model the development of two deterioration indicators of a system. We assume that there exist (corrective) thresholds s i (i = 1, 2) for each indicator, above which the system is considered to be deteriorated. The system is not continuously monitored but only inspected at time intervals, with a perfect observation of the deterioration level. When one (or both) indicator(s) is observed to be beyond its corrective threshold, an instantaneous maintenance action is undertaken, which brings the system back to a better state, not necessarily as good as new. When both indicators are observed to be below their corrective thresholds or after a maintenance action, a new inspection is planned. The time to next inspection (τ ) must ensure with a high probability that neither X

(1) t nor X

(2) t go beyond their corrective thresholds s i before the next inspection.

Let (

x 1 , x 2 ) ∈ [0, s 1 [×[0, s 2 [
be the observed deterioration level at some inspection time, say at time t = 0 with no restriction. (If x 1 ≥ s 1 or x 2 ≥ s 2 , a maintenance action is immediately undertaken).

For i = 1, 2, let T (i) be the hitting time of the threshold s i for the marginal process

X (i) t t≥0 , with T (i) = inf t > 0 : X (i) t ≥ s i . Also, let ε ∈]0, 1[ be some confidence level.
Different points of view are possible: in the first case, τ (i) , i = 1, 2 is the "book" -2009/9/22 -17:52 -page 10 -#10 (1) , τ (2) , τ U and τ B (ρmax).

time to next inspection associated to the marginal process X (i) t t≥0 , with

τ (i) = max τ ≥ 0 such that P xi T (i) > τ ≥ 1 -ε ,
where P xi stands for the conditional probability given X (i) 0 = x i . One then gets:

P xi T (i) > τ (i) = 1 -ε.
Without a bivariate model, a natural time to next inspection for the system is: (1) , τ (2) .

τ U = max τ ≥ 0 s.t. P x1 T (1) > τ ≥ 1 -ε and P x2 T (2) > τ ≥ 1 -ε , = min τ
Using a bivariate Gamma process, the time to next inspection becomes: τ B = max τ ≥ 0 such that P (x1,x2) T (1) > τ, T (2) > τ ≥ 1 -ε .

The goal is to compare τ U and τ B , and more generally, to understand the influence of the dependence between both components on τ B . Using

P xi T (i) > t = P xi X (i) t < s i = P 0 X (i) t < s i -x i = F ait,bi (s i -x i ) ,
where F ait,bi (x) is the cumulative distribution function of the distribution Γ (a i t, b i ), the real τ (i) is computed by solving the equation F aiτ (i) ,bi (s i -x i ) = 1 -ε, for i = 1, 2, and τ U = min τ (1) , τ (2) is derived. Similarly,

P (x1,x2) T (1) > t, T (2) > t = P (0,0) X (1) t < s 1 -x 1 , X (2) t < s 2 -x 2 , = s1-x1 0 s2-x2 0 g t (y 1 , y 2 ) dy 1 dy 2 , ≡ G t (s 1 -x 1 , s 2 -x 2 ) ,
where g t is the p.d.f. of X t (see [START_REF] Van Noortwijk | A survey of the application of gamma processes in maintenance[END_REF]). This provides τ B by solving With a 1 = 0.03, b 1 = 20, = 0.5 and s 1 = s 2 = 1, and different values for a 2 , b 2 , x 1 and x 2 , Table 4 gives the corresponding values for ρ max (as provided by ( 1)) and the resulting τ (1) , τ (2) , τ U and τ B (ρ max ). The value of τ B is plotted with respect to ρ in the Figures 2 for the two different cases of Table 4, and the corresponding value of τ U is indicated.

G τ B (s 1 -x 1 , s 2 -x 2 ) = 1 -ε. "book" -
In both figures, one can observe that with all other parameters fixed, the bivariate preventive time τ B is an increasing function of ρ, such that τ B ≤ τ U . Also, both τ B = τ U and τ B < τ U are possible. The theoretical proof of such results is not provided here because of the reduced size of the present paper, but will be provided in a forthcomming one.

In conclusion to this section, one can see that using a bivariate model instead of two separate univariate models generally shortens the time to next inspection (τ B ≤ τ U ). This means that taking into account the dependence between both components provides safer results. Also, the optimal time to next inspection is increasing with dependence (τ B increases with ρ), which implies that the error made when considering separate models (τ U ) is all the more important the less the components are dependent. This also implies that the safest attitude, in case of an unkown correlation, is to consider both components as independent and chose τ = τ ⊥ , where

τ ⊥ = max τ ≥ 0 such that P x1 T (1) > τ P x2 T (2) > τ ≥ 1 -ε .

APPLICATION TO TRACK MAINTENANCE

A bivariate Gamma process is now used to model the evolution of the two track indicators NL and NT (see the Introduction) and times to next inspection are computed, as described in the previous section.

Using univariate maximum likelihood and EM methods on data corresponding to the Paris-Lyon high-speed line provide the estimations â1 = 0.0355, b1 = 19.19, â2 = 0.0387, b2 = 29.72, ρ = 0.5262. Usual thresholds are "book" -2009/9/22 -17:52 -page 12 -#12 Figure 3. τ (1) , τ (2) and τ B with respect to x2 with x1 = 0.4 s 1 = 0.9 for NL and s 2 = 0.75 for NT. With these values, τ (1) , τ (2) and τ B are plotted in Figure 3 with respect of x 2 when x 1 is fixed (x 1 = 0.4). In that case τ (1) = 150.

This figure shows that taking into account the single information x 1 = 0.4 may lead to too late maintenance actions. As an example, if x 2 = 0.4, one has τ B = 134.7 (and τ (2) = 152.9). The preventive maintenance action based only on NL is consequently scheduled 15 days too lately. If x 2 = 0.5, one obtains τ B = 95.9 (τ (2) = 97.5) and the maintenance action is undertaken 54 days too late. If x 2 = 0.6, one obtains τ B = 47.1 (τ (2) = 47.2) and this is 103 days too late.

Concluding this section, one can finally observe that if x 1 is not too close to x 2 , the value τ U = min τ (1) , τ (2) seams reasonable for maintenance scheduling (see Figure 3), contrary to the currently used τ (1) , which may entail large delay in its planning (more than 100 days in our example). If x 1 is close to x 2 , the values of τ U and τ B have the same order of magnitude, however with τ U > τ B , so that the preventive maintenance action is again planned too lately (15 days in the example).

CONCLUSION

A bivariate Gamma process has been used to model the development of two deterioration indicators. Different estimation methods have been proposed for the parameters and tested on simulated data. Based on these tests, the best estimators seem provided by univariate likelihood maximization for the marginal parameters and by an EM algorithm for the correlation coefficient.

Preventive maintenance scheduling has then been studied for a system that deteriorates according to a bivariate Gamma process. In particular, "book" -2009/9/22 -17:52 -page 13 -#13 it has been shown that, given an observed bivariate deterioration level, the optimal time to maintenance is increasing with dependence. It has been proven that the optimal time to maintenance is always shorter when taking into account the dependence between the two deterioration indicators than when considering them separately (or only considering one of them).

Finally, a bivariate Gamma process has been used to study a real track maintenance problem. The application shows that when both observed deterioration indicators are close to each other, the bivariate process gives safer results for maintenance scheduling than both univariate processes considered separately or one single univariate process, with the same order of magnitude in each case however. When the observed deterioration indicators are clearly different, considering one single univariate process as it is done in current track maintenance, may lead to clearly unadaquate results. This application to real data of railway track deterioration hence shows the interest of a bivariate model for a correct definition of a maintenance strategy.

Bibliography

Figure 1 .

 1 Figure 1. Levelling defects

  .1 and ρ (0) = 0.75). For Tables

Table 1 .

 1 Results for the first series of data.

		True	Empirical	Univariate	EM algorithm
		values	estimators	max likelihood	EM1	EM2
	a1	0.0330	0.0348	0.0342	0.0347	-
	b1	13.5	14.38	14.14	-	-
	a2	0.0350	0.0362	0.0357	0.0354	-
	b2	20	20.58	20.25	-	-
	ρ	0.5296	0.5637	-	0.5231	0.5214
		True	Empirical	Univariate	EM algorithm
		values	estimators	max likelihood	EM1	EM2
	a1	0.0330	0.0315	0.0326	0.0328	-
	b1	13.5	12.80	13.16	-	-
	a2	0.0350	0.0357	0.0361	0.0365	-
	b2	20	20.25	20.54	-	-
	ρ	0.5296	0.5750	-	0.5272	0.5257

Table 2 .

 2 Results for the second series of data.

Table 3 .

 3 Results for the third series of data.

			True		Empirical		Univariate		EM algorithm
			values	estimators	max likelihood	EM1	EM2
		a1	0.0330	0.0297			0.0340	0.0343	-
		b1	13.5		11.71			13.43		-	-
		a2	0.0350	0.0340			0.0385	0.0389	-
		b2	20		18.79			21.28		-	-
		ρ	0.5296	0.5645			-	0.5060	0.5027
		a2	b2	x1	x2	ρmax	τ (1)	τ (2)		τ U	τ B (ρmax)
	case 1	0.03	30	0.2	0.2	1		341.12	558.31	341.12	341.12
	case 2	0.04	20	0.4	0.2	0.866	237.33	255.84	237.33	229.91

Table 4 .

 4 Two different combinations of values for a2, b2, x1 and x2, and the resulting ρmax, τ

  2009/9/22 -17:52 -page 11 -#11Figure 2. τ B with respect to ρ and τ U , for the four cases of Table 4

			Modelling track geometry by a bivariate Gamma wear process
	341.2									238
										236
	341.1									234
	341									232
										230
	340.9
										228
	340.8									226
	340.7									224
										222
	340.5 340.6	0	0.2	0.4	ρ	0.6	0.8	τ B τ U	1	220 0 218	0.1	0.2	0.3	0.4	ρ	0.5	0.6	0.7	τ B 0.8 τ U
				case 1									case 4


Gabriele Boschetto

Stefania Carapezzi

Aida Todri-Sanial

Quantum Mechanical Simulations of 2D Materials for Unconventional Computing and Biosensing Applications

Quantum mechanical simulations of 2D materials for unconventional computing and biosensing applications

Gabriele Boschetto 1 , Stefania Carapezzi 1 , Aida Todri-Sanial 1 1 LIRMM, Université de Montpellier, CNRS, Montpellier, France Atomically thin two-dimensional (2D) materials have been -and are still currently beingextensively studied due to their unique mechanical, electrical, and optical properties, which enable the development of innovative devices and technologies. Within the wide chemical space of transition metal dichalcogenides (TMDs), single-layer molybdenum disulphide (MoS2) is with no doubt one of the most studied materials due to its stability and its direct optical band gap of 1.8 eV, which make it the ideal candidate to be used in a wide range of nanoelectronic devices, going beyond conventional CMOS technology.

Here we look at MoS2 in the context of two different cutting-edge applications: first, as the core material for memristors for neuromorphic computing, and then as the sensing material in field-effect biosensors. We want to bridge the gap between materials' properties and device physics, and to do so we carry out atomistic computer simulations in the framework of density functional theory (DFT). We study a planar MoS2 device which, depending on its geometry, can be exploited either as a memristor device or field-effect biosensing device.

2D memristors based on MoS2 present several advantages with respect to conventional memristors based on transition metal oxides, such as good flexibility, high transparence, and the potential to work when applying low voltages (0.1 -0.2 V), thus allowing the fabrication of compact and energy-efficient devices. [START_REF] Qiao | Fabricating molybdenum disulfide memristors[END_REF] Monolayer-MoS2 memristors ("atomristors") are technologically relevant, although there is still much debate regarding the physics of their resistive switch, which is thought to depend on both surface defects and on the type of metal-semiconductor junction. [START_REF] Ge | Atomristor: nonvolatile resistance switching in atomic sheets of transition metal dichalcogenides[END_REF][START_REF] Kim | Zero-static power radio-frequency switches based on MoS atomristors[END_REF] To shed light on the memristive mechanism, we systematically investigate the effect of some of the most common defects found in experiments (for instance, sulphur mono-and di-vacancies, adatoms) on the electronic properties of MoS2 in the vacuum and also in the presence of a metal top contact. Our study constitutes the first step of a wider multi-scale modelling approach in which the aim is to construct a full atomistic-to-device level model.

Recently, MoS2 has also been studied as a sensing platform for the detection of gas and small biological molecules, such as glucose. [START_REF] Jeevanandham | Nickel oxide decorated MoS2 nanosheet-based non-enzymatic sensor for the selective detection of glucose[END_REF] Enzymatic biosensing is the most common approach, however non-enzymatic sensing can provide a higher stability of the sensor and a prompt response, potentially at a lower cost but at the expense of chemical selectivity. Here we are interested in the non-enzymatic detection of cortisol in human sweat, as a mean to monitor the risk of cardiovascular diseases. However, it is not known if such analyte interacts in a suitable way with MoS2, and if so, how. Thus, we thoroughly explore the MoS2/cortisol interaction in terms of both structural and electronic properties to assess the possible sensing mechanism. We also study the impact of some of the most used metal dopants employed in lab-scale experiments, such as Ni, Pt, Au, Ag, in order to improve the sensing platform with respect to bare MoS2.

Overall, our work ultimately aims at obtaining a deep understanding of the properties of MoS2 when used in different applications, in order to drive the design of devices towards better performance.