Oleg Egorov
email: oleg.egorov@iao.ru

Michaël Rey

Andrei V Nikitin

Dominika Viglaska

New theoretical infrared line list for the methyl radical with accurate vibrational band origins from high-level ab initio calculations

Keywords: CH3, ab initio, potential energy surface, dipole moment surface, band origins, line list, large-amplitude motion

In the present work, high-level ab initio calculations were carried out for the ground electronic state of the methyl radical (CH3). The Dunning's augmented correlation-consistent orbital basis sets were employed up to the quintuple-ζ valence quality with the core-valence electron correlation [augcc-pCV5Z] combined with the single-and double-excitation unrestricted coupled cluster approach with a perturbative treatment of triple excitations [RHF-UCCSD(T)]. The explicitly-correlated version of the coupled cluster approach [RHF-UCCSD(T)-F12x{x=a, b}] was additionally applied with the core-valence cc-pCVQZ-F12 basis set in order to study convergence with respect to the basis set size. The contributions beyond the coupled cluster level of the theory like Douglas-Kroll-Hess scalar relativistic Hamiltonian, diabatic Born-Oppenheimer corrections and high-order electronic correlations have been included to the ab initio potential energy surfaces (PESs). It is shown that the theoretical band origins of CH3 converge gradually to the experimental values when applying the ab initio PESs using the aug-cc-pCVXZ [X = T, Q, and 5] basis sets. For the first time, all available experimental band origins of the gaseous CH3 are reproduced within an accuracy of 0.2 cm -1 using a newly developed PES extrapolated to the Complete Basis Set limit [CBS(TQ5Z)]. The reached accuracy is one order of magnitude better than that of the best available calculations. A new theoretical infrared line list was generated for astrophysical applications using an ab initio dipole moment surface computed at the RHF-UCCSD(T)/aug-cc-pCVQZ level of the theory. The manifestation of a large amplitude motion in CH3 is also discussed.

INTRODUCTION

Methyl radical (CH3) is an intermediate species in chemical reactions associated with the production and destruction of the complex hydrocarbons, precursor in the formation of the carbon films during the chemical vapor deposition (see e.g. refs 1-4). Due to its very high reactivity, the CH3 lifetime cannot be long under the Earth's conditions but it increases in rarefied gaseous media which is typical for astronomical objects. There are essentially two natural ways to form CH3. The first one is based on the photolysis of methane, namely CH4 + hv → CH3 + H that can occur in planetary atmospheres. [START_REF] Mordaunt | Primary product channels in the photodissociation of methane at 121.6 nm[END_REF][START_REF] Gans | Photolysis of methane revisited at 121.6 nm and at 118.2 nm: quantum yields of the primary products, measured by mass spectrometry[END_REF][START_REF] Romanzin | Combined experimental and theoretical studies on methane photolysis at 121.6 and 248 nm-implications on a program of laboratory simulations of Titan's atmosphere[END_REF] The second one is based on a series of ion-molecular reactions in the cold low-density interstellar medium CHx + + H2 → CHx+1 + + H, where the reaction С + H3 + → CH + + H2 initiates the formation of the C-H bonds. [START_REF] Herbst | The formation and depletion of molecules in dense interstellar clouds[END_REF] Note that ethane, which is the second strongest abundant hydrocarbon in the atmospheres after CH4, can be formed via the reaction: 2CH3→C2H6. So far, CH3 was detected by the Infrared Space Observatory in Saturn, [START_REF] Bézard | Detection of methyl radicals (CH3) on Saturn[END_REF] Neptune, [START_REF] Bézard | Detection of the Methyl Radical on Neptune[END_REF] and in the interstellar medium toward the Sagittarius A* radio source. [START_REF] Feuchtgruber | Detection of Interstellar CH3[END_REF] The present work aims at first computing an accurate PES for CH3 using advanced ab initio calculations before constructing a comprehensive, room-temperature line list from variational calculations. The previously published works (see e.g. refs 12-15) were based on both the multireference configuration interaction method and the explicitly correlated coupled cluster approach but they were not able to accurately reproduce the observed energy levels. For example, in the last study by Adam et al. [START_REF] Adam | Variationally Computed IR Line List for the Methyl Radical CH3[END_REF] the errors on the band centers were of 4 and 9 cm -1 for v2 and 2v2, respectively. Such deviations are surprisingly big if we keep in mind that the ground electronic state of CH3 is quite well isolated from the first excited state located at 216 nm or ~ 46000 cm -1 . This paper is a continuation of our previous studies focused on the construction of theoretical line lists of rigid and non-rigid polyatomic molecules of different symmetries from accurate ab initio potential energy (PESs) and dipole moment surfaces (DMSs). It is structured as follows. First, the state-of-the-art in the CH3 spectroscopy is given in Section 2. The construction of the PES and DMS will be reported in Section 3. The final results including the corrections due to advanced ab initio methods and high order electronic correlation effects will be discussed in Section 4 before concluding in Section 5.

The ab initio calculations reported in this work were performed by using the MOLPRO, [START_REF] Werner | The Molpro quantum chemistry package[END_REF][START_REF] Werner | version 2019.2, a package of ab initio programs[END_REF] MRCC, [START_REF] Kállay | The MRCC program system: Accurate quantum chemistry from water to proteins[END_REF][START_REF] Kállay | a quantum chemical program suite[END_REF] and CFOUR [START_REF] Matthews | Coupled-Cluster Techniques for Computational Chemistry: the CFOUR Program Package[END_REF][START_REF] Stanton | a quantum chemical program package[END_REF] packages while the nuclear motion calculations were carried out from our home-made variational computer code TENSOR [START_REF] Rey | Ab initio ro-vibrational Hamiltonian in irreducible tensor formalism: a method for computing energy levels from potential energy surfaces for symmetric-top molecules[END_REF][START_REF] Rey | Complete nuclear motion Hamiltonian in the irreducible normal mode tensor operator formalism for the methane molecule[END_REF][START_REF] Rey | Group-theoretical formulation of an Eckartframe kinetic energy operator in curvilinear coordinates for polyatomic molecules[END_REF] applied with success in refs 25-34 for computing accurate spectroscopic line lists of semirigid molecules. As in ref 35, a hybrid nonrigid Hamiltonian model was built in this work using the Hougen-Bunker-Johns formalism [START_REF] Hougen | The vibration-rotation problem in triatomic molecules allowing for a large-amplitude bending vibration[END_REF] where the relation between the curvilinear  and effective coordinates  was explicitly established in ref 37.

STATE-OF-THE-ART IN THE CH3 SPECTROSCOPY

CH3 has one unpaired electron, therefore its ground electronic state (denoted to as 𝑋 ̃2 2 A in the D3h point group) is a doublet state. Consequently, the vibration-rotation spectra should exhibit spin-rotation splittings. Due to its planar configuration at the equilibrium geometry in the 𝑋 ̃2 2 A state, CH3 has no permanent dipole moment and consequently exhibits only weak magnetic-dipole rotational transitions (see ref 38 and reference therein for further details).

Early experimental studies by Herzberg and Shoosmith [START_REF] Herzberg | Absorption spectrum of free CH3 and CD3 radicals[END_REF][START_REF] Herzberg | The Spectra and Structures of Free Methyl and Free Methylene[END_REF] were devoted to the ultraviolet absorption spectrum where the band at 216 nm (and 214.5 nm for CD3) was observed. This band obeys Rydberg formula and corresponds to the excitation of the unpaired electron located on the molecular 2pz character orbital into the 3s one. The structure and dissociation dynamics of the vibronic levels of the 3s Rydberg state (𝐵 ̃2 1 A ) which is the lowest excited electronic state of CH3 (D3h) were further studied in refs 41-45. The vibrational feature at 212.7 nm observed by Callear et al. 41 was assigned to the hot band 𝐵 ̃2 1 A (0100)← 𝑋 ̃2 2 A (0100) by Westre et al. 43 This hot band was detected at 211.6 nm for CD3. Another absorption at 207 nm was assigned to the cold band 𝐵 ̃2 1 A (1000)← 𝑋 ̃2 2 A (0000) from analysis of the resonance Raman spectrum. 43 At higher energies, in particular at 193 nm, a broad continuum was observed due to the photodissociation of CH3. [START_REF] North | Photodissociation dynamics of the methyl radical 3s Rydberg state[END_REF] Note that the energies of the Rydberg states correlate with the excited electronic states of methylene (CH2) in the dissociation channel CH2(C2v) + H( 2 S) which leads to the strong predissociation properties. 47- 49 In particular, for the 3s Rydberg state the empirical and theoretical values of the barrier estimated in refs 42 and 50 were of 6.29 kcal/mol (~2200 cm -1 ) and 2417 cm -1 , respectively for the reaction:

