Robert R Gamache
email: robert_gamache@uml.edu

Bastien Vispoel

Michael M Rey

Vladimir Tyuterev

Alain Barbe

Andrei Nikitin

Oleg L Polyansky

Jonathan Tennyson

Sergei N Yurchenko

Attila G Császár

Michaël Rey

Tibor Furtenbacher

Valery I Perevalov

Sergei A Tashkun

Partition sums for non-local thermodynamic equilibrium conditions for nine molecules of importance in planetary atmospheres

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

Introduction

When studying planetary atmospheres via spectroscopy, the radiative transfer equations are solved by following the radiation as it passes through roughly isothermal-isobaric layers, where at each layer the incoming and outgoing radiations are determined. At each layer, the line intensity and line shape parameters must be determined for the pressure and temperature of the layer based on reference values, which can be found in molecular spectroscopic databases such as HITRAN [START_REF] Gordon | The HITRAN2020 molecular spectroscopic database[END_REF] or GEISA [START_REF] Delahaye | The 2020 edition of the GEISA spectroscopic database[END_REF]. As long as the pressure and temperature are not too low, the number of collisions per unit time will be sufficient for the populations of the quantum energy states to follow the Boltzmann law and the molecular velocity distribution to remain Maxwellian. This condition is called Local-Thermodynamic Equilibrium (LTE). However, if the conditions change such that the number of collisions per unit time is too low, one finds that the populations deviate from the Boltzmann law while the molecular velocity distribution remains Maxwellian. This occurs whenever collisions between molecules are infrequent enough that the radiative lifetimes of the important vibrational states of a molecule become comparable to the mean time between collisions. This condition is labeled Non-Local Thermodynamic Equilibrium (NLTE).

The change in conditions can be understood by considering the hydrostatic equation, which shows that atmospheric pressure and density decrease exponentially as a function of altitude. Kinetic theory states the collision frequency of a molecule is proportional to the molecule's cross-section, mean speed, and the density; the first two depending on the molecule in question and the latter on the conditions (p, T, and z). Hence, the number of collisions per unit time decreases exponentially as well, resulting in the conditions needed to maintain LTE no longer existing. At this point, the population of states is not given by the Boltzmann distribution but can be described by NLTE modeling. The population of energy states is crucially important in modeling atmospheres.

Considering the role of CO2 on the energy budget of Earth's middle atmosphere, the effects of modeling with or without NLTE can be understood. The 2 vibrational band of CO2 at 15 m is the dominant source of cooling by transferring radiation to space and the near-IR bands contribute substantially to the heating rate during the day by the absorption of incoming solar radiation. The more (less) the population of states in the 15 m band means the more (less) radiation transferred to space (cooling / less cooling); the more (less) the population of states in the near-IR bands the more (less) radiation absorbed (heating / less heating). Thus, moving molecules in or out of these states will affect the energy budget of the middle atmosphere.

Because the change from LTE to NLTE conditions affects the populations of these states, NLTE modeling must be done in order to obtain the population data to determine the correct energy budget.

Strictly speaking, perfect LTE or NLTE conditions are seldom encountered in planetary atmospheres. As shown by [START_REF] López-Puertas | Studies of Solar Heating by CO2 in the Upper Atmosphere Using a Non-LTE Model and Satellite Data[END_REF], the atmosphere is in a state that might be called partial LTE, in which the velocity distribution of molecules is Maxwellian but only some groups of levels obey the Boltzmann relation. If one follows a volume of gas and traces its state as the density decreases, then generally one would observe a series of partial LTE states, in which ever fewer levels follow the Boltzmann relation [START_REF] Kutepov | Solution of the non-LTE problem for molecular gas in planetary atmospheres: superiority of accelerated lambda iteration[END_REF]. In conditions of NLTE, collisional energy exchange alone, under conditions of lower atmospheric pressure, no longer determines the internal energy level populations because of the increased importance of radiative and chemical pumping and radiative losses. The partition of population is then characterized by a temperature different from the local kinetic temperature. For these conditions, the solution of the radiative transfer equations is preceded by solving the statistical equilibrium equations (SEE), which consider a large number of excited states of different molecular species coupled by a variety of collisional energy exchange processes and by band overlapping, together with a detailed model of atmospheric stratification.

