Patrick Dupré
email: pdupre@gmx.com

Dipole Saturated Absorption Modeling in Frequency Modulation Spectroscopy: Dealing with a Gaussian Beam, Resonance Narrowing

Keywords: Saturated Absorption, Lamb-dip, Gaussian beam, transit-time broadening, acetylene, collision rate, FMS, modeling

Quantitative, nonlinear sub-Doppler sensitive spectroscopy requires reliable formalism for optimizing the detection of weak molecular absorption patterns. Phase modulation spectroscopy (FMS) is a common technique usedto challenge the ultimate photon-shot-noise limit. In addition, nonlinear spectroscopy can exhibit surprising effects like resonance narrowing, i.e., resonance width below the transit-time rate, like it has been pointed out in the seventies. Here, by extending our previous work, we propose an analytical development associated with numerical integration to simulate the equivalent Lamb-dip under Gaussian beam conditions. Simulations of nonlinear profiles of two transitions belonging to the polyad ∆P 11 of acetylene are discussed, and analyzed by fitting with standard functions. The respective roles of the transit-time, of the collision rate, and of the Rabi frequency (including the power broadening) are carefully discussed and reviewed. It shows "super-narrowing" effects under specific low pressure and low power conditions. A reviewing introduction states previous experimental reports as well about previous modeling.

A C C E P T E D M A N U S C R I P T

• Formalism for calculating the Saturated Absorption under Frequency Modulation Spectroscopy • Gaussian Beam characterized by a Maxwellian transit-time rate distribution are targeted • Profile simulations are discussed versus the mean Rabi frequency, the collision rate and the transittime rates

• Broadening/narrowing, and shape of the FM dispersion profiles are discussed

• The study focus on 2 R(0) transitions of C 2 H 2 (polyad 11).

1

Introduction

Since the 1970's, saturation spectroscopy of gaseous species has received considerable attention with the advent of the laser sources. Nevertheless, it still requires for modeling refinements, because of the onset of new, very accurate and sensitive spectrometers. For example, the possibility to address gas under very low pressure conditions relevant in the frequency-standards context, requires taking into careful considerations the spatial shape of the interacting photon beam. Intuitively, the weak collision rate stresses the pivotal role of the "interaction duration" between the molecular (or atomic) species and the probing photons resulting from the finite spatial extent of the electromagnetic field. Highly sensitive absorption techniques allows probing the line shape beyond the Voigt profile, and even beyond the most elaborated line profiles under some specific conditions. In the infrared (IR) domain where the spontaneous emission rates are negligible, only techniques based on absorption, (and on stimulated emission) processes are experimentally conceivable. Potentially, this also allows challenging the elastic collisional processes associated with the long range interaction. Indeed, the molecular longitudinal speed selection associated with the absorption/emission of two counter propagating beams requires us to consider the velocity-changing collisions (VCC) and phase changing collision (PCC) in a context of the field saturation, i.e., when the Rabi splitting can compete with processes altering the inelastic collisional environment as suggested by some authors (vide infra).

Absorption techniques, based on a single energy step, are not intrinsically extremely sensitive, because the molecular transitions are detected over a bright background. Hence, the ultimate detection limits are supplied by the usual background photon-shot-noise. Nevertheless, targeting ultimate levels of sensitivity requires elaborated experimental design. In addition, ultra-high resolution spectroscopy techniques are based on Doppler-free molecular responses, which requires population transfers involving two-photon processes (like saturated absorption or Doppler-free twophoton absorption). A way to improve the sensitivity (i.e., the signal-to-noise ratio [SNR]) is, of course, to implement optical cavities, because these devices provide large enhancement factors (mostly the cavity finesse). On the other hand, the residual noise needs to be properly addressed. Employing radio-frequency (RF) heterodyne techniques allows eliminating the intensity fluctuations of the source (noise in 1/ f ). Phase (or Frequency) Modulation Spectroscopy (FMS) techniques [START_REF] Bjorklund | Frequency-modulation spectroscopy: a new method for measuring weak absorptions and dispersions[END_REF] have been implemented, either in the free space by partially or entirely reflecting the beam over itself, or by using optical cavity as original demonstrated by J. Hall group [START_REF] Ye | Ultrastable optical frequency reference at 1.064 µm using a C 2 HD molecular overtone transition[END_REF] (inside optical cavities, the saturation absorption conditions are de facto satisfied, because of the beam counter propagation). Heterodyne techniques of detection can also be employed to challenge dispersion profiles.

We recently reported on the Lamb-dip shape associated with the saturated absorption under different collisional regimes [START_REF] Dupré | Dipole saturated absorption modeling in gas phase: Dealing with a Gaussian beam[END_REF], pointing out the relative roles of the Rabi (flopping)frequency, and of the transit-time broadening (or timeof-flight broadening) associated with the transverse extension of Gaussian beams (as opposed to usual plane waves).

This applies to usual setup without modulation, including Cavity Ring Down Spectroscopy (CRDS). The following step is obviously to address these saturation conditions when a phase modulation is applied. This is the purpose of the current paper addressing the SA-FMS (Saturated Absorption FMS), extending the previous studies, i) by proposing a formalism based on a Gaussian electro-magnetic field (EMF), ii) by simulating some specific transitions of acetylene, iii) by extracting generic behaviors based on fitting the simulated profiles by standard functions.

A C C E P T E D M A N U S C R I P T

The role of the transit-time, has been discussed multiple times since pioneer experiments have been performed under laser saturation. Unexpectedly narrow resonances have been identified, and authors introduced the concept of super-or ultra-narrow resonances [START_REF] Barger | Pressure shift and broadening of methane line at 3.39 µm studied by laser-saturated molecular absorption[END_REF][START_REF] Bagaev | Anomalous decrease of the shift of the center of the Lamb dip in low-pressure molecular gases[END_REF][START_REF] Thomas | Transit-time effects in power-broadened Doppler-free saturation resonances[END_REF][START_REF] Chebotayev | Supernarrow saturated absorption resonances[END_REF], hypothesizing on the role of the slow molecules [START_REF] Thomas | Transit-time effects in power-broadened Doppler-free saturation resonances[END_REF][START_REF] Ch | Ultrahigh-resolution saturation spectroscopy using slow molecules in an external cell[END_REF][START_REF] Ma | Ultrasensitive frequency-modulation spectroscopy enhanced by a high-finesse optical cavity: theory and application to overtone transitions of C 2 H 2 and C 2 HD[END_REF]. Such behavior has been reported by J. Hall in an extensive study of CH 4 obtained by probing both the low pressure and power dependencies [START_REF] Hall | Saturated Absorption Line Shape[END_REF]. More recently, the Gaussian beam expansion has already been considered by O. Axner group in the Doppler-broadened regime of NICE-OHMS (RF modulation in cavity) experiments [START_REF] Foltynowicz | Wavelength-modulated noise-immune cavity-enhanced optical heterodyne molecular spectroscopy signal line shapes in the Doppler limit[END_REF][START_REF] Foltynowicz | Fiber-laser-based Noise-Immune Cavity-Enhanced Optical Heterodyne Molecular Spectrometry[END_REF]. The current paper interrogates this concept under the light of transit-time rate and of the Rabi frequency.

In this study, the density matrix formalism in the frequency domain has been performed to determine the saturated absorption induced in a quasi two-level system, and in the electric dipole moment approximation, paying attention to the degenerated Zeeman structure [START_REF] Pinard | Third-order nonlinear optical response resulting from optical pumping: Effects of atomic motion[END_REF] giving rise to the sub-transitions controlled by the impinging EMF polarization (an isotropic space is assumed). We first solved analytically the system of coupled equations for a specific Doppler shift. This approach, based on a perturbative expansion provides the complete frequential response. The relevant tensorial components of the complex third order molecular susceptibility are established for an optically thin medium (only stationary responses are eventually considered here). The "multichromaticity" of the EMF is considered. The full implementation requires for from one to three numerical integrations to take into consideration the distribution of the molecule velocity, of the beam impact parameter, and eventually of the distribution of the transit-time rate associated with a normalized 2-dimension Maxwell-Boltzmann velocity distribution. Crossover resonances possibly generated in a system involving more than two levels have been previously reported for an unmodulated EMF in the plane wave approximation [START_REF] Dupré | Saturated absorption and crossover resonances in a high-finesse cavity: Formalism and application to the hyperfine structure of jet-cooled NO 2 by saturated-absorption cavity-ring-down spectroscopy[END_REF], similarly the crossover sideband resonances resulting from the multi components of a modulated EMF are reasonably well predicted in the frame of a 2-level system. However, the complexity associated with a 3-level system in modulated EMF (yielding multiple crossover resonances) is postponed to future works.

Generic behaviors will be discussed by running simulations for two specific relatively weak transitions of acetylene (in Sec. 3). The global behavior of these simulations is examined, but, we will mainly draw our attention on the shape of the (Doppler-free) central resonance imprinting the transition observed by collecting the dispersion (optical phase shift), since it is a trivial observation in any FMS experiment. It addition, it is the resonance to analyze when targeting the transition frequency determination for metrology purposes [START_REF] Cozijn | Sub-Doppler frequency metrology in HD for test of fundamental physics[END_REF], for VCC and PCC [START_REF] Berman | Theory of collision effects on line shapes using a quantum-mechanical description of the atomic center-of-mass motion-application to lasers[END_REF][START_REF] Berman | Theory of saturation spectroscopy including collisional effects[END_REF][START_REF] Cygan | Spectral line-shapes investigation with Pound-Drever-Hall-locked frequency-stabilized cavity ring-down spectroscopy[END_REF], for determining the molecular hyperfine structure resulting from quadrupole molecular Hamiltonians, and eventually for testingthe Physics beyond the standard model [START_REF] Salumbides | Bounds on fifth forces from precision measurements on molecules[END_REF]. Currently, focusing on the resonance shape, we ignore the frequency line position which also depends on the experimental environment (i.e., collisions, EMF shape, incident power, etc..).

In a model without VCC, the collision rate does not depend on the molecule speed while the transit-time rate depends on the molecule transverse speed (relatively to the direction of the photon beam). It results that the molecules may experience a collisional regime or a collision-free regime, according to their speed, depending on the value of the product T 1 γ tr relatively to 1 (T 1 is the longitudinal relaxation time, and γ tr is the mean transit-time rate, the elastic collisions are ignored). However, both quantities contribute to the degree of saturation. Then, the contribution to the saturation of the slow molecules versus the fast molecules differs. This inhomogeneous environment give rise to complex resonance shapes, like it as been reported by several authors [START_REF] Ch | Ultrahigh-resolution saturation spectroscopy using slow molecules in an external cell[END_REF][START_REF] Ma | Ultrasensitive frequency-modulation spectroscopy enhanced by a high-finesse optical cavity: theory and application to overtone transitions of C 2 H 2 and C 2 HD[END_REF] Mostly, the current formalism exhibits some minor basic changes compared with the previous one published in recent papers [START_REF] Dupré | Saturated absorption and crossover resonances in a high-finesse cavity: Formalism and application to the hyperfine structure of jet-cooled NO 2 by saturated-absorption cavity-ring-down spectroscopy[END_REF][START_REF] Dupré | Sub-Doppler noise-immune cavity-enhanced optical heterodyne molecular spectrometry modeling: from Doppler broadening to cross-sideband resonances[END_REF][START_REF] Dupré | Dipole saturated absorption modeling in gas phase: Dealing with a Gaussian beam[END_REF] based on an intense use of the Fourier domain. These changes take into account some formal inaccuracies, and allow us to extend the initial approach to the Gaussian EMF modulation. Nevertheless, the tracking of the equations requires some approximations which are clearly documented. One is relevant in the framework of the large Doppler broadening versus the transit-time rate, while the other one alters the amplitude of the crossover sideband resonances. If our model involves a single, but multi-featured electromagnetic field pattern, the concept of pump-probe interaction is not specifically considered: the main saturation is due to the central carrier frequency for both, the back and the forth going EMFs, while a pump-probe approach could be more relevant to calculate the crossover sideband resonance contributions (barely discussed here). Hence, additional refinements may be required, but they have been postponed to further developments at the light of forthcoming experimental data.

The concepts of both transit-time rate and Rabi frequency have been redefined to match our purposes.

Generally speaking, the current formalism is an improvement of models proposed by several authors who set up a saturation coefficient (or degree of saturation) and, who ignore the individual magnetic sub-transitions, while these sub-transitions have a specific "saturation factor". This remark addresses the formalism developed to analyze CRDS experiments [START_REF] Giusfredi | Theory of saturated-absorption cavity ring-down: radiocarbon dioxide detection, a case study[END_REF], or NICE-OHMS experiments [START_REF] Axner | Sub-Doppler dispersion and noise-immune cavity-enhanced optical heterodyne molecular spectroscopy revised[END_REF]. Furthermore, the crossover sideband resonances are usually not simulated in FMS. In addition, the regime of the low pressure (i.e., less than 1 Pa) which is associated with collision rates less than the transit-time rate requires for dealing with Gaussian EMF as introduced by K. Shimoda [START_REF] Shimoda | Line Broadening and Narrowing Effetcs[END_REF][START_REF] Dupré | Dipole saturated absorption modeling in gas phase: Dealing with a Gaussian beam[END_REF], and by C. Bordé and co-authors [START_REF] Bordé | Saturated absorption line shape: Calculation of the transittime broadening by a perturbation approach[END_REF]. The mean transit-time rate, and the transit-time rate distribution have been clearly

