© The Authors 2025

1 Introduction

The formation of polycyclic aromatic hydrocarbons (PAHs) in astrophysically favorable environments – particularly in circumstellar envelopes of asymptotic giant branch (AGB) stars – has been extensively studied over several decades (Frenklach & Feigelson 1989; Cherchneff et al. 1992; Cadwell et al. 1994; Cherchneff 2012). Nevertheless, a complete understanding of the evolution of aromatic hydrocarbons remains elusive.

Acetylene (C₂H₂) has been identified as a key molecule in PAH and carbon dust formation in AGB stars (Gail & Sedlmayr 1987; Cherchneff 2012; Dhanoa & Rawlings 2014; Pentsak et al. 2024). The pioneering kinetic model by Frenklach & Feigelson (1989), based on combustion chemistry networks, introduced the hydrogen abstraction–acetylene addition (HACA) mechanism as a primary pathway for PAH growth. While this work established chemical routes for PAH formation under circumstellar conditions, it also revealed that the process is efficient only within a narrow gas temperature range (≈900–1000 K) and in stars with high mass-loss rates. These findings prompted further exploration of alternative PAH formation pathways and detailed reaction analysis.

Alternative mechanisms, such as the recombination of two propargyl radicals proposed by Miller et al. (1990) and Cherchneff et al. (1992), proved more promising in circumstellar environments. However, stationary-flow models alone could not account for the observed PAH abundance (Cherchneff et al. 1992; Helling et al. 1996). By coupling astrochemical kinetical model with envelope dynamics – including shock waves from stellar pulsations – Willacy & Cherchneff (1998); Cherchneff (2012) significantly enhanced PAH production efficiency. Their simplified estimates yielded a carbonaceous dust-to-gas mass ratio of ≈ 1.2 · 10⁻³−5.8 · 10⁻³, consistent with observational and galactic dust evolution values of 10⁻⁴−10⁻³ (Ferrarotti & Gail 2006; Zhukovska et al. 2008).

While most studies focus on gas-phase PAH formation, concurrent dust nucleation and growth processes are equally critical. Carbonaceous dust nucleation in AGB stars may occur via homogeneous or heterogeneous pathways (Gail & Sedlmayr 1988; Cadwell et al. 1994). Cadwell et al. (1994) proposed induced nucleation, with silicon carbide (SiC) particles serving as nucleation seeds due to their higher formation temperatures (Frenklach et al. 1989; Lorenz-Martins & Lefevre 1994; Kozasa et al. 1996). These particles, observed in AGB stars and meteorites (Speck et al. 2009), provide substrates for PAH and other hydrocarbons adsorption, making heterogeneous nucleation far more effective than homogeneous processes. Notably, Cadwell et al. (1994) hypothesized that reactions on the SiC surface could explain the high PAH abundance in the interstellar medium (ISM). The PAHs formed on grains might later desorb under ultraviolet (UV) radiation during post-AGB or planetary nebula phases.

Apart from providing nucleation sites, SiC surfaces catalyze aromatic formation. Studies confirm efficient acetylene chemisorption on SiC active sites (Frenklach et al. 1989; Carmer et al. 1993; Cho et al. 2001; Navarro-Ruiz et al. 2013; Pollmann et al. 2014), followed by cyclotrimerization (CTM) to benzene (C₆H₆) – a process experimentally validated by Zhao et al. (2016). Subsequent surface reactions could coat SiC particles with graphene-like layers. In its turn, Merino et al. (2014) suggested that SiC surfaces first graphitize, then fragment into PAHs via hydrogen passivation or etching. Both Zhao et al. (2016) and Merino et al. (2014) underscore SiC importance in PAH formation. Carbonaceous materials also facilitate acetylene CTM, a key carbocatalytic process (Pentsak et al. 2020; Pentsak et al. 2024). First discovered by Berthelot (1867) and then independently by Zelinsky (1923), this reaction gained renewed interest through density functional theory (DFT) studies by Gordeev et al. (2020), who demonstrated its viability on zigzag PAH edges with low activation barriers. Gordeev et al. (2020) studied several PAHs containing from 6 to 361 C atoms and concluded that most of them (except phenyl radical) can be considered as catalysts for acetylene CTM which demonstrates the versatility of this mechanism.

In this work, we quantify the impact of acetylene CTM on SiC and PAH zigzag edges in AGB envelopes. Combining the dynamical model of IRC+10216 from Cherchneff (2012) with an astrochemical model – incorporating gas-phase, surface reactions, and dust growth (adsorption and coagulation) – we assess PAH abundance enhancement under realistic circumstellar conditions.

2 The astrochemical model

2.1 Main statements

The model used in this work builds upon our previous study, Murga (2024, hereafter Paper I), in which we incorporated acetylene CTM on SiC particle surfaces into the chemical reaction network. The key assumptions of the model in Paper I are as follows:

• The physical conditions for astrochemical calculations – gas number density (n_gas) and temperature (T_gas) – are adopted from Cherchneff (2012); Willacy & Cherchneff (1998). These conditions were derived from a dynamical model of the envelope of the star IRC+10216, accounting for periodic shock wave propagation related to stellar thermal pulses. Such propagation results in a sharp increase in the gas density and temperature. Cooling of disturbed gas occurs in two steps due to: (1) dissociation of H₂ molecules; (2) adiabatic expansion. After cooling, gas moves as a stationary flow before being disturbed by another shock wave.

• SiC particles grow via SiC molecule adsorption, while coagulation is neglected. The initial number of particle seeds is set to 10⁻¹² relative to the number of H₂ molecules.

• The chemical reaction network is based on that from Cherchneff (2012). The original network consists of 377 gas-phase reactions, with benzene as the largest molecule. In Paper I, this network has been extended with surface reactions involving acetylene adsorption on SiC particles and subsequent CTM into benzene.

In the present work, we present the updated model (designated as BRAHMA), which further enhances the previous one by incorporating the following processes:

• Chemisorption of acetylene molecules on PAH zigzag edges and their CTM;

• Aromatic growth pathways from benzene to pyrene;

• Dust grain growth via hydrocarbon adsorption;

• Dust grain coagulation.

We exclude homogeneous nucleation of SiC particles, as the results of Paper I suggest that SiC dust forms predominantly close to the star before aromatic molecules and carbonaceous particles emerge. Subsequent dust growth is assumed to occur via hydrocarbon adsorption onto SiC particles, consistent with the induced nucleation mechanism studied in Frenklach et al. (1989); Cadwell et al. (1994). Homogeneous nucleation and growth of purely carbonaceous particles are omitted, as test calculations indicate their inefficiency compared to heterogeneous nucleation, while their inclusion significantly slows down computations.

Our reaction network is limited to pyrene, justified by experimental and theoretical considerations. Studies on the formation of interstellar dust analogs show that small PAHs are predominantly present in gas during the growth process (Santoro et al. 2020). This is partially due to observational constraints but it also reflects the lower abundance of larger PAHs. Similarly, combustion theories suggest that moderate-sized PAHs dominate soot formation (Wang 2011; Jacobson et al. 2020). Thus, we argue that restricting PAH size to pyrene is well founded. Intriguingly, PAHs in the experiment on the formation of dust analogs (Santoro et al. 2020) often bear aliphatic substituents (e.g., – CH₂, – C₂H₃), which may bridge aromatic units, yielding hydrogen-rich, irregular structures. Such PAHs are present in our network and also participate in the dust growth process.

2.2 Kinetic equations

