Open Access
Volume 614, June 2018
Article Number A131
Number of page(s) 12
Section Numerical methods and codes
Published online 26 June 2018

© ESO 2018

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1 Introduction

We present a Fortran 2003 program EXOCROSS to compute spectra as well as spectral properties of molecules using line lists. EXOCROSS is specifically developed to work with huge molecular line lists such as those generated as part of our ExoMol project (Tennyson & Yurchenko 2012) or similar endeavours (Rey et al. 2016). EXOCROSS takes such line lists as input and returns pressure- and temperature-dependent cross sections as well a variety of other derived molecular properties which depend on the underlying spectroscopic data. These include state-dependent lifetimes, temperature-dependent cooling functions, and thermodynamic properties such as partition functions and specific heats.

The main challenge when working with hot line lists for polyatomic molecules is their extremely large sizes. Thus, for example, there are several line lists generated as part of the ExoMol project containing in excess of 10 billion transitions (Yurchenko & Tennyson 2014; Sousa-Silva et al. 2015; Underwood et al. 2016a; Yurchenko et al. 2017a; Owens et al. 2017; Pavlyuchko et al. 2015; Al-Refaie et al. 2015a,b). The size of these datasets makes them impractical for direct use in line-by-line applications. We note that simply ignoring the billions of often very weak lines does not give reliable results (Yurchenko et al. 2014, 2017a). While there are a number of approaches to this problem such as the use of k-coefficients (see, for example, Showman et al. 2009; Amundsen et al. 2014; Malik et al. 2017; Min 2017), the most practical approach which does not involve making significant approximations is to produce cross sections for a set of predefined conditions. These cross sections are then easier to handle in, for example, radiative transfer codes than the original line lists as they can be stored on far fewer grid points than there are lines. However handling these large line lists requires care and, in particular, the generation of cross sections on an appropriate temperature-, pressure- and frequency/wavelength-dependent grid is data intensive and can become computationally highly demanding. EXOCROSS provides a computational solution to this problem; it has been extensively optimised to process huge datasets, including the introduction of an efficient algorithm for generating large numbers of Voigt profiles which is discussed below. EXOCROSS is optimized to provide high throughput via efficient parallelization and vectorization. This is especially important when working with line lists containing tens of billions lines. At different stages of development EXOCROSS was used to generate spectra by Underwood et al. (2016a,b); McKemmish et al. (2016); Wong et al. (2017); Tennyson & Yurchenko (2017); Yurchenko et al. (2017a, 2018); Owens et al. (2017); Prajapat et al. (2017) and Rutkowski et al. (2018).

EXOCROSS is designed to generate molecular cross sections (absorption or emission) on a grid for a set of temperatures and pressures using different line profiles (e.g. Doppler, Voigt, etc.) under the local thermodynamic equilibrium (LTE) as well as non-LTE (Darby-Lewis et al. 2018). Other useful functionality include computing lifetimes (Tennyson et al. 2016a), stick spectra, partition functions, cooling functions, and specific heats. The HITRAN molecular spectroscopic database (Gordon et al. 2017) is a widely-used compilation aimed at radiative transport studies of the Earth’s atmosphere. EXOCROSS is capable of working with HITRAN line lists (.par) as well as super-lines (Rey et al. 2016; Yurchenko et al. 2017a). It can be easily extended to accept other formats.

As part of this implementation, we have developed two new algorithms to perform convolution integrals needed for the Voigt line profile. The first algorithm is based on the Gauss–Hermite quadratures and is developed specifically to guarantee conservation of the Voigt line area. The second algorithm is based on exploiting the similarity of the Voigt profile at large distances from the line centre to compute the opacities quickly.

There are a number of other similar programs available which are designed to work with line lists. These include the HITRAN interface HAPI (Kochanov et al. 2016), SPECTRAPLOT.COM (Goldenstein et al. 2017) and SPECTRA (Tennyson et al. 1993). However, all of these programs would struggle to handle the huge line lists required for models of atmospheres at elevated temperatures. EXOCROSS is designed to be flexible; it takes input in both ExoMol (Tennyson et al. 2013, 2016b) and HITRAN (Rothman et al. 2005) formats. Data can be returned in a variety of formats: ExoMol, HITRAN and Phoenix (Jack et al. 2009), where Phoenix is a full non-LTE atmospheric transfer code accounting for depth-dependent abundances (cloud formation, element diffusion, etc.) using the line-by-line approach. Thus as a subsidiary function the code can be used to interconvert between ExoMol and HITRAN formats.

The paper is organised as follows. The main functionality of EXOCROSS is presented in Sect. 2. The line profile implemented in EXOCROSS are discussed in Sect. 3. Section 4 presents EXOCROSS calculation steps. The data format are described in Sect. 5. Section 7 offers some conclusions. The EXOCROSS manual provided as part of the supplementary data as well as GitHub and CCPForge repositories gives full working details of the program so the description below is restricted to outlines and examples.

2 Main functionality

2.1 Intensities and partition function

An absorption line intensity Ifi (cm molecule−1), also known as absorption coefficient, is given by (1)

where Afi is the Einstein-A coefficient (s−1), is the transition wavenumber (cm−1), Q(T) is the partition function defined as (2)

is the total degeneracy given by

is the nuclear-spin statistical weight factor, c2 = hckB is the second radiation constant (cm K), i = Eihc is the energy term value (cm−1), and T is the temperature in K.

The emissivity (erg (molecule sr)−1) is given by: (3)

Note that the isotopic abundance is not included in the definition of the line intensities (absorption or emission) in Eqs. (1) and (3). This is different from the HITRAN convention, where the absorption coefficients of an isotopologue contain the corresponding natural (terrestrial) isotopic abundances, see For such applications where the isotopic abundance is required, the intensities in Eqs. (1) and (3) can be scaled by an abundance factor specified in the input.

thumbnail Fig. 1

Lifetimes of CH4 computed using the 10to10 line list (Yurchenko & Tennyson 2014). The colors range from dark blue (J = 0) to red (J = 45). See Tennyson et al. (2016a) for a full discussion of methane lifetimes.

2.2 Radiative lifetime