established. In the current context, the treatment has required for abundantly dealing with the Faddeeva function [START_REF] Abramowitz | Handbook of mathematical functions with Formulas, Graphs, and Mathematical Tables[END_REF],

an extension of the absorption Voigt profile to the dispersion profile. This is a key function to analyze the resonance profile in SA-FMS. Nevertheless, the limit of such a profile is discussed.

To illustrate the current formalism, simulations target two specific transitions of acetylene which should be easy to probe by setting up appropriate experiments. Acetylene has drawn a high spectroscopy interest [START_REF] Amyay | Vibration-rotation spectroscopic database on acetylene X 1 Σ + g ( 12 C 2 H 2 )[END_REF]. But, it also attracted the interest of the International Committee for Weights and Measures, because of its relatively strong transition in the InfraRed whose frequency can be determined by using the modern frequency etalons which are the Cs-referenced Optical Frequency Comb (OFC). Indeed, rovibrational lines of C 2 H 2 have been used as secondary frequency standards for the important field of optical telecommunications [START_REF] De Labachelerie | Ultranarrow 13 C 2 H 2 saturated-absorption lines at 1.5 µm[END_REF][START_REF] Nakagawa | Accurate optical frequency atlas of the 1.5-µm bands of acetylene[END_REF][START_REF] Onae | Optical frequency standard at 1.5 µm based on doppler-free acetylene absorption[END_REF][START_REF] Quinn | Practical realization of the definition of the metre, including recommended radiations of other optical frequency standards[END_REF][START_REF] Ch | Absolute frequency measurement of a 1.5-µm acetylene standard by use of a combined frequency chain and femtosecond comb[END_REF][START_REF] Van Leeuwen | Measurement of pressure-broadened, ultraweak transitions with noise-immune cavity-enhanced optical heterodyne molecular spectroscopy[END_REF][START_REF] Edwards | High-accuracy frequency atlas of 13 C 2 H 2 in the 1.5 µm region[END_REF][START_REF] Twagirayezu | Frequency-comb referenced spectroscopy of v 4and v 5 -excited hot bands in the 1.5 µm spectrum of C 2 H 2[END_REF]. Furthermore, C 2 H 2 can be used as a molecular target for the spectroscopic determination of the Boltzmann constant by means of Doppler Broadening Thermometry (DBT). [START_REF] Gianfrani | Linking the thermodynamic temperature to an optical frequency: recent advances in Doppler broadening thermometry[END_REF] It is also a benchmarking molecule for VCC/PCC [START_REF] Forthomme | Application of the Hartmann-Tran profile to precise experimental data sets of 12 C 2 H 2[END_REF][START_REF] Fasci | Dual-laser absorption spectroscopy of C 2 H 2 at 1.4 µm[END_REF][START_REF] Twagirayezu | Frequency measurements and self-broadening of sub-Doppler transitions in the v 1 + v 3 band of C 2 H 2[END_REF] while F. Constantin also identified a few transitions with strong sensitivity to the electron-to-proton mass ratio on purpose of testing Physics beyond the standard model. [START_REF] Constantin | Sensitivity to electron-to-proton mass ratio variation from 12 C 2 H 2 rovibrational transitions to v 1 + v 3 and v 1 + v 2 + v 4 + v 5 interacting levels[END_REF] To validate our models (actually four models depending mainly on the number of integrations performed), we selected two transitions of acetylene belonging to two different vibrational combination modes, but, belonging to same polyad ∆P = 11, and exhibiting a different line intensity. To minimize any possible issue related to the hyperfine structure of the levels, two para transitions R e (0) have been chosen. To give more relevance to the model, two Rabi

A C C E P T E D M A N U S C R I P T
frequencies have been selected with significantly different values relevant for the collision rate dependence.

The analysis of the resonances shape has been mainly obtained by varying the value of the collision rate, an alternative to the gas pressure, because of the lack of well-established proportionality between the collision rate against the gas pressure under saturation conditions: specifically, self-pressure broadening coefficients available in the literature can exhibit a pressure dependence (factor 2 to 3) for methane [START_REF] Barger | Pressure shift and broadening of methane line at 3.39 µm studied by laser-saturated molecular absorption[END_REF][START_REF] Hall | Measurement of methane hyperfine structure using laser saturated absorption[END_REF][START_REF] Hall | Saturated Absorption Line Shape[END_REF][START_REF] Chebotayev | Nonlinear narrow optical resonances induced by laser radiation[END_REF], for carbon dioxide [START_REF] Meyer | High-resolution line broadening and collisional studies in CO 2 using nonlinear spectroscopic techniques[END_REF][START_REF] Vasilenko | Nonlinear dependence of optical resonance widths at CO 2 transitions on pressure[END_REF],

for iodine [START_REF] Goncharov | Study of relaxation of molecular iodine by the method of saturated absorption spectroscopy[END_REF], for water [START_REF] Chen | A phase-stable dual-comb interferometer[END_REF], and for acetylene [START_REF] Ma | Ultrasensitive frequency-modulation spectroscopy enhanced by a high-finesse optical cavity: theory and application to overtone transitions of C 2 H 2 and C 2 HD[END_REF][START_REF] Trivikram | Relativistic and QED effects in the fundamental vibration of T 2[END_REF][START_REF] Hua | Cavity-enhanced saturation spectroscopy of molecules with sub-kHz accuracy[END_REF]. The dependence of the resonance shape with the mean transit-time rate is also discussed.

Simulation analyses have been performed to validate our models. The final confidence in the simulated profiles depend on the number of numerical integrations considered. However, only two models involving a transit-time rate have been mainly considered through a unique mean transit-time, or a Maxwellian distribution of the transit-time rate.

To conduct this study, the dispersion shape of the central FM resonance is mainly fitted by using the dispersive part of a Faddeeva profile. This analysis clearly shows three different regimes: collision-free, collisional and intermediate. These regimes are associated with the relative values of the collision rate to the mean transit-time rate, and to a mean Rabi frequency. This allows us to define generic behaviors. In this context, the resonance narrowing is discussed as well as the limits of validity of the models versus the collision rate and the impinging power. Actually, the collision-free regime can potentially induce high saturation degrees which are beyond the scope of the current approach. Experimental data (under consideration) would help to benchmark the conclusive behaviors. Nevertheless, we show that recovering the input simulation parameters (or similarly experimental parameters), from the fitting processes is not a trivial task, despite it is one of the original goals.

Specifically, this work is an alternative approach of the resonance narrowing like previously reported by C. Chardonnet and co-authors [START_REF] Ch | Ultrahigh-resolution saturation spectroscopy using slow molecules in an external cell[END_REF]. Theoretically, it is not limited to some particular ranges of thermal velocity, of pressure broadening, or of Rabi frequency; all these parameters eventually contribute to the "saturation factor". This approach is certainly not unique, the approximations linked to the coupled system of equations and to the generic Gaussian shape of the beam considered for performing tractable simulations, certainly appeals for additional refinements of the model.

However, the "first order" perturbative approach considered, here requires caution when discussing strong saturation, cf. see Ref. [START_REF] Autler | Stark effect in rapidly varying fields[END_REF][START_REF] Thomas | Transit-time effects in power-broadened Doppler-free saturation resonances[END_REF] for a more complete Taylor expansion. Indeed, experimental data may even emphasize the current purposes, particular in the low pressure regime despite that elastic collision rate should be carefully considered.

Power mode pulling, mode splitting, mode beating [START_REF] Bennett | Hole burning effects in a He-Ne optical maser[END_REF], relativistic Doppler shift and recoil effects [START_REF] Basov | Optical frequency standards[END_REF][START_REF] Hall | Direct optical resolution of the recoil effect using saturated absorption spectroscopy[END_REF][START_REF] Bagayev | Second-order Dopplerfree spectroscopy[END_REF] which can generate asymmetric shape [START_REF] Kol'chenko | Interaction of an atom with a strong electromagnetic field with the recoil effect taken into consideration[END_REF] are disregarded from the present formalism (actually the recoil effect [doublet] can be easily set in the equations [START_REF] Kol'chenko | Interaction of an atom with a strong electromagnetic field with the recoil effect taken into consideration[END_REF][START_REF] Rautian | Saturation effects for long-lived systems in spatially bounded fields[END_REF]). The second order effects due to the molecular absorption and dispersion [START_REF] Mugino | High sensitive detection of trace gases using optical heterodyne method with a high finesse intra-cavity resonator[END_REF] are neglected, because an optically thin medium is assumed. The present simulations also ignore the possible optical pumping and N -scheme crossover resonance since R e (0) transitions preclude observing such effects if the hyperfine couplings are ignored.

The proposed formalism is applicable to non radiative two-level systems interacting with any derivable spectral shape of a classical EMF. However, it has been particularized to Gaussian-shaped beam that is the shape which is associated with the usual continuous wave EMF, or with an EMF trapped inside an optical cavity if it is mode matched

A C C E P T E D M A N U S C R I P T
against a TEM 00 cavity eigen mode. The spatial variations along the z axis are ignored through a mean value (cylindrical beam). However, the beam divergence can be readily added to the model by running an additional integral along the z axis.

The fields of application of the current formalism are typically the determination of the pressure broadening coefficients. However, this may rise new definition of these coefficients beyond the standard collisional models.

Of course, FMS without cavity (but with two equivalent counter propagating beam intensities) [START_REF] Hall | Optical heterodyne saturation spectroscopy[END_REF], the wavelength modulation spectroscopy (WMS) in cavity [START_REF] Moses | High-sensitivity laser wavelength-modulation spectroscopy[END_REF], and the NICE-OHMS [START_REF] Ye | Ultrastable optical frequency reference at 1.064 µm using a C 2 HD molecular overtone transition[END_REF] are specifically targeted.

Species trace detection by saturated absorption is possible as demonstrated by Ma et al. [START_REF] Ma | Ultrasensitive frequency-modulation spectroscopy enhanced by a high-finesse optical cavity: theory and application to overtone transitions of C 2 H 2 and C 2 HD[END_REF]. However, because the intensity of the saturated absorption signal is a nonlinear function of the pressure conditions, it is of crucial importance to properly model the entire resonance shape. Astonishingly, the sensitivity of a nonlinear detection (SA) can be comparable to that of a linear absorption detection [START_REF] Ma | Optical heterodyne spectroscopy enhanced by an external optical cavity: Toward improved working standards[END_REF].

The modeling is developed in Sec. 2 with the help of the Appendix 5. It includes the revised formalism, dealing with a Gaussian beam (comprising a full and an approximated approach). It is followed by Sec. 3 devoted to the profile simulations targeting the two specific transitions under study, while the model validation is discussed in Sec. 4.

Modeling

Background

The present modeling is largely based on previous developments proposed in recent papers [START_REF] Dupré | Saturated absorption and crossover resonances in a high-finesse cavity: Formalism and application to the hyperfine structure of jet-cooled NO 2 by saturated-absorption cavity-ring-down spectroscopy[END_REF][START_REF] Dupré | Sub-Doppler noise-immune cavity-enhanced optical heterodyne molecular spectrometry modeling: from Doppler broadening to cross-sideband resonances[END_REF][START_REF] Dupré | Dipole saturated absorption modeling in gas phase: Dealing with a Gaussian beam[END_REF]. However, slight modification has been brought to the formalism to properly deal with the polarization of the EMF. The net effect is a reduction of the saturation degree by a factor 3 (for the same value of the dipole moment of the transition), because of the redefinition of the dipole moment of the sub-transitions.

