Master Equation Modelling of Non-equilibrium Chemistry in Stellar Outflows

A current challenge in astrochemistry is to explain the formation of Fe-Mg silicate dust around evolved stars. The dust is observed to form within 2 to 3 stellar radii of oxygen-rich AGB stars, where the typical conditions are kinetic (translational) temperatures between 1200 and 1600 K, and total gas densities below 1011 cm-3. At these high temperatures, molecules with bond energies < 400 kJ mol-1 should be short-lived, and this results in kinetic bottlenecks in postulated mechanisms for converting the observed Fe, Mg, SiO and H2O into silicate. Here we show that, in the very low pressure regime of a stellar outflow, molecules can exhibit significant vibrational disequilibrium because optical transitions - both spontaneous and stimulated by the stellar radiation field - occur on a much faster timescale than collisions. As a result, relatively less stable molecules can form and survive long enough to provide building blocks to silicate formation. Here we use the molecule OSi(OH)2, formed by the recombination of SiO2 and H2O, as an example. When vibrational disequilibrium is accounted for in a master equation treatment which includes optical transitions, the quantity of metal silicates produced in a low mass loss rate evolved star (R Dor) is increased by 6 orders of magnitude.


Introduction
The chemical evolution of galaxies depends on mass loss from stars during the late stage of their evolution: asymptotic giant branch (AGB) stars lose a signicant fraction of their mass in a short period, during which they tend to be surrounded by an optically thick shell of dust. 1 The absorption and emission of radiation by this circumstellar dust creates a radiation pressure which induces a substantial acceleration in the outow velocity of the dust (and gas), allowing escape from the star's gravitational eld.The dust is also a signicant source of interstellar dust particles on which heterogeneous chemical synthesis can occur, and is central to the formation of new planets around young stars.
In oxygen-rich outows the dust is observed (by emission spectroscopy in the IR and far-IR) to be mostly composed of the Mg-Fe-silicates olivine and pyroxene (Fe 2À2x Mg 2x SiO 4 and Fe x Mg 1Àx SiO 3 , respectively, where 0 # x # 1). 2,3However, a continuing challenge in astrochemistry is to explain how these metal silicates form from the componentschemically inert Fe, Mg, SiOthat are observed in the inner wind region close to the star.6][7][8][9] Fig. 1(a) illustrates the typical kinetic (i.e.translational) temperature and total gas density (mostly H 2 ) modelled as a function of distance r (in units of stellar radius, R*) for the semi-regular variable AGB star R Dor. [10][11][12] Fig. 1(b) shows the concentration proles of H, H 2 O, OH and O as a function of r, calculated using standard outow chemistry. 6The inner edge of observable dust clouds is typically where the temperature in the outow is between 1400 and 1000 K, in this case beyond 2R*.There are several reasons why this is a challenging environment in which to form metal silicates.Firstly, the precursor molecules need to be very stable to survive thermal dissociation.Secondly, the major constituents are H 2 and H; these are chemically reducing species, which reduce metal oxides and hydroxides back to metal atoms, and SiO 2 back to SiO. 6 Indeed, we have shown recently using the Atacama Large Millimeter/submillimeter Array (ALMA) that metal oxides and hydroxides are not detectable in high sensitivity observations of R Dor; the upper limit abundances are in accord with our experimental and theoretical work on the reaction kinetics of these species with H and H 2 . 13Thirdly, the total pressure in the stellar outow at 2R* is less than $0.003Pa, which will slow down the rates of pressuredependent recombination reactions.Lastly, the rates of formation of silicates from species such as Mg, Fe and SiO will have at least a second-order dependence on the total density in the outow, and decrease rapidly as the ow expands away from the star. 6t has therefore been concluded that metal silicates are unlikely to form directly in an outow at temperatures above 1000 K, [7][8][9] so that attention has focused on more refractory alternatives in an O-rich stellar outow, such as corundum (Al 2 O 3 ), 11 titania (TiO 2 ), 14 and perovskite (CaTiO 3 ). 6If small "seed" particles of these substances were able to form in the inner wind at relatively high temperatures, then they might act as condensation nuclei for metal silicates in the cooler outer wind.However, this would still require metal silicate molecules to form at some point in the outow.
In this discussion paper we explore the role of vibrational disequilibrium, which we dene here as the internal vibrational temperature of a molecule being signicantly different from its kinetic temperature.Vibrational disequilibrium can arise when the gas density is low enough so that optical transitions occur at a faster rate than collisional energy transfer.Nuth and co-workers pointed out in the 1980s that the assumption of local thermodynamic equilibrium (LTE) in some astrophysical regions may not be valid. 15,16This idea does not appear to have been pursued further, probably because the assumption of LTE greatly simplies astrochemical modelling.However, observations in this century are clearly showing that molecules in the inner winds of AGB stars are far from LTE.For example, a very recent paper from Fonfría et al. 17 shows that at 2R* around the AGB star Y CVn the vibrational temperature of HCN is only 590 K, compared with a kinetic temperature of 1800 K (see Fig. 9 of that paper); similarly low vibrational temperatures have been observed for C 2 H 2 and C 2 H 4 . 18,19ere we consider the effect of vibrational disequilibrium on the formation and dissociation of the dehydrated form of silicilic acid, OSi(OH) 2 .This molecule is postulated to form in an outow from the recombination reaction 6 where the reaction enthalpy change is calculated at the CBS-QB3 level of theory. 20,21Metal atoms, such as Ca, Fe and Mg, can then react with OSi(OH) 2 to form pyroxene-like molecules, e.g.
Two problems have been identied with this (otherwise) straightforward route to forming silicates. 6First, the only apparent route to forming SiO 2 is oxidation of SiO by OH: This reaction has been studied at room temperature and is reasonably fast. 22owever, because the reaction is only slightly exothermic, the very large concentration of atomic H in an outow (e.g.Fig. 1(b)) rapidly destroys SiO 2 through the reverse reaction.The second problem is that reaction (1) is only exothermic by around 280 kJ mol À1 . 6Thus, assuming it is in LTE, OSi(OH) 2 should have a relatively short lifetime against thermal dissociation at the high temperatures in the inner wind: In this study we rst determine the Einstein coefficients for spontaneous emission (A) and stimulated emission/optical absorption (B) for the vibrational infrared active transitions of OSi(OH) 2 .These are then incorporated into a master equation treatment to estimate the non-LTE rates for reactions ( 1) and (À1), as a function of radial distance r from the star.The resulting rates are then used in a simple stellar outow model to show that the non-LTE treatment produces a substantial increase in the formation of metal silicates in the outow.
It is interesting to note that in the Faraday Discussions meeting of 1922 there was a vigorous debate about whether dissociation reactions were caused by collisions, the side championed by F. A. Lindemann, or by infra-red photon absorption (from the sun, or reactor walls)the "radiation theory of chemical action" championed by the French physical chemist J. Perrin. 23The story of the relatively short-lived radiation hypothesis is well described by King and Laidler. 24ur aim is not to resurrect it here as originally proposed, but rather to show that in the very unusual environment of a stellar outow optical transitions can have an important inuence on the kinetics of recombination and dissociation.

Photochemistry of OSi(OH) 2
The harmonic vibrational frequencies and infra-red intensities of OSi(OH) 2 were calculated at the B3LYP/6-311+g(2d,p) level of theory, using the Gaussian suite of programs. 20The absorption cross section spectrum is illustrated in Fig. 2, and the vibrational frequencies are listed in Table 1 (note that the 4 th mode is infra-red inactive).
The integrated cross section s int (in units of cm 2 Hz) for each vibrational mode was used to calculated the Einstein coefficients for spontaneous and stimulated emission (A ij and B ij , respectively) from the relations: 25 where v is the band frequency (Hz), c is the speed of light, and h is Planck's constant.The rates of absorption and stimulated emission were estimated by multiplying the B ij coefficients by a stellar irradiance ux from the MARCS database 26 for an evolved star with T* ¼ 2500 K, which is close to the surface temperature of R Dor. 27 The irradiance range was extended to 80 000 nm by tting a Planck distribution curve at 20 000 nm, the upper wavelength limit of the MARCS model.The values of A ij , B ij (in units of m 2 J À1 s À1 ) and the product of B ij with the irradiance I (in units of J m À2 s À1 Hz À1 h J m À2 ) at 1R* are listed in Table 1 for each vibrational mode.Note that close to the stellar surface the rates of  spontaneous and stimulated emission are fairly similar, apart from the highest frequency modes 11 and 12 where spontaneous emission is much faster.

Optical transitions in the master equation
The Einstein coefficients in Table 1 are for excitation from the ground state to the rst vibrational level of each mode.Here we assume that these coefficients remain constant for transitions between higher vibrational levels in each manifold.We also do not include overtone transitions (which would be formally correct for the harmonic oscillator approximation).The justication for these assumptions is that there is some evidence that the Einstein coefficients actually increase somewhat going up a vibrational manifold (which would enhance the role of optical transitions even more), and overtone transitions are much less probable. 28ormation or dissociation of OSi(OH) 2 is assumed to proceed via an excited adduct, which can either dissociate to SiO 2 + H 2 O or be stabilized by collision with the H 2 third body and spontaneous/stimulated emission, depending on the third body pressure and stellar radiation eld i.e., the reaction occurs in a single well.In practice, the potential energy surface connecting SiO 2 and H 2 O with OSi(OH) 2 is somewhat more complicated: an initial SiO 2 -H 2 O complex rearranges over a barrier to produce an asymmetric form of OSi(OH) 2 , which then undergoes an internal torsional rotation to produce the most stable symmetric form.The complete surface is shown in the ESI, † along with all the relevant molecular parameters for carrying out a full multi-well master equation analysis.In fact, the SiO 2 -H 2 O complex and subsequent barrier are 89 and 59 kJ mol À1 below the entrance channel, respectively; we used the Master Equation Solver for Multi-well Energy Reactions program (MESMER) 29 to show that these features are sufficiently submerged to have a very small effect (<1%) on the overall kinetics under the conditions of interest in this study.Hence, reaction (1) is treated here as a simple recombination involving a single well containing the more stable OSi(OH) 2 isomer.
The internal energy of the OSi(OH) 2 adduct was divided into a contiguous set of grains (width ¼ 40 cm À1 ) containing a bundle of rovibrational states.Each grain was then assigned a set of microcanonical rate coefficients for dissociation, which were determined using inverse Laplace transformation to link them directly to k rec,N , the high-pressure limiting recombination coefficient.This was estimated using long-range transition state theory 30 to be k rec,N ¼ 8.7 Â 10 À10 (T/ 1200 K) 1/6 cm 3 molecule À1 s À1 .The density of states of each adduct was calculated with the vibrational frequencies in Table 1, without making a correction for anharmonicity, and a classical densities of states treatment for the rotational modes (see the ESI † for all rotational constants and vibrational frequencies used, including those of SiO 2 and H 2 O).The probability of collisional transfer between grains was estimated using the exponential down model, where the average energy for downward transitions, hDE down i, was assigned a value of 200 cm À1 for H 2 at 300 K and a small temperature dependence of the form T 0.25 . 31The probabilities for upward transitions were determined by detailed balance. 31,32The collision rate of H 2 with the adduct as a function of temperature was calculated using Lennard-Jones parameters to characterise the intermolecular potential, giving Z(T) ¼ 8.4 Â 10 À10 (T/1200 K) 1/2 cm 3 molecule À1 s À1 .
The master equation, which describes the evolution with time of the adduct grain populations, is given by: Pr ij r j ðtÞ À ur i ðtÞ À k i r i ðtÞ þ X 11 s¼1 ðA i;iþd s r iþd3 s ðtÞ À A iÀd3 s ;i s r i ðtÞ þ I s ðB i;iþd3 s s r iþd3 ðtÞ À ðB where r i (t) is the grain population; u is the frequency of collisions of OSi(OH) 2 with the bath gas (i.e.Z(T) [M], where M ¼ H 2 ); Pr ij is the probability of transfer from grain j to grain i on collision with M; and k i is the microcanonical rate coefficient for dissociation of the adduct.s is the vibrational mode number of OSi(OH) 2 (the sum is over 11 because the 4 th mode is inactive -Table 1).Optical transitions for each mode can only occur between states that are separated by a xed energy, d3 s , which depends on the frequency of mode s.So for the spontaneous transitions, state i can receive density via an optical transition from a grain that is i + d3 s higher in energy (Einstein coefficient A i,i+d3 s s ), and lose density to a grain that is d3 s lower in energy (Einstein coefficient A iÀd3 s ,i s ).The stimulated transitions can go both ways, so grain i can gain population by absorption from grain i À d3 s (Einstein coefficient B i,iÀd3 s s ) and stimulated emission from grain i + d3 s (Einstein coefficient B i,i+d3 s s ), but also lose population via absorption from grain i to i + d3 s or emission from state i to i À d3.I s is the radiance.Lastly, R i is the rate of population of grain i if recombination is occurring.The master equation was then expressed in matrix form and solved via diagonalization to yield the recombination rate coefficient (k 1 ) or dissociation rate coefficient (k À1 ) at a specied pressure and temperature.

Stellar outow model
The stellar outow is modelled here using a b-velocity law, which describes the net acceleration of the wind without explicitly considering pulsations: where v N is the terminal velocity, which is 5.7 km s À1 for R Dor, 27 v 0 ¼ v(R*) ¼ 1.5 km s À1 is the initial velocity at the stellar surface r ¼ R* (¼ 2.5 Â 10 8 km), and the parameter b ¼ 1.The temperature prole is given as a power-law with an exponent a ¼ 0.6 and the density n(r) is expressed using the pressure scale height H ¼ RTR* 2 /(mM*G) (where R is the gas constant, m is the molecular mass of H 2 , M* the mass of the star, and G the gravitational constant), the ideal gas law, and parameter g (¼ 0.89): 5 and The surface pressure for R Dor is taken as n* ¼ 10 14 cm À3 and the effective temperature T* is 2400 K. 11,27 The chemistry network in the model is taken from Plane. 6Measured rate coefficients are adopted from the NIST kinetics database, 33 and those of other reactions (e.g.reactions ( 2), ( 3) and (À3)) were calculated 6 using electronic structure theory combined with MESMER; 29 these are listed in the ESI.† The model predictions of the total gas density, kinetic temperature, and major constituent concentrations are shown in Fig. 1. the kinetic temperature is 1583 K. Also shown is the effective internal temperature of the molecule, determined from the relative population of the grains in the molecular well.The solid lines show the case when optical transitions are turned on.As [M] approaches 10 14 cm À3 , both the dissociation rate and internal temperature increase to their respective thermal-only values (shown by the dashed lines which are the model run with optical transitions turned off).This is because the rate of collisional energy transfer now dominates over optical transitions.In contrast, when [M] is less than 10 9 cm À3 the internal temperature and dissociation rate become independent of [M], because energy transfer between grains is now controlled by optical transitions rather than collisions.Note that the internal temperature is now only 660 K. Fig. 3(b) shows how the internal temperature of OSi(OH) 2 and its rst-order dissociation rate vary as a function of stellar intensity ([M] is xed at 1.2 Â 10 11 cm À3 ).As the star gets brighter, the internal temperature and dissociation rate increase.The dissociation rate exceeds that with optical transitions turned off when the intensity is increased by a factor of 42, and the internal temperature exceeds the translational temperature when the radiance is increased by a factor of 520.In contrast, if the star gets very dim (radiance factor < 10 À3 ), the dissociation rate and temperature become constant because spontaneous emission keeps the molecule cool, but collisional excitation is still occurring.

Results and discussion
Fig. 4 shows the second-order rate coefficient for the recombination of SiO 2 and H 2 O (reaction (1)), as a function of [M].The dashed line is the case for optical transitions turned off, where the reaction exhibits the expected fall-off behaviour as [M] increases.However, when optical transitions are included (solid line), the rate coefficient levels off at low [M] (<10 9 cm À3 ) to a constant value, which is the radiative recombination rate.Once [M] > 10 14 , the effect of optical transitions is minimal and collisional stabilization dominates.At this point the ratio of the recombination rate to the dissociation rate equals the thermal equilibrium constant of reaction (1), consistent with the internal temperature approaching the kinetic temperature (Fig. 3(a)).disequilibrium of the molecule that is comparable to that observed recently for HCN, 17 as discussed in the Introduction.Fig. 5(b) shows that when optical transitions are included, the dissociation rate of OSi(OH) 2 becomes orders of magnitude lower because of the much colder internal temperature of the molecule.It is only close to the star that the H 2 density is high enough for collisional energy transfer to warm the molecule to the point where the dissociation rate approaches that when optical transitions are turned off in the model.The effect on the second-order recombination rate of SiO 2 and H 2 O is illustrated in Fig. 5(c).Moving away from the star, the decrease in [H 2 ] by 4 orders of magnitude between 1.5 and 5R* is only partially offset by the decrease in temperature, so that k rec decreases by $3 orders of magnitude.In contrast, when optical transitions are turned on, there is very little change in the recombination rate with r, because radiative recombination quickly becomes the dominant process, particularly as the stellar radiance is decreasing as r 2 so that spontaneous emission becomes more important than absorption.
The rate coefficients plotted in Fig. 5(b) and (c) were parameterised as a function of r (see Table S2 in the ESI †) and then input into the stellar outow model (Methods).Fig. 6 illustrates the predicted concentrations of the three metal silicates -CaSiO 3 , FeSiO 3 and MgSiO 3as well as OSi(OH) 2 as a function of r.The black and red lines show the model run with and without optical transitions, respectively.When optical transitions are included for reactions (1) and (À1), the OSi(OH) 2 concentration is on average 9.0 Â 10 5 times higher between r ¼ 1.5 and 5R*.This accounts for the metal silicate concentrations being around 10 6 times larger when the OSi(OH) 2 optical transitions are turned on.Note that CaSiO 3 is predicted to be the dominant silicate that forms in the outow: on average it is 3 Â 10 4 higher in concentration than FeSiO 3 , and 3 Â 10 7 higher than MgSiO 3 , in this part of the outow.The reason is that reaction (2) does not have an energy barrier, whereas the analogous reactions of Fe and Mg with OSi(OH) 2 have barriers of 31 and 51 kJ mol À1 , respectively (see the ESI †). 6Fe, Mg, SiO and H 2 O probably then add to CaSiO 3 , producing mixed metal pyroxenes and olivines.This chemistry will be the subject of future work.

Conclusions
This study is a rst step in incorporating optical transitions into a master equation treatment of recombination and dissociation reactions, motivated by the observation of signicant molecular vibrational disequilibrium in stellar outows.For this exploratory work we have made various assumptions and approximations, the most important being that the vibrational modes are treated as harmonic oscillators with Einstein coefficients that remain constant in each of the vibrational manifolds.Nevertheless, this model is able to reproduce the magnitude of vibrational disequilibrium that has been observed in small molecules. 17,18If a molecule is vibrationally cool then its dissociation rate will be reduced compared to the dissociation rate if the molecule were vibrationally equilibrated with the high kinetic temperatures in this environment.In addition, radiative recombination becomes more important than collisional energy transfer at these low pressures, so that the recombination rates are greatly enhanced over what would be expected by extrapolating from the higher pressures where these reactions tend to be studied in the laboratory.An important astronomical result is the substantial increase in the concentration of metal silicates that are predicted to form via this mechanism involving OSi(OH) 2 when optical transitions are included.The modelled silicate mixing ratio with respect to H 2 is 6 Â 10 À11 , which would produce a relatively small amount of dust. 6However, R Dor has a relatively low mass loss rate of only 2 Â 10 À7 M per year (where M ¼ solar mass), whereas mass loss rates range from 2 Â 10 À8 to 4 Â 10 À5 M per year for M-type AGB stars. 34Because the rates of silicate-forming reactions like reaction (2) depend on the square of the concentrations of the relevant species (metals and SiO) in the outow, 6 for a high mass loss rate star of 1 Â 10 À5 M per year the predicted silicate mixing ratio would be $2500 times larger.

Fig. 1
Fig. 1 (a) Predicted total gas density and kinetic temperature as a function of radial position (r, in units of stellar radius R*) for the b-velocity trajectory model of the SRV star R Dor.(b) Predicted concentration profiles of H, H 2 O, OH and O.

Fig. 3 (
Fig. 3(a) shows how the rst-order dissociation rate dissociation rate of OSi(OH) 2 varies as a function of [M], for the conditions of the star R Dor at r ¼ 2R*, where

Fig. 4
Fig. 4 Second-order recombination rate coefficient for reaction (1) as a function of [M].The stellar radiance is for R Dor at 2R*.

Table 1
Vibrational frequencies, Einstein coefficients and optical absorption/stimulated emission rates for OSi(OH) 2Mode numberFrequency/cmÀ1 A ij /s À1 B ij /m 2 J À1 s À1 B ij I b /s À1a Not infra-red active.b Optical absorption/stimulated emission rate at 1R* for the AGB star R Dor.