The radiative lifetime (s) can be computed as (Tennyson et al. 2016a) (4)

See an example of the lifetimes in Fig. 1 computed from the 10to10 line list for CH4. Examples of ExoMol lifetimes and cooling functions can be found in Tennyson et al. (2016a); Melnikov et al. (2016) and Mizus et al. (2017).

2.3 Cooling function

The emissivity (erg (s sr molecule)−1) can be used to produce the cooling function W(T) as the total energy emitted by a molecule (Neale et al. 1996) (5)

2.4 Stickspectra

A stick spectrum is a list of frequencies and line intensities, accompanied by the full description (quantum numbers) of the upper and lower states. When plotted, each line is represented by a “stick” with the intensity given by its height, see Table 1 where an extract from an output file containing an absorption stick spectrum of KCl (Barton et al. 2014) is shown. A stick spectra of CaO is shown in Fig. 2.

2.5 Cross sections

A cross section from a single line fi is related to the corresponding integrated absorption coefficient Ifi as (6)

where is a transitional wavenumber. By introducing a line profile the cross section (cm/(molecule cm−1)) can be defined as (7)

where is an integrable function with the area normalized to unity: (8)

Figure 3 shows an example of cross sections of H2 S at T = 300 and 2000 K using the ExoMol line list of Azzam et al. (2016).

Table 1

Extract from a stick spectrum output generated using the KCl line list of Barton et al. (2014).

thumbnail Fig. 2

Stick spectra (cm molecule−1) of CaO (Yurchenko et al. 2016a) compared the CDMS (Müller et al. 2005) rotational band at T = 298 K.

thumbnail Fig. 3

Absorption spectrum of H2S at T = 300 and 2000 K simulated using the ExoMol line list AYT2 (Azzam et al. 2016).

2.6 Grids

By default EXOCROSS uses an equidistant grid, defined by the wavenumber of wavelength range and the number of the grid points Npoints. The latter includes both the first and last bounds. The grid bin size is defined by (9)

The number of intervals is then Npoints − 1. Usually the number of points is an odd number in order to make a “round” value.

Non-equidistant wavenumbers grids can be generated either as grids of constant resolving power or equidistant wavelength grids.

2.7 Partition function and specific heat

The partition function Q(T) is given by Eq. (2). The evaluation of Q(T) requires the energy term values i and degeneracies gtot, which are usually included in molecular line lists. As part of the intensity calculations, the partition function must be either evaluated using these quantities, or directly provided as part of the input. These values can be, e.g. taken from the .pf files provided as part of the ExoMol database (Tennyson et al. 2016b) or as part of the TIPS program provided by HITRAN (Gamache et al. 2017). The direct input option is recommended as often the ExoMol or HITRAN partition functions are more accurate as they contain additional, higher energy contributions which make an important contribution, particularly at elevated temperatures.

The molar specific heat is given by (J K−1 mol−1) (10)

where R is the gas constant and the 1st and 2nd moments Q′ and Q″ are

These latter moments can be also requested from EXOCROSS. An example of Cp (T) of CH4 generated using the 10to10 line list is shown in Fig. 4.

It is often instructive to plot individual contributions to the partition function from different J states defined as (11)

This is useful to assess the convergence of the line list with respect to J and thus to estimate Tmax the line listis applicable to. Figure 5 shows the such individual QJ (T) contributions for the UYT2 line list for SO3 (Underwood et al. 2016a).

thumbnail Fig. 4

Specific heat Cp(T) of CH4 computed using the 10to10 line list Yurchenko & Tennyson (2014).

thumbnail Fig. 5

Contributions QJ(T) to the partition function of SO3 using the line list UYT2 of Underwood et al. (2016a).

2.8 Intensity thresholds

An intensity threshold can be used to speed up the cross-section calculation or to reduce the output in stick-spectra type calculations (done by simply specifying a constant intensity threshold value in cm molecule−1 in the input file). The constant intensity cut-offs are however known to cause problems at long wavelengths, where the density of lines is small and each line, even weak, can be important. A more sophisticated method is to use the dynamic HITRAN’s intensity cut-off (Rothman et al. 2013), defined as

where the HITRAN values for and Icrit are 2000 cm−1 and 10−29 cm molecule−1, respectively.These values are also default in EXOCROSS but can be changed in the input.


EXOCROSS can be used to work with the line list in the HITRAN native format .par, which covers almost all its functionality. It can also be used to convert to ExoMol to HITRAN format (see Sect. 5.2).

2.10 Phoenix

EXOCROSS has the facility to output data in Phoenix format (Jack et al. 2009). In order to speed up the line-by-line calculations Phoenix’s atomic and molecular line lists have a compact structure, where all required properties (line positions, oscillator strengths, lower state energies and broadening parameters) are stored as 4- and 2-bytes integers. For the wavelength (μm, 4 byte-integers) this is defined as:


The oscillator strength gffi for a fi transition, energy term values , and broadening parameters γ and n are mapped onto 2-byte integers according to

where p is one of these properties. The integers iλ, iγ and in are then written as unformatted records with direct access, each of which containing data for 65 536 lines (block-size). For molecules the broadening parameters include the reference Voigt line widths due to H2 () and He (γHe) and the corresponding temperature exponents and nHe (see below). It should be noted that Phoenix uses the so-called astrophysics-convention for the nuclear statistical weights, which are related to the physics convention (adopted by ExoMol and HITRAN) as follows: (12)

where i counts different nuclear statistics. For example, in case of water (HO), the nuclear statistics factors (physics convention) are 1 (para) and 3 (ortho), thus in the astrophysics convention are 1∕4 (para) and 3∕4 (ortho). Since Phoenix’s partition functions are directly affected by the astrophysics convention, in order to be consistent, the ExoMol gffi values have to be scaled by the factor 1∕4 for water, or in general.

2.11 Treating non-local thermodynamic equilibrium (non-LTE)

EXOCROSS provides a simple approach to treating non-LTE environment by differentiating between the rotational and vibrational (vibronic) temperatures when calculating intensities or partition functions (or other T-dependent properties). To this end we approximate the total energy as a sum of the vibrational (or vibronic) and rotational contributions; (13)