The astrochemical model presented above is described by the system of discrete kinetic equations, which include separate equations for the number density of species in gas (n_i), on the dust surface $\left(n_i^{\text{s}}\right)$, and dust grains themselves $\left(n_i^{\text{d}}\right)$. Kinetic equations for the number densities normalized to n_gas are $$\begin{array}{lll}\frac{d{n}_i}{dt}\hfill & =\hfill & {\displaystyle \sum _j{k}_j{n}_j+{\displaystyle \sum _{l,m}k{}_{lm}{n}_j{n}_m{n}_{gas}+{\displaystyle \sum _{o,p,t}{k}_{opt}{n}_o{n}_p{n}_t{n}_{gas}^2}}}\hfill \\ \hfill & \hfill & -n_i\left|{\displaystyle \sum _s{k}_s+{\displaystyle \sum _q{k}_{iq}{n}_q{n}_{gas}+{\displaystyle \sum _{v,w}{k}_{ivw}{n}_v{n}_w{n}_{gas}^2}}}\right|\hfill \\ \hfill & \hfill & +{\displaystyle \sum _{z=1}^{N_a}{k}_{iz}^{\text{des}}{n}_{iz}^s-{\displaystyle \sum _{z=1}^{Na}{k}_{iz}^{\text{ad}}{n}_i},}\hfill \end{array}$$(1) $$\frac{d{n}_{iz}^{\text{s}}}{dt}=k_{iz}^{\text{ad}}{n}_i-k_{iz}^{\text{des}}{n}_{iz}^{\text{s}}+k_{ijz}^{\text{sr}}{n}_{jz}^{\text{s}}{n}_{{\text{C}}_2{\text{H}}_2}{n}_{\text{gas}}-k_{liz}^{\text{sr}}{n}_{iz}^{\text{s}}{n}_{{\text{C}}_2{\text{H}}_2}{n}_{\text{gas}},$$(2) $$\frac{d{n}_{iz}^{\text{d}}}{dt}=\frac{1}{2}{\displaystyle \sum _{j=1}^{N_a}{\displaystyle \sum _{l=1}^{N_a}{C}_{ili}{\beta }_{jl}{n}_j^{\text{d}}{n}_l^{\text{d}}{n}_{\text{gas}}}-n_i^{\text{d}}}{\displaystyle \sum _{j=1}^{N_a}{\beta }_{ij}{n}_j^{\text{d}}{n}_{\text{gas}}}$$(3) where k_j and k_s, k_lm and k_iq, k_opt and k_ivw are rate coefficients of unimolecular, bimolecular and termolecular reactions, correspondingly, $k_{iz}^{\text{des\hspace{0.17em}}}$ and $k_{iz}^{\text{ad\hspace{0.17em}}}$ are desorption and adsorption rates, $k_{ijz}^{\text{sr}}$ (and $k_{liz}^{\text{sr}}$) is a rate coefficient of a surface reaction between C₂H₂ and j(i)-th molecule on the z-th grain with formation of the i(l)-th molecule. The second index in $k_{ijz}^{\text{sr}}\left(j\right)$ can correspond to sC₂H₂¹ and sC₄H₄, the first index in $k_{\text{liz\hspace{0.17em}}}^{\text{sr}}\left(l\right)$ – to sC₄H₄ or sC₆H₆. C_jli is the value related to the mass (m[index]) that was put to the i-th bin after coagulation of j-th and l-th grains, which we calculated according to Akimkin (2015) as $$C_{jli}=\{\begin{array}{ll}ϵ,\hfill & \text{\hspace{0.17em}if\hspace{0.17em}}{m}_i=\text{max}\left\{m_n\right\}\text{\hspace{0.17em}where\hspace{0.17em}}{m}_n<m_j+m_l,\hfill \\ 1-ϵ,\hfill & \text{\hspace{0.17em}if\hspace{0.17em}}{m}_i=\text{min}\left\{m_n\right\}\text{\hspace{0.17em}where\hspace{0.17em}}{m}_n>m_j+m_l,\hfill \\ 0,\hfill & \text{\hspace{0.17em}other\hspace{0.17em}cases}.\hfill \end{array}$$(4)

In Eq. (3), N_a is the number of dust grain populations, the factor β_jl is the frequency of collisions between j-th and l-th grains, which is calculated as $${\beta }_{jl}=2.2\sqrt{\frac{\pi k_{\text{B}}{T}_{\text{gas}}}{2{\mu }_{jl}}}{\left(2a_j+2a_l\right)}^2,$$(5) where k_B is the Boltzmann constant, µ_jl is the reduced mass of j-th and l-th grains, aj and a_l are their radii, and factor 2.2 is the Van der Waals enhancement factor (Harris & Kennedy 1988).

At each calculation step we performed two corrections. The first one was to take into account the change of the grain mass and radius due to surface reactions. For the i-th grain, we calculated the difference, dm_i, between masses of all surface components at the current and previous steps and added this difference to the mass of the i-th grain. The radius of the i-th grain, a_i, was recalculated from the new mass.

The second correction was to redistribute surface components between grains according to changes in dust grain number densities after coagulation process for the time interval dt. To find the corrected number density of i-th surface component on the z-th dust grain ${\left(n_{iz}^{\text{s}}\right)}^{\text{corr}}$, we modified Eq. (3). Instead of the dust grain number density, the equation was rewritten for the number density of surface components. The first term in the right-hand side of Eq. (3) was multiplied by the sum of number densities of surface components of each grains normalized to the number densities of these grains. The number density of the z-th dust grains was replaced by the uncorrected number density of the i-th component on the z-th dust grain in the second term in the right part of Eq. (3). Since we calculated the number density instead of its change, the right part obtained was multiplied by dt. The expression for ${\left(n_{iz}^{\text{s}}\right)}^{\text{corr}}$ is $${\left(n_{iz}^{\text{s}}\right)}^{\text{corr}}=dt\times \left(\frac{1}{2}{\displaystyle \sum _{j=1}^{N_a}{\displaystyle \sum _{l=1}^{N_a}{C}_{jli}{\beta }_{jl}}\left(\frac{n_{ij}^{\text{s}}}{n_j^{\text{d}}}+\frac{n_{ij}^{\text{s}}}{n_l^{\text{d}}}\right)n_j^{\text{d}}{n}_l^{\text{d}}{n}_{\text{gas}}-n_{iz}^s{\displaystyle \sum _{j=1}^{N_a}{\beta }_{jz}{n}_j^{\text{d}}{n}_{\text{gas}}}}\right).$$(6)

2.3 Rate coefficients

The rate coefficients for monomolecular, bimolecular, and termolecular reactions were calculated using the modified Arrhenius law with kinetic parameters A, n and E_a (the activation barrier), which were defined for each reaction individually (see Sect. 2.4). The adsorption rate coefficient was calculated as $$k_{iz}^{\text{ad}}=4\pi a_z^2{v}_i^{\text{th}}{\alpha }_i{n}_z^{\text{d}},$$(7) where $v_i^{\text{th}}$ is the thermal velocity of the i-th molecule and α_i is its sticking coefficient. We assume that all hydrocarbons contribute to dust grain growth proportionally to their gas-phase abundances. The most abundant species, such as acetylene (C₂H₂), play the dominant role. We adopted α_i = 1 for all molecules. We accounted for adsorption onto sites occupied by sC₂H₂ or sC₄H₄ separately, as they are intermediates in CTM (see below). Therefore, the rate coefficient for C₂H₂ adsorption was multiplied by f , the fraction of sites not occupied by these species. We consider the adsorption process to be a direct chemisorption with no activation barrier.

The desorption rate coefficient is $$k_{iz}^{\text{des}}=v_0\text{exp}\left(-\frac{E_{\text{d}}^i}{k_{\text{B}}{T}_z^{\text{d}}}\right),$$(8) where ν₀ is a characteristic frequency (see Hasegawa et al. 1992), $E_{\text{d}}^i$ is the desorption energy of the i-th molecule, and $T_z^{\text{d}}$ is the temperature of the z-th dust grain.

The value of $E_{\text{d}}^i$ is a cornerstone in our astrochemical model, so we considered its choice in more detail. Hydrocarbons first adsorb onto SiC grain surfaces. The SiC surface is chemically active for forming strong bonds with hydrocarbons. Calculations presented in Wieferink et al. (2006); Pollmann et al. (2014) provide values of 1–5 eV depending on the molecular species and the considered reconstruction of SiC surface. Such values strongly support the chemical nature of bonding between the acetylene molecule and the SiC surface. For such a scenario, the desorption process proceeds less effectively, as more energy is required for desorption of adsorbed acetylene, being in a stable configuration.

After the formation of a carbonaceous layer on SiC grains, subsequent adsorption occurs on this carbonaceous surface. While covalent bonding can occur at specific sites, adsorption is not universally chemisorptive. Bond strength critically depends on local surface structure and defect presence (e.g., radicals, zigzag edges). Whereas basal graphene planes exhibit low reactivity, PAH zigzag edges and carbon radicals are highly active sites (Radovic & Bockrath 2005). Krueger et al. (1996) argued that adsorption of C₂H₂ results in the formation of a sp³ bond, which finally leads to grains having an amorphous structure. Thus, circumstellar environments likely produce irregular carbon surfaces with numerous potential binding sites.

Determining a single desorption energy for all hydrocarbons and surfaces is challenging due to multiple crucial factors influencing it. In the absence of direct measurements, we can rely only on indirect estimations. Among them, one can mention the work of Alamdari et al. (2023), who studied ethylene adsorption on graphenes and gave $E_{\text{d}}^{{\text{C}}_2{\text{H}}_4}$ in the wide range 1–5 eV. Based on the results of Gordeev et al. (2020) who investigated the kinetics of acetylene CTM on carbene active sites, one may conclude that the dissociation of the bond between acetylene and carbene site needs around 2 eV.

From the works on laboratory formation of interstellar dust analogs, we can conclude that the structure of the dust grains corresponds to hydrogenated amorphous carbon with a mixture of sp, sp², and sp³ bonds. The strength of these bonds also varies in a wide range, ∼2–5 eV (Robertson 2002; Sugiura et al. 2020), depending on the sp²/sp³ ratio. Additionally, some estimations of binding between aromatic molecules can be borrowed from works on the soot inception. It is believed that the covalent coalescence takes place during the initial soot inception (Gentile et al. 2020). Pure van der Waals interaction would yield unstable particles. Frenklach & Mebel (2020) identify E-bridge formation (two aliphatic C–C bonds) as the optimal pyrene clustering pathway, with DFT calculations showing a bond strength of ∼2 eV.

Given the uncertainties discussed above, we performed calculations for E_d = {2, 3, 4} eV, adopting the intermediate value E_d = 3 eV as a standard one². The range 2–4 eV encompasses energies expected for desorption of surface species from SiC or carbonaceous dust grains.

We included two surface reactions of Eley-Rideal type – sC₂H₂+C₂H₂→sC₄H₄ and sC₄H₄+C₂H₂→sC₆H₆. Herewith we adopted that the formed sC₆H₆ desorbs without a barrier immediately after cyclization, which seems to be plausible based on our DFT calculations of acetylene CTM reaction on the SiC surface (Fedotova et al., in prep.). Should we consider the CTM process in more detail, then the formation of benzene from three acetylene molecules can proceed via at least two distinct reaction channels. The first involves the linear growth of an acetylene chain to form sC₆H₆, followed by cyclization. The second pathway proceeds through the initial formation and surface stabilization of sC₄H₄ ring, which a subsequent acetylene molecule then opens and extends to form cyclic sC₆H₆. In the present model, we do not computationally distinguish these two channels; their kinetics are treated as identical. A detailed discussion of the differences will be provided in a forthcoming paper by Fedotova et al. (in prep.). It should be mentioned that we assume that after forming the carbonaceous layer on SiC dust grains surface reactions proceed further with the same formalism, meaning that carbon dangling bonds act as active sites instead of Si in this case.

In the Eley–Rideal mechanism, a reaction occurs directly between a gas-phase species and an adsorbed species on a surface. The rate coefficient for such a reaction is analogous to the adsorption rate coefficient, but must be computed for a single surface site. This is achieved by dividing the adsorption rate coefficient by the total number of sites available on a grain of a specific size, $n_z^{\text{d}}{n}_{\text{sites},z}$. Thus, the rate coefficient for the Eley–Rideal reaction between the j-th adsorbed molecule and a gas-phase acetylene molecule is given by $$\begin{array}{l}{k}_{ijz}=4\pi a_z^2{v}_{{\text{C}}_2{\text{H}}_2}^{\text{th}}{\alpha }_{{\text{C}}_2{\text{H}}_2}{n}_z^{\text{d}}\frac{1}{n_z^{\text{d}}{n}_{\text{sites},z}}\hfill \\ \text{or}\hfill \\ k_{ijz}=4\pi a_z^2{v}_{{\text{C}}_2{\text{H}}_2}^{\text{th}}{\alpha }_{{\text{C}}_2{\text{H}}_2}\frac{1}{n_{\text{sites},z}}.\hfill \end{array}$$(9)

As was mentioned, apart from the acetylene CTM on the SiC surface, we also considered this reaction on zigzag edges of PAHs. According to Gordeev et al. (2020), PAHs with carbene centers can be efficient catalysts for CTM of acetylene. The approximate value of E_a is ≈0.5 eV of the first and subsequent acetylene attachments, and its free energy is ∆G ≈ −1 eV. Cyclization proceeds with E_a ≈ 0.2 eV and ∆G ≈ −2 eV. After cyclization, the formed benzene molecule can be desorbed either directly or after another acetylene becomes associated with the same site. For small PAHs, the first way is more attractive because of lower values of E_a, which are ∼0.8–0.9 eV (C₆, C₁₄) against ≈1.5 eV for C₃₇. In its turn, the second way is more favourable for large PAHs as it requires energy around 1.2 eV almost independently of the PAH size.

We adapted the results of the work of Gordeev et al. (2020) in the following way: 1) only the largest PAH in our network, pyrene (radical), acts as a catalyst; 2) we added four reactions for pyrene: ${\text{C}}_{16}{\text{H}}_9^{•}+{\text{C}}_2{\text{H}}_2\to {\text{C}}_{16}{\text{H}}_9-{\text{C}}_2{\text{H}}_2$, C₁₆H₉-C₂H₂ + C₂H₂ → C₁₆H₉-C₄H₄, C₁₆H₉-C₄H₄ + C₂H₂ → C₁₆H₉-C₆H₆, and ${\text{C}}_{16}{\text{H}}_9-{\text{C}}_6{\text{H}}_6\to {\text{C}}_{16}{\text{H}}_9^{•}+{\text{C}}_6{\text{H}}_6$ ; 3) E_a for acetylene addition was adopted to be 0.5 eV; 4) the energy of C₂H₂- and C₄H₄-loss, $E_{\text{d}}^{\text{PAH}}$, is 3 eV; 5) the energy of C₆H₆-loss is 1 eV. The expression for rate coefficients for the reactions of the formation of the intermediates C₁₆H₁₀-C₂H₂, C₁₆H₁₀-C₄H₄, and C₁₆H₁₀-C₆H₆ is analogous to Eq. (9), except that we multiplied the adsorption rate by a factor of exp(−E_a/(k_BT^PAH)), accounting for the activation barrier of the attachment. To find the rate of the reaction, the rate coefficient was multiplied by the number densities of the corresponding PAH and gaseous acetylene.

