© The Authors 2025

1 Introduction

Astrochemical models for interstellar clouds and protoplanetary disks differ in various ways from models simulating the chemistry in planetary and exoplanetary atmospheres. Using the kinetic rate networks originally developed for interstellar clouds makes sense for relatively cold and tenuous outer disks. However, the densities and temperatures in the inner disks around T Tauri stars resemble the conditions in planetary atmospheres where the equilibrium chemistry approach is often used (Marley & Robinson 2015). This raises the questions of where and how these astrochemical models and methods developed separately for interstellar clouds and disks, and planetary atmospheres connect in terms of density and temperatures.

Studies such as Lodders & Fegley (2002) and Fortney et al. (2006) discussed the departure from thermodynamic equilibrium in planetary atmospheres of hot Jupiters and demonstrated that disequilibrium chemistry better explains the observed spectra. Cooper & Showman (2006) studied the effect of disequilibrium chemistry and reported that chemical equilibrium may not exist in the exoplanetary atmosphere due to long timescales at low temperatures and pressure. Carleo et al. (2022) and Guilluy et al. (2022) used thermodynamic equilibrium to analyse the observations of warm Saturn-mass planets but suspected that disequilibrium chemistry can better interpret them. Zahnle et al. (2009a), Zahnle et al. (2009b), Moses et al. (2011), Venot et al. (2014), and Barth et al. (2021) studied (photo-) chemical kinetics. Tsai et al. (2023) used photochemically produced SO₂ to explain the observed feature at 4.05 μm in the transmission spectrum of WASP-39b, Molaverdikhani et al. (2020) demonstrate the effect of quenching in ultra-hot Jupiters, and Moses et al. (2011) and Moses (2014) reported that thermodynamic equilibrium can only prevail in deeper layers of hot planetary atmospheres. Blumenthal et al. (2018) studied whether the differences caused by equilibrium and disequilibrium chemistry are observable across the range of exoplanets and report that disequilibrium chemistry should be considered to determine the metallicity of a planet. Lodders (2003) used thermodynamic equilibrium approach to calculate condensation temperatures for a wide range of gas-phase species and condensates.

One such models developed to simulate the chemistry in planetary and exoplanetary atmospheres is ARGO (Rimmer & Helling 2016, 2019), which introduced the STAND chemical rate network. This network includes a large number of three-body reactions, it has reaction-rate coefficients in the low- and the high-pressure limits, and it uses a thermodynamical concept to reverse reactions based on the temperature-dependent Gibbs free energies of the molecules. This concept is key to ensuring that, at very high pressures, the abundances approach their thermodynamic equilibrium values.

In this paper, we aim to utilise these concepts and implement them in the radiative thermochemical disk model PROtoplanetary DIsk MOdel (PRODIMO) (Woitke et al. 2009; Kamp et al. 2017; Woitke et al. 2016). We merged this exoplanetary chemical network with the chemical network developed for disks and aimed to test the validity of thermodynamic equilibrium in the planet-forming disk environments. Section 2 explains the theoretical concepts adopted from the exoplanet community. Section 3 describes our methodology and the modelling results are shown in Sect. 4. The effect of reversing the reactions on the chemical abundances and analysing whether thermodynamic equilibrium is established in disk models are studied in Sect. 5. We present our conclusion in Sect. 6.

2 Chemical reactions

Here we briefly summarise the theoretical concepts used to determine the rate coefficients for the different types of chemical reactions. Kamp et al. (2017) included photoreactions, a few three-body reactions and chemical concepts from the interstellar medium. We used the photo-cross-sections for the photoreactions from the Leiden database (van Dishoeck et al. 2006; Heays et al. 2017). The bimolecular rates were calculated using the modified Arrhenius equation described by Arrhenius (1889). Below, we explain the concepts introduced now from exoplanet chemistry models.

2.1 Lindemann-Hinshelwood mechanism

There are a few reaction types that are regularly included in planetary chemistry, but are not typically considered in interstellar and disk chemistry because they are irrelevant for the low pressures and temperatures in interstellar clouds. These reactions proceed via an activated complex that is stabilised by another collision with an abundant third gas particle, for example

• (1)

  termolecular A + B + M → AB + M

• (2)

  thermal decomp. AB + M → A + B + M

• (3)

  3-body recomb. A⁺ + e⁻ + M → A + M

• (4)

  thermal ionization A + M → A⁺ + e⁻ + M

where A, B, and AB are atoms or molecules. The (2) and (4) reaction types, which are endothermic and have large activation energies, can be considered as the reverse of the (1) and (3) reaction types, respectively. Another example is the termolecular ion-neutral reaction where A and AB are charged molecules.

The reaction rates of these types of reactions are assumed to follow the Lindemann-Hinshelwood mechanism (Lindemann et al. 1922), which is derived in Atkin & de Paula (2006) as $$K=\frac{K_0{n}_{\mathrm{M}}}{1+K_0{n}_{\mathrm{M}}/K_{\mathrm{\infty }}},$$(1)

where K₀ and K_∞ are the rate constants in the low- and the high-pressure limits, respectively. The unit of K depends on the reaction type and the pressure regime, and n_M is the number density [cm⁻³] of the third body. $$K_{\circ }={\alpha }_0{\left(\frac{T}{300}\right)}^{{\beta }_0}\exp (-\frac{{\gamma }_0}{T})$$(2) $$K_{\mathrm{\infty }}={\alpha }_{\mathrm{\infty }}{\left(\frac{T}{300}\right)}^{{\beta }_{\mathrm{\infty }}}\exp (-\frac{{\gamma }_{\mathrm{\infty }}}{T})$$(3)

In the case of a termolecular or three-body recombination reaction, the unit of K_∘ is cm⁶ s⁻¹ and the units of K_∞ and K are cm³ s⁻¹. In PRODIMO, we assume n_M = n_H2 + n_H + n_H⁺ + n_He + n_He⁺ for the particle density of the third body. There are six constants required to calculate the effective reaction rate, K.

In some works, the reaction rate according to Eq. (1) is multiplied by another factor, F, which is a dimensionless function of temperature and pressure (Troe & Gesell 1983). This is used to obtain a more accurate transition from the low- to the high-pressure limit. We neglect these corrections in this work.

2.2 Backward reactions

We consider an example of a pair of a forward and a reverse gas phase reaction $$\mathrm{A}+\mathrm{B}\leftrightarrow \mathrm{C}+\mathrm{D}+\mathrm{E}.$$(4)

In thermodynamic equilibrium (TE), the law of mass action (see Appendix A) for this reaction reads $$\left(\frac{\stackrel{\circ }{p_{\mathrm{C}}}}{p^{\ominus }}\right)\left(\frac{\stackrel{\circ }{p_{\mathrm{D}}}}{p^{\ominus }}\right)\left(\frac{\stackrel{\circ }{p_{\mathrm{E}}}}{p^{\ominus }}\right)=\left(\frac{\stackrel{\circ }{p_{\mathrm{A}}}}{p^{\ominus }}\right)\left(\frac{\stackrel{\circ }{p_{\mathrm{B}}}}{p^{\ominus }}\right)\exp (-\frac{{\mathrm{\Delta }}_{\mathrm{R}}{G}^{\ominus }}{RT})$$(5)

where ${\stackrel{\circ }{p}}_i={\stackrel{\circ }{n}}_{i}kT$ is the partial pressure, ${\stackrel{\circ }{n}}_i$ the number density of species, i, in thermodynamic equilibrium, and p^⊖ = 1 bar is a standard pressure. Δ_RG^⊖[J/mol] is the change in Gibbs free energy per mol of reactions $${\mathrm{\Delta }}_{\mathrm{R}}{G}^{\ominus }={\mathrm{\Delta }}_{\mathrm{f}}{G}^{\ominus }(\mathrm{C})+{\mathrm{\Delta }}_{\mathrm{f}}{G}^{\ominus }(\mathrm{D})+{\mathrm{\Delta }}_{\mathrm{f}}{G}^{\ominus }(\mathrm{E})-{\mathrm{\Delta }}_{\mathrm{f}}{G}^{\ominus }(\mathrm{A})-{\mathrm{\Delta }}_{\mathrm{f}}{G}^{\ominus }(\mathrm{B})$$(6)