where v and k are generic vibrational (vibronic) and rotational quantum numbers, respectively. If the pure vibronic contributions are taken as the corresponding energy values at J = 0 (integer spin), J = 1∕2 (non-integer spin) or the lowest J allowed by the symmetry of the electronic term and the parity, corresponding to the lowest states (usually “+” or “e”). The rotational contribution is simply given by (14)

We also assume that the rotational and vibrational modes are in corresponding (Boltzmann) LTE and that the non-LTE population of a given state (used in intensity and/or partition function calculations) is given by

For this representation it is important to have all the vibrational and rotational quantum numbers defined in the line list, or at least for states accessed by non-LTE calculations.

3 Line profiles

The line broadening is important for practical applications. While temperature effects are commonly modelled by a Doppler line profile, pressure broadening is more complicated.For very high pressure regimes Lorentzian profiles can be used, while for moderate pressures Voigt profiles are generally used (see, for example, Schreier 2017).

3.1 Standard line profiles and sampling method

The most commonly used line profiles in EXOCROSS include Gaussian, Doppler, Voigt and Lorentzian.

The general Gaussian line profile is given by (Hill et al. 2013b) (15)

where is the line centre position and αD is the Gaussian half-width at half-maximum (HWHM). The Gaussian line profile is useful to model generic spectra represented by lines with constant HWHM. The Gaussian line profile can be also used to model the microturbulence broadening by choosing αfi appropriately.

The Doppler line profile is based on the Gaussian shape defined in Eq. (15), where the Doppler HWHM is given by (16)

at temperature T for a molecule of mass m.

The Lorentzian profile is given by (17)

where γL is the Lorentzian line width (HWHM), given most commonly by (18)

Here T0 and P0 are the reference temperature and pressure, respectively, γ0 and nL are broadening parameters for a given broadener, reference HWHM and temperature exponent, respectively.

The Voigt profile is a convolution of the Doppler and Lorentazian profiles: (19)

where and . The Lorentzian line width γL strongly depends on the molecule and is usually also state-dependent. The corresponding values must be given in the input including the specification of the broadeners and their mixing ratio. Each calculation can handle only one combination of broadeners.

Additionally, a simple box-type line profile given by (20)

where is the width of the box, is available.

The individual contribution from each line to the cross sections at a given frequency grid point k is evaluated by sampling the corresponding line profile (see Eq. (7)) a given by

which will be often referred to as a sampling method. This method has the disadvantage of underestimating the opacity when too coarse grids are used which can lead to lines being partially or completely left out. This is a typical problem for long wavelengths where the lines are narrow and far from each other, which is usually tackled either by re-normalizing the line area, see, for example, Sharp & Burrows (2007), or by using a random sampling (Lupu et al. 2016). Below we explore a different, more rigorous alternative.

In practical applications the cross sections are computed on a grid of frequencies (wavenumbers) . When the grid is not sufficiently dense, the line profiles lose their normalisation. This is usually not a problem, at least for most of the room temperature applications. However for high T when billions of lines are used, this leakage can lead to a significant loss of opacity. In order to prevent this effect, Hill et al. (2013b) suggested using an averaged intensity over a given frequency bin, where the corresponding cross section is integrated analytically. This method originally presented for the Gaussian (Doppler) line profile, is extended here to describing Lorentzian and Voigt profiles.

3.2 Binned Gaussian profile with analytical integrals

An averaged (integrated) cross section over a bin from a line fi is given by

where erf is the error function and (23)

are the scaled limits of the wavenumber bin centred on relative to the line centre, , and Ifi is the line intensity in units of cm−1/molecule cm−2 from Eq. (1). Here we take advantage of the fact that an analytical solution exists for the integral of the Gaussian function (24)

The total cross section at the frequency bin k is given by a sum over all contributions from individual lines fi: (25)

and can be interpreted as an average value of the cross sections from a given frequency bin k. The advantage of this approach is that in definition it always gives exact integrated cross sections independent of the number of grid points used or the integration interval. Therefore it is recommended for applications where accurate integrated cross sections or absorption coefficients on coarse grids are required. However it is known that averaged cross sections, especially on coarse grids, can lead to huge errors in integrated flux. Therefore for radiative transfer applications, the direct sampling methods are more accurate and should be used instead.

3.3 Binned Lorentzian profile with analytical integrals

Here we apply the same idea of analytical integral to the Lorentzian line profile:

where (28)

Here the following integral was used:

Again, the integration within each bin is done analytically which guarantees no loss of accuracy for any number of points.

3.4 Binned Voigt profile with analytical integrals

The two line profiles (Gaussian and Lorentzian) can be combined to produce a similar formulation for the Voigt profile, where we use the idea of Gauss–Hermite quadratures as, for example, used in Humlíček’s algorithm (Humlicek 1979). The Voigt convolution integral in Eq. (19) can be written using these quadratures as follows: (29)

where νk and are the Gauss–Hermite quadrature points and weights, respectively (k = 1, …, NG-H). In this form the computation of Voigt can be also generalised to produce the area-conserved integrals using Eq. (27): (30)

We usually take NG-H = 30 Gauss–Hermite points. This approach does not appear to have been taken previously.

thumbnail Fig. 6

Absolute relative error at 4 cm−1 for the Voigt profile at against cm−1 for H2 O with T = 5000 K and γL computed from Eq. (18) with parameters , T0 = 296.0 K, nL = 0.5, P0 = 1 and pressureP at 10−20, 100 and 101 bar.

3.5 Vectorized Voigt approximation

