C Milián

A Jarnac

Y Brelet

V Jukna

A Houard

A Mysyrowicz

A Couairon

Effect of input pulse chirp on nonlinear energy deposition and plasma excitation in water

We analyze numerically and experimentally the effect of the input pulse chirp on the nonlinear energy transfer from 5 µJ fs-pulses at 800 nm to water. Numerical results are also shown for pulses at 400 nm, where linear losses are minimized, and for different focusing geometries. Input chirp is found to have a big impact on the transmitted energy and on the plasma distribution around focus, thus providing a simple and effective mechanism to tune the electron density and energy deposition. We identify three relevant ways in which plasma features may be tuned.

I. INTRODUCTION

Laser energy deposition in condensed dielectric media has many applications ranging from micromachining of glasses [1][2][3] and medical laser surgery [4][5][6][7] to bubble formation [8,9] and sound wave generation for oceanography [8,[START_REF] Jones | Remote Intense Laser Acoustic Source[END_REF][START_REF] Jarnac | Laser-Ultrasonics LU 2013[END_REF]. The initial stage of laser energy deposition consists in the generation of a localized weakly ionized plasma with typical density of one electron per hundreds of molecules in the focal region of the laser beam (see, e.g., Refs. [START_REF] Gamaly | [END_REF]13] for reviews on laser-plasma interaction in solids and Refs. [7,[14][15][16][17][18] and references therein for works in water and other liquids).

To deposit laser energy in a well defined focal volume far from the surface, femtosecond pulses carrying energies of few µJ may be used in conjunction with tight focusing geometries to avoid nonlinear effects prior to the focus. Smaller numerical apertures trigger nonlinearities at earlier stages and are closely related to the phenomenon of laser beam filamentation [START_REF] Couairon | Self-focusing and filamentation of femtosecond pulses in air and condensed matter: simulations and experiments[END_REF][START_REF] Jarnac | [END_REF], but this allow for laser energy deposition at deeper distances from the surface (see, e.g., Refs. [21][22][23][24] for investigations in water and other liquids in the last few years).

We are interested in femtosecond laser energy deposition in water for the above-mentioned potential applications that require a well controlled localized plasma in a focal volume at depths from a few centimeters to deeper positions under the surface. From a practical point of view, liquids provide a platform where the localized plasma tracks do not lead to permanent damage since they are naturally erased via electron-ion recombination. This allows for consecutive independent material excitations at laser repetition rates ν kHz, which are much lower than the typical hydrodynamic inverse time scales for material recovery, ν hydro ∼ MHz (see recent experimental results in Ref. [8]).

The idea that input pulse chirp has a strong impact on nonlinear dynamics has been widely used during the last decade. For media exhibiting normal group velocity dispersion (GVD), an optimal negative input chirp * carles.milian@cpht.polytechnique.fr makes equal the spatial focusing and temporal compression lengths, yielding enhancement of the nonlinear effects resulting in high intensities, long plasma channels and broader spectra [25][26][27][28][29], generation of few-cycle [30,31] and ultra-short [32] pulses, and the possibility for remote spectroscopy [33]. Control of input chirp has also been reported to enhance pulse collision induced spectral broadening [34] and damage tracks in solids [35,36].

In this work, we explore numerically and experimentally the effects the input pulse chirp and focusing conditions have on nonlinear energy deposition from µJ pulses at 800 nm to water and on the electron-plasma density distribution. Comparison between numerical and experimental results for the transmission (laser energy deposition) presents an excellent agreement. Numerical simulations let us acquire a deeper understanding of the spatiotemporal dynamics and have access to experimentally inaccessible data, such as the generated plasmas, fluence distributions, intensity profiles, etc. We then identify different plasma generation regimes. In particular we find three different input pulse widths that maximize different features of the plasmas. First, the minimum of the optical transmission corresponds to negatively pre-chirped pulses that generate the plasmas with the maximum possible energy. Second, further negative pre-chirping results in a plasma volume with maximized length (and still relatively high electron density), and third, plasma densities are maximized for even larger negative values of the input chirp, at the expense of the plasma channel length. We foresee that the possibility to control the density and shape of the plasmas in the focal region is promising for developing the aforementioned applications.

The rest of this paper is organized as follows. Section II describes the experimental setup and Section III presents the theoretical model used for numerical simulations of laser energy deposition. Section IV presents the results and shows the effects on energy transmission induced by changes in the focalization geometry and chirp. Section V is devoted to the extension of the numerical results for pulses which carrier wavelength is 400 nm. Final conclusions and remarks are expounded in Section VI.

II. EXPERIMENTAL SETUP

The experimental setup used for measurements and modeled below is shown in Fig. 1. The experiment was performed by using a commercial CPA Ti:Sapphire femtosecond Laser (THALES Alpha 100) delivering t pl = 45 fs transform limited pulses with carrier vacuum wavelength λ 0 = 800 nm at a repetition rate of ν = 100 Hz. The chirped pulses are obtained by detuning the compressor stage integrated within the laser. The beam is focused inside the BK7 glass water tank with a lens of focal distance d = 7.5 cm in air placed at a = 1 cm from it. The walls of the tank are 1 cm thick and the inner length L w = 10 cm. Two power-meters are used to monitor simultaneously the input (1) and output (2) pulse energies. Their calibration with the 80 % reflectivity beam splitter (BS) provides accurate (and repeatable) measurements of the transmission as a function of input pulse chirp (see Fig. 3) and let us keep the energy of the pulses entering the tank relatively constant: E in = 5.0 ± 0.1 µJ. Results presented here are for distilled water, however we produced essentially identical data with tap water.

III. THEORETICAL MODEL

We model the propagation of the electric field envelope in Fourier space, E(ω, r, z) ≡ F [E(t, r, z)], along the z coordinate by means of a unidirectional pulse envelope propagation equation (see, e.g., Ref. [37] for details),

∂ E ∂z = i K(ω) + ∆ ⊥ 2k(ω) E + iQ(ω) P 2ǫ 0 - J n(ω)2ǫ 0 c , (1)
where k(ω) ≡ n(ω)ω/c and n(ω) is the real valued refractive index of water (for data see Ref. [38]). Linear and nonlinear dispersion functions read the carrier frequency ω 0 = 2πc/λ 0 , and β 1 is the linear attenuation coefficient. The cylindrically symmetric Laplacian ∆ ⊥ ≡ ∂ 2 r + r -1 ∂ r accounts for diffraction. In the nonlinear terms of the above equation, the instantaneous Kerr nonlinearity is accounted for through the polarization P(ω, r, z) ≡ 2ǫ 0 n 0 n 2 F [IE], where n 0 ≡ n(ω 0 ), n 2 is the nonlinear index, and I ≡ ǫ 0 cn 0 |E| 2 /2 the electric field intensity. The effects associated to quasi-free electrons are contained in the current J (ω, r, z) ≡ J MP A + J P L . J MP A ≡ ǫ 0 cn 0 β K [1 -ρ/ρ nt ]I K-1 E accounts for multiphoton absorption (MPA), where K denotes the number of absorbed photons, β K the nonlinear loss coefficient, and ρ nt the density of neutral molecules. The interaction of light with the plasma of electron density ρ(t, r, z) is described by J P L ≡ ǫ 0 cσ(ω 0 ) F [ρE]. The cross section σ(ω) may be found from the Drude model as [START_REF] Couairon | Self-focusing and filamentation of femtosecond pulses in air and condensed matter: simulations and experiments[END_REF]37,39]

