AGB Stars and Their Circumstellar Envelopes: An Operative Approach to Computing Their Atmospheres

: The study of AGB stars necessarily covers a wide range of topics, from the primary astronomical observations to their interpretation in terms of fundamental physics. All that requires proper ad hoc methodologies, among which numerical modeling of the outer layers of AGB stars plays a paramount role. In this paper, we present an iterative sequential procedure, operative and physically sound, tailored to compute extended stellar atmospheres. It will constitute the backbone of the in ﬁeri TEIDE package to be implemented into our VULCAN code. Such an improvement will allow us to compute more realistic models of the extended atmospheres of AGB stars taking into account important physical aspects that are neglected in the actual version of VULCAN.


Introduction
The modeling of AGB stars requires reliable input parameters (opacities, nuclear inputs, etc.) as well as a detailed knowledge of macroscopic phenomena (equation of state, convection, rotation, stellar pulsations, magnetic fields, etc.). The convective envelopes of AGB stars are characterized by a rich kind of nucleosynthesis: chemical elements are produced in the interiors, which are carried later to the surface by means of third dredge-up (TDU) episodes. As a matter of fact, AGB stars are the major polluters of the inter stellar medium (see, for example [1,2]). Those objects produce half of the elements heavier than iron present in the Universe. Moreover, they are the principal dust polluters of the interstellar medium. As a consequence, a detailed knowledge of their mass-loss mechanism is a necessary condition to properly model the chemical evolution history of galaxies. Depending on the amount of dredged up material, AGB stars may be O-rich (C/O<1) or C-rich (C/O>1). According to the availability of carbon seeds, different molecules form in the most external layers of AGB stars (amorphous carbon and SiC for the C-rich regime; oxides, silicates and spinel for the O-rich regime). Finally, if the temperature drops below 2000 K and the density is still high enough, dust condensation from clusters of molecules begins. Their formation is of paramount importance for AGB evolution, because they determine the ultimate stage of those objects by regulating the rate of their mass loss. Specifically, the gas is dragged along with the shock-accelerated dust grains. The last two decades have been witness to dramatic progresses in the field, as proved by a large quantity of published papers ( [3][4][5][6][7][8][9][10] to quote some of them).
For many years we have carried on the study of the physical and chemical evolution of AGB stars (the relevant models are available online in the FRUITY database web pages 1 [11][12][13]). In FRUITY models the mass-loss rate is estimated by correlating the pulsation period to the variations of surface stellar properties, such as temperature and luminosity [14]. On the other hand, external layers may become unstable against large amplitude pulsations (e.g., Mira variables) that introduce recurrent compression of the gas, which returns to its original radial structure owing to gravity [15]. As already said, these shocks lead to the formation of dust grains that absorb or scatter stellar photons with possible transfer of momentum to the gas (e.g., [6]). The details of the dust nucleation path are still largely unknown to date [16,17].
Several chemical-kinetic models have been presented so far, but they either consider steady state gas conditions [10] or take into account the solution of simplified hydrodynamic equations [17]. On the other hand, highly sophisticated hydrodynamic models including the radiative transfer equations have been published [4][5][6]18,19]. These models, however, adopt a simplified chemical equilibrium and a classical nucleation theory to treat molecules and dust. To our knowledge there are no self-consistent AGB hydrodynamic wind models including a detailed formation path for molecules and dust. Moreover, the aforementioned circumstellar envelope models have never been coupled with the evolution of the inner stellar layers and this requires a completely different model to deal with the evolution of the star. Such a coupling would definitely improve our understanding of the mass-loss mechanism and allow a better determination of AGB chemical yields. As a first step we developed the 1D lagrangian hydrodynamic VULCAN code [20] that we are now improving with the inclusion of the radiative transfer equations and a proper formation path for molecules and dust grains, with the final goal of integrating it in the framework of FRUITY models.
The complex interaction among the physical processes involved constitutes a difficult problem that is strongly non-linear in character. The transport of energy across the outer stellar layers also makes the global problem be non-local. Convective transport takes place under certain physical conditions in the atmospheres of cool stars, but radiative transfer (RT) is always in effect in stellar atmospheres. Efficient algorithms must, therefore, be envisaged for the numerical solution of the RT equation(s). Modeling stellar atmospheres unavoidably implies struggle with grids of huge dimension, both in depth and frequency. The curse of dimensionality 2 is one of the drawbacks to be faced. Another is the above mentioned entanglement of physical processes. The former dictates that in large scale numerical computation the effectiveness of the algorithm must become computationally efficient. On the other hand, brute force approaches to get rid of the latter fails and smart strategies based on a careful analysis of the physics that governs the structure of the system are required. According to our long-lasting experience, the only effective way to get rid of this difficulty is to dissect the structure of a stellar atmosphere into its constituent building blocks. The result of this ideal dissection is a set of blocks such that each of them deals with a well-defined self-consistent physical problem. According to the nature and the range of the mutual interactions, the building locks are hierarchically organized. A stellar atmosphere model designed along these lines will be an unparalleled benchmark to ascertain the quantitative effects brought about by each individual physical process and allow a direct control onto the partial results obtained within each block. In other words, following these guidelines we will have at hand a veritable Laboratory for Computational Astrophysics that allows us to perform experiments through numerical computation (see [21]). One of the most critical points is a reasonable description of the complex phenomena characterizing radiative transfer processes that regulate the exchange of energy between matter and radiation. The main goal of the present paper is to outline the main features of the iterative sequential procedure that will be the backbone of the in fieri TEIDE 3 package, designed to compute extended stellar atmosphere model by means of a novel sequential iterative procedure and that will be presented in a further paper of the present series. Our operative and physically sound iterative sequential procedure will allow us to treat properly both the fluid dynamics and the radiative transfer equations, whose solution in the actual version is obtained at the price of crude simplifications. The RT equations (in spherical geometry) are solved making use of the implicit integral method (IIM [22]). This method, that does not require the storage and inversion of huge matrices, yields an accurate solution of the latter at a computational cost that will allow for the coupling of the RT equations with those governing the inner stellar layers. We wish to stress that the exceedingly large CPU time for the numerical solution of the radiative transfer equations required by the models or codes available in the literature makes the foregoing coupling unfeasible in the practice. Such drawbacks are overcome by the IIM, which warrant the quickness and flexibility required by an effective modeling of the circumstellar envelopes of AGB stars. As a matter of fact, the preliminary results presented in Section 3.2 has been obtained with a personal computer within a reasonable lapse of time.

