Tiago Jose

Santos Moraes

Hailong Wu

Eric Semail

Ky Ngac

Duc Nguyen

Tan

Tiago José

Duc Tan Nguyen

Vu

Optimal torque/speed characteristics of a Five-Phase Synchronous Machine under Peak or RMS current control strategies

Keywords: «Multiphase drive», «Harmonic injection», «Torque optimization», «Control under constraints», «Traction drive»

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Multiphase PMSMs have been widely studied for transport applications because of their fault-tolerant ability [START_REF] Barrero | Recent Advances in the Design, Modeling, and Control of Multiphase Machines-Part I[END_REF]. One particular feature of the PMSM with (2k+1) phases is that vector control with high torque quality can be achieved even with non-sinusoidal back electromotive forces (back-EMF) in transient operation. In [START_REF] Barrero | Recent Advances in the Design, Modeling, and Control of Multiphase Machines-Part I[END_REF]- [START_REF] Wang | Optimum injected harmonics into magnet shape in multiphase surface-mounted PM machine for maximum output torque[END_REF] and [START_REF] Sui | Research on the vector control strategy of five-phase permanentmagnet synchronous machine based on third-harmonic current injection[END_REF], the first and third harmonics of currents contribute to the torque of a five-phase machine and of a nine-phase machine [START_REF] Slunjski | Control of a Symmetrical Nine-phase PMSM with Highly Non-Sinusoidal Back-Electromotive Force Using Third Harmonic Current Injection[END_REF]. It has been shown that the third harmonic current injection may lead to high torque density of a multiphase machine, depending on the harmonic distribution of the back-EMF obtained by machine design.

Besides, in order to estimate the maximum capabilities of a drive in steady states and transient operations, voltage and constraints must be considered during the control of multiphase PMSM. We will consider two current constraints for the optimization and control strategy: MTPA [START_REF] Scuiller | Maximum Reachable Torque, Power and Speed for Five-Phase SPM Machine With Low Armature Reaction[END_REF] [START_REF] Vu | Control strategies for non-sinusoidal multiphase PMSM drives in faulty modes under constraints on copper losses and peak phase voltage[END_REF] and MTPPC [START_REF] Feng | Dual three-phase PMSM torque modeling and maximum torque per peak current control through optimized harmonic current injection[END_REF] [START_REF] Vu | Torque optimization of seven-phase BLDC machines in normal and degraded modes with constraints on current and voltage[END_REF]. The first one maximizes the torque for a given RMS current, directly related to the copper losses, while the second strategy maximizes the torque for a given peak current, much important for the Voltage Source inverter sizing. When only the first harmonic of stator currents is used, the mathematical relationship between its RMS value and its peak value is fixed. Hence, both MTPA and MTPPC strategies result in the same torque and current values.

However, if the first and third harmonics of stator currents are considered, the mathematical relationship between the RMS value and the peak value can be changed by supplying different amplitudes and phases for each harmonic. Consequently, each strategy will result in different drive's behaviors.

These different strategies could have various influences on the control performance of multiphase PMSMs. However, the impact on the frontiers of the torque speed characteristics has not been studied. This paper aims to explain the impact of the MTPA and the MTPPC strategies on the maximum torque and copper losses of a 5-phase PMSM when the third current harmonic is injected.

In this paper, the improvement of torque density by adding the third current harmonics of a five-phase PMSM is introduced in the second section. Then the optimizing algorithms considering two constraints are illustrated in the third section. Section four gives and discusses the results. Section five summarizes the analysis of this paper.

Studied five-phase machine and inverter

An open-end winding (OEW) five-phase PMSM is designed and studied in this paper. This prototype has 20 stator slots and 14 poles as presented in Fig. 1(a). The rare-earth magnets are buried in the rotor in radial disposition called spoke. Besides, the inverter structure of this open-end winding machine is illustrated in Fig. 1

L M M M M M L M M M M M L M M M M M L M M M M M L         =         L ( 1 )
The cyclic inductances in dq rotating frame can be calculated by

C α α α α α α α α α α α α α α α α         =           (3) α=2π/5 cos sin 0 0 0 sin cos 0 0 0 = 0 0 c os3 sin3 0 0 0 sin3 c os3 0 θ θ θ     -     -         (4)
According to (2), the characteristics of the windings in dq rotating frame, obtained from inductances determined in stator natural frame by a finite element software, are listed in Table I. In [START_REF] Gong | Design Considerations of Five-Phase Machine with Double p/3p Polarity[END_REF], it has been shown that with MTPA strategy the optimal injection of a third harmonic of a machine leads to a torque.

= + = + - ( 5
)
Where p is the number of pole pairs, I1 and I3 are the peak amplitudes of the fundamental and the third stator current harmonics respectively; λm1 and λm3 are the peak amplitudes of the first and the third flux-linkage produced by permanent magnet respectively; θ1 and θ3 are the angles between the back-EMF harmonic and current harmonic respectively; r is the ratio between I1 and IRMS (r= I1/ [START_REF] Gong | Design Considerations of Five-Phase Machine with Double p/3p Polarity[END_REF]) and IRMS is a given RMS current value. When no constraints are imposed on current and voltage amplitudes, the optimal values for three parameters θ1, θ3 and r are calculated by dTe/dr=0.

1 1 3 3 2 2 1 1 3 3 3 1 1 3 m m opt m m E E r E E λ λ λ λ = =     + +         θ1=0 θ3= 0 ( 6 )
E1 and E3 are the first and the third harmonics of no load back-EMF respectively. In Fig. 2, the no load back-EMFs and the back-EMF/Speed harmonic distribution of the studied machine are presented. Their shapes are clearly non-sinusoidal. The amplitude of the third harmonic is close to 22.7% of the amplitude of the fundamental harmonic. Considering a 5-phase machine, it allows us to generate a constant meanvalue torque using the third harmonic. Therefore, according to the amplitude of the injected the third current harmonic should be 22% of the fundamental current harmonic.

Torque Optimization algorithm

Finite element method is an effective tool to analyze and evaluate the design of electrical machine. During the analysis stage, optimal stator currents should be obtained in advance and supplied into the machine. Therefore, a program developed on MATLAB is used to obtain the optimal current shape that generates the maximum torque considering MTPA and MTPPC strategies. The torque optimizing problem is stated as following. The objective function is ( )

1 1 3 3 1 1 3 3 1 3 ( , ,
For the RMS current constraint

I I   + ≤       ( 9 )
For the peak current constraint ( )

max 120 abcde I ≤ ( 1 0 )
The Flow chart of the optimization process is illustrated by Fig. 3. The parameters of this program are listed below:

The input inductance matrix L, resistance R and the first harmonic λ1 and the third harmonic λ3 of flux linkage without load are obtained by a no-load finite element simulation. The optimized parameters are the amplitudes (I1 and I3) and phases (θ1 and θ3) of the first and the third stator current harmonic. The optimizing algorithm is applied for the interested speed range to obtain the torque-speed characteristic. The objective function should be smaller than zero. If it is larger than zero, it means the machine is at generator mode.

The chosen constraints are associated to three components of the drive: the DC source, the Voltage Source Inverter and the machine. The value of the DC source will impact the maximum voltage which can be imposed to the machine. The OEW connection allows a maximum voltage per phase Uabcde equal to the voltage of the DC source.

The maximum RMS value, or the current density, which can be imposed in a machine is highly depending on thermal consideration. For example, we can consider that 71.4 Arms (6 A/mm²)) corresponds to steady-state operation (>30 minutes). As the thermal time constant constraints of transistor are less than 0.1s, the chosen peak value of the current Iabcde is associated to maximum allowable value 120A. If the MTPA strategy is applied, the corresponding stator currents can be optimized with the constraint [START_REF] Feng | Dual three-phase PMSM torque modeling and maximum torque per peak current control through optimized harmonic current injection[END_REF]. If the MTPPC strategy is considered, the corresponding optimal stator currents can be obtained with the constraint (10).

4.

Optimized results

a. Impact of the third harmonic injection for MTPA strategy

Firstly, the impact of the third harmonic injection is analyzed by the torque improvement algorithm with the current constraint [START_REF] Feng | Dual three-phase PMSM torque modeling and maximum torque per peak current control through optimized harmonic current injection[END_REF]. The first case (1h) is to only apply the fundamental current harmonic to optimized torque and the second case (1h+3h) injects the third current harmonic while the same maximum RMS value in both cases is considered. The results are presented in Fig. 4.

Comparing the curves of Fig. 4 (a), it is possible to conclude that the third harmonic current injection allows a torque increase for the same copper losses. The torque increase is quite small at low speed (2,4%). This low gain is expected because the amplitudes of the third harmonic of the current and of the back-EMF are equal to 22% of the fundamental ones. The obtained stator current (phase a) at low speed 750rpm is illustrated in Fig. 4 (b). The current of the second case (1h+3h) is no longer sinusoidal.

When the flux-weakening strategy is implemented, the torque increase overtakes 32% at 4500rpm as shown in Fig. 4 (a). The torque improvement at 5000rpm is over 100% because the available maximum speed is just 4800rpm for the first case (1h). It implies that injecting the third current harmonic can increase the flux-weakening region for the studied five-phase PMSM. The flux-weakening strategy is needed at high speed because the peak back-EMF exceeds the voltage limits of the drive. At high speed, over 2500rpm for this five-phase PMSM, the maximum peak voltage (60V) is a major constraint. The injection of the third harmonic increases from 2 to 4 the number of degrees of freedom (amplitudes and phases of first and third harmonics). By shifting the phase between the third and first harmonics, it is possible to significantly decrease the peak voltage. The resulting peak value of the voltage reference is less than or equal to the voltage limit. An example at 4000rpm is presented in Fig. 4 (c). The current shape is different from the sinusoidal stator current. The produced voltage and back-EMF are illustrated in Fig. 5. It is observed that back-EMF is 90V, but the voltage is reduced to 48V.

At last, the obtained optimal currents of these two cases are applied to the finite element model at speed 750rpm and 4000rpm. The torques are presented in Fig. 4 (a) by the red dots (case 1) and blue dots (case 2). At each speed, the torque is improved by the second case which is same as the conclusion summarized by analytical torque optimization algorithm. The difference is that the torques at 750rpm are smaller than the analytical results.

Because the inductances used in the algorithm are obtained at no load which means there is no saturation. But the finite element model has considered saturation.

b. MTPA and MTPPC strategies

In the previous section, the RMS current is used as a constraint for the machine characterization because this is directly related to the machine copper losses. However, when considering the whole drive, the peak-current is the limitation for the transistors. This is the reason why an analysis will be carried out by not limiting the RMS current but the peak one. This strategy is called MTPPC.

This analysis is only possible when the currents are composed by more than one harmonic, otherwise the peak and the RMS values are directly related. The chosen peak value for the MTPPC strategy is 120Â because it is the maximum peak-current value obtained with the MTPA strategy for the RMS current of 71.4Arms.

