A&A 403, 261-275 (2003)
DOI: 10.1051/0004-6361:20030357

Monte Carlo transition probabilities. II.

L. B. Lucy

Astrophysics Group, Blackett Laboratory, Imperial College London, Prince Consort Road, London SW7 2AZ, UK

Received 31 January 2003 / Accepted 6 March 2003

Abstract
The macroscopic quantizations of matter into macro-atoms and radiant and thermal energies into r- and k-energy packets initiated in Paper I is completed with the definition of transition probabilities governing energy flows to and from the thermal pool. The resulting Monte Carlo method is then applied to the problem of computing the hydrogen spectrum of a type II supernova. This test problem is used to demonstrate the scheme's consistency as the number of energy packets ${\cal N} \rightarrow \infty$, to investigate the accuracy of Monte Carlo estimators for radiative rates, and to illustrate the convergence characteristics of the geometry-independent, constrained $\Lambda$-iteration method employed to obtain the NLTE stratifications of temperature and level populations. In addition, the method's potential, when combined with analytic ionization and excitation formulae, for obtaining useful approximate NLTE solutions is emphasized.

Key words: methods: numerical - radiative transfer - stars: atmospheres - stars: supernovae: general - line: formation

1 Introduction

In an earlier paper (Lucy 2002; Paper I), a technique for imposing the constraint of statistical equilibrium on Monte Carlo (MC) transfer codes was developed. This was achieved by first replacing the rate-equation representation of statistical equilibrium by one involving macroscopic energy flows and then quantizing these flows into indivisible energy packets. These e-packets - the MC quanta - contain either radiant energy (r-packets) or kinetic energy (k-packets) and undergo changes on interaction with matter in accordance with certain probabilistic rules. Because these rules are analogous to the transition probabilities that govern a real atom's interactions with photons and electrons (or other colliding species), they can be interpreted as (MC) transition probabilities that govern a macro-atom's interactions with r- and k-packets, with the macro-atom simply representing the atoms of a particular species in a finite volume element.

In addition to providing its theoretical basis, Paper I subjected the macro-atom formalism to various numerical tests, each of which concerned a single point in a stratified atmosphere. In the first group of tests, the consistency of the scheme was confirmed by demonstrating that, when used with the exact level populations, the technique asymptotically reproduces the correct line- and continuum emissivities. But consistency, though necessary, does not guarantee that an effective NLTE code can be developed with this technique. Accordingly, Paper I also tested the sensitivity of the MC emissivities to departures from the exact level populations, since extreme sensitivity would make it unlikely that the scheme would be successful in an iterative search for a NLTE solution. Fortunately, the emissivities were found to be remarkably insensitive to the level populations, thus suggesting good convergence characteristics.

Given these successes with one-point tests, the next step is evidently to apply the macro-atom formalism to a stratified medium. Possible test problems are many and varied since the technique is quite general: potentially, it applies to multi-species plasmas, to multi-dimensional geometries, and to problems with non-radiative heating. But given the technique's novelty, step-by-step testing is advisable, starting with the simplest of problems. Accordingly, this paper treats a 1D medium comprising just one atomic species and subject only to radiative heating. Specifically, the problem is to compute the level populations and kinetic temperature throughout the outer envelope of a type II supernova, which we take to be pure hydrogen. In such extended, low density envelopes, substantial departures from LTE arise due to the strong dilution of continuum radiation and to the negligible importance of collisional excitations. The technique's ability to solve non-trivial NLTE problems is therefore tested. In particular, the convergence behaviour of the iterations required to obtain the level populations and the temperature stratification can be investigated. Moreover, this test serves to illustrate the structure of NLTE codes using the macro-atom formalism.

2 Supernova model

In this section, the basic assumptions concerning the structure of the supernova and the relevant atomic physics are stated. In contrast to previous MC codes for SNe (Lucy 1987; Mazzali & Lucy 1993; Lucy 1999b; Mazzali 2000), the interest here is to test technique rather than to construct a code immediately suitable for analysing observational data. Accordingly, state-of-the-art modelling precision with respect to the SN envelope or to the atomic physics and radiative transfer is of no concern.

2.1 Envelope stratification

In formulating this test problem, we take the computational domain D to be exterior to the layers where the energy released by radioactive decays is thermalized. Then, with the further assumptions that the radiative cooling and recombination time scales are small compared to the expansion time scale, the principles governing spectrum formation are statistical and thermal equilibrium in every local co-moving frame within D.

Because the energy sources are interior to the lower boundary (r = R) of D, this calculation does not predict the SN's luminosity L, which is therefore a parameter. For the same reason, a boundary condition at r = R on the outwardly-directed radiation field is required. This is taken to be

$$I^{+}_{\nu}(R) = B_{\nu}(T_{\rm b})$$ (1)

where $T_{\rm b}$, the black-body temperature, remains to be determined.

In Eq. (1) and throughout this paper, $\nu$ denotes frequency in the matter frame; frequency in the rest frame will be denoted by $\nu_{R}$. Similarly, the energy of an r-packet in the matter and rest frames will be denoted by  $\epsilon_{\nu}$ and  $\epsilon_{R}$, respectively.

The density distribution at time t after the explosion is obtained as usual from the assumption of homologous expansion. Thus

$$\rho(v,t) = \left( \frac{t_{1}}{t} \right)^{3} \rho_{1}(v)$$ (2)

where v = r/t is the constant velocity of the mass layer that is at radius r₁ at time t₁, and  $\rho_{1}(v)$ is the density-velocity profile at elapsed time t₁. The profile adopted is that of Arnett's (1988) explosion model for SN1987A.

With the further assumption of a pure H composition, the atomic number density is $n_{\rm H}(r,t) = \rho(r,t)/m_{\rm H}$.

2.2 Parameters

For the problem thus defined, the basic parameters are t and L. The test calculation of Sect. 7 is carried out at t = 9.8 days, and L(t) is required to be $1.68 \times 10^{41}$ erg s^-1, the value determined observationally for SN 1987A (Suntzeff & Bouchet 1990).

In addition to setting parameters, the location of the lower boundary must be specified. This is taken to be the layer moving at  v₁ = 6500 km s^-1, corresponding to  $R = 5.50 \times 10^{14}$ cm.

2.3 Discretization

The MC calculation is performed for a discretized version of the SN's envelope. Specifically, the envelope is modelled as ${\cal M}$ constant-density spherical shells, within each of which the temperature T and the level populations n_i are also constant.

The independent variable is x = R/r, in terms of which the shells have constant thickness $\Delta x = (x_{1}-x_{{\cal M}+1})/{\cal M}$, with x₁ = 1 and $x_{{\cal M}+1} = 0.25$. The mean radius of the mth shell is defined to be $\overline{r}_{m} = R/ \overline{x}_{m}$, where $(\overline{x}_{m})^{-3} = (x_{m}^{-3}+x_{m+1}^{-3})/2$. The density of the mth shell is then the value given by Eq. (2) at $v=\overline{r}_{m}/t$.

2.4 Atomic model

The model H atom is that used in the one-point consistency test in Sect. 5.2 of Paper I. Thus there are 15 levels, with level 15 being the H⁺ continuum ${\kappa }$. The 14 bound levels correspond to principal quantum numbers i = 1-14 and have consolidated statistical weights g_i = 2i².

Transition probabilities and photoionization cross sections for hydrogen can of course be computed with arbitrary accuracy. Accordingly, the oscillator strengths for bound-bound (b-b) transitions are derived from an accurate look-up table. But the continuous absorption coefficient is simplified by setting Gaunt factors = 1. Thus, the cross section for photoionization from level i is

$$a_{i \kappa}(\nu) = a_{i} \left( \frac{\nu_{i}}{\nu}\right)^{3} \;\; \mbox{for} \;\; \nu > \nu_{i}$$ (3)

where a_i is the cross section at the ionization threshold frequency $\nu_{i}$.

3 Monte Carlo preliminaries

In the following three sections, a detailed account is given of how the internal radiation field in the SN's envelope is derived using the macro-atom formalism of Paper I. This MC calculation is subsequently embedded in an outer iteration loop (Sect. 6.5) in the search for the NLTE solution.

3.1 Random numbers

A random number generator is essential for every MC code. In this investigation, the routine ran2 from Numerical Recipes (Press et al. 1992) is used, but with minor changes to convert to double precision. Independent tests confirm the high quality of sequences of random numbers created with this routine.

In the rest of this paper, the variable z denotes a random number sampling a uniform distribution in the interval (0,1) and thus indicates a call to ran2 in the MC code.

3.2 Formulation

In Paper I, three formulations of the macro-atom technique were presented. Two concerned alternative treatments of stimulated emission; and the third, moving media in the Sobolev limit. Here, in view of the large velocity gradients in a SN envelope, we treat line formation in the Sobolev approximation and so follow the formulation presented in Sect. 4.3 of Paper I.

Although accuracy is not of concern in this paper, it is nevertheless of interest that Duschinger et al. (1995) find that, in treating the formation of hydrogen lines in type II SNe, the Sobolev approximation is as accurate as the comoving frame method.

3.3 Approximations

Because of the envelope's large optical depth in the Lyman continuum, an r-packet in that continuum will propagate a negligible distance before it undergoes a bound-free (b-f) absorption. Accordingly, a useful economy of computational effort is obtained by assuming that such packets are re-absorbed at the point of emission. This is the familiar "on-the-spot'' approximation commonly used for models of photoionized nebulae (e.g., Osterbrock 1974). Moreover, as for such models, this approximation is implemented for the dominant source of r-packets with $\nu > \nu_{1}$ by setting  $\alpha_{1}$ and  $\gamma_{1}$, the ground state recombination and photoionization coefficients, to zero. Of course, Lyman continuum r-packets can in principle also be emitted by f-f emission and by recombinations to excited levels. Such packets if they occur are re-absorbed in situ by activating a macro-atom to state ${\kappa }$.

The above approximation effectively imposes an upper limit  $\nu_{\rm U} = \nu_{1}$ on the MC radiation field. In a second approximation, a lower limit  $\nu_{\rm L}$ is also imposed. The reason for this is the $\nu^{-2}$dependence of f-f absorption as $\nu \rightarrow 0$, which results in occasional jumps in estimates of  ${\cal H}^{\rm ff}$, the f-f heating rate, when an r-packet is emitted with $\lambda \ga 100~\mu$m. The limit  $\nu_{\rm L}$ eliminates such jumps at the price of a slightly biased estimate of  ${\cal H}^{\rm ff}$.

The chosen lower limit $\nu_{\rm L} = \nu_{\rm U}/196$, the ionization threshold of the highest bound level i = 14. To ensure that r-packets are not emitted with $\nu < \nu_{\rm L}$, we set oscillator strengths f_ji = 0 if $\nu_{ji} < \nu_{\rm L}$ and also exclude $\nu < \nu_{\rm L}$ when sampling the f-f emissivity function.