CH3(𝐵 ̃2 1 A ) → CH2(𝑎 ̃1 1 A ) + H( 2 S
), where 𝑎 ̃1 1 A is the first excited singlet electronic state of CH2. Next electronic transitions in CH3 are associated with the 3pxy ( 2E ) and 3pz ( 2 2 A ) Rydberg states. The energy of the 3pxy( 2E )←𝑋 ̃2 2 A vertical excitation calculated in refs 48, 51 and 52 varied within 6.95-7.17 eV (172.9-178.4 nm). The 3pz( 2 2 A )←𝑋 ̃2 2 A band was measured by the Resonance- Enhanced Multiphoton Ionization (REMPI) spectroscopy in refs 53-61. Simultaneous absorption of two photons with energies between 325 and 336 nm resulted in the detection of several absorptions with the maximum located at 333.4 nm that is about 166.7 nm (59972 cm -1 ) in terms of two-photon energy. The frequency of the v2 mode at the 3pz( 2 2 A ) Rydberg state extracted in refs 55, 57 and 61 was about two times bigger when comparing with the ground electronic state whereas other fundamentals, particularly v1 57,61 and v3 58 were red-shifted by ~74 cm -1 . Some of the higher excited Rydberg states (4p, 4f, and 5f) were also detected by the REMPI spectroscopy in refs 55, 62 and 63. As to the infrared spectra, the region of the out-of-plane bending vibration (v2 band) associated with the strongest z component of the molecular dipole of CH3 and located at 606 cm -1 (or 16.50 μm) was one of the first for which high-resolution studies were carried out. In early works, [START_REF] Milligan | Infrared and Ultraviolet Spectroscopic Study of the Products of the Vacuum-Ultraviolet Photolysis of Methane in Ar and N2 Matrices. The Infrared Spectrum of the Free Radical CH3[END_REF][START_REF] Snelson | Infrared matrix isolation spectrum of the methyl radical produced by pyrolysis of methyl iodide and dimethyl mercury[END_REF][START_REF] Jacox | Matrix isolation study of the infrared spectrum and structure of the CH3 free radical[END_REF] the spectra of CH3 were recorded by the matrix isolation technique with a certain shift in the frequency depending on the inert gas matrix applied (argon, nitrogen, neon etc.). Tan et al. [START_REF] Tan | Infrared Spectrum of Gaseous Methyl Radical by Rapid Scan Spectroscopy[END_REF] recorded the spectra of the gaseous short-lived CH3 with the help of the rapid-scanning infrared spectrometer. Yamada et al. [START_REF] Yamada | Diode laser study of the ν2 band of the methyl radical[END_REF] applied the pyrolysis and electric discharge approaches to produce CH3 and used a diode laser spectrometer. As a result, the rotational structures of the v2, 2v2-v2 and 3v2-2v2 bands were resolved in the gas phase. The spin-rotation parameters extracted from the observed spin-rotational splitting were in a good agreement with the theoretical estimation.

Although the v1 symmetric stretch band of CH3 is not infrared active, its rotational structure (especially, Q-branch) was explicitly resolved at 3004 cm -1 by the Coherent anti-Stokes Raman Spectroscopy (CARS) in refs 69-71 using the second harmonic of the YAG laser as a pumping beam and by the resonance Raman spectroscopy in ref 72. The spectra obtained in ref 71 under the nearnascent conditions show the domination of the ground vibrational state of CH3 when it is produced by the UV laser photolysis from CH3I. This agrees with the theoretical study [START_REF] Amatatsu | Ab initio potential energy surfaces and trajectory studies of A-band photodissociation dynamics: CH3I*→CH3+I and CH3+I*[END_REF] which predicts a slow relaxation from the pyramidal to the planar configurations after the bond breaking.

The x component of the molecular dipole is responsible for the infrared absorption in the regions of the v3 and v4 bands (in-plane degenerate asymmetric stretch and bend, respectively). The difference frequency laser spectrometer applied in ref 74 allowed to record several series of the Qbranch transitions of the v3 band located at 3161 cm -1 . The absorption of this band was then remeasured in ref 75 under the sub-Doppler conditions that allowed to resolve the spin-rotation and partially nuclear hyperfine splitting. To our knowledge, there is no studies of the v4 band of the gaseous CH3 to confirm results of the matrix isolation spectroscopy: its band origin varies from 1385 to 1403 cm -1 , depending on the gas matrix (see e.g. refs 65, 76-78).

When extracting the line intensities of CH3 from the measured spectra, the most complicated issue is the determination of the concentration and temperature inside the cell under the plasma conditions. Currently, the vibrational transition moments ( V V    ) of three infrared bands, namely, v2, 2v2-v2, and v3 were analyzed at different levels of accuracy. In early works, [START_REF] Yamada | The transition dipole moment of the ν2 band of the methyl radical[END_REF][START_REF] Wormhoudt | A measurement of the strength of the ν2 band of CH3[END_REF] the peak absorbances of the single transition belonging to the Q-branch of v2 was measured to calculate its V  , whereas for the hot band 2v2-v2 the V V    value is about 1.5 times stronger. [START_REF] Stancu | Measurement of the transition dipole moment of the first hot band of the v2 mode of the methyl radical by diode laser spectroscopy[END_REF] Note that the V V    values in experimental works were determined at the zero-order approximation i.e. without rotation-vibration couplings. Such approach can provide satisfactory results for the strongest bands (e.g. for v2 of CH3) but turns out not sufficient to predict weaker bands whose intensities mainly come from vibrational-rotational interactions. A new, non-perturbative approach was recently proposed in ref 87 to obtain the transition moment spectroscopic effective parameters directly from the ab initio DMS for computing both cold and hot bands.

CONSTRUCTION OF THE PES AND DMS

3.1. Analytical representation. The PES and DMS of CH3 are both expressed in this work as a Taylor power series expansion in terms of symmetry coordinates i S . This analytical representation, which converges quite rapidly despite the out-of-plane bending vibration, has been recently used to describe the NH3 PES with its 2 minima separated by a barrier of 2000 cm -1 . [START_REF] Egorov | New Ab Initio Potential Energy Surfaces for NH3 Constructed from Explicitly Correlated Coupled-Cluster Methods[END_REF] The D3h symmetry coordinates associated with the symmetric and asymmetric stretching and in-plane bending vibrations are defined using the projection technique [START_REF] Bunker | Molecular Symmetry and Spectroscopy[END_REF] as follows

         
                            (1)
where k y and k  with k = 1, 2, and 3 are functions depending on the internal coordinates. Here, we have considered the Morse-type functions

  1 exp k k e y b r r        , (2)
for the stretching part with b = 1.9 Å -1 and the cosine functions

    cos cos k k e      , (3)
for the bending part. Note that Eqs. (1-3) are the same for both the PES and DMS while the out-of-plane bending coordinate is different. Indeed, the sign of the z component of the dipole between 0 о to 90 о and between 180 о to 90 о is different whereas it is the same for the PES (see Figure 1). Consequently, a sine function

    1 2 sin sin A e S      (4)
will be used for the PES (here ee    , see ref 37) and a cosine function

  2 2 cos A S   (5)
for the DMS, where  can be related to the three in-plane bending angles as follows  

1 2 3 2 arcsin sin / 6 3             . ( 6
)
The use of symmetry coordinates which transform as the irreducible representations of the molecular symmetry or point group [START_REF] Bunker | Molecular Symmetry and Spectroscopy[END_REF][START_REF] Papoušek | Molecular Vibrational-Rotational Spectra. Theory and Applications of High Resolution Infrared[END_REF] allows to express the PES or DMS as sum-of-products of irreducible tensor operators

p j R   max 11 1 2 3 3 4 4 0
, , , , ,

p AAE E E E p a b a b j j pj V S S S S S S C        R (7)
where j C are expansion coefficients to be determined and the irreducible tensors are written as

123 4 5 6 11 1 2 3 3 4 4 1 ii i i i i n pp p p p p AA p E E E E j i a b a b i T S S S S S S                                       R . ( 8
)
The i T 's are computed from the D3h Clebsch-Gordan coupling coefficients and each term in Eq. ( 8)

is of order p. The tensor formalism used in this work shares the same coupling coefficients and phase factor conventions with the variational computer code TENSOR. Since the PES is totally symmetric, Eq. ( 7) must contain only the terms which transform as 1 A . The dipole moment operator has a similar expansion, except its z and (x,y) components will transform as 2 A and E , respectively.

Ab initio calculations.

In order to study the ground electronic state of the open-shell CH3 radical, the spin-unrestricted coupled cluster (CC) method [RHF-UCCSD(T)] implemented in MOLPRO by Knowles et al. 90,91 was employed. The unpaired electron was located on the 2pz molecular orbital during all calculations that corresponds to the 𝑋 ̃2 2 A state of CH3. The pragmatic CC approach provides accurate and size-consistent calculations of the correlation energy especially for single-reference electronic states. To estimate the multireference character of the 𝑋 ̃2 2 A state, the T1 and D1 diagnostic numbers were extracted from the output files of MOLPRO. As depicted in Contrary to previous ab initio studies, in particular refs 12 and 13, the Dunning's augmented correlation-consistent orbital basis sets [START_REF] Dunning | Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen[END_REF] were applied up to the quintuple-ζ valence quality with the core-valence electron correlation [aug-cc-pCV5Z, hereafter, ACV5Z]. In the case of the correlationconsistent basis set, the ab initio energies can be further extrapolated to the Complete Basis Set (CBS) limit using different extrapolation formula. [START_REF] Peterson | Benchmark calculations with correlated molecular wave functions. IV. The classical barrier height of the H+H2→H2+H reaction[END_REF] In this work, we used the three-point exponential functional  

exp total total X CBS E E A C X      (9)
which was recently applied with success in ref 35 for NH3. In Eq. ( 9) total X E and total CBS E are the total energies (Hartree-Fock (reference) + correlation) corresponding to the ACVXZ (X = T, Q, and 5) basis set and CBS limit, respectively.

Along with the "standard" CC method, the simplified explicitly correlated versions [RHF-UCCSD(T)-F12x{x=a, b}] of Knizia et al. [START_REF] Knizia | Simplified CCSD(T)-F12 methods: Theory and benchmarks[END_REF] were also employed in order to reach the CBS limit without extrapolation. According to refs 35 and 97, the accuracy of the F12b method combined with the correlation-consistent core-valence cc-pCVQZ-F12 basis set (hereafter, CVQZ-F12) is comparable to that of the sextuple-ζ basis set.