When solving the radiative transfer equations in conditions of NLTE, the line intensities must be corrected to account for the change in the populations of states. This corrections is done using a partition function appropriate for NLTE conditions, as shown by (D. P. [START_REF] Edwards | The non-LTE correction to the vibrational component of the internal partition sum for atmospheric calculations[END_REF]. Using these corrections Lópes-Puertas et al. (López-Puertas et al., 1986a;López-Puertas et al., 1986b) studied the 15 m bands of CO2 in the middle atmosphere and clearly demonstrated that NLTE corrections were necessary to model heating and cooling rates in the middle atmosphere.

Hummer et al. [START_REF] Däppen | The Equation of State for Stellar Envelopes III. Thermodynamic Quantities[END_REF][START_REF] Hummer | The Equation of State for Stellar Envelopes I. An Occupation Probability Formalism for the Truncation of Internal Partition Functions[END_REF][START_REF] Mihalas | The Equation of State for Stellar Envelopes II. Algorithm and Selected Results[END_REF] formulated the equation of state for materials in stellar envelopes at low densities by using an occupation probability formalism that accounts for the perturbations of bound states by both neutral and charged perturbers. The formalism was then applied to the computation of the thermodynamic properties of a partially ionized (and/or dissociated) multicomponent gas. [START_REF] Kutepov | Solution of the non-LTE problem for molecular gas in planetary atmospheres: superiority of accelerated lambda iteration[END_REF] presented a general formulation of the multi-level rotation-vibration NLTE problem for a mixture of radiating molecular gases in a planetary atmosphere, treating explicitly the coupling of molecular energy levels by collisionally-induced energy-transfer processes and by band overlap. More recently, [START_REF] Young | Non-local thermodynamic equilibrium transmission spectrum modelling of HD 209458b. A&A. 641. Table 3 Term value cutoff needed to insure convergence of vibrational partition sums for methane isotopologues. Maximum temperature (K) of converged Qvib Maximum term value in cm -1 Number of term values for[END_REF] demonstrated that exoplanet upper atmospheres are low-density environments, where radiative processes can compete with collisional ones and introduce NLTE effects into transmission spectra. By adapting the spectral synthesis code Cloudy for NLTE calculations, they produced an atmospheric structure and atomic transmission spectrum in both NLTE and LTE for the hot Jupiter HD209458b. Their results demonstrated that individual spectral lines in the NLTE spectrum exhibit up to 40% stronger absorption that the corresponding lines in the LTE spectra, demonstrating the need for NLTE modeling of spectra [START_REF] Fossati | A data-driven approach to constraining the atmospheric temperature structure of the ultra-hot Jupiter KELT-9b[END_REF] as discussed above, recognizing the need to observationally constrain the atmospheric temperature-pressure (TP) profile of exoplanets, developed a data-driven approach to constrain the TP profile of the ultra-hot Jupiter KELT-9b by fitting synthetic spectra to the observed H and H lines and identified why self-consistent planetary TP models are unable to fit the observations. They constructed 126 one-dimensional TP profiles varying the lower and upper atmospheric temperatures, as well as the location and gradient of the temperature rise. For each profile, transmission spectra were computed accounting for non-local thermodynamic equilibrium effects. Their results demonstrated that the assumption of LTE overestimates the level populations of excited hydrogen by several orders of magnitude, indicating the need for including NLTE effects when modeling exoplanet atmospheres.

Laboratory measurements are also being made to elucidate NLTE effects that are essential for determining the properties of Hot-Jupiter-type exoplanets [START_REF] Bernath | Molecular opacities for exoplanets[END_REF][START_REF] Swain | The presence of methane in the atmosphere of an extrasolar planet[END_REF][START_REF] Tinetti | Spectroscopy of planetary atmospheres in our Galaxy[END_REF]. [START_REF] Dudás | High-temperature hypersonic Laval nozzle for non-LTE cavity ringdown spectroscopy[END_REF] have developed a novel approach to produce NLTE IR spectra using a hypersonic expansion probed by cavity ringdown spectroscopy (CRDS). They have used this method to study both the CO and CH4 molecules and have extracted from the work the rotational and vibrational temperatures necessary to characterize the NLTE.

In this work, the determination of the internal vibrational and rotational partition functions needed for NLTE computations is undertaken for nine molecules that are often observed in NLTE conditions and are abundant in planetary and exoplanet atmospheres: H2O, CO2, O3, N2O, CO, CH4, NO, NO2, and OH. The methods of calculation are reviewed for each molecule/isotopologue. A code, NLTE_TIPS_2021.FOR, that can quickly recall the vibrational NLTE and rotational LTE partition sums is then discussed.

Calculations of the vibrational and rotational partition sums

The calculation of the total internal partition sum (TIPS) is discussed for both LTE and NLTE conditions. For LTE, the most accurate TIPS are made by direct summation over very large sets of rovibrational term values, which usually come from first-principles calculations.

However, for many of the molecules in question, these large sets of term values are not available, and the product approximation [START_REF] Gamache | Total Internal Partition Sums in the Temperature Range 70-3000K: Atmospheric Linear Molecules[END_REF][START_REF] Gamache | Total Internal Partition Sums for Molecules in the Terrestrial Atmosphere[END_REF] must be made. In NLTE applications, the rotational and vibrational partition sums must be determined for a number of layers in the atmosphere. When first-principles data are available for the molecules, e.g., for H2O and CO2, the calculations of the NLTE partition sums can be made. However, the summing over the energies is time-consuming and must be done for each atmospheric layer since the correct vibrational temperature must be used for each state. Note also that if the atmospheric conditions change, the sums must be redone. A faster approach is to apply the product approximation to determine Qv and Qr separately. In this study, the product approximation is used for all the molecules.

LTE conditions

Under conditions of local thermodynamic equilibrium, the line intensity for a transition from initial state i to final state f at temperature T is given by

S if T ( ) = I a A if 8p c w if 2 g i e -h c E i / k T 1-e -h c w if k T ae è ç ö ø ÷ Q T ( ) (1)
where Ia is the natural abundance of the isotopologue (note that we are following the HITRAN [START_REF] Gordon | The HITRAN2020 molecular spectroscopic database[END_REF] definition of intensity that includes the isotopic abundance, for other databases that do not include this factor Ia should be set to 1.), Aif (s -1 ) is the Einstein Acoefficient, gi is the lower state statistical weight, Ei is the energy (term value) of state i, (lower state energy in cm -1 ), h, and k are the Planck and Boltzmann constants, respectively, c is the speed of light in vacuum, if is the wavenumber of the transition, and Q(T) is the total internal partition sum at temperature T. The total internal partition sum is given by a sum over all quantum states of the molecule,

Q T ( ) = g s e -h c E s k T all quantum statess å (2).
For the LTE condition, we have the Maxwell-Boltzmann distributions in thermal equilibrium at a temperature T(K), for which the number of molecules [cm -3 ] in a rovibrational state, nvJ, labeled by v and J, of energy EvJ [cm -1 ] and of statistical weight gvJ is

n vJ T ( ) = N g vJ e -h c E vJ k T Q T ( ) (3)
where N [molecules cm -3 ] is the molecular number density. One can write the population probability as

W vJ * T ( ) = g vJ e -h c E vJ k T Q T ( ) (4)
with the condition

W vJ * T ( ) v ,J å = 1 (5).
In this study, the molecules are in the electronic ground state and rotation and vibration are the only quantized motions, i.e., tunneling, internal rotations, etc. are not considered. Assuming no rotational-vibrational interactions, the wavefunction of the state can be written as a product of the vibrational wavefunction times the rotational wavefunction and the energy is given as a sum of vibrational and rotational energies,

E vJ = E v + E J (6).
Equation ( 2) can be rewritten as

Q T ( ) = g vJ e -h c E v + E J ( ) T k v , J å = g v e -h c E v ( ) k T v å ´gJ e -h c E J ( ) k T J å (7)
where gvJ = gv gJ. This is called the product approximation (PA) and TIPS can be written as

Q T ( ) = Q v T ( ) Q r T ( ) (8)
where the vibrational partition sum is

Q v T ( ) = g v e -h c E v k T all vibrational states å (9)
and the rotational partition sum is

Q r T ( ) = g J e -h c E J k T all rotational states å (10).
Population probabilities can be written, similar to Eq. ( 4), for the vibrational and rotational states.

Using Eq. ( 3), Eq. ( 1) can be rewritten as

S if T ( ) = I a A if 8p c w if 2 n i T ( ) N 1-e -h c w if k T ae è ç ö ø ÷ (11).
Note that in most cases the induced emission (IE) term in Eqs. ( 1) or ( 11), 1-e

-h c w if k T ae è ç ö ø ÷ , is
approximately one. At 296 K the IE term for the 183 GHz (6.104 cm -1 ) line of H2O is 0.0302, however, for an IR transition, say at 2000 cm -1 , the IE term is 0.99996. The same terms at 1000 K are 0.009 and 0.9490, respectively. So, often the IE term can be neglected when the wavenumber of the transition is in the infrared or higher.

NLTE conditions

As noted first in the introduction, a gas can be said to be in a state of NLTE when the population of the levels deviate from the Boltzmann law while the molecular velocity distribution remains Maxwellian. The translational, rotational, and vibrational levels of a molecular system adjust to equilibrium at very different rates. As a consequence, for most of the atmospheric species in the upper atmosphere, the translational and rotational levels are in LTE, while the vibrational ones are not. For such conditions, the populations of the states are no longer determined by simple Boltzmann statistics, but rather by a mixture of collisional reactions and exchange of vibrational quanta. The former can involve both vibrational-thermal and vibrational-vibrational exchanges, and both can involve a multiplicity of levels and exchange partners. [START_REF] Edwards | The non-LTE correction to the vibrational component of the internal partition sum for atmospheric calculations[END_REF] has shown that this state can be described with a similar probability formalism as LTE by introducing at each point a set of vibrational temperatures, T, in place of the kinetic temperature, T. The vibrational temperatures are defined by the vibrational level populations determined by the SEE models. The model includes vibrationalthermal (V-T) and vibrational-vibrational (V-V) collisional processes, absorption of solar radiation, exchange of photons between atmospheric layers, and spontaneous emission (for details see [START_REF] Edwards | The non-LTE correction to the vibrational component of the internal partition sum for atmospheric calculations[END_REF] and references therein). Once these populations are known, the vibrational temperature for the vibrational state labeled j is defined with reference to the ground state (subscript 0)

where E j , , g j , are the energy, number of molecules in, and statistical weight of the vibrational state j, n0 and g0 are the number of molecules in and statistical weight of the ground vibrational state.

Equation ( 8) can still be used to determine the rotational partition sum, Qr(T), since the rotational states are still in LTE. Because the vibrational levels are in NLTE, the usual expression for the vibrational partition sum, see for example Herzberg's harmonic approximation (HA) [START_REF] Herzberg | Molecular Spectra and Molecular Structure II. Infrared and Raman Spectra of Polyatomic Molecules[END_REF] or including anharmonic corrections to the vibrational levels [START_REF] Gamache | Total internal partition sums for 166 isotopologues of 51 molecules important in planetary atmospheres: Application to HITRAN2016 and beyond[END_REF], cannot be used. Qv(T) in Eq. ( 9) must now be determined via a direct sum, where the appropriate vibrational temperatures must be used. Rewriting Eq. ( 9) as ( 13)

where the tilde labels a NLTE quantity, T is the kinetic temperature and T j is the vibrational temperature for state j. The sum goes up to some maximum vibrational level, vmax, that ensures convergence of the sum. Due to the advent of extensive first-principles calculations, the set of js are known for many molecules. The T j are determined for each vibrational state by chemical-dynamical models (see (D. [START_REF] Edwards | Non-local thermodynamic equilibrium studies of the 15-μm bands of CO2 for atmospheric remote sensing[END_REF][START_REF] Edwards | The non-LTE correction to the vibrational component of the internal partition sum for atmospheric calculations[END_REF][START_REF] López-Puertas | Studies of Solar Heating by CO2 in the Upper Atmosphere Using a Non-LTE Model and Satellite Data[END_REF]López-Puertas et al., 1986a;López-Puertas et al., 1986b)). It should be noted that many vibrational states are not involved in the chemistry and dynamics that shifts population, so for these states, T j is the kinetic temperature.

The Non-LTE correction to the line intensity

Under conditions of LTE, radiative transfer theory for the two levels of a transition defines a spectral line intensity in terms of the Einstein-B coefficient for absorption for a transition f  i as [START_REF] Gamache | Extension of the Hitran Database to Non-LTE Applications[END_REF])