2.4 The reaction network

There are several networks presented in the literature that are used to model the chemistry of AGB stellar envelopes. Frenklach & Feigelson (1989) in their pioneering work presented the network of reactions, which was mainly borrowed from combustion studies. This network and the calculations performed on its base become fundamental in modeling of aromatics in stellar envelopes. Later, Cherchneff et al. (1992) corrected the network, supplementing it with new studied reactions and neglecting the reactions that have little importance. Many further studies have relied on the results obtained with these networks (e.g. Allain et al. 1997; Pascoli & Polleux 2000; Cherchneff 2012), so the significance of the networks is undeniable. However, kinetic parameters for many reactions in these networks were poorly estimated or even extrapolated; namely, for the reactions of PAH growth via the HACA mechanism.

The most recent network of reactions was presented in Cherchneff (2012), but as far as aromatics are concerned it includes only reactions leading to formation of benzene. Thus, we took this network as a base for a new network and supplemented it with a number of reactions of PAH growth up to pyrene for which kinetic parameters are either available or could be estimated. We also revised the Cherchneff’s network and corrected kinetic parameters for some reactions if more accurate data are present in literature.

In our network, the first aromatic ring can be formed via following routes in the gas phase:

1. Recombination of propargyl radical resulting in the formation of benzene:

   C₃H₃^• + C₃H₃^• ⇌ a-C₆H₆³,

   for which we used the kinetic parameters given in Zhao et al. (2021), and in the formation of a phenyl radical:

   C₃H₃^•+ C₃H₃^• ⇌ a-C₆H₅^•+H, with kinetic parameters from Cherchneff (2012). The phenyl radical further forms benzene via collisions with H, H₂, or H-bearing species.

2. Isomerization of linear C₆H₆ (2-ethynyl-1,3-butadiene) to benzene directly:

   C₆H₆ ⇌ a-C₆H₆,

   or to a phenyl radical:

   C₆H₆ ⇌ a-C₆H₅^•+H.

3. The reactions of 1,3-butadien-1-yl (denoted as n-C₄H₅) and 1,3-butadien-2-yl (denoted as i-C₄H₅) with acetylene:

   n-C₄H₅ + C₂H₂ ⇌ a-C₆H₆ + H (Wang & Frenklach 1997)

   and

   i-C₄H₅ + C₂H₂ ⇌ a-C₆H₆ + H (Huang et al. 2018).

4. The reaction between 1,3-hexadien (denoted as l-C₆H₆) and H, which gives either a cycled C₆H₇ with additional hydrogen, which can dissociate later,

   l-C₆H₆+H ⇌ a-C₆H₇ (Wang & Frenklach 1997),

   or linear 1,3,5-hexatrien-1-yl, which can later cyclysize:

   l-C₆H₆+H ⇌ C₆H₇ (Wang & Frenklach 1997),

   C₆H₇ ⇌ a-C₆H₇ (Duran et al. 1988).

