# Trajectory-based simulation of far-infrared CIA profiles of CH${}_{4}-$N2
for modeling Titan’s atmosphere

Finenko A. Artem Center for Astrophysics |Harvard & Smithsonian, Atomic and
Molecular Physics Division, 60 Garden Street, Cambridge MA 02138, USA
Department of Chemistry, Lomonosov Moscow State University, GSP-1, 1-3
Leninskiye Gory, Moscow 119991, Russia  Bruno Bézard LESIA, Observatoire de
Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, place
Jules Janssen, 92195 Meudon, France  Iouli E. Gordon Center for Astrophysics
|Harvard & Smithsonian, Atomic and Molecular Physics Division, 60 Garden
Street, Cambridge MA 02138, USA  Daniil N. Chistikov  Department of Chemistry,
Lomonosov Moscow State University, GSP-1, 1-3 Leninskiye Gory, Moscow 119991,
Russia  Obukhov Institute of Atmospheric Physics, Russian Academy of Sciences,
3 Pyzhevsky per., Moscow 119017, Russia  Institute of Quantum Physics, Irkutsk
National Research Technical University, 83 Lermontov str., Irkutsk 664074,
Russia  Sergei E. Lokshtanov Department of Chemistry, Lomonosov Moscow State
University, GSP-1, 1-3 Leninskiye Gory, Moscow 119991, Russia  Sergey V.
Petrov  Department of Chemistry, Lomonosov Moscow State University, GSP-1, 1-3
Leninskiye Gory, Moscow 119991, Russia  Andrey A. Vigasin Obukhov Institute of
Atmospheric Physics, Russian Academy of Sciences, 3 Pyzhevsky per., Moscow
119017, Russia

(Accepted November 2, 2021)

###### Abstract

We report the results of the trajectory-based simulation of far-infrared
collision-induced absorption (CIA) due to CH${}_{4}-$N2 pairs at temperatures
between 70 and 400 K. Our analysis utilizes recently calculated high-level
potential energy (PES) and induced dipole surfaces (IDS) [Finenko, A. A.,
Chistikov, D. N., Kalugina, Y. N., Conway E. K., Gordon, I. E., Phys. Chem.
Chem. Phys., 2021, doi: 10.1039/d1cp02161c]. 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 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 semi-
empirical model for CH${}_{4}-$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 (CIRS) aboard the Cassini spacecraft at low and high emission
angles.

Titan, collision-induced absorption, infrared observations, radiative transfer

