© J. Leconte 2020

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 tables¹ (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 ~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², an open-source Python 3 library that can directly load precomputed k-coefficient and cross-section 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 theaccuracy 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:

• Transmission spectrum of a hot Jupiter with a 1000 K isothermal H₂/H₂O atmosphere in front of a Sun.

• Emission spectrum of a pure CO₂ Mars-like atmosphere.

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

Each time, we start by computing a reference spectrum using high-resolution monochromatic cross sections (with a typical step of 0.01 cm⁻¹). For our CO₂ atmosphere, we use spectra from Turbet & Tran (2017). For our hot-Jupiter case, we use water spectra produced by the Exomol project with H₂ 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_sp = {1000, 3000, 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_fin = {10, 20, 50, 100, 200, 500, 1000}). The transmission spectra computed in the hot-Jupiter case are shown in the left column of Fig. 2. Finally, these spectra are compared tothe properly binned reference spectrum (right column of Fig. 2). The root mean squared (RMS) error is shown for each (R_sp, R_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².

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_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_fin = 10 case for which the R_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.

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.

Fig. 1 Root mean square error on the transit spectrum of our fiducial hot Jupiter over the 1–2 μ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: Rsp = {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).

Fig. 2 Effect of sampling resolution on the transit spectrum of a hot Jupiter. The various rows show the effect for different final observing resolutions (Rfin, specified in the top left corner of the panels). Left column: transit depth in the 1–2 μ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; Rsp) shown in the legend. The curve with dots shows the reference case computed from our very high-resolution cross sections (Δσ ≈ 0.01 cm−1). The high-resolution reference spectrum before binning is shown as the shading in the top left panel. Right column: 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.

Fig. 3 Effect of sampling resolution on the emission spectrum of our Mars-like case around the 15 μm CO2 band. The rows show the effect for different final observing resolutions (Rfin). Left column: flux. The various colors correspond to the four sampling resolutions (i.e., the resolution used in radiative transfer calculation before binning; Rsp) shown inthe legend. The reference case is computed from our very high-resolution cross sections (Δσ ≈ 0.01 cm−1). The high-resolution reference spectrum before binning is shown as the shading in the top left panel. Right column: 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.

Fig. 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 opacity.

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 ν is the wavenumber, is divided into N_b spectral bands of width $\Delta {\nu }^{\text{b}}\equiv {\nu }_{+}^{\text{b}}-{\nu }_{-}^{\text{b}}$. In each band, the transmission through a slab of matter can be written as $${\mathcal{T}}^{\text{b}}(u)=\frac{1}{\Delta {\nu }^{\text{b}}}{\displaystyle {\int }_{{\nu }_{-}^{\text{b}}}^{{\nu }_{+}^{\text{b}}}{\text{e}}^{-k^{\nu }u}}\text{d}\nu ={\displaystyle {\int }_0^1{\text{e}}^{-{\widehat{k}}^{\text{b}}(g)u}}\text{d}g\approx {\sum }_{\text{i}=1}^{N_g}{w}_{\text{i}}^{\text{b}}{\text{e}}^{-{\widehat{k}}_{\text{i}}^{\text{b}}u},$$(1)

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

An intuitive way to understand that is to say that now, in band b, the opacity can be described by N_g representative values (${\widehat{k}}_{\text{i}}^{\text{b}}$) for each of the N_g points in g-space (g_i), each of these values being affected a weight (wb_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_{\text{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_{\text{in}}{}^{\nu }=F_{\text{in}}{}^{\text{b}}/\Delta {\nu }^{\text{b}}$, where $F_{\text{in}}{}^{\text{b}}\equiv {\displaystyle {\int }_{{\nu }_{-}^{\text{b}}}^{{\nu }_{+}^{\text{b}}}{F}_{\text{in}}{}^{\nu }}\text{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 $$\begin{array}{ll}{F}_{\text{out}}{}^{\text{b}}(u)\equiv \hfill & {\displaystyle {\int }_{{\nu }_{-}^{\text{b}}}^{{\nu }_{+}^{\text{b}}}{F}_{\text{in}}{}^{\nu }}{\text{e}}^{-k^{\nu }(\nu )u}\text{d}\nu \approx \frac{F_{\text{in}}{}^{\text{b}}}{\Delta {\nu }^{\text{b}}}{\displaystyle {\int }_{{\nu }_{-}^{\text{b}}}^{{\nu }_{+}^{\text{b}}}{\text{e}}^{-k^{\nu }(\nu )u}}\text{d}\nu \hfill \\ =\hfill & F_{\text{in}}{}^{\text{b}}\cdot {\mathcal{T}}^{\text{b}}(u).\hfill \end{array}$$(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 $$F={\sum }_{\text{b}=1}^{N_{\text{b}}}{F}_{\text{out}}{}^{\text{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: Can we determine the k-coefficients ${\widehat{k}}_{\text{i}}^{\overline{\text{b}}}$ of a “super-band” $\overline{\text{b}}$ that is the union of two smaller, nonoverlapping “sub-bands”³, b₁ and b₂, for which we know the k-coefficients, ${\widehat{k}}_{\text{i}}^{{\text{b}}_1}$ and ${\widehat{k}}_{\text{i}}^{{\text{b}}_2}$?

In this form, the answer is not trivial. However, when computing k-coefficients, the crucial mathematical object is notthe k-distribution (${\widehat{k}}^{\text{b}}(g)$) but its inverse, the cumulative density function (CDF) of the opacity within the band (hereafter called the g-distribution, ĝ^b (k)) defined as $${\widehat{g}}^{\text{b}}(k)\equiv \frac{1}{\Delta {\nu }^{\text{b}}}{\displaystyle {\int }_{{\nu }_{-}^{\text{b}}}^{{\nu }_{+}^{\text{b}}}H}\left(k-k^{\nu }(\nu )\right)\text{d}\nu ,$$(4)

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

The question can thus be reframed to ask whether we can compute the g-distribution of the super-band, ${\widehat{g}}^{\overline{\text{b}}}$, if we know the g-distribution within the two sub-bands (ĝ^b₁ and ĝ^b₂). This is rather straightforward, because as long as the two smaller bands do not overlap, we have $\Delta {\nu }^{\overline{\text{b}}}=\Delta {\nu }^{{\text{b}}_1}+\Delta {\nu }^{{\text{b}}_2}$ and $$\begin{array}{ll}{\widehat{g}}^{\overline{\text{b}}}(k)\hfill & \equiv \frac{1}{\Delta {\nu }^{\overline{\text{b}}}}{\displaystyle {\int }_{{\nu }_{-}^{\overline{\text{b}}}}^{{\nu }_{+}^{\overline{\text{b}}}}H}\left(k-k^{\nu }(\nu )\right)\text{d}\nu \hfill \\ \hfill & =\frac{1}{\Delta {\nu }^{\overline{\text{b}}}}\left({\displaystyle {\int }_{{\nu }_{-}^{{\text{b}}_1}}^{{\nu }_{+}^{{\text{b}}_1}}H}\left(k-k^{\nu }(\nu )\right)\text{d}\nu +{\displaystyle {\int }_{{\nu }_{-}^{{\text{b}}_2}}^{{\nu }_{+}^{{\text{b}}_2}}H}\left(k-k^{\nu }(\nu )\right)\text{d}\nu \right)\hfill \\ \hfill & =\frac{\Delta {\nu }^{{\text{b}}_1}}{\Delta {\nu }^{\overline{\text{b}}}}{\widehat{g}}^{{\text{b}}_1}(k)+\frac{\Delta {\nu }^{{\text{b}}_2}}{\Delta {\nu }^{\overline{\text{b}}}}{\widehat{g}}^{{\text{b}}_2}(k).\hfill \end{array}$$(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 ${\widehat{g}}^{\overline{\text{b}}}(k)$, that is, the probability that a monochromatic ray of light falls in a spectral region where the opacity k^ν 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 }^{\text{b}}/\Delta {\nu }^{\overline{\text{b}}}$, as long as the incoming spectral flux is constant over the super-band, multiplied by the probability that k^ν < k within this sub-band, which is ĝ^b(k). So we recover Eq. (5) 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 (${\widehat{k}}_{\text{i}}^{\text{b}}$) that are tabulated for fixed g-points (g_i) with the associated weights (wb_i). The difficulty is that although g_i can be seen as the values of ĝ^b in band b for the points located at ${\widehat{k}}_{\text{i}}^{\text{b}}$ in opacity space, the distribution ĝ^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 numberof bands. For the moment, however, we restrict ourselves to the case in which the super-band $\overline{\text{b}}$ is composed ofan integer number of sub-bands (hereafter, natural binning). Once these sub-bands have been identified, we determine $$k_{\text{min}}=\underset{\text{i},\text{b}\in \overline{\text{b}}}{\text{min}}\left\{{\widehat{k}}_{\text{i}}^{\text{b}}\right\}\text{and}{k}_{\text{max}}=\underset{\text{i},\text{b}\in \overline{\text{b}}}{\text{max}}\left\{{\widehat{k}}_{\text{i}}^{\text{b}}\right\}$$(6)

and compute an evenly log-distributed grid of N_k points (${\tilde{k}}_{\text{j}}$) between these two values. Then, for each sub-band, we interpolate $g_{\text{i}}({\widehat{k}}_{\text{i}}^{\text{b}})$ on this grid, which yields N_k values per band, ${\tilde{g}}_{\text{j}}^{\text{b}}$, representingall the ĝ^b on the same grid for all bands. It is thus now easy to compute the global g-distribution which is given by $${\tilde{g}}_{\text{j}}^{\overline{\text{b}}}={\sum }_{\text{b}\in \overline{\text{b}}}\frac{\Delta {\nu }^{\text{b}}}{\Delta {\nu }^{\overline{\text{b}}}}{\tilde{g}}_{\text{j}}^{\text{b}}.$$(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}}_{\text{j}}({\tilde{g}}_{\text{j}}^{\overline{\text{b}}})$ onto our original g-grid (g_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 longas 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}}_{\text{j}}({\tilde{g}}_{\text{j}}^{\overline{\text{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.

Fig. 5 Example of the binning process for two bands. Left: k-coefficients on Ng = 20 Gauss-Legendre quadrature points for two bands (${\widehat{k}}_{\text{i}}^{{\text{b}}_1}$ and ${\widehat{k}}_{\text{i}}^{{\text{b}}_2}$). Middle: g-distribution for the two bands resampled on ${\tilde{k}}_{\text{j}}$, a grid of Nk = 100 points between kmin and kmax (${\tilde{g}}_{\text{j}}^{{\text{b}}_1}$ and ${\tilde{g}}_{\text{j}}^{{\text{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}}_{\text{j}}^{\overline{\text{b}}}$). Right: colored lines are shown as in the left panel. The black line indicates the finely sampled k-distribution (${\tilde{k}}_{\text{j}}({\tilde{g}}_{\text{j}}^{\overline{\text{b}}})$). The gray line represents the k-distribution of the super-band re-sampled on the original g-grid with Ng = 20, i.e., the final k-coefficients (${\widehat{k}}_{\text{i}}^{\overline{\text{b}}}$).

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, Δν^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:

• Starting from very high-resolution spectra, we compute a first reference set of tables of k-coefficients with resolutions ranging from R = 1000 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.

• Starting from the R = 1000 table of k-coefficients⁴ 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.

• 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 μ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 μ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⁻³−10⁻⁴, 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.

Fig. 6 Numerical relative error introduced by the binning process as a function of the binning factor (ratio of the initial resolution, Rini, to the final one, Rfin) for four different values of the number of intermediate sampling points Nk (see Sect. 3.3; dotted: Nk = Ng; dashed: Nk = 1.5Ng; dot-dash: Nk = 2Ng; solid: Nk = 5Ng, where Ng = 20). Top: RMSrelative top-of-atmosphere flux error in our Mars-like case. Bottom: 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².

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 applications⁵ 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

All Tables

All Figures

Fig. 1 Root mean square error on the transit spectrum of our fiducial hot Jupiter over the 1–2 μ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: Rsp = {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).

In the text

Fig. 2 Effect of sampling resolution on the transit spectrum of a hot Jupiter. The various rows show the effect for different final observing resolutions (Rfin, specified in the top left corner of the panels). Left column: transit depth in the 1–2 μ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; Rsp) shown in the legend. The curve with dots shows the reference case computed from our very high-resolution cross sections (Δσ ≈ 0.01 cm−1). The high-resolution reference spectrum before binning is shown as the shading in the top left panel. Right column: 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.

In the text

Fig. 3 Effect of sampling resolution on the emission spectrum of our Mars-like case around the 15 μm CO2 band. The rows show the effect for different final observing resolutions (Rfin). Left column: flux. The various colors correspond to the four sampling resolutions (i.e., the resolution used in radiative transfer calculation before binning; Rsp) shown inthe legend. The reference case is computed from our very high-resolution cross sections (Δσ ≈ 0.01 cm−1). The high-resolution reference spectrum before binning is shown as the shading in the top left panel. Right column: 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.

In the text

	Fig. 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 opacity.
In the text

Fig. 5 Example of the binning process for two bands. Left: k-coefficients on Ng = 20 Gauss-Legendre quadrature points for two bands (${\widehat{k}}_{\text{i}}^{{\text{b}}_1}$ and ${\widehat{k}}_{\text{i}}^{{\text{b}}_2}$). Middle: g-distribution for the two bands resampled on ${\tilde{k}}_{\text{j}}$, a grid of Nk = 100 points between kmin and kmax (${\tilde{g}}_{\text{j}}^{{\text{b}}_1}$ and ${\tilde{g}}_{\text{j}}^{{\text{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}}_{\text{j}}^{\overline{\text{b}}}$). Right: colored lines are shown as in the left panel. The black line indicates the finely sampled k-distribution (${\tilde{k}}_{\text{j}}({\tilde{g}}_{\text{j}}^{\overline{\text{b}}})$). The gray line represents the k-distribution of the super-band re-sampled on the original g-grid with Ng = 20, i.e., the final k-coefficients (${\widehat{k}}_{\text{i}}^{\overline{\text{b}}}$).

In the text

Fig. 6 Numerical relative error introduced by the binning process as a function of the binning factor (ratio of the initial resolution, Rini, to the final one, Rfin) for four different values of the number of intermediate sampling points Nk (see Sect. 3.3; dotted: Nk = Ng; dashed: Nk = 1.5Ng; dot-dash: Nk = 2Ng; solid: Nk = 5Ng, where Ng = 20). Top: RMSrelative top-of-atmosphere flux error in our Mars-like case. Bottom: 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.

In the text