S if = h w if c B if n i 1- g i n f g f n i ae è ç ö ø ÷ ( 14
)
where Bif is the Einstein coefficient for absorption, ni and nf are the populations of the lower and upper states, respectively, gi and gf are the level statistical weights of the lower and upper states, respectively, and if the transition wavenumber.

The ratio of the populations of the upper and lower levels,

n f n i = N g f e -h c E f k T kin ( ) Q N g i e -h c E i k T kin ( ) Q = g f g i exp -h c E f -E i ( ) k T kin ae è ç ö ø ÷ = g f g i exp -h c w if k T kin ae è ç ö ø ÷ (15),
can be rearranged to give the last factor in the parentheses in Eq. ( 14), which is designated as the Boltzmann factor  at the local kinetic temperature Tkin,

G = g i n f g f n i = exp - h c w if k T kin ae è ç ö ø ÷ (16).
The line intensity can be rewritten, where the quantities that depend on temperature are explicitly noted,

S if T kin ( ) = h w if c B if n i T kin ( ) 1-e -h c w if k T kin ( ) (17)
In LTE, Tkin = T and given data at a reference temperature, Tref, the temperature at another T can be obtained by taking the ratio of Eq. ( 17) expressed at the two temperatures

S if T ( ) = S if T ref ( ) n i T ( ) 1-e -h c w if k T ( ) n i T ref ( ) 1-e -h c w if k T ref ( ) = S if T ref ( ) g i e -h c E i k T Q T ( ) 1-e -h c w if k T ( ) g i e -h c E i k T ref Q T ref ( ) 1-e -h c w if k T ref ( ) (18).
Cancelling like terms and assuming