During the optimization, the peak current constraint ( 10) is used and the results are obtained. Fig. 6 and Fig. 7 present the maximum torque, RMS and peak current per speed curves for the both strategies. These figures confirm that the MTPPC strategy allows the drive to attaint higher values of torque (Fig. 6) without increasing the current peak value (Fig. 7 (b)). At low speed, below 1700rpm, it is possible to have 31% of torque gain for the same peak current, despite a loss increase of 38,6% (Fig. 7 (a)). Different from the comparison of the previous section, the major gain is not in the flux-weakening operation zone. Fig. 8 shows how different the current shapes are when comparing both strategies. The major difference is the phase shift of the third harmonic. Roughly speaking, the peak-current value when the MTPA strategy is applied is the sum of the first and third harmonics amplitude. This behavior is due to the shape of the back-EMF, intrinsic to the machine. Due to the back-EMF distribution, the third harmonic generates much less torque than the first harmonic. When applying the MTPPC strategy, the third harmonic of the current will be phase-shifted of π from the relative back-EMF harmonic, generating then a negative torque. However, it allows to have a higher first harmonic, consequently a much higher torque, without overtaking the peak-current limit. Fig. 8: Current shape of the phase a for both strategies at 750rpm. This behavior is interesting because, in applications such as automotive traction, the higher torque values are usually used during the accelerations. As the duration of those accelerations are relatively short, their impact on machine losses is limited. Finally, torque/speed characteristics of Fig. 6 define which is the most adaptable strategy for each functioning point.

Lastly, these two strategies are simulated at 750rpm and 4000rpm by finite element method. The torques are presented in Fig. 6 by the red dots (MTPPC) and blue dots (MTPA). The torques are larger when MTPPC is applied which is same as the analytical result. The saturation reduces the torques obtained by finite element method. It can be observed the torque is much reduced in MTPPC (12.9% at 750rpm and 4.23% at 4000rpm) than in MTPA (9.8% at 750rpm and 0.23% at 4000rpm). Because the peak currents are same for these strategies, but the RMS current are different. So higher RMS current produces more saturation.

Conclusion

The impact of the third harmonic injection on torque generation of a five-phase open-winding PMSM with two different strategies was studied. This third harmonic injection has been optimized in order to generate the maximum torque considering the current limitations of the transistors and the machine. At first, the injected third current harmonic could enlarge the flux-weakening region for the RMS current constraint (MTPA). Besides, the analysis of two constraints with injection of current harmonic allowed us to adapt the control of the drive, in order to increase the torque density of the drive without increasing its sizing. For the same peak current constraint, MTPPC could produce higher torque than MTPA. But the copper losses could also be larger because of higher RMS current. This control was highly relevant in an industrial context in which the sizing of the machine and the transistors of the inverter are defined together. A possible control for the automobile application was to apply MTPPC for short time such as acceleration and to use MTPA for steady state.

Fig. 1 :

 1 After analysis of the results of table I, it is will considered only three cyclic inductances L1= (Ld1+Lq1)/2= 177 µH, L3= (Ld3+Lq3)/2= 97 µH, Lh=159µH. It can be remarked that the value of the minimum time constant among L1/R, L3/R, Lh/R is equal to about 6.1ms. A PWM frequency of 10 kHz will ensure a good control of the average currents of the machine even for the zero-sequence current which is present in open-winding configuration. Finite element model of the machine (a) and drive scheme (b).

2 :

 2 Back-EMF at 750rpm (a) and back-EMF/Speed harmonic distribution (RMS value) (b).

Fig. 3

 3 Fig. 3 Torque improvement algorithm

Fig. 4 :

 4 Torque per speed and torque gain curves with and without third harmonic injection for a current of 71.4Arms (a), with examples of the current shape of the phase-a obtained at 750rpm (b) and 4000rpm (c).

Fig. 5 :

 5 Fig. 5: Current, Back-EMF and Reference Voltage of phase a with third harmonic injection at 4000rpm.

Fig. 6 :

 6 Fig. 6: Torque per speed curves for MTPA for a RMS current of 71.4Arms and MTPPC for a peak value of 120Â and torque gain percentage of MTPPC in comparison to MTPA.

Fig. 7 :

 7 Fig. 7: RMS current (a) and peak-current (b) per speed curves for MTPA for a RMS current of 71.4Arms and MTPPC for a peak value of 120Â.

  (b). The inductances in natural stator frame are obtained in stator natural frame by a finite element software:

		p	a b	a c	a d	a e
		ba	p	bc	bd	be
	abcde	ca	cb	p	cd	ce
		da	db	dc	p	de
		ea	eb	ec	ed	p

Table I : Winding characteristics.

 I

	Parameter	Value	Parameter	Value
		154.8 µH		199 µH
		96.5 µH		97.1 µH
		159.

Acknowledgment

This work has been achieved within the framework of CE2I project. CE2I is co-financed by European Union with the financial support of European Regional Development Fund (ERDF), French State and the French Region of Hauts-de France.


D Fontanari

C Bray

G Dhont

G Mouret

A Cuisset

F Hindle

R Bocquet

K M Hickson

Molecules probed with a slow chirped pulse excitation: analytical model of the free induction decay signal

which is well tailored to interpret the signals obtained with very fast chirped pulses, but is not sufficient to account for the signals in the case of slower chirped pulses used in spectroscopy to increase the signal to noise ratio. A theoretical study of the polarization of molecules in a cell, uniform supersonic jet or molecular beam and subjected to a chirped pulse is presented. Three degrees of approximation for the polarization are introduced and are compared with the numerical solution of the optical Bloch equations. These expressions enter the analytic expression of the free induction decay signal which is validated against experimental data on the rotational emission spectra of OCS molecules. A relation between the pulse duration, the line position in the chirped pulse and the signal amplitude is proposed in the thermalized case. It assists in the optimization of the chirped pulse parameters and in the estimation of the error associated with the line intensity.

I. INTRODUCTION

Chirped Pulse (CP) spectrometers are now widely employed at microwave and millimeter frequencies to perform high resolution spectroscopic studies [1][2][3][4][START_REF] Hindle | Infrared Millim[END_REF]. A small number of instruments have also been set up in the submillimeter band [6]. Such systems employ a two-step measurement cycle. Firstly, the gas being studied is polarized by a CP emitted from a powerful source. Secondly, after the source extinction, the emission of the molecules corresponding to the Free Induction Decay (FID) is recorded. The signal level is proportional to the polarization reached at the end of the CP. The spectrum of the molecular emission is then recovered by the use of a Fourier transform. The broadband nature of CP spectroscopy has the advantage that many rotational transitions may be probed simultaneously.

While this technique is very well suited to the frequency measurement of rotational transitions, measurements of line strengths require more care. It has already been shown that the position of the rotational transition frequency in the pulse can affect the measured intensity, especially in cases where pulse durations are increased to obtain sufficient signal to noise ratios [START_REF] Hindle | Infrared Millim[END_REF]7]. However, most groups use a model originally developed by McGurk et al. [8] to model FID signals. It neglects re-emission effects during the pulse, together with population and coherence relaxation times, due to a very fast passage on the transition. In the model, the polarization at the end of the pulse is proportional to T c /∆ω where T c is the pulse duration and ∆ω is the spectral extension of the CP. This behaviour is in contradiction with the experi- * robin.bocquet@univ-littoral.fr mental measurements for which we observe a maximum FID signal for a particular pulse duration [START_REF] Hindle | Infrared Millim[END_REF].

This work proposes new analytical expressions for the polarization which will be helpful to analyze the experiments with slower CP. Section II deals with the polarizing step where the molecules evolve under the influence of the CP. We propose three approximations and compare them with the numerical integration of the optical Bloch equations. Both thermalized cases (cell and uniform supersonic jet) and non-thermalized cases (molecular beam) are discussed. Section III presents the second step of a CP experiment, i.e. the Free Induction Decay. Results of Section II are leveraged to propose an analytical expression of the FID recorded signal. This formula is validated through comparison with experimental data in the thermalized case using an experimental set-up described in [START_REF] Hindle | Infrared Millim[END_REF], similar to those of [2,4,6,9]. We show how to use this expression to optimize the amplitude of the signal and to take into account the dependence of the intensity on the line position.

II. ANALYTICAL EXPRESSION OF THE POLARIZATION IN A CP EXPERIMENT

A. Optical Bloch equations