and the Δ_fG^⊖ are the Gibbs free energies of formation of the reactants and products as pure substances from the standard states of the elements at room temperature. In thermodynamic equilibrium, the forward and reverse reaction rates are equal (detailed balance), hence $$K_{\mathrm{f}}{\stackrel{\circ }{n}}_{\mathrm{A}}{\stackrel{\circ }{n}}_{\mathrm{B}}=K_{\mathrm{r}}{\stackrel{\circ }{n}}_{\mathrm{C}}{\stackrel{\circ }{n}}_{\mathrm{D}}{\stackrel{\circ }{n}}_{\mathrm{E}},$$(7)

where K_f and K_r are the rate coefficients of the forward and reverse reactions. Combining Eqs. (5) and (7) results in $$\frac{K_{\mathrm{r}}}{K_{\mathrm{f}}}=\frac{{\stackrel{\circ }{n}}_{\mathrm{A}}{\stackrel{\circ }{n}}_{\mathrm{B}}}{{\stackrel{\circ }{n}}_{\mathrm{C}}{\stackrel{\circ }{n}}_{\mathrm{D}}{\stackrel{\circ }{n}}_{\mathrm{E}}}={\left(\frac{p^{\ominus }}{kT}\right)}^N\exp (+\frac{{\mathrm{\Delta }}_{\mathrm{R}}{G}^{\ominus }}{RT}),$$(8)

where N, in general, is the number of reactants minus the number of products; here N = −1.

Equation (8) implies that the reverse and forward rate coefficients are not independent and are related to each other via the reaction’s Gibbs free energy. Using this principle is standard in exoplanet chemistry, but is normally not used in interstellar cloud or disk chemistry. The concept can be applied not only to all reaction types listed in Sect. 2.1, but also to regular bimolecular neutral-neutral, ion-neutral, and ion-ion gas phase reactions that have no high-pressure limit data.

Equation (8) can be reliably applied to the case where the forward reaction is an exothermic reaction, Δ_RG^⊖ < 0, and the reverse reaction rate coefficient will be small because of the negative term in the exponential in Eq. (8). However, there are cases where the rate measurements are only available for an endothermic reaction (for example, a termolecular reaction) and the reverse exothermic reaction needs to be constructed. This is still safe as long as the measured activation energy, γ₀, in Eq. (2) is larger than the positive Δ_RG^⊖. However, if that is not the case, the constructed rate would go to infinity at low temperatures. This problem is further discussed in Sect. 3.

2.3 Thermochemical data

The thermochemical data, such as the standard formation enthalpies, Δ_fH^⊖, and the standard Gibbs free energies of formation, Δ_fG^⊖, are calculated from NASA seven-term polynomials given by Burcat & Ruscic (2005). These consist of 14 coefficients. The first seven coefficients are used to determine the thermodynamic properties for temperatures 1000−6000 K and the subsequent seven coefficients are used for 200−1000 K.

We carefully checked the standard heat of formation, Δ_fH^⊖(0 K), against the data provided by Millar et al. (1997), the National Institute of Standards and Technology (NIST) Standard Reference Database Number 69¹, and the Active Thermochemical Tables (ATcT)². In most cases, we find fine agreements between the four data sources, but also some disagreements. About 30 molecules in our selection could not be found in Burcat & Ruscic (2005) and Dick (2019). In these cases we (i) fitted the thermochemical data found in Gleich-Gewichts-Chemie (GGchem) and exported them in Burcat format, (ii) estimated the Burcat polynomials from similar molecules calibrated with the correct Δ_fH^⊖(0 K), as found in the ATcT and NIST databases, or (iii) used known ionisation potentials, or (iv) used proton affinities with respect to their mother molecules. The details are explained in Table B.1.

3 Implementation

Here, we describe the implementation of the above-mentioned concepts in the chemical network of the thermochemical disk code P_roD_IM_o (Woitke et al. 2009; Woitke et al. 2016).

3.1 Selection of species

To study whether thermodynamic equilibrium can be established in disks we made a selection of species based on the large DIANA (DIscANAlysis) chemical network (Kamp et al. 2017). However, since we discuss rather high temperatures of T ≳ 500 K, we eliminated all ice species. We have furthermore eliminated excited molecular hydrogen, ${\text{H}}_2^{\ast }$, and polycyclic aromatic hydrocarbons (PAHs), due to the lack of thermochemical data for those, resulting in the set of 169 species listed in Table 1. All but the ten doubly ionised atoms have valid thermodynamic data. It is noteworthy that our selected species included atomic and molecular ions, which is not standard in exoplanet chemistry, to study the details of photoionisations and cosmic rays, and doubly ionised atoms to study the effects of X-ray processes implemented by Aresu et al. (2011). We reintroduced the eliminated species later in Sect. 5 to study the 2D abundance structure in full disk models.

3.2 Reaction data and pairs of forward and reverse reactions

Our primary source of reaction kinetic data was the UMIST 2022 database (Millar et al. 2024). Additional gas-phase reactions were taken from the STAND network (Rimmer & Helling 2016), providing the majority of three-body reactions with low- and high-pressure limits. The STAND reactions have preference in this work: that is, if the same reaction was found in both databases, we used the kinetic reaction data from STAND. This includes the radiative association and direct recombination reactions, but not the photoionization and photodissociation reactions, which are treated in a different way in PRODIMO. Only reactions among the selected species were read from these databases. In STAND, the reactions are already paired. They come in the form of a reaction with valid Arrhenius coefficients followed by the reverse reaction with invalid data. We ignored this pairing and simply added each reverse reaction with a flag signalling that the reaction rate of this reaction needed to be derived later from its respective forward rate, using Eq. (8).

Next, each gas-phase reaction is checked whether a reverse reaction is present, in which case these two reactions are paired. If an invertible reaction is found that has no reverse counterpart, the reverse reaction is auto-created with a flag that it has no valid data, then we pair these two reactions. However, a STAND reverse reaction, which has invalid reaction data, will not overwrite the same UMIST 2022 reaction with valid data. A gas-phase reaction is considered to be invertible if their reactants and products all have valid thermodynamic data. However, gas-phase reactions with less than two reactants (spontaneous reactions) or more than three products are not inverted. Once the pairing of the reactions is completed, we check whether at least one reaction in each pair has valid kinetic data. We also check whether any unpaired invalid reactions remain, in which case the program stops with an error message.

Altogether, the reaction network merged in this way consisted of 5539 (2462+375+1155+1547) reactions. Of these, 2462 reactions were from the UMIST 2022 database and 375 were taken from a collection of additional reactions compiled in the past by various PRODIMO developers (see Aresu et al. 2011, Kamp et al. 2017 and Kanwar et al. 2024b), resulting in 2837 reactions. This compilation now also includes some additional three-body and neutral-neutral reactions of silicon hydrides, as listed in Table 2. The UMIST 2022 database only has ion-neutral reactions for silicon hydrides; see further discussion in Sect. 4. There were 1155 reactions from the STAND network and 736 existing reactions were overwritten. Finally, 1547 reverse reactions were auto-created as discussed above, resulting in 2466 reaction pairs. A total of 75 gas-phase reactions out of 5539 reactions remained unpaired.

The final network included 135 X-ray driven reactions, 252 photoreactions, and 203 cosmic-ray driven reactions. The photoreactions here also included radiative association and direct recombination reactions. The cosmic-ray reactions included cosmic-ray particle reactions and ionisation by UV photons generated by the cosmic rays. Table 3 provides an overview of the different types of reactions picked up from the STAND network. The number of forward reactions included exceeds the number of reverse reactions because for some of the STAND reverse reactions with invalid reaction data the same reaction is found in UMIST 2022 with valid reaction data.

3.3 Computation of the rate coefficients

To calculate the rate coefficients of the paired reactions, we first determined the reaction enthalpy, Δ_RH^⊖(T), from NASA seven-term polynomials to decide which one was the exothermic reaction. If that reaction had valid data, it was considered as the “forward” reaction and the paired endothermic reaction was the reverse. However, if the exothermic reaction had no valid rate coefficient data, it became the reverse reaction, and the endothermic reaction became the forward reaction. In both cases, the rate coefficient of the reverse rate was finally computed according to Eq. (8). If the exothermic reaction did not have reaction data and the endothermic reaction was identified as the forward one, we have Δ_RG^⊖ > 0 and Eq. (8) can potentially lead to K_r → ∞ for T → 0. Combining Eq. (8) with Eq. (2) for the forward rate, we have $$K_{\mathrm{r}}={\left(\frac{p^{\ominus }}{kT}\right)}^N{\alpha }_{\mathrm{f}}{\left(\frac{T}{300}\right)}^{{\beta }_{\mathrm{f}}}\exp (+\frac{{\mathrm{\Delta }}_{\mathrm{R}}{G}^{\ominus }/R-{\gamma }_{\mathrm{f}}}{T}).$$(9)