E = E v + E r to obtain the product approximation, Q T ( ) = Q v T ( ) Q r T ( ) gives S if T ( ) = S if T ref ( ) Q v T ref ( ) Q r T ref ( ) Q v T ( ) Q r T ( ) e -h c E i k 1 T -1 T ref ae è ç ö ø ÷ 1-e -h c w if k T 1-e -h c w if k T ref ae è ç ç ö ø ÷ ÷ (19).
Writing Eq. ( 15) for the NLTE case and using a tilde to signify quantities in NLTE gives (20).

where T * signifies that the vibrational and rotational temperatures may be different. Note, it was assumed that for the lower level, i, the vibrational and rotational temperatures are the same. If it is not the case, the equation must be modified.

Allowing for the final state to have different vibrational and rotational temperatures gives Ef = Ev(T) + Er(Tkin) and the ratio of populations is (21).

The line intensity for the NLTE case is

(22).
Using Eqs. ( 17) and ( 22), the NLTE intensity at some kinetic temperature can be written in terms of the LTE intensity (scaled to Tkin from HITRAN) (23).

Noting that , where again we have assumed that in the lower state the vibrational and rotational temperatures are the same. Thus, the ratio of the population of state i in NLTE to LTE is the inverse ratio of the partition sums; is written as , giving

where is evaluated using Eq. ( 13) and it should be noted that the determination of allows for a different value of T for each vibrational level. The NLTE line intensity can also be written in terms of the HITRAN line intensity at the reference temperature (25).

Calculations of Qr and Qv

The calculations of the rotational and vibrational partition sums are presented for each molecule/isotopologue of this study. In all the calculations, h, k, and c, are the 2014 CODATA values [START_REF] Mohr | CODATA recommended values of the fundamental physical constants: 2014[END_REF]. Note, the Qvs determined below are the LTE values and were used to check the calculations by comparing with the TIPS2017 values [START_REF] Gamache | Total internal partition sums for 166 isotopologues of 51 molecules important in planetary atmospheres: Application to HITRAN2016 and beyond[END_REF]. The algorithm discussed in section 4 provides the NLTE Qv.

To apply Eq. ( 13 [START_REF] Herzberg | Molecular Spectra and Molecular Structure II. Infrared and Raman Spectra of Polyatomic Molecules[END_REF]), were extracted. It was observed that the rotational term values were not available for all states. This is due to the first-principles data not having complete attribution of the quantum numbers. To remedy this situation, the missing term values for the ground vibrational state were added using the smooth variation and pairing rules of Ma et al. [START_REF] Ma | Pair identity and smooth variation rules applicable for the spectroscopic parameters of H2O transitions involving high-J states[END_REF] for H2 16 O and HD 16 O. The rules show the energies of certain states labeled by the pseudo-quantum numbers Ka, Kc, which correspond to the projection of J along the inertia a and c axis, respectively, are equal or nearly equal.

Figure 1 shows the term values for J=34 in the ground vibrational state of H2 16 O versus Ka. In the figure the Ka+Kc odd states are given by the + symbol and the Ka+Kc even states are given by the  symbol. The MARVEL [START_REF] Tennyson | IUPAC critical evaluation of the rotational-vibrational spectra of water vapor, Part III: Energy levels and transition wavenumbers for H2 16 O[END_REF] term values are given by blue symbols, the BT2 [START_REF] Barber | A high-accuracy computed water line list[END_REF] 1).

Insert Figure 1 here Insert Figure 2 here Insert Figure 3 here

The term values in the Ev files go to values large enough to ensure convergence of Qv to 5000 K so the temperatures in Table 1 reflect the temperature for which Qr converges.

Using the rotational term values and vibrational energies discussed above, the rotational and vibrational internal partition sums were computed for each isotopologue of water vapor in this study by direct summation. The calculations ranged from 1 to Tmax in 1 K steps. Note, the Qv determined is the LTE value and was used to check the NLTE calculations.

CO2

The 12 isotopologues of CO2 considered here are shown in Table 1. The vibrational states of CO2 are labeled by 12l2n, where 1, 2, 3 express the number of quanta activated for each normal mode; l2 is the l value for the degenerate 2 normal mode and its overtones; the fifth integer is the nth component of the Fermi interacting 1 and 22 vibrational states, including their overtone and combination states [START_REF] Rothman | Infrared energy levels and intensities of carbon dioxide-II[END_REF][START_REF] Toth | Spectroscopic database of CO2 line parameters: 4300-7000 cm -1[END_REF]. The vibrational energies in cm -1 , Ev(12l2n), are from the first-principles calculations of Huang et al. [START_REF] Huang | Semi-empirical 12 C 16 O2 IR line lists for simulations up to 1500 K and 20,000 cm -1[END_REF][START_REF] Huang | Reliable InfraRed Line Lists for 13 CO2 Isotopologues up to E'=18,000 cm -1 and 1500K, with Line Shape Parameters[END_REF]. Because in the work of Huang et al. the rotational states were only calculated up to J=150, the rotational term values were recalculated here for the ground vibrational state of each isotopologue up to J=226 using the constants of [START_REF] Majcherova | High-sensitivity CW-cavity ringdown spectroscopy of 12 CO2 near 1.5μm[END_REF] with updates from the Tomsk coauthors. The maximum values of the rotational term values, Fr, of the ground vibrational states and the vibrational energies for each isotopologue and the number of vibrational states are provided in Table 1. However, care must be taken as the partition sums are to be determined for Qv and Qr separately and the product approximation used for the NLTE total internal partition sum. The four isotopologues of carbon dioxide, 12 C 16 O2, 13 C 16 O2, 12 C 18 O2, 13 C 18 O2, have equivalent oxygen atoms with zero spin nuclei (i.e. Bose particles). For these isotopologues, vibrational states with 2 = 0 only even J levels are allowed when 3 is even and only odd J levels exist when 3 is odd, however, when l2 > 0 all J levels exist. The total internal partition sum can be written as