K(ω) ≡ k(ω) + iβ 1 /2 -k 0 -k ′ 0 [ω -ω 0 ] and Q(ω) = ω 2 /[k(ω)c 2 ], respec- tively, where k 0 ≡ k(ω 0 ), k ′ 0 ≡ ∂ ω k(ω)
σ(ω) ≡ ω 2 0 τ c cρ c 1 1 -iωτ c , (2)
where τ c is related to the mean collision time, ρ c ≡ ω 2 0 m e ǫ 0 /e 2 is the critical density at which the plasma becomes opaque, m e the electron mass, and e the elementary charge. The two terms in J P L associated with the real and imaginary parts of σ(ω) = σ a + iσ d account for plasma absorption and plasma defocusing, respectively. Equation ( 1) is coupled to the rate equation [16,[START_REF] Couairon | Self-focusing and filamentation of femtosecond pulses in air and condensed matter: simulations and experiments[END_REF]43],

∂ρ ∂t = β K K ω 0 I K 1 - ρ ρ nt + σ a (ω 0 ) U i ρI - ρ τ r . (3)
From left to right, the terms on the right hand side in Eq.

(3) describe multiphoton ionization (MPI), avalanche ionization (inverse Bremsstrahlung), and electron-ion recombination over the typical time scale τ r . Note validity of Eq. ( 3) is restricted to the weak plasma condition: ρ ≪ ρ nt . A list of parameters and their values is provided in table I. The model provided here disregards back reflection of light. Recent numerical studies done with Maxwell solvers on plasma generation in water micro-droplets [18] clearly show that plasma induced reflections are of the order of 1% or less. Our Gaussian pulses have a peak power P = 2/πE in /t p ≈ 3P cr when they are transform limited (P cr ≈ 30 MW is the critical power at 800 nm, see appendix) and therefore the revolution symmetry of Eq. ( 1) is justified.

The origin of the propagation coordinate, z, is chosen at the first air-glass interface of the water tank (see Fig. 1). At this position, Eq. ( 1) is initialized with the pulse

E(t, r, 0) = E 0 exp - r 2 w 2 0 1 + i k 0 w 2 0 2f 0 - t 2 t 2 p0 [1 + iC] , (4)
where w 0 , f 0 , t p0 and C denote the input beam width, curvature, pulse duration and chirp, respectively. For comparison with experimental results, we refer below to the beam width w F W HM and pulse duration, t p , at the position of the focusing lens (z = -a). Under our geometrical conditions the generation of plasma is highly localized around the nonlinear focus, which position and size are around z N L 8.4 cm and ∆z N L 0.1 cm, respectively (see e.g., Figs. 2(b) and 5), and therefore J ≈ 0 in Eq. ( 1) during most of the propagation in water. Under these conditions, the beam waist typically decreases from w 0 by two or three orders of magnitude before reaching the focal region but the beam and the pulse remain approximately Gaussian. Strong space-time reshaping of the pulse occurs mainly in the focal region. Numerical integration of Eq. ( 1) in the region prior to the focus can therefore be advantageously replaced by a less intensive numerical integration of the propagation by using the moment method [44][45][46] up to z 0 ≈ 8 cm (see appendix A for details).

This strategy allowed us to perform a parametric study leading to a good match between numerical and experimental data (see Sec. IV B), a task that would have been computationally much more expensive with the full model Eqs. (1-3) due to the fine spatiotemporal resolution required in the numerical grids to convey an input 7 mm wide beam through the focus. We checked with benchmarks that results are not significantly affected by the use of the moment method in the first propagation stage.

IV. NONLINEAR ENERGY DEPOSITION AND

PLASMA EXCITATION AT λ0 = 800 nm

In section IV A we review the numerical results regarding the effects on energy transmission induced by changes in the focusing geometry. We then show in Sec. IV B that further tuning of the transmission and plasma volumes is achieved simply by modifying the input pulse chirp, in the geometry used in experiments.

A. Influence of beam width and pulse energy