Note that, in order to avoid loss of generality, these devices for imposing frequency limits are ignored in all subsequent formulae.

A further approximation is the neglect in the transfer theory of all terms of O(v/c) except for the Doppler effect (McCrea & Mitra 1936; Lucy 1971, 1999b). But apart from time-delay, it is quite straightforward to include relativistic effects in MC codes (Mazzali & Lucy 1993).

3.4 Initiation

The NLTE solution must be obtained iteratively and so starting values are required for the basic variables. For every iteration except the first, the previous iteration's radiation field $J_{\nu}$ is available, as also are the corresponding temperatures T and level-populations n_i. For the first iteration, initial guesses are required for these quantities.

The initial mean intensity of the radiation field in the Balmer and higher continua is taken to be

$$J_{\nu}(r) = W(r) B_{\nu}(T_{\rm b})$$ (4)

where the dilution factor $W = \frac{1}{2} \left[1-\sqrt{1-x^{2}}\right]$.

The lower-boundary temperature $T_{\rm b}$, which also determines the frequency distribution of r-packets launched into D at r = R - see Eq. (1), is initiated at $T_{\rm b} = 6000$ K.

With the radiation field specified, values of T, n_i and $n_{\rm e}$ in each shell could be obtained by imposing the constraints of statistical and thermal equilibrium as described in Sect. 6. Instead, the simpler and cruder option is followed of assuming an isothermal stratification with T =5000 K and then obtaining the initial n_i and $n_{\rm e}$ from the Saha and Boltzmann formulae.

3.5 Monte Carlo transition probabilities

With the current estimates of $J_{\nu}$, T and n_i, the MC transition probabilities are computed for the next simulation of the radiation field. Following an e-packet's capture by a macro-atom, these probabilities determine stochastically the characteristics of the subsequently-emitted e-packet.

The MC transition probabilities defined in Paper I and in this paper involve ratios of the rates of various radiative and collisional processes. As in Paper I, the formulae for such rates always refer to unit volume.

Adopting the notation of Paper I, we define ${\cal R}_{ij} = R_{ij}+ C_{ij}$ to be the total rate of the transition $i \rightarrow j$, with R_ij and C_ij denoting radiative and collisional rates. Also, except where noted, the summation convention of Paper I is assumed. Thus suffixes $\ell$ and uindicate summations over all levels <i and all levels >i, respectively.

The probability that a macro-atom in state i emits an r-packet and thus de-activates is then

$$p_{i}^{r}= R_{i \ell} (\epsilon_{i} - \epsilon_{\ell}) / {\ca... ...al D}_{i} = ({\cal R}_{i\ell}+{\cal R}_{i{\rm u}})\epsilon_{i}$$ (5)

where $\epsilon_{i}$ is the excitation energy of level i. The corresponding probability for the emission of a k-packet is

$$p_{i}^{k}= C_{i \ell} (\epsilon_{i} - \epsilon_{\ell})/ {\cal D}_{i} .$$ (6)

In addition to the options of emitting an r- or a k-packet, a macro-atom in state i can spontaneously jump to any other state j. Moreover, such internal transitions occur without the absorption or emission of an e-packet. The probabilities of internal transitions to the particular upper or lower states u and $\ell$ are

$$p_{i{\rm u}}={\cal R}_{i{\rm u}}\epsilon_{i}/{\cal D}_{i} \;\... ...;\; p_{i \ell}= {\cal R}_{i \ell}\epsilon_{\ell}/{\cal D}_{i}.$$ (7)

The above options for leaving state i are complete; consequently,

$$p^{r}_{i} + p^{k}_{i} + p_{i{\rm u}} + p_{i \ell} = 1$$ (8)

as may readily be demonstrated.

In order to compute these MC transition probabilities, the rates R_ij and C_ij must be specified.

3.5.1 Collisional rates

The rates of excitation and de-excitation due to collisions by thermal electrons are

$$C_{ji} = q_{ji}n_{j}n_{\rm e} \;\; \mbox{and} \;\; C_{ij} = q_{ij}n_{i}n_{\rm e}.$$ (9)

The rates of collisional ionization and recombination have the same form, except that the latter, being a three-body process, is $\propto$ $n_{\rm e}^{2}$.

For excitations, the coefficient q_jiis computed with van Regemorter's (1962) approximate formula and the inverse coefficient q_ij by detailed balancing. For collisional ionizations, the coefficient  $q_{j \kappa}$ is evaluated with the approximate formula given by Mihalas (1978, Eq. (5-79)), with the inverse coefficient  $q_{\kappa j}$ again obtained by detailed balancing.

3.5.2 Radiative rates

The radiative bound-bound (b-b) rates are computed in the Sobolev approximation as discussed earlier. Thus, with stimulated emissions treated as negative absorptions, the radiative rates between bound level i and a lower level j are

$$R_{ij} = A_{ij}\beta_{ji} n_{i} \;\;\; \mbox{and} \;\;\; R_{ji} = (B_{ji}n_{j}-B_{ij}n_{i})\beta_{ji} J_{ji}^{\rm b} .$$ (10)

Here the mean intensity $J_{ji}^{\rm b}$ refers to the far blue wing of the transition  $j \rightarrow i$, and  $\beta_{ji}$ is the Sobolev escape probability, given by

$$\beta_{ji}= \frac{1}{\tau_{ji}}(1-{\rm e}^{-\tau_{ji}})$$ (11)

where $\tau_{ji}$ is the transition's Sobolev optical depth. In the special case of homologous expansion,

$$\tau_{ji}= (B_{ji}n_{j}-B_{ij}n_{i})\frac{hct}{4\pi}$$ (12)

with t being the elapsed time since the SN exploded.

For the first iteration, he mean intensities $J_{ji}^{\rm b}$ are from Eq. (4). For subsequent iterations, they are evaluated from the MC radiation field - see Sect. 6.2.

The radiative processes coupling bound levels to the continuum are photoionizations and radiative recombinations. The rate coefficient for spontaneous recombinations to level i is (e.g., Mihalas 1978, p. 131)