5. Isomerization of fulvene to benzene:

   a-C₅H₄-CH₂⁴ ⇌ a-C₆H₆ (Zhao et al. 2021).

   In its turn, fulvene can be formed either from phenyl radical (Zhao et al. 2021) or in the reaction between i-C₄H₅ and C₂H₂ (Huang et al. 2018).

6. The reaction between 1-buten-3-yn-lyl (denoted as n-C₄H₃) and acetylene results in formation of a phenyl radical, n-C₄H₃ + C₂H₂ ⇌ a-C₆H₅^• (Wang & Frenklach 1997), or benzyne, n-C₄H₃ + C₂H₂ ⇌ a-C₆H₄ +H (Wang & Frenklach 1997), which further may form benzene in reactions with H-bearing species as well as phenyl radical.

We consider two molecules, consisting of two aromatic rings – biphenyl and naphthalene. Biphenyl forms out of benzene and phenyl molecules and is then converted to phenanthrene. Naphthalene is formed via the HACA mechanism through subsequent steps of hydrogen abstraction and acetylene addition. On the way, two intermediates can appear – a-C₆H₅-C₂H₂ and a-C₆H₄-C₂H₃. The kinetic parameters were taken from Mebel et al. (2017). The path to naphthalene synthesis is a dead end in our network, as it stops after the formation of acenaphthylene. This limit of the HACA mechanism was found in a number of works (Kislov et al. 2005; Parker et al. 2015), although it has already been shown that naphthalene can grow further via the phenyl-addition – dehydrocyclization mechanism, which implies the addition of phenyl radical and formation of triphenylene (Zhao et al. 2019). There are also other ways to form naphthalene from benzene (Mebel et al. 2017); for example the hydrogen abstraction-vinylacetylene addition mechanism, in which the growth proceeds through addition of vinylacetylene instead of acetylene (Parker et al. 2012; Kaiser et al. 2015), but this mechanism is probably operative more efficiently at low temperatures, so we keep it outside the scope of the work.

One more benzene derivative – toluene, forming in the reaction a-C₆H₆ + CH₃ – participates in formation of anthracene and phenanthrene according to Sinha & Raj (2016); Kaiser et al. (2022). Two toluene radicals then recombine and form H-rich complexes with either phenanthrene or anthracene geometry, which subsequently lose their hydrogen atoms and transform into phenanthrene or anthracene. Finally, pyrene is formed from phenanthrene via the HACA mechanism, which has been studied in Frenklach et al. (2019).

Obviously, this reaction network is far from being complete up to pyrene and molecules of comparable size. In addition to the above-mentioned omitted pathways, we have not considered the growth reactions of PAHs with five-membered rings, such as indene, considered in Mebel et al. (2017). Besides the incompleteness, there are several inaccuracies that can be improved in the future. One of them is using a kind of a lumping method widely used in combustion theories (e.g., Frenklach & Wang 1991). We applied it when two isomers (e.g., C₄H₃) potentially had different rate coefficients for the same reactions or even participated in different reactions, but the kinetic parameters for either isomer were not available. In these cases, we used the same rate coefficients for both isomers if they were reactants and a half of the rate coefficient value if they were products.

A general uncertainty for most of the reactions is related to the difference of the atmosphere composition, at which DFT-calculations were performed – it is N₂, characteristic of the Earth, versus H₂ for stellar envelopes. This issue is discussed separately in Appendix A.

Another source of uncertainty stems from the discrepancy between the low pressures typical of stellar envelopes (below ∼10⁻³ atm) and the pressures at which the kinetic rates are typically reported. Although we selected rate coefficients corresponding to the lowest available pressure values, these rarely extend below 0.1 atm. To improve the accuracy of the reaction network, a comprehensive recalculation of rate coefficients at lower pressures would be necessary – a task that remains computationally demanding.

An additional uncertainty arises from the treatment of reverse reaction rates. When these rates are not explicitly provided, we estimated them using the corresponding forward reaction rates and equilibrium constants derived from thermochemical data given in Wang & Frenklach (1997) or GRI-Mech⁵. While this approach prevents the artificial overaccumulation of certain reaction products, it introduces greater uncertainty in our results since our physical conditions are far from thermodynamic equilibrium.

Additionally, we have excluded photoreactions from our model, as we assume that UV radiation plays a negligible role in hydrocarbon chemistry within the studied envelope regions. This assumption is justified by two factors: (1) external UV radiation is absorbed efficiently in the outer regions to be significant, and (2) stellar emission of IRC+10216 in this spectral range is insufficiently intense. While comprehensive photochemical networks exist (notably the UMIST database for stellar envelopes; Millar et al. 2000; McElroy et al. 2013; Millar 2016), these are primarily applicable to outer envelope layers where interstellar UV photons dominate the chemistry – a regime outside our current scope of investigation.

2.5 Treatment of dust evolution

The presence and properties of dust grains are fundamental to this study. As was previously noted, our model does not include homogeneous SiC nucleation. We assume that SiC seed particles first form in the stellar vicinity and present in the envelope at radial distances R ≥ 2 R_⋆, where conditions permit both their survival and subsequent growth. We adopted a size distribution for SiC seed particles following the law n^d(a) ∝ a^−3.5 (the MRN distribution, Mathis et al. 1977). While the MRN distribution was originally proposed for the ISM dust, we employed it here as it represents a well-established standard that has been applied to stellar envelope modeling in previous works (see, e.g., Speck et al. 2009). Test calculations performed under the assumption of an equal grain distribution showed that while the adopted dust size distribution influences hydrocarbon abundances, the general trends and conclusions remain unchanged. These results are provided in Appendix B.

We divided grains into N_a = 15 size bins spanning the range from 0.005 to 0.25 µm. Initially, only bins corresponding to radii a ≤ 0.1 µm were populated, with bins of larger sizes being filled gradually over time due to adsorption and coagulation. Our choice of 0.1 µm as the upper limit for initial seed particles is motivated by several factors. Firstly, theoretical models of dust growth suggest limited efficiency for particles in excess of this size (Cadwell et al. 1994; Krueger & Sedlmayr 1997). Secondly, this value represents a standard choice for grains in monodisperse models or for an upper size limit in models with multiple grain populations used for SED calculations of AGB stars (Lorenz-Martins & Lefevre 1994; Suh 2000). Finally, observations of IRC+10216 indicate that dust sizes are mostly below 0.1 µm (Martin & Rogers 1987). Unlike estimates of the circumstellar dust size, meteoritic studies reveal larger SiC grains, which are believed to have circumstellar origin. Typically, they are 0.3–0.7 µm, but some of them reach several microns (Daulton et al. 2003). However, these measurements likely reflect survival during traveling across the ISM and evolution in the Solar System rather than formation in stellar envelopes.

The dust-to-gas mass ratio in stellar envelopes is estimated to have a range from ∼10⁻³ to 5 · 10⁻³ (Groenewegen et al. 1997; Mattsson et al. 2010; Cherchneff ∼ 2012; Höfner & Olofsson 2018). Within this dust amount, SiC grains typically constitute 1–10% by mass (Speck et al. 2009; Gómez-Llanos et al. 2018), yielding an m_SiC/m_gas ratio between 10⁻⁵ and 10⁻⁴. In our model, we adopt m_SiC/m_gas = 10⁻⁴ as a representative value. Assuming a SiC bulk density of ρ_SiC = 2.5 g cm⁻³ and the MRN size distribution, the ratio m_SiC/m_gas corresponds to an initial seed number density of 5 · 10⁻¹² relative to n_gas. This value agrees well with grain number densities predicted by theoretical models of dust growth (Cadwell et al. 1994; Krueger & Sedlmayr 1997; Gail & Sedlmayr 2013). Note that we assume that homogeneous SiC grain growth terminates once carbonaceous species begin to adsorb onto their surfaces.

The dust temperature was calculated from the equation of thermal balance, assuming that dust is heated by the stellar radiation and by thermal collisions with gas particles and is cooled through infrared emission. Details on the calculations are given in Appendix C. In order to find the dust temperature, absorption and scattering efficiencies (Q_abs and Q_sca) for grains materials are required. Those were calculated from known refractive indices with the SIGMA program⁶ (Lefèvre et al. 2019; Min et al. 2005, 2016). Refractive indices are taken from Laor & Draine (1993) for SiC particles and from Zubko et al. (1996) for amorphous carbon. We distinguish purely SiC grains and grains with a SiC core and carbonaceous mantle. A proportion in which each material contributes to the cross section depends on the mass of adsorbed molecules relative to the core mass and changes at each calculation step. The high temperature of dust grains and correspondingly high desorption rate in the beginning of the pulsation phase make abundances of surface species very small, but they increase again by the end of the phase.