Thus, there is a problem when Δ_RG^⊖/R > γ_f. This case seems chemically impossible, as an endothermic reaction should have an activation energy that is at least as large as the reaction enthalpy. However, in practice this case occurs frequently whenever an endothermic reaction has been reported with low or zero activation energy. Tinacci et al. (2023) identify this problem as one of the major weaknesses of current theoretical astrochemical modelling and propose a cleaning of chemical networks by simply eliminating all these spurious reactions. Here, we followed the solution that Rimmer & Helling (2016) propose, to multiply both K_f and K_r by the same factor $\exp (-\frac{\frac{{\mathrm{\Delta }}_{\mathrm{R}}{G}^{\ominus }}{R}-{\gamma }_{\mathrm{f}}}{T})$ if Δ_RG^⊖/R > γ_f, which not only avoids the problem above, but also reduces the rate coefficients of all spurious endothermic forward reactions automatically (auto-cleaning).

3.4 At which densities do three-body reactions become important?

Computing an average over all exothermic three-body reaction rates in the new merged chemical network at 1000 K results in $⟨{\log }_{10}{K}_{3\text{-body }}^{\text{exotherm }}\left[{\mathrm{c}\mathrm{m}}^6{\mathrm{s}}^{-1}\right]⟩=-31.9\pm 4.8$, whereas the average over all exothermic two-body reaction rates is $⟨{\log }_{10}{K}_{2\text{-body }}^{\text{exotherm }}\left[{\mathrm{c}\mathrm{m}}^3{\mathrm{s}}^{-1}\right]⟩=-12.1\pm 3.4$. From this simple consideration one might be led to the conclusion that three-body reaction rates are only important when the density is as large as 10²⁰ cm⁻³. However, this conclusion is incorrect. The termolecular reactions open new chemical pathways with reactions that are not possible without a stabilising collision, and in our chemical analysis we see three-body reactions being flagged as important once the density is of the order of 10¹³ cm⁻³ or higher. The following consideration tries to explain this observation.

In order to assess the density beyond which the inclusion of the three-body (termolecular) reactions seems important, we searched in the network for pairs of three-body and two-body reactions that share the same reactants, except for the additional M in the three-body reaction. An example for such reactions is listed in Table 6: see the first two destruction reactions listed for C₂H₂. From the rate coefficients of such a pair we can calculate a critical density for the third body as $$n_{\mathrm{c}\mathrm{r}}^{\mathrm{t}\mathrm{e}\mathrm{r}}=\frac{K_{2-\mathrm{b}\mathrm{o}\mathrm{d}\mathrm{y}}}{K_{3-\mathrm{b}\mathrm{o}\mathrm{d}\mathrm{y}}}.$$(10)

Each of these pairs of reactions means a branching point. Either the activated complex is stabilised by a collision with a third body, resulting in a constructive exothermic process, or it decays spontaneously into two fragments, often endothermically. In this way, the exothermic termolecular reactions often compete with the endothermic two-body reactions which have much lower rates than their exothermic counterparts. From Eq. (1), one can derive another critical density beyond which the high-pressure limit becomes relevant: $$n_{\mathrm{c}\mathrm{r}}^{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}-\mathrm{p}}=\frac{K_{\mathrm{\infty }}}{K_0}.$$(11)

Table 4 lists the mean values we find for these two critical densities in our network, at various temperatures. The critical density for the inclusion of termolecular reactions, $n_{\text{cr}}^{\text{ter}}$, shows a broad distribution among the molecules, but values from roughly 10¹¹ cm⁻³ to 10¹⁵ cm⁻³ are often found. The values are lower at low temperatures, where the termolecular reactions are often faster than their endothermic two-body counterparts that have higher activation barriers. In contrast, the critical density for the inclusion of high-pressure data varies from about 10¹⁷ cm⁻³ to 10²² cm⁻³, depending on reaction and temperature, which seems unlikely to be important for any region in protoplanetary disks. Thus, including the three-body termolecular reactions seems important for the inner disks, but using only the low-pressure reaction data is justified. Kamp et al. (2017) and Kanwar et al. (2024b) already include a limited set of three-body reactions; however, the network developed here is more exhaustive for three-body reactions.

3.5 Disk modelling

We used the thermochemical disk code PRODIMO (Woitke et al. 2009; Woitke et al. 2016) to establish in which parts of the disk thermochemical equilibrium can be established. This code first sets up an axisymmetric physical gas and dust density structure, and then calculates the 2D continuum radiative transfer to determine the dust temperature structure, T_dust(r, z). PRODIMO then iteratively calculates the chemical and thermal structure of the gas in the disk using the enhanced chemical network described in Sect. 3.2. The gas temperature is varied at every point in these calculations until the local gas heating rate is balanced by the local gas cooling rate, taking into account heating and cooling processes (see Woitke 2015). The treatment of cosmic rays follows the approach described by Padovani et al. (2009), Padovani et al. (2013) and Rab et al. (2017, see Eq. (4)).

As a basis for all modelling steps to follow, we computed a standard T Tauri PRODIMO disk model adopted from Woitke et al. (2016) with the enhanced chemical network described above. The full list of input parameters of the disk model are provided in Appendix C. The elemental abundances were adopted from Kamp et al. (2017). The C/O elemental ratio used in this work is 0.45.

4 Results

4.1 Verification of the chemical network

Finally, we obtain a chemical network that we named CHemistry Assembled from exoplanets and dIsks for Thermodynamic EquilibriA (CHAITEA), which consists of the UMIST 2022, large DIANA, and STAND network. Our first goal was to verify whether the CHAITEA chemical network was able to converge towards thermodynamic equilibrium. We checked this by selecting a single point from the T Tauri disk model, with given hydrogen nuclei density, n_〈H〉, and only simulating the chemistry at that point either in the time-dependent or in the time-independent mode. For the purpose of this test, we eliminated all possible causes for deviations from thermodynamic equilibrium: (i) we switched off the cosmic-ray and X-ray reactions by using extremely low values for the cosmic-ray ionisation rate and the stellar X-ray luminosity, (ii) we set the dust and gas temperatures equal to a user-specified value T_gas = T_dust = T and, (iii) we set the local radiation field to a Planckian J_v = B_v(T), (iv) we switched off all photo rates including any reactions that produce photons, such as radiative associations and direct recombinations. We defined two conditions to claim that the chemistry had reached thermodynamic equilibrium: (1) the forward and reverse rates of all reaction pairs are equal (detailed balance) and (2) the calculated abundances are equal to the abundances computed with an independent code that predicts them in thermodynamic equilibrium, $n_i={\stackrel{\circ }{n}}_i$. For the latter, we used GGchem (Woitke et al. 2018) with a new option to use the same Burcat thermochemical data as used in this work. For our selection of 12 elements, GGchem finds 1450 molecules in its databases, including many isomers and molecules as complex as [(O₂N)₃C₆H₂CH]₂ (hexanitrostilbene). However, for all temperatures considered in this paper, these complex molecules only attain trace abundances.

We produced three grids of models with varying density and temperature. These models and their parameters are outlined in Table 5. All three grids explore the same density and temperature range. First, we produced grid_0, where X-rays, the photoreactions and cosmic rays are turned off. Second, we studed the deviations caused by photoreactions based on a local Planck field J_v = B_v(T) and reactions producing photons, such as radiative association and direct recombination reactions that are considered non-invertible, leading to another grid of 6 × 27 models called grid_photo. Third, to discern the impact of the cosmic rays, we again switched off the photorates and produced a third grid of models for six different cosmic-ray ionisation rates (CRIs), leading to grid_cr with 6 × 6 × 27 models.

