# Spin Effect induced Momentum Spiral and Asymmetry Degree in Pair Production

Li-Na Hu  Key Laboratory of Beam Technology of the Ministry of Education, and
College of Nuclear Science and Technology, Beijing Normal University, Beijing
100875, China  Hong-Hao Fan  Key Laboratory of Beam Technology of the Ministry
of Education, and College of Nuclear Science and Technology, Beijing Normal
University, Beijing 100875, China  Orkash Amat  Key Laboratory of Beam
Technology of the Ministry of Education, and College of Nuclear Science and
Technology, Beijing Normal University, Beijing 100875, China  Suo Tang
College of Physics and Optoelectronic Engineering, Ocean University of China,
Qingdao, Shandong 266100, China  Bai-Song Xie<EMAIL_ADDRESS>Key
Laboratory of Beam Technology of the Ministry of Education, and College of
Nuclear Science and Technology, Beijing Normal University, Beijing 100875,
China  Institute of Radiation Technology, Beijing Academy of Science and
Technology, Beijing 100875, China

###### Abstract

Spin effect on the pair production under circularly polarized fields are
investigated. Significantly different from what momentum spirals caused by two
counter-rotating fields with a time delay, we find for the first time that the
spirals can also be induced due to the particles spin effect even if in a
single field. We further examine the bichromatic combinational fields, the
inhomogeneous spiral structures can be observed in the momentum spectrum, in
particular, the spiral not only does exist in two cases of spin but also is
about two orders of magnitude amplifier than that in the single field.
Meanwhile, the spin asymmetry degree on the momentum distributions is
investigated and found that there exist the effect of spin flip with
increasing time delay between two fields. The spin asymmetry degree on the
number density can reach to $98\%$ in a certain of condition. These results
indicate that the signatures of created particles, especially the spiral
structures are strongly associated with the information of laser field as well
as the created particle spin, which can deepen the understanding of vacuum
pair production.

###### pacs:

12.20.Ds, 11.15.Tk

Introduction.–In the presence of strong background field, the quantum
electrodynamic (QED) vacuum is unstable and decays into electron-positron
pairs Sauter:1931zz ; Heisenberg:1935qt ; DiPiazza:2011tq ; Fedotov:2022 ;
Xie:2017 ; Dunne:2008kc ; Schwinger:1951nm ; Burke:1997ew . A great deal of
research has focused on the momentum distribution and number density of
created pairs for different external fields. The momentum distribution
contains rich features such as effective mass signatures Kohlfurst:2013ura ,
the self-bunching effect Hebenstreit:2011wk , the ponderomotive force effect
Kohlfurst:2017hbd , the multiple-slit interference effect Li:2014En , node
structures Li:2015Non , as well as recently discovered spiral structures
Li:2017qwd ; Hu:2023pmz ; Li:2019hzi .

In fact, spiral structures have been widely investigated in many areas in past
years, for example, in atomic and molecular ionization Ngoko:2015Ele ;
Pengel:2017Elec ; Majczak:2022 ; Macek2009 ; K2016 ; L2020 , nonlinear optics
Harris1994 , type-II superconductors Blatter1994 , plasmas physics Uby1995 ;
Shukla2009 , atomic condensates Madison2000 , and so on. Interestingly, recent
investigations indicate that the momentum spirals are discovered in electron-
positron pair production Li:2017qwd ; Hu:2023pmz ; Li:2019hzi , which are
generated by two counter-rotating fields with a time delay. It should be noted
that the spirals with even and odd arms are attributed to the fields by the
same color Li:2017qwd ; Hu:2023pmz  and by the different colors Li:2019hzi ,
respectively.

On the other hand, it is known that the spin effect of created particles plays
a critical role in vacuum pair production Strobel:2014tha ; Blinne:2015zpa ;
Li:2019 ; Kohlf2019S . For example, Strobel et al. proved that the spin
distribution of produced pairs is generally not $1:1$ in the Schwinger pair
production by the rotating electric fields depending on time Strobel:2014tha .
Blinne et al. pointed out that the pair production in different spin states is
unequal for constant rotating field Blinne:2015zpa . Interestingly, we found
for circularly polarization (CP) field that the number density of created
bosons is the geometric mean, i.e., the square root of spin-up number density
times that of spin-down Li:2019 , meanwhile, for arbitrarily polarized
multicycle fields, that the degree of spin polarization roughly increases with
$\gamma\ll 1$, where $\gamma=m\omega/eE_{0}$ is Keldysh adiabatic parameter,
$m$ and $-e$ denote electron mass and charge, respectively. Recently,
Kohlfürst showed the power of the photon absorption model by demonstrating how
to identify the imprint the different spin states leave on the electron-
positron angular distribution, and obtained the spin-dependent particle
production amplitudes Kohlf2019S .

To our knowledge, however, whether the spin of created particles can create
momentum spiral has not yet been reported. Therefore, as an important concern,
we shall investigate the problem in this Letter. The focus is on considering
the spin induced momentum spiral by the single field, and the spiral structure
as well as the spin asymmetry degree under the bichromatic combination fields.
For the single field, it is found that there is no momentum spiral when the
spin of created particles is parallel to the direction of laser field
($s=+1$). Surprisingly, when the spin is antiparallel to the field direction
($s=-1$), there exists an obvious spiral structure in the momentum spectrum.
And the spiral in the subcycle field is more remarkable than that in the
multicycle field. On the other hand, for the bichromatic combination fields,
the spiral does exist no matter what is for the $s=+1$ or/and $s=-1$. Not only
the unobserved spiral in a single field is seen in the two fields, but also
the observed spiral in a single field is amplified about two orders of
magnitude larger in the bichromatic fields. Meanwhile, it is found that the
larger the time delay between the two fields, the more pronounced the spiral
structure. Finally, by adjusting appropriate field parameters, we can also
obtain high spin asymmetry degree of created particles.

Our considered electric field form is associated with either of one single
field or their combination by two fields with a little different strength and
frequency between them. For single field, it is
$\mathbf{E}(t)=\mathbf{E}_{1}(t)=\frac{E_{1}}{\cosh(\frac{t}{\tau_{1}})}[\cos(\omega_{1}t)\mathbf{e}_{x}+\delta_{1}\sin(\omega_{1}t)\mathbf{e}_{y}]$
or
$\mathbf{E}(t)=\mathbf{E}_{2}(t)=\frac{E_{2}}{\cosh(\frac{t-T_{d}}{\tau_{2}})}[\cos(\omega_{2}(t-T_{d}))\mathbf{e}_{x}+\delta_{2}\sin(\omega_{2}(t-T_{d}))\mathbf{e}_{y}]$.
And for the combined field case, it is
$\mathbf{E}(t)=\mathbf{E}_{1}(t)+\mathbf{E}_{2}(t)$, where
$E_{1,2}=E_{01,02}/\sqrt{1+\delta_{1,2}^{2}}$ are the field strengths
($|\delta_{1,2}|=1$ are the CP degree), $\omega_{1,2}$ are the field
frequencies, $\tau_{1,2}=N2\pi/\omega_{1,2}$ denote the pulse duration ($N$
determines the number of cycles in the individual pulse), $T_{d}$ is the time
delay between the two pulses (it is in unit of $T=\tau_{1}+\tau_{2}$). In our
study, the typical field parameters are set as $E_{01}=0.1\sqrt{2}E_{cr}$,
$E_{02}=0.07\sqrt{2}E_{cr}$, $\delta_{1}=-1$, $\delta_{2}=1$,
$\omega_{1}=0.44$, $\omega_{2}=0.55$. Note that $\hslash=c=1$ are used and all
quantities are presented in terms of the electron mass $m$.

We employ the equal-time Dirac-Heisenberg-Wigner (DHW) formalism
Strobel:2014tha ; Blinne:2015zpa ; Kohlfurst:2015zxi ; Hebenstreit:2011pm ;
Bialynicki-Birula:1991jwl ; Hebenstreit:2010vz , a relativistic phase-space
quantum kinetic approach, which has been widely applied to study pair
production in arbitrary electromagnetic fields Ababekri:2019dkl ;
Kohlf2020Effect ; Blinne:2013via ; Li:2021vjf ; Hebenstreit:2011wk ;
Olugh:2018seh ; Li:2015cea ; Kohlfurst:2017git ; Olugh:2019nej . For the
spatially homogeneous and time-dependent electric field presented in our
study, the complete DHW equations of motion can be reduced to the following
form Blinne:2015zpa ; Li:2017qwd ; Olugh:2019nej

$\begin{split}\dot{f}&=\frac{e\mathbf{E}\cdot\mathbf{v}}{2\Omega},\\\
\dot{\mathbf{v}}&=\frac{2}{\Omega^{3}}[(e\mathbf{E}\cdot\mathbf{p})\mathbf{p}-e\mathbf{E}\Omega^{2}](f-1)-\frac{(e\mathbf{E}\cdot\mathbf{v})\mathbf{p}}{\Omega^{2}}\\\
&-2\mathbf{p}\times\mathbbm{a}-2m\mathbbm{t},\\\
\dot{\mathbbm{a}}&=-2\mathbf{p}\times\mathbf{v},\\\
\dot{\mathbbm{t}}&=\frac{2}{m}[m^{2}\mathbf{v}+(\mathbf{p}\cdot\mathbf{v})\mathbf{p}],\end{split}$
(1)