$$\alpha^{\rm sp}_{i} = 4 \pi \: \Phi_{i \kappa} \int_{\nu_{i... ...frac{2h \nu^{3}}{c^{2}} \; {\rm e}^{-h \nu/kT} \; {\rm d} \nu$$ (13)

with

$$\Phi_{i \kappa}(T)= \frac{n_{i}^{*}}{n_{\kappa}^{*} n_{\rm e}}$$ (14)

where an asterisk indicates the level's LTE population at the given Tand $n_{\rm e}$. The corresponding rate coefficient for stimulated recombinations  $\kappa \rightarrow i$ is

$$\alpha^{\rm st}_{i} = 4 \pi \: \Phi_{i \kappa} \int_{\nu_{i... ...nu)}{h \nu} \: J_{\nu} \: {\rm e}^{-h \nu/kT} \: {\rm d} \nu$$ (15)

and that for photoionizations  $i \rightarrow \kappa$ is

$$\gamma_{i} = 4 \pi \int_{\nu_{i}}^{\infty} \frac{a_{i \kappa}(\nu)}{h \nu} J_{\nu}~ {\rm d} \nu.$$ (16)

In terms of the above, the total recombination coefficient for level i is $\alpha_{i} = \alpha_{i}^{\rm sp}+\alpha_{i}^{\rm st}$.

If we treat stimulated recombinations as negative photoionizations, the radiative rates coupling i and ${\kappa }$ are $R_{\kappa i}= \alpha_{i}^{\rm sp} n_{\kappa} n_{\rm e}$ and $R_{i \kappa}= \gamma_{i} n_{i}- \alpha_{i}^{\rm st} n_{\kappa} n_{\rm e}$. This latter rate may be written as $R_{i \kappa}= \tilde{\gamma}_{i} n_{i}$, where the corrected photionization coefficient

$$\tilde{\gamma}_{i} = 4 \pi \int_{\nu_{i}}^{\infty} \frac{\tilde{a}_{i \kappa}(\nu)}{h \nu} J_{\nu} ~{\rm d} \nu$$ (17)

is expressed in terms an absorption cross section

$$\tilde{a}_{i \kappa}(\nu) = a_{i \kappa}(\nu) \left(1 - \fr... ... \frac{n_{i}^{*}}{n_{\kappa}^{*}} {\rm e}^{-h \nu/kT}\right)$$ (18)

which incorporates the NLTE correction factor for stimulated emissions.

3.6 Modified rate coefficients

The rates at which b-f and f-b transitions absorb and emit radiant energy are needed in the course of the calculation. These can conveniently be expressed in terms of modified photoionization and recombination coefficients. Thus the rate at which photoionizations from level i remove energy from the radiation field is $\gamma_{i}^{\rm E} \: h \nu_{i} \: n_{i}$, where

$$\gamma_{i}^{\rm E} = 4 \pi \int_{\nu_{i}}^{\infty} \frac{a_{i \kappa}(\nu)}{h \nu_{i}} J_{\nu} ~{\rm d} \nu.$$ (19)

Comparison with Eq. (16) shows that this modified coefficient is obtained from the conventional one by changing $h \nu$ to $h \nu_{i}$ in the integrand's denominator.

Similarly, the corresponding rate at which spontaneous recombinations to level i add energy to the radiation field is  $\alpha_{i}^{\rm E ,sp} \: h \nu_{i} \: n_{\kappa} n_{\rm e}$, where

$$\alpha_{i}^{\rm E ,sp} = 4 \pi \: \Phi_{i \kappa} \: \int_{... ...rac{2h \nu^{3}}{c^{2}} \; {\rm e}^{-h \nu/kT} \; {\rm d} \nu.$$ (20)

This modified coefficient is again obtained from the conventional one (Eq. (13)) by changing $h \nu$ by $h \nu_{i}$ in the integrand's denominator.

Following the same procedure, the modified coefficients  $\alpha_{i}^{\rm E, st}$ and  $\tilde{\gamma}_{i}^{\rm E}$ are also defined. Thus stimulated recombinations to level i emit radiant energy at the rate  $\alpha_{i}^{\rm E, st} \: h \nu_{i} \: n_{\kappa} n_{\rm e}$. But when treated as negative photoionizations, these stimulated recombinations reduce the rate at which photoionizations from level i absorb radiant energy to $\tilde{\gamma}_{i}^{\rm E} \: h \nu_{i} \: n_{i}$.

In terms of the above, the total rate at which recombinations to level i emit radiant energy is  $\alpha_{i}^{\rm E} \: h \nu_{i} \: n_{\kappa} n_{\rm e}$, where $\alpha_{i}^{\rm E} = \alpha_{i}^{\rm E ,sp} + \alpha_{i}^{\rm E ,st}$.

For the first iteration, the coefficients $\gamma_{i}$, $\gamma_{i}^{\rm E}$, $\alpha_{i}^{\rm st}$ and  $\alpha_{i}^{\rm E, st}$ are derived by numerical integration with $J_{\nu}$ from Eq. (4). For subsequent iterations, they are evaluated using the MC radiation field - see Sect. 6.2. In terms of these coefficients, the corrected photionization coefficient $\tilde{\gamma}_{i}= \gamma_{i} - \alpha_{i}^{\rm st} n_{\kappa} n_{\rm e}/n_{i}$and correspondingly for $\tilde{\gamma}^{\rm E}_{i}$

4 Emission and absorption of energy packets

This MC technique directly models the flow of energy through the domain D by computing the interaction histories of numerous indivisible e-packets. In consequence of these interactions, an e-packet is sometimes an r-packet and at other times a k-packet. Accordingly, in this section, the probabilistic rules governing the emission and absorption of both r- and k-packets are developed.

4.1 Initial creation of r-packets

In the absence of non-radiative heating in D, the desired solution satisfies the constraint of radiative equilibrium in the co-moving frame at every point in D. This constraint is rigorously incorporated into the macro-atom formalism by not allowing active macro-atoms to appear or disappear spontaneously within D. This then implies that every e-packet's interaction history both begins and ends as an r-packet crossing a boundary of D (Sect. 7.1, Paper I). Accordingly, in the context of a spherically-symmetrical SN envelope, the creation of r-packets corresponds to a sampling of $I^{+}_{\nu}(R)$ given by Eq. (1).

In the previous SN and stellar wind codes, the initial energies  $\epsilon_{0}$ of the packets were equal, independent of frequency. But now, in order that the weak radiation field between L$\alpha$ and the Lyman limit is well represented, r-packets are created in equal numbers in equal intervals of $\ell n \: \nu$, with energies assigned in accordance with black body emission at $T_{\rm b}$.

Despite r-packets no longer having equal initial energies, MC estimators for integrals of the radiation field can be cast in familiar form if we take  $\epsilon_{0}$ to be the unit adopted for measuring the energies of the packets. A convenient choice is defined by the equation

$${\cal N} \frac{\epsilon_{0}}{\Delta t} = 4 \pi R^{2} \times \sigma T_{\rm b}^{4}$$ (21)

where $\Delta t$ is the interval simulated by the MC experiment and ${\cal N}$ is the number of r-packets launched across the lower boundary r = R in time $\Delta t$.

Now consider an r-packet representing emission by the lower boundary in the interval  $\Delta \ell n \nu$ at $\nu$ during time $\Delta t$. Since the flux transported by  $I^{+}_{\nu}(R)$ in this interval is $\pi \nu B_{\nu}(T_{\rm b}) \Delta \ell n \nu$, this packet's energy in the co-moving frame is given by

$$\frac{\epsilon_{\nu}}{\epsilon_{0}} = {\cal N} \: \frac{\nu... ...T_{\rm b}) \Delta \ell n \nu }{\sigma T_{\rm b}^{4}/\pi} \cdot$$ (22)

A total of ${\cal N}$ packets with fractional energies given by Eq. (22) are created uniformly in $\ell n \nu$ in the frequency interval $(\nu_{\rm L},\nu_{\rm U})$- see Sect. 3.3.

Notice that the values of $\ell n \nu$ are not obtained by random sampling. This is the one point in the calculation where stratified sampling can readily and safely be used to limit MC noise.

The r-packets must also have their initial direction cosines $\mu$ specified. Assuming zero limb-darkening for emission at the lower boundary - see Eq. (1), we randomly sample this limb-darkening law by taking $\mu = \sqrt{z}$. The packet's initial rest energy is then $\epsilon_{R} = \epsilon_{\nu}/(1- \mu v_{1}/c)$.

4.2 Emission of r-packets

A macro-atom can be activated by a k-packet or an r-packet. The MC transition probabilities defined in Sect. 3.5 are then applied until the macro-atom de-activates. This occurs with the emission of a k-packet or an r-packet. If an r-packet is emitted, the next step is to assign $\nu$. The relevant procedure depends on whether de-activation occurred from the continuum state or from a discrete state.

4.2.1 De-activation from a discrete state ${1 < i < \kappa }$

In this case, the emitted r-packet belongs to an emission line arising from atomic level i. One of these lines must therefore be selected randomly, with due weight given to the lines' emissivities. Since the total emissivity of the lines  $i \rightarrow \ell$ is  $\dot{E}_{i}^{R} = R_{i \ell} (\epsilon_{i}-\epsilon_{\ell})$, where the  $R_{i \ell}$ are from Eq. (10), the selection probability for the line $i \rightarrow j$ is $R_{ij} (\epsilon_{i}-\epsilon_{j})/\dot{E}_{i}^{R}$.

If the line thus randomly selected is $i \rightarrow j$, the $\nu$of the emitted r-packet is simply the transition's rest frequency $\nu_{ji} = (\epsilon_{i}-\epsilon_{j})/h$. Note that because we are using the Sobolev approximation in the narrow-line limit, no sampling of an emission profile is required.

4.2.2 De-activation from the continuum state ${\kappa }$

In this case, the emitted r-packet's $\nu$ is obtained by sampling the energy distributions of the spontaneous recombination continua (Sect. 4.2.3, Paper I). This can be done by first selecting a continuum and then selecting a $\nu$ in that continuum.

The emissivity due to spontaneous recombinations to level i is (cf. Eq. (13))

$$4 \pi j_{i}(\nu) = 4 \pi \: \Phi_{i \kappa} \: a_{i \kappa}(... ...u^{3}}{c^{2}} \; {\rm e}^{-h \nu/kT} \:n_{\kappa}n_{\rm e} .$$ (23)

Substituting for $a_{i \kappa}(\nu)$ from Eq. (3) and integrating from $\nu_{i}$ to $\infty$, we find that the integrated emissivity of the  $\kappa \rightarrow i$ continuum is

$$\dot{E}_{\kappa i} = 4 \pi \: \Phi_{i \kappa} \: a_{i} \; \f... ... \frac{kT}{h} {\rm e}^{-h \nu_{i}/kT} \:n_{\kappa}n_{\rm e}.$$ (24)

The selection probability for emission in the $\kappa \rightarrow i$ continuum is therefore $\dot{E}_{\kappa i}/\sum_{i}\dot{E}_{\kappa i}$.

If the selected continuum is $\kappa \rightarrow i$, the next step is to sample its energy distribution. The emitted r-packet's $\nu$ is thus the solution of the equation

$$\int_{\nu}^{\infty} 4 \pi j_{i}(\nu) \: {\rm d} \nu = z \: \dot{E}_{\kappa i}.$$ (25)

Substitution from Eqs. (23) and (24) reduces this equation to the simple formula

$$\nu= \nu_{i} \: \left(1 - \frac{kT}{h \nu_{i}} \:\ell n\: z \right).$$ (26)

Note that r-packets emitted by macro-atoms de-activating from state ${\kappa }$ represent only the loss of ionization energy - see Sect. 5.2 and Fig. 4 in Paper I. But recombination emission also results in the loss of thermal energy. This loss is represented by (some of) the conversions of k-packets into r-packets - see Sect. 4.3.2 below.

4.3 Creation and elimination of k-packets

The MC transition probabilities defined in Paper I model the effect of atomic levels on the frequency redistribution of absorbed radiant energy. Another set of probabilities must now be defined to model the corresponding effect of interactions with the thermal pool. In this formalism, such interactions are represented by the emission and absorption of k-packets.

Note that in developing probabilistic rules for k-packets, we assume that b-f and f-f absorptions as well as collisional de-excitations and recombinations are all followed by an instantaneous re-adjustment of the thermal motions to a Maxwell velocity distribution. Accordingly, when considering the elimination of k-packets, we use rate coefficients that assume a Maxwell distribution for the thermal electrons.

4.3.1 Creation of k-packets

The interactions creating k-packets correspond to the mechanisms heating the gas.

Collisional de-excitations and recombinations convert excitation and ionization energy into heat. These processes are represented in this formalism as ${\cal A}^{*} \rightarrow k$ - i.e., by the emission of a k-packet by an active macro-atom - see Eq. (6).

A second heating mechanism is the conversion of radiant energy into heat by f-f absorptions. In this formalism, this corresponds to  $r \rightarrow k$ - i.e., to the direct conversion of an r-packet into a k-packet without the mediation of a macro-atom.

A third heating mechanism is the conversion of radiant energy into heat by b-f absorptions. For this mechanism, an extended discussion is necessary since not all absorbed energy appears as heat.

In terms of the modified photoionization coefficient defined in Sect. 3.6, the rate at which b-f transitions from level iabsorb radiant energy is $\tilde{\gamma}_{i}^{\rm E} h \nu_{i} n_{i}$. But each photoionization raises the ionization energy by $\epsilon_{\kappa}-\epsilon_{i} = h \nu_{i}$. Accordingly, from the total absorbed rate, the amount $\tilde{\gamma_{i}} h \nu_{i} n_{i}$is used in ionizing the atoms, with only the remainder $(\tilde{\gamma}_{i}^{\rm E}-\tilde{\gamma_{i}}) h \nu_{i} n_{i}$available to heat the gas.

The existence of these two channels might suggest that an r-packet undergoing a b-f absorption should be split. But this would lose the advantages of working with indivisible packets (Paper I, Sect. 1). Instead, therefore, the absorbed energy is assigned randomly to one or other of the channels. Thus, with probability

$$\pi_{i}^{\kappa} = \tilde{\gamma}_{i}/ \tilde{\gamma}_{i}^{\rm E}$$ (27)

the r-packet activates a macro-atom to its continuum state ${\kappa }$ - i.e., $r \rightarrow {\cal A}^{*}_{\kappa}$. Alternatively, with probability $1 - \pi_{i}^{\kappa}$, the r-packet converts directly into a k-packet - i.e., $r \rightarrow k$ . In the first case, all of the r-packet's energy is converted into ionization energy; in the second, all is contributed to the thermal pool.

Summing over all bound states, we find that the probability that a b-f absorption of an r-packet with frequency $\nu$ results in the creation of a k-packet is

$$\pi^{k}(\nu) = 1-\sum_{i} \; p^{\rm bf}_{i}(\nu) \:\pi_{i}^{\kappa}$$ (28)

where $p^{\rm bf}_{i}$ is the probability that the b-f absorption occurred from level i. Now, since the macroscopic b-f absorption coefficient  $k^{\rm bf}_{\nu}$ is given by

$$k^{\rm bf}_{\nu} \rho = \sum_{\nu_{i} < \nu} n_{i} \: \tilde{a}_{i \kappa}(\nu)$$ (29)

where the corrected cross section $\tilde{a}_{i \kappa}$ is defined by Eq. (18), the required probability  $p^{\rm bf}_{i}$ is evidently

$$p^{\rm bf}_{i}(\nu) = n_{i} \: \tilde{a}_{i \kappa}(\nu) \; / \: k^{\rm bf}_{\nu} \rho.$$ (30)

The probabilities defined in this subsection are $\geq$0 if $\tilde{a}_{i \kappa}(\nu) \geq 0$ for all $\nu$. From Eq. (18), we see that this inequality is violated as $\nu \rightarrow 0$ if  $b_{\kappa}/b_{i} > 1$, where b_i = n_i/n_i^* is the departure coefficient for level i. In the numerical tests of Sect. 7, minor violations of the condition $b_{i}/b_{\kappa} \geq 1$ occur for high levels of the model H atom, and these are probably due to the neglect of the recombinations into levels i > 14 and therefore of the subsequent cascade contributions to the included levels. To avoid negative probabilities when this happens, we set $b_{\kappa}/b_{i} = 1$ in Eq. (18) if  $b_{\kappa}/b_{i} > 1$.

4.3.2 Elimination of k-packets

On the assumption that advection and conduction of thermal energy are negligible compared with radiative transport (but note the possibilities of generalization), a k-packet is converted in situ back into an r-packet, either directly ( $k \rightarrow r$) or indirectly ( $k \rightarrow {\cal A}^{*} \rightarrow \cdots \rightarrow r$). The interactions eliminating k-packets correspond to the mechanisms cooling the gas, and their selection probabilities are determined by the mechanisms' cooling rates.

Collisional excitations and ionizations convert heat into excitation and ionization energy. The total cooling rate provided by this mechanism is

$${\cal C}^{{\rm c} \ell} = \sum_{i = 2}^{\kappa} \: {\cal C}^{... ...j = 1}^{i-1} C_{ji}. \left(\epsilon_{i}-\epsilon_{j}\right) .$$ (31)

A second cooling mechanism is the conversion of heat into radiant energy by f-f emission. The cooling rate in this case (e.g., Osterbrock 1974, p. 44) is

$${\cal C}^{\rm ff} = C\; T^{1/2}\; n_{\kappa}n_{\rm e}.$$ (32)

With the mean Gaunt factor set =1, the coefficient $C = 1.426 \times 10^{-27}$ (cgs).

A third cooling mechanism is the conversion of heat into radiant energy by f-b emission. The rate at which spontaneous recombinations $\kappa \rightarrow i$convert thermal and ionization energy into radiant energy is $\alpha^{\rm E, sp}_{i}\: h \nu_{i}\: n_{\kappa} n_{\rm e}$ - see Sect. 3.6. But each such recombination releases $h \nu_{i}$ of ionization energy, which is therefore being radiated at the rate $\alpha^{\rm sp}_{i}\; h \nu_{i}\: n_{\kappa} n_{\rm e}$. Subtracting this from the previous expression and summing over bound states, we find that the rate at which heat is radiated by spontaneous f-b transitions is

$${\cal C}^{\rm fb, sp} = n_{\kappa}n_{\rm e} \sum_{i} \left(\alpha^{\rm E, sp}_{i}-\alpha^{\rm sp}_{i}\right) \; h \nu_{i}.$$ (33)

With cooling rates thus specified, the selection probabilities are as follows: a k-packet is removed by a f-b transition ( $k \rightarrow r$) with probability

$$\pi^{\rm fb} = {\cal C}^{\rm fb, sp}/ \left({\cal C}^{\rm fb, sp} + {\cal C}^{\rm ff} + {\cal C}^{{\rm c} \ell}\right)$$ (34)

by a f-f process ( $k \rightarrow r$) with probability

$$\pi^{\rm ff} = {\cal C}^{\rm ff}/ \left({\cal C}^{\rm fb, sp} + {\cal C}^{\rm ff} + {\cal C}^{{\rm c} \ell}\right)$$ (35)

and by a collisional process with probability $1-\pi^{\rm fb}-\pi^{\rm ff}$. In this last case, the k-packet activates a macro-atom ( $k \rightarrow {\cal A}^{*}$) and so is not directly converted into an r-packet.

Note that the quantity subtracted in Eq. (33), namely the rate of loss of ionization energy, is represented in this MC procedure by macro-atoms de-activating from state ${\kappa }$ (Sect. 4.2.2). Also cooling by stimulated recombinations does not appear in the above formulae because these recombinations, treated as negative photoionization, are allowed for when creating k-packets (Sect. 4.3.1).

5 Monte Carlo calculation

The MC model of the radiation field in the SN envelope is built up from the trajectories of ${\cal N}$ individual r-packets, starting from their launch into D at r=R and ending with their exit from D, either by escaping to $\infty$ or by re-crossing the lower boundary.

As an r-packet propagates in D, it will typically undergo numerous geometrical and physical "events''. The geometrical events are simply crossings of the boundaries of the ${\cal M}$ constant-density shells. But the more important physical events constitute the MC technique's treatment of the interactions of radiation and matter.

In describing how the packets' trajectories are computed, it suffices to explain how to find the next event along one packet's trajectory. Thus, we now suppose that an r-packet has just undergone event n in its propagation history and that this occurred in shell m. The task now is first to find the location and nature of event n+1 and then to compute the frequency $\nu$ and direction of propagation of the r-packet that emerges from this event.

5.1 Absorption and scattering coefficients

Immediately following event n, the macroscopic coefficients of continuum absorption ( $k_{\nu}^{\rm bf}$, $k_{\nu}^{\rm ff}$) and scattering ($\sigma$) must be computed since the r-packet is either in a new shell or has changed its $\nu$.

The b-f coefficient, corrected for stimulated emission, is given by Eq. (18), where the n_i are the current estimates of the NLTE level populations in shell m. The f-f coefficient  $k_{\nu}^{\rm ff}$ is calculated with the standard formula but with mean Gaunt factor = 1, in accordance with the corresponding simplification for Eq. (3).

5.2 Pathlength

If s is distance from event n measured along the packet's trajectory and $\tau(s)$ is the corresponding optical depth, then the pathlength $\Delta s$ between events n and n+1 is determined by the equation

$$\tau(\Delta s) = \Delta \tau$$ (36)

where the beam's exponential decay with optical depth is randomly sampled by taking $\Delta \tau = - \ell n (z)$.

Continuum and line processes both contribute to $\tau(s)$. The continuum contibution is approximately

$$\tau_{\nu}(s) = \chi_{\nu} \rho \: s \;\;\;\;\; \mbox{with} \... ...\; \chi_{\nu} = k_{\nu}^{\rm bf} + k_{\nu}^{\rm ff} + \sigma$$ (37)

where $\nu$ is the packet's frequency immediately following event n. Notice that the small changes in the extinction coefficient seen by the packet as its $\nu$ decreases along the trajectory between events n and n+1 are neglected. This is justified since these changes  $\rightarrow 0$ as ${\cal M} \rightarrow \infty$.

The contribution of lines to $\tau(s)$ is computed with the Sobolev approximation in the narrow-line limit, according to which a line with Sobolev optical depth $\tau_{k}$ attenuates an incident beam by the factor e $^{-\tau_{k}}$ at the point where resonance occurs. Thus $\tau(s)$ undergoes discontinuous increments $\tau_{k}$ at the points s_k where the packet's frequency $\nu(s)$ equals one of the rest frequencies $\nu_{k}$ of the b-b transitions. Accordingly,

$$\tau(s) = \tau_{\nu}(s) + \sum_{k} \tau_{k}$$ (38)

where the summation extends over all s_k in the interval (0,s). Moreover, at each s_k,

$$\tau(s_{k}^{+}) = \tau(s_{k}^{-}) + \tau_{k}.$$ (39)

The above prescription for locating event n+1 applies only if it occurs within shell m. If Eq. (36) is not satisfied in shell m, event n+1 is geometrical, and the pathlength $\Delta s$ is then simply the distance from event n to the point where the trajectory crosses into a neighbouring shell.

5.3 Selecting the physical event

If event n+1 occurs within shell m then it is physical, due to either continuum extinction or a line absorption. The former is indicated if $\Delta s \neq s_{k}$; the latter otherwise.

If the event is a continuum extinction, we decide between the three contributing mechanisms using a single random number: the event is an electron scattering if $z < \sigma/\chi_{\nu}$; failing this, the event is a b-f absorption if $z < (\sigma + k_{\nu}^{\rm bf})/\chi_{\nu}$; failing this, the event is a f-f absorption.

Event n+1 is absorption in line k if $\Delta s = s_{k}$. But since $\tau(s)$ has a discontinuity at s_k, this solution of Eq. (36) is actually recognized by finding that $\tau(s_{k}^{-}) < \Delta \tau < \tau(s_{k}^{+})$.

5.4 New frequency

With the nature of event n+1 determined, the next task is to assign a frequency to the emitted r-packet.

If event n+1 is geometrical, the r-packet is simply continuing its existing flight path into the new shell, and so the $\nu$ at the beginning of the "new'' trajectory is identical to the $\nu$ at the end of the "previous'' trajectory, namely the value at the boundary crossing.

On the other hand, if event n+1 is physical, the nature of the event determines the procedure to be followed in assigning $\nu$.

5.4.1 Electron scattering

With electron scatterings treated as coherent in the comoving frame, the scattered r-packet's $\nu$ is identical to that of the incident packet.

5.4.2 Bound-bound absorption

If event n+1 is absorption by the transition $\ell \rightarrow u$, the immediate result is a macro-atom in state u - i.e. $r \rightarrow {\cal A}^{*}_{u}$. The next step therefore is the repeated application of the transition probabilities of Sect. 3.5 until the macro-atom de-activates. This was discussed at length in Paper I and specifically illustrated in Sect. 5.5 of that paper with the 15-level H atom used here.

The macro-atom de-activates with the emission of either a k-packet or an r-packet. In the former case, the k-packet is eliminated in situ as described in Sect. 4.3.2.

If the macro-atom de-activates with the emission of an r-packet, selection of $\nu$ depends on whether it was emitted from the continuum state ${\kappa }$ or from one of the bound states ${1 < i < \kappa }$. These cases were both discussed in Sect. 4.2, where the selection procedures followed in assigning $\nu$ are given.

5.4.3 Free-free absorption

If event n+1 was a f-f absorption, the immediate result is a k-packet which then in situ either converts into an r-packet or excites a macro-atom. The conversion $k \rightarrow r$ occurs by f-b emission with probability  $\pi^{\rm fb}$ or by f-f emission with probability $\pi^{\rm ff}$ (Sect. 4.3.2). For f-b emission, the selection of $\nu$ for the emergent r-packet is the two step procedure described in Sect. 4.2.2.

For f-f emission, $\nu$ is selected by sampling this mechanism's emissivity function - i.e., by solving the equation

$$\int_{\nu}^{\infty} j_{\nu}^{\rm ff} {\rm d} \nu = z \: \int_{0}^{\infty} j_{\nu}^{\rm ff} {\rm d} \nu .$$ (40)

If the $\nu$-dependence of the velocity-averaged Gaunt factor is neglected, $j_{\nu}^{\rm ff} \propto {\rm e}^{-h \nu/kT}$ and the solution is

$$\nu = -\frac{kT}{h} \: \ell n \: z\cdot$$ (41)

Besides these two direct $k \rightarrow r$ conversions, there is alternatively the indirect conversion $k \rightarrow {\cal A}^{*} \rightarrow \cdots \rightarrow r$, which occurs with probability $1-\pi^{\rm fb}-\pi^{\rm ff}$ and is due to collisional cooling. In this case, the macro-atom is activated to state i with probability ${\cal C}_{i}^{{\rm c} \ell}/ {\cal C}^{{\rm c} \ell}$ - see Eq. (31). Application of the MC transition probabilities as in Sect. 5.4.2 then leads eventually to the emission of an r-packet.

5.4.4 Bound-free absorption

If event n+1 is a b-f absorption, the immediate result is either a k-packet (Sect. 4.3.1) or a macro-atom in state ${\kappa }$. The selection probability is  $\pi^{k}(\nu)$ for the former and $1-\pi^{k}(\nu)$ for the latter, where $\nu$ is the frequency of the incident r-packet.

If this selection results in a macro-atom in state ${\kappa }$, the MC transition rules are applied as described in Sect. 5.4.2 until an r-packet is emitted. On the other hand, if the result is a k-packet, it converts in situ into an r-packet as described in Sect. 5.4.3.

5.5 New direction cosine

Re-emission following a physical event is assumed to occur isotropically in the matter frame, so the new direction cosine in every case is $\mu = 1-2z$. For electron scattering, this is an approximation. But isotropic emission is exact for thermal processes - i.e., for f-f and f-b emissions. It is also exact for the Sobolev treatment of b-b emission in a homologously expanding flow since there is then no kinematically-preferred direction.

If this technique is applied to a spherically-symmetric stellar wind, the expansion is no longer homologous and so b-b emission is angle dependent. Specifically, the pdf to be sampled for the new direction cosine is $\propto$ $I_{\nu_{k}}^{r}(\mu)/\left[\mu^{2} {\rm d}v/{\rm d}r + (1 - \mu^{2}) v/r\right]$, where $I_{\nu_{k}}^{r}(\mu)$ is the intensity in the line's far-red wing of the radiation emitted by the line (Castor 1970; Lucy 1971).

5.6 New rest energy

If $\epsilon_{R}$ is the packet's rest energy during its free flight prior to event n+1, then, at the event, $\epsilon_{\nu} = \epsilon_{R} (1 - \mu_{1} v/c)$, where v is the velocity at that location and $\mu_{1}$ is the incident direction cosine. Now, at every event, the incident and emergent r-packets have the same energy  $\epsilon_{\nu}$. Accordingly, the packet's rest energy along its free flight following event n+1is $\epsilon_{R} =\epsilon_{\nu}/(1 - \mu v/c)$, where $\mu$ is the new direction cosine.

With the assignment of $\nu$ and $\mu$ and the updating of  $\epsilon_{R}$, the treatment of event n+1 is complete. The calculation therefore now returns to Sect. 5.1 to initiate the search for event n+2.

6 Improved solution

The trajectories of ${\cal N}$ individual r-packets computed as described in Sect. 5 collectively constitute the MC estimate of the radiation field in D. The next challenge is to use it to derive improved values of n_i and T for each of the ${\cal M}$ constant-density shells.

6.1 Improved scale factor and boundary temperature

The MC estimate for the luminosity of the SN is

$$L(\infty) = \frac{\epsilon_{0}}{\Delta t} \sum \frac{\epsilon_{R}}{\epsilon_{0}}$$ (42)

where the summation is over the r-packets that escape to $\infty$.

From this equation, we derive an improved scale factor  $\epsilon_{0}/\Delta t$ by requiring that $L(\infty)$ has its specified value - see Sect. 2.2. The improved value of $T_{\rm b}$ then follows from Eq. (21).

This improved scale factor is used below in computing the radiative rates needed for the statistical and thermal equilibrium calculations. The improved $T_{\rm b}$ is used to determine the initial frequency distribution of r-packets (Sect. 4.1) at the start of the next iteration.

6.2 Monte Carlo estimators

In order to use the MC radiation field to improve the values of n_i and T, estimators are needed for the radiative rate coefficients. Fortunately, efficient estimators can be constructed with a procedure given in an earlier paper (Lucy 1999a).

In the context now of a moving atmosphere, the energy-density argument of Lucy (1999a) applied to a local co-moving frame yields the formula

$$J_{\nu}\: {\rm d} \nu = \frac{1}{4 \pi} \frac{\epsilon_{0}}{\... ... \: \sum_{d \nu} \frac{\epsilon_{\nu}}{\epsilon_{0}} {\rm d}s$$ (43)

where the summation is over all trajectory segments ds in V such that the r-packets' frequencies $\nu(s)$ fall in the interval $(\nu, \nu + {\rm d} \nu)$. This equation is the building block that allows us to construct the required estimators.

For this problem, V in Eq. (43) is the volume of one of the spherical shells into which the envelope is discretized (Sect. 3.2). But the estimators derived in this section are quite general: no particular geometry or symmetry is assumed.

If we use Eq. (43) to derive an estimator for the photoionization coefficient given by Eq. (16), the result is a summation, with the summed quantity being the integral of $\epsilon_{\nu} a_{i \kappa}(\nu) / h \nu$over the pathlength $\Delta s$ between consecutive events. But since every pathlength is confined to a single shell, $\Delta s \rightarrow 0$ as ${\cal M} \rightarrow \infty$ and so $\Delta s$ may be treated as a small quantity. This then implies that each $\Delta s$-integral can be approximated by a one-point quadrature. With this approximation, the estimator for the photoionization coefficient for level i is

$$\gamma_{i} = \frac{\epsilon_{0}}{\Delta t} \: \frac{1}{V} \:... ...\epsilon_{0}} \: \Delta s \: \frac{a_{i \kappa}(\nu)}{h \nu}$$ (44)

where $\nu$ is the packet's frequency at the mid-point of $\Delta s$ and the summation is over all r-packets in V with $\nu > \nu_{i}$.

The efficiency of this estimator derives from not being restricted to r-packets that actually cause photoionizations. Thus its precision greatly exceeds that of an estimator based simply on counting photoionizations in V. In particular, this estimator returns a non-zero estimate of  $\gamma_{i}$ even if no r-packets cause photoionizations from level iin a given shell, as commonly happens for the higher, less populated levels and in the outer, less dense shells.

Applying the same procedure now to Eq. (15), we obtain

$$\alpha^{\rm st}_{i} = \Phi_{i \kappa} \: \frac{\epsilon_{0}}{... ...lta s \: \frac{a_{i \kappa}(\nu)}{h \nu} \: e^{-h \nu / kT}$$ (45)

as the corresponding estimator for the rate coefficient of stimulated recombinations to level i, with $\nu$ again evaluated at the mid-point of $\Delta s$. This estimator's efficiency is a consequence of not being limited to r-packets that actually cause stimulated recombinations.

The above estimators refer to transitions to and from the continuum. But statistical equilibrium also depends on the radiative rates of b-b transitions - see Eq. (10). Thus we require an estimator for  $J_{ji}^{\rm b}$, the mean intensity at frequency  $\nu_{ji}^{+}$- i.e. in the far blue wing of the transition $j \rightarrow i$. This can also be derived from Eq. (43).

For homologous expansion, there is no kinematically-preferred direction, and so the rate at which the frequency $\nu(s)$ of an r-packet in free flight decreases with distance s is independent of its direction of propagation. To the accuracy of this calculation - see Sect. 3.3, this rate is ${\rm d} \nu$/ ${\rm d}s =-\nu_{R}$/ct. This equation gives the length ds of the trajectory segment within which $\nu$ is in the interval $(\nu_{ji}, \nu_{ji} + {\rm d} \nu)$. Substituting this value of ds into Eq. (43), we obtain

$$J_{ji}^{\rm b} = \frac{\lambda_{ji} t}{4 \pi} \: \frac{\epsi... ...\frac{\epsilon_{ji}}{\epsilon_{0}} \: \frac{\nu_{ji}}{\nu_{R}}$$ (46)

as an estimator for the mean intensity at  $\nu_{ji}^{+}$. In this formula, $\epsilon_{ji}$ is the packet's energy at the point where $\nu = \nu_{ji}$, and the summation is over all r-packets in V that come into resonance with the transition  $j \rightarrow i$.

The efficiency of this estimator derives from it not being restricted to r-packets that actually undergo the b-b transition  $j \rightarrow i$. Accordingly, a non-zero estimate is usually returned even when no such transitions occur in V. Evidently, the estimate is zero only if no r-packets in V happened to come into resonance during the MC calculation. To avoid such spurious zero estimates for  $J_{ji}^{\rm b}$, a minimum value of ${\cal N}$ is required (Lucy 1999b).

The above estimators refer only to the radiative rates that enter the equations of statistical equilibrium. Further quantities are required by the equation of thermal equilibrium. These include the coefficients  $\gamma_{i}^{\rm E}$ and $\alpha^{\rm E, st}_{i}$ defined in Sect. 3.6. As is evident from their definitions, estimators for these quantities are obtained from those derived above for  $\gamma_{i}$ and  $\alpha^{\rm st}_{i}$ simply by changing $h \nu$ to $h \nu_{i}$ in the factors $a_{i \kappa}(\nu)/h \nu$ - see Eqs. (44) and (45).

Finally, an estimate is also needed for the heating rate due to f-f absorptions,

$${\cal H}^{\rm ff} = 4 \pi \int_{0}^{\infty} k^{\rm ff}_{\nu} \rho \: J_{\nu} \: {\rm d} \nu.$$ (47)

If we separate out the dependence on level population by writing ${\cal H}^{\rm ff} = h_{\kappa}^{\rm ff} n_{\kappa}$ and then again apply Eq. (43), the estimator for the coefficient is

$$h_{\kappa}^{\rm ff} =n_{\rm e} \frac{\epsilon_{0}}{\Delta t} ... ...lta s \: \frac {k^{\rm ff}_{\nu} \rho}{n_{\kappa}n_{\rm e}}$$ (48)

where the summation is over all pathlengths in V and $\nu$ is evaluated at their mid-points.

The summations in the above estimators are accumulated in the course of the MC simulation and are the only quantities concerning the internal radiation field that need to be stored. When the trajectories of all ${\cal N}$ packets have been computed, the scale factor  $\epsilon_{0}/\Delta t$ is determined as described in Sect. 6.1 and then used to convert the summations into the required rate coefficients.

Although the accuracy of these various estimators remains to be demonstrated, it is clear from their derivation that they are consistent. In other words, as ${\cal N} \rightarrow \infty$ and ${\cal M} \rightarrow \infty$, they converge to the quantities being estimated.

6.3 Statistical equilibrium

The required solution is in statistical and thermal equilibrium. Accordingly, following each MC calculation, the equations representing both of these equilibria are used simultaneously to derive improved values of n_i and T throughout the envelope. Nevertheless, it is convenient to discuss these equilibria separately.

The equation of statistical equilibrium for level i can be written in the form

$$\Lambda_{\ell i} n_{\ell} - \left(\Lambda_{i \ell} + \Lambda_{i {\rm u}}\right) n_{i} + \Lambda_{{\rm u} i}n_{\rm u}=0$$ (49)

where $\Lambda_{ij}$ denotes the total rate coefficient for the transition $i \rightarrow j$. Specifically, for transitions between bound levels, we have (cf. Eqs. (9) and (10)).

$$\Lambda_{ji} = B_{ji} \beta_{ji} J_{ji}^{\rm b} + q_{ji} n_{\rm e}$$ (50)

for excitations $j \rightarrow i$, and

$$\Lambda_{ij} = A_{ij}\beta_{ji} + B_{ij} \beta_{ji} J_{ji}^{\rm b} + q_{ij}n_{\rm e}$$ (51)

for de-excitations $i \rightarrow j$. In addition, for transitions coupling bound levels to the continuum, we have

$$\Lambda_{i \kappa} = \gamma_{i} + q_{i \kappa} n_{\rm e}$$ (52)

for ionizations $i \rightarrow \kappa$, and

$$\Lambda_{\kappa i} = \left(\alpha_{i}^{\rm sp}+\alpha_{i}^{\rm st}\right)n_{\rm e} + q_{\kappa i} n_{\rm e}^{2}$$ (53)

for recombinations $\kappa \rightarrow i$.

Equation (49) applies to all levels $i = 1 \ldots \kappa$. But since every source term is a loss term for another level, the sum of the left hand sides is zero, which implies that one of the equations is redundant. This leaves  $\kappa - 1$ equations in the ${\kappa }$ unknown level populations n_i. To make the system determinate, we impose the constraint

$$n_{1}+n_{2}+ \cdots + n_{\kappa} = n_{\rm H}$$ (54)

where $n_{\rm H}$ is the known number density of H atoms - see Sect. 2.1.

The coefficients $\Lambda_{ij}$ depend on $n_{\kappa} (=$ $n_{\rm e})$ and on the n_i of bound levels through the escape probabilities  $\beta_{ji}$. Accordingly, the system of equations defined by Eqs. (49) and (54) is non-linear. But fortunately the system can be solved iteratively with a standard package for linear equations. Thus the coefficients are evaluated with the current estimates of n_i and then the system solved for improved estimates. These are then used to recompute the $\beta_{ji}$ and hence the $\Lambda_{ij}$ and the process repeated until convergence is achieved.

However, at this point a limitation in the MC approach to NLTE problems becomes apparent. The problem arises for the populations n_iof closely-spaced levels near the continuum. Monte Carlo sampling errors in the estimated rates of the processes populating such levels can give rise to level inversions, and these then prevent convergence of the above back substitutions. To overcome this, the following stabilizing device is adopted: for each level j, the n_i of a higher level i is replaced by $n_{\rm c} = n_{j}g_{i}/g_{j}$if $n_{i} > n_{\rm c}$.

Because this problem only arises for high levels in the outermost shells, its impact on the MC code's prediction of the SN's line- and continuum spectrum from the UV to the mid-IR is totally negligible. Moreover, as expected, the occurrence of these spurious inversions decreases as the size of the MC simulation increases. For example, with ${\cal N} = 2 \times 10^{4}$ in an envelope stratified into ${\cal M} = 20$shells, inversions occurred in 13 shells starting at m = 2 and affecting levels as low as i=6. With ${\cal N}$ increased to $5 \times 10^{4}$, inversions occur in only 5 shells starting at m = 8 and the lowest affected level is i = 8. With a modest further increase to ${\cal N} \ga 10^{5}$, inversions seldom occur. Nevertheless, problems requiring accurate population ratios for closely-spaced levels are perhaps best not tackled with MC methods.

6.4 Thermal equilibrium

In addition to statistical equilibrium, the required NLTE solution must also be in thermal equilibrium. This latter equilibrium is simply the condition that throughout the envelope

$${\cal G} = {\cal H} - {\cal C} = 0$$ (55)

where ${\cal H}$ and ${\cal C}$ are the total heating and cooling rates, respectively. Given the physical processes included in this calculation, we have

$${\cal H} = {\cal H}^{\rm bf} + {\cal H}^{\rm ff} + {\cal H}^{{\rm c} \ell}$$ (56)

and

$${\cal C} = {\cal C}^{\rm fb} + {\cal C}^{\rm ff} + {\cal C}^{{\rm c} \ell}.$$ (57)

In terms of the modified photoionization coefficient defined in Sect. 3.6, the rate at which b-f transitions convert radiant energy into thermal energy is

$${\cal H}^{\rm bf} = \sum_{1}^{\kappa -1} \: h_{i}^{\rm bf} n_... ...} = \left(\gamma_{i}^{\rm E} - \gamma_{i}\right) \: h \nu_{i}.$$ (58)

The corresponding rate at which f-b transitions convert thermal energy into radiant energy is

$${\cal C}^{\rm fb} = c_{\kappa}^{\rm fb} n_{\kappa} \;\;\;\; \... ...: \left(\alpha_{i}^{\rm E} - \alpha_{i}\right) \: h \nu_{i}.$$ (59)

The heating and cooling rates ${\cal H}^{\rm ff}$ and ${\cal C}^{\rm ff}$due to f-f processes have previously been defined in Eqs. (47) and (32). The cooling rate ${\cal C}^{{\rm c} \ell}$ due to collisional excitations is given by Eq. (33), and the corresponding heating rate ${\cal H}^{{\rm c} \ell}$ due to de-excitations is computed similarly.

When the MC model of the internal radiation field has been computed, the current values of a shell's level populations and temperature, $\bar{n}_{i}$ and $\bar{T}$, do not in general correspond to thermal equilibrium - i.e., $\bar{{\cal G}} \neq 0$ . To eliminate this residual, we apply corrections  $\delta n_{i}$ and $\delta T$ that satisfy the linearized version of Eq. (55),

$$\bar{{\cal G}} + \frac{\partial{\cal G}}{\partial n_{i}} \delta n_{i} +\frac{\partial{\cal G}}{\partial T} \delta T = 0.$$ (60)

But a more convenient form of this equation follows from noting that each level's contribution to the net heating rate can be expressed in terms of heating (h_i) and cooling (c_i) coefficients - see Eqs. (48), (58) and (59) - so that $G = \sum (h_{i}-c_{i})n_{i}$. With this substitution, Eq. (60) becomes

$$\sum_{1}^{\kappa} \left(h_{i}-c_{i}\right)n_{i} +\frac{\partial{\cal G}}{\partial T} \delta T = 0$$ (61)

where $\partial {\cal G}/ \partial T$ is computed by numerical differencing, and $n_{i} = \bar {n}_{i} + \delta n_{i}$ is the improved level population.

Equation (61) is one equation in $\kappa +1$ unknowns. But if we add it to Eqs. (49) and (54), we have $\kappa +1$ equations in the $\kappa +1$ unknowns n_i and $\delta T$. This combined system of equations, representing statistical and thermal equilibrium, is solved by repeated back substitution as described in Sect. 6.3. With each back substitution, the coefficients  $\Lambda_{ij}$, h_i, c_i and $\partial {\cal G}/ \partial T$ are recomputed because of their dependences on n_i, $n_{\rm e}$ and T.

This simple iteration technique for finding n_i and $\delta T$ must often start far from the final solution and can then fail due to large initial corrections. To avoid such failures, an undercorrection factor u = 0.8 is applied to the corrections  $\delta n_{i}$ and $\delta T$, and the temperature change is further limited to $\pm$$\Delta T$, with  $\Delta T = 200$ K. Moreover, if the magnitude of the correction vector  $\delta n_{i}$ increases from one iteration to the next, the factor u and the bound $\Delta T$ are both decreased. With these precautions, back-substitution proves to be a robust technique for solving the equations of statistical and thermal equilibrium.

The iteration technique also fails when sampling errors' contribution to ${\cal H} - {\cal C}$ itself implies a large correction $\delta T$. The only remedy in this case is to increase ${\cal N}$.

The iterations are deemed to have converged when $\vert\delta T\vert < 1$ K and $\overline{\vert\delta n_{i}\vert/n_{i}} < 10^{-5}$, where the averaging is over levels i=1-5 because of the occasional spurious inversions of high levels (Sect. 6.3).

When the iterations have converged to this accuracy, the typical residuals for each shell are $\vert{\cal H}-{\cal C}\vert / {\cal C} < 10^{-4}$ and similarly for the agreement between the source and sink terms, $\Lambda_{\ell i} n_{\ell} + \Lambda_{{\rm u} i}n_{\rm u}$ and $(\Lambda_{i \ell} + \Lambda_{i {\rm u}}) n_{i}$, for all levels i in Eq. (49). However, when the n_i of a high level has been stabilized to avoid a spurious inversion, the fractional departure from statistical equilibrium rises to $\sim$10^-3.

6.5 Outer iteration

During the above (inner) iterations to find n_i and T, the pathlength summations of Sect. 6.2 are of necessity held fixed. In effect, therefore, the matter in each shell is brought into equilibrium with a fixed radiation field. Accordingly, if the MC radiation field were now to be recomputed with the improved n_i and T, these summations would change, and so statistical and thermal equilibrium would again be violated. Evidently, an outer iterative loop is necessary if convergence to the NLTE solution is to be achieved.

The scheme adopted for these outer iterations is simply to repeatedly input the updated n_i and T from Sect. 6.4 into the MC calculation of Sect. 5. These outer iterations are continued until the changes in n_i and T from one outer iteration to the next can be attributed to MC sampling errors.

The novelty of this iterative technique is that, despite well-documented convergence failures in various earlier applications (see, e.g., Mihalas 1978), we are relying on the $\Lambda$-iteration device of repeatedly bringing matter into statistical and thermal equilibrium with the just- updated MC radiation field. Indeed, these $\Lambda$-iterations would also fail to converge if the radiation field were obtained by solving the Equation of Radiative Transfer (RTE). But in this MC scheme, the radiation field is subject at every iteration to constraints that other iterative schemes only aim to achieve asymptotically. Thus, even while n_i and T still depart from their equilibrium values, the constraint of radiative equilibrium in the matter frame is obeyed rigorously. Moreover, the MC transition probabilities approximately incorporate the constraints that statistical equilibrium imposes on the frequency redistribution of absorbed radiant energy. In effect, therefore, the scheme introduces a constrained $\Lambda$-operator that, for sufficiently large ${\cal N}$, generates at every iteration a radiation field that is closer to the converged solution than would be the radiation field generated by the RTE. The convergence of these constrained $\Lambda$-iterations is illustrated in Sect. 7.2.

7 Numerical results

In Sects. 3-6, a detailed description has been given of how the MC model of the internal radiation field is derived and how it is used to improve level populations and temperatures throughout the envelope. But whether such improvements will actually result in convergence to the NLTE solution remains to be demonstrated.

7.1 Accuracy of estimators

Crucial to the success or otherwise of this technique is the accuracy of the estimators developed in Sect. 6.2. Even with consistent estimators, the technique would have little current interest if adequate accuracy requires simulations with impossibly large ${\cal N}$. On the other hand, if accurate rates are achievable with feasible values of ${\cal N}$, then the coefficients in, and therefore the solutions of, the equations of statistical and thermal equilibrium will differ little from those derivable with a standard transfer calculation.

An obvious accuracy test would be to obtain an accurate solution of a NLTE problem using a conventional code and then examine the convergence of MC solutions of the same problem as ${\cal M}$ and ${\cal N}$ increase. But a simpler procedure is followed here: the special case of free-streaming r-packets allows values of the various coefficients obtained with the MC estimators to be compared with values obtained by numerical integration.

The MC code is easily modified so that r-packets propagate without interaction. Each packet is then emitted as before from the moving lower boundary but thereafter conserves its  $\epsilon_{R}$ and $\nu_{R}$, while its  $\epsilon_{\nu}$ and $\nu$ decrease along the trajectory because of the differential expansion.

A further modification made for this free-streaming test is to omit stratified sampling of the r-packets' frequencies at the lower boundary - see Sect. 4.1. Thus, the ${\cal N}$ bins in $\ell n \: \nu$ are here randomly not uniformly sampled. By thus not benefiting from stratified sampling, the results should be more typical of circumstances where the radiation field is created within D rather than emitted by an artificial boundary condition.

At radius r, the conical radiation field modelled by this free-streaming MC calculation has specific intensity $I_{\nu}(\mu) = B_{\nu'}(T_{\rm b})$, where $\nu'(\mu)$ is the frequency at r = R of a photon that reaches r with frequency $\nu$ and direction cosine $\mu$. From this expression for $I_{\nu}$, the mean intensity  $J_{\nu}(r)$ can be computed to arbitrary accuracy by numerical integration. Such calculations carried out at fequencies $\nu_{ji}$ and at the mean radii  $\bar{r}_{\rm m}$ then allow the estimator $J_{ji}^{\rm b}$ given by Eq. (46) to be tested. Similarly, $J_{\nu}(\bar{r}_{\rm m})$ can be computed for a grid of frequencies and then used to obtain accurate values for the coefficients  $\gamma_{i}$, $\gamma_{i}^{\rm E}$, $\alpha_{i}^{\rm st}$, $\alpha_{i}^{\rm E, st}$ and  $h_{\kappa}^{\rm ff}$ by numerical integration, thereby testing their MC estimators from Sect. 6.2.

These comparisons, carried out for ${\cal M} = 10{-}100$ and  ${\cal N} = 10^{4}{-}10^{6}$, support the assertion of Sect. 6.2 that these estimators are consistent. Thus numerical experiments indicate that the errors of the MC estimates  $\rightarrow$0 as ${\cal M}$ and ${\cal N} \rightarrow \infty$, with no suggestion of bias. But if ${\cal M}$ is fixed while ${\cal N} \rightarrow \infty$, the estimators' sampling errors  $\rightarrow$0 but biases due to discretization errors remain. On the other hand, if ${\cal N}$ is fixed while ${\cal M} \rightarrow \infty$, biases  $\rightarrow$0 but sampling errors become large, especially for the mean intensities  $J_{ji}^{\rm b}$.

In addition to confirming the estimators' consistency, this special case also suggests that acceptable accuracies are indeed achieved with feasible values of ${\cal N}$. For example, with ${\cal M} = 20$ and  ${\cal N} = 10^{5}$, the MC free-streaming simulation requires only about 3 min computer time and yet gives reasonably accurate coefficients. Thus, in one such simulation, the means of the fractional and of the absolute fractional errors of the coefficients  $\gamma_{i}$ for $i = 2, \ldots ,\kappa-1$ and  $m = 1, \ldots ,{\cal M}$ are  $2.8 \times 10^{-3}$ and $6.1 \times 10^{-3}$, respectively. These gratifyingly small errors are a consequence of the large number of r-packets contributing to the summation in Eq. (44) and confirm the efficiency of estimators constructed from Eq. (43).

The quantities least accurately estimated are the $J_{ji}^{\rm b}$. Thus, in the above simulation, the mean fractional and absolute fractional errors of the $J_{ji}^{\rm b}$ for the Balmer lines in all shells are $0.4 \times 10^{-2}$ and $12.6 \times 10^{-2}$, respectively. In this case, sampling errors clearly dominate over any systematic bias due to discretization errors. This is no suprise since the summation in Eq. (46) includes only those r-packets that are red-shifted into resonance with the transition $j \rightarrow i$. Accordingly, even with ${\cal N} = 10^{5}$, the number of packets contributing to the summation for a particular  $J_{ji}^{\rm b}$ is not large. For example, in the case of the Balmer lines, the number is $\sim$35 at m = 1 rising to $\sim$300 at  $m = {\cal M}$.

These relatively large sampling errors for the $J_{ji}^{\rm b}$are the main source of error for the coefficients in the equations of statistical equilibrium. Of course, for a typical level, there are several source terms, some not even subject to sampling errors; and so the fractional sampling error of the total input rate will be correpondingly reduced. Nevertheless, for a level populated predominately by a single radiative b-b transition, its n_iwill directly reflect the sampling error of the relevant  $J_{ji}^{\rm b}$. But since sampling errors differ from one iteration to the next, levels suffering this sensitivity will display fluctuating values of n_i and so should not escape notice.

7.2 Convergence of temperature profile

To investigate convergence to the NLTE solution, the values of T and n_i throughout the envelope were monitored for 20 outer iterations in a calculation with ${\cal M} = 20$ shells and ${\cal N} = 4 \times 10^{5}$ packets. The envelope's parameters are as specified in Sect. 2.2, and the calculation is initiated as described in Sect. 3.4. Thus, for the first iteration, the MC transition probabilities are computed with radiative and collisional rates corresponding to the initial analytical model of the radiation field - Eq. (4) - and to the inital guesses for n_i and T. For all subsequent iterations, the MC transition probabilities are computed from the MC radiation field via the MC estimators of Sect. 6.2, and the n_i and T are the updated values resulting from bringing the matter into statistical and thermal equilibrium with the MC radiation field.

Figure 1 shows how the envelope's temperature stratification evolves from the initial T = 5000 K as the iterations proceed. The corrections are seen to be substantial for iterations 1-4 and negligibly small thereafter. Thus the scheme appears to converge at the fifth iteration, with the changes introduced at iterations 6-20 presumably being random fluctuations about the exact NLTE solution. Such fluctuations are expected because of different sampling errors in successive realisations of the MC radiation field.

	Figure 1: Convergence of constrained $\Lambda$-iterations. Results of 20 temperature iterations starting from an isothermal stratification with T = 5000 K are plotted for a simulation with ${\cal N} = 4 \times 10^{5}$ packets. The solution for the conical radiation field is plotted as open circles.
Open with DEXTER

To support this claim that the scheme has indeed achieved rapid convergence to the close neighbourhood of the exact NLTE solution, a separate NLTE calculation has been carried out with line formation again treated in the Sobolev approximation but now with the continuum radiation field approximated by the free-streaming model described in Sect. 7.1. With the radiation field thus prescribed, the determination of n_i and T at radius r from the equations of statistical and thermal equilibrium is a one-point problem. The coefficients in these equations are calculated with accurate numerical integrations as for the accuracy tests of Sect. 7.1 and then n_i and T obtained with the back-substitution iterations of Sect. 6.4. This one-point calculation has been carried out at several radii, with the results plotted as open circles in Fig. 1. The agreement with the MC solutions for iterations 6-20 is excellent.

Given that the scheme has indeed converged and that the remaining temperature fluctations are therefore due to sampling errors, the next question is whether such fluctuations are unacceptably large. Assuming convergence at iteration 6, we derive an accurate estimate  $\overline{\rm T}$ of the exact Tfor each shell by averaging over iterations 6-20. The standard deviation of the fluctuations $T-\overline{\rm T}$ can then be computed for each shell as well as a consolidated value  $\sigma_{\rm T}$ for the entire envelope. For this simulation with ${\cal N} = 4 \times 10^{5}$, we find that $\sigma_{\rm T} = 9.0$ K. Temperature uncertainties of this magnitude are admittedly large compared to the nominal precision - often <0.01 K - achieved by conventional NLTE codes. But in view of the parametric, geometrical and kinematic uncertainties associated with specifying the structure of any celestial object, errors of $\sim$10 K are completely inconsequential.

The amplitude of the fluctuations following convergence depend of course on ${\cal N}$. Repeating the above simulation with  ${\cal N} = 10^{5}$ and again averaging over iterations 6-20, we find that  $\sigma_{\rm T}$ increases to 16.2 K, consistent with the expected scaling $\sigma_{\rm T} \propto {\cal N}^{-1/2}$.

7.3 Convergence of level populations

	Figure 2: Convergence of constrained $\Lambda$-iterations. Results of 20 level population iterations starting from Saha-Boltzmann populations at T = 5000 K are plotted for the mass shell m = 10 in a simulation with ${\cal N} = 4 \times 10^{5}$packets. The solution for the conical radiation field is plotted as open circles.
Open with DEXTER

The convergence behaviour of the n_i for the same 20 iterations of the above model with ${\cal M} = 20$ shells and ${\cal N} = 4~\times$10⁵ packets is shown in Fig. 2 for shell m = 10. Again we see substantial corrections in the first few iterations followed by fluctations of moderate amplitude about the scheme's solution in the limit  ${\cal N} \rightarrow \infty$.

As for the temperature iterations, the MC solution is compared to the accurate NLTE solution for the conical radiation field of Sect. 7.1. The agreement is excellent for low levels and for the continuum (${\kappa }$) but small discrepancies are evident for high levels. This probably reflects a small difference in the two calculations due to the device used in the MC calculation to avoid negative probabilities when $b_{\kappa} > b_{i}$- see Sect. 4.3.1.

A further test of the accuracy of the converged n_i is presented in Fig. 3. In this diagram, the variation with radius of the departure coefficients for several levels are compared to the corresponding values for the conical radiation field. The agreement is seen to be satisfactory.

The behaviour of the H level populations in this SN envelope are qualitatively similar to those of the hydrogenic ion He⁺ in W-R winds (Hillier 1987). In particular, photoionizations predominantly occur from level 2, which behaves like a ground state, becoming overpopulated in the outer layers due to the dilution of the ionizing continuum shortward of $\lambda$3646 Å and to the trapping of Lyman $\alpha$ photons.

Note that, because most H atoms are in the ground state for both the NLTE and the LTE solutions, the departure coefficient b₁ is almost exactly = 1. Accordingly, the values b₂ plotted in Fig. 3 directly indicate the decoupling of level n = 2 from the ground state.

7.4 Consistency

Given that the macro-atom formalism has now been extended to include interactions with the thermal pool, it is of interest to test this extension by repeating the consistency test of level emissivities described for H in Sect. 5.2 of Paper I.

When the outer iterations have converged, the level populations, temperature and radiative rates appropriate to statistical and thermal equilibrium are known for each mass shell. From this information, the rates at which radiant energy is being converted to other forms can be computed for each of the processes operating. Specifically, these are: conversions into excitation energy by b-b absorptions to levels $2, \ldots , \kappa-1$; conversion into ionization energy by b-f absorptions; and conversions into thermal energy by b-f and f-f absorptions. These last two processes, labelled $\kappa +1$ and $\kappa +2$, were not included in Paper I, which treated only statistical equilibrium.

	Figure 3: Test of departure coefficients bi. Results obtained at iteration 20 in Monte Carlo calculation with ${\cal N} = 4 \times 10^{5}$ packets are plotted as solid lines for the indicated levels. The solutions for the conical radiation field are plotted as filled circles (i=2), open circles (i=3) and asterisks (i=10).
Open with DEXTER

From these conversion rates, a Monte Carlo test of the MC transition probabilities is carried out for each mass shell as follows: N = 10⁶equal energy r-packets are assigned to these absorption channels in proportion to the computed conversion rates. Those assigned to processes $2, \ldots ,i, \ldots ,\kappa$ activate a macro-atom to state i, while those assigned to processes $\kappa +1$ and $\kappa +2$ result in the creation of a k-packet. Application of the MC transition probabilities then results eventually in the emission of an r-packet as described in Sect. 4. Moreover, each emitted r-packet belongs to one of the $\kappa +1$emission processes representing the inverse of the above absorption processes. Accordingly, when all N packets have been absorbed and re-emitted, we have experimental values N_i for each of the $\kappa +1$emission processes.

But the rates at which excitation, ionization and thermal energy is being converted in a given shell into radiant energy can be directly computed from the known level populations and temperature - i.e., from the rates of the various b-b, f-b and f-f processes. Thus, we can predict that N^*_i of the N packets should have been emitted by process i.

In Paper I, Fig. 4, the agreement of N_i and N^*_i was shown graphically. Here, more rigorously, we compute the values of $\chi^{2} = \sum_{i} (N_{i}-N^{*}_{i})^{2}/N^{*}_{i}$ for each mass shell. Since there are 16 emission channels, we expect the experimental values of $\chi^{2}$ to be distributed as $\chi^{2}_{\nu}$ with $\nu = 15$ degrees of freedom. Thus, 50 percent of the values should be <14.34 and 95 percent <25.00

The results obtained after the 20th iteration for the model with ${\cal M} = 20$ discussed in Sect. 7.2 are as follows: the 20 values of $\chi^{2}$ range from 7.87 to 34.22, with 6 < 14.34 and 3 > 25.00.

The outcome of this test is that the experimental values of $\chi^{2}$ are slightly biased to higher values than predicted by the $\chi^{2}_{\nu}$ distribution. However, given that the underlying model has departures from strict statistical and thermal equilibrium because of its own MC sampling errors, these results are in fact a strong confirmation that the extension of the macro-atom formalism to include interactions with the thermal pool has been carried out correctly.

As a further test of this conclusion, values of $\chi^{2}$ were computed with the summation restricted to the three channels containing r-packets emitted by f-b and f-f processes. In this case, the values for all shells are <7.81, the 95 percent confidence limit for $\nu = 3$degrees of freedom. This eliminates the possibility that the above upward bias is due to poor predictions for the continuum channels.

7.5 Monte Carlo spectrum

When the outer iteration loop has converged, the r-packets that escaped to $\infty$ during the last iteration provide a crude estimate of the SN's emergent spectrum. This is obtained by simply allocating each escaping packet's rest energy  $\epsilon_{R}$ to its appropriate bin in a grid of rest frequencies. This crude spectrum is shown in Fig. 4 for the above model with ${\cal N} = 4 \times 10^{5}$ packets. In this figure, luminosity density is plotted against vacuum wavelength in the optical domain and, despite sampling errors, the P Cygni line profiles of H$\alpha$ and H$\beta$ are clearly seen.

If this MC code were intended for use in interpreting observational data, the crude MC spectrum shown in Fig. 4 would not be compared with the observed spectrum. Instead, a far less noisy theoretical spectrum would be computed by following a procedure described earlier (Lucy 1999b). Specifically, estimates of line- and continuum source functions can be derived from the same MC simulation and used to compute a high quality spectrum from the formal integral for the emergent intensity as a function of impact parameter.

	Figure 4: Monte Carlo spectrum obtained at iteration 20 in Monte Carlo calculation with ${\cal N} = 4 \times 10^{5}$ packets. The unit of luminosity density is 1045 erg s-1 cm-1. The rest wavelengths of H$\alpha$, H$\beta$ and H$\gamma$ are indicated.
Open with DEXTER

8 Discussion and conclusion

In Paper I, the radiative and collisional interactions of an atomic species in statistical equilibrium with its environment was modelled in terms of macro-atoms being activated and de-activated by the absorption and emission of e-packets. In this paper, this macroscopic quantization has been extended to include interactions with the thermal pool on the assumption of local thermal balance - i.e., ${\cal H} = {\cal C}$. (But note that departures from thermal equilibrium due to non-radiative heating can be incorporated as described in Sect. 7.2 of Paper I.)

With this extension, the macro-atom formalism allows MC transfer codes to obtain NLTE solutions for stratified atmospheres. As a first test of this, this paper describes a MC code treating the formation of the H spectrum in a type II SN, with the emphasis being on demonstrating the consistency of the technique as ${\cal N} \rightarrow \infty$, the accuracy of the estimators of the various radiative rates for feasible values of ${\cal N}$, and the convergence behaviour of the constrained $\Lambda$-iteration scheme.

Given that consistency has been demonstrated and that the geometry-independent $\Lambda$-iteration scheme converges, applications to problems with arbitrary geometry and to multiple-species plasmas are now in principle possible. But considering the detail and resolution required to have impact on a current research topic, such problems will undoubtedly require a substantial increase in computational resources over those deployed here for a 1D envelope of pure H. Even with the efficient estimators of Sect. 6.2, each cell in the numerical grid must still be traversed by numerous r-packets if adequate accuracy is to be achieved. This requirement obviously mandates a huge increase in ${\cal N}$ for 2- and 3D problems. But also in 1D problems for static atmospheres, narrow line widths will demand large values of ${\cal N}$ in order to get accurate b-b transition rates.

A further problem arises if the lower boundary of the computational domain is at a large optical depth in the continuum. Each r-packet's interaction history will then comprise numerous events and will often end with a recrossing of the lower boundary, thus reducing the number of r-packets traversing the cells in the outer atmosphere.

This requirement for substantial computer resources suggests the use of a computer with numerous parallel processors. Fortunately, this MC technique is extremely well suited for such machines, especially those with large shared memory. A MC transfer calculation can then be equally divided between the processors, each of which carries its task to completion with no exchanges of information or packets with other processors. Following completion, the individual processors' contributions to the estimator summations are added to obtain the radiative rates for all cells, which are then equally divided among the processors to solve the equations of statistical and thermal equilibrium. These remarks suggest that close to the maximum theoretical efficiency should be achievable.

The long-term aim in devoloping the macro-atom formalism is to obtain accurate NLTE solutions for multi-dimensional problems using MC methods. But the technique is likely to be of more immediate use in obtaining approximate NLTE solutions with precision sufficient for current research problems. The reasons for optimism in this regard are that, even without the converged stratifications of temperature and level populations, the technique rigorously obeys radiative equilibrium and also provides remarkably accurate emissivites (Paper I, Sect. 6). Thus, for some problems, adequate accuracy will be achieved by using analytic ionization and excitation formulae, as in previous stellar wind and SNe codes, but then computing the radiation field with this MC technique replacing the earlier techniques that assumed coherent scattering (Abbott & Lucy 1985) or downward branching (Lucy 1999b). For yet higher accuracy, ionization could be solved for while retaining an analytic excitation formula. Compared to the full NLTE problem, this greatly reduces the demand for large ${\cal N}$ since estimates of b-b rates are no longer required. Moreover, for each atomic species, the statistical equilibrium equations for perhaps thousands of levels are replaced by the ionization equations for just $\sim$3-5 ions. Despite this enormous simplification, the anticipated loss of accuracy is slight in view of the afore-mentioned accurate MC emissivities.

Acknowledgements

I am grateful to S.Sim for a discussion on the potential application of parallel processors.

References

 
• Abbott, D. C., & Lucy, L. B. 1985, ApJ, 288, 679 In the text NASA ADS  
• Arnett, W. D. 1988, ApJ, 331, 377 In the text NASA ADS  
• Castor, J. I. 1970, MNRAS, 149, 111 In the text NASA ADS  
• Duschinger, M., Puls, J., Branch, D., Hoeflich, P., & Gabler, A. 1995, A&A, 297, 802 In the text NASA ADS  
• Hillier, D. J. 1987, ApJS, 63, 947 In the text NASA ADS  
• Lucy, L. B. 1971, ApJ, 163, 95 In the text NASA ADS  
• Lucy, L. B. 1987, in ESO Workshop on SN 1987A, ed. I. J. Danziger, 417 In the text  
• Lucy, L. B. 1999a, A&A, 344, 282 In the text NASA ADS  
• Lucy, L. B. 1999b, A&A, 345, 211 In the text NASA ADS  
• Lucy, L. B. 2002, A&A, 384, 725 In the text NASA ADS  
• Mazzali, P. A., & Lucy, L. B. 1993, A&A, 279, 447 In the text NASA ADS  
• McCrea, W. H., & Mitra, K. K. 1936, Zs.f.Ap., 11, 359 In the text  
• Mihalas, D. 1978, Stellar Atmospheres 2nd ed. (San Francisco: W. H. Freeman & Co.) In the text  
• Osterbrock; D. E. 1974, Astrophysics of Gaseous Nebulae (San Francisco: W. H. Freeman & Co.) In the text  
• Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes (Cambridge: Cambridge Univ. Press) In the text  
• Suntzeff, N. B., & Bouchet, P. 1990, AJ, 99, 650 In the text NASA ADS  
• van Regemorter, H. 1962, ApJ, 136, 906 In the text NASA ADS