Evaluation of Voigt line profile is generally one of the biggest bottlenecks in opacity calculations. Here we present a new approximate cross section algorithm for the Voigt line profile, which leads to efficient vectorization and thus fast calculations. Our approach is based on the observation that the shape of the wings of the Voigt profile (>4 cm−1 from the line centre), at least for Humlíček’s algorithm, is relatively constant over the large variation of as Lorentzian broadening is generally the largest contributor. For example, Fig. 6 shows how the wings of the Voigt profile centred at differ from the wings of other Voigt profiles centred at all other across the entire wavenumber range from 0 to 30 000 cm−1 (computed using Humlíček’s algorithm). As expected, the error grows as the Doppler HWHM (Eq. (16)) increases with transition wavenumber. However this never exceeds than 1% for even the lowest pressure. One of the most interesting observations is that at 100 bar, the relative error is almost the same as the mostly Doppler profile error at 10−20 bar. With higher pressures this error falls significantly to lower than 10−2% and lower temperatures reduces this by orders of magnitudes. It is only around the line centre, which we estimate to be within 4 cm−1, that the variation of the line shape of the Voigt profile is important.

Based on this observation, the Voigt profile can be split into two parts as follows (omitting the indexes αD, γL for simplicity): (31)

where is a reference Voigt profile centred at cm−1:

Here is a parameter that is used to prevent discontinuities at cm−1 when switching between the two profiles and is given by: (32)

This parameter is included for completeness and is generally set to β = 1 for a performance boost as the discontinuities are not visible at most scales for a single transition and invisible once the whole spectrum is considered. For a given set of pressure broadening parameters γL we pre-compute a set of points defining the wings and then simply select a relevant set. Therefore the only place where the real Voigt calculation needs to be done is around the centre. Additionally, (if used) needs to be calculated at the boundary, which completes the evaluation of the given profile.

The algorithm is based on the (Humlicek 1979) approximation for the Voigt profile in Eq. (19), which is themain method used by EXOCROSS. The Humlíček algorithm is called only for the regions within 4 cm−1 from the line centre. Using the conventionally used Lorentz cutoff of 25 cm−1, this means that only up to 8% of the calculation is computationally demanding giving a theoretical speed up of 12.5 times. This is illustrated in Fig. 7, which shows speed up using our vectorized Voigt algorithm when applied to the region of 0.0–300 cm−1 of the BT2 water linelist (Barber et al. 2006) at T = 1900 K and P = 1 bar. The speed up S for N points used to bin the wavenumber grid is defined as: (33)

where is the time required for a standard Humlíček computation on a wavenumber grid N and is the time required using the vectorized Voigt method. The speed up converges to a maximum value of about 11 times compared to the standard Humlíček calculation, close to the predicted maximum speed up.

This procedure is also efficiently vectorized. Firstly, for the inner part (top of Eq. (31)), which is symmetric, only one half is computed. The other half is then merely looped through backwards and applied to the grid, requiring only to multiply by the absorption coefficient (emissivity) and to add to the opacity grid. The second vectorization occurs when dealing with the second part of Eq. (31). Here again, only a multiply by the intensity and add to the opacity grid is required. These two loops are vectorized through the Fused-Multiply-Add (FMA) instruction.

Figure 8 presents an illustration, where both the vectorized and standard (Humlíček) Voigt methods were used to generate cross sections of water from the BT2 line list for T = 1900 K and P = 1 bar. The new algorithm captures all features with the total opacity for the range shown differing by only by 10−6 cm2 molecule−1.

Lastly, for a full opacity calculation between 0.0–30 000 cm−1, Table 2 shows that using no intensity threshold with the vectorized-Voigt method is almost 3.5 times faster than the full Humlíček method at 10−30 cm molecule−1 thresholding. Comparing like for like, the vectorized Voigt is around 10 to 12 times faster compared to the standard Humlíček method.

Future development of the algorithm will look into automatically tuning the distance from the line centre depending on the temperatureand pressure parameters given.

thumbnail Fig. 7

Speed up, Eq. (33), for the vectorized Voigt method against the standard Voigt (Humlíček) computed on varying wavenumber grid sizes (N) using the BT2 (Barber et al. 2006) water line list computed at T = 1900 K, P = 1 bar and wavenumber range between 0 and 300 cm−1.

thumbnail Fig. 8

Topplot: comparison of cross section calculations for the BT2 water line list (Barber et al. 2006) between the standard Humlíček and the vectorized Voigt method with T = 1900 K. Bottom plot: percentage difference between the Humlíček and vectorized Voigt method. The calculations used no intensity threshold and a wavenumber bin of 0.1 cm−1.

Table 2

Time taken (s) for differing methods and intensity thresholds (cm molecule−1) to compute opacities using the 500 million transitions of the BT2 water line list (Barber et al. 2006) between 0 and 30 000 cm−1 with a wavenumber binning of 0.1 cm−1.

3.6 Binned vectorized Voigt with the line area preserved

Considering the importance of preserving integrated cross section in many applications, we also provide an alternative version of the vectorized Voigt, basedon re-normalization of the line area. During the precomputation stage of the vectorized Voigt method, the total sum for all points () that lie above cm−1 is computed and stored alongside the reference Voigt profile. When computing the vectorized Voigt on a transition, the central Humlíček region is evaluated into a temporary array and its sum is added to . After which the scaled absolute intensity Ĩfi is computedas: (34)

Both the temporary Humlíček array and reference Voigt is applied to the opacity grid with the scaled intensity Ĩfi.

Whilst not a proper treatment of area conservation as that given by Eq. (30), it serves as a reasonable approximation and, as shown in Table 3, gives good results within 1% of the total summed absolute intensity for even large wavenumber bins. To our knowledge, this method does not appear to be reported before.

Table 3

Percentage relative error between the total summed absolute intensity and the total integrated intensity for BT2 water line list computed between 0 and 300 cm−1 at T = 1900 K and P = 1 bar at various wavenumber binnings.

3.7 Broadening parameters

The Voigt profile as a convolution of Doppler and Lorentzian profiles requires definition of the corresponding line widths (HWHM), αD (see Eq. (16)) and γL, given by Eq. (18). The Doppler parameter αD(T) is easy to deal with. It does not depend on the molecular states, only the line position and can be always computed on the fly. The Lorentzian (Voigt) parameters γ0(P0, T0) and nL however are very different for different molecules. Besides they show a pronounced dependence on the state quantum numbers, with the rotational (J) state dependence being the strongest.