##  1 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}-$N2 collision-induced absorption 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 N2, CH4, and
H2 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}–$N2, CH${}_{4}-$N2, CH${}_{4}–$CH4,
N${}_{2}–$H2, H${}_{2}–$H2, CH${}_{4}–$H2. Due to the large abundance of N2
and relatively small rotational constant, $B\approx 2$ cm-1, the N${}_{2}-$N2
CIA dominates the tropospheric opacity for wavenumbers shortward of $\sim 150$
cm-1, whereas the CH${}_{4}-$N2 CIA is the prevailing source of opacity in the
150–450 cm-1 region (Samuelson et al., 1997). Several previous analyses
(Anderson & Samuelson, 2011; Tomasko et al., 2008; de Kok et al., 2010)
concluded that the opacity due to CH${}_{4}-$N2 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 Gas Chromatograph Mass
Spectrometer (GCMS) 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 (1993)’s 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 (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
the 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.

Figure 1: Normal optical thicknesses of Titan’s atmosphere. Collision-induced
intensities due to N${}_{2}-$N2 (Chistikov et al., 2019), CH${}_{4}-$N2
(current work), CH${}_{4}-$CH4 (Borysow & Frommhold, 1987), and N${}_{2}-$H2
(Borysow & Frommhold, 1986a) are shown separately. The contribution from
CH${}_{4}-$N2 CIA is obtained employing semi-classical model as described in
Section 3 (see Eq. (23)). The calculations were performed utilizing atmosphere
model of Bézard & Vinatier (2020).

The most accurate laboratory measurements of collision-induced absorption in
CH${}_{4}-$N2 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 (Borysow & Frommhold, 1986a, b, c,
1987; 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 the 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}-$N2, Borysow & Tang (1993) opted to perform fitting to measurements
by Birnbaum et al. (1993).

Recently, sophisticated first-principles approaches were developed and applied
to simulate far-infrared absorption in N${}_{2}-$N2 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}-$N2 (Chistikov et al., 2019;
Gordon et al., 2021) and CO${}_{2}-$Ar (Chistikov et al., 2021), we extended
it towards linear molecules interacting with spherical top molecules. Our
consideration relies on the ab initio potential energy and induced dipole
surfaces for CH${}_{4}-$N2 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}-$N2 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 semi-empirical model based on the ab initio
simulated CH${}_{4}-$N2 CIA profiles allowed us to reproduce Titan’s spectra
between 50 and 650 cm-1 recorded by the Composite Infrared Spectrometer (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.

##  2 Methodology of CIA spectrum calculation

###  2.1 Summary

In first-order perturbation theory, the binary absorption coefficient
$\alpha(\nu,T)$ in a mixture of gases can be represented as (Frommhold, 2006)

$\displaystyle\begin{aligned}
\alpha(\nu,T)=&\frac{\tau\left(\nu,T\right)}{\rho_{1}\rho_{2}}=\\\
&\frac{\left(2\pi\right)^{4}N_{L}^{2}}{3hc}\nu\left(1-\exp\left(-\beta
hc\nu\right)\right)VG(\nu,T),\end{aligned}$ (1)

where

$\displaystyle\tau(\nu,T)=L^{-1}\ln\left(I_{0}/I\right)$ (2)

is the absorption coefficient, $\rho_{1}$ and $\rho_{2}$ are the number
densities of gas constituents; $\beta=\left(k_{\text{B}}T\right)^{-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 $\nu$ is the wavenumber in
the reciprocal centimeters. The spectral function $G(\nu,T)$ is given by the
Fourier transform of the time correlation function (TCF) $C(t,T)$

$\displaystyle
G(\nu,T)=\frac{1}{2\pi}\int\limits_{-\infty}^{+\infty}C(t,T)e^{-2\pi ic\nu
t}\mathop{}\\!\mathrm{d}{t},$ (3)

which is defined as the statistically-averaged product of the time-shifted
dipoles relevant to the compound molecular pair

$\displaystyle
C(t,T)=\frac{V}{4\pi\varepsilon_{0}}\langle\boldsymbol{\mu}(0)\cdot\boldsymbol{\mu}(t)\rangle,$
(4)

where $\varepsilon_{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)

$\displaystyle G(-\nu,T)=e^{-\beta hc\nu}G(\nu,T)$ (5)

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

$\displaystyle G_{\text{cl}}(-\nu,T)=G_{\text{cl}}(\nu,T),$ (6)

in violation of Eq. (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 the so-called desymmetrization procedures (Frommhold, 2006; Borysow et
al., 1985). These are ad hoc prescriptions for constructing spectral functions
that obey 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; Berens et al.,
1981; Bader & Berne, 1994; Egelstaff, 1962). Let us formulate several known
prescriptions in terms of the quantum correction factor $\mathcal{F}(\nu,T)$
defined by

$\displaystyle G(\nu,T)=\mathcal{F}(\nu,T)G_{\text{cl}}(\nu,T).$ (7)

Schofield’s procedure (Schofield, 1960)

$\displaystyle\mathcal{F}_{\text{Sch}}(\nu,T)=e^{\beta hc\nu/2}$ (8)

was demonstrated to provide a satisfactory semi-classical 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

$\displaystyle\mathcal{F}_{\text{Fromm}}(\nu,T)=\frac{e^{\beta
hc\nu/2}}{G_{\text{cl}}(\nu,T)}\frac{C_{\text{cl}}(0)}{C_{\text{cl}}(\beta
hc\nu/\sqrt{2})}\int\limits_{-\infty}^{+\infty}C_{\text{cl}}\left(\sqrt{t^{2}+\left(\beta\hbar/2\right)^{2}},T\right)e^{-2\pi
ic\nu t}\mathop{}\\!\mathrm{d}{t}$ (9)

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 semi-classical binary absorption coefficient

$\displaystyle\alpha_{\text{SC}}(\nu,T)=\frac{\left(2\pi\right)^{4}N_{L}^{2}}{3hc}\nu\left(1-\exp\left(-\beta
hc\nu\right)\right)\mathcal{F}(\nu,T)VG_{\text{cl}}(\nu,T),$ (10)

where quantum correction factor $\mathcal{F}(\nu,T)$ is represented either by
Schofield’s (8) or Frommhold’s procedure (9).

###  2.2 Trajectory-based approach

Invoking Hamilton 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
$\mathbf{R}=\left(R,\Phi,\Theta\right)$ be the vector connecting the monomers’
centers of mass. The Euler angles
$\zeta_{1}=\left(\varphi,\vartheta,\psi\right)$ and
$\zeta_{2}=\left(\eta,\chi\right)$ describe an active rotation of CH4 and N2
moiety, respectively, from an initial position, in which a reference frame
fixed on 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}-$N2
system breaks up into

$\displaystyle H=K_{\text{tr}}+K_{\text{rot}_{1}}+K_{\text{rot}_{2}}+U,$ (11)

where $U$ is potential energy function and the kinetic energies of
translational and rotational motions are given by

$\displaystyle\begin{aligned}
K_{\text{tr}}&=\frac{p_{R}^{2}}{2\mu}+\frac{p_{\Theta}^{2}}{2\mu
R^{2}}+\frac{p_{\Phi}^{2}}{2\mu R^{2}\sin^{2}\Theta},\\\
K_{\text{rot}_{1}}&=\frac{\left(p_{\varphi}-p_{\psi}\cos\vartheta\right)^{2}}{2I_{1}\sin^{2}\vartheta}+\frac{p_{\psi}^{2}+p_{\vartheta}^{2}}{2I_{1}},\\\
K_{\text{rot}_{2}}&=\frac{p_{\chi}^{2}}{2I_{2}}+\frac{p_{\eta}^{2}}{2I_{2}\sin^{2}\chi}.\end{aligned}$
(12)

Here $\mu$, $I_{1}$, and $I_{2}$ denote the reduced mass of the molecular pair
and moments of inertia of CH4 and N2, respectively. The momentum conjugated to
the generalized variable $\xi$ is designated via $p_{\xi}$. Hereafter, the
generalized coordinates and conjugated momenta are assembled into vectors
$\mathbf{q}$ and $\mathbf{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):

$\displaystyle\begin{aligned} \dot{R}&=\frac{p_{R}}{\mu},\\\
\dot{p}_{R}&=\frac{p_{\Theta}^{2}}{\mu R^{3}}+\frac{p_{\Phi}^{2}}{\mu
R^{3}\sin^{2}\Theta}-\frac{\partial U}{\partial R},\\\
\dot{\Phi}&=\frac{p_{\Phi}}{\mu R^{2}\sin^{2}\Theta},\\\
\dot{p}_{\Phi}&=-\frac{\partial U}{\partial\Phi},\\\
\dot{\Theta}&=\frac{p_{\Theta}}{\mu R^{2}},\\\
\dot{p}_{\Theta}&=\frac{p_{\Phi}^{2}\cos\Theta}{\mu
R^{2}\sin^{3}\Theta}-\frac{\partial U}{\partial\Theta},\\\
\dot{\varphi}&=\frac{p_{\varphi}-p_{\psi}\cos\vartheta}{I_{1}\sin^{2}\vartheta},\\\
\dot{p}_{\varphi}&=-\frac{\partial U}{\partial\varphi},\\\
\dot{\vartheta}&=\frac{p_{\vartheta}}{I_{1}},\\\
\dot{p}_{\vartheta}&=\frac{\left(p_{\varphi}-p_{\psi}\cos\vartheta\right)\left(p_{\varphi}\cos\vartheta-
p_{\psi}\right)}{I_{1}\sin^{3}\vartheta}-\frac{\partial
U}{\partial\vartheta},\\\
\dot{\psi}&=\frac{p_{\psi}-p_{\varphi}\cos\vartheta}{I_{1}\sin^{2}\vartheta},\\\
\dot{p}_{\psi}&=-\frac{\partial U}{\partial\psi},\\\
\dot{\eta}&=\frac{p_{\eta}}{I_{2}\sin^{2}\chi},\\\
\dot{p}_{\eta}&=-\frac{\partial U}{\partial\eta},\\\
\dot{\chi}&=\frac{p_{\chi}}{I_{2}},\\\
\dot{p}_{\chi}&=\frac{p_{\eta}^{2}\cos\chi}{I_{2}\sin^{3}\chi}-\frac{\partial
U}{\partial\chi}.\end{aligned}$ (13)

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 CH4 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 CH4 monomer. The orientation of the vector $\mathbf{R}$ is
described by angles $\left(\Phi^{\text{BF}},\Theta^{\text{BF}}\right)$. Frame
2, parallel to Frame 1 fixed to CH4 molecule, has its origin in the center of
mass of N2. Orientation of N2 molecule with respect to Frame 2 is described by
angles $\left(\eta^{\text{BF}},\chi^{\text{BF}}\right)$. 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 $\Phi,\Theta,\varphi,\vartheta,\psi,\eta,\chi$. The numerical
scheme of calculating the desired derivatives is described in Appendix A.
Transformation of angles.

Figure 2: Body-fixed coordinate frame of the CH${}_{4}-$N2 complex.

From a formal point of view, one can opt for direct simulation of the TCF
using classical trajectories method for the purpose of then transforming it to
spectral function according to Eq. (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 $\nu$, 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 $G_{\text{cl}}^{\text{F+Q}}$ can
be written as

$\displaystyle G_{\text{cl}}^{\text{F+Q}}(\nu,T)\approx\frac{1}{2\pi
Q(T)}\int\limits_{-\infty}^{0}\frac{-p_{R}}{\mu}\mathop{}\\!\mathrm{d}{p_{R}}\idotsint\limits_{H>0}\exp\left(-\beta
H\right)\mathop{}\\!\mathrm{d}{\mathbf{q}^{\star}}\mathop{}\\!\mathrm{d}{\mathbf{p}^{\star}}\Bigg{|}\int\limits_{-\infty}^{+\infty}\boldsymbol{\mu}(t)e^{-2\pi
ic\nu t}\mathop{}\\!\mathrm{d}{t}\Bigg{|}^{2}.$ (14)

Here, $Q(T)$ represents the partition sum

$\displaystyle Q(T)=\int e^{-\beta
H}\mathop{}\\!\mathrm{d}{\mathbf{q}}\mathop{}\\!\mathrm{d}{\mathbf{p}}.$ (15)

An instantaneous value of the dipole moment vector of the collisional complex
is designated with $\boldsymbol{\mu}(t)$, and vectors $\mathbf{q}^{\star}$ and
$\mathbf{p}^{\star}$ are obtained upon exclusion of the coordinate $R$ and
conjugated momentum $p_{R}$ from the vectors $\mathbf{q}$ and $\mathbf{p}$,
respectively.

Integral representation (14) should be understood in the following way. First,
some fixed value of intermolecular distance $R_{\text{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 fly-by
trajectories for which the intermolecular distance $R$ always exceeds
$R_{\text{max}}$ that are naturally excluded from the expression (14). The
larger value of $R_{\text{max}}$ is chosen, the more negligible 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 quasi-periodic 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_{\text{cl}}^{\text{F+Q}}$. The sampling of quasi-periodic
trajectories through representation (14) can be achieved only when the value
of $R_{\text{max}}$ is somewhat close to the potential well. Thus it becomes
clear that no value of $R_{\text{max}}$ allows for taking into account
contributions from all significant trajectories. The spectral signatures from
the quasi-periodic trajectories are restricted to comparatively low
frequencies and hence are adequately simulated through the calculation of the
TCF. Selecting a large value of $R_{\text{max}}$ for the representation (14)
permits simulation of the spectral profile, excluding the quasi-periodic
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 $\exp\left(-\beta
H\right)$ over the domain of unbound states. The sampling algorithm is
outlined in Appendix B. Initial condition generation. Given a supply of $N$
random variables
$\left\\{\mathbf{q}_{k}^{\star},p_{R,k},\mathbf{p}_{k}^{\star}\right\\}$, the
Monte Carlo estimate is given by (for a more detailed discussion see Chistikov
et al. (2019))

$\displaystyle G_{\text{cl}}^{\text{F+Q}}(\nu,T)\approx\frac{3}{4\pi
R_{\text{max}}}\frac{1}{N}\sum_{k=1}^{N}\frac{-p_{R,k}}{\mu}\Bigg{|}\int\limits_{-\infty}^{+\infty}\boldsymbol{\mu}\left(t;R_{\text{max}},\mathbf{q}_{k}^{\star},p_{R,k},\mathbf{p}_{k}^{\star}\right)e^{-2\pi
ic\nu t}\mathop{}\\!\mathrm{d}{t}\Bigg{|}^{2}.$ (16)

Here, $\boldsymbol{\mu}\left(t;\cdot\right)$ is the value of dipole moment at
the time $t$ during classical evolution of 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. Initial condition
generation. Given $N$ phase space points
$\left\\{\mathbf{q}_{k},\mathbf{p}_{k}\right\\}$ supplied by the rejection
algorithm, the TCF in the domain $\Omega=\left\\{H<0\right\\}$ or
$\left\\{H>0\right\\}$ is obtained according to

$\displaystyle\begin{aligned}
C_{\Omega}(t;T)=\frac{1}{4\pi\varepsilon_{0}}\frac{V}{N}\frac{Q_{\Omega}\left(R_{\text{min}},R_{\text{max}};T\right)}{Q\left(T\right)}\sum_{i=1}^{N}\frac{1}{2}\Big{[}\boldsymbol{\mu}\left(0;\mathbf{q}_{i},\mathbf{p}_{i}\right)\boldsymbol{\mu}\left(t;\mathbf{q}_{i},\mathbf{p}_{i}\right)+\boldsymbol{\mu}\left(0;\mathbf{q}_{i},\mathbf{p}_{i}\right)\boldsymbol{\mu}\left(-t;\mathbf{q}_{i},\mathbf{p}_{i}\right)\Big{]},\end{aligned}$
(17)

where $Q_{\Omega}\left(R_{\text{min}},R_{\text{max}};T\right)$ is the partial
partition sum defined by

$\displaystyle
Q_{\Omega}\left(R_{1},R_{2};T\right)=\int\limits_{R_{1}}^{R_{2}}\mathop{}\\!\mathrm{d}{R}\int\limits_{-\infty}^{+\infty}\mathop{}\\!\mathrm{d}{p_{R}}\idotsint\limits_{\Omega}\exp\left(-\frac{H}{k_{\text{B}}T}\right)\mathop{}\\!\mathrm{d}{\mathbf{q}^{\star}}\mathop{}\\!\mathrm{d}{\mathbf{p}^{\star}}.$
(18)

Averaging of the products of dipole moments in Eq. (17), containing either
positive or negative time arguments, ensures time-reversibility of the TCF
estimate: $C_{\Omega}(t;T)=C_{\Omega}(-t;T)$. In the case of rigid linear
molecule-spherical top complex, the total partition sum can be expressed via
the molecular complex parameters (see Chistikov et al. (2018) for related
discussion)

$\displaystyle Q(T)=32\pi^{3}V\left(2\pi
k_{\text{B}}T\right)^{4}\mu^{3/2}I_{1}^{3/2}I_{2}.$ (19)

The adaptive Monte Carlo integration method VEGAS111A C++11 template library
for Monte Carlo integration: https://github.com/cschwan/hep-mc (Lepage, 1978)
is utilized to compute the partial partition sum
$Q_{\Omega}\left(R_{\text{min}},R_{\text{max}};T\right)$ prior to the
calculation of the TCF. We opted to use $R_{\text{min}}=4.0$ $a_{0}$ and
$R_{\text{max}}=35.0$ $a_{0}$ for the calculations described above.

The dipole variations along individual trajectories in Eq. (17) can be re-
weighted 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 Monte Carlo staged multitemperature 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_{\text{max}}=1.58$ ns is reached. 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 Eq.
(3). Direct application of standard numerical procedures implementing 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 the variable-width
window to remove the residual Monte Carlo noise. Subsequently, absorption
coefficients are produced according to Eq. (10) with properly chosen quantum
correction factors.

###  2.3 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 $\nu=200$ cm-1. The stitching is
performed independently for the spectral profiles obtained using Schofield’s
and Frommhold’s quantum correction factors. 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 dipole-multipole 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 N2 molecule has an almost
identical intensity for Schofield’s and Frommhold’s quantum corrections. It
should be noted that the true 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.

Figure 3: CH${}_{4}-$N2 CIA spectral profiles corresponding to unbound states
obtained through TCFs processing (in color) and through formalism (14) (in
black) at 70–90 K. Both sets of profiles utilize Schofield’s quantum
correction factor. The TCF-based profiles longward of 200 cm-1 are shown with
dashed lines. Resultant profiles are obtained upon the stitching of two sets
at 200 cm-1.

Figure 4: CH${}_{4}-$N2 CIA spectral profiles corresponding to both bound and
unbound states (in color) from 70 to 400 K with a step of 10 K. The black
lines represent profiles corresponding to unbound states at 70–100 K. The left
panel shows profiles subjected to Schofield’s quantum correction factor,
whereas Frommhold’s correction is utilized for the profiles on the right
panel.

###  2.4 Spectral moments as convergence control parameters

The convergence of 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

$\displaystyle
M_{n}=2\int\limits_{0}^{\infty}\nu^{n-1}\left[1-\exp\left(-\frac{hc\nu}{k_{\text{B}}T}\right)\right]^{-1}\alpha\left(\nu,T\right)\mathop{}\\!\mathrm{d}{\nu}.$
(20)

On the other hand, the lower-order spectral moments can be obtained through
phase-space averaging according to (see, e.g., Frommhold (2006))

$\displaystyle\begin{aligned}
M_{0}&=\frac{\left(2\pi\right)^{4}N_{L}^{2}}{3h}\frac{V}{4\pi\varepsilon_{0}}\frac{1}{2\pi
c}\langle\boldsymbol{\mu}(0)^{2}\rangle,\\\
M_{2}&=\frac{\left(2\pi\right)^{4}N_{L}^{2}}{3h}\frac{V}{4\pi\varepsilon_{0}}\frac{1}{\left(2\pi
c\right)^{3}}\langle\dot{\boldsymbol{\mu}}(0)^{2}\rangle.\end{aligned}$ (21)

When comparing spectral moments obtained through phase-space integration with
the values derived from experimental spectra, one should take into account
that the corresponding spectra obey different detailed balance conditions
shown in Eqs. (5) and (6) (the issue is expounded upon in Finenko et al.
(2021) and Chistikov 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, Eq. (6). Consequently, prior to computing the
moments from the laboratory-measured spectra using Eq. (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

$\displaystyle\mathcal{F}_{\text{inv-Sch}}\left(\nu,T\right)=e^{-\beta
hc\nu/2}.$ (22)

In Fig. 5, classical spectral moments obtained through phase space integration
(21) are compared to the values derived from the trajectory-based spectral
profiles via Eq. (20) as well as from the laboratory-measured spectra. Perfect
agreement between the spectral moments obtained via frequency-domain and
phase-space integration validates calculated spectral profiles.

Figure 5: Calculated and experimental temperature variations of zeroth (left
panel) and second (right panel) spectral moments of the CH${}_{4}-$N2
rototranslational band. The solid line refers to the values obtained through
phase-space integration (21). The crosses show the values derived from the
trajectory-based profiles using Eq. (20). The shaded area indicates the true
bound states’ contribution. Circles and triangles represent the values
retrieved from the measurements by Dagg et al. (1986) and Birnbaum et al.
(1993), respectively.

Figures 6 and 7 compare the calculated CIA CH${}_{4}-$N2 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.

Figure 6: CH${}_{4}-$N2 collision-induced absorption spectra at selected
temperatures. Dashed yellow lines denote experimental measurements by Dagg et
al. (1986) at 126, 149, 179, and 212 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.

Figure 7: CH${}_{4}-$N2 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.

##  3 Modeling CIRS/Cassini measurements

###  3.1 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 (14 January 2005) 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 H2 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).

###  3.2 Radiative transfer model

Synthetic spectra were generated using the radiative transfer code described
in Bézard & Vinatier (2020) and including CIA from N${}_{2}-$N2,
CH${}_{4}-$N2, H${}_{2}-$N2 and CH${}_{4}-$CH4, rotational lines of CO, HCN
and CH4, rovibrational bands from C4H2 and CH3C2H, 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 C4H2 and CH3C2H 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}-$N2 and CH${}_{4}-$N2 CIA coefficients were updated with respect
to Bézard & Vinatier (2020)’s analysis. We opted for the use of the most
recent N${}_{2}-$N2 CIA coefficients calculated using a trajectory-based
approach (Chistikov et al., 2019; Gordon et al., 2021) utilizing ab initio PES
and IDS of Karman et al. (2015). In contrast, Bézard & Vinatier (2020)
modelled the N${}_{2}-$N2 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). Figure 8 compares the experimental measurements by Sung
et al. (2016) at 89.3 K against the theoretical results by Chistikov et al.
(2019), 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.

Figure 8: N${}_{2}-$N2 collision-induced absorption spectra at 89.3 K. The red
curve indicates the results of the trajectory-based calculation (Chistikov et
al. (2019)) utilizing Frommhold’s desymmetrization procedure. The experimental
measurements of Sung et al. (2016) are denoted with the blue line. The yellow
diamonds and gray curve denote the results of quantum mechanical calculations
of Karman et al. (2015) and Borysow & Frommhold (1986b), respectively. The
adjusted profile utilized in Bézard & Vinatier (2020) is indicated with the
green line.

The analysis of spectral selections at low (13∘) and high (59∘) emission
angles was performed using CH${}_{4}-$N2 CIA coefficients obtained for both
quantum correction factors. The simulated radiances are shown in Fig. 9.
Bearing in mind the comparison of the calculated CH${}_{4}-$N2 CIA
coefficients with the experimental measurements, it seems reasonable to
consider a semi-empirical model in the form of linearly-weighted desymmetrized
profiles

$\displaystyle\alpha_{\text{SE}}(\nu,T)=\Big{[}c_{1}\mathcal{F}_{\text{Sch}}\left(\nu,T\right)+c_{2}\mathcal{F}_{\text{Fromm}}\left(\nu,T\right)\Big{]}\alpha_{\text{cl}}\left(\nu,T\right),$
(23)

where coefficients $c_{1},c_{2}$ are wavenumber- and temperature-independent.
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 semi-empirical 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 so-derived 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 years and 40∘ in latitude. The synthetic spectra simulated with
our semi-empirical CH${}_{4}-$N2 CIA model and retrieved temperature profile,
shown in Figure 10, excellently reproduce the CIRS spectra over the spectral
range 50–500 cm-1.

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}-$N2 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 (CH3C2H), while the broad absorption
centered at 355 cm-1 is due to the S${}_{0}(0)$ H2 line. Figure 10: Far-
infrared low (13∘) and high (59∘) emission angle spectra. CIRS FP1 radiances
are represented with black circles. The synthetic profiles simulated using
semi-empirical CH${}_{4}-$N2 CIA model (23) are denoted with orange lines. The
radiances simulated using semi-empirical CIA CH${}_{4}-$N2, as well as
modified temperature profile, are represented with magenta lines. The high-
emission spectra are shifted by 0.5 erg s-1 cm-2 sr-1/cm-1 for clarity.

##  4 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 Earth-like
exoplanets (Catling & Zahnle, 2020; Krissansen-Totton et al., 2018b). In the
Archean, the atmosphere was devoid of oxygen, and the bulk gases, N2 and CO2,
were supplemented with reducing gases such as CH4 (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 interacting pairs of CH${}_{4}-$N2 and CH${}_{4}-$CO2 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)). Collision-induced absorption 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}-$N2 and CO${}_{2}-$CO2 absorption is strong shortward of 300 cm-1, it
is too far from the peak of black-body emission at terrestrial temperatures to
have a significant greenhouse effect. The CH${}_{4}-$N2 and CH${}_{4}-$CO2 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}-$N2, CH${}_{4}-$N2,
CO${}_{2}-$CO2, and CH${}_{4}-$CO2 collision-induced components assuming a
model Archean and modern atmospheres. As illustrated, the opacities due to
methane-containing 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}-$N2 and CH${}_{4}-$CO2 pairs is of little importance for infrared
radiative transfer in the Archean atmosphere taking the constraint for CH4
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}-$N2 CIA coefficients parameterized by Borysow & Tang
(1993)’s model are insignificant for greenhouse warming of early Earth.
Nevertheless, should the constraints for methane mixing ratio be revised, the
contribution of methane-containing intermolecular pairs may need to be
revisited.