The molecular interaction with an electromagnetic wave has been extensively described in the literature [8,[START_REF] Mcgurk | Advances in chemical physics[END_REF][START_REF] Levenson | Quantum Electronics -Principles and Applications[END_REF][START_REF] Meystre | Elements of quantum optics 3rd edition[END_REF][START_REF] Park | [END_REF]. We discuss for simplicity a two-level isolated molecular system, given by its energy levels E a and E b with E a -E b = ω 0 . We consider the interaction between this two-level system and an electromagnetic field with angular frequency ω(t). The system of N molecules with a transition dipole moment µ ab is described by the evolution of its density matrix. Introducing the relaxation mechanisms by γ 1 = 1/T 1 and γ 2 = 1/T 2 (respectively, inverse of the decay time of the population difference and inverse of the dipole dephasing time) and W eq , the population difference at thermodynamic equilibrium, we get the generic form of the optical Bloch equations:

dz dt (t) = -(γ 2 + iω 0 )z(t) -iE(t)W (t), (1)
dW dt (t) = -γ 1 (W (t) -W eq ) + 1 2i (E * (t)z(t) -E(t)z * (t)), (2)
where W (t) and z(t) = P(t)/(N µ ab ) correspond to the difference of population and the pseudo-polarization respectively, P(t) being the polarization. The function E(t) = e -iα(t) Ω 0 represents the perturbation of the system with Ω 0 = µ ab E0 the Rabi frequency and α(t) a generic function of time. Equations ( 1) and ( 2) are solved for given functions W (t) and z(t) respectively. Using the general solution of a first order differential equation, we have:

z(t) = z(0)e -(iω0+γ2)t -i t 0 e -(iω0+γ2)(t-x) E(x)W (x)dx, (3)
W (t) = W eq + (W (0) -W eq )e -γ1t + 1 2i t 0 e -γ1(t-x) (E * (x)z(x) -E(x)z * (x))dx. (4)
We suppose that z(t) and W (t) may be written as series:

z(t) = z (0) (t) + z (1) (t) + z (2) (t) + • • • , W (t) = W (0) (t) + W (1) (t) + W (2) (t) + • • • .
(5) Replacing these series in the Eq. (3) and Eq.( 4) and identifying term by term, we find the series recurrence:

z (n+1) (t) = -i t 0 e -(iω0+γ2)(t-x) E(x)W (n) (x)dx, W (n+1) (t) = 1 2i t 0 e -γ1(t-x) (E * (x)z (n) (x) -E(x)z (n) * (x))dx. (6)

B. Mathematical formulation

We consider in the paper a linear CP with duration T c and spectral extension ∆ω. The CP frequency equals the molecular resonance at time r 0 T c with r 0 ∈ [0; 1], the relative line position in the CP (for r 0 = 1 2 the chirp is centered on the resonance ω 0 ). The angular frequency is given by:

ω(t) = dα dt (t) = ω 0 -r 0 ∆ω + ∆ω T c t. (7)
Except for free induction decay or double resonance experiments, we are dealing with z(0) = 0. Only the odd terms of the z-series and the even terms of the W -series are non-zero. In particular, the first non zero terms of the series are:

z (1) (t) = -i t 0 e -(iω0+γ2)(t-x) E(x)W (0) (x)dx, W (0) (t) = W eq + (W (0) -W eq )e -γ1t . (8)
A rotating wave approximation is performed by introducing ξ (n) (t) = e iα(t) z (n) (t). Defining β 2 = ∆ω Tc , the sweep speed of the pulse, the pseudo-polarization is then given by: ξ (1) (t) = W eq g(t, γ 2 ) + (W (0) -W eq )e -γ1t g(t, γ 2 -γ 1 ), (9) where

g(t, γ) = -iΩ 0 t 0 e -γ(t-x)-icβ 2 (t-x)Tc+i β 2 2 (t 2 -x 2 ) dx, (10)
which can be integrated using the error function:

g(t, γ) = ρ(t, γ)e iθ(t,γ) σ(t, γ), (11)
where

ρ(t, γ) = √ πΩ 0 √ 2β exp (-γ (t -r 0 T c )) , θ(t, γ) = (t -r 0 T c ) 2 β 2 2 - 3 4 π - γ 2 2β 2 , σ(t, γ) =erf e i π 4 √ 2 r 0 T c β -i γ β + erf e i π 4 √ 2 (t -r 0 T c ) β + i γ β . (12)
At time t = T c (polarization at the end of the pulse), we can simplify the expressions by ignoring small terms ( γ β 1 ) and taking a first order asymptotic expansion of the error function [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF] in order to obtain: g asy (T c , γ) = ρ asy (T c , γ)e iθasy(Tc) + g osc (T c , γ), [START_REF] Park | [END_REF] where

ρ asy (T c , γ) = √ 2πΩ 0 √ ∆ω T c exp [-T c (1 -r 0 ) γ] , θ asy (T c ) = (1 -r 0 ) 2 ∆ω 2 T c - 3 4 π, g osc (T c , γ) = Ω 0 ∆ω 1 1 -r 0 + 1 r 0 exp i ∆ω 2 T c (1 -2r 0 ) × exp [-T c (1 -r 0 ) γ] . (14)
Using Eq. ( 9), ξ (1) (T c ) is then approximated by: ξ (1) asy (T c ) =W eq g asy (T c , γ 2 )

+ (W (0) -W eq )e -γ1Tc g asy (T c , γ 2 -γ 1 ).

The term g osc (T c , γ) induces an oscillation whose amplitude decreases as T c and/or ∆ω increases. We neglect it to get an approximation of the pseudo-polarization:

ξ (1) app (T c ) =e iθasy(Tc) [W eq ρ asy (T c , γ 2 ) + (W (0) -W eq ) e -γ1Tc ρ asy (T c , γ 2 -γ 1 )]. (16)
Finally, the modulus of the pseudo-polarization is ap-

proximated by |ξ(T c )| ≈ ξ (1) (T c ) ≈ ξ (1) asy (T c ) ≈ ξ (1) app (T c ) , with | ξ (1) app (T c ) |= √ 2πΩ 0 T c ∆ω × e -Tc T 2 (1-r0) × W eq + (W (0) -W eq ) e -Tc T 1 r0 , ( 17
)
which is the product of three terms. The first one corresponds to the polarization in the McGurk approximation [8]. The second one describes an exponential decay with relaxation time T 2 , for a duration of T c (1 -r 0 ) (duration between the line position and the end of the pulse). The third term represents the relaxation of the difference population towards the equilibrium with a relaxation time T 1 , for a duration r 0 T c (duration between the beginning of the pulse and the line position).

C. Numerical validation of the approximations

Figure 1 compares the three approximations ( ξ (1)

(T c ) , ξ (1) asy (T c ) , ξ (1)
app (T c ) ) of the modulus of the polarization with the numerical integration of the optical Bloch equations using a Rosenbrock method [START_REF] Maple | Maplesoft[END_REF] for three different spectral extensions and three different line positions r 0 . We consider a two-level system consisting of the J = 16 and J = 17 rotational levels of the OCS molecule with W (0) = W eq and T 2 = T 1 = 0.3 µs. We notice that except for the case ∆ν = 10 MHz and r 0 = 0.9 of Fig. 1(c), the first approximation ξ (1) (T c ) given by Eq. ( 9) (green curve) reproduces the numerical integration (blue curve). As soon as ∆ν 100 MHz, the red curves corresponding to the second approximation ξ

(1) asy (T c ) (Eq. ( 15)) are acceptable. The black curves corresponding to the third approximation ξ (1) app (T c ) (Eq. ( 16)) have no more oscillations and reproduce the global behavior of the numerical results, in particular the optimal pulse duration and the corresponding value of the polarization. Figure 1 suggests however two remarks. First, the asymptote of the modulus of the polarization for high T c value is not zero. This effect is seen for r 0 0.5 and would be visible at T c > 10 µs for r 0 = 0.9. This asymptote decreases when ∆ν increases and/or r 0 decreases. Second, for ∆ν = 1000 MHz, the amplitude of the oscillations decreases from r 0 = 0.1 to r 0 = 0.5 but increases from r 0 = 0.5 to r 0 = 0.9.

D. Discussion

The molecules interacting with a chirped pulse can be probed under different experimental conditions: in an ordinary cell [START_REF] Hindle | Infrared Millim[END_REF][START_REF] Broquier | [END_REF]17], in a uniform supersonic jet [7], or in a molecular beam [2,3,18]. We can take advantage of the typical characteristics of each experimental setup summed up in Table I to adapt the analytic approximation (17) of the modulus of the polarization obtained in Section II B.

In an ordinary cell or in a uniform supersonic jet, the molecules are thermalized. The rotational temperature may be the ambient temperature in the former situation while it is lowered to a few dozens of kelvins in the latter one. In both cases, the relaxation time T 2 is of the same order of magnitude as T 1 and the initial population difference W (0) corresponds to its equilibrium value W eq . In these conditions, the analytic approximation of Eq. ( 17) is proportional to |W eq (T rot )| which can be significantly increased in a uniform supersonic jet with respect to an ordinary cell at 300 K by decreasing the temperature. The pulse duration giving the maximal polarization is found to be T c,max = T 2 / (2 (1 -r 0 )), which depends on T 2 but not on T 1 .

The behaviour of experiments in the non-thermalized case is different. In a molecular beam experiment, the rotational temperature is lowered to a few kelvins, the difference population W (0) can be as high as 50W eq . The analytic approximation of Eq. ( 17) simplifies to:

| ξ (1) app (T c ) |∝ T c ∆ω e -Tc T 2 (1-r0) |W (0)| e -Tc T 1 r0 .
The collisional linewidth is strongly reduced, implying generally T 2 , T 1 T c [START_REF] Park | [END_REF][START_REF] Steinfeld | Molecules and Radiation[END_REF]. Under these conditions, the modulus of polarization reduces to the McGurk's model [8]. The optimal pulse duration is then no longer visible. The molecular beam chirped pulse setup is very advantageous in term of signal to noise ratio. First, the signal is enhanced by the difference populations W (0). Second, the duration of the FID signal increases due to the relaxation time T 2 (see III) to such an extent that the decay of the FID signal is often driven by the Doppler inhomogeneous broadening (see Appendix A). The Doppler width is indeed proportional to ω 0 (see Eq. (A5)) and more pronounced in the millimeter or submillimeter ranges than in the microwave range. Unfortunately, the drawback of a molecular beam experiment is a limited operating frequency range due to the peaked shape of the population of energy levels towards the low J values [START_REF] Kroto | Molecular Rotation Spectra[END_REF].

III. FREE INDUCTION DECAY SIGNAL

A. Mathematical formulation

The source extinction occurs at t = T c . The emission of the polarized molecules which evolve freely is observed with time t = t -T c . The polarization in the rotating wave approximation P(t ) is given by (see Appendix A):

P(t ) = N µ ab ξ(T c ) exp - t T 2 - t 2 ∆ω 2 D 4 ln 2 , ( 18
)
where ∆ω D is the Doppler width. The electric field is emitted at angular frequency ω 0 , the molecular resonance, with an amplitude proportional to the polarization module [2]:

E(t ) ∝ P(t ) e -i(ω0t +Φ) + cc. ( 19
)
Chirped pulse experiments use an heterodyne detection scheme [START_REF] Hindle | Infrared Millim[END_REF] where a local oscillator of angular frequency ω LO is used to shift the spectrum from ω 0 to the intermediate angular frequency ω IF = ω 0 -ω LO . If φ is a phase term at the time origin, the output signal is (see Appendix B):

S(t ) ∝ N µ ab ξ (1) app (T c ) (20) ×e -t T 2 e -t 2 ∆ω 2 D 4 ln 2 × cos(ω IF t + φ).
The exponential terms describe dampings respectively due to the dephasing polarization and Doppler broadening. The last term is an oscillation at the intermediate angular frequency ω IF .

B. Experimental results