with initial conditions $f(\mathbf{q},-\infty)=0$,
$\mathbf{v}(\mathbf{q},-\infty)=\mathbbm{a}(\mathbf{q},-\infty)=\mathbbm{t}(\mathbf{q},-\infty)=0$,
where $f$ is the single-particle momentum distribution function, $\mathbf{v}$
is associated with the current density $\mathbbm{v}$, $\mathbbm{a}$ denotes
spin density, and the projection of the $\mathbbm{t}$ in the direction of
momentum determines the rate of change of the mass density Kohlfurst:2015zxi ;
Hebenstreit:2011pm ; Kohlf2019S .
$\Omega=\sqrt{m^{2}+(\mathbf{q}-e\mathbf{A}(t))^{2}}$ represents the total
energy of particles, $\mathbf{q}$ is the canonical momentum related to the
kinetic momentum $\mathbf{p}(t)=\mathbf{q}-e\mathbf{A}(t)$, $\mathbf{A}(t)$ is
the vector potential of the electric field $\mathbf{E}(t)$. By integrating the
$f(\mathbf{q},t)$ over full momenta at $t\rightarrow+\infty$, we can obtain
the number density of created particles:

$n=\lim_{t\to+\infty}\int\frac{d^{3}q}{(2\pi)^{3}}f(\mathbf{q},t)\,.$ (2)

Note that the above $f$ and $n$ are the forms without considering the spin of
created particles. When considering the spin, the single-particle momentum
distribution function can be changed as Blinne:2015zpa ; Li:2019

$f_{s}=\frac{1}{2}\left(f+s\,\delta f_{\mathrm{sc}}\right),$ (3)

where $s=\pm 1$, $\delta f_{\mathrm{sc}}=\frac{q_{z}}{\epsilon_{\perp}}\delta
f_{\mathrm{c}}+\frac{m}{\epsilon_{\perp}}\delta f_{\mu_{z}}\,$ is associated
with the chiral asymmetry $\delta
f_{\mathrm{c}}=\frac{1}{2\Omega}\mathbf{p}\cdot\mathbbm{a}$ and the magnetic
moment asymmetry $\delta
f_{\mathrm{\mu_{z}}}=\frac{1}{2\Omega}\left(ma_{z}+(\mathbf{p}\times\mathbbm{t}\,)_{z}\right)$,
$\epsilon_{\perp}=\sqrt{m^{2}+{q_{z}}^{2}}$. Also the corresponding number
density of considering the spin could be modified as
$n_{s}=\lim_{t\to+\infty}\int\frac{d^{3}q}{(2\pi)^{3}}f_{s}(\mathbf{q},t)\,$.

Figure 1: (color online). Momentum spectra of created particles in the
($q_{x}$, $q_{y}$) plane (where $q_{z}=0$) under the single field
$\mathbf{E}_{1}(t)$ (for the first column) and the bichromatic combination
fields for $T_{d}=8T$ (for the second column) with $N=0.5$. From up to down,
each row corresponds to $f(q_{x},q_{y})$, $f_{+}(q_{x},q_{y})$ and
$f_{-}(q_{x},q_{y})$, respectively. Other parameters are
$E_{01}=0.1\sqrt{2}E_{\mathrm{cr}}$, $\delta_{1}=-1$, $\omega_{1}=0.44$,
$E_{02}=0.07\sqrt{2}E_{\mathrm{cr}}$, $\delta_{2}=1$, $\omega_{2}=0.55$.

Momentum distribution.–The momentum spectra of created particles in the single
field $\mathbf{E}_{1}(t)$ and the bichromatic combination fields with $N=0.5$
are shown in Fig. 1. It is found that under the $\mathbf{E}_{1}(t)$, for
$f(q_{x},q_{y})$, there is no obvious spiral structure but the distribution
presents a good symmetry in the $q_{x}$ direction, see Fig. 1(a). The symmetry
is related to the fact that the electric field in the $x$ direction is an even
function Li:2015Non ; Hu:2023pmz . Meanwhile, one can see that the
$f(q_{x},q_{y})$ is mainly distributed in the lower half plane, where it
dominates the contribution to the pair production, which is consistent with
the results in terms of the turning points by $\Omega(\mathbf{q},t)=0$
Dumlu1011 ; Strobel2015 . For $f_{\pm}(q_{x},q_{y})$, however, the momentum
distributions present a clear difference. We find that the
$f_{+}(q_{x},q_{y})$ shown in Fig. 1(b) is dominant. But surprisingly, there
exists a remarkable spiral structure in the $f_{-}(q_{x},q_{y})$, see Fig.
1(c), although it is about two orders of magnitude weaker than the
$f_{+}(q_{x},q_{y})$. Actually, we do the results of $\mathbf{E}_{2}(t)$
(patterns are not shown), which are similar to the results of
$\mathbf{E_{1}}(t)$, except that the momentum distribution shifts from the
lower half plane to the upper half plane, and the spiral structures appears in
the $f_{+}(q_{x},q_{y})$. These are not surprising since that the spin
directions of the fields $\mathbf{E}_{1}(t)$ and $\mathbf{E}_{2}(t)$ are
opposite with each other.

Note that when we consider the spin of created particles, there are two ways
for spin alignment in a pair of spin-1/2 particles: parallel and antiparallel
Kohlf2019S . For the parallel, the total spin is $S=1/2+1/2=1$, while for the
antiparallel, the total spin is $S=0$. Obviously, it is found that the
distribution of created pairs for $S=1$ is larger than that for $S=0$, which
results in agreement with that in Ref.Kohlf2019S .

Now let us discuss briefly the possible relation for the distribution region,
for example it is mainly focused in the lower half plane when field
$\mathbf{E}_{1}(t)$ is applied, to the conservation of total angular momentum
in pair production process. Note that the total angular momentum is composed
of orbital angular momentum (OAM) and spin angular momentum, i.e., $J=L+S$. It
is known that the electron-positron pairs are mainly produced in about a
field-period near the $t=0$ since that the field strength is larger than those
in other time regime. For one pair of electron-positron, the OAM produced in
this period has a positive trend fortunately when the momentum is located in
lower half plane. Therefore, in the probability sense, we can think the OAM is
positive, which is indeed required for guaranteeing the conservation of total
angular momentum in pair production precess. For example, if one-pair
electron-positron is produced via absorbing at least $5$ photons (since the
original central frequency is $\omega_{1}=0.44$), then the OAM of
electron/positron should be $L=2$ when $S=1$ so that the total angular
momentum is $J=5$. Certainly exact descriptions and examinations about the
total angular momentum conservation problem need to adopt the more powerful
theoretical tool such as the S-matrix treatment.

It should be emphasized that the momentum spirals in a single field in our
study mentioned above are induced by the spin effect of created particles,
which is different from the previous one Li:2017qwd ; Hu:2023pmz ; Li:2019hzi
that is caused by the opposite photons’ spin in two counter-rotating fields
with a time delay.

Under the bichromatic combination fields with $T_{d}=8T$, it is found that
there are the remarkable spiral pattern in the momentum spectra. For
$f(q_{x},q_{y})$, the spiral structure is constituted of nine arms and its
distribution in the whole plane is nonuniform with a stronger in lower half
plane, see Fig. 1(d). This phenomenon can be understood in connection with the
results of a single field. It is known that the momentum spectra
$f(q_{x},q_{y})$ of $\mathbf{E}_{1}(t)$ and $\mathbf{E}_{2}(t)$ are mainly
distributed in the lower half plane and the upper half plane, respectively.
And due to the field strength of $\mathbf{E}_{1}(t)$ is larger than that of
$\mathbf{E}_{2}(t)$, so the $f(q_{x},q_{y})$ of $\mathbf{E}_{1}(t)$ is
stronger than ones of $\mathbf{E}_{2}(t)$. When time delay is large, the two
fields $\mathbf{E}_{1}(t)$ and $\mathbf{E}_{2}(t)$ are almost independent,
which leads to the result under the combined fields being almost a
superposition of the ones under the two individual fields. Therefore, for
$T_{d}=8T$, the $f(q_{x},q_{y})$ of the lower half plane is stronger than that
of the upper half plane.

For $f_{\pm}(q_{x},q_{y})$, one can see that the spiral of
$f_{+}(q_{x},q_{y})$ is mainly distributed in the lower half plane, as shown
in Fig. 1(e), while the result of $f_{-}(q_{x},q_{y})$ is opposite, see Fig.
1(f). The results are also associated with the momentum distributions in a
single field. Specifically, let’s first understand the result of
$f_{+}(q_{x},q_{y})$ in Fig. 1(e). In the combined fields with larger time
delay, we can think that the $f_{+}(q_{x},q_{y})$ is consisted of the
$f_{+}(q_{x},q_{y})$ of $\mathbf{E}_{1}(t)$ and $\mathbf{E}_{2}(t)$. It is
known that the $f_{+}(q_{x},q_{y})$ of $\mathbf{E}_{1}(t)$ and
$\mathbf{E}_{2}(t)$ are predominantly distributed in the lower half-plane and
the upper half-plane, respectively, and there exists an obvious spiral
structure in the $f_{+}(q_{x},q_{y})$ of $\mathbf{E}_{2}(t)$. Meanwhile, the
$f_{+}(q_{x},q_{y})$ of $\mathbf{E}_{1}(t)$ is about two orders of magnitude
larger than that of $\mathbf{E}_{2}(t)$. When the $\mathbf{E}_{1}(t)$ and
$\mathbf{E}_{2}(t)$ are combined, there should be $f_{+}(q_{x},q_{y})$ in the
whole plane. But since the $f_{+}(q_{x},q_{y})$ of $\mathbf{E}_{1}(t)$ mainly
contributes to the position of distribution, while the $f_{+}(q_{x},q_{y})$ of
$\mathbf{E}_{2}(t)$ mainly contributes to the spiral structure. Therefore,
these results ultimately lead to the $f_{+}(q_{x},q_{y})$ in the combined
fields mainly distributed in the lower half plane and accompanied by a
remarkable spiral structure, where the spiral is about two orders of magnitude
amplifier than that in the single field. For the results of
$f_{-}(q_{x},q_{y})$ in the combined fields, a similar understanding can be
performed.