Figure 11: Far-infrared collision-induced absorption from selected pairs for
the Archean atmosphere (left panel) with 1 atm of N2, 0.1 atm of CO2, $5\cdot
10^{-3}$ atm of CH4, and surface temperature of 300 K, and for the modern
Earth conditions (right panel). The N2 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 CO2 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 (Blake et al., 2010; Hren et al., 2009). CIA
coefficients for CO${}_{2}-$CO2, N${}_{2}-$N2, and CH${}_{4}-$CO2 are taken
from Gruszka & Borysow (1997), Chistikov et al. (2019), and Wordsworth et al.
(2017), respectively. The semi-classical model (Eq. (23)) is utilized to
calculate the CH${}_{4}-$N2 CIA contribution. Note significantly different
scales on the panels.

##  5 Conclusions

This work presents an extensive trajectory-based examination of the
CH${}_{4}-$N2 rototranslational collision-induced absorption band. Our first-
principles 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}-$N2 CIA by Borysow & Tang
(1993) in virtue of the use of high-level ab initio anisotropic PES and IDS
(Finenko et al., 2021), 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 semi-empirical 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 semi-empirical CH${}_{4}-$N2 CIA spectra and a slightly modified
temperature profile yielded an excellent agreement with CIRS data. For that
reason, we believe that the semi-empirical CH${}_{4}-$N2 CIA spectra we
obtained are the most accurate to date and may be beneficial for future Titan
studies. The semi-empirical CH${}_{4}-$N2 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 website222https://hitran.org/cia/) are provided in the
Supplementary material. These results are proposed to be included in the
subsequent edition of the CIA section of the HITRAN database (Gordon et al.,
2021).

Since existing laboratory measurements of CH${}_{4}-$N2 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 semi-empirical 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}-$N2 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}-$N2 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}-$N2 CIA profiles
could play an important role in radiative transfer models for the climate
simulations for warm N2/CH4 exoplanetary atmospheres.

Trajectory-based calculations were carried out using HPC resources of
Smithsonian Institution High Performance Cluster
(SI/HPC)333https://doi.org/10.25572/SIHPC. This work is supported by the NASA
grants 80NSSC20K0962 and 80NSSC20K1059.

## A. Transformation of angles

The components of any vector $\mathbf{v}$ in the space-fixed system
$\mathbf{v}_{\textup{SF}}$ are linearly related to the components of the same
vector in the body-fixed reference frame $\mathbf{v}_{\textup{BF}}$ through
orthogonal transformation matrix $\mathbb{S}$ (Goldstein et al., 2001)