The Computational Laboratory
The variety of physical processes listed in the previous section and the complexity of their interactions require a stable and steady approach for this particular case of the stellar atmosphere problem. Numerical experiments performed in the computational laboratory mentioned in Section 1 will point out the way towards achieving an effective solution. As a first step we have worked out firstly a simplified case-study, but one that includes the fundamental physical processes, that will constitute the cornerstone for further more sophisticated models.

Structure and Building Blocks of AGB Stars' Model Atmospheres
Let us consider the extended atmosphere of an AGB stars whose properties do not vary with time. The structure of such a system will be the results of a steady-state equilibrium among a number of physical processes. The gravitational confinement of matter is expressed by the constraint of mechanical equilibrium and the interchange of energy between the stellar material and the radiation field is constrained by the condition of energy conservation. The two foregoing conditions are coupled by the equation of state for matter (either macroscopic or microscopic) that will take into consideration the formation of selected molecules. On the other hand, the fundamental evidence of an outward flux of radiant energy compels us to supplement the energy equilibrium and state equations with the RT equation(s). For simplicity's sake, we will reduce the mechanical equilibrium to hydrostatic equilibrium and neglect convective transport. Moreover, the local supply of a prescribed amount of non-radiative energy may be taken into consideration in order to mimic the dissipation of mechanical energy brought about by, for example, shocks or pulsations. All these physical 'ingredients' can be organized into individual building blocks, whose proper arrangement will reproduce the structure of the physical system. Inside each block a well-defined physical problem can be solved, provided that the values of the relevant input variables from other blocks are known.

