Trajectory-based Simulation of Far-infrared Collision-induced Absorption Profiles of CH4–N2 for Modeling Titan’s Atmosphere

We report the results of the trajectory-based simulation of far-infrared collision-induced absorption (CIA) due to CH4–N2 pairs at temperatures between 70 and 400 K. Our analysis utilizes recently calculated high-level potential energy and induced dipole surfaces. Treating collision partners as rigid rotors, the time evolution of interaction-induced dipole is accumulated over a vast ensemble of classical trajectories and subsequently transformed into a CIA spectrum via Fourier transform. In our calculations, both bound and unbound states are properly accounted for, and the rigorous theory of lower-order spectral moments is addressed to check the accuracy of simulated profiles. Classically derived trajectory-based profiles are subject to two approximate desymmetrization procedures so that resulting profiles conform to the quantum principle of detailed balance. The simulated profiles are compared to laboratory measurements and employed for modeling Titan’s spectra in the 50–500 cm−1 range. Based on the desymmetrized simulated profiles, a new semiempirical model for CH4–N2 CIA is proposed for modeling Titan’s infrared spectra. Synthetic spectra derived using this model yield an excellent agreement with the data recorded by the Composite Infrared Spectrometer aboard the Cassini spacecraft at low and high emission angles.


Introduction
Nitrogen and methane were well established as the primary atmospheric constituents of Titan's atmosphere after the Voyager 1 encounter (Hunten et al. 1984). With a concentration of roughly 95%, nitrogen is the most abundant species in the atmosphere of Titan. Methane, while being less abundant, is crucial to the maintenance of the thick nitrogen atmosphere; the atmosphere would gradually condense in the absence of warming resulting from the far-infrared CH 4 -N 2 collisioninduced absorption (CIA) and hydrocarbon haze (Lorenz et al. 1997). The photochemistry of nitrogen and methane leads to the formation of complex hydrocarbons and nitriles in the atmosphere of Titan. Figure 1 highlights the magnitude and wavenumber dependence of the major CIA contributions of relevant collision pairs in Titan's atmosphere. The maximum of spectral radiance from relatively cold Titan's surface falls in the vicinity of 200 cm −1 . The most abundant atmospheric species like N 2 , CH 4 , and H 2 have no permanent dipole and are therefore nearly transparent in the relevant far-infrared spectral range. As a result, the opacity of Titan's atmosphere is largely determined by continuum absorption dominated by the collision-induced component that belongs to molecular pairs formed from the main atmospheric constituents: N 2 -N 2 , CH 4 -N 2 , CH 4 -CH 4 , N 2 -H 2 , H 2 -H 2 , and CH 4 -H 2 . Due to the large abundance of N 2 and relatively small rotational constant, B ≈ 2 cm −1 , the N 2 -N 2 CIA dominates the tropospheric opacity for wavenumbers shortward of ∼150 cm −1 , whereas the CH 4 -N 2 CIA is the prevailing source of opacity in the 150-450 cm −1 region (Samuelson et al. 1997). Several previous analyses (Tomasko et al. 2008;de Kok et al. 2010;Anderson & Samuelson 2011) concluded that the opacity due to CH 4 -N 2 could not be correctly reproduced using the model of Borysow & Tang (1993). Courtin et al. (1995) and Samuelson et al. (1997) suggested that missing opacity could be attributed to a large supersaturation of methane in low-latitude regions of the upper troposphere. However, measurements performed using the Gas Chromatograph Mass Spectrometer on the Huygens probe indicated that there is no methane supersaturation in the lower atmosphere at 10°S latitude (Niemann et al. 2005). Later studies (Tomasko et al. 2008;de Kok et al. 2010) suggested that Borysow & Tang's (1993) CIA data is possibly systematically in error at Titan's temperatures. Both studies derived that a multiplicative factor of approximately 1.5 needs to be applied to Borysow & Tang's (1993) data to obtain satisfactory fits of the Cassini's Composite Infrared Spectrometer (CIRS) observations in the 150-450 cm −1 region. Anderson & Samuelson (2011) proposed a heuristic correction factor of the same magnitude. Bézard & Vinatier (2020) reported that the wavenumber-dependent correction of Anderson & Samuelson (2011) produces too much absorption between 180 and 300 cm −1 . Instead, Bézard & Vinatier (2020) chose to use a wavenumber-independent multiplicative factor, the value of which was refined to 1.52.
The most accurate laboratory measurements of CIA in CH 4 -N 2 gaseous mixtures were performed up to 550 cm −1 at 162 K and up to 670 cm −1 at 195 and 297 K by Birnbaum et al. (1993).
An earlier effort by Dagg et al. (1986) provides experimentally measured spectra at 126, 149, 179, and 212 K over the frequency region of 40-400 cm −1 . These measurements, however, are affected by a substantial noise, most notably at the edge of the covered frequency range at 350-400 cm −1 . The temperature of Titan's atmosphere varies from 70 to 94 K, close to the surface, to 120 K in the upper atmosphere. Thus, the temperatures at which the absorption spectra are required for modeling Titan's atmosphere are much lower than those at which laboratory experiments were carried out.
Despite significant efforts undertaken in the past to characterize the CIA sources of opacity for Titan, the required data remain either incomplete or not accurate enough to be used in the modeling of the atmospheric radiation. A series of papers by Borysow and Frommhold (1986a, 1986b, 1986c, 1987 and a paper by Borysow & Tang (1993) provided a database of calculated opacities based on the effective isotropic potentials combined with adjusted short-range induced dipoles. The intrinsic parameters of the induced dipole were derived through fitting to experimental spectra. This approach inherently depends on the quality of experimental data since the possible systematic errors in the fitted data would be introduced in the model parameters. Furthermore, the validity of such profiles modeled far outside the temperature range for which the intrinsic parameters have been fitted is under question. In the case of CH 4 -N 2 , Borysow & Tang (1993) opted to perform fitting to measurements by Birnbaum et al. (1993).
Recently, sophisticated first-principles approaches have been developed and applied to simulate far-infrared absorption in N 2 -N 2 using quantum calculation with proper account for intermolecular anisotropy (Karman et al. 2015), molecular dynamics (Hartmann et al. 2011), and classical trajectory-based simulation (Chistikov et al. 2019). Encouraged by the success of the trajectory-based approach in the case of N 2 -N 2 (Chistikov et al. 2019;Gordon et al. 2022) and CO 2-Ar , we extended it toward linear molecules interacting with spherical top molecules. Our consideration relies on the ab initio potential energy (PES) and induced dipole surfaces (IDS) for CH 4 -N 2 calculated at CCSD(T)-F12/ aug-cc-pVTZ and CCSD(T)/aug-cc-pVTZ levels of theory, respectively, reported in Finenko et al. (2021). Section 2 presents the details of our trajectory-based simulation of CIA in CH 4 -N 2 molecular pairs, including transitions characteristic to unbound (free) and quasi-bound molecular pairs as well as to true bound dimers. Section 3 shows that the semiempirical model based on the ab initio simulated CH 4 -N 2 CIA profiles allowed us to reproduce Titan's spectra between 50 and 650 cm −1 recorded by the CIRS aboard the Cassini spacecraft with unprecedented accuracy. Potential implications of the simulated profiles for greenhouse warming of early Earth's atmosphere are discussed in Section 4.

Summary
In first-order perturbation theory, the binary absorption coefficient α(ν, T) in a mixture of gases can be represented as (Frommhold 2006) is the absorption coefficient, ρ 1 and ρ 2 are the number densities of gas constituents, 1 is the reciprocal temperature in units of energy, k B , h, c, and N L are the Boltzmann and Planck constants, speed of light, and Loschmidt number, respectively, V designates the gas volume, and ν is the wavenumber in reciprocal centimeters. The spectral function G (ν, T) is given by the Fourier transform of the time-correlation function (TCF) C(t, T): which is defined as the statistically averaged product of the time-shifted dipoles relevant to the compound molecular pair: where ε 0 is the vacuum permittivity, and angular brackets denote phase-space averaging. The spectral function conforms to the so-called detailed balance condition (Frommhold 2006): , , 5 hc causing the striking asymmetry of the observed spectral profiles. If exact quantum formalism is employed to simulate the spectral function, the detailed balance condition (5) is satisfied naturally. However, if one chooses to address the approach utilizing classical mechanics, the resulting spectral function is symmetric: cl cl in violation of Equation (5). This is caused by the time-reversal invariance of classical dynamics, which induces the symmetry of the corresponding TCF. An appealing approach to approximately correct this major defect is through the use of  (Borysow & Frommhold 1987), and N 2 -H 2 (Borysow & Frommhold 1986a) are shown separately. The contribution from CH 4 -N 2 CIA is obtained employing a semiclassical model as described in Section 3 (see Equation (23)). The calculations were performed utilizing the atmosphere model of Bézard & Vinatier (2020).
the so-called desymmetrization procedures (Borysow et al. 1985;Frommhold 2006). These are ad hoc prescriptions for constructing spectral functions that obey a detailed balance condition from the classical counterpart. A variety of such prescriptions were previously proposed to reproduce profiles rendered from the quantum approach for model problems (Schofield 1960;Egelstaff 1962;Berens et al. 1981;Bader & Berne 1994). Let us formulate several known prescriptions in terms of the quantum correction factor ( Schofield's procedure (Schofield 1960) was demonstrated to provide a satisfactory semiclassical correction to classically simulated CIA rototranslational bands for multiple collisional complexes (Hartmann et al. 2011;Bussery-Honvault & Hartmann 2014;Chistikov et al. 2019;Odintsova et al. 2021). The more sophisticated quantum correction factor due to Egelstaff (1962) was shown to provide the most accurate description of spectra in comparison to other corrections, at least in the case of freely rotating linear molecules (Borysow & Moraldi 1994). The extension of Egelstaff's procedure to the case of interacting polyatomics is, therefore, somewhat more grounded than any other correction suggested in the literature. Frommhold (2006) suggested a slight variant of Egelstaff's procedure enhanced by the constant rescaling factor. This variant, given by provides the closest approximation to the quantal spectral profiles for the purely translational spectra (Frommhold 2006). As a final result of our consideration, we adopt the semiclassical binary absorption coefficient where quantum correction factor ( ) n T , is represented either by Schofield's (8) or Frommhold's procedure (9).

Trajectory-based Approach
Invoking the Hamiltonian formalism, let us introduce the generalized coordinates describing the position and orientation of moieties relative to the space-fixed frame placed at the center of mass of the bimolecular complex. Let describe an active rotation of the CH 4 and N 2 moiety, respectively, from an initial position, in which a reference frame fixed on the moiety is parallel to the space-fixed frame, to its present position.
In the space-fixed frame of reference, the Hamiltonian of the CH 4 -N 2 system breaks up into where U is the potential energy function, and the kinetic energies of translational and rotational motions are given by Here μ, I 1 , and I 2 denote the reduced mass of the molecular pair and moments of inertia of CH 4 and N 2 , respectively. The momentum conjugated to the generalized variable ξ is designated via p ξ . Hereafter, the generalized coordinates and conjugated momenta are assembled into vectors q and p, respectively. We consider the time evolution of the bimolecular complex in terms of the solution of the Hamilton equations arising from the Hamiltonian (11): The potential energy and induced dipole surfaces are commonly expanded in the body-fixed frame in order to reduce the number of angular degrees of freedom. Here we utilized the PES and IDS computed at a high level of ab initio theory in Finenko et al. (2021). Those were constructed in the frame in which the orientation of the CH 4 is fixed in its initial orientation. The schematic representation of the frame is given in Figure 2. The origin is placed on the carbon atom of the CH 4 monomer. The orientation of the vector R is described by angles ( ) F Q , BF BF . Frame 2, parallel to Frame 1 fixed to the CH 4 molecule, has its origin in the center of mass of N 2 . The orientation of the N 2 molecule with respect to Frame 2 is described by angles ( ) h c , BF BF . The main computational impediment that we face when solving dynamical Equations (13) is to find the derivatives of the potential energy U with respect to the space-fixed angles Φ, Θ, j, ϑ, ψ, η, χ. The numerical scheme of calculating the desired derivatives is described in Appendix A.
From a formal point of view, one can opt for direct simulation of the TCF using the classical trajectories method for the purpose of then transforming it to a spectral function according to Equation (3). The Monte Carlo approach provides means for the simulation of a TCF since the high dimensionality of a phase space forbids practical use of deterministic numerical integration schemes. The convergence issues at large frequencies ν, however, hinder the application of the direct approach for free and quasi-bound states. An alternative integral representation of the spectral function is often adapted in the classical trajectories modeling of these pair states (van Kranendonk & Gass 1973;Oparin et al. 2017;Chistikov et al. 2019). In the seminal work by van Kranendonk & Gass (1973), the transformation into an integral expression over complete trajectories is suggested, alleviating the convergence difficulties. Chistikov et al. (2019) suggested a similar integral representation suitable for a canonical set of Hamilton variables. Adopting the latter integral representation for free (F) and quasi-bound (Q) states, the corresponding classical spectral function Here, Q(T) represents the partition sum 1 5 H An instantaneous value of the dipole moment vector of the collisional complex is designated with μ(t), and vectors q å and p å are obtained upon exclusion of the coordinate R and conjugated momentum p R from the vectors q and p, respectively.
Integral representation (14) should be understood in the following way. First, some fixed value of intermolecular distance R max is chosen, and a hyperplane in the phase space corresponding to the chosen value is drawn. All trajectories passing through the selected hyperplane are properly accounted for in the representation (14). There are infinitely many far flyby trajectories for which the intermolecular distance R always exceeds R max that are naturally excluded from the expression (14). The larger the value of R max is chosen, the more negligible the contribution to the spectral function gets from the excluded trajectories since the induced dipole vanishes at extreme separation. As discussed in Chistikov et al. (2021), there is a family of quasiperiodic trajectories that, while having positive energy, do not result in the spontaneous breakdown of the molecular complex into monomers within two-body consideration. The time evolutions of dipole moment along those peculiar trajectories contribute to the spectral function + G cl F Q . The sampling of quasiperiodic trajectories through representation (14) can be achieved only when the value of R max is somewhat close to the potential well. Thus it becomes clear that no value of R max allows for taking into account contributions from all significant trajectories. The spectral signatures from the quasiperiodic trajectories are restricted to comparatively low frequencies and hence are adequately simulated through the calculation of the TCF. Selecting a large value of R max for the representation (14) permits simulation of the spectral profile, excluding the quasiperiodic contributions. Therefore we opt to combine a TCF-based spectral function at lower frequencies with the one based on the representation (14) in the spectral wing using a basic stitching procedure (see Section 2.3). Contributions from all the mentioned pair states are sufficiently accounted for in the combined spectral function, which is confirmed by the agreement between spectral function moments and their phase-space counterparts described in Section 2.4. The calculation of the bound states' spectral function is addressed through the calculation of the TCF supplemented by the proper averaging over the true bound domain in the phase space as described below.
To numerically evaluate the multivariate integral (14), we invoke the importance sampling approach with proposal distribution over the domain of unbound states. The sampling algorithm is outlined in Appendix B. Given a supply of N random variables  , the Monte Carlo estimate is given by (For a more detailed discussion, see Chistikov et al. 2019.) Here, ( ·) m t; is the value of the dipole moment at a time t during classical evolution of the trajectory starting initially at the given phase point.
The Monte Carlo simulation of the TCF is addressed by calculating an ensemble of classical trajectories emerging from Boltzmann-distributed initial points. The sampling from the Boltzmann distribution is performed in the course of a two-step rejection-based procedure summarized in Appendix B. Given N phase space points { } q p , k k supplied by the rejection algorithm, the TCF in the domain min max is the partial partition sum defined by Averaging of the products of dipole moments in Equation (17) (17) can be reweighted to produce correlation function estimates at various temperatures. This computational technique is based on the premise that an individual classical trajectory is a microcanonical object that does not depend on temperature. The staged multitemperature Monte Carlo algorithm is described in detail in Section IV.C of Chistikov et al. (2021) and used here without modifications for the temperature range 70-400 K with a 10 K step.
The correlation functions along individual trajectories were sampled at fixed intervals of 4.8 fs until = t 1.58 max ns is reached. The differential equation solver CVODE from the SUNDIALS suite (Hindmarsh et al. 2005) was employed to propagate equations of motions (13). The obtained TCFs approach the positive limit as time goes to infinity. Because of this, one should take precautions before transforming TCFs into spectral functions according to Equation (3). Direct application of standard numerical procedures implementing a Fourier transform would cause artifacts to appear. The utilized processing technique allowing us to obtain proper spectral functions from the corresponding TCFs is explained in Chistikov et al. (2021). Obtained spectral functions are subjected to smoothing using local quadratic fits with a variable-width window to remove the residual Monte Carlo noise. Subsequently, absorption coefficients are produced according to Equation (10) with properly chosen quantum correction factors.

Stitching Procedure
As mentioned above, the wings of the spectral profiles derived from the TCFs corresponding to unbound states exhibit broad unphysical structure, especially for low temperatures, demonstrated in Figure 3. This structure fades slowly with an increasing number of trajectories and thus can be considered presumably as Monte Carlo noise. To mitigate this effect, the wings of the spectral profiles derived through the formalism (14)-(16) are stitched continuously to TCF-based profiles at ν = 200 cm −1 . The stitching is performed independently for the spectral profiles obtained using Schofield's and Frommhold's quantum correction factors. The thus-obtained arrays of spectral profiles for both quantum corrections are shown in Figure 4. The overall profile can be roughly decomposed into three spectral contributions having different dipole origins (Birnbaum et al. 1993;Borysow & Tang 1993). First, absorption due to the dipole in the nitrogen molecule induced by the electric multipole field of methane manifests itself in the form of the pronounced peak. Second, spectral components due to the electric fields of nitrogen polarizing the methane molecule result in the formation of the spectral shoulder. Lastly, there are double-transition terms, usually much weaker, originating from the interaction of the gradient of the multipolar electric field of one molecule with the dipolemultipole polarizability tensor of another molecule. The intensity of the spectral shoulder greatly depends on the choice of the quantum correction factor, especially at low temperatures. Due to the rescaling factor introduced by Frommhold (2006) to modify Egelstaff's desymmetrization, the conspicuous rotational peak caused by the collision-induced rotational transitions in the N 2 molecule has an almost identical intensity for Schofield's and Frommhold's quantum corrections. It should be noted that the bound states' contribution can be seen in Figure 4 as the difference between spectral profiles corresponding to both bound and unbound states and unbound states only. It is marginal for the lowest temperatures, amounting to 4% of integral intensity at 70 K, and becoming effectively zero at temperatures higher than 120 K.

Spectral Moments as Convergence Control Parameters
The convergence of the obtained spectral profiles (17) with respect to the number of trajectories is controlled through the calculation of the lower-order spectral moments. The usefulness of spectral moments in this context arises from the fact that they can be derived from two significantly different prerequisites. On the one hand, the classical spectral moments can be expressed as the weighted integrals of the binary absorption coefficient: On the other hand, the lower-order spectral moments can be obtained through phase-space averaging according to  (5) and (6) (the issue is expounded upon in Finenko et al. 2021 andChistikov et al. 2021). To compare the moments issued from either the observations or the calculations, one must first introduce the respective correction. We chose to compare the moments that are derived utilizing classical detailed balance, Equation (6).
Consequently, prior to computing the moments from the laboratory-measured spectra using Equation (20) a symmetrization procedure should be applied to produce the spectral shape conforming to the classical detailed balance. The symmetrization procedure is procured as an inverse to the desymmetrization procedure. For its simplicity, we opted to invert Schofield's procedure so that the measured binary absorption coefficient is multiplied by the factor .
2 2 hc inv Sch 2 In Figure 5, classical spectral moments obtained through phase-space integration (Equation (21)) are compared to the values derived from the trajectory-based spectral profiles via Equation (20) as well as from the laboratory-measured spectra. Perfect agreement between the spectral moments obtained via frequency-domain and phase-space integration validates the calculated spectral profiles. Figures 6 and 7 compare the calculated CIA CH 4 -N 2 spectral profiles with the experimental measurements by Dagg et al. (1986) and Birnbaum et al. (1993), respectively. It can be seen that in the spectral shoulder and far wing, between 150 and 600 cm −1 , the spectral profiles desymmetrized using Schofield's and Frommhold's corrections act as upper and lower bounds on the experimental spectrum. However, the deviations between profiles corrected using Frommhold's procedure and experimental data are slightly less compared to those obtained using Schofield's procedure.   (21)). The crosses show the values derived from the trajectory-based profiles using Equation (20). The shaded area indicates the bound states' contribution. Circles and triangles represent the values retrieved from the measurements by Dagg et al. (1986) and Birnbaum et al. (1993), respectively.

Observations
The Cassini/CIRS instrument consists of two interferometers, sharing a common telescope and scan mechanism. The far-infrared interferometer using Focal Plane 1 (hereafter referred to as FP1) covers the spectral range 10-650 cm −1 with an adjustable apodized spectral resolution from 15.5 to as high as 0.5 cm −1 . This focal plane consists of two thermopile detectors with a 3.9 mrad circular field of view with a half-peak diameter of 2.5 mrad (Flasar et al. 2004;Jennings et al. 2017). Bearing in mind the time (2005 January 14) and location (10.3°S , 167.7°E) of the Huygens landing site, we focused on spectra recorded between 2005 and 2007 in the 20°S-20°N latitude range. Seasonal variations of temperature and gas abundances are weak in the considered latitude range (Bézard et al. 2018), allowing us to use the constraints derived from the in situ measurements conducted by the Huygens probe. We opted to use the same two selections of 15.5 cm −1 resolution spectra as those made by Bézard & Vinatier (2020) to retrieve the H 2 mole fraction and ortho-to-para ratio. The two selections have different mean emission angles, "low" (13°) and "high" (59°), allowing us to more effectively disentangle the contributions from tropospheric CIA and stratospheric haze emission (Samuelson et al. 1997).

Radiative Transfer Model
Synthetic spectra were generated using the radiative transfer code described in Bézard & Vinatier (2020) and including CIA from N 2 -N 2 , CH 4 -N 2 , H 2 -N 2 and CH 4 -CH 4 , rotational lines of CO, HCN and CH 4 , rovibrational bands from C 4 H 2 and CH 3 C 2 H, and haze particles. The atmospheric model is the same as that of Bézard & Vinatier (2020) and makes use of in situ Huygens measurements of temperature (Fulchignoni et al. 2005), tropospheric methane abundance (Niemann et al. 2010), and main haze profile (Doose et al. 2016). Additional constraints on the aerosol opacity profile (including the nitrile ice cloud described by Anderson & Samuelson 2011) as well as C 4 H 2 and CH 3 C 2 H abundance profiles were obtained from CIRS limb spectra (at viewing altitudes of 105 and 145 km) recorded in the same latitude and time ranges as the two surface-intercepting spectral selections (Bézard & Vinatier 2020).
The N 2 -N 2 and CH 4 -N 2 CIA coefficients were updated with respect to Bézard & Vinatierʼs (2020) analysis. We opted for the use of the most recent N 2 -N 2 CIA coefficients calculated Figure 7. CH 4 -N 2 collision-induced absorption spectra at selected temperatures. Markers denote experimental data by Birnbaum et al. (1993) at 162, 195, and 297 K. The blue and red lines indicate the trajectory-based results obtained using Schofield's and Frommhold's quantum correction factors. Model profiles of Borysow & Tang (1993) are represented with gray lines. using a trajectory-based approach (Chistikov et al. 2019;Gordon et al. 2022) utilizing ab initio PES and IDS of Karman et al. (2015). In contrast, Bézard & Vinatier (2020) modeled the N 2 -N 2 absorption by slightly adjusting the spectral profiles of Borysow & Frommhold (1986b). They introduced a wavenumber-and temperature-dependent multiplicative factor such that adjusted CIA coefficients in the far wing match the quantum mechanical calculations of Karman et al. (2015).  Borysow & Frommhold (1986b), and heuristic model of Bézard & Vinatier (2020). It is seen that in the vicinity of the absorption peak, trajectory-based results of Chistikov et al. (2019) demonstrate significant improvement compared to previous theoretical efforts. In the far wing, neither of the calculations demonstrate perfect agreement with the experimental measurements of Sung et al. (2016) at temperatures 78.3, 89.3, 109.6, and 129.0 K. However, for 89.3 and 129.0 K, the deviations between the trajectory-based results and laboratory-measured spectra do not exceed 15% over the rototranslational band.
The analysis of spectral selections at low (13°) and high (59°) emission angles was performed using CH 4 -N 2 CIA coefficients obtained for both quantum correction factors. The simulated radiances are shown in Figure 9. Bearing in mind the comparison of the calculated CH 4 -N 2 CIA coefficients with the experimental measurements, it seems reasonable to consider a semiempirical model in the form of linearly weighted desymmetrized profiles: a n n n a n = +  T  c  T  c  T  T  , , , where coefficients c 1 and c 2 are wavenumber-and temperatureindependent. The best-fitting values of c 1 and c 2 , which minimize the residuals between synthetic spectra and both CIRS data, were found to be 0.35 and 0.65, respectively. The fitted weights are congruent with the laboratory measurements, in particular, of Birnbaum et al. (1993), gravitating to trajectory-based profiles desymmetrized using Frommhold's procedure. Figure 10 shows the synthetic radiances obtained using the semiempirical model (23).
The slight mismatch between synthetic radiances and CIRS data in the 100-200 cm −1 range most likely indicates that the temperature profile we used (Fulchignoni et al. 2005) is slightly too warm in the lower stratosphere and/or tropopause region. To investigate this possibility, we retrieved a temperature profile by inverting the CIRS radiances of both spectra in the range 55-545 cm −1 , using the Huygens probe profile (Fulchignoni et al. 2005) as the a priori and a correlation length of 0.75 scale height in the inversion algorithm. The soderived temperature profile is colder than our initial profile by about 1 K in the 50-150 mbar region. We consider that this difference is physically acceptable given that the Huygens measurements were made at a single location and time while the CIRS spectra correspond to an average over 3 yr and 40°in latitude. The synthetic spectra simulated with our semiempirical CH 4 -N 2 CIA model and retrieved temperature profile, shown in Figure 10, excellently reproduce the CIRS spectra over the spectral range 50-500 cm −1 .

Implications for Archean Earth
Understanding the environmental conditions in the terrestrial atmosphere during the Archean eon-from 4 to 2.5 billion years ago-may shed light on the biological and geological evolution of our planet and Earthlike exoplanets (Krissansen-Totton et al. 2018b;Catling & Zahnle 2020). In the Archean, the atmosphere was devoid of oxygen, and the bulk gases, N 2 and CO 2 , were supplemented with reducing gases such as CH 4 (Catling 2014). High levels of biogenic methane, as well as hydrocarbon haze layers similar to Titan's, have been suggested to warm the early Earth (Pavlov et al. 2001). Pieces of evidence, such as pronounced fractionation of xenon isotopes (Zahnle et al. 2019) and mass-independent fractionation of sulfur isotopes (Zahnle et al. 2006), point to high mixing ratios of methane, some reaching 0.5%. In turn, it follows that  Karman et al. (2015) and Borysow & Frommhold (1986b), respectively. The adjusted profile utilized in Bézard & Vinatier (2020) is indicated with the green line. Figure 9. Far-infrared low (13°) and high (59°) emission angle spectra. CIRS FP1 radiances are represented with black circles. The synthetic profiles using CH 4 -N 2 CIA of Borysow & Tang (1993) and trajectory-based results desymmetrized with Schofield's and Frommhold's procedures are denoted with yellow and blue and red lines, respectively. The high-emission spectra are shifted by 0.5 erg s −1 cm −2 sr −1 /cm −1 for clarity. The minor emission feature near 330 cm −1 is due to methyl-acetylene (CH 3 C 2 H), while the broad absorption centered at 355 cm −1 is due to the S 0 (0) H 2 line.
interacting pairs of CH 4 -N 2 and CH 4 -CO 2 may need to be considered as potential contributors to greenhouse warming relevant in the context of the faint young Sun problem (see, e.g., a recent review by Charnay et al. 2020). CIA is the major component of the opacity of gas giants and Titan in the dense regions of the atmosphere. However, its implications for early Earth warming are not well studied. While N 2 -N 2 and CO 2-CO 2 absorption is strong shortward of 300 cm −1 , it is too far from the peak of blackbody emission at terrestrial temperatures to have a significant greenhouse effect. The CH 4 -N 2 and CH 4 -CO 2 CIA bands are broader and may be important in the critical 450-700 cm −1 range near the maximum of outgoing radiation, where the contribution of other gases to opacity is small. Figure 11 shows the N 2 -N 2 , CH 4 -N 2 , CO 2 -CO 2 , and CH 4 -CO 2 collision-induced components assuming a model Archean and modern atmospheres. As illustrated, the opacities due to methanecontaining pairs were significantly larger in the Archean atmosphere than in the modern one. Despite that, it can be concluded that the CIA due to CH 4 -N 2 and CH 4 -CO 2 pairs is of little importance for infrared radiative transfer in the Archean atmosphere taking the constraint for CH 4 partial pressure inferred from the fractionation of xenon isotopes (Zahnle et al. 2019). Our results corroborate the findings of Pavlov et al. (2000), who concluded that CH 4 -N 2 CIA coefficients parameterized by Borysow & Tangʼs (1993) model are insignificant for greenhouse warming of early Earth. Nevertheless, should the constraints for methane mixing ratio be revised, the contribution of methanecontaining intermolecular pairs may need to be revisited.

Conclusions
This work presents an extensive trajectory-based examination of the CH 4 -N 2 rototranslational CIA band. Our firstprinciples analysis relies on the approach presented in Chistikov et al. (2021), allowing us to simulate the spectral contributions corresponding to both bound and unbound pair states. Significant progress has been made with respect to previous consideration of CH 4 -N 2 CIA by Borysow & Tang (1993) in virtue of the use of high-level ab initio anisotropic PES and IDS , obviating the need for short-range empirical dipole terms. We took advantage of the possibility to represent the lower-order spectral moments as phase-space and frequency-domain integrals to utilize them as the convergence control parameters of our Monte Carlo approach.
Classically derived trajectory-based profiles were subject to Schofield's (8) and Frommhold's (9) desymmetrization procedures. These profiles were then compared to available experimental data in two ways. First, we demonstrate an overall satisfactory agreement of profiles issued from both desymmetrization procedures with the spectra measured by Birnbaum et al. (1993) and Dagg et al. (1986). Second, the simulated profiles were utilized for modeling Titan's spectra in the 50-500 cm −1 range. A semiempirical model in the form of linearly weighted desymmetrized profiles was designed to reproduce Cassini CIRS spectra at low and high emission angles. The synthetic spectra simulated with semiempirical CH 4 -N 2 CIA spectra and a slightly modified temperature profile yielded an excellent agreement with CIRS data. For that reason, we believe that the semiempirical CH 4 -N 2 CIA spectra we obtained are the most accurate to date and may be beneficial for future Titan studies. The semiempirical CH 4 -N 2 CIA profiles from 70 to 400 K with a step of 10 K tabulated in the HITRAN CIA format (see Karman et al. 2019 and HITRAN CIA website 7 ) are provided in the .tar.gz supplementary package. These results are proposed to be included in the subsequent edition of the CIA section of the HITRAN database (Gordon et al. 2022).
Since existing laboratory measurements of CH 4 -N 2 CIA are limited to temperatures above 126 K (Dagg et al. 1986), we encourage laboratory studies down to the lowest available temperature to validate our semiempirical and simulated profiles. Naturally, since our calculations were carried out only up to 400 K and there are no experimental data above that temperature, higher temperature experiments would also be very welcome to account for possible scenarios in exoplanetary atmospheres.
We explored the potential implications of the simulated CH 4 -N 2 CIA profiles for greenhouse warming of early Earth's atmosphere. Based on the current constraints for Archean levels of methane mixing ratio, we concluded that CH 4 -N 2 CIA coefficients are insignificant for far-infrared radiative transfer in the Archean atmosphere, corroborating the findings of Pavlov et al. (2000).
The highly reduced atmospheres on rocky exoplanets are theoretically supported, assuming reduced iron can be retained in the mantle (Lichtenberg 2021). The atmospheres of many super-Earth exoplanets may resemble a warm Titan composition. For instance, Turbet et al. (2018) performed climate simulations to analyze the possibility of TRAPPIST-1 outer planets (Gillon et al. 2017) being warm Titans. We believe that our CH 4 -N 2 CIA profiles could play an important role in radiative transfer models for the climate simulations for warm N 2 /CH 4 exoplanetary atmospheres.
Trajectory-based calculations were carried out using HPC resources of Smithsonian Institution High Performance Cluster (SI/HPC). 8 A.A.F. and I.E.G.'s contributions to this work were supported by NASA (grant Nos. 80NSSC20K0962 and 80NSSC20K1059). The spatial configuration of the CH 4 -N 2 complex can be given in either body-fixed or space-fixed coordinate systems. To propagate the trajectory, one should be able to express body-fixed angles in terms of the space-fixed ones for the same configuration of the moieties since PES and IDS are expressed in terms of the body-fixed coordinates. The basic idea of the procedure is that given the initial reference orientation of the moieties, an arbitrary configuration can be reached by the sequence of rotations through either body-fixed or spaced-fixed angles expressed as the product of rotation matrices. We select two vectors embedded in the complex configuration, and by equating their coordinates in body-fixed and space-fixed frames, we can derive the values of the body-fixed coordinates. For the intermolecular vector R, we have Note that the Euler angles ( ) j J y , , defining the orientation of the body-fixed frame with respect to the space-fixed frame are the same as the ones that describe the orientation of the CH 4 molecule in the space-fixed frame. Defining the auxiliary vector  Figure 11. Far-infrared collision-induced absorption from selected pairs for the Archean atmosphere (left panel) with 1 atm of N 2 , 0.1 atm of CO 2 , 5 · 10 −3 atm of CH 4 , and surface temperature of 300 K, and for the modern Earth conditions (right panel). The N 2 partial pressure is based on the constraint derived from the analysis of fluid inclusions in Archean rocks (Avice et al. 2018). The partial pressure of CO 2 is selected as the largest value predicted by the carbonate-silicate climate model of Krissansen-Totton et al. (2018a). Methane partial pressure is based on the constraint inferred from the fractionation of xenon isotopes (Zahnle et al. 2019). A mean global temperature of 300 K is assumed based on studies of the oxygen isotope ratios in sedimentary evidence (Hren et al. 2009;Blake et al. 2010). CIA coefficients for CO 2 -CO 2 , N 2 -N 2 , and CH 4 -CO 2 are taken from Gruszka & Borysow (1997), Chistikov et al. (2019), and Wordsworth et al. (2017), respectively. The semiclassical model (Equation (23)) is utilized to calculate the CH 4 -N 2 CIA contribution. Note significantly different scales on the panels.
one can obtain the body-fixed angles Θ BF and Φ BF through the following relations: where atan2 is the two-argument arctangent, which takes into account the signs of both arguments to determine the quadrant of the result. Similarly, let us consider the vector l directed along the N 2 molecule. Its coordinates in the space-fixed and body-fixed frames of references are given by The change of variables induces a Jacobian factor that can be obtained as a determinant of the   matrix: det sin sin sin . B4 3 2 1 3 2 2 2 The distribution ρ 0 in variables x splits into a product of independent Gaussian distributions: .