In the CBS limit, further corrections beyond the CC level of the theory are required to predict the band centers within 1 cm -1 or below of accuracy. There are three types of corrections that are actually analyzed when considering the ground electronic state: scalar relativistic effects (relativistic), diabatic Born-Oppenheimer correction (DBOC), and contributions to the correlation energy from the highly-excited Slater determinants. Since the convergence of the high-order CC methods is quite fast with respect to the basis set size, the cc-pVTZ basis set will be enough to completely include the contributions from the triple [CCSDT] and quadruple [CCSDT(Q)] excitations whereas the cc-pVDZ basis set will be used for the pentuple excitations. Note that the MRCC and CFOUR packages consider the unrestricted Hartree-Fock (UHF) reference wave function, making the open-shell CC method different to that implemented in MOLPRO. For this reason, the ab initio energies computed from MOLPRO were not used when making the scaling from CCSD(T) to CCSDT. All the energy differences, namely [CCSDT-CCSD(T)]/cc-pVTZ, [CCSDT(Q)-CCSDT]/cc-pVTZ, and

[CCSDT(Q)-CCSDTQP]/cc-pVDZ were calculated using the ab initio energies from MRCC.

The relativistic correction was calculated using the Douglas-Kroll-Hess (DKH) scalar relativistic Hamiltonian available in MOLPRO. A fourth order DKH Hamiltonian was combined with the ACVTZ-DK basis set. The energy contributions from higher order transformations were found to be less than 10 -4 cm -1 . The DBOC correction was calculated at the UHF-CCSD/cc-pVTZ level in CFOUR. According to extensive benchmarks from ref 98, "…the CCSD approach is a good choice for the accurate calculation of the DBOC…" because it provides the contribution of about 90 % to the electron correlation energy.

As to the ab initio dipole moment, it was calculated via the finite difference approximation as the first derivative of the potential energy with respect to the electric field. Since CH3 is a symmetric top molecule, only two components of the dipole moment were necessary. Here, a field was applied along the x and z axes with the magnitudes of 0.002 a.u. and -0.002 a.u. The RHF-UCCSD(T)/ACVQZ method was applied for the energy calculations. S  of the stretching and bending coordinates take values between 0 to 0.20 and between 0 to 0.72, respectively, while the ranges for

1 1 A S  , 3 E a S  and 4 E a
S  are the following: [-0.2, 0.25], [-0.2, 0.25], and [-0.72, 0.90], respectively. Finally, the Grid_1, which is composed of 36032 nuclear configurations up to 15000 cm -1 , was used to compute (i) the PES at the RHF-UCCSD(T)/ACVQZ and RHF-UCCSD(T)-F12x{x=a,b}/CVQZ-F12 levels and (ii) the x and z components of the DMS at the RHF-UCCSD(T)/ACVQZ level. This

Grid_1 was fitted with a Root Mean Square (RMS) deviation of 0.18 cm -1 for the PES and of 0.000627 (resp. 0.00038) Debye for the z (resp. x) component of the DMS using 303 and 61 (resp. 143) expansion coefficients in Eq. (7). The fitted sets of the parameters included the irreducible tensor operators up to 8 th and 6 th orders for the PES and DMS, respectively.