The Iterative Sequential Procedure
We may list the fundamental equations that shape the structure of a stellar atmosphere following a 'natural' order: constitutive equations (equations of motion and conservation of both mass and momentum); macroscopic equation of state; atomic occupation numbers (Boltzmann and Saha equations in local thermodynamic equilibrium (LTE), statistical equilibrium equations in non-LTE); transport equations (radiative transfer and convective transport); and conservation of energy for both matter and radiation field.
The direct solution of the resulting system of equations is not feasible, owing to its non-linear and non-local character. On the other hand, we cannot solve the above sequence of equations one by one, because at each step the number of available equations is less than that of the variables involved. The only way to get rid of this otherwise insurmountable drawback is the iterative sequential procedure conceived by Simonneau in the 1980s and applied by himself and coworkers in the computation of stellar atmosphere models. The most recent presentation of the method can be found in [23]. As a necessary preparation for the next subsection, we will recall here just the guidelines of the method.
According to the nature of their mutual interactions, we group the individual physical processes into elementary blocks, such that each of them contains the statement of a selfconsistent physical problem that cannot be further reduced. We then organize the above building blocks sequentially into a mechanical macro-block and an energy macro-block, coupled through a block that accounts for the microscopic state of matter. If an educated guess at one of the fundamentals variable-the distribution of the temperature with depth turns out to be the optimum choice for physical reasons-is introduced as an external input, the equations of the mechanical macro-block can be solved to yield the values of all the input variables required for the successive determination of the microscopic state of matter. Hence, the macroscopic transport coefficients that specify the radiative transfer (RT) equations can be computed. In the next step, within the energy macro-block, the latter are solved and the constraint of energy conservation can be checked. If the constraint is not fulfilled, the current values of the temperature must be up-to-dated using the energy conservation constraint as a transcendental equation for the temperature. The sequential solution of the RT equations followed by the temperature correction is iterated a pre-set number of times (inner iterations) keeping unaltered the transport coefficients that are the input of the energy block. Afterwards, the procedure is repeated from the beginning with the new temperature distribution, and iterated until the energy constraint is satisfied within a prescribed criterion of convergence (outer iterations).

Flowchart of the Iterative Sequential Procedure
We will implement in the TEIDE package the above iterative sequential procedure, whose flowchart is shown in Figure 1. Just a few complementary remarks may be useful for an easier interpretation of the per se self-evident figure.
(1) Because of the spherical symmetry assumed, the natural (discrete) independent variable is the outward radial direction, defined on the external input r-grid. The main external input is the initial temperature distribution T 0 (r) that will be up-to-dated at each stage of the iterative procedure by means of the temperature corrector in the energy block; (2) The fundamental stellar parameters M, R, and T e f f that individualize a star, as well as its chemical composition, constitute a further input to the mechanical block; (3) A first checkpoint can be set in the mechanical lock: as a by-product of the solution of the equation of state we get quantitative information on the presence of the molecules under study (H 2 , H 2 O, CO, SiO, and TiO). Opacity data for molecules are extracted from the EXOMOL database [24]; (4) Within the opacity block we can ascertain the effect of individual molecules on the global absorption coefficient; (5) External input: local supply at a chosen depth of an arbitrary amount of non-radiative energy that mimics the passage of a shock wave; (6) Another valuable checkpoint can be set after the computation of the variable Eddington factors. The latter are defined as the ratio of the K-moment (2nd order angular moment of the specific intensity of the radiation field) to the zero-order J-moment, which yield the closure of the system of bolometric RT equations, whose solution J(r)-necessary for correcting the temperature distribution-automatically satisfies the ER constraint.

A Test Case
In order to check its reliability we applied the forgoing iterative sequential procedure to the test-case of a steady-state and spherically symmetric extended AGB atmosphere under the simplifying hypotheses of hydrostatic equilibrium and local thermodynamic equilibrium (LTE) for matter. The constraint of energy conservation is imposed, either in the case of radiative equilibrium (RE) or when a prescribed arbitrary amount of non-radiative energy is injected at a given depth into the atmosphere. We choose the set of fundamental stellar parameters M = 2.38 × 10 33 g , R = 1.9 × 10 13 cm and T e f f = 3000 K that mimic a low-mass star at the beginning of the AGB phase. The quantities reported above have been extracted from FRUITY full stellar evolutionary models [11] and correspond to an M star (we consider initially only O-rich stars).

