ExoMol line lists -- L: High-resolution line lists of H$_3^+$, H$_2$D$^+$, D$_2$H$^+$ and D$_3^+$

New MiZo line lists are presented for the D$_2$H$^+$ and D$_3^+$ isotopologues of H$_3^+$. These line lists plus the existing H$_3^+$ MiZATeP and the Sochi H$_2$D$^+$ line lists are updated using empirical energy levels generated using the MARVEL procedure for H$_3^+$, H$_2$D$^+$ and D$_2$H$^+$, and effective Hamiltonian energies for D$_3^+$ for which there is significantly less laboratory data available. These updates allow accurate frequencies for far infrared lines for these species to be predicted. Assignments of the energy levels of H$_3^+$ and D$_3^+$ are extended using a combination of high accuracy variational calculations and analysis of transition intensities. All line lists are made available via www.exomol.com.


INTRODUCTION
H + 3 is known to form rapidly in H2 gas following an ionisation event via the strongly exothermic reaction H2 + H + 2 −→ H + 3 + H (1) which occurs at essentially every collision. As H2 is common in a variety of astronomical bodies, H + 3 is often the dominant molecular ion. The first 30 years of H + 3 astronomy has been comprehensively reviewed by Miller et al. (2020). H + 3 formation is stimulated by cosmic ray ionisation of the interstellar medium and by collisions with fast electrons and other charged particles in planetary ionospheres. H + 3 is also believed to form in the ionosphere of planets through the ionisation of H2 by extreme ultraviolet radiation (Chadney et al. 2016). In gas giant ionospheres, H + 3 acts as a coolant through efficient infrared (IR) emissions (Miller et al. 2010); indeed it is thought that H + 3 emissions are key to determining the stability limits in hot Jupiter exoplanets (Koskinen et al. 2007).
The infrared spectrum of H + 3 has been extensively observed in giant planets in our solar system such as Jupiter (Drossart et al. 1989;Ballester et al. 1994;Miller et al. 1997; Moore et al. 2017), Saturn Stallard et al. 2008a,b), Uranus (Trafton et al. 1993;Lam et al. 1997;Trafton et al. 1999;Melin et al. 2019) and is believed to be ⋆ The corresponding author: j.tennyson@ucl.ac.uk present in Neptune (Melin et al. 2011(Melin et al. , 2018 although it is yet to be detected there. Its presence can be used as an effective temperature probe here and in other astrophysical settings (Gibbs & Fitzgerald 2022). H + 3 is similarly expected to be of importance in extrasolar giant planets (Chadney et al. 2016;Khodachenko et al. 2015), such as hot-Jupiters (Lenz et al. 2016), and an even more prominent feature in the aurorae of brown dwarfs (Gibbs & Fitzgerald 2022); however, it has so far defied observation on these objects.
H + 3 has also been observed in the interstellar medium (ISM) via absorption in the infrared light of a background star (Oka 2006) where it forms through cosmic ray ionisation (Geballe & Oka 1996). Hence it is also used in this setting to trace the cosmic ray ionisation rate (Indriolo & McCall 2012;Harju et al. 2017) and similarly primarily relies on IR emissions. These IR bands lie well within the wavelength range of the instruments onboard the JWST.
In regions where the precursors to H + 3 exist in deuterated forms, namely HD and D2, equivalent reactions occur to that described in Eq. (1) resulting in the formation of the deuterated isotopologues H2D + , D2H + and D + 3 (Merkt et al. 2022). At low temperatures fractionation drives the preferential formation of isotopically substituted H + 3 (Hewitt et al. 2005); indeed, models by Walmsley et al. (2004) suggest that in certain very cold regions D + 3 may be the dominant isotopologue of H + 3 ! Spectra of H2D + (Stark et al. 1999;Caselli et al. 2003) and D2H + (Vastel et al. 2004) have been observed in interstellar space through pure rotational transitions which lie in the far infrared / THZ region. However, D + 3 remains undetected, at least in part because its higher symmetry means that, like H + 3 , its pure rotational spectrum is very weak. Elsewhere, electrons provided by H + 3 have been shown to play an important role in the atmospheres of cool white dwarfs (Bergeron et al. 1997).
Line lists for H + 3 (Kao et al. 1991;Neale et al. 1996) and the associated partition function (Neale & Tennyson 1995; have played a key role in the astronomical study of this important molecular ion. In this work we update the MiZATeP H + 3 line list of Mizus et al. (2017) and the older ST1 H2D + line list of Sochi & Tennyson (2010). We do this using updated versions of the marvel (measured active rotation-vibration energy levels) studies due to Furtenbacher et al. (2013a) and Furtenbacher et al. (2013b). We present new line lists for D2H + and D + 3 , for which we also use empirical energy levels to improve the accuracy of the transition frequencies between key levels. These line lists are produced as part of the ExoMol project (Tennyson & Yurchenko 2012).

METHOD
The triatomic discrete variable representation (DVR) nuclear motion code DVR3D ) was used previously and here to compute initial energy levels for H + 3 and its deuterated isotopologues as well as Einstein A-coefficients for each transition. This code, which is based on the use of an exact nuclear motion kinetic energy operator, has been shown to be capable of giving highly accurate results for the H + 3 system (Polyansky & Tennyson 1999;Pavanello et al. 2012a). It is important to note that in the absence of any absolute measurements of transition intensities for the H + 3 system, all models rely on computed values which are thought to be accurate (Farnik et al. 2002;Petrignani et al. 2014). It should be noted that DVR3D only provides assignments for the rigorous quantum numbers: J, rotationless parity e/f and the interchange symmetry for two identical atoms. This means that most states in the existing version of the MiZATeP H + 3 (Mizus et al. 2017) and ST1 H2D + (Sochi & Tennyson 2010) line lists do not have full ro-vibrational labels. We partially address this issue below.
Marvel (Furtenbacher et al. 2007;Furtenbacher & Császár 2012;Tóbiás et al. 2018) takes assigned, high resolution spectra and uses them to construct empirical energy levels with spectroscopic accuracy and specified uncertainties. Use of the energy levels can greatly improve the accuracy with which a line list can predict transition frequencies, see Al-Derzi et al. (2021) for a recent example. For H + 3 , H2D + and D2H + marvel spectroscopic networks were then constructed using available laboratory spectra in order to obtain empirical energy levels for each species. These empirical energy levels were then used to improve our line lists, correcting for obs.-calc. shifts in levels where empirical energy levels are available for comparison. Such a refinement allows for a subset of the energies to be provided with very high accuracy, as has been demonstrated in similar projects (Bowesman et al. 2022). This allows for high-accuracy transition frequency predictions to be made (Al-Derzi et al. 2021), making the final line list well suited for high-resolution studies (Bowesman et al. 2021;Owens et al. 2022).
The marvel process determines an uncertainty for each energy level based on the uncertainties of the input transition that define that level.

Spectroscopic Networks
Transition data for a molecule can be aggregated to construct a spectroscopic network, where the transition frequencies represent the edges of the networks and the energy levels the nodes (Furtenbacher et al. 2007;Császár & Furtenbacher 2011). The marvel procedure (Furtenbacher & Császár 2012) achieves this by inverting transition matrices, which yields a set of empirical energy levels with individual uncertainties. Spectroscopic networks have been constructed in the past for the molecules H + 3 (Furtenbacher et al. 2013a), H2D + and D2H + (Furtenbacher et al. 2013b) but new transition data has since been published. New transition data allows us to expand the level coverage of the networks and in the case of new data that remeasures existing transitions to improve the accuracy to which term energies are known. The recent high-resolution experiments have provided new THz transition data with uncertainties on the order of MHz or kHz, well within the part-per-billion regime. When all of the transitions that determine a level's term energy are consistent within their experimental uncertainties, the uncertainty on the final term energy will generally be on the same order, but not less than, the smallest uncertainty of the transitions that define the level. Hence, with the addition of these highaccuracy transition measurements we are able to significantly improve the accuracy of our final energies, in some cases by a few orders of magnitude.
The marvel procedure requires all transitions within a network to be identified by the same set of quantum numbers. For this purpose, the deuterated isotopologues can be divided into two groups: H + 3 and D + 3 are symmetric-tops belonging to the D 3h (M) molecular symmetry group and H2D + and D2H + are asymmetric-tops belonging to the C2v(M) group. These in turn dictate the set of good quantum numbers used to define the levels of the networks.
As symmetric tops, the molecules H + 3 and D + 3 are defined by two primary vibrational modes, symmetric stretching (ν1) and bending (ν2). The bending mode is degenerate however, and as such these species are also described by the vibrational angular momentum quantum number l2 = ν2, ν2 − 2, ..., −ν2 + 2, −ν2. The energy levels of these species are identified instead by L2 = |l2| in this work. The rigorous rotational angular momentum quantum number J is used to define the rotational levels of these molecules. The projection of J along the molecular symmetry axis, k can be used to determine the parity of the energy levels, such that the total parity is the sign of (−1) k (Furtenbacher et al. 2013a). The quantum number k does not offer a complete description of the system however, due to Coriolis coupling between it and l2. As such, the value G = |g| is used, where g = k − l2 (Hougen 1962). This quantum number is important because, as members of the D 3h (M) symmetry group, the A ′ 1 , A ′′ 1 , A ′ 2 , A ′′ 2 rovibronic symmetries only exist for these molecules when G = 3n, where n is an integer. In H + 3 the A ′ 1 and A ′′ 1 symmetries do not exist however, as they are determined to have 0 nuclear spin statistical weight . This difference arises from the nuclear spin I of constituent atoms, which are I = 1 2 for hydrogen and I = 1 for deuterium. Consequently H + 3 consists of two nuclear spin isomers while D + 3 has three ; these are shown in Table 1.
Levels with equivalent vibrational, J and G assignments are differentiated by the (u|l|m) U -notation of Watson (1994). u and l are used to identify the upper and lower energy levels of the same assignment, differing in their value of K = |k|, and m is used when such a distinction is irrelevant as only one value of K can exist that produces the same G. As such, the energy levels of the species H + 3 and D + 3 are identified by the quantum number set (ν1, ν2, L2, J, G, U, K, Γrve).
As asymmetric-tops, H2D + and D2H + are defined by the symmetric stretching, bending and anti-symmetric stretching vibrational quantum numbers ν1, ν2 and ν3. The quantum number J and its projections onto the C2 axis and the axis perpendicular to the C2 axis in the plane of the molecule, Ka and Kc, are the standard ones for labelling the rotational states of an asymmetric top.
For both asymmetric-tops, the total parity is the sign of (−1) Kc . A similar expression is used to identify the spin isomers of each asymmetric-top, such that they are labelled ortho when it is positive and para when negative. For H2D + this is expression is (−1) ν 3 +Ka and for D2H + it is (−1) ν 3 +Ka+Kc . The ortho and para nuclear spin isomers identify the rovibronic symmetry group of the molecule, A1, A2, B1 or B2, as shown in Table 1. Hence, the energy levels of the species H2D + and D2H + are identified by the quantum number set (ν1, ν2, ν3, J, Ka, Kc, Γrve).
Transitions between nuclear spin isomers are forbidden, meaning networks constructed from observational transitions will be split into separate components for each isomer. Given marvel determines the energies of each level relative to the lowest energy level in the network, which is defined as 0, this presents a problem as any levels of a different nuclear spin isomer to that of the zero-energy level will not have their energies defined. This is avoided by the introduction of "magic" numbers; forbidden transitions are added to the networks to connect the lowest levels of each nuclear spin isomer to each other using calculated energies. This enables the determination of all energies for connected levels within the networks, relative to the zero-energy level. For D + 3 this treatment was not needed as the energies were determined through the use of effective Hamiltonian calculations. which combined transition data from 26 sources into a single network. Since then, five new sets of high-accuracy spectra have been reported (Berg et al. 2012;Hodges et al. 2013;Perry et al. 2015;Jusko et al. 2016;Guan et al. 2018; and are detailed below: 12BeWoPe (Berg et al. 2012): 3 R-branch transitions in the visible region are provided that had not been included from other sources.
13HoPeJeSi (Hodges et al. 2013): 10 R-branch transitions in the ν2 band with MHz and sub-MHz accuracy are reported. All of these transitions had been previously reported in other sources (Oka 1981;McKellar & Watson 1998) but are presented here at higher accuracy.
15PeHoMaKo (Perry et al. 2015): A further 10 R-branch transitions in the ν2 band are reported with MHz and sub-MHz accuracy. All of the transitions from this source had been observed previously (Oka 1981;Wu et al. 2013).
16JuKoScAs (Jusko et al. 2016): 5 measurements of low-J R-branch transitions in the ν2 band are provided, all with sub-MHz accuracy. These transitions had also been observed by 13HoPeJeSi (Hodges et al. 2013).
18GuChLiPe (Guan et al. 2018): Transition frequencies for 12 R-branch and 4 Q-branch transitions in the ν2 band are presented. The majority are reported with sub-MHz accuracy. All of the transitions in this source had also been observed in other experiments (Lindsay & McCall 2001;Hodges et al. 2013;Perry et al. 2015;Jusko et al. 2016;Perry et al. 2016).
19MaMc : The largest of the new H + 3 data sources, providing 36 measurements of transitions in the ν2 ← 0 band, 15 transitions in the 2ν 2 2 ← ν2 band and 7 transitions in the 2ν 2 2 ← 0 band. The data comprises a selection of P-, Q-and R-branch transitions with uncertainties of 4 MHz. All of the reported transitions had been observed in prior studies Majewski et al. 1987;Bawendi et al. 1990;Xu et al. 1990;Uy et al. 1994;McKellar & Watson 1998;Guan et al. 2018;Markus et al. 2018).
In total there were 102 new transitions to be added to the existing set of 1 610 transitions and a further 7 transitions that had been missed in the original marvel compilation. The inclusion of the new transitions allowed us to validate 19 previously invalidated transitions. 7 transitions had their assignment corrected based on updated assignments provided in subsequent papers. A further 9 transitions were reassigned based on apparent mixing identified through comparison with the empirical marvel energies. 2 transition assignments were updated based on suggestions provided by Jaquet (2022). Sarka & Poirier (2022) provided a series of assignments for H + 3 which were used to correct the assignments of 41 transitions and to assign 22 previously unassigned transitions. A further three unassigned transitions from the literature were assigned based on comparison with the marvel energy levels. The new transition assignments are detailed in Table 2. The transitions updated in this way all involve levels with relatively high energies, in a region where accurate determination of the vibrational state can be challenging. In total, 59 transitions from the original network were reassigned and an additional 8 had their assignment corrected due to digitisation errors from the original source.
The final set of 1 719 transitions were passed through marvel for validation, with 1 656 being validated; these transitions are summarised in Table 3. The network is divided into 39 components, the largest of which comprises 1 580 tran- Table 2. The new transition assignments for transitions from the literature that previously had no assignment or were missing an upper state assignment. Assignments to transitions from  were determined through comparison between the empirical energy levels determined by marvel, while transitions from Morong et al. (2009) were assigned based on assignments produced using the code ScalIT (Sarka & Poirier 2022). Morong et al. (2009) sitions which determine 703 empirical energy levels. Of the remaining components: two consist of 6 transitions; one of 5 transitions; two of 4 transitions; three of 3 transitions; 12 of 2 transitions and the remaining 18 of single transitions. As these minor components are not connected to the primary component the energies of the levels within them are not determined relative to the zero-energy level and are hence not considered further.

H2D + network
A marvel network for H2D + was originally produced by Furtenbacher et al. (2013b) and contained transition data from 13 sources. Since then, two new sets of transition measurements have been published by Jusko et al. (2016) and Jusko et al. (2017).
16JuKoScAs (Jusko et al. 2016): Transitions frequencies for 11 ν1 band transitions are reported with MHz and sub-MHz accuracy. This source provides significantly higher accuracy measurements of transitions previously observed by Amano (1985).
These two new sources and a "magic" forbidden transition to connect the nuclear spin isomers, derived from effective Hamiltonian calculations using molecular constants from Amano (2006), brought the total number of transitions for H2D + to 210; these are summarised in Table 4. 208 transitions were validated successfully by marvel, yielding a final network comprising 7 components, the largest of which contained 200 transitions that determine 109 unique energy levels. For the remaining 6 components, one contained 3 tran-sitions while the rest are single transition components; these disconnected components are not considered further here.

D2H + network
The D2H + marvel network was originally published by Furtenbacher et al. (2013b) and was constructed using transition data from 9 sources. Four new sets of transition data have subsequently been published (Jusko et al. 2016(Jusko et al. , 2017Yu. et al. 2017; and have now been added to the existing network. 16JuKoScAs (Jusko et al. 2016): This source provides 10 transitions frequencies observed with sub-MHz accuracy. 8 of these transitions had been previously observed by Lubic & Amano (1984).
17YuPeAmMa (Yu. et al. 2017): This source reports 5 ground state pure-rotation transitions with sub-MHz accuracy. 4 of these transitions had not been reported in other sources, with the other also being observed by Jusko et al. (2017).
19MaKoMc ): Transition frequencies for 37 ν1 band transitions are provided with MHz accuracy. 10 of the transitions in this source are new and had not been reported elsewhere, while the rest had been observed by Lubic & Amano (1984) and Jusko et al. (2016).
With the addition of 55 new transition measurements, the new network contains 210 transitions. Two existing transitions that had digitisation errors in their transitions frequencies were updated to the correct values. A "magic" forbidden transition is included to connect the otherwise distinct ortho and para components of the network, using an frequency cal- Table 3. The components of the updated H + 3 marvel network. The network is broken down by the transition data sources and the vibrational bands contained within them, given in the form For each band, the transition energy and total angular momentum J ranges are provided, along with the mean and maximum uncertainties for these transitions. The total number of transitions validated and accessed by marvel (V/A) are also provided. (Oka 1980    The network is broken down by the transition data sources and the vibrational bands contained within them, given in the form For each band, the transition energy and total angular momentum J ranges are provided, along with the mean and maximum uncertainties for these transitions. The total number of transitions validated and accessed by marvel (V/A) are also provided.

Effective Hamiltonians
There have been significantly fewer observations of D + 3 spectra than of the three other isotopologues considered here. As such, there was insufficient data to build a well-connected network. In lieu of this, we used effective Hamiltonian constants from Watson et al. (1987) and Amano et al. (1994) to calculate energies for the states in the range J = 0 − 15, up to a maximum energy of 2676.387 cm −1 . These calculations were performed using the program pgopher (Western 2017) which also provides full state assignments for the levels it computes. Hence full quantum number assignments were determined via this method for 282 levels within the bands for which constants were available: 188 levels in the 000 band (ν1, ν2, L2); 2 levels in the 010 band; 105 in the 011 band; 22 in the 100 band. The majority of the levels assigned through this method have J < 10. Further assignments were done manually for states at higher energies with the aid of the Table 5. The components of the updated D 2 H + marvel network. The network is broken down by the transition data sources and the vibrational bands contained within them, given in the form For each band, the transition energy and total angular momentum J ranges are provided, along with the mean and maximum uncertainties for these transitions. The total number of transitions validated and accessed by marvel (V/A) are also provided.

Vib.
J range V/A Energy range (cm −1 ) Uncertainty Mean/Max (cm −1 ) 84LuAm (Lubic & Amano 1984) 100 -000 0 -6 34/34 2638 -2990 0.005/0.005 86FoMcWa (Foster et al. 1986b Amano et al. (1994) and the assigned hot and overtone bands published by Alijah et al. (1995). A further 1045 states were assigned this way and a breakdown of the assigned bands is given in Table 6. These assignments were added to the calculated states file, given that DVR3D only provides assignments for the rigorous quantum numbers: J, rotationless parity e/f and interchange of two of the D atoms.

LINE LIST CALCULATIONS
3.1 Updated H + 3 and H2D + While the main line lists were computed using DVR3D, we  symmetry group. Due to the very high convergence accuracy, the full G12 PI group labels (Γrve) were easily assigned unambiguously using the Γ(G12/D 3h ) ↓ G4 correlation table.
Next, calculations were repeated with the code GENIUSH with slightly lower accuracy, but still sufficient to match the energy levels with the ScalIT ones. After the vibrational states were labelled (ν1, ν2, L2), vibrational parent labels were semi-automatically assigned to rovibrational states using the rigid rotor decomposition scheme (RRD) implemented in GENIUSH (Mátyus et al. 2010). The RRD overlaps of GENIUSH also help to assign the rotational quantum numbers (J, G, U, K), as the symmetric top rigid rotor functions are labelled by K. Using this method we were able to label by the quantum number set (ν1, ν2, L2, J, G, U, K, Γrve) an additional 1 525 states in the MiZATeP H + 3 line list. The energies of the existing H + 3 and H2D + line lists were updated with empirical energies where they were known. For levels with matching assignments in both the states files and corresponding marvel network, the term energies and their uncertainties were set to the values determined by the marvel procedure. For the levels that did not exist within the molecules' marvel network, the calculated term energies were retained and estimates for their uncertainties were calculated as follows ∆Ẽ = 0.1, ifẼ < 2000 ⌊Ẽ/2000⌋/10, otherwise. (2)

New line lists
New line lists were computed for D2H + and D + 3 molecular ions. Both calculations used the highly accurate global ab initio PES, together with adiabatic and relativistic correction surfaces, computed by Pavanello et al. (2012a,b) which were used for the MiZATeP line list. To calculate transitions intensities the high-accuracy DMS obtained for the H + 3 system was used (Petrignani et al. 2014). This DMS was obtained by fitting 7 parameters to a polynomial form written in terms of effective charges, see Röhse et al. (1994). This DMS was centred on the centre-of-mass ensuring the correct treatment of the centre-of-charge -centre-of-mass displacement which leads to D2H + (and H2D + ) having a permanent dipole moment.
The DVR3D program suite ) was used to compute the final line lists. As part of this project a new module was added to this suite which converts DVR3D results to ExoMol format (Tennyson et al. 2013) comprising a states and trans file, and allowing spectra to be easily generated using ExoCross (Yurchenko et al. 2018). The module reads the energy levels and Einstein coefficients from DVR3D, requiring as input from the user the molecular symmetry, Cs or C2v, and the nuclear statistic weight for each symmetry. The program, called LINELIST, is available as part of the DVR3D program suite from the ExoMol GitHub pages.
As with H + 3 and H2D + , the calculated term energies of D2H + were updated with empirical values from the new marvel network where available. Likewise, the unchanged calculated energies of the new D2H + and D + 3 line lists had uncertainties estimated using Eq. (2).

D2H + nuclear motion calculations
DVR3D nuclear motion calculations were performed in the following way: with 31, 31, and 50 grid points for two radial and an angular scattering coordinates. The calculations used Morse-like oscillators with parameters 3.1, 0.1, and 0.006 for r1 (D -D) radial coordinate, and spherical oscillators with parameters 0, 0, and 0.016 for the r2 (H -D2 scattering coordinate). The dimension of the final vibrational Hamiltonians were set to 5000. Calculations included all levels up to 15000 cm −1 for J ⩽ 25. 1500 vibrational basis functions calculated using DVR3DRJZ were passed to ROTLEV3 for the rotational step of the calculation. Following Polyansky & Tennyson (1999), nuclear masses mH = 1.007276 Da and mD = 2.013553 Da were used for rotational motion and effective masses intermediate between nuclear and atomic masses mH = 1.007537 Da and mD = 2.013814 Da were used for the vibrational motion. This formulation has been shown to allow for non-adiabatic effects in the calculation. The resulting MiZo line list is complete up to a temperature of 2000 K.

D + 3 nuclear motion calculations
DVR3D calculations for D + 3 were performed using the same size grids and Hamiltonians as those specified for D2H + above (31, 31, and 50 grid points for two radial and an angular scattering coordinates, correspondingly, and the final Hamiltonian dimensions equal to 5000). The calculations used Morselike oscillators with parameters 2.6, 0.1, and 0.006 for the r1 coordinate, and spherical oscillators with parameters 0, 0, and 0.016 for the r2 coordinate. Nuclear masses equal to mD = 2.013553 Da were used for all calculations.
On the basis of nuclear motion calculation results and the DMS by Petrignani et al. (2014) a linelist in the frequency region up to 15 000 cm −1 for a temperature of 800 K and using the corresponding partition function value 118.7 together with a list of corresponding energy levels for D + 3 was computed. This linelist contains transitions between energy states with J values 0-15 and energies 0-15 500 cm −1 . The new MiZo D + 3 line list is complete up to a temperature of 800 K.
Statistics for all of the line lists presented here are given in Table 7.

States files
New states files are provided for each of the four species considered here and are formatted using the standard outlined by Tennyson et al. (2013). H + 3 and D + 3 are identified using the same set of quantum numbers and hence their states files contain the same columns. Excerpts from the H + 3 and D + 3 states files are provided in Table 8. Similarly, the states files for H2D + and D2H + contain the same set of quantum number columns and an excerpts for them are shown in Table 9. All entries in the new states files are marked with a source tag, indicating the method that was used to determine the final term energy for that level. Several values for the source tag can occur in ExoMol line lists (Bowesman et al. 2021) but only three are in use here: "Ma" for marvelised energies, "EH" for energies from effective Hamiltonian calculations and "Ca" for energies calculated using DVR3D. For all levels marked "Ca" in all four states files an uncertainty estimate was calculated using Eq. (2).
All new and updated states files are ordered on increasing J, Γrve and energy. The existing H2D + states file was already formatted this way, while H + 3 was not and has been changed. Hence for H + 3 the state counting number assigned to each level has been reassigned and the references to these values in the corresponding trans file have accordingly been updated. Excerpts from the trans files for the new D2H + and D + 3 line lists can be seen in Table 10.
In general all molecular states should be characterised by three rigorous quantum numbers: J, parity and the overall symmetry, which also determines the nuclear spin state (ortho/para/meta). All calculations have J and parity rigorously determined; due to the full treatment by DVR3D of the C2v symmetry, the symmetry and nuclear spin statistics are also correctly determined for H2D + and D2H + . This is not the case for H + 3 and D + 3 however, as the program does not include a full treatment for a D 3h Hamiltonian. The symmetry was initially output for all molecular species considered here under C2v symmetry, but was subsequently mapped to the appropriate D 3h representation for H + 3 and D + 3 for states with  Table 8. Excerpts from the H + 3 and D + 3 states files, using the format defined by Tennyson et al. (2013). Note that the zero-energy level of H + 3 with A ′ 1 symmetry does not exist and hence has "nan" entries for its lifetime and isomer. Any row for which the ν 1 , ν 2 , L 2 , G, U or K assignment is not known will likewise have "nan" entries in these columns. Isomer: Nuclear spin isomer; ν 1 : Symmetric stretching quantum number; ν 2 : Bending quantum number; L 2 : Vibrational angular momentum quantum number; G: Absolute value of k − l 2 ; U : U -notation of Watson (1994); K: Rotational angular momentum quantum number; Source Tag: The method used to generate the term value; "Ma" for marvelised energies, "EH" for energies from effective Hamiltonian calculations and "Ca" for energies calculated using DVR3D.
quantum number assignments. Due to the nuclear spin statistical weight of the A ′ 1 , A ′′ 1 symmetries in H + 3 being 0, it was possible to map all of the C2v symmetries to the appropriate D 3h representation (Mizus et al. 2017). For D + 3 , where D 3h symmetry and quantum number assignment was not carried out, the entries retain the C2v symmetries output by DVR3D. These C2v symmetries are given as lower-case in the states file, to avoid confusion between the A1 and A2 irreducible representations under C2v symmetry (formatted as "a1" and "a2") and the A ′ 1 , A ′′ 1 , A ′ 2 and A ′′ 2 irreducible representations under D 3h symmetry.

Partition Functions
Calculating the partition functions of the new line lists requires the determination of the total degeneracy of each state Table 9. Excerpts from the H 2 D + and D 2 H + states files, using the format defined by Tennyson et al. (2013). Any row for which the ν 1 , ν 2 , ν 3 , Ka or Kc assignment is not known will have "nan" entries in these columns. considered, as show in Eq. (3): where gtot is the total degeneracy of the ith level, Ei its energy, kB is the Boltzmann constant and T is the temperature. The total degeneracy of a level depends directly on the nuclear spin statistical weight, gns: The nuclear spin statistical weights represent the number of nuclear spin functions that yield each symmetry. We follow the ExoMol and HITRAN convention of including the full nuclear spin degeneracy in our partition functions. In the case of H + 3 , where each atom consists of one proton with I = 1 2 , Fermi statistics apply which results in no allowed configurations corresponding to the A ′ 1 or A ′′ 1 representations. Accordingly, the gns values for these representations of H + 3 are 0, as shown in Table 1. For D + 3 however, each Deuterium atom has I = 1 and is hence a boson, meaning additional representations can occur that give rise to the "meta" nuclear spin isomer.
Our mappings of the C2v representations output by DVR3D to the full D 3h symmetry are complete for all states presented in the H + 3 line list, but not for the states of D + 3 . Hence, we must also consider the appropriate statistical weights to use when determining the degeneracy of the D + 3 states left with C2v symmetries. As the D 3h representations A ′ 1 , A ′′ 1 , A ′ 2 , A ′′ 2 , E ′ , E ′′ correspond to A1, A2, B2, B1, A1 ⊕ B2, A2 ⊕ B1 respectively under C2v symmetry, we can calculate population weighted average statistical weights for each C2v representation. These are calculated knowing that approximately two thirds of all D + 3 levels will be E ′ and E ′′ , due to the constraint that A ′ 1 , A ′′ 1 , A ′ 2 and A ′′ 2 levels only occur when G = 3n, where n is an integer; see also Berblinger et al. (1992). Accordingly, for the A1 and A2 representations in C2v symmetry, one third will be A ′ 1 or A ′′ 1 under D 3h symmetry while the remaining two thirds will be E ′ and E ′′ . The equivalent population ratio is true for B1 and B2 states and the corresponding A ′ 2 and A ′′ 2 representations. These population ratios are then weighted using the weights from Table 1,   yielding the gns values shown in Table 11. These adjusted weights were used to calculate the total degeneracies of the states of the D + 3 line list for which C2v symmetries are used and are included in the new states file.
Partition functions are given in Table 12 for a series of temperatures up to the value of Tmax for each line list. Partition functions for D + 3 had previously been calculated by  and their values are in close agreement with those presented here. For H2D + however, the partition functions differ from those given by Sochi & Tennyson (2010). This is due to the different nuclear spin statistical weights used in their earlier work that did not consider the spin of the Deuterium atom. The degeneracies given in the H2D + line list have been updated to use the nuclear spin statistical weights quoted in Table 1 and are thus consistent with the partition function computed using Eq. (3) and the same convention.

Rotational Spectra
Purely rotational transitions within the vibrational ground states are important for detections of H + 3 and particularly its deuterated isotopologues in the interstellar medium. This is due to the low temperatures in such regions driving the majority of the level populations to low-energy states. Predicted transition frequencies for such transitions are given in Table 13 for H + 3 and D + 3 and in Table 14 for H2D + and D2H + , based on the term energies derived from the marvel networks presented here. While the purely rotational transitions of H + 3 are "forbidden" due to the lack of a permanent electric dipole moment, it is believed that a small, temporary dipole moment can be induced due to the distortion of the molecule from equilibrium geometry under rotation about the C2 molecular axis (Pan & Oka 1986;Miller & Tennyson 1988). The D + 3 states file contains 43 extremely long-lived states with radiative lifetimes greater than 10 10 s and can be considered meta-stable. In the extreme case, two states are calculated to have radiative lifetimes greater than 10 18 s, a duration greater than the current age of the universe. All of these meta-stable states are in the vibrational ground state and have G = J or G = J −1. The same is true for the meta-stable states of H + 3 , though only one of their radiative lifetimes is in excess of 10 10 s. In both species, this can cause molecules to become "trapped" in these states in collision free environments with consequences for both laboratory measurements ) and for possible maser action. Fig. 1 shows example stick spectra for H + 3 and its deuterated isotopologues in the far infrared region from 10 -1000 cm −1 . These spectra were generated using the program ExoCross (Yurchenko et al. 2018) and were computed for a temperature of 200 K. The predicted transitions frequencies listed in Tables 13 and 14 that are expected to be of importance in the ISM fall within this range. Likewise, spectra were computed for the mid infrared region between 1000 -5000 cm −1 for all species considered here and are shown in Fig. 2. These cover several wavelength ranges where H + 3 has been observed in Jupiter's aurorae, such as the K band region centred at 2.2 µm (Uno et al. 2014) and the L band centred on 3.5 µm (Baron et al. 1991). This also covers the wavelength regions observed using the Jupiter infrared Auroral Mapper (JIRAM) instrument onboard NASA's Juno spacecraft (Dinelli et al. 2019;Migliorini et al. 2019) and those covered by the Near infrared Camera (NIRCam) long-wavelength channel onboard JWST.

Mid and Far Infrared Spectra
Comparisons between computed spectra and real observations cannot be made for H2D + , D2H + and D + 3 , as no measured transition intensities have been published. McKellar & Watson (1998) presented infrared absorption spectra of the ν2 fundamental band of H + 3 and provided integrated intensities, which were well reproduced by the original line list (Mizus et al. 2017).

CONCLUSION
We have updated three marvel networks for H + 3 , H2D + and D2H + to include all currently published spectroscopic data for these molecules. We have also performed variational nuclear motion calculations using the program DVR3D for D2H + and D + 3 to produce the new MiZo line lists. The empirical energy levels derived from the marvel networks have been used to update the calculated levels of their respective molecules. This allows for the subset of transitions involving these marvelised energies to be determined to much higher accuracy, making the use of these line lists well suited for high-resolution spectroscopy. The D + 3 calculations have been combined with a set of energies derived from effective Hamiltonian calculations computed using experimentally determined molecular constants. Given these effective Hamiltonian constants were derived from transitions in infrared bands, the new D + 3 line list is best suited for infrared studies. Overall, the new D + 3 line list is of lower resolution than the other three, due to the lack of marvelised energy levels.
Hyperfine effects are known to be present in the spectra of H + 3 and its deuterated isotopologues, due to the nuclear spin of the component protons and deuterons (Jensen et al. 1997). These hyperfine splittings have not yet been observed however in either an astrophysical setting or in the laboratory. Were experiments to be conducted to measure the hyperfineresolved spectra of the species considered here, it would enable the construction of a hyperfine-resolved marvel network and subsequently line lists. Similarly, any further measurements of hyperfine-unresolved spectra could be added to the current marvel network to further constrain the resultant empirical energy levels and hence improve the accuracy of the line lists. If a sufficiently large number of spectra were observed for D + 3 to form a well-connected spectroscopic network, it would enable us to further update the D + 3 line list with marvel energies to a high-resolution standard.
The four line lists presented here, each consisting of a states file and transitions file, are made available via www.exomol.com. Figure 1. Far infrared spectra for H + 3 and its deuterated isotopologues, computed using the program ExoCross (Yurchenko et al. 2018) at 200 K between 10 -1000 cm −1 (1 mm -10 µm).

DATA AVAILABILITY
The updated MARVEL input files and resulting output energy levels are given as supporting material. All other data are available via the www.exomol.com website. Table 14. Predicted pure-rotation transition wavenumbers and Einstein A coefficients for transitions involving low-energy lower levels for H 2 D + and D 2 H + . For all upper and lower levels included here, ν 1 = ν 2 = ν 3 = 0. Values quoted in brackets for the transition frequencies are the uncertainty in the last digit.