The two-file format of the ExoMol database requires special structure for the broadening parameters. Instead of using the conventional line-by-line approach employed by spectroscopic databases such in HITRAN (Gordon et al. 2017) or GEISA (Jacquinet-Husson et al. 2016), where the pressure broadening are specified for the each transitions, ExoMol’s broadening parameters are stored in separate files with the extension .broad (Tennyson et al. 2016b). This structure is justified for most applications as the same parameters are usually used for a large number of different transitions. The latter is either due to the absence of broadening information on all the lines or due to the weak dependence of these parameters for different states. This structure was recently implemented for a number of molecules including H2 O, CH4 and HCN (Barton et al. 2017; Yurchenko et al. 2017b). Table 4 shows an extract from the .broad file for CS as an example. Each line in .broad has the following structure: type (a0, a1, …), γ0 (P0, T0), nL and quantum numbers defined by the type.

Currently EXOCROSS supports three following broadening schemes, constant, a0 and a1, depending on the rigorous quantum numbers J′ and J″. The simplest case is when γ0(P0, T0) and nL are constant and the .broad data is not required. The a0 type corresponds to the J-dependence only. In this case the 4th column in the .broad file contains the J values. The J quantum number is a mandatory quantity in the ExoMol format (Col. 4 in .states) and is therefore relatively straightforward to handle. A similar scenario (a1) is when the broadening depends on the upper J′ (Col. 5 in .broad) and lower J″ (Col. 4) rotational quantum numbers. All other broadening schemes involve dependence on some non-rigorous quantum numbers (‘labels’), such as vibrational v or rotational K. The non-rigorous quantum numbers and their position in the .states file are molecule dependent and thus need to be specified. This information can be found in the ExoMol’s .def (API) file. The current version of EXOCROSS supports rigorous quantum numbers only and therefore does not require interfacing with the ExoMol database.

3.8 Mixtures of broadeners

We consider different broadeners to be independent and their effect additive. Thus the total value of γL is a weighed sum of from each broadener as given by:

where ρi is the fraction portion of the ith broadener. Here we used the fact that the cross sections from each lines are additive and thus the line profile can be represented as a weighted average of lines broadened by different species.

3.9 Off-set

Even though, at least in principle, a line profile has infinite spread, in practical calculations a frequency (or wavelength) cut-offs must be applied to limit the calculation region to around the line centre only. Not only does this influence the computation time and the accuracy of cross sections, but it is also assumed in some applicationsas a point of convention. For example, water cross section are conventionally taken to have a 25 cm−1 cutoff, with far-wing contributions outside this region assumed to form part of the so-called water-continuum (Shine et al. 2012). 25 cm−1 is the default cut-off value in EXOCROSS, alternatively it is specified in the input file.

3.10 Super-lines

The super-line approach is an efficient method for describing a molecular broadened continuum originally proposed by Rey et al. (2016) and was recently studied in detail by Yurchenko et al. (2017a). The super-lines are constructed as temperature-dependent intensity histograms as follows (see also detailed discussion by Rey et al. (2016)). We divide the wavenumber range into N frequency bins, each centred around a grid point . For each the total absorption intensity Ik (T) is computed as a sum of absorption line intensities Ifi, as in Eq. (1), from all fi transitions falling into the wavenumber bin at the given temperature T. Each grid point forms a super-line of an artificial transition with an effective absorption intensity Ik (T). The super-line lists are given in a two-column format with pre-computed intensities Ik, in the same format as used to store ExoMol cross-sections (Tennyson et al. 2016b). The filename have the extension .super. The super-lines approach does not require that histograms are of the same widths and can accept non-equidistant grids as well, see below.

The histograms in EXOCROSS can be produced as cross sections using the Bin-option in the input file (see Manual), which is basically just a sum of all intensities within a given bin i. Ones the histograms are computed (in the standard cross section two-column format), they can be treated as normal line lists. In this case the .states file is not needed as all the information has been already included into the line position and intensity. Moreover, since the states-specific information is completely lost from the line characteristics, the state-dependent line profiles can not be used for temperature/pressure broadening. Doppler line profiles require no information on the upper/lower states and are not restricted. However for the Voigt pressure broadening parameters, which usually depend (at least) on J, only constant values of γ0 and nL (see Eq. (18)) can be used in conjunction with super-lines. For this reason the super-lines are recommended for description of featureless continuum produced from the weaker lines only. The stronger lines should be treated as usual, line-by-line.

Table 4

Air .broad file for 12C32S: portion of the file (upper part); field specification (lower part).

3.11 User-defined profiles

New line profiles, see Tennyson et al. (2014) for example, can be easily implemented to EXOCROSS by the user. A detailed description is provided in the manual. The HITRAN option in EXOCROSS can be used as an example.

4 Calculation protocol

The typical EXOCROSS calculation includes the following steps (see Fig. 9):

  • Read input instruction;

  • Read the .states file: energies, quantum numbers and statistical weighs;

  • Compute the partition function (if required);

  • Read N lines with upper/lower IDs and the Einstein coefficient lines from the .trans file;

  • Apply filters;

  • Compute line intensities (absorption coefficients or emissivities, if required).

  • Compute cross sections on a grid of wavenumbers (if required);

  • Compute lifetimes (if required);

  • Compute cooling functions (if required);

  • Print the cross sections (stick spectra, life times, cooling functions) into a separate file;

  • Do time and memory reporting.

thumbnail Fig. 9

EXOCROSS program work-flow.

Table 5

Extract from the states file of the 14N16O line list.

Table 6

Extract from the transitions file of the 14N16O line list.

5 Data formats

EXOCROSS currently takes in input in either ExoMol or HITRAN format. It can provide output in these formats and in the format used by the Phoenix radiative transport code (Jack et al. 2009). These formats are discussed in turn below.

5.1 ExoMol format