The CHAITEA chemical network converges towards thermodynamic equilibrium in the absence of radiation such as cosmic rays, X-ray, and UV. We arrive at this conclusion irrespective of the molecule considered. Table 6 shows the detailed balance obtained for C₂H₂ solving the rate network in the time-independent mode at T = 1200 K and n_〈H〉 = 7.56 × 10¹³ cm⁻³ as part of grid_0. It also shows that species such as H₂ O and OH play a role in the formation of C₂H₂, similarly to the findings of Kanwar et al. (2024b). Figure 1 shows the calculated abundances of various species, determined with both PRODIMO and GGchem, at a density of n_〈H〉 = 7.56 × 10¹³ cm⁻³. These results match well for temperatures >550 K; below that, they begin to deviate from each other. These deviations are due to numerical problems in finding the correct chemical solution with the time-independent solver used in PRODIMO. However, in general, our results show that the CHAITEA chemical network will indeed approach thermodynamic equilibrium when all causes to deviate from it are eliminated.

Fig. 1 Abundances of various species, ϵ(X)=nX/n〈H〉, obtained with the time-independent PRODIMO models from grid_0 for density n〈H〉= 7.56 × 1013 cm−3 (dots), compared to the corresponding GGchem models (dashed lines). The deviations at the lowest temperatures are due to numerical problems in converging the PRODIMO models.

4.2 Chemical relaxation timescale

We now analyse the chemical relaxation timescale (τ_chem) in the disk where the kinetic chemistry approaches thermodynamic equilibrium. This can tell us whether thermodynamic equilibrium can be established on disk-evolutionary timescales. We realise that there can be various definitions for chemical relaxation timescales and explore them here. We define the chemical relaxation timescale in the following two different ways: $${\tau }_{\text{chem},1}=\underset{\text{valid}i}{max}\{n_i/min(F_i,D_i)\}$$(12) $${\tau }_{\text{chem},2}=\underset{\text{valid}n}{max}\left|\mathrm{\Re }{\left\{{\lambda }_n\right\}}^{-1}\right|,$$(13)

where F_i and D_i are the total formation and destruction rates, respectively, of species i. τ_{chem, 1} can be computed from the diagonal elements in the chemical Jacobian. The slowest species determines the relaxation timescale. Equation (12) is easy to understand but lacks a proper mathematical foundation. Using perturbation theory and linear algebra, the proper mathematical formulation of the relaxation timescale involves a transformation of the Jacobian into n eigenvectors (the chemical modes) which each relax on individual timescales given by the inverse of the real part of their eigenvalues, λ_n. The slowest of these modes determines the chemical relaxation timescale, τ_{chem, 2}, according to Eq. (13), (see Woitke et al. 2009), with values listed in Table 7. Since a chemical network obeys a number of auxiliary conditions in the form of the conservation of the elements, the major element carriers must not be considered in Eq. (12), and the N_el slowest modes must be discarded in Eq. (13), where N_el is the number of elements.

Figure 2 depicts the evolving abundances of various species in two time-dependent chemical models at 1200 K and 2000 K using a density of 7.65 × 10¹³ cm⁻³. We can visually read from this figure that the PRODIMO abundances only approach the GGchem values after ∼10⁸ yr(τ_{chem, 3}) in the 1200 K model, but already approach them after ∼0.1 yr in the 2000 K model. These examples show how strongly the chemical relaxation timescale increases with decreasing temperature. While all three timescales at 2000 K are relatively similar, the timescale at 1200 K from the first definition (τ_{chem, 1}) is about four orders of magnitude lower than τ_{chem, 2} and the visually determined value from the figure, which are both similar.

Hence, τ_{chem, 1} can be misleading. For example, proton exchange reactions between N₂ and N₂H⁺ have high rates, yet these reactions do not contribute to the formation and destruction of the strong N ≡ N bonds. Therefore, studying the chemical relaxation timescale of the chemical modes (eigenvectors of the Jacobian) is more reliable than the first definition. However, even that method does not predict the chemical relaxation timescale fully satisfactorily in all cases, especially when we consider even lower temperatures where the chemical relaxation timescale can easily become greater than the age of the universe. One of the possible explanations is that there are numerical problems when decomposing the chemical Jacobian into its eigenvectors. This is due to the extremely high conditional numbers of the matrix. Thus, from this analysis, we conclude that it is not trivial to reliably evaluate the chemical relaxation timescale that can have a huge range of fast and slow modes.

During this analysis, we find that species such as SiH take extremely long to reach their equilibrium abundances based on a preliminary network. This indicates that this preliminary chemical network is incomplete for silicon-hydrides, as both UMIST 2022 and STAND only include reactions of Si-hydrides with charged species but no neutral-neutral reactions. We therefore added a few reactions from Raghunath et al. (2013) and NIST (Manion et al. 2015), as listed in Table 2, to the CHAITEA chemical network. This leads to a much faster relaxation of the Si-hydride molecules. However, once this problem is solved, the next molecules that present this issue are N-bearing species, which is further discussed in the following section. This might indicate that the network for these species is also not yet complete.

4.3 The N₂ formation problem

Figure 2 shows that N-bearing species such as HCN, NH₃, and CN are the slowest to attain thermodynamic equilibrium at 1200 K, which happens only after ∼10⁸ years. We analyse here the time-dependent behaviour at various steps for the two temperatures 1200 K and 2000 K. We started our simulations by putting all elements into their highest ionisation state as an initial condition, except for neutral hydrogen, which has an initial abundance of 10⁻⁵. The gas recombines within 10⁻⁸ yr (i.e. in less than one second) at both temperatures. After about 10⁻⁶ yr (∼30 s), neutral molecules form by consuming the neutral atoms via barrier-free, exothermic neutral-neutral and termolecular reactions such as $$\mathrm{H}+\mathrm{H}+\mathrm{H}\to {\mathrm{H}}_2+\mathrm{H}$$(14) $${\mathrm{H}}_2+\mathrm{O}\to \mathrm{O}\mathrm{H}+\mathrm{H}$$(15) $$\mathrm{O}\mathrm{H}+\mathrm{N}\to \mathrm{N}\mathrm{O}+\mathrm{H}$$(16) $$\mathrm{N}\mathrm{O}+\mathrm{N}\to {\mathrm{N}}_2+\mathrm{O}$$(17)

This phase results in the creation of many simple neutral molecules that all attain similar abundances, as they are all produced by almost equally fast reactions. Most of the atoms are consumed after about 0.1 yr (one month) at 1200 K and about 0.0001 yr (one hour) at 2000 K. At this time, the chemistry slows down substantially. Some parts of the network, such as the hydrocarbon chemistry, have already reached internal detailed balance, but as these parts are only weakly connected to other parts of the network, such as the N- and O-chemistry, the connections have not yet reached equilibrium. Eventually, after about 0.1 yr in the 2000 K model, each molecule has reached thermodynamic equilibrium. However, in the 1200 K model, after 10³ yr all molecules have reached thermodynamic equilibrium except for the N-bearing species, which take as long as 10⁸ yr to reach equilibrium. This final phase is characterised by difficulties in forming the thermodynamically favoured N₂ from the other overabundant N molecules. All the reactions stated above become extremely inefficient due to the lack of free atoms at 1200 K. The following reaction is found to be most effective in reducing the N-bearing species and forming N ≡ N bonds: $$\mathrm{N}\mathrm{O}+{\mathrm{N}\mathrm{H}}_2\to {\mathrm{N}}_2+{\mathrm{H}}_2\mathrm{O},$$(18)

but since both NO and NH₂ are not very abundant molecules in thermodynamic equilibrium, this reaction channel takes as long as 10⁸ yr in this example.

This analysis indicates that we may miss some other neutral-neutral reaction pathways to form or destroy N ≡ N bonds, in order to reduce the chemical relaxation timescale. We experimented with adding N₂H, NCN, and HNCO to our network, which indeed leads to a slight decrease in the relaxation timescales, but only by a factor of 2−3, and we are unable to solve the general problem this way. One possibility might be to include surface reactions to form N₂ on much shorter timescales. At temperatures ⩽ 1100 K, the chemical relaxation timescale derived from Eq. (13) becomes longer than the age of the universe. This issue was also reported by Le Gal et al. (2014). The neutral-neutral reactions forming N₂ mentioned by Pineau des Forets et al. (1990) and Hily-Blant et al. (2010) are included in the chemical network but are not found to be dominant.

Fig. 2 Abundances, ϵ(X)=nX/n〈H〉, of various species (X) as a function of time obtained with Pro DI Mo at 1200 K and 2000 K at a density of n〈H〉 = 7.56 × 1013 cm−3. The large coloured dots represent the abundances in thermodynamic equilibrium obtained with GGchem.

4.4 Deviations from thermodynamic equilibrium