CP experiments have been performed in order to check the validity of the FID signal model of Eq. ( 20) in the case of a thermalized sample (W (0) = W eq ). We recorded the FID signal of the OCS molecule around 206.745 GHz corresponding to the rotational transition J = 17 → 16. For the line position r 0 as parameter, P = 100 µbar and ∆ν = 1000 MHz. Blue points are experimental data and black curve is the fit using Eq. ( 22). The green curve is a fit where the polarization used in Eq. ( 18) is given by Eq. ( 9) rather than Eq. ( 17 a given pressure (i.e. T 2 fixed) and a given relative position r 0 , we measured the FID signal for different pulse durations T c . After a Fast Fourier Transform, the amplitude at the intermediate frequency ω IF is proportional to:

F (ω = ω IF , T c ) = N T c ∆ω e -Tc T 2 (1-r0) . (21)
These amplitudes F (ω = ω IF , T c ) measured for different T c values have been compared with two models. The simplest is given by:

B P T c ∆ω exp [-T c (1 -r 0 )2πγ OCS P ] , (22)
which is Eq. ( 21) reformulated with experimental parameters. In the second model, the polarization used in Eq. ( 18) is given by Eq. ( 9) rather than Eq. ( 17), corresponding to ξ (1) (T c ) . The number of molecules N is assumed proportional to the gas pressure

P , T 2 = 1 2πγOCSP
is determined by the OCS self pressure broadening coefficient γ OCS = 4.8 GHz/atm of the J = 17 → 16 transition, measured in Ref. [START_REF] Matton | [END_REF]. The parameter B corresponds to the instrument sensitivity of a given experimental configuration. It accounts for variations in mixer conversion losses and amplifier gains. The value of B is constant for a given experimental configuration. This parameter and the pressure are fitted against the experimental data sets in the second model. These two values are also used to plot the black curve corresponding to the simplest model. Figure 2 superposes the experimental data (blue points) with the fitted model of Eq. ( 22) (black curves) and the model corresponding to ξ (1) (T c ) (green curves). A good agreement between the experimental and fitted data is observed for the three different line positions r 0 at P = 100 µbar (corresponding to T 1 ≈ T 2 = 0.3 µs) and ∆ν = 1000 MHz. For all experimental curves, the relative difference between measured and fitted amplitudes with the simplest model never exceeds 30%. Such differences, especially around T c ≈ T c,max for r 0 = 0.1 and 0.5, are attributed to the oscillations mentioned above in Section II C and seen again on the green curves. For r 0 = 0.9 the oscillations are present over the entire T c range. Oscillations have also been discussed as "edge effects" by Park et al. [22] for short pulse durations due to the "windowing" of the perturbating electric field. This effect is reduced in our experiment by programming raised cosine edges (5 ns rise time) rather than a rectangular time window in the Arbitrary Wave Generator. The edge effects might explain the underestimated amplitude of the oscillations, given by our model (green curves), visible for r 0 = 0.5 at low T c values. The signal does not fall to zero for large T c values which is consistent with Fig. 1, where the first approximation coincides exactly with the numerical simulation for ∆ν = 1000 MHz. The amplitude of the experimental data in this region is nonetheless higher than the fitted curve. This might suggest a shortcoming in the isolated two-level model itself. However, the results exhibited in Fig. 2 demonstrate that our simplest model with its approximations is sufficient to reproduce the behaviour of the emitted intensities in the chirped pulse experiment. This conclusion is supported by the analyses of experimental data at P = 20 µbar and P = 50 µbar which show strictly identical patterns.

C. Discussion

1. Optimization of the CP parameters at a given pressure Eq. [START_REF] Kroto | Molecular Rotation Spectra[END_REF] in the case of a thermalized sample is relatively simple but depends on parameters T c , r 0 , T 2 and N . The signal amplitude at t = 0 can be rewritten in the form of Eq. ( 22) to show the influence of the gas pressure or equivalently in the following form to show the influence of T c /T 2 :

S r0 Tc/T2 ∝ 1 √ T 2 ∆ω T c T 2 e -Tc T 2
(1-r0) .

Experimentally, T 2 (adjusted with the gas pressure) is fixed by considerations of Fourier transform resolution. For a given T 2 , the optimal pulse duration is given by T c,max = T2 2(1-r0) . Figure 3 is a three dimensional view of the amplitude of the signal from Eq. ( 23), versus T c /T 2 and r 0 at T 2 fixed. The dashed line gives the signal amplitude of S 0 2 . It shows that as long as Tc T2 < 0.2, the signal amplitude is quasi independent of the line position r 0 . Inside this region, the output signal is proportional to T c /∆ω and corresponds to the approximation of a fast chirped pulse obtained by McGurk et al. [8]. Outside this region, the signal amplitude can be increased significantly but at the expense of a dependence on the line position r 0 within the pulse. However, we observe that up to Tc T2 = 2, the amplitude of the signal for r 0 = 0 is always greater than the amplitude for

T # c T2 ≈ 0.04. In other words, working at T c /T 2 = 2 is more interesting than T c /T 2 T # c
T2 in term of signal to noise ratio: the gain in amplitude for r 0 = 1 is an order of magnitude while the gain for r 0 = 0 is at least one.

Dependence of the line strength on the line position

Abeysekera et al. experimentally discovered the line position effect on the line intensity [7]. They considered two related experiments to mitigate the problem: in the first experiment, the frequencies are swept to higher values, whereas the frequencies are swept to lower values in the second experiment.. They proposed to average the two CP spectra to compensate for the line position effect: S ave = (S up + S down ) /2. They worked with a uniform supersonic jet for which Eq. ( 23 r0 where A is a proportionality constant. The average of the two signals,

S ave = Ae -Tc 2T 2 cosh T c T 2 r 0 - 1 2 ,
is quasi-independent of the line position r 0 if the argument of the cosh is small, i.e. for a chirp duration shorter than T 2 . The series expansion of the hyperbolic cosine gives:

S ave = Ae -Tc 2T 2 1 + 1 2 T c T 2 r 0 - 1 2 2 + ... .
indicating that the error done is always lower than 12.5% if T c /T 2 1.

IV. CONCLUSION

The FID signal of a two-level system interacting with a linear CP depends on the polarization at the end of the CP. The value of this polarization, either in thermalized cases (cell or uniform supersonic jet) or non-thermalized cases (molecular beam), was approached through three consecutive analytical approximations which improve the model proposed by McGurk [8]. In particular, the last approximation neglected the oscillations in the polarization and a very simple formula for the optimal pulse duration (for which the modulus of the polarization is maximum) was obtained.

The simple expression between the pulse duration, the line position in the CP and the signal amplitude will certainly be useful in spectroscopy, where it is convenient to use slower CPs to increase the signal over a fast CP. Moreover, this relation describing the connection between the position of the lines in the CP and their intensity will allow for a correction of the intensities required for high resolution molecular spectroscopic studies.

∆ω D = ω 0 c 2RT ln 2 M , ( A5
)
where T is the temperature and M is the molecular weight. Different molecules have different frequency resonances due to the Doppler effect and thus different detuning δ. The polarization is given by the integral of Eq. (A4) over all the detunings: (A6)

U (t ) V (t ) =     +∞ -∞ e -

Figure 1 :

 1 Figure 1: Polarization modulus at the end of the CP against the chirp duration T c for three spectral extensions ∆ν ∈ {10, 100, 1000} MHz and three line positions r 0 ∈ {0.1, 0.5, 0.9}. Blue curves: numerical integration of optical Bloch equations. Green curves: first non zero term ξ (1) (T c ) of the series. Red curves: asymptotic expansion ξ (1) asy (T c ) of the polarization. Black curves: polarization ξ (1) app (T c ) after the last approximation. Parameters are: J = 17 → 16 rotational transition of OCS at 206.745 GHz, Ω 0 = 2.4 rd/µs, W (0) = W eq = 6 10 -4 , T 2 = T 1 = 0, 3 µs.

Figure 2 :

 2 Figure 2: Comparison between experiment and models,the line position r 0 as parameter, P = 100 µbar and ∆ν = 1000 MHz. Blue points are experimental data and black curve is the fit using Eq. (22). The green curve is a fit where the polarization used in Eq. (18) is given by Eq. (9) rather than Eq. (17). The red vertical dotted lines indicate the positions of T c,max . (a) r 0 = 0.1. (b) r 0 = 0.5. (c) r 0 = 0.9.

  Figure 2: Comparison between experiment and models,the line position r 0 as parameter, P = 100 µbar and ∆ν = 1000 MHz. Blue points are experimental data and black curve is the fit using Eq. (22). The green curve is a fit where the polarization used in Eq. (18) is given by Eq. (9) rather than Eq. (17). The red vertical dotted lines indicate the positions of T c,max . (a) r 0 = 0.1. (b) r 0 = 0.5. (c) r 0 = 0.9.

Figure 3 :

 3 Figure 3: Three dimensional view of the amplitude of the signal S in arbitrary unit from Eq. (23) versus T c /T 2 and r 0 at T 2 fixed.

2 ( 1

 21 ) applies. The signals are the same in the two experiments except that r 0 in the first one is replaced by 1 -r 0 in the second one. The corresponding signals are respectively S up = A e -Tc T -r0) and S down = A e -Tc T 2

t T 2

 2 [cos (δt ) U 0 -V 0 sin (δt )] ln 2 π e ln 2 δ 2 ∆ω 2 D ∆ω D dδ +∞ -∞ e -t T 2 [sin (δt ) U 0 + V 0 cos (δt )] ln 2

ACKNOWLEDGMENTS

We would like to acknowledge the French Agence Nationale de la Recherche for their financial support of this work via funding of the project Original Submillimetre Chirped pulse instrumentation for Astrochemical Reactivity (OSCAR) under contract number ANR-15-CE29-0017. The authors are also grateful for the financial support received from the Région Hauts-de-France, the Ministère de l'Enseignement Supérieur et de la Recherche (CPER Climibio), and the European Fund for Regional Economic Development through the TER-AFOOD project (INTERREG V FR-WA-VL 1.2.11).

Appendix A: Doppler broadening in the Free Induction Decay

The FID corresponds to the extinction of the source at t = T c and the reemission of the polarized molecules which freely evolve. We perform a change of variable in this section: t = t -T c .

System evolution

Since the source is switched off, the Rabi frequency is set to zero and the optical Bloch equations in the rotating frame [START_REF] Meystre | Elements of quantum optics 3rd edition[END_REF]23] reduce to:

where δ is the detuning from the resonance, null if the inhomogenous broadening is negligible and W eq is the population difference at the thermal equilibrium. (U, V, W ) is the Bloch vector and the polarization is

] N µ ab . This is a simple coupled system of differential equations. The third equation is directly integrable giving:

) where W (T c ) = W (t = 0). The remaining part of the system is then:

which may be diagonalizable to find eigenvalues and eigenvectors and finally the solution:

, (A4) where U 0 = U (t = 0) and V 0 = V (t = 0).

Inhomogeneous broadening

We must take into account all the contributions due to inhomogeneous broadening characterized by ∆ω D , the Half Width at Half Maximum (HWHM) given by: These integrals are tabulated or can be computed with Maple [START_REF] Maple | Maplesoft[END_REF] giving:

Appendix B: Detection of the FID signal