A line list is defined as a catalogue of transition frequencies and intensities (Tennyson & Yurchenko 2012). In the basic ExoMol format (Hill et al. 2013b), adopted by EXOCROSS, a line list has a compact structure consisting of two files: “States” and “Transitions”; an example for the list NOname line list for 14N16O (Wong et al. 2017) is given in Tables 5 and 6. The “States” (.states) file contains energy term values supplemented by the running number n, total degeneracy , rotational quantum number Jn (all obligatory fields), other quantum numbers and labels (both rigorous and not rigorous), lifetimes and Landé g-factors. For example for a generic open-shell diatomic molecule, the quantum numbers include υ, Λ, parity (±), Σ, Ω and the electronic state label (e.g. X2Sigma+) (Yurchenko et al. 2016b). The “Transitions” (.trans) file contains three obligatory columns, the upper and lower state indexes nf and ni which are running numbers from the “State” file, and the Einstein coefficient Afi. For the convenience it also sometimes provides the wavenumbers as the Col. 4. The line list in the ExoMol format can be used to simulate absorption or emission spectra for any temperature in a general way.


The current “HITRAN format” is fully specified in Table 1 of the 2004 edition of HITRAN (Rothman et al. 2005). This format, which is also used for the current release of the related high-temperature database HITEMP (Rothman et al. 2010), has been implemented here.

Although the HITRAN format is widely adopted as a de facto standard, we advise some caution before adopting it. The format is rather verbose and can become extremely unwieldy as a means of representing large line lists. The format is highly tuned towards Earth atmosphere application (e.g. in its choice pressure broadening parameters and temperature ranges) and is therefore rather inflexible for other applications. HITRAN themselves have recognised these issues and have introduced their own web-based interface HAPI (Kochanov et al. 2016) to act as front end and to perform data compression. The database itself has moved to an online-version which provides much more flexibility than the 2004 format (Hill et al. 2013a).

thumbnail Fig. 10

Overview of the absorption line intensities of NaH at T = 2000 K computed using the line list of Rivlin et al. (2015).

5.3 Improving data processing

Both the cross-section and intensity steps (see Fig. 9) are OpenMP parallelized. Users can specify the number of processors requested, which is otherwise set to 1 (no parallelization). In order to make reading and processing data from the .trans file more efficient, EXOCROSS reads line transitions in chunks of N lines, not line-by-line. “Caching” these records into RAM allows for the parallelization for both the transition filtering and of the computation of line-profiles. Each thread is given their own version of the opacity grid to perform work independently without the usage of atomic operations or mutex locks. The total opacity grid can be retrieved at the end of the program run combining all threads opacity arrays. This number N is either specified in the input file or estimated based on the memory available on the system (default). The number of processors must be specified in the input as well (see below on the memory handling).

5.4 Filters

EXOCROSS allows the selection of specific bands/states when computing intensities using the “filter” option. The filters are based on the column-numbers containing the corresponding quantum labels of the upper and lower states. For example, the vibrational quantum number v in the NOname line list is given in the column number 10 (see Table 5), which can be used to generate absorption cross section of NO for the overtone band v = 5, i.e. for transitions between v′ = 5 and v″ = 0 of NO, by referring in the input the corresponding values from the Col. 10 (see Manual for details). Another typical example is to generate cross sections for specific electronic bands, see Fig. 10, where an overview of three absorption electronic bands XX, AX, AA of NaH is shown (Rivlin et al. 2015).

The filter-feature will work even if not all states are assigned. According to the ExoMol convention, the string NaN (with any combination of upper an lower cases) is used for missing quantum labels. Thus “NaN” in this case will be effectively used by EXOCROSS’ filter as a quantum label.

5.5 Units

The defaultunits of EXOCROSS are listed in Table 7. Microns (μm) can be optionally used for wavelength as alternative to wavenumber (default). Pressure does not have designated units; it is assumed to have the same units as of the parameter P0 defining the broadening parameter γ, see Eq. (18).

Table 7

Units used by EXOCROSS.

5.6 Memory handling

The program records and controls the memory used at all processors. For proper control, the user is requested to specify the memoryavailable on the machine in Gb or Mb. This number is used, for example, to estimate the number of transitionlines from .trans processed simultaneously. At the end of the program a memory usage report is given.

6 Program repository

The up-to-date version of the EXOCROSS code together with manual and input examples are freely available from the ExoMol website1, CCPForge2 or GitHub3.

7 Conclusion

We present a new Fortran program EXOCROSS to compute different spectroscopic properties of molecules using spectral line lists. The program has being actively used by ExoMol to generate absorption cross sections using the ExoMol line lists available at In order to works with huge sizes of some line lists, EXOCROSS is optimized for efficient usage of parallelism and vectorization. Our new Voigt algorithm (vectorized Voigt) is designed to be fast and accurate.

The program can be easily extended by users with their profiles or other functionality.

We are planning to provide production of k-coefficients as part of EXOCROSS in the future; integrate the API via the ExoMol .def file; reading the partition function from an ExoMol .pf file; implement a non-LTE model, which does not require definition of non-rigorous quantum numbers (see Sect. 2.11).