The Temperature Correction
The constraint of energy conservation is given by where a ν (r) is the monochromatic absorption coefficient and Q(r) = Q 0 f (r − r 0 ) is the distribution with depth around r 0 of the local supply of non-radiative energy. The dimension of the latter, like that of the integrals, is erg · cm −3 . The function f (r − r 0 ) is a narrow profile peaked at r 0 . If J(r) is given by the quadrature of the J ν terms, the solution of the system of monochromatic RT equations that does not warrant for the fulfillment of the ER constraint, the iterative correction of the temperature distribution via Equation (1) is equivalent to the so-called Λ-iteration, whose drawbacks are well known and widely discussed in the literature on stellar atmosphere theory. An effective approach of addressing the shortcomings of the Λ-iteration is given by the iterative factors method, introduced by Simonneau and Crivellari [25].
If we introduce the standard mean Planck absorption and the mean J absorption, defined as we can define the factor α ≡ a J / a P that allows us to recast Equation (1) into the form The factor α turns out to be quasi-invariant throughout the steps of iteration because it is the ratio of two homologous quantities.
To get rid of the drawbacks of the Λ-iteration, the crucial step is to compute J(r) from the solution of the system of the two first angular moments of the monochromatic RT equation for the specific intensity, integrated over the whole frequency range. These two equations involve the first three moments J, H and K; the system can be closed by mean of the above-mentioned variable Eddington factor f (r) ≡ K(r)/J(r). Eventually one obtains the bolometric equation whose solution automatically satisfies the ER constraint. The second term on the RHS of Equation (5) is a known output of the preceding elementary blocks. The quantity r 2 0 H 0 is the constant amount of radiative energy carried layer by layer and χ H (r) is the flux-averaged extinction coefficient. The latter and the variable Eddington factors are also the ratio of homologous quantities. Together with α they constitute a set of quasi-invariant iteration factors that warrant a fast convergence to the correct physical solution. Equation (4), where J(r) is now the solution of Equation (5), is a transcendental equation in T, whose solution yields the updated temperature distribution.

Preliminary Results
The test case specified at the beginning of this section has been the benchmark for ascertaining the effects on the structure of the outermost layers brought about on the one hand by the inclusion of molecules in the opacity of the stellar material, by a prescribed supply of non-radiative energy on the other. We present in the following a set of preliminary results that must be taken as illustrative examples only. We adopted case by case an ad hoc criterion for stopping the procedure (i.e., the percentage correction of the temperature). In the future an homogeneous criterion shall be fixed, such that accurate numerical results are achieved after a number of external iterations not exceedingly large. Moreover, when the departure from radiative equilibrium is taken into consideration, the choice of the parameter Q 0 is just a rough estimate. We proceeded step by step. To start with we computed a reference model atmosphere in RE, whose opacity takes into account the first five bound-free transitions of the H atom, the bound-free and free-free transitions of the ion H − and the main bound-bound transitions of 30 strong lines that predominate in the spectrum of cool stars. In this exercise the iterative sequential procedure quickly converges (nine outer iterations) to a physically sound solution. The structure of the reference model is shown in panel (a) of Figure 2. The contributions to the opacity due to H 2 O, CO, SiO, and TiO have been successively taken into account. The frequency profile at selected depths of the relevant absorption coefficient a ν is compared in Figure 3 with that of the reference model. The new contributors, mainly CO, H 2 O, and SiO, enhance a ν at wavelengths longer than λ = 5300 Å; moreover the profile is markedly different as a sequence of abrupt maxima and minima show up. CO molecules form at T < 3900 K, those of SiO and H 2 O at T < 2600 K and 2300 K, respectively. As shown in Figure 3, the contributions to opacity by CO and SiO predominate (see the four outstanding peaks at T of order of 2300 K) till the H 2 O molecule takes over at T < 1500 K.
The effects of molecular opacity on the structure of the outer atmosphere (i.e., layers with temperature lower than 5000 K) are visible in the left panels of Figure 2. In panel (a) the temperature profile of the reference model shows a plateau for 13.28 < LogR < 13.40 (between the points marked with 3 and 4 in Figure 2) and a successive mild descent to a second plateau (from point 2). The behavior of the temperature is different when molecular opacity is introduced (see panel (c) of Figure 2): the drop at the end of the plateau (point 3) is here steep; successively (point 2) the temperature decreases steadily towards a much lower surface value. Both behaviors are consistent with theory, which predicts that the temperature distribution of an LTE-RE atmosphere is a monotonic decreasing function outwards, characterized by plateaux where previously dominant transitions become optically thin at increasing radial distance. We can monitor the corresponding depths thanks to the check point (4) of the flowchart. In the reference model the plateau extends from the depth where the b-f transitions of H − become optically thin to that where the escape of photons in the Balmer continuum starts. When the molecules are taken into account the plateau, now shorter, is due to the escape of photons in the H − b-f and Paschen continua; beyond point 2 all the molecules become optically thin.
The foregoing solutions correspond to models in RE. The effects due to a local supply of non-radiative energy (which mimics the passage of shocks within the AGB circumstellar envelope) are shown in the right panels of Figure 2. In both cases, the depth point for the injection of extra energy has been chosen at Log R = 13.34, i.e., a depth corresponding to the temperature plateau, and the parameter Q 0 that fixes the amount of non-radiative energy supplied set equal to 0.025 erg · cm −3 . This quantity corresponds to some 5% of the radiant energy, J a , absorbed at the same depth in the case without molecules, to about the 10% when the molecules are taken into account.
The rate of convergence when the contribution of the molecules to the total opacity is taken into account is lower, even worse when the local supply of non-radiative energy is considered. In order to achieve quickly illustrative results we relaxed the exit criterion on the temperature correction imposed in the reference model and set ∆T/T ≤ 0.50% almost everywhere instead of ∆T/T ≤ 0.25%. We had, however, to pay a price for that: in correspondence to the steep variation of the temperature profile in cases (b), (c), and (d) the departure from energy conservation is higher than that of the reference model, where departures from zero is observed only in the outer layers and is always less than 1.5 %.
Close examination of the results shows a clear correlation between the percent correction of the temperature and the percent departure from energy conservation. The temperature correction is much slower in correspondence to sudden variations of the temperature profile. Consequently, if we stop prematurely the cycle of outer iterations, we cannot achieve a good fulfillment of the energy conservation constraint.