Q T ( ) = ge -E vr /kT all states å = ge -E 00001,J ( )/kT + ge -E 01101,J ( )/kT + ge -E 10001,J ( )/kT + ge -E 00011,J ( )/kT + ge -E 10002,J ( )/kT + ge -E 11101,J ( )/kT + ge -E 01111,J ( )/kT + ge -E 02201,J ( )/kT + iii (26)
This expression can be rewritten by collecting terms with l2 = 0 and with l2 > 0

20 Q T ( ) = ge -E 00001,J ( )/kT + ge -E 10001,J ( )/kT + ge -E 00011,J ( )/kT + ge -E 10002,J ( )/kT + iii { } + ge -E 01101,J ( )/kT + ge -E 11101,J ( )/kT + ge -E 01111,J ( )/kT + ge -E 02201,J ( )/kT + iii { } = ge -E v,J ( )/kT allstates with ℓ 2 =0 å + ge -E v,J ( )/kT allstates with ℓ 2 >0 å (27)
Taking the ro-vibrational energy as a sum of the vibrational and rotational parts, Evr = Ev + Er and realizing that when l2 = 0 only even J levels are allowed when 3 is even and only odd J levels exist when 3 is odd and when l2 > 0 all J levels, J≧l2, are allowed, gives or

Q(T) = ge -E v kT
Q T ( ) = Q v ℓ 2 =0 n 3 even T ( )´Q r J even ( ) + Q v ℓ 2 =0 n 3 odd T ( )´Q r J odd ( ) + Q v ℓ 2 >0 T ( )´Q r all J ( ) (29).
This expression is labeled the triple product approximation (TPA). Equation ( 29) must be used for the symmetric Boze system CO2 isotopologues rather than the product approximation, Eq. ( 8). In order to make these calculations, the rotational terms values for the (01101) vibrational state were extracted from the work of Huang et al., fit to determine the rotational constants, and calculations were made up to J=226. Note, for the rotational terms values for (01101), the vibrational energy, Ev(01101), must be subtracted from the energies reported by Huang et al.

(see Eq. ( 29)). The comparison of the TPA and the PA to the TIPS determined from a direct sum in the work of [START_REF] Gamache | Total internal partition sums for 166 isotopologues of 51 molecules important in planetary atmospheres: Application to HITRAN2016 and beyond[END_REF] is presented in Fig. 4. Shown are the TIPS determined by direct summation (red dashed line), the TIPS determeined via the TPA (blue solid line with + symbols), and the TIPS determined via the PA (black dot-dashed line) versus temperature. The PA TIPS is in error by 4% at 296 K, 29% at 1000 K, and 46% at 5000 K.

Insert Figure 4 here

Using the term values and vibrational energies discussed above, the rotational and vibrational internal partition sums were computed for each isotopologue of carbon dioxide in this study by direct summation. While the vibrational partition sums are converged to 4000 K, the rotational partition sums are only converged to 3000 K, so Tmax was set to 3000 K (see Table 1).

The calculations ranged from 1 to 3000 K in 1 K steps. Note, the Qv determined is the LTE value and was used to check the NLTE calculations.

O3

The vibrational and rotational constants used in this work are from the "Ozone Spectroscopy and Molecular Properties" (S&MPO) databank [START_REF] Babikov | S&MPO -An information system for ozone spectroscopy on the WEB[END_REF]. updated in 2020. For rare isotopic species these data were obtained from the analyses of experimental spectra recorded in the Groupe de Spectrométrie Moléculaire et Atmosphérique (GMSA, Université de Reims) with the help of first-principles predictions carried out in Tomsk State University, Russia. The vibrational fundamentals, 1, 2, and 3 were used to generate vibrational energies, Ev1, 2, 3, with up to 9 quanta within the harmonic approximation of [START_REF] Herzberg | Molecular Spectra and Molecular Structure II. Infrared and Raman Spectra of Polyatomic Molecules[END_REF]. The rotational partition sums used the analytical formula of [START_REF] Watson | The asymptotic asymmetric-top rotational partition function[END_REF].

Calculations of Qv and Qr were made from 1 to 1000 K in 1 K steps.

N2O

The vibrational fundamentals and rotational constants of [START_REF] Toth | Frequencies of N2O in the 1100-to 1440-cm -1 region 1440-cm -1 region[END_REF] were used for the five isotopologues of nitrous oxide of this study: 14 N2 16 O, 14 N 15 N 16 O, 15 N 14 N 16 O, 14 N2 18 O, and 14 N2 17 O. The Ev1, 2, 2, 3 values were generated using Eq. (II-284) of [START_REF] Herzberg | Molecular Spectra and Molecular Structure II. Infrared and Raman Spectra of Polyatomic Molecules[END_REF] and

Ev data from HITRAN2020 [START_REF] Gordon | The HITRAN2020 molecular spectroscopic database[END_REF] to determine the constants. The Ev1, 2, 2, 3 data were then computed with up to 9 quanta of each fundamental frequency. The Qr data were generated using the McDowell formula for linear molecules [START_REF] Mcdowell | Rotational partition fuctions for linear molecules[END_REF].

Calculations were made from 1-5000 K in 1 K steps and are converged at all temperatures. were extracted from the work of [START_REF] Li | Rovibrational Line Lists for Nine Isotopologues of the CO Molecule in the X 1 Σ + Ground Electronic State[END_REF]. The calculations were made for each isotopologue in the electronic ground state with v41 and J150. For all isotopologues, the Ev values are complete to 63 802 cm -1 or greater and the rotational term values (on the ground vibrational state) are complete to 35 138 cm -1 or greater. Qv and Qr were computed at temperatures from 1 to 5000 K, in 1 K steps using the Ev values and term values discussed above and the resulting partition sums are converged at all temperatures for each isotopologue.

CH4

Four isotopologues were considered in this work: 12 CH4, 13 CH4, 12 CH3D, and 13 CH3D.