In a CP experiment we are interested by the electric field reemitted during the FID sequence with the polarization given by equation (A7). We calculate the electric field in the Slowly Varying Amplitude and Phase (SVAP) approximation [8,[START_REF] Meystre | Elements of quantum optics 3rd edition[END_REF]. The general forms for electric field and polarization that propagates along the y axis are:

The SVAP approximation neglects the amplitude and phase variations during temporal and spatial periods:

If we neglect also 1 c ∂E ∂t compared to ∂E ∂y , we get simplified propagation equations for the real and imaginary parts:

equations that we can easily integrate between y = 0 and y = L corresponding to the length of the gas cell:

Finally, taking Eq. (A7) into account, we get:

In the experiment, we access the electric field by an heterodyne scheme tuned in the way to observe S(t ) = A E(t ) cos(ν IF t + φ) where A is a constant depending on mixer efficiency, amplifier gain etc.., ν IF is the intermediate frequency (ν IF = ν 0 -ν LO ) and φ the phase at the origin of time.

S(t ) ∝ |P(T

Table I: Characteristic parameters for three experimental setups. The simplified modulus of polarization (see text) and the optimal duration T c,max for the maximum signal are reported. T rot is the molecular rotational temperature.


H Peng

C Riconda

M Grech

C.-T Zhou

S Weber

Dynamical aspects of plasma gratings driven by a static ponderomotive potential

Keywords: plasma grating, laser manipulation, nonlinear plasma dynamics, laser-plasma interaction

A quasi neutral plasma grating is considered which is generated by two counter-propagating identical laser beams. Typical times scales associated with such grating are given by tunit = mi/2Zme(kv0) -1 ,where v0 is the electron quiver in the laser field of wavevector k. In most situations the behavior of the grating can be characterized by a single parameter, µ, which is proportional to the ratio of background electron temperature to the square of the electron quiver velocity. Indeed even for quasi-neutral gratings the electron pressure might play an important role, with adiabatic electron heating. The influence on grating formation and lifetime of other parameters that imply the inclusion of dissipative effects or kinetic effects is also examined in detail. Finally, an approximated analytical solution to the fluid model is found that shows good agreement with first-principle simulations.

I. INTRODUCTION

In the last few years there have been worldwide efforts to construct multi-petawatt laser facilities [1][2][START_REF] Gaul | Design of the texas petawatt laser[END_REF][START_REF] Yu | [END_REF][START_REF] Lureau | 10 petawatt lasers for extreme light applications[END_REF]. These kinds of power levels prove to be very demanding on optical devices. The intrinsic fluence limitations determine that the optical elements can only sustain ultrashort laser of few hundred mJ/cm 2 [START_REF] Uteza | [END_REF][7][8][9]. Thus large-scale optical elements are required which are difficult to manufacture and very expensive.

By contrast, plasmas can withstand very high energy fluxes, up to 6 orders of magnitude higher than solid-state based materials. Also a plasma supports many linear and nonlinear effects which can be employed to control, manipulate and shape ultra-intense laser pulses [10][11][12][13][14][15][16][17][18][19]. Plasma optics, although still in its infancy, has therefore attracted a lot of attention in recent years and many successful proof-of-principle experiments have been carried out. A plasma grating is a plasma device that can work as photonic crystal, capable of compressing and stretching chirped pulses, creating optical solitons [20,21], acting as an optical switch based on the laser frequency or injection angle [19]. A plasma grating can also change the polarization and phase of a laser pulse. They have been used as laser polarizer and waveplate for both, long and moderately intense laser pulses [16][17][18], as well as short and ultraintense laser pulses [22]. On the surface of a solid target, such a generated grating acts as a plasma hologram for ultraintense lasers [23]. Plasma gratings can also be generated in underdense plasma as volume hologram for laser focusing and mode conversion [24].

Moreover, it can be used to couple the laser into surface plasma waves instead of solid grating [25].

Forslund et al. first observed that in the late phase of strongly-coupled stimulated Brillouin scattering, a large amplitude, strongly-modulated and quasi-neutral plasma grating is generated [26,27]. The formation of the plasma grating is explained by the steepening of the ion velocity profile and the accumulation of particles in the troughs of the ponderomotive potential created by the beating of counter-propagating laser beams. Its collapse was explained by X-type breaking of the grating [26][27][28]. Forslund et al. studied the behavior of the ion acoustic waves with initial sinusoidal velocity distribution, and compared the ion kinetic energy with the maximum possible thermal potential, to explain the X-type breaking [29]. However, neither the dynamical effect of the thermal pressure nor the effect of a source term was included. Sheng et al. proposed a linearized fluid model to explain the initial phase of the formation of the plasma grating in the ponderomotive potential of two counter-propagating lasers [30], which quickly fails when the nonlinearity becomes important.

In the present work, we propose a nonlinear model of a plasma grating growing in the static ponderomotive potential generated by two identical counter-propagating lasers. In the considered regime, electrons are localized in the plasma grating and adiabatic electron heating should be considered. This configuration and the model we propose is adapted to give correct predictions for pulses on the picosecond/femtosecond scale with relatively high intensity. These high intensities allow to reach high peak values of the grating close to or even higher than n c , necessary for many applications. The ratio of the adiabatic thermal potential to the ponderomotive potential (denoted by µ) is crucial to determine the dynamics of the plasma grating [31]. In this paper, the details of the growth, transient state and collapse are considered us-ing both fluid and PIC (particle-in-cell) simulations, a detailed comparison is also performed. Good agreement is found between the two methods within the limits of the fluid approximation. Particular emphasis are put on investigating kinetic effects with PIC simulations in the long timescales, which is beyond the limits of the fluid approximation. The role of parameters affecting the dynamics of the gratings such as plasma density, ion temperature and particle collisions are also discussed.

This paper is organized as follows: in Sec.II, the equations describing the model are presented. Sec.II A presents the macroscopic approach and the following section, Sec.II B, compares the results of macroscopic and microscopic simulation. In Sec.III, depending on the value of a characteristic parameter, µ, three regimes are found, in which the plasma grating displays very different behavior: complete reflection, partially reflecting and partially crossing, and crossing. In Sec.IV, the various time scales of importance to the dynamics are briefly discussed. The following section, Sec.V, it is discussed how the use of the interaction parameters allow to customize gratings for experiments. Sec.VI discusses how kinetic effects such as a finite ion temperature and particle collisions affect the grating dynamics. Finally, in Sec.VII, conclusions are presented.

II. MODELLING THE GRATING DYNAMICS

The details of the growth, transient state and collapse are considered using both fluid and PIC (particle-in-cell) simulations, which describe the grating dynamics from the macroscopic and the microscopic point of view, respectively.

As is shown subsequently, the two methods agree very well as far as the essential aspects of grating formation are concerned. The macroscopic fluid approach is more intuitive and allows better to understand simple scaling and establish governing parameters. The kinetic (PIC) approach is a first-principle-approach which includes by definition all physically relevant processes, in particular such processes as wavebreaking which can not be described in the fluid picture. The kinetic approach therefore allows to follow the dynamics all the way.

As far as a full understanding is concerned the two approaches can be considered complementarily and are both needed for a full understanding. In the following the macroscopic equations are presented in Sec.II A and the comparison with kinetic simulations is presented in Sec.II B.

A. Governing macroscopic equations

The nonlinear two-fluid equations including the ponderomotive potential and neglecting the electron momen-tum are:

0 = e ∂φ ∂x - 1 n e ∂p ∂x -e ∂φ p ∂x , (1a)
∂ 2 φ ∂x 2 = 4πe(Zn e -n i ), ( 1b
)
∂n i ∂t + ∂(n i v i ) ∂x = 0, ( 1c
)
∂v i ∂t + v i ∂v i ∂x = - Ze m i ∂φ ∂x . (1d)
The ponderomotive potential

φ p = 1 2 mec 2
e a 2 0 cos(2kx) is generated by two identical counter-propagating lasers in the plasma. Here k is the laser wave vector in the plasma and a 0 is the amplitude of lasers. k = k 0 1n 0 /n c , where k 0 , n 0 are the laser wavevector in vacuum and unperturbed plasma density, respectively. Adiabatic heating is considered. Thus:

pn -γ e = const, (2a) p = n e T e . (2b)
γ = (2 + D)/D, D are the adiabatic index and number of freedom, respectively. The grating is one-dimensional, so D = 1 and γ = 3. Then the thermal potential is

φ th = 3 2 T0 e ( ne n0 ) 2
, where T 0 is the initial electron temperature. Then Eq(1a) leads to:

φ = φ p + φ th . (3)
The ponderomotive force is f p = -e∂φ p /∂x. One can further normalized these equations with:

x unit = 1 2k , t unit = 1 2 m i Zm e (kv 0 ) -1 , (4a)
v unit = 1 2 Zm e m i v 0 , n unit = n 0 . (4b)
Then Eq(1) can be transformed to:

µ ∂ 2 n 2 e ∂x 2 -cos x = ν(Zn e -n i ), ( 5a
)
∂n i ∂t + ∂(n i v i ) ∂x = 0, ( 5b
)
∂v i ∂t + v i ∂v i ∂x = sin x -µ ∂n 2 e ∂x . ( 5c
)
where µ = 3T 0 /m e v 2 0 , ν = (ω p0 / √ 2kv 0 ) 2 , v 0 is the electron oscillation velocity for nonrelativistic laser intensities and ω p0 is the initial electron plasma frequency. Note that ν < 1 is used as a criterion of the superradiant regime [32], where the electron dynamics are very important and the electron momentum cannot be neglected. In this regime, the grating dynamics are much faster(usually on the order of tens of femtosecond). While we are focusing on ion-dominant plasma grating which is slowly evolving and quasistationary. This kind of plasma grating is suitable for manipulation of laser pulses with duration going up to 100 fs or even few ps. So this regime(ν < 1) should be avoided in the discussion here. Note also that the strong coupling stimulated Brillouin scattering takes place when (v 0 /v th ) 2 > 4k 0 c s ω p /ω 2 p0 [33], which corresponds to µ < µ th = [2a 2 0 me mi ( nc n0 ) 2 (1 -n0 nc )] -1/3 . One finds that in the case ν 1, Eq(5a) can be simply replaced with Zn e = n i .

B. Macroscopic and microscopic comparison of the process

To validate the fluid model formulated above, we first numerically solve it and compare the results to two PIC simulations with the same µ and ν: µ = 1.5, ν = 59.5(Z = 1 is used throughout the paper unless specified). The simulation results of the fluid model are shown in Fig( 1). Because the ponderomotive potential is static and periodic along the x axis, only one period is shown and periodic boundary conditions are applied. There are obviously three different stages of the plasma grating, which are shown in (a), (b), (c), respectively.

At first, the ion velocity grows and steepens under the drive of the ponderomotive potential, leading to the gathering of ions and electrons towards the potential trough. But the thermal potential φ th = 3 2 T0 e ( ne n0 ) 2 also grows quickly as the electron density grows quickly in the grating. The combined potential of φ p and φ th has a plateau, which then stop the ions moving towards the grating. A velocity plateau is established.

Then in the second stage, as the ion kinetic energy is not big enough for the ions to climb up the potential barrier and pass it(P < e∆φ, P is the kinetic energy of the ions), the ions keep accumulating at each edge of the grating. And the grating grows "two spikes" of ions. The following ions are completely reflected, leading to a X-like ion phase space that can only be observed in PIC simulations(see following discussion on PIC simulation results). This corresponds to the X-type breaking that has been widely discussed [27][28][29]34].

In the last stage, the reflected ions lead the grating to expand. The ions co-propagating with the edges of the potential barrier are accelerated to high velocity. The plasma density in the grating starts to decrease as the grating gets stretched. As more and more ions accumulates at each edge, the ion density is deviated from the electron density significantly at each edge. 1), respectively. tunit 320 fs here. From upper row to lower row are the ion phase space, electron phase space and local electron temperature(green line, ticked on the right), ion density(red line) and electron density(blue line) at different times in the center of the simulation box, respectively. The ion velocity and plasma density in normalized units are ticked in the right. Note that a Savitzky-Golay filter is applied to smooth the electron temperature to increase the precision without distorting the trend.