Thorough modeling of dust evolution requires a more sophisticated model. In particular, our model includes dust coagulation but we do not consider the reverse process, i.e., fragmentation due to uncertainties in parameters of this process. To test the possible significance of fragmentation, we performed additional calculations assuming that dust size distribution reverses back to the initial one at the beginning of each shock wave. It was found that our conclusions do not differ substantially in this case, which we show in Appendix B.

In addition to fragmentation, we neglected structural transformations of the material that may occur at elevated temperatures during shock wave passage. While amorphous carbonaceous particles undergo dehydrogenation and graphitization at dust temperatures of T^d ≳ 2000 K (Jones 2022), followed by sublimation of carbon atoms or molecules from the graphitized surface (a process distinct from our adopted desorption treatment), we excluded these effects from our model. This simplification is justified because: 1) the duration of T^d > 2000 K conditions is extremely brief (∼1 s or less); 2) such high temperatures occur only in the immediate stellar vicinity (R ≈ 2.2 R_⋆); 3) graphitization significantly affects only the smallest dust particles and only partially.

Fig. 1 Left: gas-phase abundances of aromatic molecules versus distance from the star calculated with CTM (solid lines) and without it (dotted lines). Right: gas-phase abundances of aromatic molecules at R = 2.7 R⋆. The red color indicates the model with CTM. The blue color shows the model without it.

3 Results

3.1 Abundances of aromatic molecules

We performed calculations of chemical abundances in the envelope with conditions corresponding to IRC+10216 within the region between 1.2 and 5 R_⋆. Models with surface reactions (i.e., with CTM on the dust surface) and without them are considered.

In Fig. 1 (left), we demonstrate gas-phase abundances of aromatic molecules versus distance from the star. The abundances correspond to the end of the adiabatic phase at each position. Both models with CTM (solid lines) and without it (dotted lines) are presented. It is seen that in both cases abundances behave in a similar way. Their maxima fall around 2–3 R_⋆ which is consistent with the conclusion of Cherchneff (2012) that benzene is most efficiently formed at this distance. Abundances mostly are higher in the model with CTM, although some molecules (such as fulvene, naphthalene, and even benzene⁷ at R = 3 R_⋆) can be more abundant in the purely gas-phase model. At R ≳ 3 R_⋆ abundances gradually decrease for two reasons: (1) the rate of consumption of molecules in reactions is higher than their formation rate; (2) the molecule adsorb onto dust surfaces, and their desorption is ineffective at decreasing temperature. Abundances of benzene, fulvene, naphthalene, and biphenyl fall steeper than abundances of other molecules, because the difference between the rates of their consumption and formation is bigger.

At R = 3 R_⋆, the abundance of benzene predicted by the model without CTM is higher than the value predicted by the model with CTM. According to Cherchneff (2012), it is at this point that the abundance of benzene formed in gas-phase reactions reaches its maximum. As our network of reactions up to benzene is very close to Cherchneff’s network, this statement is valid for our work in the model without CTM. However, if CTM is included, then part of acetylene had been already consumed in it or locked in sC₄H₄. Thus, the amount of available acetylene turns out to be lower. At larger distances, when efficiency of gas-phase reactions decreases, the abundance of benzene predicted by the model with CTM exceeds the value predicted by the model without CTM.

To demonstrate difference between models more clearly, we present the abundances of aromatic molecules at R = 2.7 R_⋆ as a histogram in Fig. 1 (right). The most noticeable difference (around an order of magnitude) is seen for benzene, biphenyl, phenanthrene, and pyrene. Thus, the addition of surface CTM accelerates the formation of these molecules. In its turn, naphthalene and its precursors (a-C₆H₅-C₂H₂, a-C₆H₅-C₂H₃) are slightly more abundant in the model without CTM because in this case acetylene is not consumed in CTM and consequently is available for other reactions.

In Fig. 2, we demonstrate the abundance of aromatic molecules on surfaces of dust grains. We remind the reader that it is assumed that there is no dust at R < 2 R_⋆; therefore, surface abundances of molecules are zero. From the left panel of Fig. 2, it is seen that the abundances of some molecules (toluene (a-C₆H₅-CH₃), phenanthrene, pyrene) achieve their maxima at R ≈ 3 R_⋆ and remain constant up to R = 5 R_⋆. Abundances of benzene, fulvene, naphthalene, and biphenyl reach the maxima at R = 2.7 R_⋆, which approximately coincides with the position of their maxima in the gas phase. Further from the star their abundances on the dust surface decrease, which is consistent with that their gas-phase abundances fall steeper than for other molecules at these distances. The histogram on the right panel of Fig. 2 is similar to that in Fig. 1, but abundances of molecules on the surface are lower than in the gas phase at R = 2.7 R_⋆. Note that a detailed analysis of the dust size distribution evolution is beyond the scope of this study; however, a related figure is provided in Appendix D.

Fig. 2 Same as in Fig. 1 but abundances of aromatic molecules on surfaces of dust grains are presented.

3.2 Total abundance of aromatics and its sensitivity to the desorption energy

In order to analyze the results of our modeling and to estimate the contribution of the surface CTM, we calculated the total abundance of aromatics (n_PAH) in gas and on dust at each position. The abundances on dust are presented to show the total abundance of formed aromatics independently of their current phase as the state can change when some influencing factors (such as shock waves and UV radiation) appear. We illustrate the variation of the total abundance with the distance at different E_d in Fig. 3 (on the left). Additionally, we give the values in Table E.1. As in previous figures, the abundances correspond to the end of the adiabatic phase at each position.

Cherchneff (2012) performed simple estimations of the total PAH and carbonaceous dust mass assuming that all benzene convert to coronene molecules and that all these coronene molecules are involved in dust grain growth. She did not perform detailed calculations of all thermal pulses that parcels may undergo. Instead she multiplied the maximum values of PAH and dust masses by an approximate number of thermal pulses, implying that each one leads to formation of the same amount of benzene and, consequently, coronene and more complex carbonaceous dust. The obtained dust-to-gas ratio was found to be in the range 1.2 · 10⁻³–5.8 · 10⁻³, which is consistent with observational values and theoretical predictions of ≈10⁻³–5 10⁻³ (Groenewegen et al. (1997), etc.). We did not repeat the procedure of Cherchneff (2012) because we do not aim to reproduce the observed dust mass. Instead, we focus mainly on the PAH formation process and its variation with distance. Therefore, we discuss the calculated abundances without extrapolation of these values to other pulsations. Thus, it should be kept in mind that our total abundances can be underestimated due to the lower number of pulses included in the model (it is 7) in comparison with the real number of pulses (from 9 to 43 in Cherchneff 2012).

In Fig. 3, n_PAH at E_d = 3 eV reaches its maximum at R ∼ 2 - 3 R_⋆. This region is predicted to be the most favorable for dust formation (Krueger & Sedlmayr 1997; Gail & Sedlmayr 2013). At R = 2.2 R_⋆, aromatic molecules are mainly in gas, but they adsorb intensively on the dust surface when moving away from the star. At R = 5 R_⋆, the gas phase n_PAH is only 2.80 ・ 10⁻¹² relative to n_gas. Certainly, when the parcel undergoes another shock wave, part of aromatics is ejected into gas again. Nevertheless, under the conditions corresponding to R = 5 R_⋆ and at E_d = 3 eV gas components tend to be solidified. Overall, n_PAH reaches 9.86 · 10⁻⁶ at 2.7 R_⋆, which exceeds the maximum abundance of benzene predicted by Cherchneff (2012). Cherchneff (2012)’s value of $n_{{\text{C}}_6{\text{H}}_6}$ is 1.31 10⁻⁶ relative to $n_{{\text{H}}_2}$ ($n_{{\text{H}}_2}$ is generally very close to n_gas), which is closer to the maximum n_PAH predicted by the model without CTM.

It is not possible to compare the predicted PAH abundance with observations as PAHs in stellar envelopes are not excited and consequently are not observed. However, we may compare the predictions for $\hspace{0.17em}{n}_{{\text{C}}_2{\text{H}}_2}$ with observations, although such a comparison is also not straightforward. The observations provide the summed abundances over the phase when shock waves pass through the parcel and most of the surface components evaporate into gas and the stage when gas components adsorb back on the surface. As $\hspace{0.17em}{n}_{{\text{C}}_2{\text{H}}_2}$ at the beginning and at the end of penetration of a shock wave into a parcel are different, we suppose that it is the mean abundance predicted by the model that should be compared with observations. We provide both abundances of C₂H₂ in Table E.1.