$\displaystyle\mathbb{S}\mathbf{v}_{\textup{SF}}=\mathbf{v}_{\textup{BF}}.$
(A.1)

The elements of the matrix $\mathbb{S}$ can be expressed in terms of three
Euler angles $\left(\varphi,\vartheta,\psi\right)$. In the present paper, we
adhere to the definition of Euler angles given by Goldstein et al. (2001)
leading to the following form of the transformation matrix $\mathbb{S}$

$\displaystyle\mathbb{S}(\varphi,\vartheta,\psi)=\mathbb{S}_{\textup{Z}}\left(\psi\right)\mathbb{S}_{\textup{X}}\left(\vartheta\right)\mathbb{S}_{\textup{Z}}\left(\varphi\right),$
(A.2)

where a matrix $\mathbb{S}_{\xi}(\zeta)$ implements a counterclockwise
rotation on angle $\zeta$ around $\xi$-axis. Rotation matrices about the $X$-,
$Y$\- and $Z$-axes are defined as

$\displaystyle\mathbb{S}_{\textup{X}}(\zeta)=\begin{bmatrix}1&0&0\\\
0&\cos\zeta&\sin\zeta\\\
0&-\sin\zeta&\cos\zeta\end{bmatrix},\quad\mathbb{S}_{\textup{Y}}(\zeta)=\begin{bmatrix}\cos\zeta&0&-\sin\zeta\\\
0&1&0\\\
\sin\zeta&0&\cos\zeta\end{bmatrix},\quad\mathbb{S}_{\textup{Z}}(\zeta)=\begin{bmatrix}\cos\zeta&\sin\zeta&0\\\
-\sin\zeta&\cos\zeta&0\\\ 0&0&1\end{bmatrix}.$ (A.3)

The spatial configuration of the CH${}_{4}-$N2 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
$\mathbf{R}$ we have

$\displaystyle\mathbf{R}_{\textup{SF}}=R\begin{bmatrix}\cos\Phi\sin\Theta\\\
\sin\Phi\sin\Theta\\\ \cos\Theta\end{bmatrix},$ (A.4)

and

$\displaystyle\mathbf{R}_{\textup{BF}}=R\begin{bmatrix}\cos\Phi^{\text{BF}}\sin\Theta^{\text{BF}}\\\
\sin\Phi^{\text{BF}}\sin\Theta^{\text{BF}}\\\
\cos\Theta^{\text{BF}}\end{bmatrix}.$ (A.5)

Bearing in mind Eq. (A.1) and rewriting the spherical coordinates in terms of
the rotation matrices, we obtain

$\displaystyle\mathbb{S}_{\textup{Z}}^{+}(\varphi)\mathbb{S}_{\textup{X}}^{+}(\vartheta)\mathbb{S}_{\textup{Z}}^{+}(\psi)\mathbb{S}_{\textup{Z}}^{+}(\Phi^{\text{BF}})\mathbb{S}_{\textup{Y}}^{+}(\Theta^{\text{BF}})\begin{bmatrix}0\\\
0\\\
1\end{bmatrix}=\mathbb{S}_{\textup{Z}}^{+}(\Phi)\mathbb{S}_{\textup{Y}}^{+}(\Theta)\begin{bmatrix}0\\\
0\\\ 1\end{bmatrix}.$ (A.6)

Note that the Euler angles $\left(\varphi,\vartheta,\psi\right)$ 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 CH4 molecule in the
space-fixed frame. Defining the auxiliary vector

$\displaystyle\mathbf{e}_{1}=\mathbb{S}_{\textup{Z}}(\psi)\mathbb{S}_{\textup{X}}(\vartheta)\mathbb{S}_{\textup{Z}}(\varphi)\mathbb{S}_{\textup{Z}}^{+}(\Phi),\mathbb{S}_{\textup{Y}}^{+}(\Theta)\begin{bmatrix}0\\\
0\\\ 1\end{bmatrix},$ (A.7)

one can obtain the body-fixed angles $\Theta^{\text{BF}}$ and
$\Phi^{\text{BF}}$ through the following relations

$\displaystyle\Theta^{\text{BF}}=\arccos\left(\mathbf{e}_{1}\right)_{Z},\quad\Phi^{\text{BF}}=\operatorname{atan2}\left(\frac{\left(\mathbf{e}_{1}\right)_{Y}}{\sin\Theta^{\text{BF}}},\frac{\left(\mathbf{e}_{1}\right)_{X}}{\sin\Theta^{\text{BF}}}\right),$
(A.8)

where $\operatorname{atan2}$ is the 2-argument arctangent, which takes into
account the signs of both arguments to determine the quadrant of the result.

Similarly, let us consider the vector $\mathbf{l}$ directed along the N2
molecule. Its coordinates in the space-fixed and body-fixed frames of
references are given by

$\displaystyle\mathbf{l}_{\text{SF}}=l\begin{bmatrix}\cos\eta\sin\chi\\\
\sin\eta\sin\chi\\\ \cos\chi\end{bmatrix}$ (A.9)

and

$\displaystyle\mathbf{l}_{\text{BF}}=l\begin{bmatrix}\cos\eta^{\text{BF}}\sin\chi^{\text{BF}}\\\
\sin\eta^{\text{BF}}\sin\chi^{\text{BF}}\\\
\cos\chi^{\text{BF}}\end{bmatrix},$ (A.10)

respectively, where $l$ is the length of the N2 molecule. Keeping in mind
(A.1), we have

$\displaystyle\mathbb{S}_{\textup{Z}}^{+}(\varphi)\mathbb{S}_{\textup{X}}^{+}(\vartheta)\mathbb{S}_{\textup{Z}}^{+}(\psi)\mathbb{S}_{\textup{Z}}^{+}(\eta^{\text{BF}})\mathbb{S}_{\textup{Y}}^{+}(\chi^{\text{BF}})\begin{bmatrix}0\\\
0\\\
1\end{bmatrix}=\mathbb{S}_{\textup{Z}}^{+}(\eta)\mathbb{S}_{\textup{Y}}^{+}(\chi)\begin{bmatrix}0\\\
0\\\ 1\end{bmatrix}.$ (A.11)

The body-fixed angles $\chi^{\text{BF}}$ and $\eta^{\text{BF}}$ angles can be
obtained according to

$\displaystyle\chi^{\text{BF}}=\arccos\left(\mathbf{e}_{2}\right)_{Z},\quad\eta^{\text{BF}}=\operatorname{atan2}\left(\frac{\left(\mathbf{e}_{2}\right)_{Y}}{\sin\chi^{\text{BF}}},\frac{\left(\mathbf{e}_{2}\right)_{X}}{\sin\chi^{\text{BF}}}\right),$
(A.12)

where the auxiliary vector $\mathbf{e}_{2}$ is given by

$\displaystyle\mathbf{e}_{2}=\mathbb{S}_{\textup{Z}}(\psi)\mathbb{S}_{\textup{X}}(\vartheta)\mathbb{S}_{\textup{Z}}(\varphi)\mathbb{S}_{\textup{Z}}^{+}(\eta)\mathbb{S}_{\textup{Y}}^{+}(\chi)\begin{bmatrix}0\\\
0\\\ 1\end{bmatrix}.$ (A.13)

To compute the derivatives of the PES with respect to the space-fixed angles
entering dynamical equations (13), we opt to employ the chain rule

$\displaystyle\frac{\partial
U}{\partial\xi}=\sum_{\zeta=\Phi^{\text{BF}},\Theta^{\text{BF}},\eta^{\text{BF}},\chi^{\text{BF}}}\frac{\partial
U}{\partial\zeta}\frac{\partial\zeta}{\partial\xi},\quad\xi=\Phi,\Theta,\varphi,\vartheta,\psi,\eta,\chi.$
(A.14)

The required derivatives of the space-fixed angles with respect to the body-
fixed can be readily obtained utilizing the same approach described above by
addressing matrix calculus. Let us introduce additional auxiliary vectors

$\displaystyle\mathbf{e}_{1}^{\zeta}=\frac{\partial}{\partial\zeta}\Bigg{[}\mathbb{S}_{\textup{Z}}(\psi)\mathbb{S}_{\textup{X}}(\vartheta)\mathbb{S}_{\textup{Z}}(\varphi)\mathbb{S}_{\textup{Z}}^{+}(\Phi)\mathbb{S}_{\textup{Y}}^{+}(\Theta)\Bigg{]}\begin{bmatrix}0\\\
0\\\ 1\end{bmatrix},$ (A.15)
$\displaystyle\mathbf{e}_{2}^{\zeta}=\frac{\partial}{\partial\zeta}\Bigg{[}\mathbb{S}_{\textup{Z}}(\psi)\mathbb{S}_{\textup{X}}(\vartheta)\mathbb{S}_{\textup{Z}}(\varphi)\mathbb{S}_{\textup{Z}}^{+}(\eta)\mathbb{S}_{\textup{Y}}^{+}(\chi)\Bigg{]}\begin{bmatrix}0\\\
0\\\ 1\end{bmatrix},$ (A.16)

where $\zeta=\\{\Phi,\Theta,\varphi,\vartheta,\psi,\eta,\chi\\}$ and
$\partial/\partial\zeta$ means a standard matrix derivative (Lax, 2007). Then
the corresponding derivatives can be computed via

$\displaystyle\frac{\partial\Theta^{\text{BF}}}{\partial\zeta}=-\frac{1}{\sin\Theta^{\text{BF}}}\left(\mathbf{e}_{1}^{d}\right)_{Z},\quad\frac{\partial\Phi^{\text{BF}}}{\partial\zeta}=\frac{1}{\sin\Theta^{\text{BF}}}\Bigg{[}\left(\mathbf{e}_{1}^{d}\right)_{Y}\cos\Phi^{\text{BF}}-\left(\mathbf{e}_{1}^{d}\right)_{X}\sin\Phi^{\text{BF}}\Bigg{]},$
(A.17)
$\displaystyle\frac{\partial\chi^{\text{BF}}}{\partial\zeta}=-\frac{1}{\sin\chi^{\text{BF}}}\left(\mathbf{e}_{2}^{d}\right)_{Z},\quad\frac{\partial\eta^{\text{BF}}}{\partial\zeta}=\frac{1}{\sin\chi^{\text{BF}}}\Bigg{[}\left(\mathbf{e}_{2}^{d}\right)_{Y}\cos\eta^{\text{BF}}-\left(\mathbf{e}_{2}^{d}\right)_{X}\sin\eta^{\text{BF}}\Bigg{]}.$
(A.18)