For the principal isotopologue, an effective Hamiltonian calculation of the ground vibrational state term values was done up to J=80, giving values complete to 30 424 cm -1 ensuring convergence at all temperatures. For 13 CH4, the rotational constants are from [START_REF] Dang-Nhu | Spectral intensities in the ν3 bands of 12 CH4 and 13 CH4[END_REF] and were used to compute Qr values using McDowell's spherical-top formula [START_REF] Mcdowell | Rotational partition functions for spherical-top molecules[END_REF]. The rotational constants for 12 CH3D are from [START_REF] Tarrago | Absorption of 12 CH3D at 6-10 μm: Triad ν3, ν5, ν6[END_REF], and those for 13 CH3D are from [START_REF] Ulenikov | Rotational Analysis of the Ground State and the Lowest Fundamentals ν3, ν5, and ν6 of 13 CH3D[END_REF], and were used to compute Qr values using McDowell's symmetric top formula [START_REF] Mcdowell | Rotational partition functions for symmetric-top molecules[END_REF].

Calculations were made for the Ev values (J=0) for 12 CH4, 13 CH4, 12 CH3D and 13 CH3D.

For semirigid molecules without large amplitude vibrations, like methane and its isotopologues, the normal mode representation turns out to be an excellent choice because it provides an adequate description of the nuclear motions. In this context, the normal mode kinetic energy operator was built in the frame of the Eckart-Watson (EW) formalism [START_REF] Watson | Simplification of the Molecular vibration-rotation Hamiltonian[END_REF]). The curvilinear potential energy surface (PES) reported by [START_REF] Nikitin | Rotational and vibrational energy levels of methane calculated from a new potential energy surface[END_REF] was transformed and

Taylor-expanded in mass-dependent normal mode coordinates q adapted to the Td and C3v point groups. Both the kinetic and potential parts have been expanded at order 14 in q and reduced at order 7 using the technique described in [START_REF] Rey | Understanding global infrared opacity and hot bands of greenhouse molecules with low vibrational modes from first-principles calculations: the case of CF4[END_REF]. For a full account of symmetry, the use of the irreducible tensor operator technique was used to build the EW Hamiltonian and rovibrational basis functions following Rey et al. [START_REF] Rey | Complete nuclear motion Hamiltonian in the irreducible normal mode tensor operator formalism for the methane molecule[END_REF][START_REF] Rey | Accurate first-principles calculations for 12 CH3D infrared spectra from isotopic and symmetry transformations[END_REF] while energy levels were computed by a variational method.

The J=0 vibrational problem was first solved using a direct product of harmonic oscillator functions where a "pruned" vibrational basis was considered by selecting a limited number of functions through the criterion

F n ( ) = k i v i i =1 M å £n (29).
with M=4 for Td species and M=6 for C3V species. Here, vi=0,…,n and ki are weighting coefficients. For this work, n=15 and 𝑘 𝑖 = {1; 1; 1.1; 1} for 12 CH4 and 13 CH4 and n=14 and 𝑘 𝑖 = {1.2; 1.1; 1; 1.2,1,1} for 12 CH3D and 13 CH3D, resulting in 36 889, 32 903, 69 751, 101 912 and 105 888 functions for the symmetry blocks and 70 552, 65 855 and 136 374 functions for the symmetry blocks A 1 , A 2 , E { } , respectively. In contrast to the effective polyad models, variational calculations require the diagonalization of very large matrices for J > 0, even for computing pure rotational energy levels. To this end, a set of vibrational reduced eigenfunctions obtained from the projection technique [START_REF] Rey | Understanding global infrared opacity and hot bands of greenhouse molecules with low vibrational modes from first-principles calculations: the case of CF4[END_REF] was introduced to drastically reduce dimensionality of the full problem, making thus calculations feasible even for high J values. These vibrational energies were used to compute the Qv values and were also put into files to compute from the code of section 4.

A 1 , A 2 , E, F 1 , F 2 { }
The number of vibrational levels for all isotopologues of methane considered here is large. Because of the density of states and the polyad structure of the vibrational levels, Table 2 should not be used to determine the convergence of Qv. Convergence plots were made for temperatures from 200 to 2500 K. Figure 5 shows the calculation of Qv as a function of the vibrational energy for a summation up to 10 000 cm -1 . Shown are Qv values (blue lines) for the temperatures of 500, 750, 1000, 1250, and 1500 K versus the vibrational energy in cm -1 . The red line indicates the final Qv value after summing over all the energies and hence gives an indication of convergence. Also shown on the 500 K curve are the values of the vibrational energies (red dots) showing the polyad structure and the density of states. Looking at the curves, it is clear that summing energies up to 10 000 cm -1 gives Qvs that are only converged up to 1000 K. The erratic structure in the curves at low energy and for higher temperatures is due to the polyad structure and the density of states. Table 3 should be used to determine the convergence of Qv for the methane isotopologues.

The Qv and Qr values were computed from 1 to 2500 k in 1 K steps.

Insert Figure 5 Assuming that the vibrational levels are the same for the 2 P 1 2 and 2 P 3 2 electronic states, as

shown by [START_REF] Amiot | Infrared study of the X2Π ν = 0, 1, 2 levels of 14[END_REF], the energy can be written as (32) which is of the form of Eq. ( 8).

E a 2 P 1 2 ( ) = E v + E r 2 P 1 2 ( ) and E b 2 P 3 2 ( ) = E v + E
For the 14 N 16 O, 15 N 16 O, and 14 N 18 O isotopologues of nitric oxide, the rotational term values were calculated for the 2 1/2e, 2 1/2f, 2 3/2e, and 2 3/2f, states using the molecular constants of [START_REF] Amiot | Infrared study of the X2Π ν = 0, 1, 2 levels of 14[END_REF] and up to J=149.5, F=150.5 and are complete to 35 391, 34 233, and 33 693 cm -1 , respectively. These term values were used to compute Qr by direct sum.

The fundamental frequency for each isotopologue ,i, were taken from [START_REF] Amiot | Infrared study of the X2Π ν = 0, 1, 2 levels of 14[END_REF] and [START_REF] Meerts | A theoretical reinvestigation of the rotational and hyperfine lambda doubling spectra of diatomic molecules with a 2Π state: the spectrum of no[END_REF] were used to determine Qv via the HA and to generated the Ev files.

However, Aimot et al. report a value for v=2 for 14 N 16 O, which replaced the HA value in the Ev file. Both partition sums, Qv and Qr, were computed from 1 to 5000 k in 1 K steps for each isotopologue.

NO2