Fonfría et al. (2008) reported acetylene abundances of 8 · 10⁻⁶ (relative to H₂) in the inner envelope and 8 10⁻⁵ in outer regions. In the region of efficient PAH and dust formation (R ∼ 2–3 R_⋆), the model with E_d = 3 eV and CTM predicts $\hspace{0.17em}{n}_{{\text{C}}_2{\text{H}}_2}$ in the range 2.92 · 10⁻⁷–10⁻⁴ at the beginning of a shock wave and 2.93 · 10⁻⁹–1.24 · 10⁻⁵ in its end. The range of values is broad; nevertheless, one may conclude that they satisfy the observations. However, at larger distances $\hspace{0.17em}{n}_{{\text{C}}_2{\text{H}}_2}$ becomes less than 10⁻¹⁰, which is quite low. If CTM is not included, predicted $\hspace{0.17em}{n}_{{\text{C}}_2{\text{H}}_2}$ is higher, but at large distances $\hspace{0.17em}{n}_{{\text{C}}_2{\text{H}}_2}$ is also low (up to ≈10⁻⁹). We suppose that such a high acetylene depletion would be corrected if an activation barrier for chemisorption on surfaces were to be added but we leave investigation of this issue for future work.

Let us consider the model with E_d = 2 eV. Such a low E_d extends the region of efficient PAH formation and dust growth from 2–3 R_⋆ to 5 R_⋆. Thus, the clearly defined PAH and dust formation region is absent. Most theoretical models predict that dust formation region centers at 2–3 R_⋆ and ends to 5 R_⋆ (Krueger & Sedlmayr 1997; Gail & Sedlmayr 2013). Our model with E_d = 2 eV may be in slight contradiction with this statement. Nevertheless, these results represent a model with the lower-bound estimate for E_d that remains physically plausible. In this case, maximum n_PAH equals to 7.24 · 10⁻⁶, which is somewhat lower than the maximum value predicted by the model with E_d = 3 eV. However, at 5 R_⋆ n_PAH is 4.19 10⁻⁶, which is ∼5 times higher than n_PAH in the case of E_d = 3·eV. Also, unlike the model with E_d = 3 eV, most molecules are in the gas phase at 5 R_⋆.

The abundance of acetylene predicted by the model with E_d = 2 eV is quite high – ∼7 · 10⁻⁵–10⁻⁴, i.e., it exceeds the observational estimates given by Fonfría et al. (2008). Such a high abundance is consistent with the values of Cherchneff (2012), who did not consider participation acetylene in dust growth. But dust growth is included in our model, and acetylene is believed to be one of the main building block in this process. Thus, the model with E_d = 2 eV does not satisfy the necessary rates of dust growth although it provides a reasonable gas-phase amount of PAHs.

The model with E_d = 4 eV represents an upper-bound estimate for the desorption energy. The results obtained with this model are outstanding in comparison with the models considered above. As soon as dust grains appear in the calculations, i.e., at 2.2 R_⋆, aromatics adsorb on their surface. Consequently, the abundance of molecules in gas (aromatic and non-aromatic hydrocarbons) drops, and all aromatics further present only on the surface. The lack of acetylene and other hydrocarbons in the gas phase results in ceasing of PAH formation. CTM also does not work for the same reason. Therefore, there is no PAH formation and dust growth as most hydrocarbons formed in the locations closest to the star are depleted on SiC grains.

Thus, the high E_d = 4 eV case appears to be implausible. Moreover, the results obtained with this model are inconsistent with observed acetylene abundance: $\hspace{0.17em}{n}_{{\text{C}}_2{\text{H}}_2}$ is ∼10⁻¹⁰ and lower at R ≳ 3 R_⋆.

Summarizing the results presented in Fig. 3, we can argue that the total abundance predicted by different models varies significantly with the adopted value of E_d, establishing it as a critical parameter governing the CTM efficiency. While the E_d = 4 eV case provides an upper bound with extreme results, the models with E_d = 2 and 3 eV model better describe the actual hydrocarbon chemistry evolution. However, we suppose that a realistic compromise falls between these two models.

The uncertainty in observational estimates of the amount of PAHs formed during the AGB phase persists during the post-AGB phase, where the increasing UV radiation from the emerging white dwarf excites PAHs and other small carbonaceous species. This excitation produces diverse mid-infrared spectra, which can be broadly classified into several types (Sloan et al. 2014), though many objects defy these categories. Crucially, the “classical” PAH emission spectra ubiquitous in the ISM (Peeters et al. 2002) are almost not observed in these objects. A possible explanation for this absence is that the newly formed PAHs are not isolated in the gas phase. Instead, they are likely bound to the surfaces of dust grains or agglomerated with each other. They also can be substitutes of grains with mixed and irregular structure such as “aromers” (Homann 1998) or aromatic-aliphatic linked hydrocarbons (see Bockhorn et al. 2009). This state would prevent PAHs from producing their characteristic emission bands. Therefore, accurately modeling PAH evolution requires one to explicitly account for their interaction with solid surfaces.

Based on our results, we conclude that it is crucial to model the chemistry of hydrocarbons and the dust growth consistently because they are tightly related to each other. Consumption of gaseous hydrocarbons in the dust growth process reduces their availability for PAH formation and growth. While in this work all hydrocarbons are consumed equally, the real picture can be more complicated, i.e., some hydrocarbons are consumed more than others. Therefore, it can affect the chemistry substantially. In turn, the results of chemical evolution of hydrocarbons determine the composition and structure of forming dust grains. However, in order to model chemistry of hydrocarbons more accurately, it is necessary to know the desorption energy and also the activation barrier of chemisorption for all considered hydrocarbons and different surface structures, which are not available so far.

Fig. 3 Left: total abundance of all aromatic molecules obtained in the model with CTM. Gas-phase abundances are shown by dashed lines, solid-phase abundances are shown by dotted lines, and the sum abundance is shown by solid lines. Red, black, and blue colors illustrate the models with Ed = 2, 3 and 4 eV, respectively. Right: ratio between the total abundance of all aromatic molecules obtained in the models with the surface acetylene CTM and without it. The designation “GR+SR” means the model with CTM, while the designation “GR” indicates the model without it. The ratio of gas-phase abundances is illustrated by solid lines, while the ratio of solid-phase abundances is shown by dotted lines. The colors correspond to the same models as in the left panel.

3.3 Contribution of cyclotrimerization to the total abundance of aromatics

To quantitatively assess the contribution of the surface acetylene CTM, we compared the total abundance of aromatic molecules for models with and without CTM by computing their ratio at each radial position (n_PAH(GR + SR)/n_PAH(GR)). The variation in this ratio with distance is presented in Fig. 3 (right), where gas and dust abundances are displayed separately. For the model with E_d = 3 eV, the ratio peaks at ≈10 near R = 2.2 R_⋆ before decreasing to ≈0.5 both for gas and dust at 3 R_⋆. At the next computed position, R = 3.5 R_⋆, the ratio again rises up to ≈ 4 and remains constant farther.

The steep fall of the ratio at R = 3 R_⋆ is related to the peak efficiency of benzene production in gas-phase reactions that was mentioned above and seen in Fig. 1 (left). Such a behavior means that benzene is formed in the CTM reaction most efficiently relative to gas-phase reactions at the enhanced temperature corresponding to R = 2.2 R_⋆ and farther, at R > 3 R_⋆. Specifically, at R = 3 R_⋆ the CTM reaction still plays a crucial role in the benzene formation but gas-phase reactions become the most efficient and provide the major contribution to the abundance of aromatics.

As was discussed earlier, the parameter E_d determines the results of our study. At E_d = 2 eV, the CTM contribution becomes negligible at R ≲ 3 R_⋆, primarily due to inefficient sticking of molecules to grains: acetylene molecules do not remain surface-bound long enough to meet another two molecules. At larger distances (R > 3 R_⋆), at which gas and dust temperatures decrease, surface reactions gain importance. The reduced desorption efficiency enhances the probability of successive acetylene collisions on grain surfaces. Concurrently, gas-phase benzene formation becomes less efficient, further amplifying the relative significance of surface CTM. At R = 5 R_⋆, the ratio n_PAH(GR + SR)/n_PAH(GR) reaches a value of ≈9.

In the case in which E_d = 4 eV, CTM plays a negligible role. At R = 2.2 R_⋆, the total aromatic abundance is identical in both models. The abundance of acetylene in gas is reduced by its rapid adsorption, which dominates over desorption. It leads to complete surface coverage by various hydrocarbons along with acetylene. Consequently, when an acetylene molecule adsorbs, the probability of finding sites occupied by sC₂H₂ or sC₄H₄ for CTM is reduced. The sites with sC₂H₂ or sC₄H₄ become rapidly covered by other molecules, terminating the CTM sequence. At R = 2.7 R_⋆, the abundance of aromatics predicted with the model with CTM exceeds the abundance predicted by the model without CTM by a factor of 2. However, this enhancement is negligible indeed and is around computational precision as gas-phase abundances of aromatics are insignificant.