## B. Initial condition generation

In what follows, we outline the sampling algorithm from the Boltzmann
distribution derived under the assumption of zero intermolecular potential

$\displaystyle\rho_{0}\left(\mathbf{q},\mathbf{p}\right)\sim\exp\left(-\beta\left(K_{\text{tr}}+K_{\text{rot}_{1}}+K_{\text{rot}_{2}}\right)\right).$
(B.1)

Let us recall that the kinetic energy function is a quadratic form and can
therefore be diagonalized (see, e. g., Chistikov et al. (2018))

$\displaystyle
K_{\text{tr}}+K_{\text{rot}_{1}}+K_{\text{rot}_{2}}=\frac{1}{2}\mathbf{x}^{\top}\mathbf{x}=\frac{1}{2}\sum_{k=1}^{8}x_{k}^{2},$
(B.2)

where the relationship between generalized momenta $\mathbf{p}$ and variables
$\mathbf{x}$ is given by

$\displaystyle\mathbf{p}=\mathbb{T}^{\top}\mathbf{x}=\begin{bmatrix}\sqrt{\mu}&0&0&0&0&0&0&0\\\
0&\sqrt{\mu}R\sin\Theta&0&0&0&0&0&0\\\ 0&0&\sqrt{\mu}R&0&0&0&0&0\\\
0&0&0&\sqrt{I_{1}}\sin\vartheta&0&\sqrt{I_{1}}\cos\vartheta&0&0\\\
0&0&0&0&\sqrt{I_{1}}&0&0&0\\\ 0&0&0&0&0&\sqrt{I_{1}}&0&0\\\
0&0&0&0&0&0&\sqrt{I_{2}}\sin\chi&0\\\
0&0&0&0&0&0&0&\sqrt{I_{2}}\end{bmatrix}\mathbf{x}.$ (B.3)

The change of variables induces a Jacobian factor that can be obtained as a
determinant of the $\mathbb{T}^{\top}$ matrix

$\displaystyle\Bigg{|}\frac{\partial\left[\mathbf{q},\mathbf{x}\right]}{\partial\left[\mathbf{q},\mathbf{p}\right]}\Bigg{|}=\text{det}\left(\mathbb{T}^{\top}\right)=\mu^{3/2}I_{1}^{3/2}I_{2}R^{2}\sin\Theta\sin\vartheta\sin\chi.$
(B.4)

The distribution $\rho_{0}$ in variables $\mathbf{x}$ splits into a product of
independent Gaussian distributions

$\displaystyle\rho_{0}\left(\mathbf{q},\mathbf{x}\right)\sim
R^{2}\sin\Theta\sin\vartheta\sin\chi\prod_{k=1}^{8}\exp\left(-\frac{1}{2}\beta
x_{k}^{2}\right).$ (B.5)

It follows that the vector components $x_{k}$ are distributed normally
$x_{k}\sim\mathcal{N}(0,1)$, cosines of each of the angles
$\left\\{\Theta,\vartheta,\chi\right\\}$ are distributed uniformly in the
range of $\left[-1,1\right]$, and the angles
$\left\\{\Phi,\varphi,\psi,\eta\right\\}$ are distributed uniformly in the
range of $\left[0,2\pi\right]$. The samples of variable $R$ can be generated
through selecting $R^{3}$ uniformly on
$\left[R^{3}_{\text{min}},R^{3}_{\text{max}}\right]$.

The sampling from the full Boltzmann distribution

$\displaystyle\rho\left(\mathbf{q},\mathbf{p}\right)=\exp\left(-\frac{H}{k_{\text{B}}T}\right)$
(B.6)

is performed by applying the rejection criterion to the samples generated from
the $\rho_{0}$-density (see details in Section IV.B of Chistikov et al.
(2021)).

## References

  * Anderson & Samuelson (2011) Anderson, C. M., & Samuelson, R. E. 2011, Icarus, 212, 762, doi: 10.1016/j.icarus.2011.01.024
  * Avice et al. (2018) Avice, G., Marty, B., Burgess, R., et al. 2018, Geochimica et Cosmochimica Acta, 232, 82, doi: 10.1016/j.gca.2018.04.018
  * Bader & Berne (1994) Bader, J. S., & Berne, B. 1994, The Journal of chemical physics, 100, 8359
  * Berens et al. (1981) Berens, P. H., White, S. R., & Wilson, K. R. 1981, The Journal of Chemical Physics, 75, 515
  * Bézard & Vinatier (2020) Bézard, B., & Vinatier, S. 2020, Icarus, 344, 113261, doi: 10.1016/j.icarus.2019.03.038
  * Bézard et al. (2018) Bézard, B., Vinatier, S., & Achterberg, R. K. 2018, Icarus, 302, 437, doi: 10.1016/j.icarus.2017.11.034
  * Birnbaum et al. (1993) Birnbaum, G., Borysow, A., & Buechele, A. 1993, Journal of Chemical Physics, 99, 3234, doi: 10.1063/1.465132
  * Blake et al. (2010) Blake, R. E., Chang, S. J., & Lepland, A. 2010, Nature, 464, 1029, doi: 10.1038/nature08952
  * Borysow & Frommhold (1986a) Borysow, A., & Frommhold, L. 1986a, Astrophysical Journal, 303, 495, doi: 10.1086/164096
  * Borysow & Frommhold (1986b) —. 1986b, Astrophysical Journal, 311, 1043, doi: 10.1086/164841
  * Borysow & Frommhold (1986c) —. 1986c, Astrophysical Journal, 304, 849, doi: 10.1086/164221
  * Borysow & Frommhold (1987) —. 1987, Astrophysical Journal, 318, 940, doi: 10.1086/165426
  * Borysow & Moraldi (1994) Borysow, A., & Moraldi, M. 1994, Molecular Physics, 82, 1277, doi: 10.1080/00268979400100904
  * Borysow & Tang (1993) Borysow, A., & Tang, C. 1993, Icarus, 105, 175, doi: 10.1006/icar.1993.1117
  * Borysow et al. (1985) Borysow, J., Moraldi, M., & Frommhold, L. 1985, Molecular Physics, 56, 913, doi: 10.1080/00268978500102801
  * Bussery-Honvault & Hartmann (2014) Bussery-Honvault, B., & Hartmann, J.-M. 2014, Journal of Chemical Physics, 140, 054309, doi: 10.1063/1.4863636
  * Catling (2014) Catling, D. C. 2014, in Treatise on Geochemistry, ed. H. D. Holland & K. K. Turekian, Vol. 8 (Elsevier), 191–233
  * Catling & Zahnle (2020) Catling, D. C., & Zahnle, K. J. 2020, Science Advances, 6, eaax1420, doi: 10.1126/sciadv.aax1420
  * Charnay et al. (2020) Charnay, B., Wolf, E. T., Marty, B., & Forget, F. 2020, Space Science Reviews, 216, 90, doi: 10.1007/s11214-020-00711-9
  * Chistikov et al. (2021) Chistikov, D. N., Finenko, A. A., Kalugina, Y. N., et al. 2021, Journal of Chemical Physics, doi: 10.1063/5.0060779
  * Chistikov et al. (2018) Chistikov, D. N., Finenko, A. A., Lokshtanov, S. E., Petrov, S. V., & Vigasin, A. A. 2018, Journal of Chemical Physics, 149, 194304, doi: 10.1063/1.5054125
  * Chistikov et al. (2019) —. 2019, Journal of Chemical Physics, 151, 194106, doi: 10.1063/1.5125756
  * Courtin et al. (1995) Courtin, R., Gautier, D., & McKay, C. P. 1995, Icarus, 114, 144, doi: 10.1006/icar.1995.1050
  * Dagg et al. (1986) Dagg, I. R., Anderson, A., Yan, S., et al. 1986, Canadian Journal of Physics, 64, 1467, doi: 10.1139/p86-260
  * de Kok et al. (2010) de Kok, R., Irwin, P. G. J., & Teanby, N. A. 2010, Icarus, 209, 854, doi: 10.1016/j.icarus.2010.06.035
  * Doose et al. (2016) Doose, L. R., Karkoschka, E., Tomasko, M. G., & Anderson, C. M. 2016, Icarus, 270, 355, doi: 10.1016/j.icarus.2015.09.039
  * Egelstaff (1962) Egelstaff, P. A. 1962, Advances in Physics, 11, 203
  * Finenko et al. (2021) Finenko, A. A., Chistikov, D. N., Kalugina, Y. N., Conway, E. K., & Gordon, I. E. 2021, Physical Chemistry Chemical Physics, doi: 10.1039/d1cp02161c
  * Flasar et al. (2004) Flasar, F. M., Kunde, V. G., Abbas, M. M., et al. 2004, Space Science Reviews, 115, 169, doi: 10.1007/s11214-004-1454-9
  * Frommhold (2006) Frommhold, L. 2006, Collision Induced Absorption in Gases (Cambridge University Press)
  * Fulchignoni et al. (2005) Fulchignoni, M., Ferri, F., Angrilli, F., et al. 2005, Nature, 438, 785, doi: 10.1038/nature04314
  * Gillon et al. (2017) Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542, 456, doi: 10.1038/nature21360
  * Goldstein et al. (2001) Goldstein, H., Poole, C. P., & Safko, J. L. 2001, Classical Mechanics (Addison Wesley Publishing Company, San Francisco)
  * Gordon et al. (2021) Gordon, I. E., Rothman, L. S., Hargreaves, R. J., et al. 2021, Journal of Quantitative Spectroscopy and Radiative Transfer, 107949, doi: 10.1016/j.jqsrt.2021.107949
  * Gruszka & Borysow (1997) Gruszka, M., & Borysow, A. 1997, Icarus, 129, 172, doi: 10.1006/icar.1997.5773
  * Hartmann et al. (2011) Hartmann, J. M., Boulet, C., & Jacquemart, D. 2011, Journal of Chemical Physics, 134, 094316, doi: 10.1063/1.3557681
  * Hindmarsh et al. (2005) Hindmarsh, A. C., Brown, P. N., Grant, K. E., et al. 2005, ACM Trans. Math. Softw., 31, 363–396, doi: 10.1145/1089014.1089020
  * Hren et al. (2009) Hren, M. T., Tice, M. M., & Chamberlain, C. P. 2009, Nature, 462, 205, doi: 10.1038/nature08518
  * Hunten et al. (1984) Hunten, D. M., Tomasko, M. G., Flasar, F. M., et al. 1984, in Saturn, ed. T. Gehrels & N. Shapley Matthews (Univ. Arizona Press, Tucson), 671–759
  * Jennings et al. (2017) Jennings, D. E., Flasar, F. M., Kunde, V. G., et al. 2017, Applied Optics, 56, 5274, doi: 10.1364/ao.56.005274
  * Karman et al. (2015) Karman, T., Miliordos, E., Hunt, K. L. C., Groenenboom, G. C., & van der Avoird, A. 2015, Journal of Chemical Physics, 142, 084306, doi: 10.1063/1.4907917
  * Karman et al. (2019) Karman, T., Gordon, I. E., van der Avoird, A., et al. 2019, Icarus, 328, 160, doi: 10.1016/j.icarus.2019.02.034
  * Krissansen-Totton et al. (2018a) Krissansen-Totton, J., Arney, G. N., & Catling, D. C. 2018a, Proceedings of the National Academy of Science, 115, 4105, doi: 10.1073/pnas.1721296115
  * Krissansen-Totton et al. (2018b) Krissansen-Totton, J., Olson, S., & Catling, D. C. 2018b, Science Advances, 4, eaao5747, doi: 10.1126/sciadv.aao5747
  * Lax (2007) Lax, P. D. 2007, Linear Algebra and Its Applications, 2nd edn. (John Wiley & Sons Inc, Hoboken, New Jersey)
  * Lepage (1978) Lepage, G. P. 1978, Journal of Computational Physics, 27, 192, doi: 10.1016/0021-9991(78)90004-9
  * Lichtenberg (2021) Lichtenberg, T. 2021, Astrophysical Journal Letters, 914, L4, doi: 10.3847/2041-8213/ac0146
  * Lorenz et al. (1997) Lorenz, R. D., McKay, C. P., & Lunine, J. I. 1997, Science, 275, 642, doi: 10.1126/science.275.5300.642
  * Niemann et al. (2005) Niemann, H. B., Atreya, S. K., Bauer, S. J., et al. 2005, Nature, 438, 779, doi: 10.1038/nature04122
  * Niemann et al. (2010) Niemann, H. B., Atreya, S. K., Demick, J. E., et al. 2010, Journal of Geophysical Research (Planets), 115, E12006, doi: 10.1029/2010JE003659
  * Odintsova et al. (2021) Odintsova, T. A., Serov, E. A., Balashov, A. A., et al. 2021, Journal of Quantitative Spectroscopy and Radiative Transfer, 258, 107400, doi: 10.1016/j.jqsrt.2020.107400
  * Oparin et al. (2017) Oparin, D. V., Filippov, N. N., Grigoriev, I. M., & Kouzov, A. P. 2017, Journal of Quantitative Spectroscopy and Radiative Transfer, 196, 87, doi: 10.1016/j.jqsrt.2017.04.002
  * Pavlov et al. (2001) Pavlov, A. A., Brown, L. L., & Kasting, J. F. 2001, Journal of Geophysical Research, 106, 23267, doi: 10.1029/2000JE001448
  * Pavlov et al. (2000) Pavlov, A. A., Kasting, J. F., Brown, L. L., Rages, K. A., & Freedman, R. 2000, Journal of Geophysical Research, 105, 11981, doi: 10.1029/1999JE001134
  * Samuelson et al. (1997) Samuelson, R. E., Nath, N. R., & Borysow, A. 1997, Planetary and Space Science, 45, 959, doi: https://doi.org/10.1016/S0032-0633(97)00090-1
  * Schofield (1960) Schofield, P. 1960, Physical Review Letters, 4, 239, doi: 10.1103/PhysRevLett.4.239
  * Sung et al. (2016) Sung, K., Wishnow, E., Venkataraman, M., et al. 2016, in AAS/Division for Planetary Sciences Meeting Abstracts, Vol. 48, AAS/Division for Planetary Sciences Meeting Abstracts #48, 424.11
  * Tomasko et al. (2008) Tomasko, M. G., Bézard, B., Doose, L., et al. 2008, Planetary and Space Science, 56, 648, doi: 10.1016/j.pss.2007.10.012
  * Turbet et al. (2018) Turbet, M., Bolmont, E., Leconte, J., et al. 2018, Astronomy & Astrophysics, 612, A86, doi: 10.1051/0004-6361/201731620
  * van Kranendonk & Gass (1973) van Kranendonk, J., & Gass, D. M. 1973, Canadian Journal of Physics, 51, 2428, doi: 10.1139/p73-317
  * Wordsworth et al. (2017) Wordsworth, R., Kalugina, Y., Lokshtanov, S., et al. 2017, Geophysical Research Letters, 44, 665, doi: 10.1002/2016GL071766
  * Zahnle et al. (2006) Zahnle, K. J., Claire, M. W., & Catling, D. C. 2006, Geobiology, 4, 271, doi: https://doi.org/10.1111/j.1472-4669.2006.00085.x
  * Zahnle et al. (2019) Zahnle, K. J., Gacesa, M., & Catling, D. C. 2019, Geochimica et Cosmochimica Acta, 244, 56, doi: 10.1016/j.gca.2018.09.017