The temporal evolution of the plasma gratings in the center of the simulation box obtained from the PIC simulations are shown in Fig (2) and Fig( 3). In the first PIC simulation showing in Fig( 2), the laser wavelength, laser amplitudes, electron temperature, unperturbed plasma density are λ 0 = 1 µm, I = 5 × 10 15 W/cm 2 (a 0 0.06), T e = 0.92 keV, n 0 = 0.3n c , respectively. In the second PIC simulation showing in Fig( 3), the laser wavelength, laser amplitudes, electron temperature, unperturbed plasma density are λ 0 = 1 µm, I = 1 × 10 16 W/cm 2 (a 0 0.0855), T e = 1.84 keV, n 0 = 0.46n c , respectively. In both PIC simulations, µ = 1.5 and ν = 59.5. All the lasers have slowly linearly growing fronts of 10πω -1 p0 to prevent the triggering of the Langmuir oscillation of the electrons and the remainder of the laser pulse is constant. Plasma densities above quarter-critical are used throughout the paper, in order to avoid the noise-induced Raman scattering and the interference of the resultant Raman heating. The ion temperatures in both PIC simulations are T i = 1 eV. Note that here the strong coupling threshold is µ th 32. Thus we are in the strong coupling regime, which means at the beginning the thermal potential is much smaller than the ponderomotive potential and can be neglected. The one-dimensional SMILEI code [35] is used for the PIC simulations. Cell size of λ 0 /256, 50 particles per cell, m i /m e = 1836 and Z = 1 are used. We use a 6λ 0 plasma plateau with 2λ 0 vacuum at each side of the simulation boxes. The same three stages are found easily in the PIC simulations. They agree well with the fluid simulation results. The electron temperature in the grating does increase as the electrons are compressed in the grating from the initial 920 eV to about 2 keV in Fig( 2), from 1.84 keV to 3.08 keV in Fig( 3), and decrease as the grating gets stretched. Fig( d) of both figures show the start of another round of formation of plasma grating. The fluid simulation fails at this time, because that the ion velocity is double-valued in the spikes. The ion spikes collide with the counter-propagating ones and coalesce into a bigger one and separate from each other after that. The plasma gratings are compressed to a higher amplitude in the second round. The possible reason is that in the first round, the lasers have a slow rising front but not in the following rounds. During this whole process, compression and stretch of the plasma, a small part of ions are ejected to each direction.

Notice that even if the inital laser parameters are the same in the PIC and fluid simulation, we expect a small systematic discrepancy in the grating parameters, due to the following effect : once the plasma gratings are established, the grating's dispersion starts to have a band gap structure [19,20,36], and the driving lasers get reflected, which is not included in fluid simulation. This effect is discussed in detail in Appendix(A). Note that it is applicable here to neglect the energy exchange between the lasers, because that the growth length of the strong coupling c/γ sc (about 24.4 µm and 18 µm in Fig (2) and Fig( 3) respectively) is much larger than the plasma length. The consistency between the fluid simulation and PIC simulation approves our theoretical model and it shows that indeed the electrons are adiabatically com-pressed in the grating, µ and ν are a pair of important parameters to determine the grating behavior. The analytical solution of the fluid model in a simplified case is discussed in Appendix(B).

III. THE ROLE OF THE PARAMETER µ AS FAR AS THE ION DYNAMICS IS CONCERNED

Apparently the second stage, i.e. the stage when the velocity plateau forms, is crucial to determine the subsequent behavior of the grating. In Fig( 4), velocity plateaus with different µ are shown. At the time of this transient stage(t = t plateau ), one can compare the maximum kinetic energy of the ions P m with the amplitude of the potential barrier e∆φ = eφ maxeφ(P m )(φ max is the maximum total potential and φ(P m ) is obtained at the x where ions have the peak velocity) to know how the ions act: if P m < e∆φ, the ions will be reflected and the grating expand, as the case shown above; if P m > e∆φ, then the ions cross the potential barrier. As shown in the simulation, with a big µ(e.g. µ = 1.5), P m -e∆φ is negative. And obviously in the cold limit(µ = 0), P m -e∆φ is positive. Fig( 4) can be explained as follows: as µ decreases, the thermal potential φ th = µ(n e ) 2 in the grating needs higher grating density to stop the ions moving towards the grating. Higher grating density is obtained closer to the center as the wave steepens. Thus the grating has smaller width when the velocity plateau forms(note that the grating width and amplitude are defined as the width and electron density of the density plateau, respectively). But the increase of the grating density doesn't compensate the decrease of µ, so the potential barrier e∆φ decreases. One the other side, the ions can get closer to the center and are accelerated to higher kinetic energy. One can expect that at a certain µ, P m -e∆φ flips its sign from negative to positive and the ions starts to pass through the potential barrier. Finally, at a very small µ, the bulk of the ions can pass through the potential barrier. The grating is in the crossing regime. In both the reflection case and the crossing case, the grating starts to expand after the second stage. Thus the minimum grating width d, maximum grating amplitude (expressed by n m e ) is obtained at this time(t = t plateau ). How these parameters and P m -e∆φ evolve with µ is shown in Fig( 5) with both fluid simulation results and PIC simulation results. Fluid simulations agree with the PIC simulations quite well.

In the fluid simulations, at µ = 0.1(ν = 95.5) and µ = 0.085(ν = 19.5), P m -e∆φ flips its sign. In PIC simulation, complete reflection is found at µ > 0.25. But for 0.25 > µ > 0.001, a transition regime where the ions are firstly reflected then crossing is found, which is discussed in details in the next section. For µ < 0.001, we see directly crossing. Fluid simulations fail to describe the reflection and crossing of the ions in the transition regime and crossing regime because the fluid simulation doesn't allow the variables to be double-valued. The large deviation of negative P m -e∆φ between the fluid simulations and the PIC simulations can be observed in the last figure of Fig [START_REF] Lureau | 10 petawatt lasers for extreme light applications[END_REF]. It can be explained by the fact that slowly-growing front is applied to the lasers in PIC simulations, which leads to a smaller maximum ion velocity thus a smaller maximum ion kinetic energy P m in PIC simulations.

In the case of ν 1, it is found that ν has little effect on the main results obtained above. This can be readily known from Eq(5a) and from the simulation results with different ν in Fig [START_REF] Lureau | 10 petawatt lasers for extreme light applications[END_REF]. But note the value of ν af- fects the amplitude of the ion spikes at the edges of the grating in the complete reflecting regime. With smaller ν(still much larger than 1), the spike amplitude is larger and vice versa. But when ν is close to 1, the electron dynamics become important and can no longer be neglected, which has been proven by simulations using cold nonlinear two-fluid equations and PIC simulation (not shown). The case of ν < 1, as mentioned above, corresponds to the superradiant regime [32]. The electrons form a grating in a very short timescale, during which the ions barely move. But this paper is devoted to study the plasma grating with long time scale, which is suitable to manipulate femtosecond laser pulses. The case ν < 1 is beyond this paper's scope and left out of discussion.

IV. CHARACTERISTIC TIME SCALES OF THE NONLINEAR GRATING DYNAMICS

PIC simulations are conducted to show long-term plasma density evolution in time for different µ. For this section the parameters of the PIC simulations, except the temperature, are the same as in Fig( 2). For relatively large µ(µ ≥ 0.25), the grating forms and expands, the ions concentrate and get reflected periodically, as shown in Fig [START_REF] Uteza | [END_REF], in which µ = 1. The periodicity time T grating is roughly inversely proportional to µ, as is shown in Fig (11). As seen in Fig( 6) T grating = t 2t 1 (t 1 is the establish time of the grating, note that t 1 = t plateau . t 2 is the time that the grating stretches to its largest extent, respectively) is very regular and sustained over many periods. However, this changes if µ is decreased as the periodicity becomes irregular as discussed below. If µ is dragged back by the ponderomotive potential, and cross the center at t 4 again. After two more rounds of oscillation of particles in the ponderomotive potential, the particles are phase-mixed and no more gratings occur. The expansion of the grating due to the reflection, the subsequent shrinking and crossing at t =900 fs, 1500 fs, 1750 fs are shown in Fig( 8). The ion density evolution in time at µ = 0.02 and µ = 0 are shown in Fig( 9) and Fig (10). For smaller µ, except that at the time the velocity plateau forms, t plateau , the grating width is smaller and the grating amplitude is larger, the life time of the grating, T lif e = t 3t 1 , is shorter, the time of crossing and being dragged back by the ponderomotive force(T bounce = t 4t 3 ) is longer, as shown in Fig (11). At µ = 0, the grating only exists for a very short time and then dies out quickly.

V. PARAMETER WINDOW FOR EXPERIMENTS

The analysis of the fluid equations and the PIC simulations allows to identify two values of the µ parameter that correspond to a transition between different regimes. The threshold value below which the system goes from complete reflection to partial reflection (transition regime) is µ = 0.25. As µ decreases and reaches the threshold value µ = 0.001, the system enter the crossing regime. Another threshold µ = 2 is identified in PIC simulation (not shown), above which the grating is in the low-intensity regime, which means that the laser intensity is not high enough to produce a functional plasma grating. In this regime, as µ increases, the system slowly transitions to the weakly-coupled SBS. These three lines and the relation between µ and laser intensity I and T e are shown in Fig (12).

In the same figure we superimpose some dots that corresponds to parameters taken from the simulations in Refs [19] and [36] for comparison. For a given value of µ, and roughly independent on ν [31] the characteristics of the plasma gratings are fixed, and in particular the amplitude of the plasma grating is inversely proportional to µ. Thus for a given laser amplitude a 0 = 0.05 increasing the temperature corresponds to increasing µ and thus to obtain smaller plasma grating. Analogously, for a given initial electron temperature T 0 fixed at 25 eV, increasing the laser amplitude corresponds to decreasing µ and obtaining plasma gratings at higher amplitude. This is exactly what it is obtained in the aforementioned papers. The orange circle and the cross in the Fig (12) have very close values of µ(equal to 0.163 and 0.176, respectively),and the corresponding gratings in PIC simulation have very similar dynamics (as shown in Fig( 4) of Ref [19]). This also happens to the upward triangle and the plus sign (µ = 0.09 and 0.12), the rightward triangle and the leftward triangle(µ = 0.024 and 0.04). [36] and [19], respectively.

The result from Ref [36] (red dot) cannot be compared directly with the simulations in Ref [19] even if the value of µ is the same since in those simulations the pump lasers have finite durations and the driving ponderomotive potential exists for only 250 fs. As a consequence, the grating are not driven to grow to the maximum amplitude they can attain. In this paper we solve the model equations for a long-standing ponderomotive potential, and the results are summarized in Fig [START_REF] Lureau | 10 petawatt lasers for extreme light applications[END_REF] . We can thus predict the peak grating amplitude and formation time in Ref [19], while the grating amplitude in Ref [36] is overestimated. However the fluid model is still valid, so to obtain the value of the density grating in a general situation, it is enough to solve the model equation with a ponderomotive driver of finite duration, or as an initial condition problem. PIC simulations are necessary to obtain the subsequent evolution of the grating, however as we have shown we can identify typical behaviours associated to a given value of µ. This allows in a straightforward way to identify the plasma and laser conditions necessary in order to produce experimentally a specific grating. In the following, we examine how other parameters can also affect the the grating formation and should be considered in some particular cases.

VI. THE ROLE OF KINETIC EFFECTS FOR THE GRATING GENERATION

A. Influence of a finite ion temperature