Figure 2(a) shows the optical transmission obtained for different focusing geometries, i.e., varying w F W HM from 3 to 7 mm and E in from 0.1 to 5 µJ (pulses are initially un-chirped: C = 0, t p = t pl = 45 fs). The drop in transmission when increasing either E in or w F W HM is due to the enhancement of multiphoton and/or plasma absorption. This happens because the localization of light imposed by the lens and further enhanced by self-focusing on the transverse plane, is more efficient than all other processes that yield decrease of intensity, such as pulse dispersion [47] and plasma defocusing [48,49]. water (B, D) present a high intensity region showing one (or several) flat plateaus. In the latter cases the beam does not diffract as fast as in the former ones due to the prolongated effect of self-focusing, which keeps the beam waist relatively narrow despite the presence of plasma. Indeed, the range of parameters in the case D suggests it falls in the filamentation regime [50].

B. Influence of the input pulse chirp

The effect of input pule chirp, C = ±[{t p /t pl } 2 -1] 1/2 , was studied for E in = 5 µJ and w F W HM = 7 mm, corresponding to the experimental conditions (the relatively wide beam is required to keep intensity below the BK7 glass damage threshold [51]). The input pulse width was varied from t p = 45 fs (C = 0) to t p ≈ 4 ps, for positive and negative chirp. Figure 3 shows the transmission as a function of the input pulse width obtained numerically, by integrating Eqs. ( 1)-(3) (stars), and experimentally, in the setup of Fig. 1 (circles). Note all input pulses have exactly the same bandwidth. Moreover, because we are in the deep normal GVD regime of water self phase modulation (SPM) is the main frequency conversion effect which does not widen substantially the spectra at the energy levels and focusing conditions used here. Therefore, all results presented below are interpreted solely in terms of the spatiotemporal dynamics of the different pulses and spectral broadening effects are dismissed in all discussions (see, e.g., Ref. [52] for recent results on supercontinuum generation in water by pumping close to the zero GVD, and Refs. [53,54] for reviews on the topic in one-dimensional systems). Results in Fig. 3 present a minimum transmission for negatively chirped pulses (t p ≈ -100, -150 fs, the minus sign stands for negative chirp). This is expected because water exhibits normal GVD at λ ∼ 800 nm, and therefore the maximum intensity levels at focus (maximizing nonlinear losses) will be achieved when the pulses reach this region with nearly compensated chirp C(z N L ) ≈ 0. For highly chirped pulses the levels of transmission tend ∼ 0.82, in excellent agreement with the predicted value in the linear regime, exp(-β 1 L w ) ≈ 0.82.

The good agreement simulations present with experiments was obtained via a parametric study which involved varying β K and τ c . An important aspect to consider before a parametric study is the relative importance of MPA and plasma absorption in the overall transmission. In the slowly varying envelope model Eqs. ( 1)-(3), the instantaneous optical kinetic power density transferred to the medium due to MPA and plasma absorption is given by (see, e.g., [START_REF] Haus | Electromagnetic noise and quantum optical measurements[END_REF]) and the loss ratio l pl /l MP A , shown in Fig. 4(a), is calculated from:

W(t, r, z) ≡ 1 2 Re {JE * } = β K I K 1 - ρ ρ nt + σ a ρI, (5)
l pl ≡ rdrdzdt W βK =0 , l MP A ≡ rdrdzdt W σa=0 . ( 6
)
Whilst avalanche ionization governs nonlinear losses for large values of chirp, MPA acquires its maximum importance around t p ≈ -100, -150 fs, where transmission is minimized (see Fig. 3). Note from Fig. 4(a) that this is rather general for the relatively wide range of β K , τ c values. Therefore, in first place τ c was chosen as to match experimental data (within the error bars) in the transmission asymptotes, |t p | 3000 fs. Figure 4(b) shows the typical effect the collision time, τ c , has on the transmission for large input chirp, C. Large values of collision time, τ c ≫ 10 fs, tend to remove plasma absorption (see Fig. 4(d)) and the transmission tends to the linear one. This highlights the fact that M P A is not relevant for large values of input pulse chirp. In second place, β K was adjusted to match experiments around the minimum of transmission, t p ≈ -100, -150 fs. The resulting pair of parameters (β K , τ c ) obtained by this method (see table I or Fig. 3) provided a good quantitative agreement with experimental results.

The reason why β K had to be decreased by an order of magnitude, with respect to the value predicted by Keldysh ( βK in Table 1), in order to match experiments is due to the fact that the ionization rate does not entirely fall in the MPI regime for the typical intensity levels involved here, as it can be seen from Figs. 4(c) and 5(b). In Fig. 5(b), the intensity levels above the first Keldysh rate fall (at I ≈ 1.6 × 10 13 W/cm 2 ) are marked by the thick trace. For these high intensities, the linear trend of log W vs log I implies that using an effective MPA coefficient, β K , in Eqs. (1-3) is equivalent as to the full Keldysh rate. It could be argued that this effective value of β K does not apply for the lower intensity situations, for which I max < 1.6 × 10 13 W/cm 2 (E out /E in 0.7 in Fig. 3) because it underestimates MPA losses. However, as discussed above, nonlinear losses are typically dominated by plasma absorption in this region by one order of magnitude or more and MPA is responsible for only a few percent of the total losses. Note that we are not fitting here the Keldysh rate, but an effective MPA coefficient that mimics the Keldysh rate at the intensity levels of interest, i.e., those for which MPA induced losses acquire their maximum importance.