Figure 3 depicts the deviations from thermodynamic equilibrium at different temperatures and densities when photoprocesses, X-rays, and cosmic rays are switched off (grid_0). We define the deviation as $$\sigma =\sqrt{\frac{\sum _{i=1}^N{({\log }_{10}{n}_i-{\log }_{10}{\stackrel{\circ }{n}}_i)}^2{w}_i}{\sum _{i=1}^N{w}_i}},$$(19)

where n_i is the PRODIMO particle density and ${\stackrel{\circ }{n}}_i$ is the GGchem particle density. N is the number of species found in both models. To avoid picking up differences between particle densities at extremely low abundances such as 10⁻⁸⁰, which PRODIMO does not predict properly, we used a summation weight of $w_i={(max\{n_i,{\stackrel{\circ }{n}}_i\}/n_{⟨\mathrm{H}⟩})}^{0.2}$. In this way, a deviation of one order of magnitude at an abundance of 10⁻²⁰ is evaluated as as poor as a deviation of 10⁻⁴ at unit abundance. The mean deviation, σ, varies across temperature and density. We excluded Mg, Mg⁺, Fe, Fe⁺, Na, and Na⁺in Eq. (19) because the selection of Mg-, Na- and Fe- bearing species is incomplete in PRODIMO, causing large deviations at low temperatures where GGchem predicts these metals to be present in the form of hydroxides such as Mg(OH)₂, Fe(OH)₂, and (NaOH)₂, respectively.

As shown in Fig. 3, the deviations are <1% across all densities and temperatures considered, implying that our new CHAITEA chemical network is capable of relaxing towards thermodynamic equilibrium, as expected. These deviations are caused by small differences in the thermochemical data between PRODIMO and GGchem, and due to the different selections of molecular species.

The species that contributes the most to σ is marked by the white labels in Fig. 3. We do not observe any ions because these models are in thermodynamic equilibrium and the deviations are caused mostly by the abundant neutrals. At the lowest temperatures, there are a few cases where the PRODIMO models do not converge, as indicated by the blank rectangles. The deviations are much larger there, only for numerical reasons. In all the cases where the PRODIMO models converge, the deviations from thermodynamic equilibrium are found to be small (<1%). Hence, we consider our first test to be passed: the CHAITEA network relaxes towards thermodynamic equilibrium down to ∼600 K in the absence of any processes causing deviations from it.

Fig. 3 Deviation σ (see Eq. (19)), between particle densities obtained with time-independent PRODIMO (ni) and GGchem $({\stackrel{\circ }{n}}_i)$ at various temperatures and densities from grid_0. The blank rectangles at the bottom mark the models which do not converge. The numbers on the x-axis correspond to the ascending list of densities in Table 5. The central white labels mark the species that contribute most to σ.

4.5 Effect of photoreactions

We relaxed condition number (iv) (see Sect. 4.1), leading to grid_photo, and perform a similar analysis to study the effect of photoreactions on thermodynamic equilibrium. As shown in Fig. 4, the maximum deviation, σ, is then of the order of 10⁻¹. The radiative association reactions become relatively more significant at low densities, causing larger differences in the particle densities between GGchem and PRODIMO. The deviations at higher densities are caused by ions such as H₃ O⁺ and H⁻, at all temperatures. As the ions are more abundant in PRODIMO relative to GGchem, this results in larger deviations from thermodynamic equilibrium.

There is a gradual increase in σ going towards low densities. This is because the photoreactions scale with n_〈H〉. Higher temperatures are affected more as the rates of these photoreactions increase with temperature (Planck radiation field). Reactions such as $$\mathrm{H}+{\mathrm{e}}^{-}\to {\mathrm{H}}^{-}+h\nu$$(20)

cause species such as e⁻ and H⁻to deviate from detailed balance. H₃O⁺ is not in perfect detailed balance due to $${\mathrm{S}\mathrm{i}\mathrm{O}\mathrm{H}}^{+}+{\mathrm{H}}_2\mathrm{O}\to {\mathrm{H}}_3{\mathrm{O}}^{+}+\mathrm{S}\mathrm{i}\mathrm{O}$$(21)

which is an auto-generated reaction (see Sect. 3.2).

Fig. 4 Same as Fig. 3 but from grid_photo which includes photo processes based on a local Planck field Jv = Bv(T), radiative association, and direct recombination reactions.

4.6 Effect of cosmic rays

In this section, photo rates were again turned off. The deviations, σ, are presented in Fig. 5 for various CRIs at different densities and temperatures. The deviations increase with increasing CRI.

Across all densities and CRIs, σ increases strongly with decreasing temperature. This is because the CRIs are independent of temperature, whereas the neutral-neutral gas-phase reaction rates decrease strongly with temperature because of the activation barriers. A secondary trend is that σ increases for lower densities. This effect is because the CRI scale with n_〈H〉, whereas the gas-phase chemical reactions scale with $n_{⟨\mathrm{H}⟩}^2$ or $n_{⟨\mathrm{H}⟩}^3$.

As the CRI increases, the ions are affected first as they are the primary species created by the cosmic rays. Eventually, the neutrals are also affected at all temperatures and densities. High levels of cosmic rays create more reactive radicals such as OH, and CH₃, which promote neutral-neutral chemistry. This results in neutrals deviating from thermodynamic equilibrium. For instance, at 1000 K and a density of n_〈H〉 = 7.562 × 10¹³ cm⁻³, we observe that with increasing CRI, the molecule that causes the largest deviations shifts from He⁺ to NH₃. At a CRI of 1.7 × 10⁻¹⁷ s⁻¹, which is typically used in disk models, σ is of the order of 1 at high temperatures (>1100 K), but as large as 30 at low temperatures; i.e. the abundances in kinetic chemical equilibrium deviate from thermodynamic equilibrium by 30 orders of magnitude on average.

Fig. 5 Impact of cosmic-ray ionisation rate on deviations from thermodynamic equilibrium at different densities and temperatures from grid_cr. The adopted values of CRIs are listed at the top of each panel. The x- and y-axes are the same as in Fig. 3.

4.7 The chemical deviations in a full disk model

In this section, we study the abundances in a 2D disk model around a T Tauri star, and evaluate the differences between the chemical solutions obtained with PRODIMO in kinetic chemical equilibrium and with GGchem in thermodynamic equilibrium. The disk model is outlined in Sect. 3.5 and our choices of stellar, disk shape, dust size and opacity parameters are listed in Table C.1. We excluded ices, PAHs and excited H₂^* for this comparison, because these species are not included in GGchem. The model has a grid resolution of 200 × 300 and includes UV, cosmic rays and X-ray chemistry. The second disk model, called GGchem, uses GGchem and the densities and gas temperatures calculated in the PRODIMO model.

Figures 6 and 7 show the abundances of various molecules from these models, demonstrating how different these solutions are. As Woitke et al. (2021) discuss, the molecular abundances in thermodynamic equilibrium at low temperatures (T < 500 K) have a binary character: either a molecule is very abundant and is the major carrier of an element or it only attains a trace abundance. This is because at low temperatures the contribution of entropy term is low and the code maximises the binding energy of the molecules. In the case of the elemental abundances considered in this paper, these molecules are H₂, H₂O, CH₄, NH₃, H₂S, SiH₄, Mg(OH)₂, Fe(OH)₂, and (NaOH)₂. All other molecules have extremely low abundances in the GGchem model in the outer disk regions where the temperature is low. Molecules such as CO, CO₂, and N₂ are only present in the GGchem model in a relatively thin intermittent surface layer where the temperatures are ∼1000 K due to the heating by UV processes. Temperatures exceed 5000 K above this layer, which is too high for molecules to exist in high abundances.

The abundances calculated with chemical kinetics by the PRODIMO model are different (Figs. 6, 7, left column). The cosmic rays constantly create small amounts of ions, atoms, and neutral radicals, that drive all the important gas-phase reactions and prevent the chemistry from relaxing towards thermodynamic equilibrium.

In particular, oxygen is stored in H₂ O in thermodynamic equilibrium, whereas it is distributed among H₂O, CO₂, and CO in kinetic equilibrium. Carbon is locked in CH₄ in thermodynamic equilibrium throughout the cold disk, whereas it is found in the form of CO₂ and CO in kinetic equilibrium, similarly to in Bergin (2009). The main carrier of N in thermodynamic equilibrium is NH₃ in the entire cold disk, but switches to N₂ in the warm disk surface. In kinetic chemical equilibrium, N mostly forms N₂ with HCN and NH₃ in the inner disk (≲0.5 au), as seen in Fig. 7.

