| The code SUMA |
4. RELEVANT ISSUES
4.1 Geometrical thickness of the slabs
In a code accounting for shocks the choice of the slab thickness cannot be determined a priori as an input parameter.
The first slab
The choice of the first slab geometrical thickness is the most crucial: the thinner the better, but it is not always true because you must take into account:
(a) the equipartition lifetime between ions and electrons;
(b) the sputtering critical distance against deceleration and stop;
(c) numerical precision.
High temperature region
Once crossed the shock front the gas is thermalized to a temperature
T ~1.5 105 (Vs/100 kms-1)2 K.
At high temperatures (≤106 K) recombination coefficients are very low leading to a low cooling rate. The cooling rate being proportional to n2 (n is the gas density), the thickness of the first slab depends strongly on the density of the gas.
Intermediate temperatures
At T between 104 and 105 K the UV lines and the coronal lines in the IR are strong and lead to rapid cooling and compression of the gas. If the cooling rate is so high to drastically reduce the temperature eluding intermediate ionization-level lines, the calculated spectrum will be wrong (Fig. 2). Therefore, the slab width must be reduced and all the physical quantities recalculated. The choice of the slab width is determined by the gradient of the temperature. This process is iterated until the thickness of the slab is such as to lead to an acceptable gradient of the temperature (T(i-1) - T(i))/T(i-1) = 0.1, where T(i) is the temperature of slab i.
Beyond the temperature drop
In the radiation-bound case, the gas recombines completely before reaching the edge of the nebula opposite to the shock front (if str=0), or the radiation dominated edge (if str=1). Thus, as the temperature drops, a large portion of gas is maintained at a temperature of the order of 104 K by the diffuse (or secondary) radiation. Due to a lower temperature gradient, calculations in this zone proceed throughout wider slabs. This is critical to self-absorption of free-free radiation in the radio domain, in particular for symbiotic stars, where the densities are very high and the nebula sizes relatively small. The optical thickness of the slab τ is actually proportional to n2dr, where dr is the geometrical thickness of the slab.
The radiation dominated edge
In case str = 1 mentioned in the previous section, the edge opposite to the shock front is reached by the photoionizing flux from an external source (AGN center, hot star, starburst etc.). The geometrical thickness of the slabs cannot be as large as those calculated in the shock dominated case (str=1 and U=0) or in the radiation dominated case (str=0 and U≠0). The effect of the radiation flux on the gas physical conditions has to be accounted for in the subsequent iterations (see below), but the width of the photoionized slabs is evaluated in the first iteration, at the start of the calculations. Thus, to improve the precision, it is assumed that the slab width in the radiation dominated edge is a fraction of the total geometrical thickness of the photoionized gas, which is estimated from the ionization parameter and the density. Usually a fraction 1/50 to 1/100 is sufficient to reach a smooth distribution of the physical parameters and a good accuracy for the emission-lines.
4.2 The compression equation
Very different line and continuum spectra result adopting a pure photoionization code which assumes a constant density throughout the nebula, and a code which accounts for shocks, because compression downstream leads to a different distribution of the physical conditions.
After some algebra the Rankine-Hugoniot equations for the conservation of mass and of momentum result in the compression equation [32]:
[(B02 /8π) x2 + n0 kT(1 + ne/n)x - P0]x = -ρ v02 ,
where x=n/n0, P0= ρ v02+B02/8π + n0(1 + f0/1.16)kT, f0 being the initial hydrogen ionization fraction, and ρ0 the initial mass density.
This is a third degree equation which is solved in each slab by a subroutine that finds the roots of a real polynomial using Muller's method of inverse parabolic interpolation, choosing the adequate solution, since some solutions may be imaginary. Compression strongly affects the cooling rate and consequently, the distribution of the physical conditions downstream, as well as that of the elemental fractional abundances.
4.3 The electron temperature
The temperature in each slab depends on energy gains (G) and losses (L) of the gas. Close to the shock front downstream, collisional mechanisms prevail and the temperature is calculated from the energy equation in terms of the enthalpy change [32]. In the slabs with temperature ≤ 2 104 K, photoionization and heating by both the primary and the secondary radiation, dominate and the temperature is calculated by thermal balance (G=L).
Gains
The gains are calculated by the rate at which energy is given to the electrons by the radiation field [38]. The energy of suprathermal electrons created by photoionization is rapidly distributed among the thermal electrons through collisions, heating the gas.
Losses: the cooling rate
Several processes contribute to gas cooling. The cooling components are given by:
L = Lff + Lfb + Llines + Ldust
where Lff corresponds to losses due to free-free transitions, particularly strong at high temperatures and high frequencies, Lfb to free-bound transitions, being high at temperatures near 105 K; Llines is due to line emission mainly between 104 and 105 K.; Ldust represents the energy lost by the gas in collisional heating of dust grains, being higher the higher d/g and agr. Self-absorption of free-free radiation is included in the calculations.
Ionization equilibrium
The ionization equations are solved for each element up to the maximum ionization level, taking into account all the processes previously listed (Sect. 2).
4.4 Lifetimes
Although the code is built for a hydrostatic regime, the most significant lifetimes must be checked:
(a) tev corresponds to the evolution in one slab,
(b) trad to the radiative lifetime,
(c) teq is the equipartition time between ions and electrons in the first slab downstream, and
(d) tsput is the lifetime against sputtering.
They must satisfy the following relations: tev ≥ trad , tev ≥teq and tev ≤ tsput
4.5 Dust
Dust reprocessed radiation and sputtering of the grains (Fig. 3) is calculated consistently with gas emission. A detailed presentation of the equations is given in sect. 3 of [22].
4.6 Iterations
The calculations start at the shock front.
In the infalling case (shock and photoionization on the same edge) one iteration is enough.
In the ejected case (shock and photoionization on opposite edges) the code proceeds through some iterations until convergence of the results (e.g. the distribution of the temperature and of the density throughout the nebula have not changed from the results in the previous iteration ).
I-iteration: the calculations start at the shock front and are performed in a shock dominated case (U=0) in which the primary radiation from the external source is neglected. The diffuse radiation from each slab is stored
II-iteration: the calculations start at the edge reached by the primary radiation. The calculations account for diffuse radiation from both the sides of the nebula. The secondary radiation and the primary radiation fluxes calculated by radiation transfer in each slab in the photoionized zone are stored. In the region close to the shock front the temperature and the density from the I-iteration are used.
III-iteration: the calculations start at the shock front and account for the primary radiation and secondary radiation from the ionized edge which were stored in the previous iteration. The process proceeds by recalculating the physical parameters and the secondary radiation from the slabs heated by the shock front. The next iterations follow the procedures used in II and III. From II iteration on, the results are compared to those obtained in the previous iterations to check for convergence. Generally three iterations are enough to reach the desired accuracy.