Figure 5 shows maxima of the electron plasma density and intensity around the nonlinear focus for several input pulse widths, t p . Comparison with Fig. 3 reveals several interesting features. First, we note that the minimum of absorption in Fig. 3 corresponds exactly with the case in which the maximum intensity is achieved (see Fig. 5(b)). This might thought to be expected because indeed MPA seeds more free electrons at higher intensities and avalanche effects induce further losses to the optical field. However, the overall absorption is a spa- tiotemporal integrated quantity whereas the maximum intensity represents only a local feature. The fact that there is a correspondence in between the local and integrated quantities highlights the importance of strong MPA in the overall losses for locally high intensity levels, since it is the highest order nonlinear effect (∼ I 5 ). For this value of the input chirp, t p = -150 fs, the plasma volume having absorbed the biggest possible energy is generated (note this is relatively easy to identify in experiments). By further negatively pre-chirping the input pulses to t p = -300 fs we obtain the longest possible plasma volume, shown in Fig. 5(a). This is physically understood from Figs. 6(a) and(c) showing several temporal profiles around the focal region, maxima of intensity, and fluence. The pulse reaches the focal volume with t p (z N L ) ≈ -150 fs (see the solid black temporal profile in Fig. 6(a), corresponding to z -z 0 = 0.15 cm in Fig. 6(c)). In this case, intensity is high enough for the Kerr effect to produce the refocusing cycles observed in the fluence, Fig. 6(c), that keep the on axis intensity relatively high, 10 13 W/cm 2 , for ∼ 1 mm. Indeed, refocusing cycles are linked to the local increase of intensity in the temporal profiles and to the long high intensity and high plasma density region. Note, interestingly, that the maximum plasma channel length does not coincide exactly with the t p minimizing transmission (t p ≈ -100, -150 fs), but, instead, the input pulse has to be slightly further negatively chirped (t p ≈ -300 fs). This is simply because in the latter case, the pulse still undergoes temporal compression on its way through the focus, what helps keeping higher intensities.

Transmission values in Fig. 3 rapidly increase towards the linear transmission limit for pulse widths t p -700, 45 fs. In these cases the plasma tracks are dramatically shortened down to ∼ 100 µm. However, on the negative chirp side, the maximum plasma density achieved presents an increasing trend as input duration increases up to, t p ∼ -2000 fs, as shown in Fig. 5(b). For longer input pulses, maximum of plasma density drops, as expected for very long pulses with finite energy. The fact that the maximum plasma density is located for rela-tively strongly chirped pulses (t p ≈ 2 ps) highlights the increasing importance of avalanche effects with t p (see Fig. 4(a)). Interestingly, such density increase observed in Fig. 5(b) is associated to the existence of a hump in intensity at t p = -1500, -2000, shown also in Fig 5(b). This hump is originated from small refocusing cycles that are manifested as intensity spikes after the first focusing event. This is caused by a pronounced and peaked trailing part of the pulse (after the leading one is attenuated) reaching the maximum intensity of the whole pulse propagation. This intense trailing part further accelerates electrons and is therefore able to further increase plasma density. Note that this feature of the temporal profile is shown in Fig. 6(b): the dashed black profile there is precisely the one corresponding to the intensity spike seen in the top part of Fig 6(d). This effect is possible only if the back part of the pulse is time-shifted from the front part by a delay smaller than recombination times. Outside the intensity hump, the maximum intensity is always reached by the leading part of the pulse. This explains the often experimentally reported observation that at a constant pulse energy, the longer the incident pulse is, the more permanent damage is produced in solids (see Ref. [35] for similar experiments in synthesized silica). When these long pulse reach the focus their low peak intensity has the double drawback that less plasma is excited and self-focusing is weaker, therefore high intensity region is shorter and less intense, as shown in Figs. 6(b) and 6(d).

To summarizing the results presented in this section, depending on the desired application, plasmas may be substantially long for the slightly negative chirped pulses at focus, or plasmas may be shortened and highly localized with higher densities for pulses pre-chirped up to a few ps.

V. RESULTS FOR λ0 = 400 nm

A motivation to study the effects presented in the above sections at the shorter λ 0 = 400 nm is that linear losses of water are close to their minimum for this wavelength (see β 1 in table I), so the focal point can be placed, in principle, at much longer distances from the laser source (∼ 50 times farther with same attenuation). However, in this case ionization will set up at lower intensities than for 800 nm (see Fig. 4(c)), providing seed electrons for avalanche effects at earlier stages relative to the focus. Since plasma cross section is almost constant for these wavelengths the total nonlinear losses, rdrdtdzW, will be enhanced (intensity levels were observed to stay relatively similar in both cases). This is in agreement with numerical findings, Fig. 7(a), where the minimum of transmission requires bigger pre-chirping than in Fig. 3 simply because normal GVD is stronger at 400 nm. The increase in the ionization rate is such that nonlinear losses are now dominated by MPA, as shown in Fig. 7(a). Intensity levels achieved here (∼ 10 13 W/cm 2 ) suggest that β K would not need to be significantly modified from the Keldysh value to match experimental results since the rate is dominated by MPA in this case (see Fig. 4(c)), oppositely to the 800 nm case.

Plasmas obtained here are shown in Fig. 7(b) and share all qualitative features with those observed at 800 nm, i.e., elongation close to the minimum of the optical transmission and shortening otherwise, suggesting these properties are rather general and wavelength independent.

VI. CONCLUSIONS

We showed that the control on focusing conditions and input pulse chirp provides a simple and effective mechanism to modify at will the electron-plasma density distribution generated by high energy (∼ 5 µJ) fs-pulses in water. Pulses reaching the focal volume with the shortest temporal profiles generate elongated plasma regions with relatively constant high densities. Plasma densities can be made higher by further pre-chirping the pulses (up to few ps), at the expense of the plasma channel length. Nonlinear losses, dominated by plasma absorption at 800 nm and by MPA at 400 nm, are found to be maximized for the cases in which MPA acquires its maximum efficiency. These remarks are independent of wavelength across the visible spectrum. The different plasma shapes indicate where laser energy is deposited by the pulse. This knowledge might be useful for developing applications of laser energy deposition to medical therapies and surgeries as well as to the controlled generation of sound waves in water.

VII. ACKNOWLEDGMENTS

Authors acknowledge financial support from the French Direction Générale de l'Armement (DGA). The need of a model that describes spatiotemporal propagation far from nonlinear focus is easily justified as follows. Under our experimental conditions, an order of magnitude estimate of the beam waist at focus, w f , may be obtained from initial (at position of the lens, z = -a) waist, w i = w F W HM / √ 2 ln 2, curvature, f i , and a, d (lens separation from water tank and focal length), with the laws of Gaussian (linear) optics [START_REF] Svelto | Principles of Lasers[END_REF]