Figure 2: (color online). The same as in Fig. 1 except that the distributions
are in the ($q_{x}$, $q_{z}$) plane (where $q_{y}=0$). Figure 3: (color
online). Momentum spectra $f_{-}(q_{x},q_{y})$ of created particles in the
($q_{x}$, $q_{y}$) plane (where $q_{z}=0$) for the single field
$\mathbf{E}_{1}(t)$ with different $N$. (a) and (b) correspond to $N=1$ and
$N=2$, respectively. (c) is the amplified result of (b).

Now let us to see the momentum distributions in the ($q_{x}$, $q_{z}$) plane
(where $q_{y}=0$) under the single or/and the bichromatic combined fields with
$N=0.5$, the results are shown in Fig. 2. Due to the distributions are
symmetric about both of $q_{x}=0$ and $q_{z}=0$, we only show a quarter of the
results here. Compared to the results in the previous work Blinne:2015zpa  for
the single field, our $f_{-}(q_{x},q_{z})$ is similar to that $f_{+}$ of Fig.5
in Blinne:2015zpa  since the field is opposite to ours so that our $f_{-}$
corresponds to $f_{+}$ of them. However, our $f(q_{x},q_{z})$ and
$f_{+}(q_{x},q_{z})$ are different from them. We think the differences may be
attributed to that the field of Blinne:2015zpa  is for the multicycle case,
and ours is for the subcycle case. In Blinne:2015zpa  they emphasized that the
pattern is a typical interference effect, however, in our case here we have
observed a delicate relation of interference to the weak spiral. For the
combination fields, see the right column of Fig. 2, similar to the results in
the ($q_{x}$, $q_{y}$) plane of Fig. 1, all momentum distributions have
exhibited the interference or/and weak spiral including the magnitude
amplifier effect (Fig. 2f) compared to that in the single field (Fig. 2c).

As $N$ increases from $N=0.5$ to $N=1$ and $N=2$, we take the results of
$f_{-}(q_{x},q_{y})$ in the $\mathbf{E}_{1}(t)$ as an example, see Fig. 3, and
find that the distribution range of the spiral becomes larger and the spiral
arms become thinner as well as longer, while the symmetry of spiral becomes
better in the both of $q_{x}$ and $q_{y}$ directions (To see more clearly, we
amplify the result of $N=2$ and show a quarter of the ones.) These phenomena
are associated with the fact that the pair production tends more and more to
be the multiphoton absorption process with increasing $N$.

Spin asymmetry degree.–The spin asymmetry degree on the momentum distributions
$f(q_{x},q_{y})$ in the combined fields with different time delays are shown
in Fig. 4, in which there are the following three cases in general. (i) For
$T_{d}=0$, the $f_{+}(q_{x},q_{y})$ and $f_{-}(q_{x},q_{y})$ are uniformly
inter-nested in the whole plane, and the fraction of two distributions is
almost $1:1$. (ii) For the small time delay $T_{d}=T$, there is a preliminary
separation between the $f_{+}(q_{x},q_{y})$ and $f_{-}(q_{x},q_{y})$, but the
$f_{-}(q_{x},q_{y})$ is dominant. (iii) For the large time delay $T_{d}=8T$,
the $f_{+}(q_{x},q_{y})$ and $f_{-}(q_{x},q_{y})$ are completely separated, in
which the upper half plane and the lower half plane are dominated by the
$f_{-}(q_{x},q_{y})$ and $f_{+}(q_{x},q_{y})$, respectively.

Figure 4: (color online). Spin asymmetry degree on the momentum spectra in the
($q_{x}$, $q_{y}$) plane (where $q_{z}=0$) under the combination fields for
different time delays with $N=0.5$. (a), (b) and (c) correspond to $T_{d}=0$,
$T_{d}=T$ and $T_{d}=8T$, respectively. The gray line denotes $\kappa_{f}=0$.

.

For convenience, we define the spin asymmetry degree as

$\displaystyle\kappa_{f}=\frac{f_{+}(q_{x},q_{y})-f_{-}(q_{x},q_{y})}{f_{+}(q_{x},q_{y})+f_{-}(q_{x},q_{y})}.$
(4)

First, for $T_{d}=0$, one can see that the distributions of $\kappa_{f}$ has a
good symmetry in the $q_{x}$ direction, see Fig. 4(a), which is related to the
symmetry of corresponding momentum distributions. Moreover, it is found that
the $f_{+}(q_{x},q_{y})$ and $f_{-}(q_{x},q_{y})$ are uniformly inter-nested
in Fig. 4(a). This is because the the fraction of $f_{+}(q_{x},q_{y})$ and
$f_{-}(q_{x},q_{y})$ is almost $1:1$. Second, as the time delay increases to
$T_{d}=T$, the symmetry of $\kappa_{f}$ is destroyed, at the same time, a wide
and nonuniform spiral structure can be observed, see Fig. 4(b). In addition,
there is preliminary separation between the $f_{+}(q_{x},q_{y})$ and
$f_{-}(q_{x},q_{y})$, where the lower half plane is mainly dominated by the
$f_{+}(q_{x},q_{y})$, while the results in the upper half plane are opposite.
But in the whole plan, the $f_{-}(q_{x},q_{y})$ is dominant.

Third, as $T_{d}$ further increases to $T_{d}=8T$, the distribution of
$\kappa_{f}$ is approximately symmetric in the $q_{x}$ direction, and we can
observe the thin and uniform spiral fringes, see Fig. 4(c). Moreover, it is
found that there exists a clear dividing line in this distribution, where the
lower half plane and the upper half plane are primarily dominated by the
$f_{+}(q_{x},q_{y})$ and $f_{-}(q_{x},q_{y})$, respectively, which are
associated with the ones of Figs. 1(e) and (f). This clear dividing phenomenon
may be attributed to the fact that when the $\mathbf{E}_{1}(t)$ and
$\mathbf{E}_{2}(t)$ are combined, as we described above, the result with a
large time delay is almost contributed by the superposition of the results in
two individual fields.

Figure 5: (color online). Spin asymmetry degree on the number density of
created particles dependence on time delays by the combinational fields. (a)
for $N=0.5$. (b) for $N=2$, $E_{02}=0.05\sqrt{2}E_{\mathrm{cr}}$ and
$\omega_{2}=1.2$. Other parameters are the same as in Fig. 1.

.

The spin asymmetry degree on the number density of created particles in the
combined fields with different time delays are shown in Fig. 5, where we
define the spin asymmetry degree as $\kappa_{n}=\kappa^{1}_{n}-\kappa^{2}_{n}$
with $\kappa^{1}_{n}=n_{+}/(n_{+}+n_{-})$ and
$\kappa^{2}_{n}=n_{-}/(n_{+}+n_{-})$. From Fig. 5(a), we find that for
$T_{d}=0$, $\kappa^{1}_{n}$ is smaller than $\kappa^{2}_{n}$, i.e., the
$n_{-}$ is dominant. As $T_{d}$ increases, a crossing point between
$\kappa^{1}_{n}$ and $\kappa^{2}_{n}$ occurs, which is located in the range of
$T_{d}=1\sim 2T$. It indicates that there exists a spin-flip effect in the
range. When $T_{d}\geq 2T$, $\kappa^{1}_{n}$ is larger than $\kappa^{2}_{n}$,
i.e., $n_{+}$ is dominant. Based on the above results, we can expect to
control the momentum spectrum via changing the time delay, for example, if we
want to obtain a momentum spectrum dominated by $f_{-}(q_{x},q_{y})$, we can
set the time delay $T_{d}$ to be smaller. While if we want to get the spectrum
dominated by $f_{+}(q_{x},q_{y})$, the $T_{d}$ needs to be extended. Moreover,
one can see that with the increase of $T_{d}$, the $\kappa_{n}$ first
increases and then remains unchanged as $12\%$. In order to obtain a higher
spin asymmetry degree, we have performed some other simulations, for example
the result shown in Fig. 5(b). It is found that the $\kappa_{n}$ monotonically
decreases with increasing $T_{d}$, specifically, it decreases from $-50\%$ at
$T_{d}=0$ and eventually reaches to $-98\%$ at $T_{d}=8T$.

Discussion.– Different from the previous studies Li:2017qwd ; Hu:2023pmz , we
found that the spin can induce the momentum spirals in a single field. The
bichromatic combinational fields lead to an inhomogeneous spiral structures of
the momentum, especially, the spiral pattern does exist in both of $f_{+}$ and
$f_{-}$, and it is about two orders of magnitude larger than the single field.
It seems that one field provides to create spiral ”seed” for a certain spin
type, while the other field plays a role as the ”amplifier” to this spiral.
For the opposite spin case, the role of ”seed” and ”amplifier” played by the
fields is exchangeable.

It is known that the purely time-dependent electric field actually consists of
two counter-propagating waves Kohlfurst:2022edl  and they are focused on the
narrow regime of original electron Compton scales, i.e., a dipole
approximation is met. For our single field, it is equivalent to a combination
of a left-handed wave and a right-handed wave, and the two counter-propagating
waves are attributing for most of the interaction time at the focusing spatial
scale. So this is very similar to the case of combined fields that are
consisted of two counter-rotating pulses. This is maybe a reason why the
spiral can exist even if the single field is present.

These results suggest that the spin effect of created particles plays an
important role in momentum spiral on pair production. It can not only induce
the momentum spiral in a single field, but also contribute to the formation of
the spiral in the bichromatic combinational fields as well as the appearance
of spin flip effect. Moreover, the momentum spiral due to the spin is
sensitive to the field parameters. It is found that either the larger the time
delay or the smaller the number of cycles, the more pronounced the momentum
spiral, which provides a flexible way to control the spiral features. While we
have only considered the simple case of CP fields, we speculate that the more
delicate spin effects exist in different cases of fields and further features
associated to the momentum patterns on pair production are expected.