Conclusions
A journey of a ten thousand miles starts beneath one's foot. (Lao-Tsu, 6th century BC) As shown in Section 2, from its very nature the above iterative sequential procedure makes it possible to set up a computational laboratory for ascertaining the effects brought about by different physical conditions. The iterative sequential procedure detailed above is the master tool of such a laboratory. Vis-à-vis the difficulties due to the entanglement of different physical processes, following Lao-Tsu's teaching we started our far-reaching investigations with the test-case specified in Section 3. In spite of the simplifying hypotheses adopted, it takes into account all the essential physics that shapes the structure of the atmosphere of the stars under study.
The preliminary results achieved so far pave the way to further investigations that will provide new insights into the physics at work in AGB stars. The future realization of the TEIDE package and its implementation into the VULCAN code will allow us to handle more realistic physical scenarios. On the one hand in the mechanical block the assumption of hydrostatic equilibrium will be replaced by the equations of motion and the occupation numbers of the atomic populations will be determined by the system of statistical equilibrium equations to remove the actual LTE hypothesis. Moreover, we will take into consideration convective transport whenever brought about by the existing physical conditions, as in [26]. In future work, we shall also implement dust growth in VULCAN using a non-equilibrium kinetic network linking molecules seeds to dust grains. On the other hand, the analysis of the results presented in Section 3.2 highlights some weak points in the numerical treatment of radiative transfer, both the monochromatic RT transfer and the bolometric equation for J(r). At first glance it seems that the blame should be laid at the door of the numerical computation of the first and second derivatives of the source function of each specific RT equation in the energy block, performed by means of a piece-wise parabolic approximation that yielded fairly good results in previous applications. However, our previous experience showed that numerical instabilities can arise when the frequency grid is not fine enough. Therefore, following [27], we will replace the piece-wise parabolic approximation with a cubic spline model for each specific source function. This delicate step requires a careful fine-tuning, but it will improve both the computational accuracy of the solution of the monochromatic RT equations and the successive computation of the coefficients of the bolometric equation for J. These numerical improvements are expected to bring about a better fulfillment of the energy conservation constraints also at the critical depth points pointed out in Section 3.2.
Author Contributions: Conceptualization, L.C., S.C. and L.P.; software, L.C.; validation, L.C., S.C. and L.P. All authors have read and agreed to the published version of the manuscript.

Funding:
We acknowledge the support by the Instituto de Astrofísica de Canarias, the INFN Section of Perugia and the INAF-Osservatorio Astronomico d'Abruzzo.

Conflicts of Interest:
The authors declare no conflict of interest. Notes 1 fruity.oa-abruzzo.inaf.it, accessed on 10 June 2021. 2 The expression 'curse of dimensionality' was coined by R.E. Bellman (Dynamic programming, 1957, Princeton University Press).