© The Authors 2024

1 Introduction

Hydrogen deuteride (HD), the second most abundant isotopolog of molecular hydrogen, is an important tracer of deuterium in the Universe. The small and constant primordial fraction of deuterium to hydrogen, D/H ((2.8 ± 0.2) × 10⁻⁵; Pettini et al. 2008), is one of the key arguments supporting the Big Bang theory. Measurements of the D/H ratio in the Solar System provide information about planetary formation and evolution. The standard ratio on Earth is considered to be Vienna Standard Mean Ocean Water (VSMOW; D/H = 1.5576 × 10⁻⁴; Araguás-Araguás et al. 1998). Employing it, Donahue et al. (1982) found the D/H ratio in the Venusian atmosphere to be (1.6 ± 0.2) × 10⁻² (i.e., two orders of magnitude higher than VSMOW). This higher ratio is attributed to the evaporation of oceans and the subsequent photodissociation of H₂O in the upper parts of the atmosphere (Donahue & Pollack 1983). A comparison of the D/H ratios on Earth and Mars with those determined for various comets indicates the role of this ratio in the volatile accretion on these two planets (Drake & Righter 2002; Hartogh et al. 2011). The two largest gas giants, Jupiter and Saturn, are incapable of nuclear fusion; thus, planetary models predict that the D/H ratio in their atmospheres should be close to the primordial value for the Solar System (Lellouch et al. 2001) or slightly larger due to the accretion of deuterium-rich icy grains and planetesimals (Guillot 1999). The abundance of deuterium is larger by a factor of 2.5 in the atmospheres of Uranus and Neptune (Feuchtgruber et al. 2013), owing to their ice-rich interiors (Guillot 1999). Additionally, the determination of the D/H ratio in comets and moons provides information about the formation of ices in the early Solar System (Hersant et al. 2001; Gautier & Hersant 2005; Horner et al. 2006).

The D/H ratio in planets, moons, or comets can be derived from either in situ mass spectrometry (Eberhardt et al. 1995; Mahaffy et al. 1998; Niemann et al. 2005; Altwegg et al. 2015) or spectroscopic observations of various molecules and their deuterated isotopologs in the millimeter and infrared range (Bockelée-Morvan et al. 1998; Meier et al. 1998; Crovisier et al. 2004; Fletcher et al. 2009; Pierel et al. 2017; Krasnopolsky et al. 2013; Blake et al. 2021). Interestingly, as pointed out by Krasnopolsky et al. (2013), the ratios determined from spectra of different molecules can differ substantially. Moreover, even on Earth, the HDO/H₂O ratio in the atmosphere can vary substantially across the globe (Araguás-Araguás et al. 1998). In this context, the observation of isotopologs of molecular hydrogen can be a reliable benchmark for determining D/H ratios. The most accurate values of the D/H ratio in gas giants stem from analyses of the pure rotational R(0), R(1), R(2), and R(3) lines of HD and S(0) and S(1) lines in H₂. Lellouch et al. (2001) determined the D/H ratio in the Jovian atmosphere using the Short Wavelength Spectrometer (SWS) on board the Infrared Space Observatory to be (2.25 ± 0.35) × 10⁻⁵. Pierel et al. (2017) analyzed the far-infrared spectra of Saturn’s atmosphere gathered by Cassini's Composite Infrared Spectrometer. Interestingly, the Pierel et al. result suggests that the D/H ratio on Saturn is lower than that of Jupiter, which contradicts predictions based on interior models (Guillot 1999; Owen & Encrenaz 2006) and points to an unknown mechanism of deuterium fractionating in Saturn’s atmosphere. The D/H ratio in the atmospheres of Uranus and Neptune was determined from measurements of the pure rotational R(0), R(1), and R(2) lines in HD using the Photoconductor Array Camera and Spectrometer (PACS) on board the Herschel space observatory (Feuchtgruber et al. 2013). The analysis revealed similar values for these two giants - (4.4 ± 0.4) × 10⁻⁵ and (4.1 ± 0.4) × 10⁻⁵, respectively – confirming the expected deuterium enrichment with respect to the protosolar value.

Atmospheric models of the Solar System’s gas giants, from which the relative abundance of HD with respect to H₂ is retrieved (and consequently, the D/H ratio), require knowledge about the line parameters of HD and H₂. The temperature profiles considered in the models cover seven orders of magnitude of pressure (1 µbar−10 bar; see for instance Pierel et al. 2017). Thus, in addition to line position and intensity, knowledge about collisional effects that perturb the shape of observed lines is crucial for the accurate determination of the relative HD abundance. Furthermore, the incorporation of non-Voigt line-shape effects (such as Dicke narrowing and speed-dependent effects) allows for reducing the systematic errors in atmospheric models, as shown for Jupiter (Smith 1989) and for Uranus and Neptune (Baines et al. 1995). Indeed, although the spectral lines of most molecules observed in planetary atmospheres are not sensitive to non-Voigt effects, considering the resolving power of available telescopes, the lines of molecular hydrogen and its isotopologs are (Smith 1989; Baines et al. 1995).

In this article we report accurate collision-induced line-shape parameters for the rotational transitions R(0), R(1), and R(2) within the ground vibrational state. These HD lines are most frequently used for studies of the atmospheres of the giant planets. The results cover all thermodynamic conditions relevant to the atmospheres of giant planets in the Solar System, that is, all relevant temperatures, pressures, and H₂−He perturbing gas compositions (the HD-H₂ data are provided in this work, while the HD-He data are taken from Stankiewicz et al. 2020, 2021). Besides the basic Voigt-profile collisional parameters (i.e., the broadening and shift parameters), we also report their speed dependences and the complex Dicke parameter, which, as we will show, can influence the effective width and height of the HD lines up to almost a factor of 2 for giant planet conditions.

For the R(0) line, the non-Voigt regime coincides with the maximum of the monochromatic contribution function (see Fig. 1 in Feuchtgruber et al. 2013), which has a direct influence on the abundance of HD inferred from observations.

We performed quantum-scattering calculations on our new, highly accurate ab initio 6D potential energy surface (PES), and we used the scattering S matrices obtained in this way to determine the collision-induced line-shape parameters. We used the cavity ring-down spectroscopy to validate our theoretical methodology, demonstrating a sub-percent-level accuracy that considerably surpasses the accuracy of any previous theoretical (Schaefer & Monchick 1992) or experimental study of line-shape parameters of pure rotational lines in HD (Ulivi et al. 1989; Lu et al. 1993; Sung et al. 2022) and offers valuable input for the HITRAN database (Gordon et al. 2022). This work represents significant methodological and computational progress; calculations at this level of theory and accurate experimental validation have already been performed for a molecule-atom system (3D PES; Słowiñski et al. 2020, 2022), but in this work we extend it to a molecule-molecule system (6D PES). All the reported parameters (and their temperature dependences) are consistent with the HITRAN database format, hence allowing for the use of the HITRAN Application Programming Interface (HAPI; Kochanov et al. 2016) for generating the beyond-Voigt spectra of HD for any H₂−He perturbing gas composition and thermodynamic condition.

2 Ab initio calculations of the line-shape parameters