Acknowledgements. We are grateful to ZL Li and C Kohlfürst for helpful
discussions. This work was supported by the National Natural Science
Foundation of China (NSFC) under Grant No. 12375240, No. 11935008 and No.
12104428. The computation was carried out at the HSCC of the Beijing Normal
University.

## References

  * (1) F. Sauter, Uber das Verhalten eines Elektrons im homogenen elektrischen Feld nach der relativistischen Theorie Diracs, Z. Phys. 69, 742 (1931).
  * (2) W. Heisenberg and H. Euler, Consequences of Dirac’s theory of the positron, Z. Phys. 98, 714 (1936).
  * (3) A. Di Piazza, C. Muller, K. Z. Hatsagortsyan, and C. H. Keitel, Extremely high-intensity laser interactions with fundamental quantum systems, Rev. Mod. Phys. 84, 1177 (2012).
  * (4) A. Fedotov, A. Ilderton, F. Karbstein, B. King, D. Seipt, H. Taya, and G. Torgrimsson, Advances in QED with intense background fields, Phys. Rept. 1010, 1 (2023).
  * (5) B. S. Xie, Z. L. Li, and S. Tang, Electron-positron pair production in ultrastrong laser fields, Matter Radiation. Extremes 2, 225 (2017).
  * (6) G. V. Dunne, New strong-field QED effects at ELI: nonperturbative vacuum pair production, Eur. Phys. J. D 55, 327 (2009).
  * (7) J. S. Schwinger, On gauge invariance and vacuum polarization, Phys. Rev. 82, 664 (1951).
  * (8) D. L. Burke, R. C. Field, G. Horton-Smith, J. E. Spencer, D. Walz, S. C. Berridge, W. M. Bugg, K. Shmakov, A. W. Weidemann, C. Bula, K. T. McDonald, and E. J. Prebys, Positron production in multi-photon light by light scattering, Phys. Rev. Lett. 79, 1626 (1997); C. Bamber, et al., Studies of nonlinear QED in collisions of 46.6-GeV electrons with intense laser pulses, Phys. Rev. D 60, 092004 (1999).
  * (9) C. Kohlfürst, H. Gies, and R. Alkofer, Effective mass signatures in multiphoton pair production, Phys. Rev. Lett. 112, 050402 (2014).
  * (10) F. Hebenstreit, R. Alkofer, and H. Gies, Particle self-bunching in the Schwinger effect in spacetime-dependent electric fields, Phys. Rev. Lett. 107, 180403 (2011).
  * (11) C. Kohlfürst and R. Alkofer, Ponderomotive effects in multiphoton pair production, Phys. Rev. D 97, 036026 (2018).
  * (12) Z. L. Li, D. Lu, B. S. Xie, L. B. Fu, J. Liu and B. F. Shen Enhanced pair production in strong fields by multiple-slit interference effectwith dynamically assisted Schwinger mechanism, Phys. Rev. D 89, 093011 (2014).
  * (13) Z. L. Li, D. Lu, B. S. Xie, B. F. Shen, L. B. Fu, and J. Liu, Nonperturbative signatures in pair production for general elliptic polarization fields, Europhys. Lett. 110, 51001 (2015).
  * (14) Z. L. Li, Y. J. Li, and B. S. Xie, Momentum vortices on pairs production by two counter-rotating fields, Phys. Rev. D 96, 076010 (2017).
  * (15) L. N. Hu, O. Amat, L. Wang, A. Sawut, H. H. Fan, and B. S. Xie, Momentum spirals in multiphoton pair production revisited, Phys. Rev. D 107, 116010 (2023).
  * (16) Z. L. Li, B. S. Xie, and Y. J. Li, Vortices in multiphoton pair production by two-color rotating laser fields, J. Phys. B 52, 025601 (2019).
  * (17) J. M. Ngoko Djiokap, S. X. Hu, L. B. Madsen, N. L. Manakov, A. V. Meremianin, and Anthony F. Starace1, Electron vortices in photoionization by circularly polarized attosecond pulses, Phys. Rev. Lett. 115, 113004 (2015).
  * (18) D. Pengel, S. Kerbstadt, D. Johannmeyer, L. Englert, T. Bayer, and M. Wollenhaupt, Electron vortices in femtosecond multiphoton ionization, Phys. Rev. Lett. 118, 053003 (2017).
  * (19) M. M. Majczak, F. Cajiao Vélez, J. Z. Kamiński, and K. Krajewska, Carrier-envelope-phase and helicity control of electron vortices and spirals in photodetachment, Opt. Express 30, 43330 (2022).
  * (20) J. H. Macek, J. B. Sternberg, S. Y. Ovchinnikov, T. G Lee, and D. R. Schultz, Origin, evolution, and imaging of vortices in atomic processes, Phys. Rev. Lett. 102, 143201 (2009).
  * (21) K. J. Yuan, S. Chelkowski, and A. D. Bandrauk, Photoelectron momentum distributions of molecules in bichromatic circularly polarized attosecond uv laser fields, Phys. Rev. A 93, 053425 (2016).
  * (22) L. Geng, F. C. Vélez, J. Z. Kamiński, L. Y. Peng, and K. Krajewska, Vortex structures in photodetachment by few-cycle circularly polarized pulses, Phys. Rev. A 102, 043117 (2020).
  * (23) M. Harris, C. A. Hill, and J. M. Vaughan, Optical helices and spiral interference fringes, Opt. Commun. 106, 161 (1994).
  * (24) G. Blatter, M. V. Feigelman, V. B. Geshkenbein, A. I. Larkin, and V. M. Vinokur, Vortices in high-temperature superconductors, Rev. Mod. Phys. 66, 1125 (1994).
  * (25) L. Uby, M. B. Isichenko, and V. V. Yankov, Vortex filament dynamics in plasmas and superconductors, Phys. Rev. E 52, 932 (1995).
  * (26) P. K. Shukla and B. Eliasson, Colloquium: Fundamentals of dust-plasma interactions, Rev. Mod. Phys. 81, 25 (2009).
  * (27) K. W. Madison, F. Chevy, W. Wohlleben, and J. Dalibard, Vortex formation in a stirred Bose-Einstein Condensate, Phys. Rev. Lett. 84, 806 (2000).
  * (28) E. Strobel and S. S. Xue, Semiclassical pair production rate for rotating electric fields, Phys. Rev. D 91, 045016 (2015).
  * (29) A. Blinne and E. Strobel, Comparison of semiclassical and Wigner function methods in pair production in rotating fields, Phys. Rev. D 93, 025014 (2016).
  * (30) Z. L. Li, B. S. Xie, and Y. J. Li, Boson pair production in arbitrarily polarized electric fields, Phys. Rev. D 100, 076018 (2019).
  * (31) C. Kohlfürst, Spin states in multiphoton pair production for circularly polarized light, Phys. Rev. D 99, 096017 (2019).
  * (32) C. Kohlfürst, Electron-positron pair production in inhomogeneous electromagnetic fields, Ph.D. thesis, arXiv:1512.06082 (2015).
  * (33) F. Hebenstreit, Schwinger effect in inhomogeneous electric fields, Ph.D. thesis, arXiv:1106.5965 (2011).
  * (34) I. Bialynicki-Birula, P. Gornicki and J. Rafelski, Phase space structure of the Dirac vacuum, Phys. Rev. D 44, 1825 (1991).
  * (35) F. Hebenstreit, R. Alkofer, and H. Gies, Schwinger pair production in space and time-dependent electric fields: Relating the Wigner formalism to quantum kinetic theory, Phys. Rev. D 82, 105026 (2010).
  * (36) M. Ababekri, B. S. Xie, and J. Zhang, Effects of finite spatial extent on Schwinger pair production, Phys. Rev. D 100, 016003 (2019).
  * (37) C. Kohlfürst, Effect of time-dependent inhomogeneous magnetic fields on the particle momentum spectrum in electron-positron pair production, Phys. Rev. D 101, 096003 (2020).
  * (38) A. Blinne and H. Gies, Pair production in rotating electric fields, Phys. Rev. D 89, 085001 (2014).
  * (39) L. J. Li, M. Mohamedsedik, and B. S. Xie, Enhanced dynamically assisted pair production in spatial inhomogeneous electric fields with the frequency chirping, Phys. Rev. D 104, 036015 (2021).
  * (40) O. Olugh, Z. L. Li, B. S. Xie, and R. Alkofer, Pair production in differently polarized electric fields with frequency chirps, Phys. Rev. D 99, 036003 (2019).
  * (41) Z. L. Li, D. Lu, and B. S. Xie, Effects of electric field polarizations on pair production, Phys. Rev. D 92, 085001 (2015).
  * (42) C. Kohlfürst, Phase-space analysis of the Schwinger effect in inhomogeneous electromagnetic fields, Eur. Phys. J. Plus 133, 191 (2018).
  * (43) O. Olugh, Z. L. Li, and B. S. Xie, Dynamically assisted pair production for various polarizations, Phys. Lett. B 802, 135259 (2020).
  * (44) C. K. Dumlu and G. V. Dunne, The stokes phenomenon and Schwinger vacuum pair production in time-dependent laser pulses, Phys. Rev. Lett. 104, 250402 (2010).
  * (45) J. Oertel and R. Schützhold, WKB approach to pair creation in spacetime-dependent fields: The case of a spacetime-dependent mass, Phys. Rev. D 99, 125014 (2019).
  * (46) C. Kohlfürst, Pair Production in Circularly Polarized Waves, arXiv:2212.03180.