# Hybrid electroluminescent devices composed of (In,Ga)N micro-LEDs and
monolayers of transition metal dichalcogenides

K. Oreszczuk1, J. Slawinska2, A. Rodek1, M. Potemski1,3, C. Skierbiszewski2,
and P. Kossacki1 1 Institute of Experimental Physics, Faculty of Physics,
University of Warsaw, Pasteura 5, 02-093 Warsaw, Poland  2 Institute of High
Pressure Physics, Polish Academy of Sciences, Sokołowska 29/37, 01-142 Warsaw,
Poland  3 Laboratoire National des Champs Magnétiques Intenses, CNRS-UGA-UPS-
INSA-EMFL, 25 rue des Martyrs, 38042 Grenoble, France

###### Abstract

We demonstrate a novel electro-luminescence device in which GaN-based
$\upmu$-LEDs are used to trigger the emission spectra of monolayers of
transition metal dichalcogenides, which are deposited directly on $\upmu$-LED
surface. A special $\upmu$-LED design enables the operation of our structures
even in the limit of low temperatures. A device equipped with a selected WSe2
monolayer flake is shown to act as a stand-alone, electrically-driven single-
photon source.

##  I Introduction

Atomically thin semiconductors, such as monolayers of transition metal
dichalcogenides (TMDs) and their homo- and hetero-structures, exhibit
interesting optical properties. The list of prominent attributes of these two-
dimensional (2D) structures includes: strong excitonic bindingYou et al.
(2015); He et al. (2014); Ross et al. (2013) which favors room temperature
capabilitiesLorchat et al. (2018), possibilities to manipulate the valley
degree freedom Mak et al. (2010); Splendiani et al. (2010); Zeng et al.
(2012); Mak et al. (2012), and emergence of new class of quantum
emittersKoperski et al. (2015); Tonndorf et al. (2015); Brooks and Burkard
(2018); Kern et al. (2016); Kumar et al. (2015); Branny et al. (2017); Luo et
al. (2019). In parallel to the intense fundamental research, pertinent efforts
are also undertaken to exploit the TMD semiconductors in designing and
fabricating of new optoelectronic devices of possible practical useBlauth et
al. (2018); Peyskens et al. (2019). The operation of such devices is desired
to be electrically driven. Controlled doping of TMD materials, although
essential for the operation of conventional electroluminescent p-n diodes,
still remains a largely unresolved challenge Wang et al. (2018); Qiu et al.
(2013); Wang et al. (2016).

Therefore, alternative methods are being developed, relying most often on
tunneling mechanisms of carrier injection into the active 2D component of the
device. In fact, the electrically driven emission from TMD semiconductors has
been successfully demonstrated in a number of differently designed structures
Palacios-Berraquero et al. (2016); Sundaram et al. (2013); Ye et al. (2014);
Cheng et al. (2014); Withers et al. (2015); Wang et al. (2017); Li et al.
(2015). However, the architecture of such structures implies rather complex
manufacturing processes, and the performance of the reported devices remains
not yet fully satisfying. This calls for further efforts to improve/optimize
the already proposed device schemes and/or to search for optional solutions.

Taking advantage of the advanced technology of nitridesWu et al. (2018), we
put here forward a new concept of a compact pseudo-electroluminescent device
by integrating the TMD semiconductors into (In,Ga)N light-emitting diodes
(LEDs). Our devices are composed of micro-(In,Ga)N-LEDs ($\upmu$-LEDs), on top
of which we directly deposit the TMD layers: the electrically driven (In,Ga)N
$\upmu$-LEDs serve as the excitation sources to generate the PL emission from
TMDs. A special design of our (In,Ga)N diodes ensures their functionality at
cryogenic temperatures, which are preferential conditions for generating light
from TMD monolayers and an unavoidable requirement when tracing the quantum
emission centres from these layers. Low-temperature operation of our
electroluminescent devices comprised of different: MoS2, MoSe2, WSe2 and WS2
monolayers, is demonstrated. Importantly, the hybrid
(In,Ga)N-$\upmu$-LED/WSe2-monolayer device is also shown to function as an
electrically driven source of single photons. While the fabrication of the
present devices implies a laborious exfoliation manufacturing, we speculate
that the industrial scalability of such devices could be approached in the
future, taking into account the progress in the growth of TMDs by epitaxy
methods and thus evoking a possibility to grow the TMDs directly on (In,Ga)N
LEDs, particularly with molecular beam epitaxy techniquesPacuski et al.
(2020).

The technology of nitride LEDs has undergone intense development in recent
yearsWu et al. (2018). Relevant efforts have been made to respond to pertinent
demands for low energy consumption, that is, to achieve high luminous efficacy
of LEDs, but also to improve their brightness and contrast, their shock
resistance, and degradation time Wu et al. (2018). Although many nitride LEDs
already find relevant applications in display technology Day et al. (2011);
Chen et al. (2021), the new possibilities for nitride optoelectronics have
been raised more recently by employing a scheme of a tunnel junction (TJ)
Skierbiszewski et al. (2018a); Kowsz et al. (2017); Skierbiszewski et al.
(2018b). Such architecture of (In,Ga)N LEDs is also crucial for the present
work. The buried TJ enables the nitride LEDs to operate at cryogenic
temperaturesChlipala et al. (2020). Equally important is that the use of the
TJ allows us to replace the top p-type layer (of a conventional LED) with a
highly conductive n-type layer. This enhances the current spreading and
enables the n-type application of the side contacts to the device. Thus, TMD
flakes can be directly deposited on the GaN surface, and the experiments can
be carried out at low temperatures.