In recent years, the methodology for accurate ab initio calculations of the line-shape parameters (including the beyond-Voigt parameters) was developed and experimentally tested for He-perturbed H₂ and HD rovibrational lines, starting from accurate ab initio H₂-He PES calculations (Bakr et al. 2013; Thibault et al. 2017), through the state-of-the-art quantum scattering calculations and line-shape parameter determination (Thibault et al. 2017; Jóźwiak et al. 2018), up to accurate experimental validation (Słowiñski et al. 2020, 2022) and using the results for populating the HITRAN database (Wcisło et al. 2021; Stankiewicz et al. 2021).

In this work, we extend the entire methodology to a much more complex system of a diatomic molecule colliding with another diatomic molecule. First, we calculated an accurate 6D H₂-H₂ PES. Second, we performed state-of-the-art quantum-scattering calculations. Third, we calculated the full set of the six line-shape parameters in a wide temperature range. Finally, we report the results in a format consistent with the HITRAN database.

The 6D PES was obtained using the supermolecular approach based on the level of theory similar to that used to calculate the 4D H₂-H₂ surface (Patkowski et al. 2008). The crucial contributions involve: (1) interaction energy calculated at the Hartree–Fock (HF) level, (2) the correlation contribution to the interaction energy calculated using the coupled-cluster method with up to perturbative triple excitations, CCSD(T), with the results extrapolated to the complete basis set limit (Halkier et al. 1999), (3) electron correlation effects beyond CCSD(T) up to full configuration interaction (FCI), and (4) the diagonal Born– Oppenheimer correction (DBOC; Handy et al. 1986). The details regarding the basis sets used in calculations of each contribution and the analytical fit to the interaction energies are given in Appendix A. The PES is expected to be valid for intramolecular distances r_i ∈ [0.85, 2.25] a₀.

For the purpose of performing quantum scattering calculations, the 6D PES was expanded over a set of appropriate angular functions and the resulting 3D numerical function in radial coordinates was then expanded in terms of rovibrational wave functions of isolated molecules (see Appendix B for details). The close-coupling equations were solved in the body-fixed frame using a renormalized Numerov’s algorithm, for the total number of 3 014 energies ($E_{\text{T}}=E_{\text{kin }}+E_{v_1{j}_1}+E_{v_2{j}_j}$, where E_kin is the relative kinetic energy of the colliding pair, and $E_{v_1{j}_1}$ and $E_{v_2{j}_2}$ are the rovibrational energies of the two molecules at large separations) in a range from 10⁻³ cm⁻¹ to 4000 cm⁻¹. Calculations were performed using the quantum scattering code from the BIGOS package developed in our group (Jóźwiak et al. 2024; Jóźwiak 2024). The scattering S matrix elements were obtained from boundary conditions imposed on the radial scattering function. Convergence of the calculated S matrix elements was ensured by a proper choice of the integration range, propagator step, and the size of the rovibrational basis (see Appendix C for details).

Next, we calculated the generalized spectroscopic cross-sections, ${\sigma }_{\lambda }^q$ (Monchick & Hunter 1986; Schaefer & Monchick 1992), which describe how collisions perturb the shape of molecular resonance. Contrary to the state-to-state cross-sections, which give the rate coefficients (see, for instance, Wan et al. 2019), the ${\sigma }_{\lambda }^q$ cross-sections are complex. For λ = 0, real and imaginary parts of this cross-section correspond to the pressure broadening and shift cross-section, respectively. For λ = 1, the complex cross-section describes the collisional perturbation of the translational motion and is crucial for the proper description of the Dicke effect. The index q is the tensor rank of the spectral transition operator and equals 1 for electric dipole lines considered here.

We used the ${\sigma }_0^1$ and ${\sigma }_1^1$ cross-sections to calculate the six line-shape parameters relevant to collision-perturbed HD spectra, the collisional broadening, and the shift,

$${\gamma }_0-i{\delta }_0=\frac{1}{2\pi c}\frac{1}{k_{\text{B}}T}\langle {\upsilon }_{\text{r}}\rangle {\displaystyle {\int }_0^{\infty }x}{\text{e}}^{-x}{\sigma }_0^1(x)\text{d}x,$$(1)

the speed dependences of collisional broadening and shift,

$$\begin{array}{l}{\gamma }_2-i{\delta }_2=\frac{1}{2\pi c}\frac{1}{k_{\text{B}}T}\frac{\langle {\upsilon }_{\text{r}}\rangle \sqrt{M_{\text{a}}}}{2}{\text{e}}^{-y^2}\hfill \\ \times {\displaystyle \underset{0}{\overset{\infty }{\int }}\left(2\overline{x}\cosh (2\overline{x}y)-\left(\frac{1}{y}+2y\right)\sinh (2\overline{x}y)\right)}{\overline{x}}^2{\text{e}}^{-{\overline{x}}^2}{\sigma }_0^1\left(\overline{x}{\overline{\upsilon }}_{\text{p}}\right)\text{d}\overline{x},\hfill \end{array}$$(2)

and the real and imaginary parts of the complex Dicke parameter,

$${\tilde{v}}_{\text{opt}}^r-i{\tilde{v}}_{\text{opt}}^i=\frac{1}{2\pi c}\frac{\langle {\upsilon }_{\text{r}}\rangle M_{\text{a}}}{k_{\text{B}}T}{\displaystyle \underset{0}{\overset{\infty }{\int }}x}{\text{e}}^{-x}\left[\frac{2}{3}x{\sigma }_1^1\left(x{k}_{\text{B}}T\right)-{\sigma }_0^1\left(x{k}_{\text{B}}T\right)\right]\text{d}x,$$(3)

where υ_r is the relative (absorber to perturber) speed of the colliding molecules, 〈υ_r〉 is its mean value at temperature T, ${\overline{v}}_{\text{p}}$ is the most probable speed of the perturbed distribution, $M_{\text{a}}=\frac{m_{\text{a}}}{m_{\text{a}}+m_{\text{p}}}$, $x=\frac{E_{\text{kin}}}{k_{\text{B}}T}$, $\overline{x}=\frac{2v_{\text{r}}}{\sqrt{\pi M_{\text{a}}}{v}_{\text{r}}}$, $y=\sqrt{\frac{m_{\text{p}}}{m_{\text{a}}}}$ and m_a and m_p are the masses of the active and perturbing molecules, respectively (Wcisło et al. 2021). We estimate the uncertainty of the calculated line-shape parameters in Appendix C. The six line-shape parameters define the modified Hartmann–Tran (mHT) profile (see Appendix D), which encapsulates the relevant beyond-Voigt effects. To make the outcome of this work consistent with the HITRAN database (Gordon et al. 2022), we provide temperature dependences of the calculated line-shape parameters within the double-power-law (DPL) format (Gamache & Vispoel 2018; Stolarczyk et al. 2020):