4 Discussion

Let us discuss key physical quantities related to PAHs predicted by our model and compare them with observational estimates from the ISM. Specifically, we shall discuss the PAH number density, the fraction of carbon locked in PAHs, and the fraction of large PAHs.

Observational studies indicate the average PAH abundance in the ISM of approximately n_PAH/n_H ~ 3 · 10⁻⁷ (Tielens 2008; Draine & Li 2007). The total number densities of aromatic species produced in our model are generally consistent with this value. For instance, models with desorption energies of E_d = 2 and 3 eV predict n_PAH/n_gas values on the order of ∼10⁻⁶–10⁻⁵ within the dust formation region. The inclusion of the CTM reaction further enhances n_PAH by a factor of a few, depending on the specific model parameters.

However, a direct comparison of number densities requires caution due to a significant discrepancy in PAH sizes. The PAHs in the ISM are substantially larger than those in our current model. Studies on PAH photodissociation under UV radiation suggest that only species with more than 40–50 carbon atoms (N_C > 40–50) can survive in the harsh ISM environment (Allain et al. 1996; Murga et al. 2019). The limit size is even more drastic if we also take into account interaction of PAHs with shock waves (Micelotta et al. 2010). Our chemical network does not include pathways to PAHs of such a size, nor does it contain all possible reactions leading to even our largest considered species (pyrene, N_C = 16). Nevertheless, preliminary conclusions can be drawn. The fraction of PAHs with N_C ≥ 10 (i.e., up to pyrene) among all formed aromatics is too small, rarely exceeding 10⁻⁴. Note that the fraction increases by several times if CTM is included in case of E_d = 3 eV.

This indicates that the gas-phase growth mechanisms considered here are inefficient. We admit that we underestimate the growth rate. It could be higher by an order of magnitude or so if parcels undergo more thermal pulses than are simulated and more reactions are included up to PAHs of a size comparable to pyrene. But even in this case, the final fraction of large PAHs would be insufficient. Therefore, it is necessary to incorporate alternative growth pathways, such as chemistry on dust grains (Merino et al. 2014; Zhao et al. 2016), PAH dimerization (Webster et al. 2022), polymerization of polyynes (Krestinin 2000), or the partial aromatization of amorphous carbon materials followed by detachment during restructuring events (Marks et al. 2006; Jones et al. 2013). The conclusion about the necessity to include alternative mechanisms will be tested in the future by expansions of our chemical network.

Given the size disparity, a more relevant comparison may be for the fraction of elemental carbon locked in PAHs $\left(f_{\text{C}}^{\text{PAH}}\right)$. Tielens (2008) estimated this value to be ∼3.5 · 10⁻² for the ISM. In our work, $f_{\text{C}}^{\text{PAH}}$ – considering PAHs in both gas and dust phases – varies from 9.20 · 10⁻⁴ to 8.76 10⁻² across different models. The values obtained in the dust·formation region for models with E_d = 2 and 3 eV are comparable to the interstellar value, while the model with E_d = 4 eV yields a value approximately two orders of magnitude lower. Thus, assuming the carbonaceous material formed can be efficiently incorporated into the ISM and that the proportional allocation of carbon is preserved (though not necessarily the specific molecular structures), our model can account for the observed fraction of carbon locked in PAHs. The small aromatics formed in our simulations would be unlikely to survive travel to the ISM in their original form but could serve as building blocks for larger PAHs or other carbonaceous grains.

5 Conclusions

We present a model for the chemical evolution of hydrocarbons in the inner circumstellar envelope of the C-rich AGB star IRC+10216. The primary goal of this work is to incorporate the formation of the first aromatic ring on the surface of dust grains and to evaluate its influence on the resulting abundances of aromatic species. We assume the initial grain material is SiC, which can catalyze the CTM of acetylene to form benzene. Furthermore, we propose that the resulting carbonaceous mantle and growing PAHs continue to catalyze CTM on their carbene-like edges.

This model, named BRAHMA, is based on the framework established by Cherchneff (2012). We adopted their chemical reaction network and the physical conditions (gas density and temperature) derived from a hydrodynamical model that includes periodic shock waves from stellar pulsations. This network was expanded to include (or update, when possible) reactions with hydrocarbons up to pyrene (C₁₆H₁₀). The system of kinetic equations for gas-phase species was extended to include equations for surface species and dust grain populations.

We computed the abundances of aromatic molecules in the gas and dust phases at seven radial positions from 1.2 to 5 R_⋆. Several scenarios were investigated: models with and without the surface CTM reaction, each with varying desorption energies of 2, 3, and 4 eV. Our modeling results lead to the following conclusions:

1. A coherent model that couples gas-phase hydrocarbon kinetics with surface chemistry and carbonaceous dust growth significantly alters the predicted abundances of aromatic hydrocarbons. Consequently, this self-consistent approach is essential for obtaining accurate results;

2. Including surface CTM of acetylene increases the total abundance of aromatics by up to an order of magnitude. This process specifically accelerates the formation of biphenyl, phenanthrene, and pyrene. In contrast, the abundances of naphthalene and acenaphthylene decrease, as their formation pathways require acetylene, which is consumed by the efficient surface reactions;

3. The desorption energy of hydrocarbons from the grain surface is a critical parameter. It directly controls the efficiency of CTM by determining the surface residence time of reactants, which in turn regulates the availability of hydrocarbons for subsequent gas-phase reactions. The activation barrier for chemisorption of hydrocarbons onto the surface is equally crucial, as it controls the rate of dust growth and mantle formation.

In summary, our results indicate that the catalytic CTM of acetylene on both SiC and carbonaceous surfaces is a viable pathway in the circumstellar envelopes of C-rich AGB stars, in particular, IRC+10216. However, its precise contribution depends on several poorly constrained parameters, most notably desorption and chemisorption energies, which we aim to refine in the future.

Data availability

The developed network is publicly available at https://zenodo.org/records/17347157. In total, it includes 194 reagents and 636 reactions.

Acknowledgements

This work was supported by the grant of Russian Science Foundation 24-22-20104, https://rscf.ru/project/24-22-20104.

Reference

Appendix A Difference between rate coefficients at N₂ and H₂ atmospheres

The observed trends in pressure-dependent rate coefficients can be understood through the lens of third-body energy transfer theory, particularly through the work of Jasper (2020). Jasper systematically studied third-body collision parameters across over 300 organic molecules including alkanes, alcohols, and hydroperoxides with several bath gases (H₂, N₂, He, Ar) over temperatures ranging from 300 to 2000 K. His trajectory-based approach yielded generalized expressions for collision cross sections (σ), Lennard-Jones well depths (ϵ), and energy transfer efficiency, often quantified via an effective rotor formalism (N_eff). Notably, Jasper found that for hydrocarbons, energy transfer is consistently 25–40% higher for N₂ than for H₂, depending on temperature and molecular structure. This difference arises because nitrogen, being heavier and more polarizable than hydrogen, has a larger collision cross section and deeper interaction potentials, enabling more effective vibrational energy dissipation during collisions. In contrast, hydrogen’s low mass and weak polarizability result in poorer energy removal efficiency. These findings provide a theoretical benchmark for understanding the role of bath gases in radical recombination reactions.

Figure A.1 illustrates the reaction between benzyl (C₇H₇) and propargyl (C₃H₃) radicals, which serves as a critical pathway for naphthalene (C₁₀H₈) formation, a key step in PAH growth. The initial collision of the radicals may also lead to the formation of reactive intermediates such as methyleneindanyl and methyleneindene species. At higher pressures, frequent collisions with the bath gas (H₂ or N₂) stabilize these intermediates, slowing their progression toward naphthalene. In contrast, at very low pressures, where collisions are rare, the intermediates are less likely to be stabilized, favoring the unimolecular pathway that yields naphthalene and H₂ as the final products. The absolute rate coefficients for intermediate formation are up to two orders of magnitude higher than those for naphthalene across the pressure range, but their pressure-dependent trends mirror those of the naphthalene pathway, making the latter a representative case for analysis.