For very low ion temperature, that is when the velocity spread in the ion phase space δv = T i0 /m i is much smaller than the maximum ions velocity in the ponderomotive potential v m : δv v m , the fluid model proposed so far gives the correct description of the system. We can deduce the condition of validity of the fluid model by considering that the maximum kinetic energy of the ions P m = m i v 2 m /2 is smaller but on the same order of the maximum ponderomotive potential(eφ m p = m e v 2 0 /2). Then the inequality above becomes T i0 /m e c 2 a 0 , usually satisfied in experiments. This is verified by PIC simulations with T i = 50 eV. The other parameters are the same as in Fig( 2), and in particular a 0 = 0.06, so that eφ m p = 0.92 keV. The results are shown in Fig( 13). The plasma grating has the same behaviour as before, only in this simulations the ion spikes at the edges of the grating are smoothed out. To include finite ion temperature a kinetic treatment that allows to describe the most energetic ions has to be considered :

∂f ∂t + v ∂f ∂x -e ∂φ ∂x ∂f ∂v = 0, ( 6a
) eφ = eφ p + eφ th = cos x + µn 2 e , (6b)
n i = f dv, (6c)
∂ 2 φ ∂x 2 = ν(n e -n i ). ( 6d
)
Finite ion temperature will give rise to two different effects. First of all, the maximum attainable amplitude of the grating will be reduced. Indeed the study of a nonlinear ion acoustic wave with large-amplitude sinusoidal velocity perturbation and no driver [37] with the waterbag model results in a limiting value of the density inversely proportional to the initial ion temperature n max 6 1/3 (v m /δv) 2/3 n 0 . An analogous effect will be present in our case. By using the estimate above we observe that the normalized thermal potential is reduced :

eφ th = µn 2 max ∝ T e0 /T 2/3
i0 . The different regimes identified in the cold ion case will still exist: reflection, transition and crossing, but the thresholds between these three regimes will imply different values of the µ parameter.

When the ion temperature is high enough a second effet is expected. Namely a steady state solution with a population of trapped ions. To understand this effect let us neglect the thermal potential and the deviation between ion and electron density, then Eq(6) becomes simply :

∂f ∂t + v ∂f ∂x + sin x ∂f ∂v = 0, ( 7
)
which describes the ions' dynamic driving by sinusoidal ponderomotive potential and is the analogous of the superradiant regime for electron waves with finite initial electron temperature [38]. This equation was studied for an initial "water-bag" velocity distribution and a static ponderomotive potential in the Hamiltonian mean-field model [39] and it was shown that the initial distribution evolves into a steadystate structure with trapped particles, characterised by the energy ε h :