Figure 1:  (a) Schematic of the hybrid (In,Ga)N $\upmu$-LED/TMD-monolayer
device. (b) Sketch of the nitride $\upmu$-LED structure grown by plasma
assisted molecular beam epitaxy, with a characteristic sequence of epitaxial
sections. (c) Ion-implantation profile for the implemented He+ energy of $100$
keV and dose of $1.5\cdot 10^{14}\,$cm-2.

##  II Experimental details

The epitaxial structure, schematically shown in Figure 1, was grown by plasma-
assisted molecular beam epitaxy (PAMBE) on (0001) bulk GaN crystals. The PAMBE
grown TJ ($90$ nm) and $\upmu$-LED structure that emits light at a wavelength
of $450$ nm was capped with the n-type GaN of $150$ nm (Figure 1b). The
emission area of $\upmu$-LEDs was defined by a photolithography mask and He+
implantation ($100$ keV, $D=1.5\cdot 10^{14}\,$cm-2) of the $\upmu$-LED
structure. Ion energy and dose were adjusted to penetrate the region at the
location of the TJ (Figure 1c). The vacancies produced during the ion
implantation of the TJ deteriorate the p-type and n-type conductivity and
significantly increase the resistance of the TJ in the areas outside the
$\upmu$-LEDs. After the ion implantation, the regrowth of $200$ nm of high
conductive n-type GaN enhances the current spreading at the top of
$\upmu$-LEDs and allows the creation of the side contacts to the device. After
regrowth, the $150\times 300\,\mu$m devices containing arrays of individual
microdevices were separated by reactive-ion etching. Next, on the Ga-polar
side of the samples, the Ti/Al/Ni/Au metal contact was deposited with a
photolithographic mask, followed by lift-off. The same metal contacts were
used for bottom contact on the nitrogen polar side.

The TMD monolayers have been placed on top of (In,Ga)N $\upmu$-LEDs (see
Figure 1a) using deterministic transfer procedures. Monolayer flakes were
mechanically exfoliated from bulk crystals with chemically pure backgrinding
tape and transferred to the $\upmu$-LED surface using Gel-Pak DGL-X4
elastomeric films. A number of devices have been fabricated, comprising
different, MoS2, MoSe2, WSe2 and WS2 monolayers.

The optical measurements were carried out in the two experimental
configurations. Time-integrated studies were performed with the samples placed
in the cold-finger cryostat. The emission of TMD monolayers has been measured
in two different configurations. Primarily, the measurements have been carried
out using the external laser excitation at 532 nm ($2330$ meV). Secondly, the
luminescence was triggered from underneath with (In,Ga)N $\upmu$-LEDs.

The spectra composed of sharp emission lines characteristic of a WSe2
monolayer have also been examined with photon correlation experiments. These
experiments have been carried out in the Hanbury-Brown-Twiss configuration,
under (In,Ga)N LED excitation, and required the device to be immersed in
superfluid helium to minimize the undesirable effects of heating (when large
currents are driven to trigger the (In,Ga)N LEDs).

Figure 2:  (a-d) Images of the four different monolayers of TMDs. The flakes
are exfoliated on $\upmu$-LEDs, which span the whole photo area. The
$\upmu$-LED current density is equal to $0.3\,\mu$A$\,\mu$m-2. Red circles
mark the area from which the luminescence was acquired. (e-h) Spectra of the
monolayers acquired under $\upmu$-LED excitation (red) and $532\,$nm
($2330\,$meV) external laser excitation (blue). Dashed line: $\upmu$-LED
luminescence acquired at the off-flake spot on the sample.

##  III Results and discussion

The optical microscopy images of four different, WSe2, MoSe2, WS2 and MoS2
monolayers deposited on top of individual (In,Ga)N $\upmu$-LEDs are
illustrated in Figure 2a-d. The characteristic, low temperature ($10\,$K)
microluminescence spectra of these layers are shown, correspondingly, in
Figure 2e-h. The spectra shown with blue traces were measured under excitation
with the external laser ($532$ nm), whereas those shown in red have been
activated by the underneath (In,Ga)N $\upmu$–LEDs. The photoluminescence
signals measured under external excitation are typical (with respect to
spectral shape and characteristic energies and intensities) of TMD monolayers,
as often studied when they are directly deposited on conventional Si/SiO2
substrates. Characteristic spectra of our WSe2, MoSe2 and WS2 monolayers are
fairly well reproduced when they are excited with (In,Ga)N $\upmu$-LEDs. On
the other hand, the apparent background signal (see dashed traces in Figure
2e-f), which is due to an unfiltered low-energy tail of the (In,Ga)N emission,
largely distorts the spectrum of the MoS2 monolayer. This is, in fact, because
of a rather weak luminescence efficiency of MoS2 monolayers, as commonly
reported when this particular TMD monolayer is not protected by hexagonal
boron nitride.

Figure 3:  (a) Map of the luminescence of a long and narrow flake of a
monolayer WSe2. The signal is integrated over $1625$ meV - $1640$ meV range,
at the $\upmu$-LED current density of $0.5\,\mu$A$\,\mu$m-2. Blue circles mark
the outlines of the two circular $\upmu$-LEDs. Symbols A-D mark the spots
presented in panels (b-c). (b) $\upmu$-LED excited PL of the homogeneous part
of the WSe2 flake (black) and the background luminescence of the $\upmu$-LED
(grey). (c) $\upmu$-LED excited PL from the two spots with the narrow-line
emitters. (d) The dependence of the luminescence intensity of line D on the
$\upmu$-LED current density. The excitation sublinearity is caused by the
elevated sample surface temperature at the highest current densities.

Having established the overall functionality of our pseudo-electroluminescence
devices, we now turn our attention to structures displaying sharp emission
lines, possibly comprising the quantum emitter centers. Such centers
preferentially appear in WSe2 layers. Several
$\upmu$-(In,Ga)N-LED/WSe2-monolayer devices have been fabricated and tested.
The irregular shape of WSe2 flakes has been found to favour the observation of
sharp emission lines. One may speculate this is due to the increased
probability of the formation of the specific centers trapping the photoexcited
carries at the edges of TMD flakesKoperski et al. (2015); Tonndorf et al.
(2015) or induced by inhomogenous strain distributionBrooks and Burkard
(2018); Kern et al. (2016); Kumar et al. (2015); Branny et al. (2017); Luo et
al. (2019).

The representative example of a device displaying sharp emission lines
contains a narrow, elongated flake of a WSe2 monolayer. The device has been
mounted on a cold finger of the cryostat, cooled down to $10$ K and primarily
tested under such conditions. After turning on the (In,Ga)N diodes, the
emission spectra, in the spectral range $1590-1775$ meV, were measured when
scanning the surface of the WSe2 flake and its surrounding, with a spatial
resolution of $2\,\mu$m. A map of the intensity of the PL integrated over the
spectral range of $1625-1640$ meV is shown in Figure 3a. This image
illustrates well the shape of the flake, which extends over the active areas
of two neighbouring (In,Ga)N $\upmu$-LEDs. The background signal due to the
emission tail of the (In,Ga)N diode is fairly weak - see the spectrum measured
outside the flake area (Figure 3b). The spectra measured in the areas in the
middle of the flake, away from its edges, exhibit a broad feature, centered at
$1660$ meV – a typical response of WSe2 ML with significant spatial
inhomogeneity (see Figure 3b). However, near the edges and narrow parts of the
flake, the luminescence spectra often display sharp lines; see Figure 3c.

Notably, the low activation current (below $1\,$nA$\,\mu$m-2) of our device
permits the measurements of the emission spectra of the WSe2 monolayer over a
broad range of current densities driving the (In,Ga)N diodes and thus enables
the straightforward tuning of the excitation power to match the saturation
point of a particular emitter. To illustrate this effect, we present a typical
PL intensity of the emitter as a function of the current density (Figure 3d).
As the luminescence intensity of narrow lines in WSe2 monolayers is known to
drop down upon increase of temperatureOreszczuk et al. (2020) we conclude that
the observed effect is likely due to a low cooling efficiency of our device
under the current arrangement (limited by the heat transfer effectiveness
between our device and the mounting base in the cryostat). However, the heat
generation process can be suppressed by reducing the size of the $\upmu$-LED
and the precise positioning of the TMD flakes.

To demonstrate the operation of our device as a single-photon source, the
device was immersed in a superfluid helium bath and cooled down to a
temperature of $1.65$ K. This resulted in an improved cooling efficiency of
the device, enabling the application of higher diode currents, and at the same
time gaining higher emission intensities of narrow emission lines. This
allowed for the correlation measurements of the photon statistics emitted from
the particular centre.

In Figure 4a we present the spectrum of the narrow line emitter in the WSe2
monolayer under the (In,Ga)N $\upmu$-LED excitation that was powered with
$1.8\,\mu$A$\,\mu$m-2 current. We note the particularly low ratio of the
background $\upmu$-LED spectral tail to the narrow emission line intensity.
Any uncorrelated photons detected in the Hanbury-Brown-Twiss experiment would
proportionally reduce the amplitude of coincidence dip as well as the
feasibility of our device as a single-photon source.

Figure 4:  (a) Luminescence spectrum of the narrow emission line chosen for
the single-photon correlation experiment. (b) Coincidence correlation counts
as a function of the time interval between photons. The photoluminescence is
excited with $\upmu$-LED. The broad background luminescence from the
$\upmu$-LED was not subtracted in the post-processing.

The acquired photon coincidence correlation function, $g_{2}(\tau)$, is
presented in Figure 4b. The photon coincidence correlation function
$g_{2}(\tau)$ denotes the relative probability density of an event with a time
separation $\tau$. The decrease in the $g_{2}$ value around $\tau=0$ evidences
the single-photon component in the emission. The emission has predominantly
single-photon characteristics when $g_{2}(0)<0.5$.

The amplitude of the correlation function at zero delay is significantly low
($g_{2}(0)=0.22$), which attests to the single-photon character of the
emission. The photon coincidence correlation can also be used to determine the
lifetime of the emitter state (Figure 4b). The lifetime is equal to the
characteristic decay time of the antibunching feature in the $g_{2}(\tau)$
profile at the low photon rate limit. The lifetime obtained in our experiment
($t_{decay}=4.8\,$ns) is in agreement with the results obtained in other
studies Srivastava et al. (2015); Kumar et al. (2015); He et al. (2015);
Koperski et al. (2015); Chakraborty et al. (2015).

We also used the photon coincidence correlation experiment to verify the
temporal stability of the emitter. Fluctuations of either the emitter
intensity or the intensity of its excitation source would result in the effect
of photon bunching. We observe no decrease of $g_{2}$ at large $\tau$,
confirming the good stability of the emission center and the $\upmu$-LED. We
also observed good long-term stability of our device. During several days of
measurements, we have observed no significant variation in the emission
efficiency of the selected center.