$$\begin{array}{l}{\gamma }_0(T)=g_0{\left(T_{\text{ref }}/T\right)}^n+{g^{\prime }}_0{\left(T_{\text{ref }}/T\right)}^{n^{\prime }},\hfill \\ {\delta }_0(T)=d_0{\left(T_{\text{ref}}/T\right)}^m+{d^{\prime }}_0{\left(T_{\text{ref}}/T\right)}^{m^{\prime }}\text{, }\hfill \\ {\gamma }_2(T)=g_2{\left(T_{\text{ref }}/T\right)}^j+{g^{\prime }}_2{\left(T_{\text{ref }}/T\right)}^{j^{\prime }},\hfill \\ {\delta }_2(T)=d_2{\left(T_{\text{ref }}/T\right)}^k+{d^{\prime }}_2{\left(T_{\text{ref }}/T\right)}^{k^{\prime }},\hfill \\ {\tilde{v}}_{\text{opt}}^r(T)=r{\left(T_{\text{ref}}/T\right)}^p+r^{\prime }{\left(T_{\text{ref}}/T\right)}^{p^{\prime }}\hfill \\ {\tilde{v}}_{\text{opt}}^i(T)=i{\left(T_{\text{ref }}/T\right)}^q+i^{\prime }{\left(T_{\text{ref }}/T\right)}^{q^{\prime }},\hfill \end{array}$$(4)

where T_ref = 296 K.

3 Results: Line-shape parameters for the R(0), R(1), and R(2) lines in HD perturbed by a mixture of H₂ and He

In Fig. 1, we show the main result of this work: all six line-shape parameters for the R(0), R(1), and R(2) 0–0 lines in H₂-perturbed HD calculated as a function of temperature. Figure 1 covers the temperature range relevant for the giant planets, 50 to 200 K (Lellouch et al. 2001; Feuchtgruber et al. 2013; Pierel et al. 2017; for our full temperature range, 20 to 1000 K, see the data available at CDS). In Fig. 1, we also recall the corresponding He-perturbed data calculated with the same methodology at the same accuracy level (Stankiewicz et al. 2021). The difference between the two perturbers is not negligible, and for many cases the line-shape parameters differ by a factor of 2 or even more.