Figure 8 shows the average chemical deviation, σ, between the two models as a function of disk radius, r, and height over the midplane, z. The minimum σ is ∼1, indicating that thermodynamic equilibrium likely does not exist in any disk region. The inner midplane regions of the disk, directly behind the inner rim inside about 0.1 au, are closest to thermodynamic equilibrium, because radiation such as cosmic rays is shielded here. As we go either vertically up or radially outwards in the disk, σ becomes larger. These large σ are due to the effect of irradiation by UV, X-rays, and cosmic rays. The outer midplane is affected by cosmic rays whereas UV and X-rays affect the chemistry in warm surface layers and prevent relaxation towards thermodynamic equilibrium. There is, however, the warm molecular intermittent surface layer where σ again reaches values as low as 5. Here, the gas-phase chemical reactions can compete with the cosmic-ray, X-ray, UV-induced chemistry, because temperatures are high enough to overcome the activation energies of the neutral-neutral reactions. However, the chemistry is still not in thermodynamic equilibrium. We also tested another criterion where σ is only affected when $(max\{n_i,{\stackrel{\circ }{n}}_i\}/n_{⟨\mathrm{H}⟩})>{10}^{-13}$. However, we still observe a similar deviation between the kinetic chemical equilibrium and the thermodynamic equilibrium.

Fig. 6 Abundances of selected species in the 2D disk model calculated using chemical kinetics (left) and thermodynamic equilibrium (right). The white contours correspond to the tick values on the colorbars.

Fig. 7 Continued from Fig. 6. We note the different scaling of the GGchem results for HCN and C2H2.

Fig. 8 Deviation (σ) between kinetic chemical equilibrium and thermodynamic equilibrium in the 2D disk model. σ = 1 denotes a deviation of one order of magnitude between the PRODIMO and GGchem abundances.

5 Implications for disk kinetic chemistry

To assess the implications of the above work for kinetic chemical networks, we studied the effect of extending the large DIANA standard network (Kamp et al. 2017) to the new CHAITEA chemical network on the 2D abundances in thermochemical disk models. We investigated this question in two steps; see Table 8. First, we added the reverse reactions based on the same Gibbs free energies as discussed in this paper, and second, we added the STAND reactions which have preference over the UMIST 2022 and DIANA reactions.

In all these models, the same 239 species were selected. We added excited molecular hydrogen ${\mathrm{H}}_2^{\star }$, five charging states of PAHs, and the ices of all neutral species to the list of species in Table 1. Compared to the original large DIANA (Kamp et al. 2017) standard network consisting of 235 species, there was also N₂O, N₂O⁺, HNO₂⁺, and N₂O#, as we needed N₂ O for the Si hydride reactions (Table 2). Following the rules laid out by Kamp et al. (2017) for developing a network, we also added N₂O⁺ and HNO₂⁺. This was our final CHAITEA network. The disk model was the same as in the previous section, see Table C.1. The gas temperature was consistently calculated with the chemistry in each of the following models. UV-photo, X-ray and cosmic rays processes were included.

5.1 Effect of adding the reverse reactions

Table 8 shows that the inclusion of the reverse reactions mostly increases the number of endothermic reactions, as most of these are derived from their exothermic counterparts. It also “auto-cleans” the network by correcting spurious endothermic reactions with zero or extremely low activation energies, which are unrealistic, as discussed in Sect. 3.3.

Our results are presented in Figs. D.1, D.2, and D.3. The left columns in these figures show the model with the large DIANA and UMIST 2022 reactions, whereas the middle column also includes the reverse of these reactions. By comparing the left and middle columns in these figures, we observe that the abundances of most molecules are only slightly affected (for example, H₂O). However, certain species show significant differences, particularly atomic H close to the outer midplane, and C₂H₂, CO, CO₂, and N₂ close to the inner midplane. For these species, the abundances decrease significantly when reverse reactions are included. Additionally, we observe an increase in N₂H⁺ abundance in several disk regions when reverse reactions are considered.

The relatively high abundance of atomic H in the outer midplane of the DIANA-standard model is an example of artefacts caused by spurious endothermic reactions with zero activation energy. One such reaction is $$\mathrm{A}\mathrm{r}+{\mathrm{H}}_3^{+}\to {\mathrm{A}\mathrm{r}}^{+}+\mathrm{H}+{\mathrm{H}}_2.(\text{erroneous})$$(22)

This reaction was included as one of the additional reactions collected by Kamp et al. (2017), referencing Anicich (1993). However, there was an error in its inclusion; the correct reaction is $$\mathrm{A}\mathrm{r}+{\mathrm{H}}_3^{+}\to {\mathrm{A}\mathrm{r}\mathrm{H}}^{+}+{\mathrm{H}}_2.(\text{correct})$$(23)

The erroneous reaction has a reaction enthalpy of about +6.5 eV. While the correct reaction is also endothermic, with +0.43 eV, it is reported as barrier-free with zero activation energy (see Johnsen 2005 for details). In the outer midplane, where C-, N-, O-, and S-bearing species are locked in their respective ices, H₃⁺ is the most abundant molecular cation, maintaining charge balance with free electrons (Balduin et al. 2023). As a result, H₃⁺ is relatively abundant, as is Ar, which leads the model with the erroneous reaction to produce a significant amount of neutral H in the midplane. However, when reverse reactions are included, the forward rate of this reaction is automatically corrected, as it would be for the correct reaction, and the artefact is eliminated. We identify 45 such spurious reactions in the large DIANA-standard network, and 68 when the STAND reactions are added, demonstrating once again that this is a common issue in astrochemistry (Tinacci et al. 2023).

We identify the following reasons for the increase in N₂H⁺ abundance in various disk regions when reverse reactions are considered. First, the following constructed endothermic formation reactions become significant at very high temperatures (several thousand Kelvin): $$\mathrm{N}+\mathrm{N}\mathrm{H}\to {\mathrm{N}}_2{\mathrm{H}}^{+}+{\mathrm{e}}^{-}$$(24) $${\mathrm{N}}_2+\mathrm{H}\to {\mathrm{N}}_2{\mathrm{H}}^{+}+{\mathrm{e}}^{-}.$$(25)

Second, one of the dominant N₂H⁺ destruction reactions: $${\mathrm{N}}_2{\mathrm{H}}^{+}+\mathrm{O}\to {\mathrm{N}}_2+{\mathrm{O}\mathrm{H}}^{+}$$(26)

is found to be endothermic (+0.12 eV) but is reported as barrier free. This means that at temperatures as low as 20 K, the rate of this reaction is automatically corrected to very low values, leading to an increase in the N₂H⁺abundance.

The abundance structures of H₂O, CO, and NH₃ are similar between the two models. The abundance of c−C₃H₂ has three reservoirs. The innermost reservoir within 0.1 au decreases in abundance when the reverse reactions are taken into account, whereas the reservoir around 1 au increases. A similar trend is observed for CH₃ at 1 au. The inner reservoir of HCN is smoother in abundance when including the reverse reactions.

The abundance of C₂H₂ drops in the inner disk (≲0.1 au) between DIANA-standard and DIANA-standard-rev. The reactions mostly contributing to the formation of C₂H₂ in DIANA-standard are $${\mathrm{C}}_2{\mathrm{H}}_3+\mathrm{H}\to {\mathrm{C}}_2{\mathrm{H}}_2+{\mathrm{H}}_2$$(27) $${\mathrm{C}}_2\mathrm{H}+{\mathrm{H}}_2\to {\mathrm{C}}_2{\mathrm{H}}_2+\mathrm{H}.$$(28)

These reactions are also reported as dominant by Kanwar et al. (2024b) (grid point 3, fiducial model). However, when reverse reactions are included, each of these major formation reactions is balanced by the reverse of these reactions, thus resulting in no net contribution in the inner region (similar to detailed balance). Subsequently, the destruction of C₂H₂ by cosmic rays becomes dominant, leading to a decrease in abundance. These changes caused by the addition of reverse reactions can be crucial for the interpretation of JWST observations that can probe these deeper inner disk regions.

5.2 Effect of expanding to the CHAITEA network