f = ξ + Z 2 min /ξ, w 2 = w 2 min [1 + ξ 2 /Z 2 min ]
, where w min is the minimum waist and Z min ≡ πw 2 min /λ 0 the associated Rayleigh length in vacuum:

w 2 f = λ 0 2πn 0 Za 1 -1 - 4[d -a] 2 Z2 a , ( A1
)
Za ≡ n 0 Z a = n 0 πw 2 a λ 0 = n 0 Z f + {d -a} 2 Z f , Z f = dλ 0 πw 2 i f i , f i = Z 2 i 2d 1 -1 - 4d 2 Z 2 i , Z i ≡ πw 2 i λ 0 ,
where w a and Za are the beam waist and Rayleigh length at z = 0 + cm, i.e., just inside the water tank, respectively. In this situation, for the w F W HM = 7 mm beam we typically have w a /w f ≈ 10 3 and Z f ≈ 30 µm (numerical aperture N A ≡ w F W HM /[2n 0 d] ∼ 0.03). This implies that the computational grids discretizing r, z should contain ∼ 10 5 , 10 4 points, if the window widths are to be of the order 2∆r ∼ 5w a , 2∆z ∼ 10L w /Z f , to properly resolve the dynamics through the focal point, computation times taking up to weeks, exceeding any practical time scale, particularly in the frame of a parametric study demanding many simulation runs.

In our simulations, the typical input intensities at the entrance of the water tank, I ∼ 10 8 W/cm 2 , are well below ionization thresholds (J ≈ 0 in Eq. ( 1)). Pulse powers, however, are larger than the critical power for self-focusing, hence, the Kerr term dominates nonlinear effects until the increase of intensity due to beam focusing is such that typical Kerr and MPA lengths satisfy -2] ∼ 10 12 W/cm 2 , in our range of wavelengths. Below this intensity level, it is a reasonable approximation to assume that the propagation of the input beam (Ψ ≡ E/ ǫ 0 cn 0 /2: I = |Ψ| 2 ) Ψ(t, r, z = 0) = Ψ 0 exp -