11institutetext:  Laboratoire d’astrophysique de Bordeaux, Univ. Bordeaux,
CNRS, B18N, allée Geoffroy Saint-Hilaire, 33615 Pessac, France

#  Spectral binning of precomputed correlated-$k$ coefficients ††thanks: I
wish to dedicated this paper to the memory of Adam P. Showman and Franck
Hersant. You both were, each in your own way, an inspiration to me. Godspeed
on your journey among the stars.

Jérémy Leconte 11

With the major increase in the volume of the spectroscopic line lists needed
to perform accurate radiative transfer calculations, disseminating accurate
radiative data has become almost as much a challenge as computing it.
Considering that many planetary science applications are only looking for
heating rates or mid-to-low resolution spectra, any approach enabling such
computations in an accurate and flexible way at a fraction of the computing
and storage costs is highly valuable. For many of these reasons, the
correlated-$k$ approach has become very popular. Its major weakness has been
the lack of ways to adapt the spectral grid/resolution of precomputed
$k$-coefficients, making it difficult to distribute a generic database suited
for many different applications. Currently, most users still need to have
access to a line-by-line transfer code with the relevant line lists or high-
resolution cross sections to compute $k$-coefficient tables at the desired
resolution. In this work, we demonstrate that precomputed $k$-coefficients can
be binned to a lower spectral resolution without any additional assumptions,
and show how this can be done in practice. We then show that this binning
procedure does not introduce any significant loss in accuracy. Along the way,
we quantify how such an approach compares very favorably with the sampled
cross section approach. This opens up a new avenue to deliver accurate
radiative transfer data by providing mid-resolution $k$-coefficient tables to
users who can later tailor those tables to their needs on the fly. To help
with this final step, we briefly present Exo_k, an open-access, open-source
Python library designed to handle, tailor, and use many different formats of
$k$-coefficient and cross-section tables in an easy and computationally
efficient way.

###### Key Words.:

Radiative transfer, correlated-k method, Atmosphere

††offprints<EMAIL_ADDRESS>
##  1 Introduction

Despite the considerable increase in computing power, line-by-line radiative
transfer is still considered a computationally intensive task in many cases of
interest. One of the reasons for this state of affairs is that recent progress
in spectroscopy has multiplied the size of molecular line lists by more than
three orders of magnitude (Tennyson & Yurchenko, 2012). We are at a point at
which line list themselves, and a fortiori, high-resolution cross sections,
represent a data volume that is not trivially handled and exchanged.

To face this while maintaining flexibility, some have opted for sampled or
binned down cross-section tables111Here, sampled cross sections specifically
refer to cross sections that are computed monochromatically but at a
resolution that does not necessarily resolve the lines or their shape. Binned-
down cross sections are computed on a finely sampled spectral grid and then
integrated over a coarser grid, thereby preserving area. (Line et al., 2015;
Waldmann et al., 2015). But Garland & Irwin (2019) showed that this approach
is inaccurate if the resolution used is not sufficiently high, and too
computationally expensive if the resolution is high enough (we briefly revisit
this question in Sect. 2 whose results are summarized in Figs. 1 and 2). It
seems that for typical resolutions below $\sim$1000, methods designed to group
absorptions, such as the so-called correlated-$k$ method (Liou, 1980; Lacis &
Oinas, 1991), are more efficient while remaining accurate (Irwin et al., 2008;
Amundsen et al., 2017; Garland & Irwin, 2019; Zhang et al., 2020). It is
especially true for multidimensional models that need a very fast radiative
transfer and use only 10-30 bins for the whole spectrum (Showman et al., 2009;
Wordsworth et al., 2011; Leconte et al., 2013; Amundsen et al., 2017).

For these reasons, it is interesting to directly distribute reference
$k$-coefficients for each molecule instead of high-resolution cross sections.
This is currently done for the ExoMol project (Chubb et al., 2020). This
approach decreases the data volume necessary to handle while keeping an
optimal accuracy. A current drawback with precomputed $k$-coefficients,
however, is that the person computing the tables needs to choose a spectral
resolution (or wavenumber grid) without knowing exactly what the table will be
used for, while $k$-coefficients are used at their highest potential when
tailored to the needed resolution.

To circumvent this problem we present a way to accurately bin down precomputed
$k$-coefficients to an arbitrary spectral grid without any significant
additional loss in accuracy. This method is faster than computing the
$k$-coefficients directly from line-by-line data. It therefore opens up a new
way of disseminating accurate, reference molecular opacities: Spectroscopists
only need to compute $k$-coefficients once with a resolution that is higher
than the highest resolution they envision their users need (e.g., the
resolution of the instrument you want to compare your model to). Users can
then easily bin down the data to their spectral grid of choice, possibly even
on the fly.

To help with the latter, we developed `Exo_k`222See the online documentation
for tutorials, example scripts and notebooks, and an extensive and up-to-date
description of the (evolving) features of the library:
http://perso.astrophy.u-bordeaux.fr/~jleconte/exo_k-doc/index.html , an open-
source Python 3 library that can directly load precomputed $k$-coefficient
tables, change their spectral resolution as well as pressure-temperature
grids, and save them back on disk. The library can handle many different
formats used in various codes: for example, TauREx (Al-Refaie et al. 2019;
both pickle and HDF5), petitRADTRANS (Mollière et al., 2019), Nemesis (Irwin
et al., 2008), ARCIS (Min et al., 2020), Exo transmit (Kempton et al., 2017),
LMD Generic global climate model (Wordsworth et al., 2011; Leconte et al.,
2013). Several of these formats are provided by the latest ExoMol release
(Chubb et al., 2020). Other formats are being implemented and new ones can be
added on request. It is also possible to compute $k$-coefficient tables from
high-resolution spectra. One can also combine the opacities of several
molecules using the random overlap method (Lacis & Oinas, 1991) to create a
table for a given atmospheric mix, which is very useful for climate models,
for example. Finally, the library includes methods to directly interpolate and
combine the opacities, and even compute transmission and emission spectra for
1D planetary atmospheres. Leveraging on the `Numba` library, these operations
are performed in a very computationally efficient way. So `Exo_k` can be
imported directly into any radiative transfer code to handle molecular
opacities easily and efficiently.

In this article, we first briefly compare the relative accuracy and efficiency
of the sampled cross section and correlated-$k$ approaches in Sect. 2. After
introducing the necessary concepts and notations, Sect. 3 demonstrates
rigorously why and how to bin down precomputed $k$-coefficients to an
arbitrary resolution. Then, in Sect. 4, we validate the numerical algorithm
presented in this work and implemented in `Exo_k`.

##  2 Sampled cross sections versus correlated-$k$ tables

To accurately compare the efficiency of various radiative transfer methods
when a given resolution and precision are needed, we briefly estimate the
accuracy of atmospheric modeling when using the following types of input
radiative data: sampled and binned cross section and $k$-coefficient tables.

We compare this for two different atmospheric configurations and observation
geometries:

  * $\bullet$

Transmission spectrum of a hot Jupiter with a 1000 K isothermal H2/H2O
atmosphere in front of a Sun.

  * $\bullet$

Emission spectrum of a pure CO2 Mars-like atmosphere.

The parameters for these two cases are given in Table 1.

Table 1: Standard parameters used in the our two fiducial cases. A hyphen stands for parameters that are not needed for a given case. Case |  | Hot Jupiter | Mars
---|---|---|---
Planet radius | $R_{\mathrm{Jup}}$ | 1. | -
Surface gravity | m s-2 | 10. | 3.72
Surface temperature | K | 1000. | 200.
Stratospheric temperature | K | 1000. | 100.
Adiabatic index |  | - | 0.22
Surface pressure | Pa | 106 | 640.
Volume mixing ratio of H2O |  | 10-3 | -
Stellar radius | $R_{\odot}$ | 1. | -
Figure 1: Root mean square error on the transit spectrum of our fiducial hot
Jupiter over the 1-2 $\mu$m window as a function of the final observed
resolution. Dashed curves represent the error obtained with sampled cross-
section tables at various sampling resolution (from top to bottom:
$R_{\mathrm{sp}}=\\{1000,\,3000,\,10\,000,\,30\,000\\}$). The spectra are
shown in Fig. 2. The black curve with dots represents the error using
$k$-coefficients computed from high-resolution spectra directly at the final
resolution. The gray curves show results when using binned down
$k$-coefficients using natural (solid with stars) and non-integer (dotted)
binning (See Sect. 3.4 for details).

Each time, we start by computing a reference spectrum using high-resolution
monochromatic cross sections (with a typical step of 0.01cm-1). For our CO2
atmosphere, we use spectra from Turbet & Tran (2017). For our hot-Jupiter
case, we use water spectra produced by the Exomol project with H2 broadening
(K. Chubb, private communication) based on the line list produced by Polyansky
et al. (2018). Collision-induced-absorptions (CIA) are taken from the HITRAN
database (Richard et al., 2011).

Then, the high-resolution cross sections are sampled at various resolutions
($R_{\mathrm{sp}}=\\{1\,000,\,3\,000,\,10\,000,\,30\,000\\}$) and a new set of
spectra are computed at each of these resolutions.

The spectra obtained from the sampled cross sections are then binned
(preserving area) at selected final resolutions
($R_{\mathrm{fin}}=\\{10,\,20,\,50,\,100,\,200,\,500,\,1\,000\\}$). The
transmission spectra computed in the hot-Jupiter case are shown in the left
column of Fig. 2. Finally, these spectra are compared to the properly binned
reference spectrum (right column of Fig. 2). The root mean squared (RMS) error
is shown for each ($R_{\mathrm{sp}}$, $R_{\mathrm{fin}}$) in Fig. 1. The
emission spectra computed in the Mars-like case are shown in Fig. 3.