When comparing the DIANA-standard-rev model to the CHAITEA model (the middle and right columns of Figs. D.1, D.2, D.3), we observe that many molecular abundances are affected in various ways. This outcome is expected, as we replaced the UMIST 2022 reaction data with the STAND reaction set, leading to significant modifications in most of the important chemical reactions. However, the overall impact is not substantial. For example, the abundances of important species such as H₂O, CO, and e⁻ remain largely similar between the two models. In light of the mid-infrared and (sub-)millimetre observations, we discuss below a few species that show large differences in their abundance structures between the DIANA-standard-rev and CHAITEA models in more detail.

Between the two models, there are significant differences in the inner midplane that is probed by JWST. For example, the abundance of CO₂ increases within 1 au when using the CHAITEA chemical network. This is due to more formation of CO₂ in CHAITEA relative to the DIANA-standard-rev model. The dominant formation reaction in CHAITEA model is $$\mathrm{C}\mathrm{O}+\mathrm{O}\mathrm{H}\to {\mathrm{C}\mathrm{O}}_2+\mathrm{H}.$$(29)

This reaction is barrier free according to the STAND network. However, Frost et al. (1993) and Millar et al. (2024) report an activation barrier of 176 K. We investigated whether incorporating this barrier would alter the CO₂ abundance structure, but find no significant changes. As a result, the final CHAITEA network retains the reaction as barrier free, consistent with the STAND network. This can be further revisited in future work. The second notable change is a 2 orders of magnitude decrease in C₂H₂ abundance in the inner and central midplane when using the CHAITEA network. Kanwar et al. (2024b) also report a reduction in C₂H₂ abundance in the inner midplane when larger hydrocarbons were included in an extended hydrocarbon chemical network. These larger hydrocarbons are not included in the current network and are beyond the scope of this study. Integrating the two networks could thus lead to substantially lower C₂H₂ abundances but eventually provide a more accurate abundance structure for mid-infrared observations (Tabone et al. 2023; Arabhavi et al. 2024; Xie et al. 2023; Kanwar et al. 2024a; Schwarz et al. 2024). The change in HCN abundance is more complex. In the DIANA-standard-rev model, the thin, higher-abundance surface layer extends farther radially than in the CHAITEA model (see Fig. D.2). There is less formation of HCN in the latter. The major formation pathway in the surface layer in both the models is $$\mathrm{C}\mathrm{N}+{\mathrm{H}}_2\to \mathrm{H}\mathrm{C}\mathrm{N}+\mathrm{H}.$$(30)

In the DIANA-standard-rev model, this reaction has a very low energy barrier (6.3 · 10⁻⁴ K). However, in the CHAITEA model, the reaction has a higher barrier of 2.3 · 10³ K, resulting in a decrease in HCN abundance in the surface layer by an order of magnitude. The abundance of some ions such as CH⁺increases in the surface layers of the outer disk in the CHAITEA model and this may affect the interpretation of the Herschel Space Telescope observations (Thi et al. 2011).

There are also significant differences in the outer disk that may affect the inferences from ALMA observations. For example, the abundance of HCN decreases in the CHAITEA model in the surface layers of the outer disk relative to DIANA-standardrev by at least 1.5 orders of magnitude. This could lead to lower line fluxes at millimetre wavelengths when using the CHAITEA chemical network. Other neutral molecules such as NH₃, C₂H, CH₃, and c−C₃H₂ also show substantially lower abundances in the outer disk when using the CHAITEA chemical network. C₂H plays a key role in determining the carbon budget in the outer disk. The abundance of the key molecular ion N₂H⁺ on the other hand increases in the outer disk when using the CHAITEA chemical network. Qi et al. (2015) and, van ’t Hoff et al. (2017) used N₂H⁺to infer the location of the CO iceline in the disk as it is difficult to observe directly. There is also a decrease in the abundances of CH₃ and c-C₃H₂ in the inner disk by 2−3 orders of magnitude in the CHAITEA model compared to DIANA-standard-rev. This can help to interpret JWST (Arabhavi et al. 2024; Kanwar et al. 2024a) and ALMA (Qi et al. 2013; Bergin et al. 2016) observations better.

6 Summary and conclusions

We merged the STAND chemical network developed for planetary and exoplanetary atmospheres (Rimmer & Helling 2016, 2019), which includes pressure-dependent termolecular and reverse reactions based on Gibbs free energy data, with the UMIST 2022 chemical network and our large DIANA standard chemical network developed for disks (Kamp et al. 2017), which was later updated in Kanwar et al. (2024b). A cornerstone in developing the new network was compiling a comprehensive set of Gibbs free energy data from various sources. For some molecules, data had to be estimated using ionisation potentials and proton affinities to ensure all gas-phase reactions could be reversed. This new, integrated chemical network is named CHAITEA.

We then studied under what conditions thermodynamic equilibrium can be established in planet-forming disks. Finally, we explored how the chemical composition of the disk changed when we used the new CHAITEA chemical network instead of the large DIANA standard network. We used a 2D T Tauri disk model for this investigation, including UV, X-ray, and cosmic-ray chemistry:

• In the absence of any UV, X-ray or cosmic-ray irradiation, we verify that the new CHAITEA network converges towards thermodynamic equilibrium by benchmarking it against the chemical equilibrium code GGchem (Woitke et al. 2018);

• However, at a gas density of 7.56 × 10¹³ cm⁻³ (∼10 dyn/cm²), the chemical relaxation timescale τ_{chem, 3}, towards thermodynamic equilibrium can be as short as a month at 2000 K, which increases to ∼10⁸ yr at 1200 K, and exceeds the age of the universe at 1100 K;

• The chemical relaxation proceeds simultaneously via fast and slow chemical modes with characteristic timescales that differ by many tens of orders of magnitude. We find that one of the slowest modes is the production of N₂ from other N-bearing molecules after the neutral atoms are consumed;

• The inclusion of radiative association, direct recombination reactions (considered non-reversible), and photo rates based on a local Planck radiation field, J_v = B_v(T), causes deviations from thermodynamic equilibrium by 2−20% across all densities and temperatures studied;

• The inclusion of cosmic rays leads to significant deviations from thermodynamic equilibrium, initially affecting the charged molecules and eventually impacting the abundant neutral species. At low temperatures (T ≲ 750 K), even a tiny CRI of 10⁻⁴⁰ s⁻¹ alters the neutral chemistry. At the standard CRI of 1.7 × 10⁻¹⁷ s⁻¹, below ∼1200 K, deviations from thermodynamic equilibrium can exceed 10−30 orders of magnitude;

• In a 2D disk model, no region is found to be in thermodynamic equilibrium. We find that a small, warm, dense region in the midplane directly behind the inner rim (r ≲ 0.1 au), where the cosmic rays are mostly shielded, is closest to thermodynamic equilibrium, with σ ∼ 1. There is a warm intermittent molecular layer where σ ∼5. In all other regions σ is much larger.

• In kinetic equilibrium, nitrogen is primarily stored in N₂, carbon in CO and CO₂, and oxygen in H₂ O and CO, whereas the thermodynamically favoured molecules at low temperatures are NH₃ for nitrogen, CH₄ for carbon, and H₂ O for oxygen;

• When using the CHAITEA chemical network, the change in the abundances of species such as CO, H₂O and free electrons is not major (an order of magnitude). The abundance of species such as C₂H₂, HCN and CH₄ decreases in the inner midplane, while the abundance of CO₂ increases within 1 au. Some molecular cations such as N₂H⁺and CH⁺ attain larger abundances in the intermittent warm molecular layer and in the outer disk regions. The reasons for these differences are complex. Differences in the rate coefficients of the leading reactions can have profound effects on the resulting abundances. Defining reaction pairs based on Gibbs free energy allows for the identification and automatic correction of spurious endothermic reactions reported with zero or minimal activation energies, which are chemically implausible. Such reactions can significantly impact the chemistry of the outer disk.

The techniques implemented in this work, and the new chemical network developed, can be used to interpret both the (sub-) millimetre molecular lines observable with ALMA and the midinfrared observations with JWST. The latter probes the inner, dense, and warm regions where this network is more complete. The CHAITEA chemical network can also be used to quantify the significance of the termolecular and reverse reactions using the observations.

Acknowledgements

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie SklodowskaCurie grant agreement No. 860470. We thank the anonymous referee for their constructive comments. We acknowledge discussions with W.F. Thi. J.K. acknowledges productive discussions with Till Kaeufer, Nidhi Bangera, David Lewis, and Ashwani Rajan.

Appendix A Derivation of the law of mass action