f (ε) = η 0 [Θ(ε F -ε)+χΘ(ε h -ε)Θ(ε-ε F ).
Here ε is the particle energy, ε F the initial maximum particle energy in the water-bag model, and η 0 and χ parameters that depend on the fraction of trapped particles. We can expect a similar evolution in our system for a large ion temperature and an initially gaussian distribution function, leading to static plasma grating. This was also found in kinetic simulations Ref [38] for electron plasma gratings.

The trapping of the ions and the equilibrium with a steady-state ion grating implies that the average ion kinetic energy(T i0 ) should be comparable to the maximum of the ponderomotive potential(m e v 2 0 /2) as seen by the ions via the electrons. In general the contribution of the electron thermal potential, even if reduced because of the ions temperature (eφ th ∝ T e0 /T 2/3 i0 ) will not be negligible. In laboratory experiments, especially on short time scales, the electron temperature is usually higher than the ion temperature T e0 ≥ T i0 , thus both effects will be present for different populations : ion reflection and ion trapping.

PIC simulations are conducted to show the two expected effects. The simulation parameters here are the same as those in Fig( 2), except for the electron temperature and ion temperature. In Fig (14), T e0 = T i0 = 920 eV, which corresponds to µ = 1.5 and T i0 = m e v 2 0 /2. In this case we expect mainly ion trapping a small part of the ions is reflected. Indeed compared with Fig( 2), the grating amplitude is smaller and a quasi-static plasma grating appears. The time evolution of the phase space (not shown here) confirms that there is a population of trapped ions in the potential well that conserves the structure shown in the top panel on the right in Fig( 14). By keeping the same value of µ and decreasing only the ion temperature, T i0 = 200 eV , so that T i0 < m e v 2 0 /2, we do not obtain anymore a steady state solution. As shown in Fig( 15) the grating amplitude is higher than in Fig( 14), but less than that of the case with negligible T i0 , and it oscillates in time. The population of ions trapped in the potential well is not present anymore.

If we maintain this value of the ion temperature but the electron temperature is decreased as well, i.e. for smaller values of µ, we expect only a reduction of the grating amplitude if compared with a case with the same value of µ and completely negligible ion temperature, but no steady state solution. Simulations (not shown here) with T e0 = T i0 = 200 eV show indeed this behaviour: an oscillating grating is formed with maximum values of the amplitude larger than in Fig( 15) but reduced with respect to the case of smaller T i0 . The electron-ion and electron-electron collision fre-quency increase with increasing density and decreasing electron temperature to the power 3/2 times a function weakly dependent on density and temperature, the socalled Coulomb logarithm ln Λ [START_REF] Kruer | The Physics of Laser Plasma Interactions[END_REF][START_REF] Huba | [END_REF] .

B. Consideration of particle collisions

During the adiabatic compression the temperature increases and the collision frequency will thus decrease as the compression takes place and the grating forms. For example, for Fig (2), at the beginning when the plasma is not perturbed, ln Λ ≈ 7, ν ei ≈ 2.48 × 10 11 /s. When the plasma grating reaches its maximum amplitude, i.e. at time t = 500 fs (Fig( 2a)), the density in the grating increases to n e = 0.55n c , but the temperature also increases to T e 2 keV, so we have ln Λ ≈ 7.7, ν ei ≈ 1.49 × 10 11 /s. Both ν -1 ei and ν ei -1 are on the order of several ps, much longer than the typical dynamic time of the plasma grating, so we expect they will have little effect on the plasma grating formation and evolution.

The PIC simulation results are shown in Fig( 16), the electron-ion collisions are included and there are no significant differences with Fig( 2) except that the ion spikes at the edges of the grating are smoothed. The situation is different when the initial electron temperature is very low, e.g. for Fig (9), T e = 12.3 eV, µ = 0.02. In that case we have ln Λ ≈ 3, ν ei ≈ 6.88 × 10 13 /s, ν -1 ei ≈ 14.5 fs. In this case, the electrons are heated via inverse bremsstrahlung by the electromagnetic wave through the electron-ion collisions [START_REF] Kruer | The Physics of Laser Plasma Interactions[END_REF], since ν ei t unit 1. The overall temperature increased is higher than that of the case only considering the adiabatic heating, it results in a slightly reduced peak value of the grating and longer lifetime. This is shown in Fig (17): when the collisions are not considered in Fig( 17a), the grating dynamic is in the transition regime and the peak value of the grating is twice n c . When collisions are included Fig( 17b), the electron spread in the phase space is wider than in the previous case, and the electrons are heated by electron-ion collisions. As a consequence the grating maximum amplitude is slightly smaller and the time evolution of the grating dynamic increases, even if it is still in the transition regime. The typical ion crossing time becomes about 1750 fs instead of 800 fs without collisions(not shown).

The e-e collision frequency is of the same order of the e-i collisions, ν ee = ν ei /Z and Z = 1 here. These collisions will induce energy transfer among the electrons during the adiabatic compression, reducing the adiabatic temperature increase. This effect is shown by artificially switching-off of the e-i collisions. The result is shown in Fig( 17c).

To summarize, the electron-ion collisions leads to additional heating of the electrons besides the adiabatic heating, but whether the collisions are important or not depends on the initial state. If in the unperturbed plasma the reciprocal of the collision frequency is already on the same order of or even higher than the dynamic time of the grating it will have a strong effect on the grating dynamic. In the opposite case in which the collision time is longer than the dynamic time of the grating, the effect of collisions can be neglected.

Moreover we wish to point out that the model we propose is not limited to lasers in the visible domain interacting with under-dense plasmas. It can be applied to X-ray laser-generated plasma grating, laser-generated plasma grating on solid surface, and infrared-laser-generated plasma grating. In the two former cases, the collisions are often important, in the last one they are negligible.

VII. CONCLUSIONS

The laser-generated plasma grating induced by strongly coupled stimulated Brillouin scattering is characterized by large-amplitude and long-lasting ion density modulation [31], which can last much longer than the electron grating [38,[START_REF] Chen | [END_REF]. Combined with its high damage threshold, this kind of plasma grating is useful to manipulate ultraintense and ultrashort lasers. Many important applications are proposed theoretically. To experimentally realise these applications, it is necessary to understand the mechanism of the formation, collapse process and how experimental parameters affect the properties of the plasma grating.

In this paper, a fluid model and adiabatic compression assumption is proposed to explain the nonlinear dynamics of plasma grating. Good agreements are found with the corresponding PIC simulation in specific cases, which validate the fluid model and adiabatic compression assumption. Proofs of the adiabatic heating are also found in PIC simulation. Further simulation with different µ in both fluid and PIC shows that as µ decreases, maximum grating amplitude is higher and the minimum grating width is smaller. Three different regimes of dynamics of the plasma grating are also found: complete reflection, transition regime and crossing. The complex dynamics of the grating in the transition regime are discussed in details. Analytical solutions of the simplified model are also presented, which help explain the grating's behavior. The parameter window on how µ depends on laser intensity and electron temperature is also shown, which facilitates choosing suitable experiment parameters to obtain plasma grating with desired properties. The finite ion temperature is found to have two effects on the grating dynamics. First, as the ion temperature increases, the attainable grating amplitude is smaller. Second, when the ion kinetic energy is comparable to the ponderomotive potential, the ions will be trapped steadily in the potential. One may think that the localization of particles in the grating can lead to severe collisions, which can not be ignored in our model. But we show that as the grating forms, the adiabatic compression heats the electron besides localizing the ions, which mitigates the collisions. Whether the collisions need to be considered or not depends on the initial state of the plasma.

This model also helps to explain the formation of plasma density modulation in some special cases. For example, this may be able to explain the periodic density grating driven by the beating of laser projected on the surface of a solid target and the reflected laser [28]. The periodic formation and collapse of grating has been found in the case of high-intensity laser interacting with plasma [34,43], which can explain the periodic changes of the reflectivity of SBS in inertial confinement fusion.

This work is focused on the mechanism of the grating dynamics and one-dimensional situation. In the case of small-scale plasma and long laser pulses, in which the effects of energy transfer between lasers and the laser intensity distribution can be neglected, our model provides important instructions for practically building a plasma grating device. But in real experiments, these approximations usually do not hold, then it is necessary to further investigate how the energy transfer between lasers and the laser distribution in 3-dimensional space affect the grating formation.

VIII. ACKNOWLEDGEMENTS

authors wish to acknowledge useful discussions with F. Amiranoff and J. Wurtele. H.P. acknowledges the discussions with R. Li. This work has been done within the LABEX Plas@par project, and received financial state aid managed by the Agence Nationale de la Recherche, as part of the program "Investissements d'avenir" under the reference ANR- Data from PIC simulation of Fig( 2) is analysed. The results at t = 500 fs are shown in Fig (18), in the first row one finds that the plasma grating has been fully established at this time.

The second row shows the ponderomotive potential characterized by the mean squared transverse electrical field |E y | 2 . Note that in the free space, is smaller in the center which should be explained by that the grating reflects the following injected lasers, the lasers' transmittance drops as the lasers get closer to the center. If with enough layers of plasma grating, one could see the ponderomotive potential drops to zero in the center.

The third row shows the lasers propagating to each direction: E right and E lef t , note that

E right = (E y + v p c B z )/2, (A3a)
E lef t = (E y - v p c B z )/2, ( A3b
)
v p = c 1 -n e /n c , (A3c)
where v p is the phase velocity of the lasers, n e is the local electron density. E right and E lef t are attenuated by the grating because of the reflection in Fig (18). It is proven here that the ponderomotive potential actually decreases in the grating when the grating is fully established because the grating reflects the pump lasers. The ponderomotive restores to it's original level after the grating collapses. This effect is not included in Sec(??), it partly attributes to the differences between the fluid model and the PIC simulation. But the PIC simulations already supports the model well, and give enough information of the grating dynamics. While for the further study of plasma grating dynamics in a large-scale plasma, both the energy exchange between pump lasers and the reflection by the grating should be considered.

(-1) n 2J n (nt), and in the center of the grating center(at x = π) n i = 1/(1t), ∂v i /∂x = -t/(1t), noting that cos(nπ) = (-1) n and ∞ n=1 J n (nt) = t/2(1t). This means that as the ions keep accumulating to the grating center(x = π) and time goes from t = 0 to t = 1, the grating density goes to infinity and the slope of the ion velocity goes to minus infinity. At t = 1 the grating finally breaks due to that the ion sheets cross each other. When µ > 0, the results given above are quite complex to analyses. But one can cut the infinite series to the first several terms which gives very good approximations, and the numerical results can be provided readily by scientific softwares, e.g. Mathmatica. Fig (19) shows the ion velocity and density at different times with µ = 1.0 as given by Eq(B6) cut at n = 10. One can see clearly the formation of velocity plateau and density grating at t = 0.7, the ion reflection and grating stretch at t = 1.5. Fig (20) shows the temporal evolution of the maximum grating amplitude along x with different µ by cutting Eq(B6c) at n = 10. At finite µ, the ion density doesn't go to infinity, and as µ gets bigger then peak grating amplitude decreases, which agrees well with the simulations in SectionIII.For the limit that µ is very big, Eq(B1) can be linearized with n i = 1 + ñ, v i = ṽ and has solutions:

n i = 1 -sin(x - √ 2µt)/ √ 2µ, v i = sin(x - √ 2µt
). The velocity and density perturbation are quite small and linear.

In this appendix, the analytical solutions of homogeneous fluid model(no driving ponderomotive potential but the ion velocity has an initial sinusoidal perturbation) are which supports the simulation results in previous sections. They can also be used to predict the grating behavior conveniently.

FIG. 1 :

 1 FIG. 1: Simulation results of the fluid model with µ = 1.5 and ν = 59.5. Shown are ion velocity, plasma density, thermal potential, ponderomotive potential and total potential at different times, respectively. All the potentials are normalized to the maximum ponderomotive potential: φ m p = me/2ev 2 0 .

FIG. 2 :

 2 FIG. 2: The results of the first PIC simulation at different times. Note that Fig(a),(b),(c) are at the same normalized times as Fig(a),(b),(c) of Fig(1), respectively. tunit 320 fs here. From upper row to lower row are the ion phase space, electron phase space and local electron temperature(green line, ticked on the right), ion density(red line) and electron density(blue line) at different times in the center of the simulation box, respectively. The ion velocity and plasma density in normalized units are ticked in the right. Note that a Savitzky-Golay filter is applied to smooth the electron temperature to increase the precision without distorting the trend.

FIG. 3 :

 3 FIG. 3: The results of the second PIC simulation at different times. Note that Fig(a),(b),(c) are at the same normalized times as Fig(a),(b),(c) of Fig(1), respectively. tunit 256 fs here. The layout is the same as that of Fig(2).

FIG. 4 :

 4 FIG.4: Fluid simulation of the second stages with different µ. ν is fixed at 59.5. Normalized units are used.

FIG. 5 :

 5 FIG. 5: From top to bottom are: the minimum grating width d, the time to form the velocity plateau t plateau , the maximum electron density n m e , Pm -e∆φ(normalized by the maximum ponderomotive potential eφ m p = mev 2 0 /2) at t plateau , respectively. Fluid simulation results are shown with blue lines(ν = 59.5) and black broken line(ν = 19.5). The PIC simulation with varying temperature are shown with green triangles(I = 5 × 10 15 W/cm 2 , ν = 59.5) and orange circles (1.5 × 10 16 W/cm 2 , ν = 19.5), both with n0 = 0.3nc and other parameters the same as those of Fig(2).

8 FIG. 6 :

 86 FIG. 6: Ion density evolution in time at the center of the simulation box. The periodic formation and stretching of the plasma grating at µ = 1.0 is seen. The two insets on the right show the zooming in of the main figure at different times.Note that t1 = 540 fs, t2 = 1415 fs.

6 FIG. 7 :

 67 FIG. 7: Ion density evolution in time for µ = 0.1. The three insets on the right show the zooming in of the main figure at different times. Note that t1 = 500 fs, t2 = 900 fs, t3 = 1500 fs, t4 = 2000 fs.

3 FIG. 8 :FIG. 9 :

 389 FIG.8: Shown are the ion phase space, electron phase space and plasma density at different times when µ = 0.1.

5 FIG. 10 :

 510 FIG. 10: Ion density evolution in time at µ = 0.0. The three insets on the right show the zooming in of the main figure at different times. Note that t1 = 560 fs, t3 = 635 fs, t4 = 1860 fs.

FIG. 11 :

 11 FIG. 11: Left: the time period of grating formation and collapse Tgrating for relatively big µ. Right: the life time of the grating T lif e and the following bouncing time of ions in the ponderomotive potential T bounce for relatively small µ.

FIG. 12 :

 12 FIG.12:The relation between µ and laser intensity and Te. The upper broken line: µ = 2. The middle broken line: µ = 0.25. The lower broken line: µ = 0.001. The orange dots and red dot are from G. Lehmann's simulation[36] and[19], respectively.

FIG. 13 :

 13 FIG. 13: PIC simulation with Ti = 50 eV, others are the same as Fig(2). t = 500 fs, 800 fs, 1300 fs in (a), (b), (c), respectively.

FIG. 14 :

 14 FIG. 14: Left: the amplitude of the grating(denoted by ni/n0) versus time. Right: from up to bottom: ion phase space, electron phase space, ion density(orange line) and electron density(blue line) at t = 1500 fs. In this simulation, Te0 = ZTi0 = 920 eV, Z = 1.

FIG. 15 :

 15 FIG.15:The same as Fig(14), except that Ti0 = 200 eV.

FIG. 16 :

 16 FIG. 16: PIC simulation with the same parameters as Fig(2), except that electron-ion collisions are included. Times shown are (a)t = 500 fs, (b)t = 800 fs, (c)t = 1300 fs.

3 FIG. 17 :

 317 FIG.17: PIC simulation at t = 500 fs with the parameters same as Fig(9), i.e. µ = 0.02, except (a) is without collisions, (b) is with all the collisions(electron-ion,electron-electron,ionion), (c) is only with electron-electron collision.

|E y | 2 = 2 =

 22 |E y1 + E y2 | |E 0 cos(kxωt) + E 0 cos(-kxωt + ∆φ)| 2 = E 2 0 cos 2 (kxωt) + E 2 0 cos 2 (-kxωt + ∆φ) + E 2 0 [cos(-2ωt + ∆φ) + cos(2kx -∆φ)],(A1)where E y1 and E y2 are the transverse electrical fields of two pump lasers, E 0 is their amplitude, ω and k are the pump frequency and wavevector, ∆φ is their phase discrepancy, respectively. After averaged in one optical cycle, i.e. T = 2π/ω:|E y | 2 = E 2 0 [1 + cos(2kx -∆φ)]. (A2) It is well known the spatial derivative along x of |E y | 2 is the ponderomotive force generated by the beating of two pump lasers. |E y | 2 is shown in normalized unit. E unit = m e cω/e, for pump lasers with a 0 = 0.06, E 0 = 0.06E unit . So for laser without reflection, |E y | 2 should range from 0 to 7.2e -3. While in the case the grating is fully established as shown in Fig(18), |E y | 2

5 FIG. 20 :

 520 FIG.20:In the grating center, the ion density evolution in time with different µ according to the analytical solutions.

  The results of the Project LQ1606 were obtained with the financial support of the Ministry of Education, Youth and Sports as part of targeted support from the National Programme of Sustainability II (SW).

	pean Regional Development Fund. APPENDIX A: REFLECTION OF PUMP LASERS
		BY ESTABLISHED PLASMA GRATING
	c n/n	0.00 0.25 0.50 0.75 1.00		electron ion				t= 500 fs
		[E 2 unit ]	4 8	×10 -3	|Ey| 2
				0
	[E unit ]	0 -0.10 -0.05 0.00 0.10 0.05	1	2	3	4	5 x/λ 0	6	7	8	9	10
	11-IDEX-0004-
	02. H.P acknowledges the funding from China Schol-
	arship Council, and is partially founded by the Nat-
	ural Science Foundation of China and 11875240 and
	Key lab foundation of CAEP 6142A0403010417. S.W.
	was supported by the project Advanced research us-
	ing high intensity laser produced photons and par-
	ticles (ADONIS) (CZ.02.1.01/0.0/0.0/16 019/0000789)
	and by the project High Field Initiative (HiFI)
	(CZ.02.1.01/0.0/0.0/15 003/0000449), both from Euro-

Eright Eleft

FIG.

18

: From up to bottom are the plasma densities, mean squared transverse electrical field, transverse electrical field transporting to the right and the left at t = 500 fs, respectively.

In Section(III), we have shown that ν has little effect on the main results as long as ν 1. So in the model function Eq [START_REF] Lureau | 10 petawatt lasers for extreme light applications[END_REF], Eq(5a) can be neglected and n e can be replaced with n i in Eq(5c). The model is still strongly nonlinear and with a spatially inhomogeneous term sin x, one can simplify the model and the solution by assuming that the velocity has an initial sinusoidal distribution and no spatially driving term. Then the model is simplified to:

with the initial conditions:

Then multiply Eq(B1a) with ± √ 2µ and add to Eq(B1b), we obtain:

with u = v i ± √ 2µn i and the initial condition u| t=0 = sin x ± √ 2µ. Transforming from Eulerian variables (x, t)

to Lagrangian variables (x 0 , τ ):

Eq(B2a) is just:

or x ∓ √ 2µτ = x 0 + τ sin x 0 . We then know that sin x 0 is a periodic function of x ∓ √ 2µτ with period 2π, and it can be transformed back to Eulerian variables with the following Fourier-series [44]:

where J n is the nth Bessel function of the first kind. Then u has two solutions:

) and:

c n (t) = (-1) n 2 sin √ 2µnt √ 2µnt

J n (nt). (B6d)

When µ = 0, the solutions converge to the cold limit, we have: b n (t) = (-1) n+1 (2/nt)J n (nt), c n (t) =