The C₁₀H₁₀ potential energy surface (PES) accessed by the benzyl–propargyl radical reaction has been described in detail in previous publications (Matsugi & Miyoshi 2012; Krasnoukhov et al. 2022; He et al. 2023). The kinetics and mechanisms of the recombination reaction between benzyl (C₇H₇) and propargyl (C₃H₃) radicals have been theoretically investigated using B3LYP, CBS-QB3, and CASPT2 quantum chemical methods. Additionally, a steady-state unimolecular master equation analysis was conducted based on Rice-Ramsperger-Kassel-Marcus (RRKM) theory, with N₂ as the bath gas. In this study, we report rate coefficients results using H₂ as the bath gas in comparison to N₂, employing the same kinetic methodologies used previously. Figure A.2 presents the pressure-dependent rate coefficients for naphthalene formation at temperatures ranging from 800 to 1250 K. At the lowest pressures (p < 10⁻³ Torr), the rate coefficients for both H₂ and N₂ converge, as the infrequency of collisions becomes the limiting factor, rendering the identity of the bath gas irrelevant. This behavior is typical in astrophysical environments like the ISM, where particle densities are extremely low. As pressure increases into the intermediate range (10⁻² to 100 Torr), the influence of the third-body properties becomes apparent: N₂ consistently yields higher rate coefficients than H₂ due to its superior energy transfer efficiency. At even higher pressures, the rate of naphthalene formation decreases further, reflecting the reduced efficiency of the ring-closure pathway under conditions favoring collisional stabilization of intermediates.

The choice of bath gas is context-dependent. In combustion and atmospheric modeling, nitrogen is the standard third body due to its abundance and moderate efficiency. In astrochemical environments, however, H₂ dominates as both a collider and a potential reactant, particularly at low pressures where its chemical reactivity can open secondary pathways. This underscores the need for kinetic models that account for both collisional stabilization and chemical reactivity, especially in extreme low-pressure regimes where classical third-body formalisms may fall short.

The pressure-dependent kinetics of the C₇H₇ + C₃H₃ reaction not only align with Jasper’s framework for third-body efficiencies but also extend its applicability to a new mechanistic class involving reactive radical intermediates and aromatic ring closure. By corroborating Jasper’s findings, specifically, the superior energy transfer efficiency of N₂ over H₂, this study bridges his generalized collision parameters to complex PAH-forming reactions. The results highlight how bath gas properties influence not just simple stabilization but also the branching between collisional quenching and chemically driven pathways in radical-aromatic systems, a critical consideration for modeling PAH growth in both combustion and interstellar environments.

Appendix B Variations in a-C₆H₆ abundance with different treatments of dust description

In Fig. B.1 we show the benzene abundance predicted by different versions of the model. Solid lines represent various dust size distributions and demonstrate that the choice of size distribution influences the resulting benzene abundance. Specifically, using a uniform dust size distribution instead of the standard MRN distribution yields a higher benzene abundance. Conversely, if dust grains are assumed to fragment back to the original MRN distribution at the start of each shock wave, the benzene abundance is lower than in models where dust is not fragmented.

Despite these differences, the main conclusion of the paper remains robust: the addition of surface CMT to the reaction network increases the abundance of benzene. This is evidenced by the fact that the solid lines (with CMT) of a certain color are located above the corresponding dotted lines (without CMT), at least within the dust formation region (2–3 R_⋆). Therefore, our findings are not sensitive to the adopted treatment of the dust size distribution.

Appendix C Dust temperature

The dust temperature is calculated from the equation of thermal balance between heating by stellar radiation and due to collisions with gas particles and cooling due to thermal emission (see Sect. 24.1 in Draine 2011 and Sect. 9.2 in Bovino & Grassi 2024): $$\begin{array}{l}2\pi {\displaystyle \underset{0}{\overset{\infty }{\int }}{Q}_{\text{abs}}}\left(v,a\right)J_{*}\left(v,R\right)dv+n_{gas}{v}_{\text{gas}}^{\text{th}}\gamma k_{\text{B}}\left(T_{\text{gas}}-T^{\text{d}}\right)=\hfill \\ 2\pi {\displaystyle \underset{0}{\overset{\infty }{\int }}{Q}_{\text{abs}}}\left(v,a\right)B\left(v,T^{\text{d}}\right)dv\hfill \end{array}$$(C.1)

Fig. A.1 Naphthalene formation scheme via reaction of benzyl and propargyl radicals, highlighting initial adducts and final product.

Fig. A.2 Pressure-dependent naphthalene formation rate coefficients at various temperatures, showing convergence at low pressures and divergence at higher pressures.

Fig. B.1 Gas-phase benzene abundances at various distances from the star. The black solid line is the result of the “standard” model with CMT (Fig. 1). The red solid and dotted lines illustrate the results obtained for the uniform initial size distribution. The blue solid and dashed lines correspond to the cases when dust grains fragment back to the initial size distribution (MRN) at the start of each shock wave.

Fig. D.1 Dust size distribution at different positions. The initial distribution is illustrated by black color, colors from brown to orange illustrate the positions indicated in the legend.

where Q_abs(ν, a) is the absorption efficiency for the grain of radius a at frequency ν, J_⋆(ν, R) is the intensity of the stellar radiation at frequency ν and radial distance R, B(ν, T^d) is the Planck function, γ is the thermal accommodation coefficient (adopted to be 0.5 following Draine 2011). It is assumed that dust grains are mostly heated by collisions with hydrogen molecules, so we do not integrate heating rates over all gas-phase species.

We find J_⋆(ν, R) via $$J_{\star }\left(v,R\right)=\frac{1}{2}B\left(v,T_{\star }\right)\left(1-\sqrt{1-{\left(\frac{R_{\star }}{R}\right)}^2}\right)\text{exp}\left(-\tau \left(v,R\right)\right),$$(C.2) where τ(ν, R) is the optical depth which is calculated as $$\tau \left(v,R\right)={\displaystyle \sum _{i=1}^{N_a}{\displaystyle \underset{R_{*}}{\overset{R}{\int }}\pi a_j^2{n}_i^d\left(R\right)n_{gas}{Q}_{\text{ext}}\left(a_i\right)dr,}}$$(C.3) where Q_ext is the extinction efficiency, that is, the sum of absorption and scattering efficiencies.

Appendix D Dust size distribution

In Fig. D.1 we demonstrate the change of the dust size distribution during evolution of the envelope. Initially there are only grains smaller than 0.1 µm, and they are distributed by the MRN distribution. Further, dust grains grow via adsorption and coagulation. Consequently, the dust size distribution extends to larger grains, herewith the fraction of small grains decreases.

Appendix E Calculation results: Extension

In this section we present Table E.1 with PAH abundances calculated with different versions of our model. The table includes the following columns: 1 – the distance from the star; 2 – the model version; 3-5 – the total PAH abundances in gas, on dust surfaces and summed, respectively; 6-7 – the gas-phase abundance of acetylene at the beginning and in the end of propagation of a shock wave; 8-10 – the fraction of carbon in PAHs in gas, on dust surface and summed, respectively; 11 – the fraction of PAHs with N_C > 10 among all PAHs. All abundances are given relative to n_gas.

All Tables

All Figures

	Fig. 1 Left: gas-phase abundances of aromatic molecules versus distance from the star calculated with CTM (solid lines) and without it (dotted lines). Right: gas-phase abundances of aromatic molecules at R = 2.7 R⋆. The red color indicates the model with CTM. The blue color shows the model without it.
In the text

	Fig. 2 Same as in Fig. 1 but abundances of aromatic molecules on surfaces of dust grains are presented.
In the text

Fig. 3 Left: total abundance of all aromatic molecules obtained in the model with CTM. Gas-phase abundances are shown by dashed lines, solid-phase abundances are shown by dotted lines, and the sum abundance is shown by solid lines. Red, black, and blue colors illustrate the models with Ed = 2, 3 and 4 eV, respectively. Right: ratio between the total abundance of all aromatic molecules obtained in the models with the surface acetylene CTM and without it. The designation “GR+SR” means the model with CTM, while the designation “GR” indicates the model without it. The ratio of gas-phase abundances is illustrated by solid lines, while the ratio of solid-phase abundances is shown by dotted lines. The colors correspond to the same models as in the left panel.

In the text

	Fig. A.1 Naphthalene formation scheme via reaction of benzyl and propargyl radicals, highlighting initial adducts and final product.
In the text

	Fig. A.2 Pressure-dependent naphthalene formation rate coefficients at various temperatures, showing convergence at low pressures and divergence at higher pressures.
In the text

Fig. B.1 Gas-phase benzene abundances at various distances from the star. The black solid line is the result of the “standard” model with CMT (Fig. 1). The red solid and dotted lines illustrate the results obtained for the uniform initial size distribution. The blue solid and dashed lines correspond to the cases when dust grains fragment back to the initial size distribution (MRN) at the start of each shock wave.

In the text

	Fig. D.1 Dust size distribution at different positions. The initial distribution is illustrated by black color, colors from brown to orange illustrate the positions indicated in the legend.
In the text