L Kerr ≡ [n 2 k 0 I] -1 ≪ L MP A ≡ [2β K I K-1 ] -1 if I ≪ [n 2 k 0 / {2β K }] 1/[K
r 2 w 2 i 1 + i k 0 w 2 i 2f i - t 2 t 2 p0 [1 + iC] ,
(A2) can be described accurately by the method of pulse characteristics (see e.g., [46]), which describes the evolution in terms of the pulse moments. We define pulse energy, square beam radius and pulse duration as

FIG. 1 .

 1 FIG. 1. Sketch of the experimental setup. Pulses are incident from the left. Convergent lenses with focal lengths of 100, 75, 60 mm are represented by double arrows. The other elements are: power-meters (1, 2), beam splitter (BS), and water tank.

FIG. 2 .

 2 FIG. 2. (a) Transmission as a function of input energy for several wF W HM (values in legend). The dashed horizontal line marks the transmission in the linear limit, exp(-β1Lw). (b) Evolution across the nonlinear focus of (top) normalized energy, (center) maximum peak intensity (solid) and maximum plasma density (dashed), and (bottom) maximum beam waist (FWHM of the fluence), corresponding to the points A -D in (a). Vertical lines mark the linear focus, z f oc ≈ n0[da] ≈ 8.64 cm.

Figure 2 (FIG. 3 .

 23 Figure2(a) shows the optical transmission obtained for different focusing geometries, i.e., varying w F W HM from 3 to 7 mm and E in from 0.1 to 5 µJ (pulses are initially un-chirped: C = 0, t p = t pl = 45 fs). The drop in transmission when increasing either E in or w F W HM is due to the enhancement of multiphoton and/or plasma absorption. This happens because the localization of light imposed by the lens and further enhanced by self-focusing on the transverse plane, is more efficient than all other processes that yield decrease of intensity, such as pulse dispersion[47] and plasma defocusing[48,49].Figure2(b) shows the evolution of energy, peak intensity and plasma, and beam waist across the nonlinear focus for four different situations (A, B, C, and D in Fig.2(a)). Whilst low energy transfer (columns A, C) is associated with a short focal region of high intensity and a nearly symmetric diffraction from the focal plane, the situations in which high energy is transmitted to the

FIG. 4 .

 4 FIG. 4. (a) Plasma absorption to MPA ratio for several values of βK and τc as a function of input pulse width: (1) βK = 3.6× 10 -51 cm 7 /W 4 , τc = 3 fs; (2) βK = 6.2×10 -51 cm 7 /W 4 , τc = 3 fs; (3) βK = 3.6 × 10 -50 cm 7 /W 4 , τc = 3 fs; (4) βK = 8.3 × 10 -52 cm 7 /W 4 , τc = 10 fs. Vertical line marks the minimum pulse duration. (b) Energy transmission as a function of τc for fixed tp (βK = 8.3 × 10 -52 cm 7 W -4 ). Solid horizontal lines delimit the experimental interval of the measurement (see Fig. 3) and the dashed one the linear losses. (c) Keldysh ionization rate, W , (solid) versus I for water at λ0 = 800 and 400 nm. Dashed lines account for MPI only, neglecting tunnel ionization. (d) plasma absorption and defocusing dependence on collision time. Vertical line marks the maximum of σa.

FIG. 5 .FIG. 6 .

 56 FIG. 5. (a) Plasma volumes around the focus for several initial pulse widths (marked by the labels). Densities are given by ρ(r, z) = 10 x cm -3 , where the power x is given by the color bar. (b) absolute maxima of the electron density and intensity achieved in the focal volumes as a function of the chirped input pulse widths.

FIG. 7 .

 7 FIG. 7. (a) Transmission (left) and ratio of plasma absorption to MPA (right) as a function of the input pulse duration. Dashed vertical line marks the minimum pulse width (C = 0). (b) Plasma volumes for several input pulse widths (see labels) and electron densities ρ = 10 x cm -3 , where x is given by the color bar.

  Appendix A: Numerical model far from NL focus

TABLE I .

 I | ω0 are evaluated at Parameters used in Eqs. (1)-(3). All numerical results are produced with these values, with the exception of Fig. 4(a) in which several values of βK are used.

	magnitude	symbol (units)	value	Ref.
	vacuum wavelength	λ0 (nm)	800	400	-
	linear absorption	β1 (cm -1 )	0.0196	5 × 10 -4	[40]
	MPA order MPA (Keldysh)	K = U i ω0 + 1 βK (cm 2K-3 W 1-K ) 3.6 × 10 -50 5.4 × 10 -24 [41] 5 3 -
	MPA (this paper)	βK (cm 2K-3 W 1-K ) 8.3 × 10 -52 5.4 × 10 -24	-
	critical electron density	ρc (cm -3 )	1.7 × 10 21	7 × 10 21	-
	plasma absorption	σa (cm 2 )	6.3 × 10 -18 1.6 × 10 -18 Eq. (2)
	plasma defocusing	σd (cm 2 )	4.4 × 10 -17 2.2 × 10 -17 Eq. (2)
	linear refractive index	n0	1.33	[38]
	nonlinear index	n2 (cm 2 /W)	1.9 × 10 -16	[42]
	ionization potential	Ui (eV)	6.5		[43]
	e --ion collision time	τc (fs)	3		[43]
	recombination time	τr (fs)	100	[16]
	neutral atom density	ρnt (cm -3 )	6.7 × 10 22	[43]
	speed of light in vacuum	c (m/s)	299792458	-
	vacuum permittivity	ǫ0 (F m -1 )	8.85 × 10 -12	-

We use Eqs. (A6) to estimate the initial conditions for Eq. (1): w(z 0 ), f (z 0 ), and t p (z 0 ) inside the water tank for which I ≈ 10 11 W/cm 2 . For our input energies E in ≤ 5 µJ, this method lets us approach significantly the focal region, in the sense that w(z 0 ) ∼ 10w f .

where

Evolution equations for the above quantities are obtained through moment [45] or variational [44] methods that reduce the dimensionality of the problem by transforming Eq. ( 1) with J = 0 to (dots denote z-derivatives: ẋ ≡ dx/dz):

where P cr ≡ αλ 2 0 /πn 0 n 2 is the critical power for beam collapse (α = 4 in this analytical approach) and the pulse evolution is assumed to remain Gaussian:

Inserting Eq. (A5) into Eqs. (A4a), we eliminate the pulse intensity from the relation

which is to be solved with the initial conditions

p /t 2 pl -1] 1/2 , t pl = 45 fs being the width of the transform limited pulse as emitted by the laser source.


T Hoang
email: tien.hoang@enpc.fr

D Duhamel
email: denis.duhamel@enpc.fr

G Foret

H Yin

P Joyez

R Caby

DYNAMIC RESPONSE OF RAILWAY TRACKS IN TUNNEL

Keywords: Railway, Periodically Supported Beam, Dirac Comb, Fourier Transform

. More complex models can be constructed by including supports with damping or non-linear stiffness elements. This study deals with the dynamical modeling of non-ballasted railways, especially railways in tunnels. The model is developed as a dynamical system of multi-degree of freedom. Under the periodic assumption on the reaction force of the supports, the equation of motion for a periodically supported beam subjected to a moving load has been written. Then the Fourier transform has been used to solve this equation in case of damped supports. Analytical solutions have been established for the motion of the wheel and rail and also for the reaction force of the supports. The analytical solutions have been compared with in situ experimental measurements. The comparison shows that the theoretical results agree well with the measured results if damped supports are included in the model.

INTRODUCTION

In dynamic analysis of railway tracks, the analytical models of periodically supported beams have been often used to approximate the response of the system. In 2002, V.H. Nguyen [START_REF] Nguyen | Comportement dynamique de structures non-linaires soumises à des charges mobiles[END_REF] developed an analytical model by considering the railway as a continuously supported beam and by neglecting the mass of supports. In order to take account the vertical displacement of non-continued supports, X.Sheng and all. [START_REF] Sheng | Using the Fourier-series approach to study interactions between moving wheels and a periodically supported rail[END_REF] used Fourier series while G. Bonnet and R. Lassoused [START_REF] Bonnet | Response of a periodically supported beam to a moving load from the superposition of Bloch waves[END_REF] applied the superposition of Bloch waves. These methods are difficult to use if the mechanical behaviour of railway supports is more 1 T. Hoang, D. Duhamel, G. Foret, H. Yin, P. Joyez, R. Caby R n (t)

Q j Q 1 L x z D j v
Figure 1: Forces applied on a periodically supported beam complex (ex. hyperelastic or non-linear). In this paper, we will establish a general relation between displacements and forces of railway supports. This relation holds for periodically supported beams with any kind of supports and platform behaviour (ballast or nonballast). Then we will apply this relation to the case of linear supports to find out a response of the system and a validation with experiment results.

AN ANALYTICAL MODEL FOR PERIODICALLY SUPPORTED BEAM

In this study, we consider a general case of periodically supported beam. It means that the supports have the same mechanical constitutive law and they are distributed periodically, separated by a same length L as shown in figure 1.

The beam is subjected to moving forces Q j characterized by the distance to the first moving force D j (j = 1..K with K is the number of moving forces). The reaction force of support at the coordinate x = nL (with n = -∞..∞) is denoted by R n (t). The total force is calculated as following:

F (x, t) = ∞ n=-∞ R n (t)δ(x -nL) - K j=1 Q j δ(x + D j -vt) (1)
For the Euler-Bernoulli beam, the vertical displacement w r (x, t) of the beam under the force F (x, t) is solution of the following dynamic equation:

EI ∂ 4 w r (x, t) ∂x 4 + ρS ∂ 2 w r (x, t) ∂t 2 -F (x, t) = 0 (2)
where ρ and E are density and Young's modulus of the beam. The solution of equations ( 1) and (2) with initial conditions describe the dynamics of the beam with multi-degree of freedom. In general, these equations can not be solved analytically. However, we can find a general solution when the response of the system is stationary by using the following hypothesis:

Hypothesis: When the dynamic response of a periodically supported beam is stationary, the reaction forces of all supports are described by a same function but with a delay equal to the time for a moving load from a support to an other.

In other words, the reaction force repeats when a moving force passes from one to an other support: R n (t) = R(t -nL v ) where R(t) is the reaction force of the support at the origin. Thus:

F (x, t) = ∞ n=-∞ R(t - x v )δ(x -nL) - K j=1 Q j δ(x + D j -vt) (3)
By injecting this expression of the total force in equation ( 2) and by taking the Fourier transform, we obtain:

EI ∂ 4 ŵr (x, ω) ∂x 4 -ρSω 2 ŵr (x, ω) + K j=1 Q j v e -i ω v (x+D j ) -R(ω) ∞ n=-∞ e -i ω v x δ(x -nL) = 0 (4)
where ŵr (x, ω) and R(ω) are Fourier transforms of w r (x, t) and R(t) respectively. By taking the Fourier transform of ( 4) with regard to the variable x, we obtain:

(EIλ 4 -ρSω 2 )Π(λ, ω) + 2πδ(λ + ω v ) K j=1 Q j v e -i ω v D j -R(ω) ∞ n=-∞ e -i(λ+ ω v )nL = 0 (5)
where Π(λ, ω) is the Fourier transform of ŵr (x, ω) with regard to the variable x. The last term of ( 5) is a Dirac comb which has a following propriety:

∞ n=-∞ e -i(λ+ ω v )nL = 2π L ∞ n=-∞ δ λ + ω v + 2π L n (6)
Then from equation ( 5), the expression of Π(λ, ω) can be obtained:

Π(λ, ω) = 2π EI(λ 4 -k 4 b ) R(ω) L ∞ n=-∞ δ λ + ω v + 2π L n -δ λ + ω v K j=1 Q j v e -i ω v D j (7)
where k b = 4 ρSω 2 EI . Finally, the expression of ŵr (x, ω) is deduced by making the inverse Fourier transform of Π(λ, ω):

ŵr (x, ω) = R(ω) LEI ∞ n=-∞ e -i( ω v + 2πn L )x ω v + 2πn L 4 -k 4 b - K j=1 Q j e -i ω v (x+D j ) vEI ω v 4 -k 4 b (8)
For instance, the vertical displacement on the top of a support (x = 0) is:

ŵr (0, ω) = R(ω)η E (ω) - K j=1 Q j e -iω D j v vEI ω v 4 -k 4 b (9) R n (t), w r (nL, t) M 2 1 k 1 , η 1 k 2 , η 2 w n (t)
η E (ω) = 1 LEI ∞ n=-∞ 1 ω v + 2πn L 4 -k 4 b (10)
It is also possible to write η E (ω) in the following way:

η E (ω) = 1 4k 3 b EI sin(Lk b ) cos(Lk b ) -cos(L ω v ) - sinh(Lk b ) cosh(Lk b ) -cos(L ω v ) (11)
Equation ( 9) is a general relationship between the vertical displacement and the reaction force on the top of a support for a periodically supported beam under moving forces. The final solution of the displacement and the reaction force can be found by combining this relationship with the constitutive law of the support itself. In the next section, we will solve this equation in case of support with linear behaviour and apply it to a non-ballasted railway track.

APPLICATION IN CASE OF LINEAR BEHAVIOUR

In this section, supports with linear behavior will be considered. Damping will be taken into account and Kelvin-Voigt viscoelastic model will used as an example (see figure 2).

Let w n (t) denote the vertical displacement of block of the support n. The reaction force of support on the rail is given by:

R n (t) = -η 1 d(w r (nL, t) -w n (t)) dt -k 1 (w r (nL, t) -w n (t)) (12)
where η 1 , k 1 are the damping and spring coefficients of rail pad. The displacement of a block is governed by the following equation:

M d 2 w n (t) dt 2 + (η 1 + η 2 ) dw n (t) dt + (k 1 + k 2 )w n = η 1 dw r (nL, t) dt + k 1 w r (nL, t) (13)
where η 2 , k 2 are damping and spring coefficients of the equivalent elastic pad under the block and M is mass of the block.

Making the Fourier transform of the last equation leads to the following results:

ŵn (ω) = A(ω) A(ω) + B(ω)
ŵr (nL, ω)

Rn (ω) = -θ(ω) ŵr (nL, ω) (14)
where:

A(ω) = iωη 1 + k 1 B(ω) = -M ω 2 + iωη 2 + k 2 θ(ω) = A(ω)B(ω) A(ω) + B(ω) (15)
We see that the hypothesis R n (t) = R(t -nL v ) is equivalent to the following expression:

Rn (ω) = e -iω nL v R(ω) ∀n = -∞..∞ (16)
From equation ( 14), we show the last equation is satisfied if and only if:

ŵr (nL, ω) = e -iω nL v ŵr (0, ω) ∀n = -∞..∞ (17)
In the next section, we will prove that the response of the system under the hypothesis satisfies equation (17).

Response of railway track and verification of hypothesis

The following solution is deduced from equations ( 9) and (14) (with x = 0):

ŵr (0, ω) = -K j=1 Q j e -i ω v D j vEI [1 + θ(ω)η E (ω)] ω v 4 -k 4 b (18)
By injecting equations ( 9) and ( 14) into equation ( 8), we obtain:

ŵr (x, ω) = ŵr (0, ω)e -i ω v x 1 + θ(ω)η E (ω) - θ(ω) LEI ∞ n=-∞ e -i 2πn L x ω v + 2πn L 4 -k 4 b ( 19
)
By using equation (10), the following result can be obtained from the last equation:

ŵr (x, ω) = ŵr (0, ω)e -i ω v x 1 + θ(ω) LEI ∞ n=-∞ 1 -e -i 2πn L x ω v + 2πn L 4 -k 4 b ( 20
)
We see that equation (20) satisfies equation ( 17) when x = nL. Thus, the hypothesis is validated.

T. Hoang, D. Duhamel, G. Foret, H. Yin, P. Joyez, R. Caby

The motion of the rail can be calculated by taking the inverse Fourier transform of the last equation:

w r (x, t) = w r 0, x v -t + ∞ n=-∞ 1 -e -i 2πn L x 2πLEI ∞ -∞ θ(ω) ŵr (0, ω)e -iω x-vt v ω v + 2πn L 4 -k 4 b dω (21)
Then the reaction force and the displacement of the block can be calculated via equation ( 14). In a similar way, we could find solutions for railways with other linear supports.

Vertical vibration of wheels

From equation ( 21), the displacement of wheel -rail contact point at x = vt can be deduced:

w w (t) = w r (0, 0) + ∞ n=-∞ 1 -e -i 2πn L vt 2πvEI ∞ -∞ θ/LEI 1 + θη E K j=1 Q j e -iω D j v dω ω v 4 -k 4 b ω v + 2πn L 4 -k 4 b (22)
This displacement is a periodical motion of frequency f = v L described by the second term of equation ( 22). The amplitude A 0 of this motion can be obtained at the middle point between two supports (t = L 2v ).

A 0 = 1 2πvEI ∞ -∞ ∞ n=-∞ 1 -e -iπn ω v + 2πn L 4 -k 4 b K j=1 Q j e -iω D j v (1 + θη E ) ω v 4 -k 4 b dω (23)
Let ηE (ω) define by:

ηE (ω) = ∞ n=-∞ 1 -e -iπn ω v + 2πn L 4 -k 4 b
The amplitude A 0 of the motion of wheel can be reduced from equation (23) as following:

A 0 = 1 2πvEI ∞ -∞ θη E 1 + θη E K j=1 Q j e -iω D j v ω v 4 -k 4 b dω (24)
where ηE (ω) can be also written in the following form:

ηE (ω) = 1 4k 3 b EI sin( Lk b 2 ) cos( Lk b 2 ) + cos( Lω 2v ) - sinh( Lk b 2 ) cosh( Lk b 2 ) + cos( Lω 2v ) (25)

Approximation of the railway tracks response under the load of a train

The load by a train can be considered as a series of point load. The charges of each wagon on the train are often equal to the limit charge. It can be approximately by a series of identical charges (Q j = Q) at a distance from the first wheel D j = jH (for front wheels) and D j = jH + D (for back wheels) where H and D are the length of each wagon and distance between the font and the back wheel of a boogie. Using the propriety of Dirac comb, one can write:

∞ j=-∞ Q j e -iω D j v = 2πv H ∞ n=-∞ Qδ ω + 2πv H j (26)
By injecting the last equation into equation ( 18) and making the inverse Fourier transthe following analytical solution can be found for the response of the railway track:

w r (0, t) = -2QLk 1 k 2 H(k 1 + k 2 ) - j =0 e i2πj v H t    Q HEI 1 + e -iω D v ω v 4 -k 4 b [1 + θη E ]    ω=j 2πv H (27)
In the same way, we can get the displacement and the reaction force of a block by taking the inverse Fourier transform of ŵ0 (ω) and R(ω) which are deduced from equations ( 14) and (18) as following:

w 0 (t) = -2QL Hk 2 - j =0 e i2πj v H t    Q HEI 1 + e -iω D v A ω v 4 -k 4 b [A + B + ABη E ]    ω=j 2πv H (28) R(t) = 2QL H + j =0 e i2πj v H t    Q HEI 1 + e -iω D v θ ω v 4 -k 4 b [1 + θη E ]    ω=j 2πv H (29)
where A(ω), B(ω) in equation ( 15). These expressions show that the response of railway is a periodical motion with the frequency f = v H .

Example

In this section we will compare the analytical solution with experimental results for a non-ballasted railway. The transverse force and rotating motion are ignored. The experiments were conducted in 2005 in the Channel tunnel. The parameters used in the calculation are given in table 1.

The displacement of rail track is computed by using equation (21). Figure 3 shows this displacement of the track in the length of wagon H when the first wheel is on the top of

CONCLUSION

The railway response under moving trains has been studied by using a periodically supported beam model. First a general relationship between the rail track response, the support reaction forces and the moving loads has been established. Then the response of displacement and reaction forces have been calculated out for a railway with linear supports including damping. The analytical solution agrees well with the experimental results measured in the Channel tunnel. The model can be extended to Timoshenko beam to take into account shear forces. Extensions may also be made for including longitudinal forces or rotating motions of rail tracks and supports.

Figure 2 :

 2 Figure 2: Model of non-ballasted railway support where (1) is rail pad and (2) is support pad

T

  . Hoang, D. Duhamel, G. Foret, H. Yin, P. Joyez, R

Figure 4 :

 4 Figure 4: Displacement and reaction force at rail -support contact point

MN/m 18 a support. In other positions of the wheels, the displacement curves have similar form but they are below this curve (the amplitude of this change calculated by equation ( 24) is A 0 = 0.0065 mm). It is noticeable that the curve has a symmetrical form if damping is not included and the damping reduces the displacement at the first wheel position. The dynamic response of rail track was measured in 2005. Sensors were positioned under the rail at the rail-support contact point. The analytical solution and experimental results are compared in figure 4 and show good agreement for both displacement and reaction force if damping is included in the model.