The data shown in Fig. 1 are given in Table 1 in a numerical form within the HITRAN DPL format (see Eq. (4)). The accuracy of our ab initio line-shape parameters is within 1% of the magnitude of each parameter (see Appendix C for details). The DPL approximation of the temperature dependences introduces additional errors. For γ₀(T) and ${\tilde{v}}_{\text{opt}}^r(T)$, the DPL error is negligible, for (δ₀(T) the DPL error is at the 1 % level, and for other line-shape parameters it can be even higher, but their impact on the final line profile is much smaller (see Appendix E for details). For applications that require the full accuracy of our ab initio line-shape parameters, we provide the line-shape parameter values explicitly on a dense temperature grid in the data available at CDS.

The set of parameters in Table 1 contains all the information necessary to simulate the collision-perturbed shapes of the three HD lines at a high level of accuracy at any conditions relevant to the atmospheres of giant planets (pressure, temperature, and He/H₂ relative concentration). In Fig. 2, we show an example of simulated spectra based on the data from Table 1. It should be emphasized that at the conditions relevant to giant planets, the shapes of the HD lines may considerably deviate from the simple Voigt profile. In Fig. 2a, we show the difference between the Voigt profile and a more physical profile, which includes the relevant beyond-Voigt effects such as Dicke narrowing and the dependence on the speed of the broadening and shift (the mHT profile; see Appendix D). For the moderate pressures, the error introduced by the Voigt-profile approximation can reach almost 70%. The orange, yellow, red, and black lines in Fig. 2a are the temperature-pressure profiles for Jupiter (Seiff et al. 1998), Saturn (Lindal et al. 1985), Uranus (Lindal et al. 1987), and Neptune (Lindal et al. 1990), respectively. We illustrate this with spectra simulations for the case of Neptune’s atmosphere (panels b–d corresponds to the black dots in Fig. 2a). Panels b–d illustrate three different regimes. Panel b is the low-pressure regime with a small collisional contribution in which the line shape collapses to Gaussian, and hence the beyond-Voigt effects are small. Panel d is the opposite: the shape of a resonance is dominated by the collisional effects, but the simple Lorentzian broadening dominates over other collisional effects and again the beyond-Voigt effects are small. Panel c illustrates the nontrivial situation in which the shape of resonance is dominated by the beyond-Voigt effects (see the green horizontal ridges in Fig. 2a). The discrepancies between the beyond-Voigt line-shape model and the Voigt profile are also clearly seen as a difference between the blue lines and blue shadows in panel c. In the context of the giant planet studies, it should be noted that the beyond-Voigt regions marked in panel a (the green horizontal ridges) coincide well with the maxima of the monochromatic contribution functions for these three HD lines (see Fig. 1 in Feuchtgruber et al. 2013 for the case of the atmospheres of Neptune and Uranus).

Panels b–d of Fig. 2 also illustrate the influence of atmosphere composition on the collision-induced shapes of the HD lines for the example of Neptune atmosphere. The differences between the He- and H₂-perturbers are negligible in the low-pressure regime (panel b) since at these conditions the line shape is mainly determined by thermal Doppler broadening. In moderate- and high-pressure ranges (panels c and d), the profiles differ at the peak by a factor of 2. Hence, including both perturbing species is important for spectra analyses of the atmospheres of giant planets, especially for the R(0) line, whose contribution function dominates at moderate and high pressures (Feuchtgruber et al. 2013).

The data reported in this article (see Table 1) account for three factors that are necessary for reaching sub-percent-level accuracy: (1) separate ab initio data for both perturbers (which allows one to simulate perturbation by any H₂-He mixture), (2) accurate representation of temperature dependences, and (3) parametrization of the beyond-Voigt line-shape effects. In general, simulating the beyond-Voigt line-shape profiles is a complex task (see Appendix D). In this work, we used HAPI (Kochanov et al. 2016) to generate the beyond-Voigt spectra shown in Fig. 2 (based on the DPL parameters from Table 1):

from hapi import *
db_begin (’hitran_data’)
nu,coef = absorptionCoefficient_mHT(
       SourceTables=’HD’,
       Diluent={’He’:0.2,’H2’:0.8},
       WavenumberRange=[xmin,xmax],
       WavenumberStep=step,
       Environment={’p’:press,’T’:temp},
       HITRAN_units=True)

The combination of the data reported in Table 1 and the Python-based HAPI constitutes a powerful tool that allows one to efficiently generate accurate HD spectra (based on advanced beyond-Voigt model, the mHT profile) for arbitrary temperature, pressure, and mixture composition.

At low temperatures, relevant for studies of giant planet atmospheres and the chemistry and dynamics of the interstellar medium and protoplanetary disks, the spin isomer (para/ortho) concentration ratio of H₂ at thermal equilibrium (eq-H₂) deviates from 1:3 (the ratio of so-called normal H₂, n-H₂). Moreover, various processes, such as diffusion between atmospheric layers in gas giants, might result in the sub-equilibrium distribution of H₂. These nontrivial para/ortho distributions play a key role in atmospheric models that involve collision-induced absorption (Karman et al. 2019) and spectral features originating from hydrogen dimers (Fletcher et al. 2018), as well as in isotope chemistry of the interstellar medium, where para/ortho ratio controls the deuterium fractionation process (Flower et al. 2006; Nomura et al. 2022). In Fig. 3, we show the influence of the spin isomer concentration on the line-shape parameters. Spin isomer concentration has a large impact at low temperatures. All the line-shape parameters reported in this work are calculated for the thermal equilibrium spin isomer concentration.

Fig. 1 Ab initio temperature dependences of the collisional line-shape parameters (in units of 10−3 cm−1 atm−1 ) of the first three electric dipole lines of HD perturbed by H2 (red curves) and He (black curves).

Fig. 2 Impact of the beyond-Voigt effects and bath mixture composition on collision-perturbed spectra of HD at conditions relevant for giant planet atmospheres. Panel a: relative error of the Voigt-profile approximation as a function of pressure and temperature, shown as the relative difference between the Voigt and mHT profiles at profile maximum. The panels show, from top to bottom, the R(0), R(1), and R(2) lines. Panels b–d: simulations of the HD spectra (blue lines) at conditions relevant for the Neptune atmosphere (the perturbing bath is 80% H2 and 20% He). The spectra are generated with the mHT profile using HAPI based on the DPL temperature parametrization. As a reference, we show the same lines for the cases of pure H2 and pure He perturbers (see the red and black lines, respectively). The blue shadows show the same simulations as the blue lines but generated with the simple Voigt profile. Panels b–d correspond to points b–d shown in the temperature-pressure maps in panel a (the three selected points lie on the Neptune temperature-pressure line). The three cases illustrate three different line-shape regimes. The first one, (b), is the low-pressure case in which the lines are broadened mainly by the Doppler effect, and the pressure-induced collisional effects do not dominate the line shapes. The intermediate-pressure case, panel c, illustrates the extreme non-Voigt regime (the differences between the blue curves and blue shadows reach almost a factor of 2; see also the green ridge in the maps in the bottom panel). The third case, panel d, illustrates the high-pressure regime at which the HD lines are well described by a simple Voigt profile (the blue shadows almost overlap with the blue lines), but setting a proper composition of the perturber gas components plays an important role.

4 Experimental validation

In Fig. 4, we show a comparison between our ab initio calculations (black lines) and the experimental data available in the literature. Fourier-transformed scans from the Michelson interferometer were used to obtain the high-pressure spectra reported in the works of Ulivi et al. (1989) and Lu et al. (1993); the spectra were collected in a temperature range from 77 to 296 K. Recently, the same lines were measured at low pressures (< 1 bar) with the Fourier transform spectrometer coupled to the Soleil-synchrotron far-infrared source (Sung et al. 2022) in a temperature range from 98 to 296 K (see the olive lines in Fig. 4). The discrepancy between these experimental data is by far too large to test our theoretical results at the one percent level.

To validate our ab initio calculations at the estimated accuracy level, we performed accurate measurements using a frequency-stabilized cavity ring-down spectrometer linked to an optical frequency comb, referenced to a primary frequency standard (Cygan et al. (2016, 2019); Zaborowski et al. (2020)). Our 73.5-cm-long ultrahigh finesse (ℱ = 637 000) optical cavity operates in the frequency-agile rapid scanning spectroscopy mode (Truong et al. 2013; Cygan et al. 2016, 2019; see Zaborowski et al. 2020 for details regarding the experimental setup). Since our spectrometer operates at 1.6 µm, we chose the S (2) 2−0 line in the H₂-perturbed D₂ (we repeated all the ab initio calculations for this case). From the perspective of theoretical methodology, the H₂-perturbed D₂ and H₂-perturbed HD are equivalent and either can be used for validating the theoretical methodology (for both cases two distinguishable diatomic molecules are considered and the PES is the same except for the almost negligible DBOC term; see Appendix A). We used a sample of 2% D₂ and 98% of H₂ mixture and collected the spectra at four pressures (0.5, 1, 1.5, and 2 atm) and two temperatures (296 and 330 K; see the black dots in Fig. 5). The corresponding theoretical spectra are the red curves. The methodology for simulating the collision-perturbed shapes of molecular lines (based on the line-shape parameters calculated from Eqs. (1)–(3)) is described in our previous works (Wcisło et al. 2018; Słowiñski et al. 2020, 2022). The two sets of residuals depicted in Fig. 5 show comparisons with two line-shape models, the speed-dependent billiard-ball (SDBB) profile and the mHT profile. The SDBB profile (Shapiro et al. 2002; Ciuryło et al. 2002) is the state-of-the-art approach that gives the most realistic description of the underlying collisional processes. As expected it gives the best agreement with experimental spectra (the mean residuals are 0.65%; see Fig. 5), but it is computationally very expensive (Wcisło et al. 2013). The mHT profile is slightly less accurate (the mean residuals are 1.23%) but it is highly efficient from a computational perspective and, hence, well suited for practical spectroscopic applications. In conclusion, the ab initio line-shape parameters reported in this work (Fig. 1 and Table 1) lead to profiles that are in excellent agreement with accurate experimental spectra, and the theory-experiment comparison is limited by a choice of a line-shape model used to simulate the experimental spectra.

Fig. 3 Dependence of pressure broadening, γ0 (solid lines), and shift, δ0 (dashed lines), on the spin isomer concentration ratio (x = np/no) for the H2-perturbed R(0) line in HD at different temperatures. The vertical gray lines correspond to the value of x for normal H2, x = 1/3, and xeq(T), as determined by the Boltzmann distribution at T = 50, 195, and 296 K.

5 Conclusion

We have computed accurate collision-induced line-shape parameters for the three pure rotational HD lines (R(0), R(l), and R(2)) that are currently employed for the analysis of the giant planets’ atmospheres. To this end, we investigated HD–H₂ collisions using coupled channel quantum scattering calculations on a new, highly accurate ab initio PES. Scattering S matrices determined from these calculations allowed us to obtain the collisional width and shift, as well as their speed dependences and the complex Dicke parameter of H₂-perturbed HD lines. By integrating data from our previous work on the HD-He system (Stankiewicz et al. 2020, 2021), we provide comprehensive results that cover a wide range of thermodynamic conditions, including temperature, pressure, and H₂–He concentration, relevant to the atmospheres of giant planets. We validated our theoretical methodology using cavity ring-down spectroscopy, demonstrating a sub-percent-level accuracy that surpasses the accuracy of previous theoretical and experimental studies of line-shape parameters in HD.

All the reported line-shape parameters and their temperature dependences are consistent with the HITRAN database format. Utilizing HAPI, we demonstrated how our results can be applied to simulate HD spectra under various conditions pertinent to giant planets in the Solar System.

Until now, the analysis of observed collision-perturbed spectra in astrophysical studies has predominantly relied on the simple Voigt profile. We have introduced a methodology and provided a comprehensive dataset that enables the simulation of beyond-Voigt shapes for HD in H₂–He atmospheres. Our work demonstrates that accounting for the speed dependence of collisional width and shift, along with the complex Dicke parameter, is crucial. These factors can alter the effective width and height of HD lines by up to a factor of 2. To quantitatively assess the impact of these results on the D/H ratio on giant planets, radiative transfer modeling is necessary. Given that the beyond-Voigt effects result in a narrower line width (Fig. 2c), the previously inferred HD abundance from astrophysical observations may be underestimated.

Fig. 4 Comparison of the experimental and theoretical values of the pressure broadening and shift parameters, γ0 and δ0 (in units of 10−3 cm−1 atm−1 ). Black curves correspond to the ab initio calculations performed in this work, while green and red points report the experimental measurements from Ulivi et al. (1989) and Lu et al. (1993), respectively. The olive curves are the single-power-law (for γ0) and linear (for δ0) temperature dependences retrieved from the measurements of Sung et al. (2022).

Fig. 5 Direct validation of the ab initio quantum-scattering calculations on the accurate experimental spectra of the S(2) 2−0 line of D2 perturbed by collisions with H2 molecules (see the main text for details). The black dots are the experimental spectra, and the red lines are the ab initio profiles. Below each profile, we show the absolute residuals of two models: the SDBB profile and the mHT profile. To quantify how well theory agrees with experiments, we report the relative (with respect to the profile peak value) root mean square errors (rRMSEs) of the experiment-theory differences calculated within the ±FWHM range around the line center (see the percentages below the residuals). The mean rRMSEs are also summarized for each of the models (see the numbers on the right side of the figure).

Acknowledgements

H.J. and N.S. were supported by the National Science Centre in Poland through Project No. 2019/35/B/ST2/01118. H.J. was supported by the Foundation for Polish Science (FNP). K.S contribution is supported by budgetary funds within the Minister of Education and Science program “Perły Nauki”, Project No. PN/01/0196/2022. D.L. was supported by the National Science Centre in Poland through Project No. 2020/39/B/ST2/00719. S.W. was supported by the National Science Centre in Poland through Project No. 2021/42/E/ST2/00152. P.J. was supported by the National Science Centre in Poland through Project No. 2017/25/B/ST4/01300. K.P. was supported by the U.S. National Science Foundation award CHE-1955328. K.Sz. was supported by the US NSF award CHE-2313826. I.E.G.’s contribution was supported through NASA grant 80NSSC24K0080. P.W. was supported by the National Science Centre in Poland through Project No. 2022/46/E/ST2/00282. For the purpose of Open Access, the author has applied a CC-BY public copyright licence to any Author Accepted Manuscript (AAM) version arising from this submission. We gratefully acknowledge Polish high-performance computing infrastructure PLGrid (HPC Centers: ACK Cyfronet AGH, CI TASK) for providing computer facilities and support within the computational grant, Grant No. PLG/2023/016409. Calculations have been carried out using resources provided by the Wroclaw Centre for Networking and Supercomputing (http://wcss.pl), Grant No. 546. The research is a part of the program of the National Laboratory FAMO in Toruń, Poland.

Appendix A Details of the PES calculations

The 6D PES for the H₂-H₂ system was calculated at the following level of theory;

$$E_{\text{int}}=E_{\text{int}}^{\text{HF}}[5]+\delta E_{\text{int}}^{\text{CCSD}(\text{T})}[\text{Q}5]+\delta E_{\text{int}}^{\text{T}(\text{Q})}[\text{Q}]+\delta E_{\text{int}}^{\text{FCl}}[\text{T}]+\delta E_{\text{int}}^{\text{DBOC}}[\text{T}].$$(A.1)

In all cases, the aug-cc-pVXZ basis sets (Kendall et al. 1992) were employed with the cardinal number X taking values 2 (D), 3 (T), 4 (Q), and 5. The consecutive terms are defined in the following way: $E_{\text{int}}^{\text{HF}}$[5] is the interaction energy calculated at the HF level using the aug-cc-pV5Z basis set, $\delta E_{\text{int}}^{\text{CCSD}(\text{T})}$[Q5] is the correlation contribution to the interaction energy calculated using the coupled-cluster method with up to perturbative triple excitations, CCSD(T), with the results extrapolated to complete basis set limits using the 1/X³ formula (Halkier et al. 1999) from the calculations in the aug-cc-pVQZ and aug-cc-pV5Z basis sets. The next contribution, $\delta E_{\text{int }}^{\text{T}(\text{Q})}[\text{Q}]=E_{\text{int}}^{\text{CCSDT}(\text{Q})}[\text{Q}]-E_{\text{int}}^{\text{CCSD}(\text{T})}[\text{Q}]$, accounts for the electron correlation effects beyond CCSD(T) included in the CC method with up to perturbative quadruple excitations, CCSDT(Q), computed with the aug-cc-pVQZ basis set, whereas $\delta E_{\text{int}}^{\text{FCI}}\left[\text{T}\right]=E_{\text{int}}^{\text{FCl}}[\text{T}]-E_{\text{int}}^{\text{CCSDT}(Q)}[\text{T}]$ describes the electron correlation effects beyond CCSDT(Q), calculated with the aug-cc-pVTZ basis set. The DBOC (Handy et al. 1986), $\delta E_{\text{int}}^{\text{DBOC}}[\text{T}]$, was calculated with the masses of ¹H with the CCSD densities (Valeev & Sherrill 2003; Gauss et al. 2006) obtained with the aug-cc-pVTZ basis set. The $\delta E_{\text{int}}^{\text{DBOC}}$ term is the only one that depends on masses. However, it is small compared to other terms and would be even smaller if calculated for HD–H₂ instead of H₂–H₂. Thus, the resulting surface can be applied to the former system with full confidence.

The interaction energies were fitted by an analytic function that consisted of short- and long-range parts, with a smooth switching at R values between 9 and 10 a₀, using the switching function from Babin et al. (2013). The short-range part was taken in the form of a sum of products of exponentials ${\text{e}}^{-\alpha r_{ab}}$, where r_ab are atom-atom distances (Fernández et al. 1999; Braams & Bowman 2009). In contrast to most published work that uses the same α for all terms, we used four different optimized values. The form of the long-range part was taken from Patkowski et al. (2008), but the parameters were multiplied by linear combinations of symmetry-invariant polynomials of r₁, r₂, and of magnitudes of their differences. The linear coefficients are determined from the fit to the ab initio energies obtained at the same level of theory as for the short-range part. The PES is expected to be valid for r_i ∈ [0.85, 2.25] a₀. Zuo et al. (2021) recently published a 6D PES for the H₂ dimer obtained from the complete active space self-consistent field calculations combined with the multi-reference configuration interaction calculations. According to Zuo et al. (2021), the potential is valid up to r_i = 3.45 a₀, beyond the upper limit of our PES. However, in the region of validity of both PESs, our surface should be more accurate due to the higher level of theory used.

Appendix B Quantum scattering calculation details

The dependence of the 6D PES on the three Jacobi angles was separated from the radial and intramolecular distances by the expansion of the PES over the bi-spherical harmonics, $I_{l_1{l}_{2}l}\left({\theta }_1,{\theta }_2,\varphi \right)$:

$$V\left(R,r_1,r_2,{\theta }_1,{\theta }_2,\varphi \right)={\displaystyle \sum _{l_1,l_2,l}{A}_{l_1{l}_{2}l}}\left(R,r_1,r_2\right)I_{l_1{l}_{2}l}\left({\theta }_1,{\theta }_2,\varphi \right),$$(B.1)

where the bi-spherical harmonics are defined as

$$I_{l_1{l}_{2}l}\left({\theta }_1,{\theta }_2,\varphi ={\varphi }_1-{\varphi }_2\right)=\sqrt{\frac{2l+1}{4\pi }}{\displaystyle \sum _m\left(l_{1}m{l}_2-m\mid l_1{l}_{2}l0\right)}{Y}_{l_{1}m}\left({\theta }_1,{\varphi }_1\right)Y_{l_2-m}\left({\theta }_2,{\varphi }_2\right).$$(B.2)

If the i-th (i = 1, 2) molecule is homo-nuclear, the l_i index in the expansion in Eq. (B.1) takes only even values. The $A_{l_1{l}_{2}l}\left(R,r_1,r_2\right)$ expansion coefficients were obtained by integrating the product of the PES and the corresponding bi-spherical harmonic, over Jacobi angles (see, for instance, Eq. (3) in Zadrożny et al. (2022) and the discussion therein). We employed a 19-point Gauss-Legendre quadrature to integrate over θ₁ and θ₂ and a 19-point Simpson's rule for the integral over ϕ. The integration resulted in a tabular representation of the $A_{l_1{l}_{2}l}\left(R,r_1,r_2\right)$ expansion coefficients, calculated for R in the range of 2.5 to 200 a₀ with a step of 0.1 a₀, and for intramolecular distances ranging from 0.85 to 2.25 a₀ with a step of 0.1 a₀.

In this work we used terms up to the I₄₄₈(θ₁, θ₂, ϕ) bi-spherical harmonic, which corresponds to a total of 19 and 32 terms in the D₂-H₂ and HD-H₂ case, respectively. Such numbers of terms represent an intermediate complexity of the problem – the number of terms in the HD-H₂ case is larger by a factor of 4 in comparison to the HD-He case (Stankiewicz et al. (2020, 2021)), but significantly smaller in comparison to more anisotropic PESs, studied in our previous works (85 for O₂-N₂ (Gancewski et al. (2021)), and 205 for CO-N₂ and CO-O₂ (Jóźwiak et al. (2021); Zadrożny et al. (2022))). The error introduced by the truncation of the PES expansion is discussed in Appendix C.

The dependence of the expansion coefficients on r₁ and r₂ was reduced by averaging $A_{l_1{l}_{2}l}\left(R,r_1,r_2\right)$ over rovibrational wave functions of isolated molecules, $\left({\chi }_{{\eta }_i}\left(r_i\right)\right)$:

$$A_{l_1{l}_{2}l,{\eta }_1,{{\eta }^{\prime }}_1,{\eta }_2,{{\eta }^{\prime }}_2}(R)={\displaystyle \int \text{d}}{r}_2{\chi }_{{\eta }_2}\left(r_2\right)\left({\displaystyle \int \text{d}}{r}_1{\chi }_{{\eta }_1}\left(r_1\right)A_{l_1{l}_{2}l}\left(R,r_1,r_2\right){\chi }_{{{\eta }^{\prime }}_1}\left(r_1\right)\right){\chi }_{{}_{{{\eta }^{\prime }}_2}}\left(r_2\right),$$(B.3)

where η_i = (v_i, j_i) denotes the quantum numbers of a rovibrational state of the i-th molecule. The wave functions of H₂, HD, and D₂ were obtained by solving the nuclear Schrödinger equation for isolated molecules using the potential energy curve of Schwenke (1988). We used the standard trapezoidal rule to perform the integration in Eq. (B.3), and we obtain $A_{l_1{l}_2,{\eta }_1,{\eta }_1^{\prime },{\eta }_2,{\eta }_2^{\prime }}(R)$ coupling terms for R within the range 2.5 to 200 a₀, with a step size of 0.1 a₀.

The average in Eq. (B.3) provides a large number of possible coupling terms. In the HD-H₂ case, we considered pure rotational transitions (up to 1000 K); thus, we neglected the terms that couple excited vibrational states ($v_{\text{HD}}^{\prime }\ne v_{\text{HD}},v_{{\text{H}}_2}^{\prime }\ne v_{{\text{H}}_2}$). In the D₂-H₂ case, we observe that the terms that couple different vibrational levels are three orders of magnitude smaller than terms diagonal in v. Since we performed quantum scattering calculations in the $v_{{\text{D}}_2}=0$ and $v_{{\text{D}}_2}=2$ states separately (while maintaining $v_{{\text{H}}_2}=0$), we neglected radial coupling terms off-diagonal in vibrational quantum numbers. This approximation is additionally justified by the fact that the 2-0 S(2) line in H₂-perturbed deuterium is measured at room temperature, where the population of H₂ in vibrationally excited states is negligible.

The dependence of the coupling terms in Eq. (B.3) on rotational quantum numbers (usually at the level of a few percent) is one of the key factors that affect theoretical predictions of the pressure shift of pure rotational lines in light molecules (Shafer & Gordon 1973; Dubernet & Tuckey 1999; Thibault et al. 2016; Jóźwiak et al. 2018). This is due to the fact that the line shift is sensitive to the difference in the scattering amplitude in the two rotational states that participate in an optical transition. Some authors neglect the j dependence of the radial coupling terms (the centrifugal distortion of the PES) in scattering calculations for rovibrational transitions (Green et al. (1989); Thibault et al. (2017)) and average the expansion coefficients for a given vibrational v over the rovibrational wave function v, j = 0. This approximation is invoked either to save computational resources or due to a lack of information about the dependence of the PES on the stretching coordinates but works well for Q(j) lines. We have shown that taking the centrifugal distortion of the PES into account is crucial for achieving a sub-percent agreement with the experimental spectra in He-perturbed vibrational lines in H₂ (Słowiñski et al. 2022) and HD (Stankiewicz et al. 2020). Thus, in both the HD-H₂ and D₂-H₂ cases, we included centrifugal distortion of the PES in the scattering calculations. This leads to a large number of coupling terms (22 960 for v = 0 state in HD, 16 359 in v = 0 of D₂ and 17 157 for v = 2), which is an order of magnitude more than in the case of HD-He (1 029). Similar to the HD-He case, this effect is crucial for achieving a sub-percent agreement with the cavity-enhanced spectra.

Appendix C Convergence of quantum scattering calculations and uncertainty budget for line-shape parameters

In this appendix we provide a detailed analysis of the convergence of generalized spectroscopic cross-sections with respect to specific parameters and how these parameters influence the uncertainty of the line-shape parameters. They include the range of propagation, the propagator step, the number of terms in the expansion of the PES (B.1), the size of the rotational basis set, and the number of partial waves. The uncertainties were estimated by calculating six line-shape parameters for the R(0) line at temperatures ranging from 20 to 1000 K. These values were then compared with those obtained from the generalized spectroscopic cross-sections calculated with the reference values of each parameter, which are significantly larger than the ones used in the final calculations. The stated uncertainties represent the maximum error observed within this temperature range. We assume that a similar level of accuracy is maintained for all other transitions considered in this paper. The results are summarized in Table C.1.

The range of propagation is defined by the starting and ending points (R_min and R_max). The smallest R value that can be reliably calculated using quantum chemistry methods was employed as R_min, which in this case (taking into account the coordinate transformation from the H₂-H₂ to HD-H₂ system) was 3 a₀. On the other hand, R_max should be large enough to apply boundary conditions to the scattering equations (i.e., in the range of R where the PES becomes negligible compared to the centrifugal barrier). We tested the sensitivity of our results to R_max by performing calculations with R_max = 50 a₀, 75 a₀, 100 a₀, 150 a₀, and 200 a₀. After assessing the tradeoff between computational cost and accuracy, we chose R_max = 100 a₀ for the final calculations. Table C.1 provides uncertainties for six line-shape parameters impacted by this choice, estimated with respect to the reference calculations with R_max = 200 a₀.

The step size of the propagation directly affects the precision and the computational cost of the quantum scattering calculations. We performed tests with a varying number of steps per half-de Broglie wavelength (N_steps), including 10, 20, 30, 50, 100, 200, and 500. Based on these tests, we chose a step size of 50 for E_kin > 3 cm⁻¹, 100 for E_kin ∈ (1.5, 3〉 cm⁻¹, and 200 for E_kin ≤ 1.5 cm⁻¹. Uncertainties introduced by the choice of the number of steps, estimated with respect to the reference calculations with 500 steps per half-de Broglie wavelength, are gathered in Table C.1.

We tested the convergence of the results with respect to the number of terms in the PES expansion (Eq. (B.1)), comparing line-shape parameters for the R(0) line obtained from cross-sections calculated using a truncated expansion of the PES (with terms up to the I_{4 4 8}(θ₁, θ₂, ϕ) bi-spherical harmonic) and an expanded set of expansion coefficients describing higher anisotropies of the system (up to the I_{6 6 12}(θ₁, θ₂, ϕ) term). The results are gathered in Table C.1.

The number of partial waves (or equivalently, blocks with given total angular momentum J) necessary to converge the scattering equations was determined based on a criterion of stability in the calculated cross-sections. We solved the coupled equations for an increasing number of J-blocks until four consecutive J-blocks contributed to the largest elastic and inelastic state-to-state cross-sections by less than 10⁻⁴ Ų. The convergence criterion ensured that the estimated error introduced by the number of partial waves was smaller than the smallest uncertainty attributable to the other parameters in our study. This implies that the uncertainty in the number of partial waves did not significantly contribute to the overall uncertainty in our results; thus, we do not consider this factor in Table C.1.

The size of the rotational basis set is a critical factor in quantum scattering calculations, and it was chosen with great care to ensure a consistent level of accuracy across different rotational states. For each calculation, we checked that the basis set included all energetically accessible (open) levels of the colliding pair, as well as a certain number of asymptotically energetically inaccessible (closed) levels. We gradually increased the size of the basis set until the calculated cross-sections did not show appreciable differences, identifying a fully converged basis set. We then determined the smallest basis set that ensured convergence to better than 1% with respect to the fully converged basis. This was done for each initial state of the HD-H₂ system in a way that the estimated error for all transitions (R(j_HD), j_HD = 0, 1, 2), including those involving rotationally excited states, remained within the specified limit.

The tests were conducted separately for collisions with para-H₂ (which involves only even rotational quantum numbers) and ortho-H₂ (which involves only odd rotational quantum numbers). In the case of para-H₂, all rotational levels of HD and H₂ with j ≤ j^max = 4 were consistently included in the calculations. For specific calculations with H₂ initially in the $j_{{\text{H}}_2}=4$ state, or HD initially in the j_HD = 3 state, the basis set was expanded to incorporate all rotational levels of HD and H₂ with j ≤ 6. For ortho-H₂, the basis set consistently included all rotational levels of HD and H₂ with j ≤ j^max = 5. We extended the basis set to cover $j_{\text{HD}}^{\max }=6$ and $j_{{\text{H}}_2}^{\max }=5$ for cases where HD and H₂ were initially in the (j_HD, $j_{H_2}$) = (0, 3), (1, 3), or (2, 1) states. The largest basis set, with $j_{\text{HD}}^{\max }=j_{{\text{H}}_2}=7$, was employed in all calculations involving HD and H₂ in (j_HD, $j_{{\text{H}}_2}$) = (0, 5), (1, 5), (2, 3), (2, 5), (3, 1), (3, 3), and (3, 5) states.

Assuming the uncertainties associated with each parameter as independent, we estimated the maximum total uncertainty of each line-shape parameter using the root-sum-square method. The results are gathered in the last line of Table C.1.

Appendix D Modified Hartmann-Tran profile

This appendix describes the mHT profile, which we used to simulate the spectra. We considered the mHT profile to be the best compromise between the accurate but computationally demanding SDBB profile and the simple Voigt profile. The mHT profile can be expressed as a quotient of two quadratic speed-dependent Voigt (qSDV) profiles,

$${\tilde{I}}_{\text{mHT}}(f)=\frac{{\tilde{I}}_{\text{qSDV}}^{*}(f)}{1-\left(v_{\text{opt}}^r+i{v}_{\text{opt}}^i\right)\pi {\tilde{I}}_{\text{qSDV}}^{*}(f)},$$(D.1)

which are directly linked to the spectral line-shape parameters from Eqs (1)–(3),

$${\tilde{I}}_{\text{qSDV}}^{*}(f)=\frac{1}{\pi }{\displaystyle \int d^3}\upsilon f_m(\upsilon )\frac{1}{{\text{Γ}}_0+i{\text{Δ}}_0+\left({\text{Γ}}_2+i{\text{Δ}}_2\right)\left({\upsilon }^2/{\upsilon }_m^2-3/2\right)+v_{\text{opt}}^r+i{v}_{\text{opt}}^i-i\left(f-f_0-f_D{\upsilon }_z/{\upsilon }_m\right)}.$$(D.2)

The line-shape profile uses the parameters in the pressure-dependent form:

$$\begin{array}{l}{\text{Γ}}_0={\gamma }_0\cdot p,{\text{Δ}}_0={\delta }_0\cdot p\hfill \\ {\text{Γ}}_2={\gamma }_2\cdot p,{\text{Δ}}_2={\delta }_2\cdot p\hfill \\ v_{\text{opt}}^r={\tilde{v}}_{\text{opt}}^r\cdot p,v_{\text{opt}}^i={\tilde{v}}_{\text{opt}}^i\cdot p.\hfill \end{array}$$(D.3)

The f_m(υ) is the Maxwell-Boltzmann distribution of the active molecule velocity, υ_m is its most probable speed, and υ_z is one of the three Cartesian components of the υ vector. The f, f₀, and f_D are the frequency of light, the central frequency of the transition, and the Doppler frequency, respectively.

The hard-collision model of the velocity-changing collisions, which is used in the mHT profile, suffices to describe the velocity-changing line-shape effects (such as the Dicke narrowing) in the majority of the molecular species. However, in the cases with a significant Dicke narrowing, such as molecular hydrogen transitions, the hard-collision model does not reproduce the line shapes at the required accuracy level. To overcome this problem, a simple analytical correction (the β correction function) was introduced (Wcisło et al. 2016; Konefał et al. 2020), which mimics the behavior of the billiard ball model and, hence, considerably improves the accuracy of the mHT profile for hydrogen, at negligible numerical cost. The correction is made by replacing the $v_{\text{opt }}^r$ with ${\beta }_{\alpha }(\chi )v_{\text{opt}}^r$, where a is the perturber-to-absorber mass ratio and $\chi =v_{\text{opt}}^r/{\text{Γ}}_D$ (where Γ_D is the Doppler width; see Konefał et al. (2020) for details). It should be emphasized that the β correction does not require any additional transition-specific parameters (it depends only on the perturber-to-absorber mass ratio α). The β correction was applied every time the mHT profile was used in this work.

Appendix E DPL representation of the temperature dependences

In this appendix, we discuss the details of the DPL representation of the temperature dependences of the spectral line-shape parameters. The DPL function is used to convert the exact temperature dependence of the line-shape parameters into a simple, analytical expression, suitable for storing in spectroscopic databases (Stolarczyk et al. (2020)). This conversion is done by fitting the DPL function to the actual ab initio temperature dependence data.

Within this work, we calculated the ab initio values at temperatures ranging from 50 to 1000 K. Due to its relevance to the atmospheres of giant planets, we chose to prioritize the 50-200 K temperature range and use the data only from this range to generate the DPL coefficients (Stolarczyk et al. 2020). Projection of the ab initio data on the DPL representation is performed by fitting the DPL function in the selected temperature range. For the most faithful reconstruction of the temperature dependence, we performed several different fitting procedures (i.e., Newton, Quasi-Newton, Levenberg-Marquardt, global optimization, and gradient methods) and selected the one that gives the best result (i.e., the lowest relative root mean square error; see the next paragraph). This selection was done separately for each of the line-shape parameters and molecular transitions. Mathematically, the two terms of the DPL functions are identical. Thus, to avoid swapping them, we followed the convention that the first base coefficients should always be greater than second base coefficients. Furthermore, to reduce the number of significant digits, if the two base coefficients have opposite signs, we required that their absolute values differ at least by 1%o.

The DPL coefficients listed in Table 1 can be used to retrieve the temperature dependence of the line-shape parameters through Eq. (4). The efficiency of the DPL representation is depicted in Fig. E.1. The red curves show the results of the ab initio calculations, while the black curves (covered by the red ones in some cases) are the values reconstructed from the DPL coefficients from Table 1. The corresponding residuals are presented under each of the plots. Table E.1 quantifies the accuracy of the DPL representation by presenting the values of the relative root mean square error (rRMSE) of the differences between the ab initio data and the DPL fit. The values of the rRMSE are normalized with respect to the value of the corresponding line-shape parameter at 296 K. Even though the rRMSE values for some parameters are on the level of several percent, the overall error of the shape of the line is much smaller because the two parameters that impose the highest impact on the line width, γ₀ and ${\tilde{v}}_{\text{opt}}^r$, are reproduced with high accuracy. A detailed discussion of the propagation of errors from the parameters to the final shape of the line is provided in Słowiñski et al. (2022).

Fig. E.1 Temperature dependences of the six collisional line-shape parameters, γ0, δ0, γ2, δ2, ${\tilde{v}}_{\text{opt}}^r$, and ${\tilde{v}}_{\text{opt }}^i$, of the R(0) 0-0, R(1) 0-0, and R(2) 0-0 lines of HD perturbed by H2. The red and black curves are the ab initio results and DPL approximations, respectively. The small panels show the residuals from the DPL fits. The vertical axes for all the panels (including residuals) are in [10−3 cm−1 atm−1.

References

All Tables

All Figures

	Fig. 1 Ab initio temperature dependences of the collisional line-shape parameters (in units of 10−3 cm−1 atm−1 ) of the first three electric dipole lines of HD perturbed by H2 (red curves) and He (black curves).
In the text

Fig. 2 Impact of the beyond-Voigt effects and bath mixture composition on collision-perturbed spectra of HD at conditions relevant for giant planet atmospheres. Panel a: relative error of the Voigt-profile approximation as a function of pressure and temperature, shown as the relative difference between the Voigt and mHT profiles at profile maximum. The panels show, from top to bottom, the R(0), R(1), and R(2) lines. Panels b–d: simulations of the HD spectra (blue lines) at conditions relevant for the Neptune atmosphere (the perturbing bath is 80% H2 and 20% He). The spectra are generated with the mHT profile using HAPI based on the DPL temperature parametrization. As a reference, we show the same lines for the cases of pure H2 and pure He perturbers (see the red and black lines, respectively). The blue shadows show the same simulations as the blue lines but generated with the simple Voigt profile. Panels b–d correspond to points b–d shown in the temperature-pressure maps in panel a (the three selected points lie on the Neptune temperature-pressure line). The three cases illustrate three different line-shape regimes. The first one, (b), is the low-pressure case in which the lines are broadened mainly by the Doppler effect, and the pressure-induced collisional effects do not dominate the line shapes. The intermediate-pressure case, panel c, illustrates the extreme non-Voigt regime (the differences between the blue curves and blue shadows reach almost a factor of 2; see also the green ridge in the maps in the bottom panel). The third case, panel d, illustrates the high-pressure regime at which the HD lines are well described by a simple Voigt profile (the blue shadows almost overlap with the blue lines), but setting a proper composition of the perturber gas components plays an important role.

In the text

	Fig. 3 Dependence of pressure broadening, γ0 (solid lines), and shift, δ0 (dashed lines), on the spin isomer concentration ratio (x = np/no) for the H2-perturbed R(0) line in HD at different temperatures. The vertical gray lines correspond to the value of x for normal H2, x = 1/3, and xeq(T), as determined by the Boltzmann distribution at T = 50, 195, and 296 K.
In the text

Fig. 4 Comparison of the experimental and theoretical values of the pressure broadening and shift parameters, γ0 and δ0 (in units of 10−3 cm−1 atm−1 ). Black curves correspond to the ab initio calculations performed in this work, while green and red points report the experimental measurements from Ulivi et al. (1989) and Lu et al. (1993), respectively. The olive curves are the single-power-law (for γ0) and linear (for δ0) temperature dependences retrieved from the measurements of Sung et al. (2022).

In the text

Fig. 5 Direct validation of the ab initio quantum-scattering calculations on the accurate experimental spectra of the S(2) 2−0 line of D2 perturbed by collisions with H2 molecules (see the main text for details). The black dots are the experimental spectra, and the red lines are the ab initio profiles. Below each profile, we show the absolute residuals of two models: the SDBB profile and the mHT profile. To quantify how well theory agrees with experiments, we report the relative (with respect to the profile peak value) root mean square errors (rRMSEs) of the experiment-theory differences calculated within the ±FWHM range around the line center (see the percentages below the residuals). The mean rRMSEs are also summarized for each of the models (see the numbers on the right side of the figure).

In the text

Fig. E.1 Temperature dependences of the six collisional line-shape parameters, γ0, δ0, γ2, δ2, ${\tilde{v}}_{\text{opt}}^r$, and ${\tilde{v}}_{\text{opt }}^i$, of the R(0) 0-0, R(1) 0-0, and R(2) 0-0 lines of HD perturbed by H2. The red and black curves are the ab initio results and DPL approximations, respectively. The small panels show the residuals from the DPL fits. The vertical axes for all the panels (including residuals) are in [10−3 cm−1 atm−1.

In the text