For the correlated-$k$ approach, $k$-coefficient tables are computed from the
high-resolution spectra directly at the final desired resolution (computed
with 20 Gauss-Legendre quadrature points in $g$-space) and compared to the
reference spectrum binned at the same resolution. This gives the black dashed
curves in Fig. 2. The RMS error is shown by the black curve in Fig. 1.

For all the computations described above we used a standard algorithm to
describe the transmission and emission spectra of a 1D atmosphere that have
been implemented in the `Exo_k` library${}^{\ref{exok_foot}}$.

Figure 2: Effect of sampling resolution on the transit spectrum of a hot
Jupiter. The various rows show the effect for different final observing
resolutions ($R_{\mathrm{fin}}$, specified in the top left corner of the
panels). The left column shows the transit depth in the 1-2 $\mu m$ window
(with an offset of 11500 ppm). The various colors correspond to the four
sampling resolutions (i.e., the resolution used in radiative transfer
calculation before binning; $R_{\mathrm{sp}}$) shown in the legend. The curve
with dots shows the reference case computed from our very high-resolution
cross sections ($\Delta\sigma\approx 0.01$cm-1). The high-resolution reference
spectrum before binning is shown as the shading in the top left panel. The
right column shows the difference (in ppm) between the spectra of the left
panel and the reference case. The difference between the calculation with
$k$-coefficients and the reference case is shown with the black dashed curve;
the two spectra would be indistinguishable in the left panels. The RMS
standard deviation over the spectral region is shown in Fig. 1.

We show that the correlated-$k$ approach reaches an accuracy that is one to
two orders of magnitude better than sampled cross sections, especially for
resolutions bigger than 100. The relative drop in accuracy of the
correlated-$k$ method at $R_{\mathrm{fin}}=10$ is due to the low sampling of
the CIA absorption, which is treated as constant in each spectral bin. This
could be alleviated by applying the random overlap method to the CIA as well.
To compare numerical efficiency, we focus on the $R_{\mathrm{fin}}$=10 case
for which the $R_{\mathrm{sp}}=30\,000$ roughly manages to reach an accuracy
equivalent to the correlated-$k$ method. Then, considering that our
$k$-coefficients were computed using 20 quadrature points, the computation is
still 150 times faster with the correlated-$k$ method compared to the sampled
cross sections.

Figure 3: Effect of sampling resolution on the emission spectrum of our Mars-
like case around the 15$\mu$m CO2 band. The rows show the effect for different
final observing resolutions ($R_{\mathrm{fin}}$). The left column shows the
flux. The various colors correspond to the four sampling resolutions (i.e.,
the resolution used in radiative transfer calculation before binning;
$R_{\mathrm{sp}}$) shown in the legend. The reference case is computed from
our very high-resolution cross sections ($\Delta\sigma\approx 0.01$cm-1). The
high-resolution reference spectrum before binning is shown as the shading in
the top left panel. The right column shows the relative error between the
spectra of the left panel and the reference case. The error between the
calculation with $k$-coefficients and the reference case is shown with the
black dashed curve; the two spectra would be indistinguishable in the left
panels. Figure 4: Effect of sampling resolution on the transit spectrum of a
hot Jupiter using cross sections that have been binned by preserving the total
area. The panels can be compared to the results for sampled cross sections in
Fig. 2 (see the caption for details). The binning systematically overestimates
the opacities.

We conclude this comparison by mentioning binned-down cross sections – in
which the finely sampled cross sections are integrated over a coarser grid,
thereby preserving area. Although it seems the most natural way to bin cross
sections, this method performs even worse than the sampled cross sections (by
about one order of magnitude; Fig. 4). This rather counterintuitive result
stems from the nonlinear nature of radiative transfer: the opacity dilution
due to the binning procedure is not sufficient to increase the transparency of
the gas near optically thick lines, while decreasing it in the line wings.
Overall, this systematically over estimates the opacity of the atmosphere and
creates a systematic bias – on the order of 100 ppm in our hot-Jupiter test
case.

##  3 Method to bin down $k$-coefficients

###  3.1 Correlated-$k$ formalism

Before demonstrating how to bin down $k$-coefficients, we briefly introduce
the basic concepts we need. In the correlated-$k$ method (see Liou 1980 for
details), the wavenumber space, where $\nu$ is the wavenumber, is divided into
$N_{\mathrm{b}}$ spectral bands of width
$\Delta\nu^{\mathrm{b}}\equiv\nu^{\mathrm{b}}_{+}-\nu^{\mathrm{b}}_{-}$. In
each band, the transmission through a slab of matter can be written as

$\displaystyle\mathcal{T}^{\mathrm{b}}(u)=\frac{1}{\Delta\nu^{\mathrm{b}}}\int_{\nu^{\mathrm{b}}_{-}}^{\nu^{\mathrm{b}}_{+}}\textsl{e}^{-k^{\nu}u}\mathrm{d}\nu=\int_{0}^{1}\textsl{e}^{-\hat{k}^{\mathrm{b}}(g)u}\mathrm{d}g\approx\sum_{\mathrm{i}=1}^{N_{g}}w^{\mathrm{b}}_{\mathrm{i}}\textsl{e}^{-\hat{k}^{\mathrm{b}}_{\mathrm{i}}u},$
(1)

where $N_{g}$ is the number of points (or abscissas) used to discretize the
$g$-space and $w^{\mathrm{b}}$ are the associated set of weights for each
band. The symbol $\hat{k}^{\mathrm{b}}(g)$ denotes the inverse of the
cumulative density function of the opacity within band $\mathrm{b}$ – the so-
called $k$-distribution – and should not be confused with $k^{\nu}$, the
absorption coefficient, which is a function of the wavenumber. The so-called
correlated-$k$ coefficients (or simply $k$-coefficients) are the
$\hat{k}^{\mathrm{b}}_{\mathrm{i}}$, which are the discretized version of the
$k$-distribution ($\hat{k}^{\mathrm{b}}(g)$).

An intuitive way to understand that is to say that now, in band $\mathrm{b}$,
the opacity can be described by $N_{g}$ representative values
($\hat{k}^{\mathrm{b}}_{\mathrm{i}}$) for each of the $N_{g}$ points in
$g$-space ($g_{\mathrm{i}}$), each of these values being affected a weight
($w^{\mathrm{b}}_{\mathrm{i}}$). Then any radiative transfer calculation can
be computed separately in the $N_{g}$ bins to be later summed up with the
proper weights.

An important, although often forgotten assumption of the correlated-$k$
formalism is that both the spectral incoming flux impinging on our medium
($F_{\mathrm{in}}^{\nu}$) and the source function of the latter should be
nearly constant within each spectral band (Lacis & Oinas, 1991) so that it can
be written
$F_{\mathrm{in}}^{\nu}=F_{\mathrm{in}}^{\mathrm{b}}/\Delta\nu^{\mathrm{b}}$,
where
$F_{\mathrm{in}}^{\mathrm{b}}\equiv\int_{\nu^{\mathrm{b}}_{-}}^{\nu^{\mathrm{b}}_{+}}F_{\mathrm{in}}^{\nu}\mathrm{d}\nu$
is the total flux within the band. Simply put in the context of a purely
absorbing medium, which can be generalized, this ensures that the outgoing
flux after a path $u$ is given by

$\displaystyle
F_{\mathrm{out}}^{\mathrm{b}}(u)\equiv\int_{\nu^{\mathrm{b}}_{-}}^{\nu^{\mathrm{b}}_{+}}F_{\mathrm{in}}^{\nu}\textsl{e}^{-k^{\nu}(\nu)u}\mathrm{d}\nu$
$\displaystyle\approx\frac{F_{\mathrm{in}}^{\mathrm{b}}}{\Delta\nu^{\mathrm{b}}}\int_{\nu^{\mathrm{b}}_{-}}^{\nu^{\mathrm{b}}_{+}}\textsl{e}^{-k^{\nu}(\nu)u}\mathrm{d}\nu$
$\displaystyle=F_{\mathrm{in}}^{\mathrm{b}}\cdot\mathcal{T}^{\mathrm{b}}(u).$
(2)

This assumption is important in our context because it is the only one that we
need to make to be able to combine the $k$-coefficients of various bands.
Finally, the bolometric flux can be obtained with

$\displaystyle
F=\sum_{\mathrm{b}=1}^{N_{\mathrm{b}}}F_{\mathrm{out}}^{\mathrm{b}}.$ (3)

###  3.2 Combining $k$-coefficients of various bands

To reduce the resolution of precomputed $k$-coefficients, we essentially need
to answer the following question:

  * $\bullet$

Can we determine the $k$-coefficients
$\hat{k}^{\bar{\mathrm{b}}}_{\mathrm{i}}$ of a ”super-band” $\bar{\mathrm{b}}$
that is the union of two smaller, nonoverlapping ”sub-bands”333We restrict
ourselves to the case of the combination of two bands. The generalization to
an arbitrary number of bands is trivial., $\mathrm{b}_{1}$ and
$\mathrm{b}_{2}$, for which we know the $k$-coefficients,
$\hat{k}^{\mathrm{b}_{1}}_{\mathrm{i}}$ and
$\hat{k}^{\mathrm{b}_{2}}_{\mathrm{i}}$?

In this form, the answer is not trivial. However, when computing
$k$-coefficients, the crucial mathematical object is not the $k$-distribution
($\hat{k}^{\mathrm{b}}(g)$) but its inverse, the cumulative density function
(CDF) of the opacity within the band (hereafter called the $g$-distribution,
$\hat{g}^{\mathrm{b}}(k)$) defined as