##  IV Conclusion

In conclusion, we have demonstrated that light emission from TMD monolayers
can be efficiently generated when they are directly deposited on, and excited
by the (In,Ga)N $\upmu$-LEDs. Our hybrid devices, $100\,\mu$m in size, can be
seen as an optional route towards the effective, electrically driven
generation of light from TMD semiconductors. The devices equipped with
different MoS2, WSe2, WS2 and MoSe2 monolayers have been shown to be
operational at low temperatures thanks to the unique design of our (In,Ga)N
micro-LEDs. A proof of concept for a new, electrically driven source of single
photons has been demonstrated with one of our (In,Ga)N
$\upmu$-LEDs/WSe2-monolayer devices.

## Acknowledgements

This work is supported by TEAM-TECH POIR.04.04.00-00-210C/16-00 and TEAM
POIR.04.04.00-00-1A18/16-01 (ATOMOPTO) projects of the Foundation for Polish
Science co-financed by the European Union under the European Regional
Development Fund. This work was supported under project
2021/41/N/ST3/04240 by National Science Centre, Poland. A.R. acknowledges
support from the "Diamentowy Grant" under decision DI2017/008347 of MEiN of
Poland. The work has been also supported by EU Graphene Flagship project.

## References

  * You et al. (2015) Y. You, X.-X. Zhang, T. C. Berkelbach, M. S. Hybertsen, D. R. Reichman, and T. F. Heinz, Nature Physics 11, 477 (2015).
  * He et al. (2014) K. He, N. Kumar, L. Zhao, Z. Wang, K. F. Mak, H. Zhao, and J. Shan, Physical review letters 113, 026803 (2014).
  * Ross et al. (2013) J. S. Ross, S. Wu, H. Yu, N. J. Ghimire, A. M. Jones, G. Aivazian, J. Yan, D. G. Mandrus, D. Xiao, W. Yao, et al., Nature communications 4, 1474 (2013).
  * Lorchat et al. (2018) E. Lorchat, S. Azzini, T. Chervy, T. Taniguchi, K. Watanabe, T. W. Ebbesen, C. Genet, and S. Berciaud, ACS Photonics 5, 5047 (2018).
  * Mak et al. (2010) K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, Physical review letters 105, 136805 (2010).
  * Splendiani et al. (2010) A. Splendiani, L. Sun, Y. Zhang, T. Li, J. Kim, C.-Y. Chim, G. Galli, and F. Wang, Nano letters 10, 1271 (2010).
  * Zeng et al. (2012) H. Zeng, J. Dai, W. Yao, D. Xiao, and X. Cui, Nature nanotechnology 7, 490 (2012).
  * Mak et al. (2012) K. F. Mak, K. He, J. Shan, and T. F. Heinz, Nature nanotechnology 7, 494 (2012).
  * Koperski et al. (2015) M. Koperski, K. Nogajewski, A. Arora, V. Cherkez, P. Mallet, J.-Y. Veuillen, J. Marcus, P. Kossacki, and M. Potemski, Nature nanotechnology 10, 503 (2015).
  * Tonndorf et al. (2015) P. Tonndorf, R. Schmidt, R. Schneider, J. Kern, M. Buscema, G. A. Steele, A. Castellanos-Gomez, H. S. van der Zant, S. M. de Vasconcellos, and R. Bratschitsch, Optica 2, 347 (2015).
  * Brooks and Burkard (2018) M. Brooks and G. Burkard, Physical Review B 97, 195454 (2018).
  * Kern et al. (2016) J. Kern, I. Niehues, P. Tonndorf, R. Schmidt, D. Wigger, R. Schneider, T. Stiehm, S. Michaelis de Vasconcellos, D. E. Reiter, T. Kuhn, et al., Advanced Materials 28, 7101 (2016).
  * Kumar et al. (2015) S. Kumar, A. Kaczmarczyk, and B. D. Gerardot, Nano letters 15, 7567 (2015).
  * Branny et al. (2017) A. Branny, S. Kumar, R. Proux, and B. D. Gerardot, Nature communications 8, 15053 (2017).
  * Luo et al. (2019) Y. Luo, N. Liu, X. Li, J. C. Hone, and S. Strauf, 2D Materials 6, 035017 (2019).
  * Blauth et al. (2018) M. Blauth, M. Jürgensen, G. Vest, O. Hartwig, M. Prechtl, J. Cerne, J. J. Finley, and M. Kaniber, Nano Letters 18, 6812 (2018), ISSN 1530-6984.
  * Peyskens et al. (2019) F. Peyskens, C. Chakraborty, M. Muneeb, D. Van Thourhout, and D. Englund, Nature Communications 10, 4435 (2019), ISSN 2041-1723.
  * Wang et al. (2018) J. Wang, I. Verzhbitskiy, and G. Eda, Advanced Materials 30, 1802687 (2018), ISSN 1521-4095.
  * Qiu et al. (2013) H. Qiu, T. Xu, Z. Wang, W. Ren, H. Nan, Z. Ni, Q. Chen, S. Yuan, F. Miao, F. Song, et al., Nature Communications 4, 2642 (2013), ISSN 2041-1723.
  * Wang et al. (2016) S. Wang, W. Zhao, F. Giustiniano, and G. Eda, Physical Chemistry Chemical Physics 18, 4304 (2016), ISSN 1463-9084.
  * Palacios-Berraquero et al. (2016) C. Palacios-Berraquero, M. Barbone, D. M. Kara, X. Chen, I. Goykhman, D. Yoon, A. K. Ott, J. Beitner, K. Watanabe, T. Taniguchi, et al., Nature Communications 7, 12978 (2016), ISSN 2041-1723.
  * Sundaram et al. (2013) R. S. Sundaram, M. Engel, A. Lombardo, R. Krupke, A. C. Ferrari, P. Avouris, and M. Steiner, Nano Letters 13, 1416 (2013), ISSN 1530-6984.
  * Ye et al. (2014) Y. Ye, Z. Ye, M. Gharghi, H. Zhu, M. Zhao, Y. Wang, X. Yin, and X. Zhang, Applied Physics Letters 104, 193508 (2014), ISSN 0003-6951.
  * Cheng et al. (2014) R. Cheng, D. Li, H. Zhou, C. Wang, A. Yin, S. Jiang, Y. Liu, Y. Chen, Y. Huang, and X. Duan, Nano Letters 14, 5590 (2014), ISSN 1530-6984.
  * Withers et al. (2015) F. Withers, O. Del Pozo-Zamudio, A. Mishchenko, A. P. Rooney, A. Gholinia, K. Watanabe, T. Taniguchi, S. J. Haigh, A. K. Geim, A. I. Tartakovskii, et al., Nature Materials 14, 301 (2015), ISSN 1476-4660.
  * Wang et al. (2017) S. Wang, J. Wang, W. Zhao, F. Giustiniano, L. Chu, I. Verzhbitskiy, J. Zhou Yong, and G. Eda, Nano Letters 17, 5156 (2017), ISSN 1530-6984.
  * Li et al. (2015) D. Li, R. Cheng, H. Zhou, C. Wang, A. Yin, Y. Chen, N. O. Weiss, Y. Huang, and X. Duan, Nature Communications 6, 7509 (2015), ISSN 2041-1723.
  * Wu et al. (2018) T. Wu, C.-W. Sher, Y. Lin, C.-F. Lee, S. Liang, Y. Lu, S.-W. Huang Chen, W. Guo, H.-C. Kuo, and Z. Chen, Applied Sciences 8, 1557 (2018).
  * Pacuski et al. (2020) W. Pacuski, M. Grzeszczyk, K. Nogajewski, A. Bogucki, K. Oreszczuk, J. Kucharek, K. E. Połczyńska, B. Seredyński, A. Rodek, R. Bożek, et al., Nano Letters 20, 3058 (2020), ISSN 1530-6984.
  * Day et al. (2011) J. Day, J. Li, D. Y. C. Lie, C. Bradford, J. Y. Lin, and H. X. Jiang, Applied Physics Letters 99, 031116 (2011), ISSN 0003-6951.
  * Chen et al. (2021) Z. Chen, S. Yan, and C. Danesh, Journal of Physics D: Applied Physics 54, 123001 (2021), ISSN 0022-3727.
  * Skierbiszewski et al. (2018a) C. Skierbiszewski, H. Turski, M. Zak, K. Nowakowski-Szkudlarek, G. Muziol, M. Siekacz, A. Feduniewicz-Zmuda, and M. Sawicka, in _2018 76th Device Research Conference (DRC)_ (2018a), pp. 1–2.
  * Kowsz et al. (2017) S. J. Kowsz, E. C. Young, B. P. Yonkee, C. D. Pynn, R. M. Farrell, J. S. Speck, S. P. DenBaars, and S. Nakamura, Optics Express 25, 3841 (2017), ISSN 1094-4087.
  * Skierbiszewski et al. (2018b) C. Skierbiszewski, G. Muziol, K. Nowakowski-Szkudlarek, H. Turski, M. Siekacz, A. Feduniewicz-Zmuda, A. Nowakowska-Szkudlarek, M. Sawicka, and P. Perlin, Applied Physics Express 11, 034103 (2018b), ISSN 1882-0786.
  * Chlipala et al. (2020) M. Chlipala, H. Turski, M. Siekacz, K. Pieniak, K. Nowakowski-Szkudlarek, T. Suski, and C. Skierbiszewski, Optics Express 28, 30299 (2020), ISSN 1094-4087.
  * Oreszczuk et al. (2020) K. Oreszczuk, T. Kazimierczuk, T. Smoleńński, K. Nogajewski, M. Grzeszczyk, A. Łopion, M. Potemski, and P. Kossacki, Physical Review B 102, 245409 (2020).
  * Srivastava et al. (2015) A. Srivastava, M. Sidler, A. V. Allain, D. S. Lembke, A. Kis, and A. Imamoğlu, Nature nanotechnology 10, 491 (2015).
  * He et al. (2015) Y.-M. He, G. Clark, J. R. Schaibley, Y. He, M.-C. Chen, Y.-J. Wei, X. Ding, Q. Zhang, W. Yao, X. Xu, et al., Nature nanotechnology 10, 497 (2015).
  * Chakraborty et al. (2015) C. Chakraborty, L. Kinnischtzke, K. M. Goodfellow, R. Beams, and A. N. Vamivakas, Nature nanotechnology 10, 507 (2015).