We would like to acknowledge help of Derek Homeier and France Allard with Phoenix format. We thank Ingo Waldmann and Marco Rocchetto for testing EXOCROSS. EXOCROSS uses the Fortran 90 input parsing module input.f90 supplied by Anthony J. Stone, which is gratefully acknowledged. This work was supported by the UK Science and Technology Research Council (STFC) No. ST/M001334/1 and the COST action MOLIM No. CM1405. This work made extensive use of UCL’s Legion and STFC’s DiRAC@Darwin high performance computing facilities.


  1. Al-Refaie, A. F., Ovsyannikov, R. I., Polyansky, O. L., Yurchenko, S. N., & Tennyson, J. 2015a, J. Mol. Spectr., 318, 84 [NASA ADS] [CrossRef] [Google Scholar]
  2. Al-Refaie, A. F., Yachmenev, A., Tennyson, J., & Yurchenko, S. N. 2015b, MNRAS, 448, 1704 [NASA ADS] [CrossRef] [Google Scholar]
  3. Amundsen, D. S., Baraffe, I., Tremblin, P., et al. 2014, A&A, 564, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Azzam, A. A. A., Tennyson, J., Yurchenko, S. N., & Naumenko, O. V. 2016, MNRAS, 460, 4063 [NASA ADS] [CrossRef] [Google Scholar]
  5. Barber, R. J., Tennyson, J., Harris, G. J., & Tolchenov, R. N. 2006, MNRAS, 368, 1087 [NASA ADS] [CrossRef] [Google Scholar]
  6. Barton, E. J., Chiu, C., Golpayegani, S., et al. 2014, MNRAS, 442, 1821 [NASA ADS] [CrossRef] [Google Scholar]
  7. Barton, E. J., Hill, C., Yurchenko, S. N., et al. 2017, J. Quant. Spectr. Rad. Transf., 187, 453 [NASA ADS] [CrossRef] [Google Scholar]
  8. Darby-Lewis, D., Tennyson, J., Lawson, K. D., et al. 2018, J. Phys. B At. Mol. Opt. Phys., submitted [Google Scholar]
  9. Gamache, R. R., Roller, C., Lopes, E., et al. 2017, J. Quant. Spectr. Rad. Transf., 203, 70 [NASA ADS] [CrossRef] [Google Scholar]
  10. Goldenstein, C. S., Miller, V. A., Spearrin, R. M., & Strand, C. L. 2017, J. Quant. Spectr. Rad. Transf., 200, 249 [CrossRef] [Google Scholar]
  11. Gordon, I. E., Rothman, L. S., Hill, C., et al. 2017, J. Quant. Spectr. Rad. Transf., 203, 3 [Google Scholar]
  12. Hill, C., Gordon, I. E., Rothman, L. S., & Tennyson, J. 2013a, J. Quant. Spectr. Rad. Transf., 130, 51 [CrossRef] [Google Scholar]
  13. Hill, C., Yurchenko, S. N., & Tennyson, J. 2013b, Icarus, 226, 1673 [NASA ADS] [CrossRef] [Google Scholar]
  14. Humlicek, J. 1979, J. Quant. Spectr. Rad. Transf., 21, 309 [CrossRef] [Google Scholar]
  15. Jack, D., Hauschildt, P. H., & Baron, E. 2009, A&A, 502, 1043 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Jacquinet-Husson, N., Armante, R., Scott, N. A., et al. 2016, J. Mol. Spectr., 327, 31 [NASA ADS] [CrossRef] [Google Scholar]
  17. Kochanov, R., Gordon, I., Rothman, L., et al. 2016, in XVIIIth Symposium on High Resolution Molecular Spectroscopy (HighRus-2015), Tomsk, Russia, J. Quant. Spectr. Rad. Transf., 177, 15 [NASA ADS] [CrossRef] [Google Scholar]
  18. Lupu, R. E., Marley, M. S., Lewis, N., et al. 2016, ApJ, 152, 217 [Google Scholar]
  19. Malik, M., Grosheintz, L., Mendonça, J. M., et al. 2017, ApJ, 153, 56 [NASA ADS] [CrossRef] [Google Scholar]
  20. McKemmish, L. K., Yurchenko, S. N., & Tennyson, J. 2016, MNRAS, 463, 771 [NASA ADS] [CrossRef] [Google Scholar]
  21. Melnikov, V. V., Yurchenko, S. N., Tennyson, J., & Jensen, P. 2016, Phys. Chem. Chem. Phys., 18, 26268 [CrossRef] [Google Scholar]
  22. Min, M. 2017, A&A, 607, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Mizus, I. I., Alijah, A., Zobov, N. F., et al. 2017, MNRAS, 468, 1717 [CrossRef] [Google Scholar]
  24. Müller, H. S. P., Schlöder, F., Stutzki, J., & Winnewisser, G. 2005, J. Mol. Struct., 742, 215 [NASA ADS] [CrossRef] [Google Scholar]
  25. Neale, L., Miller, S., & Tennyson, J. 1996, ApJ, 464, 516 [NASA ADS] [CrossRef] [Google Scholar]
  26. Owens, A., Yurchenko, S. N., Yachmenev, A., Thiel, W., & Tennyson, J. 2017, MNRAS, 471, 5025 [CrossRef] [Google Scholar]
  27. Pavlyuchko, A. I., Yurchenko, S. N., & Tennyson, J. 2015, MNRAS, 452, 1702 [CrossRef] [Google Scholar]
  28. Prajapat, L., Jagoda, P., Lodi, L., et al. 2017, MNRAS, 472, 3648 [CrossRef] [Google Scholar]
  29. Rey, M., Nikitin, A. V., Babikov, Y. L., & Tyuterev, V. G. 2016, J. Mol. Spectr., 327, 138 [NASA ADS] [CrossRef] [Google Scholar]
  30. Rivlin, T., Lodi, L., Yurchenko, S. N., Tennyson, J., & Le Roy R. J. 2015, MNRAS, 451, 634 [NASA ADS] [CrossRef] [Google Scholar]
  31. Rothman, L. S., Jacquemart, D., Barbe, A., et al. 2005, J. Quant. Spectr. Rad. Transf., 96, 139 [NASA ADS] [CrossRef] [Google Scholar]
  32. Rothman, L. S., Gordon, I. E., Barber, R. J., et al. 2010, J. Quant. Spectr. Rad. Transf., 111, 2139 [NASA ADS] [CrossRef] [Google Scholar]
  33. Rothman, L. S., Gordon, I. E., Babikov, Y., et al. 2013, J. Quant. Spectr. Rad. Transf., 130, 4 [Google Scholar]
  34. Rutkowski, L., Foltynowicz, A., Johansson, A. C., et al. 2018, J. Quant. Spectr. Rad. Transf., 205, 213 [NASA ADS] [CrossRef] [Google Scholar]
  35. Schreier, F. 2017, J. Quant. Spectr. Rad. Transf., 187, 44 [NASA ADS] [CrossRef] [Google Scholar]
  36. Sharp, C. M., & Burrows, A. 2007, ApJS, 168, 140 [NASA ADS] [CrossRef] [Google Scholar]
  37. Shine, K. P., Ptashnik, I. V., & Rädel, G. 2012, Surv. Geophys., 33, 535 [NASA ADS] [CrossRef] [Google Scholar]
  38. Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009, ApJ, 699, 564 [NASA ADS] [CrossRef] [Google Scholar]
  39. Sousa-Silva, C., Al-Refaie, A. F., Tennyson, J., & Yurchenko, S. N. 2015, MNRAS, 446, 2337 [Google Scholar]
  40. Tennyson, J., & Yurchenko, S. N. 2012, MNRAS, 425, 21 [Google Scholar]
  41. Tennyson, J., & Yurchenko, S. N. 2017, Mol. Astrophys., 8, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Tennyson, J., Miller, S., & Le Sueur C. R. 1993, Comput. Phys. Commun., 75, 339 [NASA ADS] [CrossRef] [Google Scholar]
  43. Tennyson, J., Hill, C., & Yurchenko, S. N. 2013, in 6th International Conference on Atomic and Molecular Data and Their Applications ICAMDATA-2012 (New York: AIP), AIP Conf. Proc., 1545, 186 [NASA ADS] [CrossRef] [Google Scholar]
  44. Tennyson, J., Bernath, P. F., Campargue, A., et al. 2014, Pure Appl. Chem., 86, 1931 [Google Scholar]
  45. Tennyson, J., Hulme, K., Naim, O. K., & Yurchenko, S. N. 2016a, J. Phys. B At. Mol. Opt. Phys., 49, 044002 [NASA ADS] [CrossRef] [Google Scholar]
  46. Tennyson, J., Yurchenko, S. N., Al-Refaie, A. F., et al. 2016b, J. Mol. Spectr., 327, 73 [NASA ADS] [CrossRef] [Google Scholar]
  47. Underwood, D. S., Tennyson, J., Yurchenko, S. N., Clausen, S., & Fateev, A. 2016a, MNRAS, 462, 4300 [NASA ADS] [CrossRef] [Google Scholar]
  48. Underwood, D. S., Tennyson, J., Yurchenko, S. N., et al. 2016b, MNRAS, 459, 3890 [NASA ADS] [CrossRef] [Google Scholar]
  49. Wong, A., Yurchenko, S. N., Bernath, P., et al. 2017, MNRAS, 470, 882 [NASA ADS] [CrossRef] [Google Scholar]
  50. Yurchenko, S. N., & Tennyson, J. 2014, MNRAS, 440, 1649 [NASA ADS] [CrossRef] [Google Scholar]
  51. Yurchenko, S. N., Tennyson, J., Bailey, J., Hollis, M. D. J., & Tinetti, G. 2014, PNAS, 111, 9379 [NASA ADS] [CrossRef] [Google Scholar]
  52. Yurchenko, S. N., Blissett, A., Asari, U., et al. 2016a, MNRAS, 456, 4524 [NASA ADS] [CrossRef] [Google Scholar]
  53. Yurchenko, S. N., Lodi, L., Tennyson, J., & Stolyarov, A. V. 2016b, Comput. Phys. Commun., 202, 262 [NASA ADS] [CrossRef] [Google Scholar]
  54. Yurchenko, S. N., Amundsen, D. S., Tennyson, J., & Waldmann, I. P. 2017a, A&A, 605, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Yurchenko, S. N., Tennyson, J., & Barton, E. J. 2017b, J. Phys. Conf. Ser., 810, 012010 [CrossRef] [Google Scholar]
  56. Yurchenko, S. N., Sinden, F., Lodi, L., et al. 2018, MNRAS, 473, 5324 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1