$\displaystyle\hat{g}^{\mathrm{b}}(k)\equiv\frac{1}{\Delta\nu^{\mathrm{b}}}\int_{\nu^{\mathrm{b}}_{-}}^{\nu^{\mathrm{b}}_{+}}H\left(k-k^{\nu}(\nu)\right)\mathrm{d}\nu,$
(4)

where $H(x)$ is equal to 0 where $x$ is positive and 1 elsewhere.

Figure 5: Example of the binning process for two bands. Left: $k$-coefficients
on $N_{g}=$ 20 Gauss-Legendre quadrature points for two bands
($\hat{k}^{\mathrm{b}_{1}}_{\mathrm{i}}$ and
$\hat{k}^{\mathrm{b}_{2}}_{\mathrm{i}}$). Middle: $g$-distribution for the two
bands resampled on $\tilde{k}_{\mathrm{j}}$, a grid of $N_{k}=100$ points
between $k_{\mathrm{min}}$ and $k_{\mathrm{max}}$
($\tilde{g}_{\mathrm{j}}^{\mathrm{b}_{1}}$ and
$\tilde{g}_{\mathrm{j}}^{\mathrm{b}_{2}}$). The black line indicates the
weighted sum of these two $g$-distributions, i.e., the $g$-distribution of the
super-band ($\tilde{g}_{\mathrm{j}}^{\bar{\mathrm{b}}}$). Right: The colored
lines are shown as in the left panel. The black line indicates the finely
sampled $k$-distribution
($\tilde{k}_{\mathrm{j}}(\tilde{g}_{\mathrm{j}}^{\bar{\mathrm{b}}})$). The
gray line represents the $k$-distribution of the super-band re-sampled on the
original g-grid with $N_{g}=$ 20, i.e., the final $k$-coefficients
($\hat{k}^{\bar{\mathrm{b}}}_{\mathrm{i}}$).

The question can thus be reframed to ask whether we can compute the
$g$-distribution of the super-band, $\hat{g}^{\bar{\mathrm{b}}}$, if we know
the $g$-distribution within the two sub-bands ($\hat{g}^{\mathrm{b}_{1}}$ and
$\hat{g}^{\mathrm{b}_{2}}$). This is rather straightforward, because as long
as the two smaller bands do not overlap, we have
$\Delta\nu^{\bar{\mathrm{b}}}=\Delta\nu^{\mathrm{b}_{1}}+\Delta\nu^{\mathrm{b}_{2}}$
and

$\displaystyle\hat{g}^{\bar{\mathrm{b}}}(k)$
$\displaystyle\equiv\frac{1}{\Delta\nu^{\bar{\mathrm{b}}}}\int_{\nu^{\bar{\mathrm{b}}}_{-}}^{\nu^{\bar{\mathrm{b}}}_{+}}H\left(k-k^{\nu}(\nu)\right)\mathrm{d}\nu$
$\displaystyle=\frac{1}{\Delta\nu^{\bar{\mathrm{b}}}}\left(\int_{\nu^{\mathrm{b}_{1}}_{-}}^{\nu^{\mathrm{b}_{1}}_{+}}H\left(k-k^{\nu}(\nu)\right)\mathrm{d}\nu+\int_{\nu^{\mathrm{b}_{2}}_{-}}^{\nu^{\mathrm{b}_{2}}_{+}}H\left(k-k^{\nu}(\nu)\right)\mathrm{d}\nu\right)$
$\displaystyle=\frac{\Delta\nu^{\mathrm{b}_{1}}}{\Delta\nu^{\bar{\mathrm{b}}}}\hat{g}^{\mathrm{b}_{1}}(k)+\frac{\Delta\nu^{\mathrm{b}_{2}}}{\Delta\nu^{\bar{\mathrm{b}}}}\hat{g}^{\mathrm{b}_{2}}(k).$
(5)

Therefore the $g$-distribution in the super-band is just the sum of the
$g$-distributions in the sub-bands weighted by their spectral extent.

This can be understood more intuitively when discussing in terms of
probability. We want to know if $\hat{g}^{\bar{\mathrm{b}}}(k)$, that is, the
probability that a monochromatic ray of light falls in a spectral region where
the opacity $k^{\nu}$ is lower than $k$. To compute this, we have to sum, for
each band, the probability of the ray falling in the given sub-band; this
probability is $\Delta\nu^{\mathrm{b}}/\Delta\nu^{\bar{\mathrm{b}}}$, as long
as the incoming spectral flux is constant over the super-band, multiplied by
the probability that $k^{\nu}<k$ within this sub-band, which is
$\hat{g}^{\mathrm{b}}(k)$. So we recover Eq. (3.2) as long as we keep the
essential assumption of the correlated-$k$ method that the incoming flux and
sources functions do not significantly vary over the wavelength range that we
wish to describe by a single distribution.

###  3.3 Numerical algorithm

As shown above, combining the $g$-distributions of various bands in a single
$g$-distribution for a large band is straightforward. However, to use this in
a radiative transfer code, we usually need the $k$-coefficients
($\hat{k}^{\mathrm{b}}_{\mathrm{i}}$) that are tabulated for fixed $g$-points
($g_{\mathrm{i}}$) with the associated weights
($w^{\mathrm{b}}_{\mathrm{i}}$). The difficulty is that although
$g_{\mathrm{i}}$ can be seen as the values of $\hat{g}^{\mathrm{b}}$ in band
$\mathrm{b}$ for the points located at $\hat{k}^{\mathrm{b}}_{\mathrm{i}}$ in
opacity space, the distribution $\hat{g}^{\mathrm{b}}$ is sampled at different
locations in opacity space in each sub-band.

To circumvent this problem, we recommend the following approach, which is the
one implemented in `Exo_k`. The whole process is illustrated in Fig. 5 with
only two bands, but can be carried out for any number of bands. For the
moment, however, we restrict ourselves to the case in which the super-band
$\bar{\mathrm{b}}$ is composed of an integer number of sub-bands (hereafter,
natural binning). Once these sub-bands have been identified, we determine

$\displaystyle
k_{\mathrm{min}}=\underset{\mathrm{i},\,\mathrm{b}\in\bar{\mathrm{b}}}{\mathrm{min}}\left\\{\hat{k}^{\mathrm{b}}_{\mathrm{i}}\right\\}\
\ \mathrm{and}\ \
k_{\mathrm{max}}=\underset{\mathrm{i},\,\mathrm{b}\in\bar{\mathrm{b}}}{\mathrm{max}}\left\\{\hat{k}^{\mathrm{b}}_{\mathrm{i}}\right\\}$
(6)

and compute an evenly log-distributed grid of $N_{k}$ points
($\tilde{k}_{\mathrm{j}}$) between these two values. Then, for each sub-band,
we interpolate $g_{\mathrm{i}}(\hat{k}^{\mathrm{b}}_{\mathrm{i}})$ on this
grid, which yields $N_{k}$ values per band,
$\tilde{g}^{\mathrm{b}}_{\mathrm{j}}$, representing all the
$\hat{g}^{\mathrm{b}}$ on the same grid for all bands. It is thus now easy to
compute the global $g$-distribution which is given by

$\displaystyle\tilde{g}_{\mathrm{j}}^{\bar{\mathrm{b}}}=\sum_{\mathrm{b}\in\bar{\mathrm{b}}}\frac{\Delta\nu^{\mathrm{b}}}{\Delta\nu^{\bar{\mathrm{b}}}}\tilde{g}^{\mathrm{b}}_{\mathrm{j}}.$
(7)

This generally yields an oversampled version of the distribution because we
choose $N_{k}>N_{g}$ (the black curve in the middle panel of Fig. 5).
Recomputing the $k$-coefficients is thus now simply a matter of resampling
$\tilde{k}_{\mathrm{j}}(\tilde{g}_{\mathrm{j}}^{\bar{\mathrm{b}}})$ onto our
original $g$-grid ($g_{\mathrm{i}}$) – going from the black to the gray curve
in the right panel of Fig. 5 – as discussed in Lacis & Oinas (1991) and
Amundsen et al. (2017). However, unlike those authors, we think that, to be
consistent with the usual quadrature rules, the finely sampled
$k$-distribution should not be weight-averaged over the final $g$-bins, but
sampled at the precise abscissa of the $g$-point determined by the quadrature
rule used, following Irwin et al. (2008).

As long as we resample on the original grid, we keep the original quadrature
with $N_{g}$ points, and the weights are left unchanged. However, it is still
possible at this point to use different $g$-points, especially if we need to
use less $g$-points for computational reasons (as is often the case with 3D
climate models). For this, we can simply resample
$\tilde{k}_{\mathrm{j}}(\tilde{g}_{\mathrm{j}}^{\bar{\mathrm{b}}})$ on a
different $g$-grid and use the relevant weights, which is an option of
`Exo_k`. In general, `Exo_k` works for any choice of $g$-space sampling
specified by the user.

###  3.4 Non-integer binning

Up to now, we only considered the case in which the final super-band
encompasses an integer number of sub-bands. This is required for our
demonstration to be exact without any further assumptions.

In practice, it is still possible to bin $k$-coefficients on a spectral grid
that does not exactly respect the limits of the sub-bands (hereafter, non-
integer binning). To do that, we replace the spectral extent,
$\Delta\nu^{\mathrm{b}}$, of the two bands at the boundaries of the super-band
in Eq. (7) by their spectral extent inside the super-band. This is equivalent
to binning a high-resolution spectrum to a low-resolution grid using a
weighted average. In theory, this would remain perfectly accurate if the
statistical distribution of the absorption coefficient were constant within
the sub-band. As this is not true in general, we test the accuracy of this
approach in the next section.

##  4 Validation

To validate our approach, we compute spectra in the two cases detailed in
Sect. 2. As we do not want to test the accuracy of the correlated-$k$ approach
itself, which has already been done before, but whether the binning introduces
additional errors, the comparison proceeds as follows:

  * $\bullet$