We assume a mixture of ideal gases with partial pressures p_i = N_iRT/V, where N_i are the mol numbers of the gas species i, V is the common volume, and R [J/K/mol] is the ideal gas constant. P = ∑_ip_i is the total gas pressure. (T^⊖, p^⊖) denotes a thermodynamical reference state at room temperature T^⊖ = 298 K and at standard pressure p^⊖ = 1 bar, to which some of the thermodynamic quantities refer to. The thermodynamic description of such mixtures, see Atkin & de Paula (2006), is $$G=U+PV-TS$$(A.1) $$U=\sum _i{N}_i(u_i^{\ominus }+{\int }_{T^{\ominus }}^T{c}_{V,i}dT)$$(A.2) $$PV=\sum _i{N}_{i}RT$$(A.3) $$S=\sum _i{N}_i(\underset{s_i(T,p^{\ominus })}{\underbrace{s_i^{\ominus }+{\int }_{T^{\ominus },p^{\ominus }}^{T,p^{\ominus }}(c_{V,i}+R)\frac{dT}{T}}}-\underset{R\ln \frac{p_i}{p^{\ominus }}}{\underbrace{{\int }_{T,p^{\ominus }}^{T,p_i}R\frac{dP}{P}}}),$$(A.4)

where G is the system Gibbs free energy, U the internal energy, and S the entropy. The small symbols refer to the molar properties of the pure gases $i:u_i^{\ominus }$ is the internal energy in the standard state (T^⊖, p^⊖), c_{V, i} the heat capacity, and $s_i^{\ominus }$ is the entropy in the standard state (T^⊖, p^⊖). s_i(T, p^⊖) is the entropy of the pure gas i at temperature T and standard pressure p^⊖. The terms in the integrals in Eq. (A.4) can be derived from the first law of thermodynamics in the form dU = −PdV + TdS for processes of an ideal gas at constant temperature or pressure. The term R ln (p_i/p^⊖) is the mixing entropy. Altogether, the system Gibbs free energy is $$G=\sum _i{N}_i(\underset{g_i(T,p^{\ominus })}{\underbrace{u_i^{\ominus }+{\int }_{T^{\ominus }}^T{c}_{V,i}dT+RT-T{s}_i(T,p^{\ominus })}}+RT\ln \frac{p_i}{p^{\ominus }}),$$

where g_i(T, p^⊖) is the molar Gibbs free energy of the pure gas i at temperature T and standard pressure p^⊖.

Next, we consider a hypothetical chemical reaction N₂ + 3 H₂ → 2 NH₃, where k = 3, {1 → N₂, 2 → H₂, 3 → NH₃} and stoichiometric coefficients {v₁ = −1, v₂ = −3, v₃ = + 2}. We also consider a small number of mols of such reactions Δ to take place at constant temperature and pressure. The principle of minimisation of the system Gibbs free energy in this case requires that $$\mathrm{\Delta }G=G(T,{\stackrel{\circ }{N}}_1+v_1\mathrm{\Delta },\dots ,{\stackrel{\circ }{N}}_k+v_k\mathrm{\Delta })-G(T,{\stackrel{\circ }{N}}_1,\dots ,{\stackrel{\circ }{N}}_k)=0$$(A.5)

in the thermodynamical equilibrium state that is characterised by the mol numbers ${\stackrel{\circ }{n}}_i$. From Eq. (A.5) it follows that $$\mathrm{\Delta }G=\mathrm{\Delta }\sum _i{v}_i(g_i(T,p^{\ominus })+RT\ln \frac{{\stackrel{\circ }{p}}_i}{p^{\ominus }})=0$$(A.6)

for any small Δ. The Gibbs free energy data for the pure gases i can be found in various thermochemical databases, here we use Burcat & Ruscic (2005). These data are given as the Gibbs free energy of formation ${\mathrm{\Delta }}_{\mathrm{f}}{G}_i^{\ominus }$(T)[J/mol] of a molecule from the standard states of the elements at room temperature. Since any chemical reaction obeys the principle of element conservation, these constants cancel in Eq. (A.6), so we can write this equation also as $$\sum _i{v}_i{\mathrm{\Delta }}_{\mathrm{f}}{G}_i^{\ominus }=-RT\sum _i{v}_i\ln \frac{\stackrel{\circ }{p_i}}{p^{\ominus }}.$$(A.7)

Introducing the reaction Gibbs free energy as $${\mathrm{\Delta }}_{\mathrm{R}}{G}^{\ominus }=\sum _{i=1}^k{v}_i{\mathrm{\Delta }}_{\mathrm{f}}{G}_i^{\ominus }$$(A.8)

we find the law of mass action $$\prod _{i=1}^k{\left(\frac{\stackrel{\circ }{p_i}}{p^{\ominus }}\right)}^{v_i}=\exp (-\frac{{\mathrm{\Delta }}_{\mathrm{R}}{G}^{\ominus }}{RT}).$$(A.9)

Equation (A.9) is valid for all hypothetical chemical reactions in thermodynamical equilibrium, even if they cannot proceed kinetically. We note that ΔG = 0 in thermodynamical equilibrium, whereas Δ_RG^⊖ ≠ 0 because of the mixing entropy. In the main text in Sect. 2.2 we consider the reaction A + B → C + D + E, where k = 5 and {v₁ = −1, v₂ = −1, v₃ = + 1, v₄ = + 1, v₅ = + 1}.

Appendix B Thermochemical data of missing species in Burcat format

Table B.1 summarises our efforts to collect, convert, estimate or construct the thermochemical in Burcat format (NASA 7-polynomials) of those species missing in Burcat & Ruscic (2005).

Appendix C Disk model parameters

We provide the list of input parameters used to model the disk in Table C.1.

Appendix D Abundances of various species

Figures D.1, D.2 and D.3 show the abundances of various species discussed in Sect. 5. These abundances are calculated using different chemical networks and rate databases that are outlined in Table 8.

Fig. D.1 The abundances of various species in the three models calculated using different chemical networks and rate databases outlined in Table 8. The white dashed contours correspond to the tick values on the colourbar.

Fig. D.2 Continued from Fig. D.1

Fig. D.3 Continued from Fig. D.2

References

All Tables

All Figures

	Fig. 1 Abundances of various species, ϵ(X)=nX/n〈H〉, obtained with the time-independent PRODIMO models from grid_0 for density n〈H〉= 7.56 × 1013 cm−3 (dots), compared to the corresponding GGchem models (dashed lines). The deviations at the lowest temperatures are due to numerical problems in converging the PRODIMO models.
In the text

	Fig. 2 Abundances, ϵ(X)=nX/n〈H〉, of various species (X) as a function of time obtained with Pro DI Mo at 1200 K and 2000 K at a density of n〈H〉 = 7.56 × 1013 cm−3. The large coloured dots represent the abundances in thermodynamic equilibrium obtained with GGchem.
In the text

Fig. 3 Deviation σ (see Eq. (19)), between particle densities obtained with time-independent PRODIMO (ni) and GGchem $({\stackrel{\circ }{n}}_i)$ at various temperatures and densities from grid_0. The blank rectangles at the bottom mark the models which do not converge. The numbers on the x-axis correspond to the ascending list of densities in Table 5. The central white labels mark the species that contribute most to σ.

In the text

	Fig. 4 Same as Fig. 3 but from grid_photo which includes photo processes based on a local Planck field Jv = Bv(T), radiative association, and direct recombination reactions.
In the text

	Fig. 5 Impact of cosmic-ray ionisation rate on deviations from thermodynamic equilibrium at different densities and temperatures from grid_cr. The adopted values of CRIs are listed at the top of each panel. The x- and y-axes are the same as in Fig. 3.
In the text

	Fig. 6 Abundances of selected species in the 2D disk model calculated using chemical kinetics (left) and thermodynamic equilibrium (right). The white contours correspond to the tick values on the colorbars.
In the text

	Fig. 7 Continued from Fig. 6. We note the different scaling of the GGchem results for HCN and C2H2.
In the text

	Fig. 8 Deviation (σ) between kinetic chemical equilibrium and thermodynamic equilibrium in the 2D disk model. σ = 1 denotes a deviation of one order of magnitude between the PRODIMO and GGchem abundances.
In the text

	Fig. D.1 The abundances of various species in the three models calculated using different chemical networks and rate databases outlined in Table 8. The white dashed contours correspond to the tick values on the colourbar.
In the text

	Fig. D.2 Continued from Fig. D.1
In the text

	Fig. D.3 Continued from Fig. D.2
In the text