Calculations of Qv and Qr were made for the 14 N 16 O2 and 15 N 16 O2 isotopologues of nitrogen dioxide. The rotational and vibrational constants were taken from [START_REF] Perrin | The far infrared spectrum of 14 N 16 O2 Electron spin-rotation and hyperfine Fermi contact resonances in the ground state[END_REF] for 14 N 16 O2 and from [START_REF] Perrin | Detection of atmospheric 15 NO2 in the ν3 spectral region (6.3 μm)[END_REF] for 15 N 16 O2. The analytical formula of [START_REF] Watson | The asymptotic asymmetric-top rotational partition function[END_REF] was used to determine Qr. Note, the state independent statistical factor is determined by considering that both isotopologues are Bose particles, so half of the states are missing.

However, lambda doubling adds a factor of 2. These factors must be multiplied by the nuclear spin factor of the N atom giving: ½ ⋅2⋅(2⋅1+1) = 3 for 14 N 16 O2 and ½ ⋅2⋅(2⋅½+1) = 2 for 15 N 16 O2.

The vibrational fundamentals, 1, 2, and 3 were used to generate vibrational energies, Ev1, 2, 3, with up to 9 quanta within the harmonic approximation of [START_REF] Herzberg | Molecular Spectra and Molecular Structure II. Infrared and Raman Spectra of Polyatomic Molecules[END_REF]. The final partition functions were calculated from 1 to 5000 K in 1 K steps.

3.9 OH Three isotopologues of the hydroxyl radical were considered in this work: 16 OH, 18 OH, and 16 OD. The rotational term values of the principal isotopologue were taken from [START_REF] Brooke | Line strengths of rovibrational and rotational transitions in the ground state of OH[END_REF] and are complete to J=58.5 E= 46 065.460 3 cm -1 . For the 18 OH and 16 OD isotopologues, the rotational term values for the 2 1/2, 2 3/2, and 2 1/2, states were calculated using the formalism of [START_REF] Beaudet | Mirowave Spectra of Molecules of Astrophysical Interest. XII. Hydroxyl Radical[END_REF], which includes the fine structure interaction and lambda doubling. This procedure generated energies complete to 42 000 cm -1 for 18 OH, and 39 000 cm -1 for 16 OD. The rotational partition sums were calculated by direct summation over the rotational term values. The vibrational fundamentals were taken from [START_REF] Beaudet | Mirowave Spectra of Molecules of Astrophysical Interest. XII. Hydroxyl Radical[END_REF] and used to generate harmonic approximation vibrational energies, Ev, with up to v=13

giving maximum energies of 33 086, 31 586, 32 643, for 16 OH, 18 OH, and 16 OD, respectively.

The final partition functions were calculated from 1 to 5000 K in 1 K steps.

Recall of the NLTE partition sums

Equations ( 24) and ( 25) show that to correctly determine the intensity of a spectral line in NLTE conditions both Qr(T) and are needed. The Qr(T) values can easily be precalculated as a function of temperature, stored in the code for fast recall. Computing is more difficult in that a file of vibrational temperatures is needed for each layer of the A note about the T files, most of the calculations have a Tmax equal to 5000 K. To ensure convergence of , the vibrational states were determined to high energies. For applications at lower temperatures the files can be truncated appropriately.