The second grid (called Grid_2) was generated using totally symmetric ( 1 A ) irreducible tensors up to order 8 and is composed of 5849 points for the PES up to 10000 cm -1 . The Grid_2 was used to calculate the energies at the RHF-UCCSD(T)/ACVXZ (X=T, Q, and 5) levels as well as all corrections that are quite demanding (see section 3. a Details of the calculation can be found in section 3.2.

All geometries were obtained from the fitted PESs.

The biggest contribution to re of 0.000148 Ang. was obtained using the CCSDT correction.

The contributions from CCSDT(Q) and CCSDTQP of 0.000063 Ang. and 0.000007 Ang. clearly indicate the negligible effect of the higher excited determinants (e.g. sextuple order) to the correlation energy of CH3. We also tested the convergence of the highly-correlated CC methods with respect to the basis set size. The (T)→T scaling using the lower level cc-pVDZ basis set resulted in a difference of +0.000055 Ang. in re compared to the cc-pVTZ one. On the other hand, the scaling from T to (Q)

gave the same contribution to re when using the cc-pVDZ and cc-pVTZ basis sets. This simple test confirms the fast convergence of the high-order CC methods with respect to the basis set size and gives the guaranty that the cc-pVDZ basis set is more than enough for the (Q)→P scaling. Currently, the CCSDTQP/cc-pVTZ calculations are too costly. [START_REF] Amano | Difference frequency laser spectroscopy of the ν3 band of the CH3 radical[END_REF] and Davis et al. [START_REF] Davis | Jet-cooled molecular radicals in slit supersonic discharges: Sub-Doppler infrared studies of methyl radical[END_REF] For variational calculations the hybrid model implemented in our homemade computer code TENSOR was applied. The kinetic energy and potential parts were Taylor expanded in the five normal mode coordinates describing the `small' vibrations (associated here with the ν1, ν3, and ν4 modes) for each point of a numerical grid of the nonrigid coordinate  describing the large amplitude motion (associated here with the ν2 mode). To this end, the correspondence between the coordinate  involved in the kinetic energy operator of the Hougen-Bunker-Johns Hamiltonian [START_REF] Hougen | The vibration-rotation problem in triatomic molecules allowing for a large-amplitude bending vibration[END_REF] known that the F12a method combined with the CVQZ-F12 basis set provides an overestimated values of the correlation energy whereas the results based on the F12b method are comparable with the sextuple-ζ basis set. This is corroborated by our calculation: the frequency of the v2 band is 0.5 lower than the experimental value when using the F12b method and it is about 1 cm -1 bigger when the F12a method is applied. Therefore, the values of the converged vibrational energies with respect to the basis set size must fall between those provided by the F12a and F12b methods. In others words, they fall in the region of the sextuple-ζ basis set or slightly lower if we remember the approximate nature of the F12x{x=a, b} methods.

Although the ACV6Z basis set was not considered here because of its cost of calculation, the results based on the extrapolated to the complete basis set potential energies [RHF-UCCSD(T)/CBS(TQ5Z)] were enough to reproduce for the very first time all available experimental band origins of CH3 with an RMS error of 0.177 cm -1 . Such deviation is close to that obtained from the weighted fit of the PES (0.183 cm -1 ). Note that the accuracy of the faster method F12b/CVQZ-F12 is of 1 cm -1 which can be thus seen as a good alternative to the standard CC approaches for CH3. A state of CH3 using the default settings in MOLPRO and the CVQZ-F12 basis set. The following re values (in Ang.) were obtained: 1.075829 and 1.075823 using RCCSD(T)-F12x{x=a,b}, and 1.075878 and 1.075874 using RHF-UCCSD(T)-F12x{x=a,b}. As it can be seen, the difference of 0.00005 Ang. is small and can be explained by the default energy threshold of 10 -6 a.u. during the geometry optimizations in MOLPRO. Thus, we assume two main reasons why the predictions given in ref 15 are less accurate than expected. 1) the equilibrium geometry re=1.077369 Ang. used in ref 15 for the construction of the PES is not accurate enough and corresponds to a value obtained between ACVTZ and ACVQZ (see Table 1). 2) the corrections due to relativistic, DBOC, and high-order electronic correlation effects were not considered though they do not contribute significantly for the v2 and 2v2 bands (see Section 4.2).

Corrections beyond RHF-UCCSD(T).

In this section, the contributions to the band origins of CH3 obtained from various corrections are analyzed and discussed. Note that in Tables 2 and 3 all corrections have been already included to our PESs. It is also worth mentioning that according to Table 4 the various contributions on the CH3 vibrational levels are closer to those of the CH4 molecule than to those of NH3 belonging to the same symmetry. Indeed, in NH3, [START_REF] Egorov | New Ab Initio Potential Energy Surfaces for NH3 Constructed from Explicitly Correlated Coupled-Cluster Methods[END_REF] the relativistic and DBOC corrections had big contributions (up to 9.3 cm -1 and -4.8 cm -1 , respectively) especially for the vibrational levels where v2 is involved. In CH4, [START_REF] Nikitin | First fully ab initio potential energy surface of methane with a spectroscopic accuracy[END_REF] on the contrary, their effects were within ±1 cm -1 for most of the vibrational bands whereas the high-order electronic correlations gave the main contribution. a First 50 levels are displayed. Total list up to 7000 cm -1 is given in Supporting Information I. Given values are the differences in the band origins (in cm -1 ) when the corresponding correction is taken into account.

For example, the "DBOC" column means the difference in the band origins calculated at the [RHF-UCCSD(T)/CBS(TQ5Z)+Relativistic+DBOC] and [RHF-UCCSD(T)/CBS(TQ5Z)+Relativistic] levels of the theory

In the case of CH3, the relativistic corrections give rather small contributions both for the stretching and out-of-plane bending vibrations and contribute more significantly on the bands where v4 is involved. Conversely, the contributions from DBOC is quite small for v4. In terms of absolute values, the relativistic and DBOC corrections were of 0.38 cm -1 and 0.45 cm -1 in average. As expected, the inclusion of the CCSDT correction resulted in the significant shifts (-2.73 cm -1 in average) for all band centers, in particular for the stretching vibrations. The last two corrections, namely CCSDT(Q) and CCSDTQP, gave contributions which gradually decrease: 1.27 cm -1 and 0.09 cm -1 in average, respectively. Finally, among the four fundamental vibrations, the sum of these corrections was the biggest for the v1 and v3 stretching bands (-4.221 cm -1 and -4.061 cm -1 , respectively) while we have a shift of -1.184 cm -1 for v4 and of 0.093 cm -1 for v2. In this work, the CH3 PES was constructed by following the same approach as for NH3 in ref 35. Both even and odd powers were allowed for the 2 S symmetry coordinate which transforms as the totally symmetric representation (see Eq. ( 4)). As expected, there is a strong linear term in the PES with a value of 0.10323 a.u associated with the S2 coordinate. For NH3 the linear terms were removed from the PES expansion after a geometry optimization. As an alternative, a fit using  

22

SA  has been performed, leading to an RMS error of 0.54 cm -1 against 0.18 cm -1 when using  

21

SA  . The deviations of the fit of the 1D cut along the  coordinate are displayed in Figure 3. After several trials, the PES had to be expanded up to order 10 when using  

22

SA  in order to obtain an RMS error similar to that when using   21 SA  . Moreover, high order terms in the PES generally require computing additional points in the grid to obtain well-defined fitted high-order parameters. All the results presented in Tables 2-4 have been obtained using the new hybrid version of the TENSOR dealing with a large-amplitude coordinate.

New theoretical CH3 line list in the infrared.

The new ab initio PES and DMS presented in the previous section were used to generate the variationally-computed infrared line list of CH3 that is compared below with the semi-empirical line list of Adam et al. 15 In this latter reference the ab initio DMS corresponds to the RCCSD(T)/AVTZ theoretical level of the theory. Here, "semiempirical" means that the ab initio band origins of certain bands (namely, v2, 2v2, v1, v4, and v3) were replaced by the experimental values while the remaining bands are displayed in Table 3 in column "III".

Since the ground electronic state of CH3 is a doublet there is a fine structure due to the spinrotation coupling. According to ref 68, the spin-rotation splitting has the magnitude of 0.01 cm -1 that is lower than the expected accuracy of our theoretical line list. This was also corroborated in Tab. VII S are y components of the rovibronic and electron spin angular momenta, respectively. [START_REF] Bunker | Molecular Symmetry and Spectroscopy[END_REF] In this work, the spin-rotation coupling term was not considered because of its quite weak contribution. Therefore, our spectral line positions are simply the transition frequencies between the vibration-rotation energy levels:

-1 cm fi f i v E E    .
The intensity of the spectral line [in cm -1 (10) where the subscripts f and i refer to the final and initial (lower) states;   Our theoretical line list has been computed using the TENSOR computer code in the spectral range [0-7000] cm -1 up to J = 20 ( JN  , S = 0). The integrated intensities of the strongest cold and hot bands are displayed in Table 5. The absorbance spectra simulated at P = 0.1 Torr and L = 20 cm with a Gaussian profile and at resolution of 0.1 cm -1 are given in Figure 4. According to and4c) where empirical corrections were applied. In the regions where no observation is available, the difference is much bigger. For example, the spectra in the v4( E ) and v3+v4( E ) bands are shifted by +5 cm -1 and -6 cm -1 (see Figures 4b and4d). It is also worth mentioning that our line intensities for v4( E ) are two times stronger compared to those in ref 15. Our calculation is corroborated by the measurements of Snelson et al. [START_REF] Snelson | Infrared matrix isolation spectrum of the methyl radical produced by pyrolysis of methyl iodide and dimethyl mercury[END_REF] who obtained an intensity ratio of 10 between v2( 2 A ) and v4( E ) (see Table 5).

    3 8 exp 1 exp 3 fi fi i f i i BB vv E I g S f i h c Q T k T k T                                

SUMMARY AND CONCLUSION

In this article, high-level ab initio calculations for CH3 have been performed in order to reproduce the vibrational band centers of the ground electronic state within an accuracy of 1 cm -1 .

The augmented correlation-consistent orbital basis sets were tested up to quintuple-ζ valence quality combined with the core-valence electron correlation [aug-cc-pCV5Z]. The quality of the calculated levels obtained from the "standard" open-shell coupled cluster approach [RHF-UCCSD(T)] was improved by including the following contributions: Douglas-Kroll-Hess scalar relativistic Hamiltonian, diabatic Born-Oppenheimer correction, and high-order electronic correlations. It has been thus demonstrated that the variationally computed band origins of CH3 converge gradually to the experimental ones when applying the ab initio PESs using the aug-cc-pCVXZ [X = T, Q, and 5] basis sets. Finally, an RMS error of 0.18 cm -1 with respect to all the available experimental band origins of the gaseous CH3 was obtained for the first time by extrapolating the ab initio potential energies to the Complete Basis Set limit [CBS(TQ5Z)]. Note that the reached accuracy is one order of magnitude better than that of the best available calculations. Another ab initio PES based on the less demanding explicitly correlated RHF-UCCSD(T)-F12b/cc-pCVQZ-F12 method allowed to predict the band origins with an RMS error of 0.66 cm -1 .

Among the fundamental vibrations, the corrections beyond the [RHF-UCCSD(T)] level of the theory had the biggest contributions for the stretching ones (up to -4.221 cm -1 ) while the in-plane and out-of-plane bending vibrations have been corrected by -1.184 and 0.093 cm -1 , respectively. This work also confirmed the presence of the large-amplitude behavior for the out-of-plane v2 mode of CH3. Since the  

21

SA  linear term cannot be neglected, the shape of the PES is not totally quadratic near the equilibrium configuration making very limited the description of v2 in terms of the "standard" normal mode approach.

Our newly-developed ab initio PES was applied to calculate all the band centers falling in the region [0-7000] cm -1 . In this work, we report the most accurate theoretical values for the fundamental v4 band at 1390.66 cm -1 as well as for all other combination cold and hot bands with strong integrated intensities. The main advantage of using "pure" ab initio calculations for the construction of the line list has been clearly demonstrated when compared with the previously published empirically corrected line list. [START_REF] Adam | Variationally Computed IR Line List for the Methyl Radical CH3[END_REF] Our theoretical line list is able to predict quite accurately spectral regions which were not studied so far. Undoubtedly, the results presented in this work will be useful for future highresolution experimental studies of CH3 as well as for astrophysical applications.

Our line list as well as the ab initio PESs and DMS are made available in Supporting Information.

Supporting Information

  et al.[START_REF] Stancu | Line strengths and transition dipole moment of the ν2 fundamental band of the methyl[END_REF] examined the line intensities of ten transitions of the v2 band by including line profiles to determine accurately the temperature. The obtained value of V V    was in a good agreement with the ab initio estimations, particularly, with ref 82. According to refs 83-85, the v3 band is about 6.5 times weaker than v2 in terms of V V  

Figure 1 .

 1 Figure 1. Behavior of the potential energy (a) and z component of the dipole (b) of CH3 along the out-ofplane bending angle (calculations were made at the RHF-UCCSD(T)/ACVQZ level; all non-scanned coordinates were fixed to their equilibrium values).

Figure 2 ,

 2 Figure 2, we can see that the maximum values of T1 and D1 are about two times less than their critical values of 0.02 and 0.05, respectively, established in refs 92 and 93 for all the nuclear configurations of the grid.

Figure 2 .

 2 Figure 2. T1 (a) and D1 (b) diagnostic numbers for the 𝑋 ̃2 2 A state of CH3 produced by MOLPRO when calculating the potential energy at the RHF-UCCSD(T)/ACVQZ level of the theory for "Grid_1" composed of 36032 nuclear configurations (for more details, see section 3.3). Note that the averaged T1 and D1 values are one third of their critical values of 0.02 92 and 0.05, 93 respectively.

  2). As in refs 35 and 99, the energy differences, namely, RHF-UCCSD(T)/[ACVTZ-ACVQZ], RHF-UCCSD(T)/[ACV5Z-ACVQZ], and RHF-UCCSD(T)/[CBS(TQ5Z)-ACVQZ] have a smooth behavior and were fitted with an RMS error less than 0.1 cm -1 . The predicted energy differences were then used to correct the energies RHF-UCCSD(T)/ACVQZ calculated for the whole Grid_1 in order to obtain extrapolated ACVTZ, ACV5Z, and CBS(TQ5Z) PESs. The same strategy was employed to include the sum of the five corrections (relativistic + DBOC + CCSDT + CCSDT(Q) + CCSDTQP) fitted with an RMS error of 0.045 cm -1 to the whole Grid_1.

4. 3 .

 3 Large amplitude motion in CH3. The choice of the irreducible representation for the 2 S symmetry coordinate has to be discussed. On the one hand, the CH3 PES is symmetric with respect to the equilibrium configuration (contrary to the NH3 molecule), therefore the out-of-plane bending vibration should obey the "standard" rules for the bending mode, namely allowed only even powers in the PES expansion and thus implying a quadratic shape near the equilibrium configuration.[START_REF] Schwenke | A theoretical study of the ro-vibrational spectrum of the X state of CH3[END_REF] This quadratic shape is crucial if, for example, the normal mode coordinates are used.On the other hand, an anomalously large anharmonicity was revealed for the v2 mode by Yamada et al.,68 proving that the shape of the PES is not totally quadratic but quite flat instead, even for vibrations with small amplitudes. This feature was further examined by Hirota et al., 101 who confirmed that higher order terms are needed in the PES for an accurate description of the v2 mode because of the slow convergence of the normal mode series. According to ref 101, the out-of-plane bend strongly interacts with the 1 A symmetric stretch. In particular, the negative sign of the F122 anharmonic potential constant (defined in terms of the curvilinear coordinates) can be explained by the elongation of the C-H bonds when the out-of-plane angle is excited.

Figure 3 .

 3 Figure 3. Fit of the ab initio 1D cut of the CH3 PES along the out-of-plane bending angle. Blue and red scatters correspond to the fits based on the 1 A and 2 A irreducible tensors operators up to the 8 th order. The ab initio energies are calculated at the RHF-UCCSD(T)/ACVQZ level of the theory; all non-scanned coordinates were fixed to their equilibrium values.

  of ref 102. In this latter reference, the spin-rotation Hamiltonian was fitted as

Figure 4 .

 4 Figure 4. Absorbance of CH3 simulated using the ab initio line list developed in this work compared to that based on the semi-emprirical line list of Adam et al. 15 Labels a, b, c, and d correspond to the spectral regions of the following bands: a-v2( 2 A ), b-v4( E ), c-v3( E ), and d-v3+v4( E ).

3.3. Definition of the two grids for the nuclear configurations.

  Because of the high-quality but costly basis sets and corrections beyond the RHF-UCCSD(T) level employed in this work to construct the PES, two different grids of nuclear configurations have been considered. A similar strategy was already applied for other molecules (see e.g. CH4, 99 H2CO, 100 and NH3[START_REF] Egorov | New Ab Initio Potential Energy Surfaces for NH3 Constructed from Explicitly Correlated Coupled-Cluster Methods[END_REF] ).

	The first grid (called Grid_1) was generated by varying simultaneously all the six coordinates
	defined by Eqs. (1) and (6). The  coordinate takes the following values (here in degrees): 90, 94,
	98, 102, 106, 108, 110, and 112. The two components 3 E b S  and 4 E b

4. RESULTS AND DISCUSSION 4.1. Equilibrium geometry and vibrational band origins.

  The equilibrium geometry (re) of CH3 was obtained as the stationary point of the fitted PES and allowed us to analyze the contributions of various corrections to re. According to Table1, there is a good agreement (error of 0.00006 Ang.) between the results at the CBS(TQ5Z) limit and those provided by the explicitly correlated F12x{x=a, b} methods. Through a series of corrections, the equilibrium geometry became closer to the value re = 1.076298 Ang. obtained in ref 15 from an empirical optimization of the ab initio PES.It is necessary to stress, however, that our theoretical re cannot be directly compared to the experimental data derived from the rotational constants because the ground vibrational state lies higher (by Zero-Point-Energy) than the global minimum of the PES. Our ab initio A=B=9.5820 cm -1 and C=4.7396 cm -1 rotational constants of the ground vibrational state agree within 0.0041 cm -1 and 0.0024 cm -1 with the empirical values from ref 68.

Table 1 . Equilibrium geometry (re in Ang.) of CH3 (D3h) at the RHF-UCCSD(T) level of the theory and after adding various corrections a

 1

	Basis	RHF-UCCSD(T)	+Rel.	+DBOC	+CCSDT +CCSDT(Q) +CCSDTQP
	ACVTZ	1.078220	1.078103 1.078243 1.078394	1.078458	1.078465
	ACVQZ	1.076552	1.076433 1.076574 1.076723	1.076787	1.076794
	ACV5Z	1.076025	1.075907 1.076047 1.076196	1.076259	1.076266
	CBS(TQ5Z)	1.075788	1.075670 1.075810 1.075958	1.076022	1.076029
	-F12a/CVQZ-F12	1.075849	1.075731 1.075872 1.076020	1.076083	1.076090
	-F12b/CVQZ-F12	1.075842	1.075724 1.075865 1.076013	1.076076	1.076083

Table 2 . Theoretical band origins of СH3 predicted from the ab initio PESs constructed using different orbital basis sets a

 2 All our PESs were constructed using the RHF-UCCSD(T) and RHF-UCCSD(T)-F12x{x=a,b} methods with the basis set displayed in Table2and including the following corrections: relativistic, DBOC, and high-order electronic correlations [CCSDT, CCSDT(Q), and CCSDTQP]. The details can be found in section 3.2.

	Band
	Exp. V1 V2 V3 V4 Γvib	ACVTZ ACVQZ ACV5Z CBS(TQ5Z) CVQZ-F12b CVQZ-F12a
	0 1 0 0	2 A 606.4531(6) c	8.690	2.377	0.538	-0.223	0.575	-1.143
	0 2 0 0	1 A 1288.0900(10) c 14.772	4.156	1.086	-0.167	0.850	-1.757
	0 3 0 0	2 A 2019.1657(11) c 19.894	5.701	1.646	0.016	1.050	-2.181
	1 0 0 0	1 A 3004.436(12) d, e 9.065	3.704	1.265	-0.194	-0.046	0.094
	0 0 1 0 E 3160.82118(6) f, h 10.430	4.125	1.398	-0.202	-0.224	-0.170
	Total RMS	13.270	4.151	1.244	0.177	0.664	1.356

b Exp. b -Calc. a b Only gas-phase IR data are used for validation: c -Yamada et al.; 68 d and e -Triggs et al. 70 and Zahedi et al.; 71 f and h -Amano et al.

Table 3 . Ab initio band origins of CH3 calculated in this work in comparison with previous studies a

 3 First 50 levels are displayed. Total list up to 7000 cm -1 is given in Supporting Information I. This work PES has Zero-Point Energy of 6482.697 cm -1 . The sources of the experimental data are the same as in Table2. Iic-MRCI(Q)/AVTZ by Schwenke; 12 IIic-MRCI(Q)/AVTZ by Medvedev et al.;[START_REF] Medvedev | Methyl radical: ab initio global potential surface, vibrational levels and partition function[END_REF] III -RHF-RCCSD(T)-F12b/CVQZ-F12 by Adam et al.[START_REF] Adam | Variationally Computed IR Line List for the Methyl Radical CH3[END_REF] Since the RHF-UCCSD(T)/CBS(TQ5Z) PES provided the most accurate results, it was used for a comparison with the previous published studies. Our vibrational labels are the four vibrational quantum numbers as well as the vibrational symmetry in the D3h point group. According to Table3, the calculations based on the ic-MRCI(Q)/AVTZ PES from ref 12 are comparable with our results obtained at the RHF-UCCSD(T)/ACVTZ level (see Table2). The best available band centers have been published in ref 15 using the partially spin restricted open-shell explicitly correlated CC method

	Band			Energy (cm -1 )			Exp. b -Calc.
	V1 V2 V3 V4 Γvib	I	II	III	This work	I	II	III This work
	0 1 0 0	2 A	596.32	591.70	602.43	606.68	10.13 14.75 4.02	-0.22
	0 2 0 0	1 A	1278.89	1266.20	1281.24	1288.26	9.20 21.89 6.85	-0.17
	0 0 0 1 E 1387.53 1388.39	1387.26	1390.66
	0 1 0 1 E 1996.99 1991.79	2002.22	2008.63
	0 3 0 0	2 A	2025.6	1994.20	2010.09	2019.15	-6.43 24.97 9.08	0.02
	0 2 0 1 E 2690.55 2673.95	2691.61	2700.01
	0 0 0 2	1 A	2748.17	2750.70	2739.64	2756.82
	0 0 0 2 E 2766.09 2767.81	2762.05	2773.08
	0 4 0 0	1 A	2829.01	2763.20	2776.86	2787.63
	1 0 0 0	1 A	2991.45	2988.50	3002.76	3004.63	12.99 15.94 1.68	-0.19
	0 0 1 0 E 3144.36 3142.62	3158.83	3161.02	16.46 18.21 1.99	-0.20
	0 1 0 2	2 A	3371.05	3367.10	3371.59	3386.62
	0 1 0 2 E 3388.51 3382.97	3391.11	3402.09
	0 3 0 1 E 3447.82 3407.54	3430.06	3440.16
	0 5 0 0	2 A	3686.04	3557.70	3569.95	3581.61
	1 1 0 0	2 A 3575.47	3572.80	3596.3	3603.44
	0 1 1 0 E 3715.98	3710.21	3736.97	3743.45

a b [RCCSD(T)-F12b/CVQZ-F12] and correspond to our RHF-UCCSD(T)/ACVQZ PES. According to the trial calculations in refs 90 and 91, the difference between RHF-UCCSD and RCCSD is generally smaller than the basis set error. To demonstrate this, a geometry optimization was performed for the ground 𝑋 ̃2 2

Table 4 . The contributions of the ab initio corrections beyond the RHF-UCCSD(T) level of the theory into band origins of CH3 a

 4

	Band
	V1 V2 V3 V4 Γvib	Relativistic	DBOC	CCSDT	CCSDT(Q) CCSDTQP	Total
	0 1 0	2 A	-0.095	0.413	-0.472	0.184	0.064	0.093
	0 2 0	1 A	-0.091	0.686	-0.967	0.176	0.095	-0.101
	0 0 0 E	0.294	0.004	-0.750	-0.672	-0.059	-1.184
	0 1 0 E	0.188	0.407	-1.233	-0.498	0.003	-1.133
	0 3 0	2 A	-0.038	0.888	-1.484	0.076	0.112	-0.445
	0 2 0 E	0.184	0.678	-1.740	-0.517	0.033	-1.362
	0 0 0	1 A	0.548	0.030	-1.516	-1.324	-0.116	-2.379
	0 0 0 E	0.573	0.015	-1.507	-1.347	-0.119	-2.385
	0 4 0	1 A	0.053	1.024	-2.008	-0.062	0.117	-0.876
	1 0 0	1 A	-0.040	-0.485	-2.486	-1.125	-0.084	-4.221
	0 0 1 E	-0.068	-0.397	-2.485	-1.011	-0.101	-4.061
	0 1 0	2 A	0.430	0.408	-2.009	-1.161	-0.056	-2.388
	0 1 0 E	0.454	0.406	-2.003	-1.180	-0.057	-2.380
	0 3 0 E	0.232	0.883	-2.269	-0.629	0.047	-1.735
	0 5 0	2 A	0.085	0.897	-2.605	-0.004	0.095	-1.531
	1 1 0	2 A	-0.106	0.126	-2.900	-0.714	0.010	-3.584
	0 1 1 E	-0.201	0.018	-2.975	-0.800	-0.032	-3.989
	0 2 0	1 A	0.412	0.668	-2.535	-1.182	-0.026	-2.664
	0 2 0 E	0.452	0.650	-2.515	-1.239	-0.035	-2.687
	0 0 0 E	0.795	0.067	-2.287	-1.978	-0.172	-3.574
	0 0 0	1 A	0.836	0.033	-2.273	-2.021	-0.178	-3.603
	0 0 0	2 A	0.838	0.029	-2.272	-2.025	-0.179	-3.608
	0 4 0 E	0.310	1.038	-2.812	-0.693	0.052	-2.104
	1 2 0	1 A	-0.192	0.176	-3.469	-0.890	0.022	-4.354
	1 0 0 E	0.262	-0.461	-3.272	-1.849	-0.147	-5.467
	0 2 1 E	-0.218	0.283	-3.490	-0.806	0.000	-4.229
	0 6 0	1 A	0.267	1.041	-3.013	0.779	0.115	-0.810
	0 0 1	2 A	0.247	-0.331	-3.286	-1.740	-0.166	-5.276
	0 0 1 E	0.224	-0.385	-3.264	-1.702	-0.161	-5.289
	0 0 1	1 A	0.228	-0.403	-3.261	-1.720	-0.162	-5.318
	0 1 0 E	0.667	0.415	-2.789	-1.836	-0.115	-3.657

  / (molecule•cm -2 ) units] is calculated by the well-

known expression (see e.g. ref 103):

Table 5 . Integrated intensities (Sv) of the strongest bands of CH3 at T = 296 K a

 5 is a statistical weight of the lower state; i E is the lower state energy; T is the temperature; A , and E irreducible representations in the D3h point group. Finally, the partition function calculated at T = 296 K is 715.96. All data were based on the developed ab initio line list. The integrated intensities are the sum of the intensities of the individual ro-vibrational lines of the certain band. The lines which intensities were bigger than 10 -27 cm/molecule were included into the sum. The integrated intensities of the displayed bands are bigger than 10 -20 cm/molecule. There is one more band centered at 5910.589 with the same vibrational assignment.

	QT is a partition function;

i g a b

Table 5 ,

 5 there are three hot bands which fall in the region of the fundamental band v2( 2 A ) hot bands. Above, the combination bands v1+v4( E ) and v3+v4( E ) are the strongest ones. The fundamental band v1( 1 A ) is not infrared active and is thus very weak.Our ab initio line list is in good agreement, both for line positions and line intensities, with the semi-empirical list,[START_REF] Adam | Variationally Computed IR Line List for the Methyl Radical CH3[END_REF] in particular in the spectral regions of v2(2

	A ). Two bands, namely, 2v2(

A ) and v3( E ) (see Figures. 4a

ACKNOWLEDGEMENTS

This work was supported by the joint "TEMMEX" project between the Russian Scientific Foundation (RSF, No. 22-42-09022) and Agence Nationale de la Recherche (ANR, French National Research Agency, Grant No. 21-CE30-0053-01). The support from the ROMEO computer center of Reims-Champagne-Ardenne as well as from the French-Russian collaboration program IRP SAMIA2 is also acknowledged.

Supporting Information I, total list of the ab initio band origins of CH3 calculated in this work up to 7000 cm -1 (Table S1); the effect of the ab initio corrections on the band origins of CH3 (Table S2); Supporting Information II, the C++ code for computing potential energies and dipole moment from the constructed PESs and DMS; the ab initio line list of CH3 with the complete quantum assignment up to 7000 cm -1 (ZIP); Supporting Information III, the ab initio points of the potential energies and dipole moment for the nuclear configurations of "Grid_1" and "Grid_2" calculated in this work (ZIP).

TOC graphic


Ben Elyes

Guillaume Hamida
email: guillaume.chelius@inria.fr

Jean-Marie Chelius
email: jean-marie.gorce@insa-lyon.fr

Gorce

Scalable versus Accurate Physical Layer Modeling in Wireless Network Simulations

In wireless networking, due to the high complexity of analytical and theoretical models, simulations are generally considered as the most convenient methodology for performance evaluation. Nonetheless, the physical complexity of the wireless medium induces a clear tradeoff between accuracy and scalability in a wireless network simulator design. In this paper, we focus on this tradeoff and study the impact of physical layer (PHY) modeling accuracy in the computational cost of simulations. We first discuss the main aspects of the wireless medium and briefly show how they have been handled in existing simulators. Then, we introduce a flexible and modular PHY simulation framework to analyze in more details their influence on the scalability of simulations.

Introduction

In wireless multi-hop networking, due to the highcomplexity of analytical and theoretical models, simulations are generally considered as the most convenient methodology to explore the behavior of protocols and distributed applications. Nonetheless, the complexity of the physical phenomena constituting the radio medium introduces a tradeoff between accuracy and computational cost in wireless network simulation. Several wireless network simulators have been proposed in the last years. Examples are NS-2, GloMoSim [START_REF] Bajaj | Simulation of large-scale heterogeneous communication systems[END_REF], JiST/SWANS [START_REF] Barr | Scalable Wireless Ad hoc Network Simulation[END_REF], GTSNetS [START_REF] Ould-Ahmed-Vall | Simulation of large-scale sensor networks using gtsnets[END_REF], etc. They all provide an advanced and complete simulation environment to investigate and evaluate networking protocols and wireless systems. However, the complexity of the wireless physical layer (PHY) leads to implementation choices during the simulators design. As a consequence, the PHY simulation accuracy varies drastically from one simulator to another.

As it has been highlighted in previous publications [START_REF] Cavin | On the accuracy of manet simulators[END_REF][START_REF] Heidemann | Effects of detail in wireless network simulation[END_REF][START_REF] Takai | Effects of wireless physical layer modeling in mobile ad hoc networks[END_REF], these variations largely impact on the results of a simulation. A correct modeling of the PHY layer is then crucial for confidence in the simulation results. Nonetheless, most of wireless network simulators are still based on inaccurate PHY models. The reason that generally prevails to justify this low-accuracy is scalability.

In this paper, we study the impact of PHY modeling accuracy on the computational cost of simulations. The question we raise is what is the real cost of PHY simulation accuracy ? We deliberately keep aside optimizations and scalability of the node and protocol aspects of simulations which have been the subject of other studies [START_REF]Simulation scalability issues in wireless sensor networks[END_REF]. To evaluate the PHY tradeoff, we introduce a flexible and modular PHY simulation framework, called WSNet 1 . We chose not to analyze the PHY tradeoff through the study of an existing simulator as none of them offer a sufficient diversity in PHY models. Comparing several existing simulators would not have helped, as they differ in many other aspects than PHY modeling.

The remainder of this paper is organized as follows. In section 2, we discuss the main aspects of wireless communications. Then, in Section 3, we review some common wireless network simulators and briefly compare their PHY modeling. Next, the WSNet simulation framework is presented in Section 4. In Section 5 we investigate the impact of the PHY modeling on the simulation speedup. Finally, the tradeoff involving accuracy and scalability is illustrated in Section 6.

Toward a precise physical layer modeling

For the sake of realism and confidence in simulation results, an accurate PHY modeling is a key point. In analytical studies as in simulations, the disk model has long prevailed. It relies on a set of strong assumptions: time stationarity

l ij (t) = l ij , independence l ij = f (x i , x j ), switched link (on/off) l ij ∈ {0, 1}, symmetry l ij = l ji , isotropy l ij = f (x i , d ij ) and homogeneity l ij = f (d ij )
, where l ij refers to the radio link between nodes i and j and d ij to the geometric distance between i and j.

The disk model provides the radio network with three axioms: the radio range is constant, the radio link is switched, and the network is interference free. If the asset of this model holds in its simplicity for both theoretical studies and simulations, the resulting simulations are far from realistic. Nonetheless, improving this model is not a trivial task as a hard tradeoff between complexity and realism holds. Basically, this model can be improved by relaxing either of the three previously stated axioms, as discussed below.

Radio range modeling

The range of a radio system is based upon the definition of a signal to noise ratio (SNR) threshold noted γlim . If the system is interference free, the range is a constant and the radio link is defined by:

l ij : Ω 2 → B = {0, 1} (x i , x j ) → l(x i , x j ) = 1 if γij ≥ γlim 0 else (1)
where the SNR γij is given by: γij = h ij • Pi Nj , where h ij is the path-loss and P i and N j are the transmission power and the noise level respectively.

The transceiver properties. They are the transmission power P i , the noise level N i , the antenna gain and its radiation pattern g i (θ, φ). Variations of P i , N i or antenna gains affect the spatial homogeneity assumption (and so far the symmetry), which means that all nodes not further have the same range. Note that a non-uniform noise level N j is highly probable for low cost small radio systems. The isotropy which is not statistically affected by these parameters does not further hold if the radiation patterns of the antennas are introduced according to:

h ij = g(x i , x j ) • g i (θ ij , φ ij ) • g j (θ ji , φ ji ),
where g(x i , x j ) is the propagation path-loss. It should be noted that 3D radiation pattern and 3D distribution of radio systems is pertinent for the simulation of small indoor environments for example.

Propagation models. The simplest model refers to the line of sight (LOS) scenario but in urban and indoor environments, more complicated scenarios occur due to shadowing and multiple paths. Two complementary approaches can be used to deal with propagation.

The former approach relies on a deterministic modeling of the wave propagation and provides fine simulations of any environment. The most usual algorithms are ray-tracing based [START_REF] Almers | Survey of channel and radio propagation models for wirelessmimo systems[END_REF] but discrete methods have been also proposed [START_REF] Gorce | A deterministic approach for fast simulations of indoor radio wave propagation[END_REF]. The high accuracy of these methods is definitely balanced by their high computational cost. Another limitation of purely deterministic models is that simulating one real environment is often too specific.Thus, the later approach relies on a statistical description complementing the deterministic model. A stochastic variable s ij is then introduced in the propagation path-loss to handle shadowing:

g ij = g ij • s ij .
The most usual model is the log-normal shadowing. A very challenging issue with shadowing is to introduce a spatial correlation between radio links which is not found in current simulators. Thus, s ij should be a spatial correlated stochastic process constant in time; it is not the case in some works as shadowing is often confused with fading. Fading refers to SNR time variations due to multi-path interference. It has a leading role in wireless communications and is introduced also as a stochastic variable

f ij : g ij (t) = g ij • s ij • f ij (t).
This late variable is not spatially correlated as it relies on small scale phenomenon. Meanwhile, it is a time variant parameter. Considering its temporal correlation may be highly relevant [START_REF] Pham | New cross-layer design approach to ad hoc networks under Rayleigh fading[END_REF] and represents another challenge for wireless simulators.

Radio link modeling

A frame error rate (FER) as a function of the mean SNR can substitutes for the SNR threshold of Equation 1. It derives from the bit error rate (BER) function which itself relies on the radio interface properties. l ij then relates to the probability of a successful transmission. Theoretical asymptotic expressions are well-known for various modulation techniques [START_REF] Wang | A simple and general parametrization quantifying performance in fading channels[END_REF].

Interference modeling

Interference disturbs the packet reception at the physical layer. It appears as a crucial point in PHY simulations as final results can be strongly influenced by the interference model [START_REF] Cavin | On the accuracy of manet simulators[END_REF][START_REF] Heidemann | Effects of detail in wireless network simulation[END_REF][START_REF] Takai | Effects of wireless physical layer modeling in mobile ad hoc networks[END_REF]. As we will see in Section 3, interference management is probably the point where current simulators differ the most largely. Sources of interference include nodes operating in the same frequency band or in different frequencies. The first type of interference is known as co-channel interference, while the latter is termed adjacent channel interference.

The most efficient approach for introducing interference consists in replacing the SNR by a signal to interference plus noise ratio, SINR, which can be derived according to:

γij = h ij • P i N j + k =i,j h kj • P k (2)
The proper derivation of the SINR requires the knowledge, at a given time, of all the signals which are concurrently received at a given receiver.

To be exhaustive, it should be noted that non linear receivers (with multi-user detection for instance) [START_REF] Andrews | Interference cancellation for cellular systems: a contemporary overview[END_REF] can outperform classical receivers in the presence of interference. In this case, it is necessary to compute the FER not from the general SINR, but rather from the vector of received powers at each node.

MIMO systems

Multiple input multiple output (MIMO) interfaces are very promising and their simulation can be assessed in two ways. The former but simplest approach is referred to as the node-based approach and exploits the single antenna framework described above by adjusting the FER function according to a MIMO specific model. This approach is however not relevant for realistic systems where the radiation patterns of the multiple antennas differ from each other as well as for correlated channels. The later, referred to as the antenna-based approach, is more powerful and simulates the channel state for each antenna-to-antenna link, separately. For instance, in a 2 × 2 MIMO system, four channels are simulated. It is very important to take care about the correlation between these channels to have accurate simulations. In this context, modeling spatially correlated shadowing and time correlated fading is of great importance.

State of the art

Numerous wireless network simulators have been developed and are concurrently used in the wireless network research domain. In this section, the PHY modeling of four widely used simulators is presented and briefly investigated. The simulators are NS-2, GloMoSim [START_REF] Bajaj | Simulation of large-scale heterogeneous communication systems[END_REF], JiST/SWANS [START_REF] Barr | Scalable Wireless Ad hoc Network Simulation[END_REF] and GTSNetS [START_REF] Riley | Large-scale network simulations with gtnets[END_REF]. This review is far from being exhaustive and has been deliberately restricted for space restriction reasons. Other available simulators with wireless support are OPNET, OMNeT++, etc.

Let us consider the PHY modeling as implemented in these simulators. Due to space limitations, we restrict our focus to one of the most crucial aspects for cross-layer and protocol studies, namely interference modeling. The first step toward interference evaluation is to identify which signals are interfering with each other in order to assess the terms in the denominator of eq. 2, on the basis of timing considerations only. This set of interfering signals can be very large for large scale simulations. As a consequence, various simulators rather limit the range at which any signal can propagate and thus can interfere. In other words, disregarding the radio range model effectively used for the received signal strength computation, the simulator does not generate receptions at nodes further than a given range from the source. Consequently, the considered source cannot induce interference at nodes further than this range. This optimization, called Limited interference model as opposed to the Full interference model, is depicted in Figure 1. It privileges scalability at the cost of accuracy. This optimization is implemented in JiST/SWANS, GloMoSim, etc.

Regarding the SINR computation, several strategies have been investigated and implemented in existing simulators. They are all variations of eq. 2 regarding timing granularity. They induce a varying level of realism, precision but also complexity. They are summed up in Figure 2: It is obvious that these interference models offer different levels of complexity, accuracy or realism. In order to better evaluate the computational cost of interference modeling, we introduce the WSNet PHY simulation framework.

• Fig 2(
Note that none of these simulators supports multichannel systems. For this purpose, the interference model in (2) should be replaced by:

γij = h ij • P i N j + k =i,j α ik • h kj • P k , ( 3
)
where α ik stands for the rejection factor between the channels associated with signals i and k.

The WSNet PHY simulation framework

The WSNet PHY simulation framework is well-designed for the purpose of evaluating the computational cost of the PHY modeling accuracy. The design of WSNet has been mainly guided by two constraints: modularity and flexibility. Thanks to modularity, WSNet can offer a wide range of PHY models, from ideal ones to pseudo-realistic ones.

WSNet operates either as a complete event-driven wireless network simulator, including the simulation of both the wireless nodes and the radio medium, or as a physical layer simulator with nodes being simulated by external simulators in a distributed simulation framework [START_REF] Chelius | Worldsens: development and prototyping tools for application specific wireless sensors networks[END_REF].

PHY layer simulation

In WSNet, a radio signal is characterized by the following values: an emission power, a symbol rate, a modulation scheme and a channel number. The PHY layer is abstracted as a combination of independent blocs which represent radio medium properties or hardware components. The radio simulation is built upon the following blocs: (i) a propagation bloc dealing with propagation aspects including pathloss, shadowing and fading according to the definitions provided in section 2; (ii) an interference bloc which implements eq. 3; (iii) modulation blocs which provide SINR-to-BER conversions; and (iv) antenna blocs which implement antenna properties: loss, gain, radiation pattern, orientation and position. Note that WSNet offers a 3D spatial representation.

These blocs are black-boxes with well-defined interfaces. At simulation runtime, they are linked to dynamic libraries which effectively implement the interface. As an example, for the propagation bloc, a dynamic library implements the free space propagation model while another one provides a Rayleigh channel and a third one inputs pathloss values from an external file or an external propagation tool. More generally, this architecture is an opportunity to offer a wide range of radio medium models, from a basic ideal physical layer with no interference and no path-loss to a more realistic one including a Rayleigh channel, multiple frequencies, complex modulation schemes and smart antennas. It is thus a good opportunity to study the computational cost of using pseudo-realistic radio models in term of simulation complexity.

Radio range simulation

When considering a transmission initiated by an antenna i, WSNet generates two events at each receiving antenna. The first one notifies the apparition of the signal while the second one notifies its end. Given a receiving antenna j, the two events are scheduled according to the distance between the emitting and receiving antennas, the radio propagation speed and the signal symbol rate. By default, the set of receiving antennas is extended to all antennas of the network. Thus, by default, WSNet operates in a full interference mode. However, the simulation configuration may indicate a limited interference mode together with a maximum range after which receptions and interference are not propagated. This range is independent from the radio range model implemented in the propagation bloc. For each reception, i.e. each receiving antenna j, the received signal strength is computed according to the operations and properties depicted in Figure 3: (1) the emitting antenna noise N i , (2) the emitting antenna gain g i (θ ij , φ ij ), (3) the signal attenuation including shadowing and fading effects g ij (t), (4) the receiving antenna gain g j (θ ji , φ ji ) and ( 5) the receiving antenna noise N j .

Interference simulation

Upon signal reception, after the received signal strength computation, interference is considered: SINR values are evaluated according to eq. 3. WSNet divides the signal into slots and computes a SINR value for each of the slots. In order to trade accuracy for scalability and vice-versa, the number of slots can be tuned from 0 -no interference -up to the packet length -one SINR value for each signal byte. In each slot, WSNet computes a SINR according to the weighted cumulative interference model described in Channels have been introduced in WSNet to support multi-channel and MIMO systems. Channels divide the radio resource in sub-resources which mutually interact. As an example, channels can be used to model frequencies in a FDMA (Frequency Division Multiple Access) network as well as codes in a CDMA (Code Division Multiple Access) one. We do not explicit what channels are but we rather explicit the correlation, i.e. how they interfere with each others. While computing the interference on channel i, the contribution of a signal emitted on channel j is weighted by a correlation factor α ji as given in eq. 3. α ji values are provided by the interference bloc.

Radio link simulation

Upon signal reception and after each slot SINR computation, the impact of interference on the radio signal is evaluated by computing the BER values. The SINR to BER conversion is provided by the modulation bloc that characterizes the radio signal. As for the SINR, a BER value is associated with each slot of the radio signal. Based on the BER values and the slot lengths, a FER is finally computed. Contrarily to GloMoSim, WSNet forwards the radio signal to the receiving antenna whatever the FER value and whenever the first error arises. Indeed, in the majority of cases, errors are not detected until the signal is completely received or lost.

An optional feature of WSNet is to randomly introduce errors in the radio signal with a probability equal to the computed BER. In this case, an erroneous signal is transmitted to the receiving antenna. Given this feature, it becomes possible to study precisely the performance of channel coding and error correction algorithms through their implementation in the node simulation part. More generally, CDMA schemes can be either simulated statistically using channels and correlation factors or studied precisely through their implementation and the error introduction feature.

Impact of the physical layer modeling

We now study the cost of PHY modeling accuracy using WSNet, focusing on the three aspects that have been developed in Section 2: radio range modeling, interference modeling and radio link modeling. As an estimation of the com-putational cost, we use the speedup metric. The speedup of a simulation is the ratio between the logical simulated time and the effective simulation time. Unless specified otherwise, the simulations consist in 100 nodes deployed randomly in a 100 × 100 × 0 area. Each node emits a 100B/s Constant Bit Rate (CBR) broadcast traffic through an IEEE 802.15.4 868Mhz compliant radio. No MAC protocol is executed in order to focus on the PHY simulation. The presented speedups are averaged over 100 simulations with the confidence intervals being omitted as they have been found to be negligible.

Radio range models

To assess the impact of radio range models on the speedup of simulations, we compare four common propagation models: a disk model (no path-loss) with transmission ranges of 60m and 200m, the free-space model, the two-ray ground model and a Rayleigh channel with a pathloss exponent of 2 for the last three cases. Interference and modulation support are disabled. The results are drawn on Figure 4. Quite obviously, a better accuracy in the radio range model induces a decrease in the simulation speedup. There are two reasons for this overhead. The first one is the computational cost associated with the radio range model. It takes more time to compute a complex path-loss with random variables than a single distance. The second one is an increase in the number of receptions generated by a single emission. Consider for example the Range(60) and Range(200) models. They have the same computational cost but the number of receivers is larger with a range of 200 meters. In a realistic modeling, this overhead is even higher as the signal reaches all nodes. If the first overhead can hardly be reduced, the second one can be reduced using the limited interference model.

Interference models

Keeping the Friis propagation model, we now analyze the impact of interference simulation on the scalability of simulations. In a first step, we make the number of SINR slots vary from 0, no interference simulation, to n -one SINR per byte. Results are depicted on Figure 5. In a second step, we evaluate the cost of multi-channel support and adjacent channel interference computation. We vary the number of simultaneously simulated channels from 1 to 16. Results are depicted on Figure 6. Support of a cumulative interference model induces a strong overhead as the speedup is roughly reduced by 50% from no interference simulation to a single SINR computation. However, in interference modeling, it seems that the first step is the one that costs the most. Indeed, from a single SINR to a per-byte computation, the extra overhead remains quasi constant for each refinement: from 4% for two SINR values to 45% for n values. In consequence, once the interference price has been paid, there is no strong reason to decline paying for a better accuracy. This same conclusion holds for multi-channel support as its cost is negligible compared to the one associated with an accurate interference modeling.

Radio link models

In addition to interference, we now evaluate the cost of an accurate radio link modeling by integrating three modulation models in the simulations: SINR threshold, bpsk and oqpsk. The results are shown in Figure 7. As a BER value is derived from each computed SINR, the overhead induced by the radio link model is function of the interference modeling accuracy. As in Section 5.1, a realistic model induces a higher overhead: the computational cost of the erf c function used in the bpsk model is much higher than a simple comparison between a SINR and a threshold. However, from a SINR threshold to a bpsk function, the speedup only decreases by 4 to 17% depending of the number of SINR slots. This price may clearly be worth the gain.

Limited interference model

We finally study the speedup gain that can be achieved through the use of a limited interference model. If this model decreases the accuracy of the PHY layer simulation, it is supposed to increase the scalability of simulators, justifying its use in GloMoSim, JiST/SWANS, etc. For this evaluation, we consider a network of 2000 nodes randomly deployed in a network of size 1000m × 1000m × 0m. No interference nor modulation is simulated. We compare the speedup achieved using the limited interference model for various maximum interference/reception ranges to the one achieved using the full interference model. Results are presented in Figure 8. In the same Figure, we also highlight the interest of using an optimized structure to maintain the node positions and to search for nodes in reception/interference range. Indeed, the complexity of range searching using a list is O(n) and thus does not scale well in the number of nodes. In consequence, several optimized data structures have been proposed for the problem of range searching [START_REF] Bentley | Data structures for range searching[END_REF]. Two structures have been implemented in WSNet and are compared here: a db-tree (dynamic balanced tree) and a grid. A more complete study of optimized data structures for the range search problem can be found in [START_REF] Bentley | Data structures for range searching[END_REF]. As it was already pointed out in Section 5.1, a reduction in the number of receivers induces a large gain in the simulation speedup. For a limited interference range of 30 meters, with the correct data structure, the speedup can be increased by a factor of 332. In large-scale networks, this gain remains high even for larger ranges. With a limited range of 120 meters, a gain of 37 is still achieved. It is quite obvious that the limited interference model is a must to ensure a high scalability. However, its use raises the issue of determining a correct range (i.e., the distance limit for the propagation). As a solution, [START_REF] Ji | Improving scalability of wireless network simulation with bounded inaccuracies[END_REF] proposes an empirical method to derive a range with a limited impact on the simulation accuracy.

Scalability versus accuracy

Finally, we propose a short case study to clearly summary the impact of PHY modeling on the scalability and accuracy of simulations. We consider a varying number of nodes randomly deployed in a 200×200×0 area. Each node emits periodically a hello packet (100B/s) through an IEEE 802.15.4 868Mhz compliant radio. We consider three metrics: the speedup, the average number of discovered neighbors and the corresponding average number of connex components. The first metric allows to assess the impact of radio models on the scalability while the two latter metrics show the impact of the PHY layer modeling on the evaluation of higher level protocols in terms of network connectivity. We performed the same set of simulations with various PHY models. Starting from (i) a basic disk model (range = 50) with no interference nor modulation, we slightly increased the PHY simulation accuracy through the introduction of (ii) path-loss and (iii) log-normal shadowing with a transmission power of 0dBm. Then, we considered cumulative interference with a limited interference model (range of 50) and BPSK modulation for (iv) one slot per packet and (v) n slots per packet. We finished with the (vi) full interference model. Simulation results are shown on Figure 9.

As pointed out in Section 5.1, we first notice that the addition of path-loss and fading to the range model induces a low impact on the speedup. Still considering the speedup, the main gap occurs with the introduction of interference and modulation support. Next, as previously observed in Section 5.2, considering n SINR slots instead of 1 induces a regular extra overhead. Finally, the full interference model does not induce too much overhead as the network area is small in this case. Regarding the behavior of higher level protocols under variable PHY models, it can be observed that the average number of discovered neighbors varies systematically. If the largest gap also occurs when interference and modulation are introduced, the results still degrade when the accuracy of the interference modeling is increased. This degradation is more important with a high number of nodes. As the number of nodes increases, so does the level of interference and its impact on higher protocols. Sensibility to PHY accuracy increases with the network size.

Considering the average number of connex components, we observe that the logical network connectivity varies with the PHY modeling. As expected, the use of a path-loss or a log-normal shadowing with no interference increases the network connectivity (i.e, low number of connex components). However, this connectivity decreases significantly (i.e, high number of connex components) when increasing the interference model accuracy, from 1 to n slots, from a limited to a full model with modulation support. This underlines the importance of using accurate radio medium models when evaluating protocols, especially in large-scale networks.

Conclusion

In this paper, through the introduction of the WSNet PHY simulation framework, we have studied the impact of PHY modeling on the scalability and accuracy of wireless network simulations. To the question raised in the introduction: what is the real cost of PHY simulation accuracy ?, we can give the following conclusions. Accuracy has a variable cost depending on the considered PHY aspect. On one hand, for radio range modeling, the cost remains low, especially compared to the gain that can be achieved using realistic models. On the other hand, interference simulation induces a real high overhead. However, the extra cost remains regular when the simulation accuracy gets refined. Given the overhead of interference modeling, the cost added by accurate radio link models remains low. So as the support of multiple channels in the simulation. As a consequence, it is prejudicial to simulate interference but no realistic modulation schemes.

Given these trends, it remains up to the user choice to trade PHY accuracy for a desired scalability. However, we must keep in mind that, as highlighted in this work, the validity of simulation results highly depends on the choices made for the PHY modeling. In particular, if the use of optimizations as a limited interference model drastically increases the scalability of simulators, it does not solve the accuracy versus scalability tradeoff. The optimization performance and correctness now depends on the chosen limited range. In the future, we plan to investigate the limited interference model in more details. In parallel to the empirical model proposed by [START_REF] Ji | Improving scalability of wireless network simulation with bounded inaccuracies[END_REF], we think that a stochastic theoretical analysis can help in determining an adequate range. The goal here is still to enhance the scalability of simulators without affecting the confidence in the simulation results.

Figure 1 .

 1 Figure 1. Limited versus full interference model (T: transmitter, R: receiver and I: interfering nodes).

Figure 2 .

 2 Figure 2. SINR computation strategies at the node receiving the packet B.

Figure 3 .

 3 Figure 3. Signal reception

  Figure 2(b). As a consequence, given the configured number of slots, WSNet accuracy ranges from Fig. 2(b) to Fig. 2(d).

Figure 4 .

 4 Figure 4. Impact of propagation models.

Figure 5 .

 5 Figure 5. Impact of co-channel interference.

Figure 6 .

 6 Figure 6. Impact of multi-channel interference.

Figure 7 .

 7 Figure 7. Impact of modulation on speedup.

Figure 8 .

 8 Figure 8. Space partitioning methods.

Figure 9 .

 9 Figure 9. Impact of the PHY modeling accuracy on simulation speedup and results.

  a): only one SINR value is computed for the whole packet. The noise value considered for the SINR computation is the cumulative power of the signals that interfere with the considered signal at reception time. This method is the one implemented in JiST/SWANS.

• Fig 2(b): only one SINR value is computed for the whole packet. The noise value considered for the SINR computation is the cumulative power of all interfering signals weighted by the interference duration. • Fig 2(c): only one SINR value is computed for the whole packet. The noise value considered for the SINR computation is the power of the strongest interfering signal. This method is implemented in NS-2 and is close to the collision model of GTSNetS. • Fig 2(d): several SINR values are computed for the packet; whenever the set of interfering signals changes, a new SINR value is computed. The noise value considered for each SINR computation is the cumulative power of all concurrent interfering signals. This method is implemented in GloMoSim.

WSNet is available at: http://wsnet.gforge.inria.fr.