Starting from very high-resolution spectra, we compute a first reference set
of tables of $k$-coefficients with resolutions ranging from $R=1\,000$ to
$R=5$ and a $g$-space sampled using 20 Gauss-Legendre quadrature points. The
exact $R$ values used do not matter. However, we make sure that the wavenumber
points in the low-resolution grids can be found in the high-resolution grid to
use only natural binning.

  * $\bullet$

Starting from the $R=1\,000$ table of $k$-coefficients444This choice is
motivated by the resolution chosen by the ExoMol project for their
$k$-coefficient tables (Chubb et al., 2020) we just computed, we use the
binning approach described above to compute a new set of ten tables of
$k$-coefficients with the exact same resolutions and wavenumber grids as the
reference set.

  * $\bullet$

The emitted or transmitted flux of our fiducial atmosphere is computed at each
resolution with both the reference and the binned $k$-coefficients and the
results are compared two by two.

The results of the comparison are shown in Fig. 6. For the Mars emission case,
we compare the RMS standard deviation of the relative difference between the
two spectra in the 1-200$\,\mu$m range. For the Hot-jupiter in transmission
test case, the comparison metric used is the RMS standard deviation of the
relative difference between the two spectra in the 1-2$\,\mu$m window.

All the $k$-coefficients computed in this work use 20 Gauss-Legendre
quadrature points in $g$-space. Because the number of intermediate points used
to resample the $g$-distribution in opacity space ($N_{k}$; see Sect. 3.3) is
a free parameter, we tried different values of this parameter. We see that an
insufficient intermediate sampling results in significant numerical errors, as
could be anticipated. However, when $N_{k}$ is a factor of a few larger than
the initial number of $g$-points, the accuracy of the method asymptotes below
$10^{-3}-10^{-4}$, which is itself much lower than the error introduced by the
correlated-$k$ method in the first place (see for example Amundsen et al.
2017). We note that although choosing a larger $N_{k}$ may slightly increase
the binning time, this does not affect the efficiency of the subsequent
radiative transfer calculations in any way. We thus used results with
$N_{k}=5N_{g}$ as a baseline. We also found that using initial data with 40
Gauss-Legendre points did not significantly affect the performances of the
binning procedure. Our algorithm is thus robust to the initial choice of
$g$-space sampling.

The spectra obtained with this set of binned-down $k$-coefficients are also
compared directly to the reference high-resolution spectrum as described in
Sect. 2 to give the absolute difference shown in gray with stars in Fig. 1.
Although not strictly equal because they are computed using a slightly
different reference spectrum, the absolute difference shown in Fig. 1 and the
relative difference from Fig. 6 can be shown to be consistent by multiplying
the latter by the average transit depth as shown on the right axis of Fig. 6.
In this particular case, we further wanted to test the impact of non-integer
binning. For this a new set of $k$-coefficients were computed in the exact
same way, but with a spectral grid that is slightly shifted in wavenumber to
force non-integer binning. As shown by the dotted gray curve in Fig. 1, the
accuracy is lower, although this method remains more accurate than sampled
cross sections. In addition, the error decreases with the binning factor
because the errors at the boundaries of our spectral bins are more diluted.

We therefore conclude that our binning method does not introduce any
significant additional errors.

Figure 6: Numerical relative error introduced by the binning process as a
function of the binning factor (ratio of the initial resolution,
$R_{\mathrm{ini}}$, to the final one, $R_{\mathrm{fin}}$) for four different
values of the number of intermediate sampling points $N_{k}$ (See Sect. 3.3;
dotted: $N_{k}=N_{g}$; dashed: $N_{k}=1.5N_{g}$; dot-dash: $N_{k}=2N_{g}$;
solid: $N_{k}=5N_{g}$, where $N_{g}=20$). Top: The RMS relative top-of-
atmosphere flux error in our Mars-like case. Bottom: The RMS relative transit
depth error in our hot-Jupiter case. The right axis gives the absolute error
to be compared to the results of Figs. 1 and 2. With a sufficient sampling,
the errors are under the $10^{-3}-10^{-4}$ level depending on the geometry,
which is below the error introduced by the correlated-$k$ method in the first
place.

##  5 Conclusions

We have demonstrated mathematically how the $k$-distribution of the opacity in
a large band can be obtained from the $k$-distributions inside each of the
smaller bands that it encompasses. Then we have presented a numerical
algorithm to use this property to bin down precomputed tables of
$k$-coefficients to any arbitrary resolution. Finally, we showed, in two
concrete cases, that the numerical error added by this procedure is relatively
small compared to the errors due to the more general use of the correlated-$k$
approach.

To facilitate the implementation of this new method, we developed `Exo_k`, an
open-source Python library that has been designed to handle many different
sources of radiative data (including, of course, $k$-coefficient tables, but
also cross section and collision-induced-absorption tables). Because this
library is still rapidly evolving, see our extensive online documentation for
tutorials, example notebooks, and a complete description of the library’s
features${}^{\ref{exok_foot}}$.

We think the flexibility that this approach brings to the correlated-$k$
method opens up a completely new way of disseminating reference radiative data
for use in atmospheric models: precomputed $k$-coefficients can now be
directly distributed at a fraction of the data volume while providing fast and
accurate results for any user that does not need a higher resolution. For
example, $R=1000$ $k$-coefficients computed on 20 gauss points such as that
distributed by the ExoMol project take about 50 times less space on disk than
the high-resolution spectra used to produce them. This is a considerable
compression factor, and that does not include the hundreds of CPU hours that
would be needed to compute the spectra.

This of course makes us wonder what would be the optimum resolution and
quadrature for a flexible and versatile correlated-$k$ database. Our test show
that $R=1000$ $k$-coefficients computed on 20 gauss points, when combined with
the flexibility brought by a tool like `Exo_k`, seem to be sufficient for most
1D-3D atmospheric model applications555Very cold and low pressure atmospheres,
such as Pluto’s, may of course require special care because of their very
narrow lines, which usually call for tailored quadrature rules. as well as
modeling synthetic spectroscopic observations of planetary atmospheres up to
that resolution. This method can thus be used to model HST, Spitzer, Ariel
(Tinetti et al., 2018), and low-resolution JWST data. For high-resolution JWST
observation modes, higher resolution radiative data will be needed.

###### Acknowledgements.

I thank F. Selsis for being, as always, an outstanding scientific sparring
partner and the TauREx team at UCL for pointing out many important Python
tricks and libraries. I am also indebted to K. Chubb and M. Turbet for making
some radiative data available before their publication. This project has
received funding from the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation programme (grant agreement n∘
679030/WHIPLASH).

## References

  * Al-Refaie et al. (2019) Al-Refaie, A. F., Changeat, Q., Waldmann, I. P., & Tinetti, G. 2019, arXiv e-prints, arXiv:1912.07759
  * Amundsen et al. (2017) Amundsen, D. S., Tremblin, P., Manners, J., Baraffe, I., & Mayne, N. J. 2017, A&A, 598, A97
  * Chubb et al. (2020) Chubb, K. L., Rocchetto, M., Yurchenko, S. N., et al. 2020, arXiv e-prints, arXiv:2009.00687
  * Garland & Irwin (2019) Garland, R. & Irwin, P. G. J. 2019, arXiv e-prints, arXiv:1903.03997
  * Irwin et al. (2008) Irwin, P. G. J., Teanby, N. A., de Kok, R., et al. 2008, J. Quant. Spec. Radiat. Transf., 109, 1136
  * Kempton et al. (2017) Kempton, E. M. R., Lupu, R., Owusu-Asare, A., Slough, P., & Cale, B. 2017, PASP, 129, 044402
  * Lacis & Oinas (1991) Lacis, A. A. & Oinas, V. 1991, J. Geophys. Res., 96, 9027
  * Leconte et al. (2013) Leconte, J., Forget, F., Charnay, B., et al. 2013, A&A, 554, A69
  * Line et al. (2015) Line, M. R., Teske, J., Burningham, B., Fortney, J. J., & Marley, M. S. 2015, ApJ, 807, 183
  * Liou (1980) Liou, K. N. 1980, An introduction to atmospheric radiation. (Academic Press)
  * Min et al. (2020) Min, M., Ormel, C. W., Chubb, K., Helling, C., & Kawashima, Y. 2020, A&A, 642, A28
  * Mollière et al. (2019) Mollière, P., Wardenier, J. P., van Boekel, R., et al. 2019, A&A, 627, A67
  * Polyansky et al. (2018) Polyansky, O. L., Kyuberis, A. A., Zobov, N. F., et al. 2018, MNRAS, 480, 2597
  * Richard et al. (2011) Richard, C., Gordon, I., Rothman, L., et al. 2011, Journal of Quantitative Spectroscopy and Radiative Transfer,
  * Showman et al. (2009) Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009, ApJ, 699, 564
  * Tennyson & Yurchenko (2012) Tennyson, J. & Yurchenko, S. N. 2012, MNRAS, 425, 21
  * Tinetti et al. (2018) Tinetti, G., Drossart, P., Eccleston, P., et al. 2018, Experimental Astronomy, 46, 135
  * Turbet & Tran (2017) Turbet, M. & Tran, H. 2017, Journal of Geophysical Research (Planets), 122, 2362
  * Waldmann et al. (2015) Waldmann, I. P., Tinetti, G., Rocchetto, M., et al. 2015, ApJ, 802, 107
  * Wordsworth et al. (2011) Wordsworth, R. D., Forget, F., Selsis, F., et al. 2011, ApJ, 733, L48
  * Zhang et al. (2020) Zhang, M., Chachan, Y., Kempton, E. M. R., Knutson, H. A., & Chang, W. H. 2020, ApJ, 899, 27