3. 5

 5 CO Nine isotopologues of CO are studied here: 12 C 16 O, 13 C 16 O, 12 C 18 O, 12 C 17 O, 13 C 18 O, 13 C 17 O, 14 C 16 O, 14 C 18 O, and 14 C 17 O. The rotational term values (v=0) and the Ev values (J=0)

  rise to e and f energy states. The states are labeled by J and F=J±1,0 and e or f. The total internal partition sum can be written as

  odd states are given by the + symbol and the Ka+Kc even states are given by the  symbol. The MARVEL term values are given by blue symbols, the BT2 term values are given by red symbols, and the term values obtained using the rules ofMa et al. are given by black symbols. 2 The percent difference between the H2 16 O and H2 18 O term values versus an energy ordered index, J(J+1)+Ka-Kc+1 3 The rotational internal partition sums for H2 16 O versus term value for the following six temperatures: 500, 1000, 2000, 3000, 4000, and 5000 K 4 The comparison of the TPA and the PA to the TIPS determined by direct sum. Shown are the TIPS determined by direct summation (red dashed line), the TIPS determeined via the TPA (blue solid line with + symbols), and the TIPS determined via the PA (black dotdashed line) versus temperature 5 The calculation of Qv for 12 CH4 as a function of the vibrational energy for a summation up to 10 000 cm -1 . Shown are Qv values (blue lines) for the temperatures of 500, 750, 1000, 1250, and 1500 K versus the vibrational energy in cm -1 . The red line indicates the final Qv value after summing over all the energies and hence gives an indication of convergence. Also shown on the 500 K curve are the values of the vibrational energies (red dots) showing the polyad structure and the density of states.

  Figure 1 812

  ) to determine , files of the vibrational term values, Ev, for each isotopologue of each molecule of this study are needed. These Ev files have been prepared and are available with the FORTRAN code discussed in section 4. The number of vibrational states in each file is given in Table 1; listed are the molecular isotopologue, the maximum rotational Nine isotopologues of water vapor were considered in this work: H2 16 O, H2 18 O, H2 17 O, HD 16 O, HD 18 O, HD 17 O, D2 16 O, D2 18 O, and D2 17 O. The initial sets of term values are from[START_REF] Barber | A high-accuracy computed water line list[END_REF][START_REF] Tennyson | IUPAC critical evaluation of the rotational-vibrational spectra of water vapor, Part III: Energy levels and transition wavenumbers for H2 16 O[END_REF] for H2 16 O, from[START_REF] Tennyson | IUPAC critical evaluation of the rotational-vibrational spectra of water vapor. Part I-Energy levels and transition wavenumbers for H2 17 O and H2 18[END_REF] for H2 18 O and H2 17 O, from Refs.[START_REF] Tennyson | IUPAC critical evaluation of the rotational-vibrational spectra of water vapor. Part II: Energy levels and transition wavenumbers for HD 16 O, HD 17 O, and HD 18 O[END_REF][START_REF] Voronin | A high accuracy computed line list for the HDO molecule[END_REF] for HD16 O, from[START_REF] Tennyson | IUPAC critical evaluation of the rotational-vibrational spectra of water vapor. Part II: Energy levels and transition wavenumbers for HD 16 O, HD 17 O, and HD 18 O[END_REF] for HD18 O and HD 17 O, and from[START_REF] Simkó | Recommended Ideal-Gas Thermochemical Functions for Heavy Water and its Substituent Isotopologues[END_REF][START_REF] Tennyson | IUPAC critical evaluation of the rotational-vibrational spectra of water vapor. Part IV. Energy levels and transition wavenumbers for D2 16 O, D2 17 O, and D2 18 O[END_REF] for D216 O, D2 18 O, and D2 17 O. These term values were taken and the rotational states of the ground vibrational state were extracted and the vibrational states with J=0 states, designated Ev (note that these are called G(1, 2, 3) in Herzberg's notation

	Insert Table 1 here
	Insert Table 2 here
	3.1 H2O

term value, the Qr convergence (RC), the maximum vibrational term value, the Qv convergence (VC), the number of vibrational term values, and the maximum temperature for which Qv Qr has converged. The convergence of the partition sums (Qr and Qv) is determined by the percent difference between Q determined by summing over all the term values and Q determined by summing over three-fourths of the term values. For some molecular isotopologues the number of term values is quite large. Users may want to truncate the list of Ev values; however, the convergence of the vibrational partition sum is dependent on temperature. Table

2

presents a rough guide relating the maximum term value to the maximum temperature for which Qv is converged. Table

3

gives similar information as Table

2

for methane (see discussion in section 4). These cutoffs are discussed below when appropriate for each molecule/isotopologue combination.

  term values are given by red symbols, and the term values obtained using the rules ofMa et al. are given by black symbols. For the lesser abundant non-deuterated and singly deuterated isotopologues, H2 18 O, H2 17 O, HD 18 O, and HD 17 O, there are not enough MARVEL or labeled ab initio term values to apply the rules ofMa et al. The term value differences between the lesser abundant species and H2 16 O or HD 16 O were considered. Figure2shows the percent difference between the H2 16 O and H218 O term values versus an energy ordered index, J(J+1)+Ka-Kc+1. The singly deuterated isotopologues show similar features. Since there is not a great difference in term values, when term values were missing in the files for H2 18 O, H2 17 O, HD 18 O, and HD 17 O they were added using the H2 16 O or HD 16 O data. This procedure allows the rotational partition sums to be converged to 5000 K and 4000 K for the H2O and HDO isotopologues, respectively. The convergence of Qr for H216 O is shown in Fig.3. Shown are the rotational internal partition sums versus term value for the following six temperatures: 500, 1000, 2000, 3000, 4000, and 5000 K. At 5000 K, Qr is converged to better than 0.4%. Note that for the doubly deuterated isotopologues, the available terms values allow convergence of Qr to 3000 K for D2 16 O, D2 18 O, and D2 17 O (see Table

Table 2

 2 where the maximum temperatures are determined for each isotopologue such that Qv and Qr are converged, see Table1for details. These data are presented in Table1. The calculations of Qv and Qr are described for each molecule/isotopologue combination. A FORTRAN code, NLTE_TIPS_2021.for, to compute the NLTE , the recall of Qr, and the needed vibrational energy files are provided. Truncating the list of vibrational energies is possible for some lower temperature applications. To aid users, tables of the maximum vibrational energies needed for convergence of as a function of temperature are presented for the molecules in general (Table2) and separately for methane (Table3). The code uses the MIT license and it and the needed vibrational energy files and can be downloaded at http://faculty.uml.edu/Robert_Gamache , at www.HITRAN.org, or from zenodo.org(TBD).Tables 1 Isotopologue, maximum rotational term value, maximum vibrational term value, number of vibrational states, maximum temperature of reported partition sums.2 Term value versus maximum temperature for convergence of vibrational partition sum.Note that this table excludes methane.3 Term value cutoff needed to insure convergence of vibrational partition sums for methane isotopologues. Table 1 Isotopologue, maximum rotational term value, maximum vibrational term value, 770 number of vibrational states, maximum temperature of reported partition sums. 771 The term values for J=34 in the ground vibrational state of H2 16 O versus Ka. The Ka+Kc

	provides a guide to the

a convergence of Qrot in percent; see text. b convergence of Qvib in percent; see text.c analytical formula used

Acknowledgements B. Vispoel is pleased to acknowledge le Fonds de la Recherche Scientifique -FNRS for post-doctoral financial support. The GSMA Reims group acknowledges the support from CNRS (France) in the frame of the International Research Project SAMIA with IAO-Tomsk. V.T. acknowledges the support from the Russian Science Foundation (project RNF-19-12-00171). JT and SY thank ERC's Advanced Investigator Awards 267219 (2011-16) and 883830 225 (2020-25), as well as STFC as part of the UCL consolidated grant ST/R000476/1. The work in Tomsk was supported by Ministry of Science and Higher Education of the Russian Federation. The work in Budapest was supported by NKFIH under grant no. K138233.

atmosphere. To accomplish this task, a FORTRAN code was written to store Qr(T) and Qv (T) data and recall the values as a function of temperature using Lagrange interpolation and using files of vibrational wavenumbers and temperatures to calculate for the atmospheric layer in question. The code, NLTE_TIPS_2021.FOR, stores Qr(T) and Qv(T) from 1-20 K in 1 K steps and then the values at every 10 K to the maximum temperature of the computations for a given molecule. The user is prompted for the molecule, isotopologue, and kinetic temperature to determine Qr(T) and Qv(T). The code uses a 4-point Lagrange interpolation scheme, similar to earlier versions of TIPS [START_REF] Fischer | Total Internal Partition Sums for Molecular Species on the 2000 Edition of the HITRAN Database[END_REF][START_REF] Laraia | Total internal partition sums to support planetary remote sensing[END_REF]. The value of Qv(T) is only computed for reference, it is not used in Eqs. ( 24) or (25). The computation for is made by summing over the vibrational energies with the appropriate T values for the atmospheric layer in question. There are files for each isotopologue of this study provided with the code.

These T files are named by molecule_isotopologue#_layer.dat; for example, for