Extract from a stick spectrum output generated using the KCl line list of Barton et al. (2014).

Table 2

Time taken (s) for differing methods and intensity thresholds (cm molecule−1) to compute opacities using the 500 million transitions of the BT2 water line list (Barber et al. 2006) between 0 and 30 000 cm−1 with a wavenumber binning of 0.1 cm−1.

Table 3

Percentage relative error between the total summed absolute intensity and the total integrated intensity for BT2 water line list computed between 0 and 300 cm−1 at T = 1900 K and P = 1 bar at various wavenumber binnings.

Table 4

Air .broad file for 12C32S: portion of the file (upper part); field specification (lower part).

Table 5

Extract from the states file of the 14N16O line list.

Table 6

Extract from the transitions file of the 14N16O line list.

Table 7

Units used by EXOCROSS.

All Figures

thumbnail Fig. 1

Lifetimes of CH4 computed using the 10to10 line list (Yurchenko & Tennyson 2014). The colors range from dark blue (J = 0) to red (J = 45). See Tennyson et al. (2016a) for a full discussion of methane lifetimes.

In the text
thumbnail Fig. 2

Stick spectra (cm molecule−1) of CaO (Yurchenko et al. 2016a) compared the CDMS (Müller et al. 2005) rotational band at T = 298 K.

In the text
thumbnail Fig. 3

Absorption spectrum of H2S at T = 300 and 2000 K simulated using the ExoMol line list AYT2 (Azzam et al. 2016).

In the text
thumbnail Fig. 4

Specific heat Cp(T) of CH4 computed using the 10to10 line list Yurchenko & Tennyson (2014).

In the text
thumbnail Fig. 5

Contributions QJ(T) to the partition function of SO3 using the line list UYT2 of Underwood et al. (2016a).

In the text
thumbnail Fig. 6

Absolute relative error at 4 cm−1 for the Voigt profile at against cm−1 for H2 O with T = 5000 K and γL computed from Eq. (18) with parameters , T0 = 296.0 K, nL = 0.5, P0 = 1 and pressureP at 10−20, 100 and 101 bar.

In the text
thumbnail Fig. 7

Speed up, Eq. (33), for the vectorized Voigt method against the standard Voigt (Humlíček) computed on varying wavenumber grid sizes (N) using the BT2 (Barber et al. 2006) water line list computed at T = 1900 K, P = 1 bar and wavenumber range between 0 and 300 cm−1.

In the text
thumbnail Fig. 8

Topplot: comparison of cross section calculations for the BT2 water line list (Barber et al. 2006) between the standard Humlíček and the vectorized Voigt method with T = 1900 K. Bottom plot: percentage difference between the Humlíček and vectorized Voigt method. The calculations used no intensity threshold and a wavenumber bin of 0.1 cm−1.

In the text
thumbnail Fig. 9

EXOCROSS program work-flow.

In the text
thumbnail Fig. 10

Overview of the absorption line intensities of NaH at T = 2000 K computed using the line list of Rivlin et al. (2015).

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.