In an optically thin medium, the scalar complex amplitude absorption coefficient (i.e., including the absorption and the dispersion contributions) can be defined as the ratio of the infinitesimal change of the total EMF irradiance along a distance dz (i.e., inside the infinitesimal volume dV = dS dz), over the total EMF irradiance [START_REF] Weissbluth | Photon-Atom Interactions[END_REF][START_REF] Bordé | Saturated absorption line shape: Calculation of the transittime broadening by a perturbation approach[END_REF][START_REF] Alekseev | Maximum resolving power of laser spectroscopy with use of frequency resonances[END_REF]. In the frequency domain (using the inverse Fourier transform convention, and the angular frequency ω [START_REF] Kammler | A First Course in Fourier Analysis[END_REF], assuming a cylindrical EMF (i.e., exhibiting an envelope independent of z [the direction of propagation], but with a radial dependence), it becomes for molecules moving at the velocity v z along the z axis, (ignoring the relativistic Doppler and the recoil effects [START_REF] Basov | Optical frequency standards[END_REF][START_REF] Hall | Direct optical resolution of the recoil effect using saturated absorption spectroscopy[END_REF][START_REF] Bagayev | Second-order Dopplerfree spectroscopy[END_REF]), by following Ref. [START_REF] Dupré | Dipole saturated absorption modeling in gas phase: Dealing with a Gaussian beam[END_REF] and by paying attention to the unidirectionality of the absorbed dipole EMF (in other words, by taking the mean value [factor 1/3])

dα m u m l (ω, δ) = d dz i ω dP * m u m l (ω, δ, r, z) ⊗ E (ω, r, z) dV c 0 E * (ω, r, z) ⊗ E (ω, r, z) dS dω , (1)
where

dP m u m l (ω, δ, r, z) = 0 dχ m u m l (ω, δ, r ) E (ω, r, z) , (2)
is the macroscopic polarization (each molecule is considered as a microscopic target of radiation) of the medium, dχ m u m l (ω, δ, r ) is the nonlinear susceptibility tensor of a blocked 2-level system (u, l ), E (ω, r, z) is the amplitude of the incident EMF, c is the usual speed of the light, δ = v z k is the Doppler shift (k is the EMF wavevector:

k = k z = ω/c, A C C E P T E D M A N U S C R I P T k x = k y = 0)
, 0 is the vacuum permittivity, ⊗ represents a convolution product, and represent a product operator which depends on the rank of the development ( * denotes a complex conjugated quantity). The molecular energy levels can be characterized by their energy ω u and ω l (with ω l < ω u ), and by the total angular momentum F u and F l , respectively. Actually, each of the energy level has a Zeeman multiplicity, 2F u + 1 and 2F l + 1, which is associated with the degenerated magnetic sub-levels in the absence of static magnetic and electric field, denoted by the quantum numbers m u and m l , respectively.

In the linear regime of absorption, Eq. ( 1) can be greatly simplified if dχ m u m l (ω, δ, r ) is independent of the spacial coordinates (i.e., independent of the radial variable r , here) since the susceptibility becomes independent of the EMF in an homogeneous medium. Similar conditions are recovered inside an infinitesimal volume associated with a given value of r . Hence, the total absorption can be obtained by summing up all the "infinitesimal" absorption volumes with the relevant weight. Under theses conditions, Eq. ( 1) becomes:

dα m u m l (ω, δ) = i ω dP * m u m l (ω, δ, r, z) ⊗ E (ω, r, z) 2πr dr 0 c E * (ω, r ) ⊗ E (ω, r ) dω 2πr dr , (3)
where E (ω, r ) denotes the frequency and radial dependence of E (ω, r, z). dP * m u m l (ω, δ, r ) will be obtained by solving the density matrix equations.

Liouville Equation, Density Matrix

For a blocked 2-level molecular system (|u, F u , m u 〉, |l , F l , m l 〉), the Liouville Equation [START_REF] Fano | Description of states in quantum mechanics by density matrix and operator techniques[END_REF][START_REF] Meystre | Elements of Quantum Optics[END_REF] can be derived in the frequency domain. We obtain the master equation

iωρ (ω) = - i H (ω) ⊗ ρ (ω) -ρ (ω) ⊗ H (ω) + Λ (ω) -Γρ (ω) , (4)
where the convolution product has been extended to the matrix product (i.e., to the matrix element products), and where a dipole interaction V dip (ω) is assumed (actually, this formalism can be easily extended to quadrupole interaction):

H (ω) = 2πH 0 δ (ω) + V dip (ω) , (5)
where δ (ω) is the usual Dirac function, H 0 is the (static) molecular Hamiltonian in the absence of external EMF, and where Λ (ω) is the source term (at the first order of the development, this term is constant in the absence of population modulation), and Γρ (ω) is the relaxation matrix.

Assuming that the total coherence results from the sum of two counter-propagating coherences (with δ ≥ 0), after a first factorization, the matrix elements become

ρ m u m u (ω, δ) = -1 -ω+iγ u ρ m l m u (ω, δ) ⊗ Ω m u m l (ω) -ρ m u m l (ω, δ) ⊗ Ω m l m u (ω) + Λ u mu (ω,δ) γ u -iω ρ m l m l (ω, δ) = 1 -ω+iγ l ρ m u m l (ω, δ) ⊗ Ω m l m u (ω) -ρ m l m u (ω, δ) ⊗ Ω m u m l (ω) + Λ l m l (ω,δ) γ l -iω ρ m u m l (ω, δ) = ρ (+) m u m l (ω, δ) + ρ (-) m u m l (ω, δ) = L ul (ω, δ) ρ m l m l (ω, δ) -ρ m u m u (ω, δ) ⊗ Ω m u m l (ω) , (6)
with

L ul (ω, δ) = L (+) ul (ω) + L (-) ul (ω) , (7)

A C C E P T E D M A N U S C R I P T

and where the population relaxation rates γ u and γ l are assumed identical for all the sub-levels |u, F u , m u 〉 and |l , F l , m l 〉, respectively, while the polarization decay rate γ ul is independent of the sub-levels and of the molecule velocity. The source terms Λ u m u (ω, δ) and Λ l m l (ω, δ) are also assumed identical for each magnetic sub-level (isotropic space).

The term L (±) ul (ω) (independent of the Zeeman structure) characterizes the molecular transition (i.e., a complex

[not normalized] Lorentzian function, or CLF) frequency shifted by the Doppler effect δ (ω ul > 0)

L (±) ul (ω) = 1 -ω -ω ul ± δ + iγ ul = -L (∓) * l u (-ω) . ( 8
)
The dipole angular Rabi function Ω m u m l (ω) has been introduced by extension of the usual angular Rabi characteristic frequency. It has the same r dependence as E (ω), i.e.,

Ω m u m l (ω, r ) = 1 u, F u , m u - → µ • * l , F l , m l E (ω, r ) = µ m u m l E (ω, r ) , (9)
where the sub-transition dipole moment operators µ m u m l between two sub-levels have been introduced. They are linked by the reduced matrix element (Wigner-Eckart theorem [START_REF] Messiah | Quantum Mechanics[END_REF][START_REF] Zare | Angular Momentum[END_REF]) such as

µ m u m l 2 =    F u 1 F l -m u ∆m m l    2 ψ u F u - → µ ψ l F l 2 =    F u 1 F l -m u ∆m m l    2 µ ul 2 , (10)
where

µ ul 2 = S ul µ 2
band (S ul is the usual line strength factor, and µ 2 band is the squared molecular dipole moment of the vibration band considered, µ ul is a symmetry signed quantity [START_REF] Bernath | Spectra of Atoms and Molecules[END_REF][START_REF] Bunker | Molecular Symmetry and Spectroscopy[END_REF]), denotes the direction of the dipole EMF in the space-fixed referential frame, "." denotes the scalar (or inner) product, and where the Wigner three-j coefficients have been brought in. Moving from the Cartesian coordinates to the irreducible spherical components [START_REF] Fano | Irreducible Tensorial Sets[END_REF] provides the sub-transition selection rule value through ∆m (∆m = m u -m l = 0, ±1) associated with the EMF polarization component .

For a given polarization of the EMF, the Wigner symbol sum up according to

m u m l    F u 1 F l -m u ∆m m l    2 = 1 3 , (11)
which allows linking to the usual line strength S ul . Indeed, S ul [Eq. ( 10)] depends on the molecular Hamiltonian and on the internal couplings.

The degeneracy of the magnetic sub-levels, the absence of coupling and coherence between the sub-levels belonging to the same level, the electric transition selection rule between the magnetic components ∆m = 0 or ±1, and the assumption of a unique relaxation rate for all the sub-levels belonging to a given energy level, allow us to deal with a quasi two-level system, i.e., to ignore the structure of a multiple (2F u + 1) × (2F l + 1)-level systems.

Nonlinear Susceptibility for two counter-propagating Beams

As demonstrated previously, assuming dP

* (ω, δ) = Trace ρ * (ω, δ) µ , [with ρ * m u m l (ω, δ) = ρ m l m u (-ω, δ)],
the complex absorption can be obtained by calculating

A C C E P T E D M A N U S C R I P T dα m u m l (ω, δ) = -i µ m u m l 2 c 0 ω ∞ 0 L ul (ω, δ) × ∆ρ m l m u (ω, δ, r ) ⊗ E * (ω, r ) ⊗ E (ω, r ) 2πr dr E * (ω, r ) ⊗ E (ω, r ) dω 2πr dr dδ, (12)
where, at the first order of the limited development

∆ρ m l m u (ω, δ, r ) =∆ρ (0) ul (δ) 2πδ (ω) + L ul (ω) L ul (ω, δ) Ω * m u m l (ω, r ) ⊗ Ω m l m u (ω, r ) -cc , ( 13
)
where we can define the quantities

∆ρ (0) m l m u (ω, δ) = Λ l m l (ω, δ) γ l + iω - Λ u m u (ω, δ) γ u + iω , (14)
and ∆ρ (1)

m l m u (ω, δ) = L ul (ω) L ul (ω, δ) ∆ρ (0) m l m u (ω, δ) ⊗ Ω * m u m l (ω) ⊗ Ω m l m u (ω) -cc , (15)
and where

L ul (ω) = 1 2 1 -ω + iγ l + 1 -ω + iγ u . ( 16
)
The complex Lorentzian shape L ul (ω, δ) is not the Fourier transform of a real function. Hence, to recover a function which infers the "reality" of the molecular system, we introduce the sum function

L ul (ω, δ) = L l u (ω, δ) = L ul (ω, δ) + L * l u (ω, δ) . ( 17
)
The double integral at the denominator of Eq. ( 12) is easy to obtain, it is the total power of the EMF (disregarding the factor c 0 ).

Following Eq. 14, for a static, isotropic and uniform medium (i.e., independent of r and z), the zero order population difference can be set to

∆ρ (0) m l m u (ω, δ) = 2π∆ρ (0) m l m u (δ) δ (ω) , (18)
where ∆ρ (0) m l m u (δ) is provided by the usual normalized Maxwell-Boltzmann distribution

∆ρ (0) m l m u (δ) = g ∆ D (δ) g ns e - ω l k B T -e - ωu k B T Q (T ) N 0 T, p dδ = g ∆ D (δ) N m u m l T, p dδ, ( 19
)
where g ∆ D (δ) is the usual normalized velocity distribution (independent of the collision rate)

g ∆ D (δ) = ln 2 π 1 ∆ D e -ln 2 δ ∆ D 2 , (20)
and where

Q(T ) = i g ns i (2F i + 1) e -E i kT , (21)
is the partition function, N 0 T, p is the number density of absorbing molecules, and ∆ D is the Doppler half-with at

A C C E P T E D M A N U S C R I P T
half-maximum (HWHM). The nuclear spin statistic (degeneracy) g ns has been introduced. It only includes possible degeneracy due to the internal molecular structure (like due to the degenerated nuclear spin statistic, actually, this factor may require more attention because the degeneracy may be easily lifted).

Finally, the coefficient of absorption in the frequency domain is obtained by summing-up the individual contributions

α ul (ω) = m u m l α m u m l (ω) = m u m l dα m u m l (ω, δ) = m u m l ∞ 0 dα m u m l (ω, δ, r ) 2π r dr. ( 22
)

Application to a Modulated EMF

By keep using the notation which has been previously introduced in the frame of the slowly varying envelope approximation (SVEA), the EMF resulting from two counter propagating waves can be defined by

E (ω, r, z) = κ E κ (ω, r ) e ik κ z + e -ik κ z = 1 2 κ E 0 κ (ω -ω κ , r ) + E 0 κ (ω + ω κ , r ) e ik κ z + e -ik κ z =2πE 0 (ω, r ) ⊗ κ J κ (ξ) δ ω -ω 0 -κω rf + δ ω + ω 0 + κω rf cos k 0 + κk rf z ( 23
)
where E 0 (ω, r ) is a symmetric EMF profile (it is the Fourier transform of a real function) which only spreads over low frequencies. This allows dealing with the rotating wave approximation (RWA) if this spectral extension of the EMF is less than the modulation frequency (here at the radio-frequency ω rf ). In the case of a phase modulated EMF, the quantities J κ (ξ) are the first kind Bessel coefficients, parametrized by the modulation index ξ such as κ J 2 κ (ξ) = 1, and

J -κ (ξ) = (-1) κ J κ (ξ)
. ω 0 ω rf is the carrier (optical) frequency, and

k rf = ω rf /c.
The analytical solutions of the complex absorption are obtained by calculating the key quantity [ωL ul (ω, δ) × E * (ω, r )]⊗ E (ω, r ) for different shapes of the EMF of interest (the calculation of this quantity is available in the Appendix section).

We will perform the calculation for a modulated monochromatic EMF before focusing on a Gaussian EMF, as well as we will first examine the linear responses before dealing with the nonlinear responses.

It is worth commenting onthe product L ul (ω, δ) × E * (ω, r ) under FM conditions. The carrier frequency and the sidebands of the EMF comprises a pseudo optical frequency comb (i.e., limited to a few teeth). It acts as a sampler of the absorbing line shape L ul (ω, δ) at frequencies ω 0 ± κω rf constituting of a new pseudo optical frequency comb at the repetition rate ω rf .

The convolution of the two pseudo optical frequency combs with the same repetition rate gives arise to a new frequency comb, still with the same repetition rate, but centered at the origin of the frequencies (i.e., around ω = 0).

We will so-call DC components the spectrum generated around the frequency origin, and harmonic components, the spectra generated at the harmonic frequencies ∆κ ω rf (∆κ > 0), keeping a generic formalism. Actually, it is easy to deal with the DC terms, since they are similar to the case of a non modulated EMF: the unique contribution being replaced by a sum of all the contributions sampling the absorbing line shape (i.e., at each value of ∆κ ω rf ) weighted by the Bessel function coefficients.

A C C E P T E D M A N U S C R I P T

The harmonic beating notes, around ∆κ ω rf = 0, can be experimentally detected by employing a photodetector with a suitable bandwidth, and by frequency shifting the signal to the DC range by "multiplying" with the local oscillator frequency used to modulate the EMF. This provides quasi stationary responses.

Linear absorption

As it has been previously demonstrated, the linear absorption car be obtained by calculation Eq. ( 12) with the help of Eq. (B7).

m u m l (ω) = -i π µ m u m l 2 N m u m l T, p c 0 1 π ln 2 π ω 0 ∆ D ∞ -∞ κ≷0 J 2 κ (ξ) ω 0 + κω rf e -ln 2 δ ∆ D 2 ω ul -ω 0 -κω rf + δ -iγ ul dδ 2πδ (ω) , (24)
The spectral extension of α

(mono 0 )
m u m l (ω) is strictly limited to ω = 0. Hence, the temporal response is trivial: a strict stationary behavior which can be experimentally obtained through a standard DC transimpedance photodetector. By inserting the usual Faddeeva function (or plasma dispersion function) w erf [START_REF] Wells | Rapid approximation to the Voigt/Faddeeva function and its derivatives[END_REF][START_REF] Abramowitz | Handbook of mathematical functions with Formulas, Graphs, and Mathematical Tables[END_REF] (area [real part] normalized by the factor ∆ D π/ ln 2 1 ), it becomes the dependence versus the source carrier frequency ω 0 α

DC mono 0 m u m l (ω 0 ) = π µ m u m l 2 N m u m l T, p c 0 ln 2 π 1 ∆ D κ≷0 J 2 κ (ξ) ω 0 + κω rf w erf ω 0 + κω rf -ω ul -iγ ul ∆ D / ln 2 . ( 25
)
Now, considering a Gaussian EMF with a beam waist size w 0 . At a given impact parameter b [START_REF] Hall | The Linewidth Problem in Laser-Saturated Molecular Absorption[END_REF][START_REF] Hall | Sub Doppler spectroscopy, methane hyperfine spectroscopy, and ultimate resolution limits[END_REF], and for a mean transit-time rate

γ tr = ∆ D c
2 πw 0 ω 0 (see Sec. II.D.2.b of Ref. [START_REF] Dupré | Dipole saturated absorption modeling in gas phase: Dealing with a Gaussian beam[END_REF])

E 0 ω, b, γ tr = 2 π ln 2 γ tr E b e -ln 2 ω γ tr 2 = 2 π ln 2 γ tr E 0 e -b w 0 2 e -ln 2 ω γ tr 2 , (26)
one gathers

α DC G 0 m u m l ω, b, γ tr = π µ m u m l 2 N m u m l T, p c 0 e -2 b w 0 2 ln 2 γ 2 tr e -ln 2 ω 2γ tr 2 × ∞ -∞ g ∆ D (δ) κ≷0 J 2 κ (ξ) ω 0 + κω rf w * erf ω/2 -δ ul κ -δ + iγ ul γ tr / 2 ln 2 + w erf ω/2 + δ ul κ + δ -iγ ul γ tr / 2 ln 2 dδ, ( 27
)
where

δ ul κ = ω ul -ω κ = ω ul -ω 0 -κω rf = δ ul -κω rf . ( 28
)
It is worth noting that α

DC G 0 m u m l (ω, b)
is the sum of products comprising two functions of ω, i.e., a Gaussian function with a HWHM 2 γ tr centered at ω = 0 multiplied by a more complex function which is actually the convolution

1 w erf (z) = i π ∞ -∞ e -t 2 dt z-t , with Im (z) > 0 A C C E P T E D M A N U S C
Hence, by inverse Fourier transforming, we carry out

α DC G 0 m u m l ω 0 , b, γ tr = π µ m u m l 2 N m u m l T, p c 0 e -2 b w 0 2 1 γ tr 2 ln 2 π × ∞ 0 g ∆ D (δ) κ J 2 κ (ξ) ω 0 + κω rf w erf -δ ul κ + δ + iγ ul γ tr / 2 ln 2 + w * erf δ ul κ + δ -iγ ul γ tr / 2 ln 2 dδ. ( 30
)
In the approximation ω 0 + κω rf ω 0κω rf ω 0 , the frequency term ω 0 + κω rf can be factorized, and one obtains

α DC G 0 m u m l ω 0 , b, γ tr = π ω 0 c 0 g ns N 0 T, p e - ω l k B T -e - ωu k B T Q (T ) S ul    F u k F l -m u ∆m m l    2 e -2 b w 0 2 1 γ tr 2 ln 2 π × ∞ 0 g ∆ D (δ) κ J 2 κ (ξ) w erf -δ ul κ + δ + iγ ul γ tr / 2 ln 2 + w * erf δ ul κ + δ -iγ ul γ tr / 2 ln 2 dδ. (31)
From this formula it is easy to obtain the usual complex coefficient of absorption α DC G 0 ul ω 0 , b, γ tr , because the sum rule applies on the three-j coefficients.

The full linear response is obtained by calculating the integral

α DC G 0 m u m l ω 0 , γ tr = 2 π 2 w 0 ∞ 0 e -2 b w 0 2 α DC G 0 m u m l ω 0 , b, γ tr db. ( 32
)
Harmonic Components

. By using the same approach than that used in the previous paragraph (i.e., for a Gaussian beam identified by a mean transit-time rate, see 2.4.1), but now considering the frequency cross terms, the different harmonics can be obtained by combining and summing the relevant terms of Eq. (B7). Each harmonic term, or beating frequency (∆κ ω r f = 0), is a sum over the indexes κ. For a given value of ∆κ > 0, (a single value of δ is used for sake of compactness), four coupling terms, (i.e., κ, κ ± ∆κ, and κ ± ∆κ, κ), can be first derived under the generic following shape

A C C E P T E D M A N U S C R I P T ωL ul (ω, δ) × E * (ω, b) ⊗ E (ω, b) = 4π 2 ln 2 E 2 b γ 2 tr κ≥∆κ-1 J κ (ξ) J κ+∆κ (ξ) cos ∆κ k rf z ×      ω κ e -ln 2 ω γ tr 2 -ω -δ ul κ + δ -iγ ul ⊗ e -ln 2 ω γ tr 2 ⊗ δ ω -∆κ ω rf - ω κ e -ln 2 ω γ tr 2 -ω + δ ul κ + δ + iγ ul ⊗ e -ln 2 ω γ tr 2 ⊗ δ ω + ∆κ ω rf + ω κ+∆κ e -ln 2 ω γ tr 2 -ω -δ ul κ+∆κ + δ -iγ ul ⊗ e -ln 2 ω γ tr 2 ⊗ δ ω + ∆κ ω rf - ω κ+∆κ e -ln 2 ω γ tr 2 -ω + δ ul κ+∆κ + δ + iγ ul ⊗ e -ln 2 ω γ tr 2 ⊗ δ ω -∆κ ω rf + (-1) ∆κ ω κ e -ln 2 ω γ tr 2 -ω -δ ul κ + δ -iγ ul ⊗ e -ln 2 ω γ tr 2 ⊗ δ ω + ∆κ ω rf -(-1) ∆κ ω κ e -ln 2 ω γ tr 2 -ω + δ ul κ + δ + iγ ul ⊗ e -ln 2 ω γ tr 2 ⊗ δ ω -∆κ ω rf + (-1) ∆κ ω κ-∆κ e -ln 2 ω γ tr 2 -ω -δ ul κ-∆κ + δ -iγ ul ⊗ e -ln 2 ω γ tr 2 ⊗ δ ω -∆κ ω rf -(-1) ∆κ ω κ-∆κ e -ln 2 ω γ tr 2 -ω + δ ul κ-∆κ + δ + iγ ul ⊗ e -ln 2 ω γ tr 2 ⊗ δ ω + ∆κ ω rf      . ( 33
)
The different convolution products involve a CLF, and two Gaussian functions. There are associated with a propagation term driven by cos ∆κ k rf z . These integrals can be calculated by introducing the Faddeeva function as follows, (noting that theses function are not area normalized)

ω κ+∆κ e -ln 2 ω γ tr 2 -ω -δ ul κ+∆κ ± δ -iγ ul ⊗ e -ln 2 ω γ tr 2 = - i 2 ω κ+∆κ e -ln 2 ω 2γ tr 2 w erf ω/2 + δ ul κ+∆ ± δ + iγ ul γ tr / 2 ln 2 = W ± ω 2 + δ ul κ+∆κ + iγ ul . (34)
It is worth pointing that the frequency extension of these Faddeeva profile integrals are centered around ω = ±2 ω ulω 0 -∆κ ω rf , i.e., around twice the frequency detunings [see Eq. ( 28)]. The contribution due to the Doppler shifts at ±δ can be considered by setting

W ω 2 + δ ul κ+∆κ + iγ ul = W + ω 2 + δ ul κ+∆κ + iγ ul +W - ω 2 + δ ul κ+∆κ + iγ ul = W ω 2 + δ ul κ+∆κ + iγ ul +iW ω 2 + δ ul κ+∆κ + iγ ul , (35)
where the real W and imaginary W parts of the complex integrals have been introduced, and noting that W (z * ) = W * (-z).

It becomes for a given value of ∆κ

ωL ul (ω, δ) × E * (ω, b) ⊗ E (ω, b) = - 2iπ ln 2 E 2 b γ 2 tr e -ln 2 ω γ tr 2 κ≥∆κ-1 J κ (ξ) J κ+∆κ (ξ) cos ∆κ k rf z × W ω 2 + δ ul κ + iγ ul ⊗ δ ω -∆κ ω rf + W * ω 2 -δ ul κ + iγ ul ⊗ δ ω + ∆κ ω rf + W ω 2 + δ ul κ+∆κ + iγ ul ⊗ δ ω + ∆κ ω rf + W * ω 2 -δ ul κ+∆κ + iγ ul ⊗ δ ω -∆κ ω rf (-1) ∆κ W ω 2 + δ ul κ + iγ ul ⊗ δ ω + ∆κ ω rf + (-1) ∆κ W * ω 2 -δ ul κ + iγ ul ⊗ δ ω -∆κ ω rf (-1) ∆κ W ω 2 + δ ul κ-∆κ + iγ ul ⊗ δ ω -∆κ ω rf + (-1) ∆κ W * ω 2 -δ ul κ-∆κ + iγ ul ⊗ δ ω + ∆κ ω rf . ( 36
)
Experimentally, the amplitude of each harmonic term, beating at the frequency ∆κ ω rf , can be gathered by het-

A C C E P T E D M A N U S C R I P T
-J κ (ξ) [J κ-1 (ξ) + J κ+1 (ξ)] J κ (ξ) [J κ-1 (ξ) -J κ+1 (ξ)] 0 -2J 0 (ξ) J 1 (ξ) J κ (ξ) [J κ-1 (ξ) + J κ+1 (ξ)] J κ (ξ) [J κ-1 (ξ) -J κ+1 (ξ)] 2 -J κ (ξ) [J κ-2 (ξ) -J κ+2 (ξ)] J κ (ξ) [J κ-2 (ξ) + J κ+2 (ξ)] 2J 0 (ξ) J 2 (ξ) 0 -J κ (ξ) [J κ-2 (ξ) -J κ+2 (ξ)] -J κ (ξ) [J κ-2 (ξ) + J κ+2 (ξ)]
erodyning the collected signal with the desired harmonic of the purely monochromatic local oscillator (i.e., at the frequency ∆κ ω rf ). Formally, each harmonic can be obtained by a convolution (i.e., a multiplication in the time domain) by a voltage at the frequency of interest, which provides the desired DC (or quasi stationary) responses, and by filtering out the high frequencies. A way to calculate these components is to identify the terms proportional to π δ ω -∆κ ω rf + δ ω + ∆κ ω rf (for the terms in-phase signal), and to iπ δ ω -∆κ ω rfδ ω + ∆κ ω rf (for the signal in-quadrature). Hence, we obtain by factorization

ωL ul (ω, δ) × E * (ω, b) ⊗ E (ω, b) = - iπ ln 2 E 2 b γ 2 tr e -ln 2 ω γ tr 2 κ≥∆κ-1 J κ (ξ) J κ+∆κ (ξ) cos ∆κ k rf z × 2 W ω 2 + δ ul κ + iγ ul + (-1) ∆κ W * ω 2 -δ ul κ + iγ ul ⊗ δ ω -∆κ ω rf + (-1) ∆κ δ ω + ∆κ ω rf + W ω 2 + δ ul κ+∆κ + iγ ul + W * ω 2 -δ ul κ+∆κ + iγ ul + (-1) ∆κ W ω 2 + δ ul κ-∆κ + iγ ul + W * ω 2 -δ ul κ-∆κ + iγ ul ⊗ δ ω + ∆κ ω rf + (-1) ∆κ δ ω -∆κ ω rf + W ω 2 + δ ul κ+∆κ + iγ ul -W * ω 2 -δ ul κ+∆κ + iγ ul -(-1) ∆κ W ω 2 + δ ul κ-∆κ + iγ ul -W * ω 2 -δ ul κ-∆κ + iγ ul ⊗ δ ω + ∆κ ω rf -(-1) ∆κ δ ω -∆κ ω rf , ( 37
)
where we note the sign alternation with the parity of ∆κ.

Equation 36 comprises of 3 sum terms. For ∆κ = 1, the two first terms provide the in-quadrature FM signal involving only dispersion components (because it only involves the dispersion parts of the CLF), while the last term provides the in-phase FM signal involving only purely absorptive components of the CLF.

Considering the impact parameter (b) and the mean transit-time rate (γ tr ), Eq. ( 37) can be inserted into Eq. 12.

By weighting the Doppler shifts by the Maxwell-Boltzmann speed distribution, we obtain the linear absorption, i.e., the absorption (in-phase detection) and the dispersion (in-quadrature detection), where the coefficients X κ (ξ) (see Table 1) have been factorized

α FM G 0 m u m l ω, b, γ tr = ω 0 π µ m u m l 2 N m u m l T, p c 0 e -2 b w 0 2 ln 2 γ 2 tr e -ln 2 ω 2γ tr 2 ∞ -∞ g ∆ D (δ) κ≥∆κ-1 ω 0 + κω rf X κ (ξ) P w erf ω 2 -δ ul κ + δ + iγ ul γ tr / 2 ln 2 + P w erf ω 2 -δ ul κ -δ + iγ ul γ tr / 2 ln 2 + ω 0 -κω rf X -κ (ξ) P w erf ω 2 -δ ul κ + δ + iγ ul γ tr / 2 ln 2 + P w erf ω 2 -δ ul κ -δ + iγ ul γ tr / 2 ln 2 dδ, ( 38
)
and where P w erf (x) represents either the real part, or the imaginary of the Faddeeva function associated with, either the in-phase, or with the in-quadrature detection, respectively.

A C C E P T E D M A N U S C R I P T

The behavior of the terms w erf

α FM G 0 m u m l ω 0 , b, γ tr = π µ m u m l 2 N m u m l T, p c 0 0 e -2 b w 0 2 2 ln 2/π γ tr × ∞ 0 g ∆ D (δ) κ≥0 ω 0 + κω rf X κ (ξ) P w erf ω 0 -ω ul + κω rf + δ + iγ ul γ tr / 2 ln 2 + P w erf ω 0 -ω ul + κω rf -δ + iγ ul γ tr / 2 ln 2 + ω 0 -κω rf X -κ (ξ) P w erf ω 0 -ω ul + κω rf + δ + iγ ul γ tr / 2 ln 2 + P w erf ω 0 -ω ul + κω rf -δ + iγ ul γ tr / 2 ln 2 dδ. (39)
Naturally, by applying the integration previously mentioned [see Eq. ( 32)] to the relevant quantities, the coefficients

α FM G 0
m u m l ω 0 , γ tr become independent of the beam impact parameter.

As previously indicated, in the approximation ω 0 + κω rf ω 0κω rf ω 0 , the resonant frequencies can be easily factorized [see for example Eq. ( 31)] in the numeric implementation.

Saturated Absorption

The saturated absorption is obtained by calculating the term ∆ρ (1) m l m u (ω, δ) [see Eq. ( 15)], and then by plugging it into Eq. ( 12). It becomes

α (S A) m u m l ω, b, γ tr = -2i µ m u m l 2 N m u m l T, p c 0 0 × ∞ 0 g ∆ D (δ) ∞ 0 ω L ul (ω, δ) × E * ω, b, γ tr ⊗ E ω, b, γ tr + S m u m l (ω, δ, b) 2πr dr E * (ω, r ) ⊗ E (ω, r ) dω 2πr dr dδ, (40)
where

S m u m l ω, δ, b, γ tr = ω L ul (ω, δ) A ω, δ, b, γ tr ⊗ E * ω, b, γ tr ⊗ E ω, b, γ tr , (41)
with

A ω, δ, b, γ tr = 2 µ m u m l 2 L ul (ω, δ) × E * ω, b, γ tr ⊗ E ω, b, γ tr -cc L ul (ω) . (42)
A quick analysis shows that the quantity A ω, δ, b, γ tr is centered around ω = ∆κω rf , its spectral extension also finger prints the term S m u m l ω, δ, b, γ tr . Actually, A ω, δ, b, γ tr is the key quantity to calculate the saturation, and it requires special attention. As we have discussed about the linear absorption, the quantity L ul (ω, δ) × E * ω, b, γ tr ⊗ E ω, b, γ tr is a sum of cross terms, and of diagonal terms. Theoretically, we should deal with all the terms, but the full development becomes too cumbersome to be fully treated here. Indeed, if we focus only on the cross terms, the spectral extension is limited around ω ∼ 0, i.e., there are no beating terms around ω ∼ ±κω rf . Plugging the beating terms into S m u m l ω, δ, b, γ tr generates harmonics of a higher order [i.e., at ± (∆κ ± 1, ...) ω rf ]. At the opposite, ignoring these beating terms does not generate higher harmonics than those discussed in the frame of the linear absorption approximation. The net result of this approximation is the missing of some "interfering" term. Typically, the amplitude

A C C E P T E D M A N U S C R I P T

of the crossover resonances cannot be expected to be accurately predicted, specidically in the case of strong saturation.

Nevertheless, summing up the terms associated with the two opposite values of δ provides the desired resonances and crossover sideband resonances.

Before considering the Gaussian EMF, we will first deal with a pseudo-monochromatic EMF, i.e., a Dirac EMF exhibiting a radial dependence in r (instead of the impact parameter b previously considered). This r dependence can be employed to weight the EMF (i.e., to weight the saturation parameter, vide infra).

Case of a pseudo-monochromatic EMF

. This case is easily mediated by setting r = b in the previous equations (40 to ( 42)), and by ignoring the transit-time rate concept. The weighting factor becomes exp -2r 2 /w 2 0 , while E (ω) = 2π E 0 δ (ω). By plugging this latter expression into Eq. ( 42), and then, into the different constituting expressions, it becomes for A(ω,δ,r ) [with the use of Eq. (A1)]

A (mono) (ω, δ, r ) = κ J 2 κ (ξ) Ω (r ) 2 m u m l L 0 -δ ul κ , δ -cc L ul (0) 2 cos 2 κk rf z 2πδ (ω) , (43)
where the Rabi frequency follows the usual definition

Ω (r ) m u m l = E r µ m u m l = E 0 µ m u m l e -r w 0 2 . ( 44
)
Hence, we can rewrite Eq. ( 43) as

A (mono) (ω, δ, r ) = -2 κ J 2 κ (ξ) Ω (r ) 2 m u m l T 1 γ ul    1 δ ulκ +δ γ ul 2 + 1 + 1 δ ulκ -δ γ ul 2 + 1    2 cos 2 κk rf z 2πδ (ω) = -2 κ J 2 κ (ξ) s m u m l δ ul κ , δ, r 2πδ (ω) , (45)
where

T 1 = γ u +γ l 2 γ u γ l
is the usual longitudinal relaxation time.

By plugging A (mono) (ω, δ, r ) into Eq. ( 41), and by averaging the z dependence, the nonlinear contribution becomes

S (mono) m u m l (ω, δ, r ) = 1 -2 κ J 2 κ (ξ) s m u m l δ ul κ , δ ω L ul (ω, δ) × E * (ω, r ) ⊗ E (ω, r ) . (46)
By inserting Eq. ( 46) into Eq. ( 40), the linear term can be factorized. Assuming that a similar development at higher orders can be extrapolated, we can rearrange the expression to recover the usual form of the saturation absorption.

As previously stated, the stationary complex absorption can be obtained thanks to the multiplication by the Dirac function δ (ω). Thus, one obtains the stationary response versus the source carrier frequency ω 0 . It becomes for the DC component of the complex absorption

α (DCSA mono ) m u m l (ω 0 , r ) = -2i µ m u m l 2 N ul 3 0 c 2 ln 2 1 ∆ D 1 π e -2 r w 0 2 w 0 × ∞ 0 κ 1 ω 0 -ω ul -κω rf -δ+iγ ul + 1 ω 0 -ω ul -κω rf +δ+iγ ul J 2 κ (ξ) ω 0 + κω rf e -ln 2 δ ∆ D 2 dδ 1 + 2 κ s m u m l ω ul -κω rf -ω 0 , δ e -2 r w 0 2 , ( 47
)
where the absorption (imaginary part) and the dispersion (real part) contributions are clearly identifiable in the CLF.

A C C E P T E D M A N U S C R I P T

Likewise, the FM harmonic components can be deduced

α (FMSA mono ) m u m l (ω 0 , r ) = 2 2 π 1 w 0 e -2 r w 0 2 ∞ 0 g ∆ D (δ) α FM mono 0 m u m l (ω 0 , δ) dδ 1 + 2 κ s m u m l ω ul -κω rf -ω 0 , δ e -2 r w 0 2 , (48)
where

α FM mono 0 m u m l (ω 0 , δ) = π µ m u m l 2 N m u m l T, p c 0 0 × 1 π κ≥0 ω 0 + κω rf X κ (ξ) P L (+) ω 0 + κω rf + P L (-) ω 0 + κω rf ( 49
)
+ ω 0 -κω rf X -κ (ξ) P L (+) ω 0 -κω rf + P L (-) ω 0 -κω rf ,
where L (±) ul (ω) are defined by Eq. 8, and where the coefficients X κ (ξ) are defined in Table 1. A similar formula is available in Ref. [START_REF] Dupré | Sub-Doppler noise-immune cavity-enhanced optical heterodyne molecular spectrometry modeling: from Doppler broadening to cross-sideband resonances[END_REF] It is worth noting that the expression of the saturation term (i.e., the denominator of the fraction) is the same for all the components, i.e., the DC and the FM terms. This results from the approximation mentioned at the end of paragraph 2.4.2.

Case of a Gaussian EMF

. A similar approach to that previously carried out can be handled to calculate the nonlinear contribution A (G) ω, δ, b, γ tr by carefully taking into account the beam impact parameter. By plugging Eq. ( 26) into Eq. ( 42) [with the help of Eqs. (A1) and (A2)], one obtains (for a given value of the Doppler shift)

A (G) ω, δ, b, γ tr = -2iL ul (ω) ln 2 γ 2 tr 2π e -ln 2 ω 2γ tr 2 κ J 2 κ (ξ) Ω (b) 2 m u m l 2 cos 2 κ k rf z × Voigt ω/2 -δ ul κ + δ + iγ ul γ tr / 2 ln 2 + Voigt ω/2 + δ ul κ + δ + iγ ul γ tr / 2 ln 2 + Voigt ω/2 -δ ul κ -δ + iγ ul γ tr / 2 ln 2 + Voigt ω/2 + δ ul κ -δ + iγ ul γ tr / 2 ln 2 , ( 50
)
where the Voigt function is the real part of the Faddeeva function, and where Ω (b) m u m l is the Rabi frequency depending on the impact parameter b following the previous definition [see Eq. ( 44)].

A (G) ω, δ, b, γ tr is a sum of combinations of direct and convolution products of three functions of ω. Indeed, it is a sum of peaked functions, all centered around ω ∼ 0 (again, because of the approximation discussed paragraph 2.4.2).

As explained previously, we are incline to replace each term of A

m u m l ω, δ, b, γ tr = -4 T 1 γ tr 2π ln 2 × κ J 2 κ (ξ) Ω (b) 2 m u m l × Voigt δ ul κ + δ + iγ ul γ tr / 2 ln 2 + Voigt δ ul κ -δ + iγ ul γ tr / 2 ln 2 2π δ (ω) , (51)
which can be inserted in Eq. ( 41) to provide

S (G) m u m l ω, δ, b, γ tr = -4 T 1 γ tr 2π ln 2 L ul (ω, δ) E * (ω, z) ⊗ E (ω, z) × κ ω 0 + κω rf J 2 κ (ξ) Ω (b) 2 m u m l Voigt δ ul κ + δ + iγ ul γ tr / 2 ln 2 + Voigt δ ul κ -δ + iγ ul γ tr / 2 ln 2 . ( 52
)
This correction term acts as a saturation factor. It can be easily inserted into Eq. ( 40) as previously shown, to supply the complex saturated absorption spectrum. The stationary behavior (ω = 0) is obtained as usually for a given value of

ω 0 . It becomes α FM SA G m u m l ω 0 , b, γ tr = ∞ 0 g ∆ D (δ) α FM G 0 m u m l ω 0 , b, γ tr , δ 1 + 2 T 1 γ tr 2π ln 2 κ J 2 κ (ξ) Ω (b) 2 m u m l Voigt δ ulκ +δ+iγ ul γ tr / 2 ln 2 + Voigt δ ulκ -δ+iγ ul γ tr / 2 ln 2 dδ, (53)
where α FM G 0 m u m l ω 0 , b, γ tr , δ has be generalized to set for the linear components of the DC absorption, and for the absorption and dispersion FM terms (see Sec. 2.4.1). The DC linear component is similar to Eq. . [START_REF] Quinn | Practical realization of the definition of the metre, including recommended radiations of other optical frequency standards[END_REF] without the δ integration, while the FM terms are similar to Eqs. (39) (without the δ integration).

Finally, the last calculations require for summing the individual transitions all over the m u and m l Zeeman sublevels.

The intensity of the sub-transition depends on the EMF polarization., Integrating only over the impact parameter, b is a valid approximation if a mean transit-time is considered. An additional integral is required for taking into account the distribution of the transit-time rate according to

α FM SA G ul (ω 0 ) = m u m l ∞ 0 g γ t r γ tr ∞ 0 α FM SA G m u m l ω 0 , b, γ t r db dγ tr , ( 54
)
with the weighting function

g γ tr γ tr = π 2 γ tr γ 2 tr e -π γ tr 2γ tr 2 . ( 55
)
Actually, these sums upand integrations apply to all the frequency harmonics as well as to the DC component (i.e., to the harmonics 0). In the frame of the current approximation, the analytical shape of the saturation factor is the same for all the harmonic terms. A careful exam shows that the (crossover) resonances potentially arise around δ = 0, δ = ω r f /2, etc.. The crossover sidebands resonances exhibit amplitudes proportional to the factor J 2 0 (ξ) + J 2 1 (ξ) for the resonances located at ω r f /2, while intuitively, the amplitude of these resonances are expected to be proportional to the product J 0 (ξ) J 1 (ξ) since they result from two Doppler shifted counter-propagating EMFs with amplitude J 0 (ξ)

A C C E P T E D M A N U S C R I P T

and J 1 (ξ). Nevertheless, the Doppler envelop is really proportional to the product J 0 (ξ) J 1 (ξ). The correct crossover resonance amplitude could be calculated by using a pump-probe model. On the another hand, we know that the amplitude of the resonances is barely sensitive to the saturation factor [START_REF] Dupré | Sub-Doppler noise-immune cavity-enhanced optical heterodyne molecular spectrometry modeling: from Doppler broadening to cross-sideband resonances[END_REF] (see also Sec. 3.1.1) which minimizes any resonance amplitude inaccuracy of the current model.

Implementation

The implementation of the different models of the saturated absorption has been performed on a personal computer by using the gcc language, and on a graphical interface based on perl and Tk languages. The integration routines are the simulations involving only two integrations, it is relevant to compare these latter two "models" (actually the CPU demand strongly depends on the saturation factor). The fitting procedures required for discussing the simulations are based on gnuplot [START_REF]Gnuplot: A portable, multi-platform, command-line driven graphing utility[END_REF], which provides the Faddeeva function (libcerf).

Simulations

Aiming to validate the modeling described in section 2, numerous simulations have been performed. We targeted the R e (0) para-transition belonging to two perpendicular combination bands of the polyad ∆P 11 of acetylene C 2 H 2 (transition symmetry Π u ← Σ + g ) [START_REF] Herman | Vibrational spectroscopic database on acetylene, X 1 Σ + g ( 12 C 2 H 2 , 12 C 2 D 2 , and 13 C 2 H 2 )[END_REF]; one is located at ∼ 7143.8289 cm -1 and the other one is located at ∼ 7230.61313 cm -1 .

Regardless of the intra-polyad coupling, effective electric dipole moments equal to ∼ 0.912 mD and to ∼ 0.515 mD are assumed, respectively [START_REF] Jacquemart | The IR spectrum of 12 C 2 H 2 : Line intensity measurements in the 1.4 µm region and update of the databases[END_REF]. These dipole moments remain relatively large, and the saturation conditions can be reached at rather low incident power. Two different Rabi frequencies have been initially chosen to compare the resonance shape for a large range of collision rates. Actually, following the above definition, the mean Rabi frequency is defined by

Ω Rabi = µ ul h 1 w 0 P c 0 J 0 (ξ) , ( 56
)
where P is the unidirectional beam power, and µ ul is the dipole moment of the transition derived from the line intensity.

Finally, in Sec. 3.2, the power broadening of one of the transitions is shown for a given collision rate. However, the discussion applies straightforward to the other transition.

We compare the four types of simulations: type #1 involving a strictly monochromatic EMF, or plane wave (i.e., involving a single integration over the Doppler shift), type #2 involving a monochromatic EMF weighted by a Gaussian radial extension (two integrations), type #3 involving a Gaussian EMF controlled by the impact parameter and the mean transit-time (two integrations), and type #4 involving a Gaussian EMF with a Maxwellian distribution of the transit-time (three integrations). We could assume that the inelastic collisional self broadening is controlled by the gas pressure featured by a self-pressure broadening coefficient equal to 0.197 cm -1 /atm. at room temperature as suggested by HITRAN [START_REF] Jacquemart | The spectrum of acetylene in the 5-µm region from new line-parameter measurements[END_REF][START_REF] Jacquemart | The IR spectrum of 12 C 2 H 2 : Line intensity measurements in the 1.4 µm region and update of the databases[END_REF][START_REF] Rothman | The HITRAN 2012 molecular spectroscopic database[END_REF], but the value of this coefficient is probably incorrect in the low pressure regime probed by sub-Doppler spectroscopy (see Introduction). Hence, we decided to directly parametrize the current simulations by the collision rate, and not by the gas pressure. We will also mainly focus on the FM-dispersion signal because the FM-absorption signal does not exhibit any resonance at the center of the transition. All the broadening values are provided at Half-Width at Half-Maximum (HWHM). The EMF power can reach relatively high values which is relevant when optical cavities are used. Nevertheless, the pump and probe beams are assumed to be identical, with a cylindrical shape, i.e., without any dependence of the size along the z propagation direction. If the Rayleigh length of the EMF is large enough, the EMF curvature effects are marginal, minimizing the broadening by the front-wave [START_REF] Shimoda | Line Broadening and Narrowing Effetcs[END_REF]. Actually, away from the waist of a Gaussian beam, the transit-time increase compensates the broadening due to the curvature (any adiabatic rapid passages are ignored [START_REF] Adam | Rabi oscillations and rapid-passage effects in the molecularbeam CO 2 -laser stark spectroscopy of CH 3 F[END_REF]).

The calculated values of the absorption or dispersion are provided in cm -1 , the weak "absorbance" can be obtained by multiplying the y-scale of the plots by the equivalent absorption length.

3.1.

Transition R e (0), ν 1 + ν 3 + ν 1 4

Profile Analysis

The first series of simulations have been obtained for an estimated pressure of ∼ 1 Pa (see above), a regime where the resonance broadening due to the transit-time rate dominates the collision rate, and where the mean Rabi frequency is intermediate between the transit-time rate and the collision rate. Figure 1

Residual

Figure 2: Zoom on the central resonance of the FM dispersion signal of the transition R e (0) belonging to the combination band ν 1 + ν 3 + ν 1 4 (polyad ∆P = 11) of C 2 H 2 (see Fig. 1). The caption is identical to that of Fig. 1. The upper panel shows three different simulations (here the sole Doppler shift integration [in cyan] and the Doppler shift supplemented by a r integration overlap) while the lower panel shows the residual difference after fitting with a linear frequency slope and the imaginary part of a Faddeeva function (or the dispersive part of a area-normalized CLF, see Table 2). The fitted values are available in Table 2. The linear behavior (here a slope) holding the resonance is indicated by a dashed line. resonances at ±350 MHz (due to the two closest sidebands), and the crossover sideband resonances at ±175 MHz are clearly identifiable. To evaluate the effects of the different models which have been previously documented, the central resonance is targeted. Figure 2 shows the FM dispersion resonance, i.e., a zoom of Fig. 1. Since the Doppler-broadened profile is barely sensitive to the saturation conditions, the Doppler profile is almost independent of the model used under the current simulation conditions (i.e., ω rf ∆ D ).

To supply quantitative information, the profiles provided by the different models of simulations have been fitted, see Tab. Nevertheless, this slope is proportional to the number density.

The amplitude of the DC Doppler-broadened profiles depends on the simulation modeling. We observe a similar behavior which has already been reported in the previous study [START_REF] Dupré | Dipole saturated absorption modeling in gas phase: Dealing with a Gaussian beam[END_REF]: when the number of integration is increased, a desaturation of the transition can be observed in the following order: the integration over r (model #2) desaturating less than the integration over the impact parameter (model #3). The effect of the third integration (i.e., taking into account Standard errors are three times the standard deviation provided by gnuplot [START_REF]Gnuplot: A portable, multi-platform, command-line driven graphing utility[END_REF]. All the broadenings are given at half width at half maximum (HWHM). The data have not been weighted.

A C C E P T E D M A N U S C R I P T

A C C E P T E D M

A N U S C R I P T the distribution of transit-time rates) is not trivial on the Doppler-broadened profile for the simulation parameters of Figure 2. We previously reported that performing this integration brings an additional slight desaturation, but this effect cannot be observed here. Indeed, the transit-time rate integration increases the amplitude of the Lamb-dip while only the r integration decreases the width of the Lamb-dip. An accurate evaluation of the saturation is available in Table 2.

It is worthy to note that only the simulation involving a broadening by the transit-time rate can be reliably fitted by a Voigt profile. The Lamb-dips (DC) associated with the integrations of type #1 and #2 can only be satisfactory fitted by a standard Lorentzian function. However, the χ 2 coefficient (least square sum) attached to the triple integration (type #4) is clearly larger than that attached to the model involving only a mean transit-time rate. This reveals that the equivalent Lamb-dip profile can slightly deviate from a Voigt profile, at least under the current conditions. Actually, such mismatches are visible on the residual plots.

The analysis of the dispersion FM resonances is based on a Faddeeva profile. It shows a behavior slightly different from that just described (see Fig. 2). First of all, the single integration (#1) and the r integration (#2) show resonance profiles which cannot be distinguished, they can only be properly fitted by using the imaginary part of a Lorentzian function, i.e., they do not require for considering the Gaussian component of the imaginary part of the Faddeeva function which indeed cannot be determined. Again, considering the transit-time rate slightly reduces the amplitude of the resonance, and the Gaussian component of the Faddeeva profile could not be determined for the integrations which do not comprise a transit-time rate. The residual least square sum (χ 2 coefficient, see Table 2) usually show a very good match for the integrations of type #3 when fitting with a Faddeeva function (plot in green).

A comparison of the Lamb-dip profiles with the FM dispersion profiles shows that the fitted equivalent resonance linewidths are close (differences less than 12%). This demonstrates a good coherence between the different approaches.

Again, for the FM components, as previously mentioned for the DC component, the χ 2 coefficients show a slightly better match for the simulations based on a double integration (#3) than for those based on a triple integration (#4)

which may indicate that the proposed fitting profile (based on a Faddeeva function) is not optimum since we think that the model type #4 is the most accurate model in the frame of this present study.

Collision Rate Effects

In order to compare the behavior of the sub-Doppler resonances with the behavior of the Lamb-dip (DC), we simulated the FM dispersion resonances under a large range of collision rates (from 0.2 to 1000 kHz, i.e., an equivalent pressure range estimated from ∼ 0.003 to ∼ 20 Pa). For the transition considered here, the FM dispersion resonance is washed out below the considered pressure range (see sub-Fig. 5.b). We will only discuss simulations involving the transit-time broadening, either by assuming a mean transit-time rate, or by assuming a Maxwell-Boltzmann distribution of the transit-time rate. This analysis has been possible by fitting the resonance profiles as it is explained in Sec. 3.1.1.

We report, on the width of the two components (i.e., on the Lorentzian and Gaussian widths, see Fig. 3 and4), and on the integral value (see Fig. 5) of the imaginary part of the Faddeeva function obtained by fitting.

For the two simulation models, i.e. #3 and #4, the first noticeable behavior is the convergence of the component widths towards an identical asymptotic slope in the "high" pressure range, i.e., a pressure regime where the collision rates are higher than the mean transit-time rate. Actually, in this regime, both HWHMs converge approximately towards The intermediate pressureregime is a little bit more tricky to characterize: both components of the imaginary Faddeeva function (for both simulation types) exhibit a minimum which is controlled by the mean Rabi frequency, even if the convergence is only partial. In the collisional regime, when the amplitude of the resonance remains not

A C C E P T E D M

A N U S C R I P T extremely small (i.e., here, for collision rates less than 500 kHz), the Lorentzian component converges, more or less, towards the collision rate (inclined asymptote on both pictures), as already observed on the DC Lamb-dip [START_REF] Dupré | Dipole saturated absorption modeling in gas phase: Dealing with a Gaussian beam[END_REF]. The width of the Gaussian component remains below this slope, i.e., this component barely contributes to the effective resonance width (this behavior has already been noted). Actually, at "high" pressure, the resonance shapes incline to follow a derivative Gaussian function (i.e., a Dawson integral [START_REF] Press | Numerical Recipes in C[END_REF] 3 ).

There are multiple ways to analyze the amplitude of the resonances. Here, we adopted to discuss the integral of the normalized Faddeeva function, and not the peak-to-peak resonance amplitude, because the integral value is the relevant parameter, indeed. This integral value is proportional to the number density, or to the gas pressure if these two quantities are proportional. Figure 5 shows the number density normalized Faddeeva profile integral (panel a), and the non-normalized Faddeeva profile integral (panel b). First of all, it clearly appears that both simulation types (i.e., #3 and #4) provide very similar integral values, while this similarity is less trivial when considering the peak-to-peak amplitude because the resonance widths might differ. To summarize, the normalized integral is almost constant in the collision-free regime, while it drops off in the collisional regime. This demonstrates that there is an optimum collision regime to detect for saturation resonances, actually, around the mean transit-time rate. This optimum varies with the value of the mean Rabi frequency (i.e., with the value of the dipole moment, and of the power density of the impinging EMF).

Transit-time Broadening

Despite that no trivial experiment can be easily performed for providing evidence tothe behavior of the resonance width versus the transit-time rate (at the opposite of the collision rate which is easy to control), the transit-time rate will be considered as a possible experimental adjustable parameter. Actually, the beam diameter can be relatively well managed in the free space by arranging a set of optical lenses. However, the use of a cavity limits the possible adjustments of the transit-time rate. Furthermore, the value of transit-time rate is directly driven by the Doppler broadening in a bulk vessel (i.e., it depends on the molecular mass, and on the temperature of the species). For the current simulations, the beam power and the beam diameter are kept constant while the mean transit-time rate is widely varied (here, from 1 to 2000 kHz, see Fig. 6), in a range, which is of course impossible to cover with a single experimental setup. The collision rate is kept constant at 60 kHz.

The low transit-time rate regime is well-defined since it coincides with mean transit-time rates less than the collision rate. Both integration types of interest converge towards a dispersion Lorentzian shape of the resonances. The HWHM of the different components considered converge towards a value close to the mean Rabi frequency (∼ 141 MHz for ∼ 100 MHz), i.e., towards a value less than the sum of the collision rate and of the mean Rabi frequency. Above the chosen collision rate, the behavior is more problematic to describe, because the two integration types do not converge towards a close behavior. For the model of type #3, the Lorentzian component vanishes while the resonance shape is dominated by a Dawson integral. On the other hand, for the model type #4, both contributions remain similar, while the fitted width associated with the Dawson integral never passes over the width associated with the dispersion Lorentzian component. In summary, at low power, both models show an asymptotic behavior with a convergence value of the most relevant HWHM equal to approximately half the mean transit-time rate. On the hand, the narrowest component of the Faddeeva profile converges towards the collision rate.

A C C E P T E D

Discussion about the Simulations

In a first step, we mainly compared our predictive models developed for the FM dispersion resonance to two specific R e (0) transitions (exhibiting no Zeeman, neither hyperfine sub-structure) belonging to two different vibrational combination band of the ground electronic potential energy surface of C 2 H 2 , under two different amplitudes of the EMF. We primary wanted to validate the two types of simulation involving a transit-time associated with the Gaussian shape of the counter-propagating identical "saturating" EMFs. One model involves only a mean transit-time rate while the more elaborated model (and more CPU demanding) involves an integration over all the transit-time rates available in a Maxwellian gas. It appears that the two models behave slightly differently in the three dominant collision regimes:

the collision regime, the collision-free regime and the intermediate regime. frequency is lower than the mean transit-time rate which is the case for the two transitions analyzed here. This regime evidences the crucial role of the Rabi frequency when targeting the narrowest resonances. This narrowing effect appears more trivial when the two dispersion components of a Faddeeva profile can be fitted.

The collisional regime is easy to identify (the collision rate is larger than the mean transit-time rate), and it is readily to analyze and interpret, since it is also a regime where the collisions unsaturate the absorption. In this regime the Lorentzian component dominates the broadening with a width close to the collision rate (the fitted data approach an asymptote behavior, in purple color on the different figures) until it vanishes with the washing out of the resonances.

The collision-free regime is more complex to describe, because the two types of model involving two or three integrations provide different relative values of the two width components. This behavior can also be identified on the plot showing the fitted component widths versus the mean transit-time rate (see Fig. 6). Actually, this mismatch between the two integration types clearly demonstrates the role of a (Maxwellian) distribution of the transit-time rate.

Nevertheless, at very low collision rate, the saturation becomes very large, inducing an increase of the resonance width.

However, the high saturation degree induced by the weak collision relaxation points out the possible limit of the current model which may become irrelevant under these extreme conditions. Indeed, in the absence of solid experimental data, and despite of some strong suggestion, the elastic collisions have been ignored. This must be reconsidered, because in the low pressure regime, the dephasing rate (γ ul ) may dominates the inelastic collision rate (i.e., T 2 < T 1 ). Under these circumstances, the low pressure saturation effects are less pronounced than what we assumed here: the net result might be an absence of a power broadening at low pressure.

All over, the behavior of the FM dispersion resonance is not very different to that observed on the DC Lamb-dips as we discussed in a recent paper [START_REF] Dupré | Dipole saturated absorption modeling in gas phase: Dealing with a Gaussian beam[END_REF] for a single Rabi frequency, or in Table 2. However, here we did not consider an

A C C E P T E D M A N U S C R I P T
additional Lorentzian profile like we did in the previous analysis of the Lamb-dips to improve the residual least square sum, typically when the mean transit-time rate was not the dominant broadening factor. Nevertheless, on the analysis of the residual least square sum, and of the standard deviation, the simulated data issued from the integration type #3 exhibit a better agreement with the dispersion Faddeeva function than those issued from the integration type #4. This indicates that the fitting Faddeeva function chosen is not the optimum choice, because we think that the model involving three integrations must describe more accurately the resonance shape. Nevertheless, the unique mean transit-time rate approximation may remain a reasonable approximation if qualitative behaviors and resonance shapes are requested.

O. Axner and co-workers [START_REF] Axner | Sub-Doppler dispersion and noise-immune cavity-enhanced optical heterodyne molecular spectroscopy revised[END_REF] reported on the Lorentzian dispersion of the saturation resonances while here, we show that the resonance shape is not that simple, at least, the resonance shape exhibits a pressure and a transit-time rate dependence. Furthermore, the same group reported on the peak-to-peak amplitude of the dispersion resonances which is in agreement with the current simulation. Nonetheless, we were able to benchmark the relative amplitude of the resonance to the Doppler-broadened profile of our simulations with the experimental plots available in Fig. 4 of Ref. [START_REF] Axner | Sub-Doppler dispersion and noise-immune cavity-enhanced optical heterodyne molecular spectroscopy revised[END_REF].

It is also worth pointing out that the lack of discrimination between the contribution of each Zeeman sub-transitions, i.e., when these sub-transitions are just ignored, does not change very much the peak-to-peak amplitude of the FM resonances if the saturation factor is strong enough whatever the simulation type considered, the resonance width is mostly affected.

Conclusions

With the aim of optimizing the performances of setup implementing frequency modulated spectroscopy (FMS) in the saturation regime (rising Doppler-free resonance), we have developed a formalism based on the resolution of the Liouville equation (here limited to a 2-level system) in the frequency space, by keeping in mind feasible experiments.

Actually, this new development is a combination of the two formalisms recently published [START_REF] Dupré | Sub-Doppler noise-immune cavity-enhanced optical heterodyne molecular spectrometry modeling: from Doppler broadening to cross-sideband resonances[END_REF][START_REF] Dupré | Dipole saturated absorption modeling in gas phase: Dealing with a Gaussian beam[END_REF]. This improvement has been first performed by supplying analytical solutions (associated with a given Doppler shift), then by calculating numerical integrations, and finally by summing over the individual Zeeman sub-transitions. However, dealing with a Gaussian EMF introduces additional complexity. To avoid having to calculate too complex combination of direct and convolution products, we approximated a key product function by a Faddeeva function. This approach does not include any specific saturation coefficient since this one would result from multiple parameters (like the Rabi frequency or the active part of the laser beam). In the low pressure regime, we observed that the transit-time rate, defined in the framework of a Gaussian beam, substitutes more or less to the relaxation longitudinal rate 1/T 1 (as strictly defined for monochromatic beams).

The saturation effects on the energy level population difference have been supplied under a fractional form, following a usual limited development. Nevertheless, the number of interferences between the frequency components of the EMF (sidebands) have been limited by the use of an approximation preventing the propagation of the harmonics to higher orders. It results that the like-"saturation factor" is identical for all the harmonics components, i.e., for the detection at ω rf = 0 (DC), and at ω rf = 0, for both, the in-phase (absorption) and the in-quadrature (dispersion) FM Under the simulation conditions, we were able to establish three collision regimes, and to show the crucial role of the Rabi frequency. For example, if the value of the Rabi frequency is less than the mean transit-time rate, a narrowing effect is observed, i.e., the resonance width becomes less than the mean transit-time rate (approximately by a factor two). Despite that it may hurt the intuition, this super narrowing of the resonances have been experimentally observed in the early analyses of the saturated sub-Doppler absorption, and has been discussed in several works. This behavior which can be attributed to the dominant contribution (weighting) of the slow molecules (exhibiting the lowest transittime rates) at low pressure. This effect is slightly more pronounced for the model involving three integrations (i.e., an integration over the transit-time distribution). It is associated with the weighting of the saturation regime by the molecule velocity distribution.

We had the opportunity to compare the two most relevant models, and thus, to discuss the slight discrepancies between them. These discrepancies mainly concern the relative contribution of the two components of the Faddeeva profile (characterized by a Lorentzian and a Gaussian component) to the resonance broadening, in regimes where the mean transit-time rate dominates the collision rate, and the mean Rabi frequency.

Moreover, we discussed the respective roles of the mean transit-time, of the Rabi frequency (through the power broadening), and of the collision rate versus the modeling approaches. Nevertheless, we anticipate that this study is not a definitive contribution to the long history of the saturation resonances. For example, the domain of validity of the present model in term of equivalent coefficient of saturation, is not entire surrounded. This is peculiarly right in the low pressure regime where the predicted behaviors may not necessary match some earlier results obtained on methane, for example. Here, we clearly shown that the (very) weak pressures exhibit a broadening which has never been observed to the best knowledge of the author. The elastic collisions would counteract this effect.

Dealing with the collision rate can be approached by performing experiments under a large range of pressure.

Velocity changing collision phenomenons, as there are considered in Hartman/Tran model for example, have been ignored here, but they should be explored in the framework of the saturated absorption and of the low pressures.

  : i.e., considering only the transverse component of the molecule speed, the fast molecule prevail to the wing of the resonances while the slow molecules A C C E P T E D M A N U S C R I P T prevail to the central part of the resonance.

.

  Let us first consider the case of a "monochromatic" EMF [E 0 (ω) = 2πE 0 δ (ω)]. With the help of Eqs. 20 and (B6), it becomes [setting cos 2 (k κ z) = 1/2] α (mono 0 )

-ln 2 ω 2γ tr 2 = γ tr 2π ln 2 δ

 22 R I P Tproduct between a Doppler profile, and a w erf function. A simple analysis shows that in the Doppler dominant regime of interest (∆ D > γ tr , ∆ D > γ ul ), this convolution product is close to a Gaussian function centered on the frequency detunings ω 0κω rfω ul (i.e., it is maximal for ω 0 + κω rf = ω ul ) with the HWHM ∆ D . Hence, the overall resultant frequency spectrum spreads over ∼ 2 γ tr (centered at ω = 0) with an amplitude controlled by the detunings. The stationary response (ω = 0) can be obtained by multiplying by δ (ω). After rearranging the multiplicative term the Gaussian contribution becomes e (ω) .

ω 2 -δ ulκ +δ+iγ ul γ tr / 2 ln 2 ,

 22 which are resonant for ω 2 ω 0ω ul + κω rf have already been discussed in paragraph 2.4.1. Likewise, the stationary response can be deduced by multiplying by δ (ω) and inverse Fourier transforming

ω 2γ tr 2 /

 2 (G) ω, δ, b, γ tr by a sum of functions proportional to a Dirac function to calculate S (G)m u m l ω, δ, b, γ tr . This allows factorizing the term [ω L ul (ω, δ) × E * (ω, b)] ⊗ E (ω, b) inEq.[START_REF] Hall | Measurement of methane hyperfine structure using laser saturated absorption[END_REF], which offers a large convenience. Hence, we replace the normalized term expln 2 γ tr / ln 2/2π by a Dirac function δ (ω) in Eqs.[START_REF] De Labachelerie | Ultranarrow 13 C 2 H 2 saturated-absorption lines at 1.5 µm[END_REF] and[START_REF] Hall | Measurement of methane hyperfine structure using laser saturated absorption[END_REF]. Thus [setting 2 cos 2 (k 0 z) = 1], it becomesA C C E P T E D M A N U S C R I P TA(G)

  issued from the open source gsl library (implementing adaptive methods). The parallelization of the code is based on the OpenMP API[72]. The calculation of the Faddeeva function is based on the original algorithm provided by C. Benner[START_REF] Letchworth | Rapid and accurate calculation of the Voigt function[END_REF], acclimated to the C language2 . Because the simulations involving three integrations (excluding the calculation of the Faddeeva function) is approximately 100 times longer (i.e. ∼ 1 hour for 1000 sample point) than

Figure 1 :

 1 Figure 1: Doppler-broadened FM dispersion signal of the transition R e (0) belonging to the combination band ν 1 + ν 3 + ν 1 4 (polyad ∆P = 11) of C 2 H 2 (I = 0). Electric dipole moment: 0.912 mD, linear polarization, Doppler HWHM: 0.00863 cm -1 , collision rate: 60 kHz, number density: 2.5 × 10 14 cm -3 , line intensity: 5.985 × 10 -24 cm/mol. (approximate pressure: 1 Pa), directional beam power: 200 mW, beam waist size: 0.46 mm, mean transit-time broadening: 220 kHz, mean Rabi frequency: 97.8 kHz, frequency modulation: 350 MHz, modulation index: 0.4. The central and the two weak sideband resonances (at ±350 MHz), as well as the cross sideband resonances (at ±175 MHz), are visible.

  shows a large part of the FM dispersion Doppler-broadened profile, while the different narrow resonances, i.e., the central resonance, the A C C E P T E

2 .

 2 The (unmodulated) Lamb-dip profiles (so-called DC component) have been fitted by a Voigt profile (real part of a Faddeeva function) superimposed to a Gaussian profile accounting for the Doppler-broadened profile (not shown). The FM dispersion resonances (centered on the null detuning) have been fitted by using the imaginary part of a Faddeeva function superimposed to a linear slope accounting for the Doppler profile; on this range of detuning, it is a very good approximation of evaluating the amplitude of the Doppler absorption. As anticipated, the values of the slope do not change as a result of the quasi-invariance of the Doppler-broadened profile with the saturation degree.

10 - 21 †

 1021 The Doppler contribution is evaluated from the Doppler-broadened of the DC signal, or from the slope of the "dispersion" profile for the FM signal. ‡ Resonance integrals are deduced either from a standard Lorentzian function, from a Voigt profile, from a CLF, or from a Faddeeva function (all area normalized). * Lorentzian component of the dispersive contribution of a Faddeeva function (see note 1), or the dispersive component of a normalized CLF ( 1π 1 γ L +iω). * * The lack of Gaussian width results from the absence of convergence of the Voigt or Faddeeva profile (only a Lorentzian component could be deduced).

Figure 3 :Figure 5 :

 35 Figure 3: Analysis of the FM dispersion central resonance width of the transitionR e (0) belonging to the combination band ν 1 + ν 3 + ν 1 4 (polyad ∆P = 11) of C 2 H 2 located at ∼ 7143.8289 cm -1 : Lorentzian component (HWHM) of the dispersion (imaginary part) Faddeeva function versus the collision rate. Simulation conditions are available in the caption of Fig. 1. Here, two models (involving double and triple integrations) are shown. The inclined asymptote (dashed) indicates a HWHM equal to the collision rate. The error bars are equal to three times the standard deviation provided by the fit algorithm. The mean transit-time rate and the mean Rabi frequency are indicated.

Figure 7 :Figure 8 :

 78 Figure 7: Analysis of the FM dispersion central resonance width of the transition R e (0) belonging to the combination band ν 1 + ν 2 + (2ν 4 + ν 5 ) 1 of C 2 H 2 located at 7230.61313 cm 1 (polyad ∆P = 11): Lorentzian component (HWHM) of the dispersion Faddeeva profile. Simulation conditions: dipole moment: 0.515 mD, number density: 2.5 × 10 14 cm -3 , line intensity: 1.901 × 10 -24 cm/mol., mean transit-time rate: 220 kHz, directional beam power (linear polarization): 10 mW, beam waist size: 0.46 mm, frequency modulation: 350 MHz, modulation index: 0.3, mean Rabi frequency: 27 kHz.

  least, one numerical integral over the Doppler shift (exhibiting a Maxwell-Boltzmann distribution) is required. But, to be able of reproducing the EMF radial extension, we have three additional types of integration models. The first one consists of simply weighting the usual standard Dirac EMF (considered in the frequency domain) by a Gaussian radial extension. The two other approaches, which are mostly discussed in the paper, consider a more accurate description involving the beam impact parameter, and the molecular transit-time rate. This requires a single additional integration associated with a mean transit-time rate for the less CPU demanding model, or two additional integrations for the most sophisticated model which includes the distribution of transit-time rate associated with the usual Maxwell-Boltzmann velocity distribution relative to the standard gases.Simulations have been performed to mainly address the behavior of the FM signal in dispersion (or optical phase shift) of two transitions R e (0) exhibiting different line intensity and belonging to two different combination bands of the polyad 11 of acetylene C 2 H 2 . Realistic experimental parameters have been considered. Collision rates reflecting the gas pressure, has been used as the main variable parameter. However, the dependence to the mean transit-time rate has also been considered.

Table 1 :

 1 Coefficients X κ (ξ) associated with the FM detection on the first harmonics

	Sidebands Harmonic	in-phase	-κ	in-quadrature	0 in-phase in-quadrature	in-phase	κ	in-quadrature
	1

Table 2 :

 2 Fitted parameters extracted from the simulation of the transition R

	+ ν 1 4 of C 2 H 2 associated with Fig. 2.
	+ ν 3
	1

e (0) of the band ν

  This different behaviors prevent us of using the model involving a single mean transit-time rate cautionless. The intermediate collision regime is approximately bounded by the mean Rabi frequency and by the mean transittime rate. It is in this intermediate regime that the narrowest resonances are observed, at least when the mean Rabi

Available on request.

e -z 2 z 0 e t 2 dt

ACKNOWLEDGMENTS

The LPCA is "laboratoire conventionné au CNRS". P. Dupré is affiliated with CNRS.

the collision rate. Focusing on the resonance width components, we already notice that in the collision-free regime the width of the Lamb-dip can reach the Rabi frequency of the transition (under a power of 200 mW). The overall behavior is not very different from that observed for the Lamb-dip (DC) analysis [START_REF] Dupré | Dipole saturated absorption modeling in gas phase: Dealing with a Gaussian beam[END_REF]. However, in the collision-free regime, there is a clear divergence of the width of the Gaussian and of the Lorentzian components. In the simulations involving a double integration (in red on the figures), the Lorentzian component approaches the mean Rabi frequency while for those involving a triple integration (in green), this component approaches the mean transit-time rate. This is the result of the integration over the transit-time. The relative values of the Gaussian component to the Lorentzian component is different in the two approaches: in the simulation type #3 (in red), the Gaussian component dominates the profile width, while in the simulation type #4 (in green), both components are close (Fig. 4).

A C C E P T E D M

A N U S C R I P T Anyway, in the collisional regime, all the resonance widths are proportional to the collision rate. For the integration type #3, the Gaussian component is almost equal to the mean transit-time rate divided by 2 which is the expected value [START_REF] Dupré | Dipole saturated absorption modeling in gas phase: Dealing with a Gaussian beam[END_REF]. For the model type #4, this behavior is blurred by the large value of the Lorentzian component. Indeed, in this regime, the model type #4 provides a resonance profile which is never purely "Gaussian" (at the opposite of model type #3), but, in fact, the Lorentzian component remains larger than the Gaussian component. This can be seen on Fig. 6 where a plot obtained by fitting the resonances with a pure Lorentzian profile is also shown.

The passage between the two collisional regimes is clearly observable, it appears for a mean transit-time rate equal to the mean Rabi frequency.

Transition R e (0),

Collisional Effects

. In this section, we wish to compare the collisional effects under slightly different conditions, i.e., when the mean Rabi frequency is 4 times weaker than it is in the previous section (3.1), by keeping the same mean transit-time rate. To deal with lower Rabi frequencies, we preferred to consider a weaker transition. However, the nature of this transition makes that the simulations can be equally permuted on purpose. Under these new conditions, the ratio between the mean transit-time rate γ tr , and the mean Rabi frequency Ω Rabi [as defined by Eq. 56] reaches almost one order of magnitude.

Again, both types of modeling (i.e., involving two or three integrations) are considered. of the dispersion FM Faddeeva function, we can also observe a permutation of the Lorentzian component width versus the Gaussian component width for both models in the collision-free regime. This dispersion of the behaviors contrasts with the clear convergence of the behavior observed in the collisional regime as previous discussed. The double integration model (#3) exhibits a Lorentzian component which reaches the mean Rabi frequency (it can go even below around γ ul = 10 kHz). At very low collision rates, all the fitted width components exhibit, more or less the same behavior than that previously observed: the convergence towards the mean transit-time rate might be artificial, because of the limited range of collision rate explored. However, it is a strong saturation regime where the model may fail.

The width of each component of a dispersion Faddeeva profile be difficult to extract from experimental data potentially suffering from a poor signal-to-noise ratio. Hence, we plotted the fitted HWHM of the FM dispersion resonance by using a single dispersion Lorentzian function, see Fig. 9. This approach offers the advantage of providing data less dispersed than those obtained with a profile involving a dual width. Thus, we observe that the narrowing effects of the resonances are less pronounced. Typically, both models (involving a double and a triple integration) show a very similar behavior: the minimum resonance width is obtained for a collision rate close to the mean Rabi frequency, while the value of the HWHM is only a factor two less than the mean transit-time rate, i.e., it does not reach the mean Rabi frequency like previously shown. At very low pressure, the asymptotic regime is not accurately determined, but the high saturation degree which can be reached generates a power broadening, bringing the resonance width above the mean transit-time rate. In summary, when a single width is used to fit the dispersion resonances, the net effect is an attenuation of the difference between the two models, despite that the least-square sum (χ 2 ) can significantly increase up to two orders of magnitude between both models.

Finally, we have compared the equivalent integral of the resonances versus the collision rate. Similar plots to those shown on Fig. 5 are observed. Hence, they are not reported here. Both R e (0) transitions under study exhibit a central FM resonance whose position of the maximum of the integral is close to the mean transit-time rate. Actually, this maximum is slightly shifted towards the low collision rates for the simulations involving the lower mean Rabi frequency.

Power Broadening. In purpose of evaluating the power broadening, a set of simulations has been performed for the transition considered under low pressure conditions (collision rate: 60 kHz). The basic idea is to compare, or to predict experimental behaviors since the ultra narrowing of the resonances has been previously reported under such pressure conditions [START_REF] Baklanov | Narrow resonances of two-photon absorption of super-narrow pulses in a gas[END_REF][START_REF] Chebotayev | Supernarrow saturated absorption resonances[END_REF][START_REF] Ch | Ultrahigh-resolution saturation spectroscopy using slow molecules in an external cell[END_REF]. To analyze the simulations, the same procedure as that previously described has been used for two types of dispersion profile: a normalized CLF and a Faddeeva profile (see Fig. 10). If it is easy to extract the width of the two components of a Faddeeva profile from simulated data, this extraction may be more at risk when fitting experimental data exhibiting limited signal-to-noise ratio. As anticipated, the CLF HWHM is larger than the width of the Lorentzian component fitted from the Faddeeva profile (round symbols versus square symbols). Of course, both values converge towards the same one as soon as the Gaussian component becomes negligible. Actually, at high power, the same behavior is observed for both models: the Lorentzian component has a dominant contribution over the Gaussian component, while at low power, the "incomplete" model (#3) exhibits a Gaussian component dominating the Lorentzian contribution (under the current simulation conditions). On the other hand, at the same low power, the "full" Assuming that E (ω, r, z) is defined by Eq. 23

Taking into account the properties of the Bessel functions J κ (ξ) = (-1) κ J -κ (ξ), κ J 2 κ (ξ) = 1, and ignoring the terms at ω ≥ 2ω rf , it becomes

whose amplitude is centered at ω = 0.

Calculation of

By using Eqs. [START_REF] Shimoda | Line Broadening and Narrowing Effetcs[END_REF] and Eq. ( 17), by frequency shifting the Lorentz functions, and by ignoring the non-resonant terms (and setting ω κ = ω 0 + κω rf ), we obtain the first product in the RWA approximation (γ ul ω κ ) by temporally ignoring the r dependence

where δ ul κ is the detuning [see Eq. ( 28)], and where k κ = k 0 + κk rf , and with

By considering the expansion based on the Bessel coefficients

Finally, the terms of the convoluted product can be calculated by distinguishing (the terms around 2ω 0 , etc.. are ignored), i) for the diagonal terms (∆κ = 0)