Derivation of Methods

For a two-level system in slab geometry that includes collisional pumping between atomic levels, the radiation transport equation is

$${\frac{{{\partial f}}}{{{\partial t}}}} \mu {\frac{{{\partial f}}}{{{\partial x}}}} {\frac{{{n_2 }}}{{2}}} \phi \left(\vphantom{ {K_{12} n_1 - K_{21} n_2 } }\right. \left.\vphantom{ {K_{12} n_1 - K_{21} n_2 } }\right) \phi$$ (1)

where c is the speed of light, x is the position in the slab, $\mu$ is the direction cosine of the radiation, $\nu$ is the frequency of the radiation, f ($\mu$,$\nu$, x, t) is the photon number density distribution per unit atom density, n₂(x, t) is the upper level population fraction, n₁(x, t) is the lower level population fraction, A₂₁ is the spontaneous emission rate, $\phi$(a,$\nu$) is the Voigt line profile normalized to unit integral [7], and K₁₂ = $\kappa$N where $\kappa$ is the lower state absorption cross section and N is the atom number density. The coefficient K₂₁ is defined by

$${\frac{{{g_1 }}}{{{g_2 }}}}$$ (2)

where g₁ and g₂ are the statistical weights for levels 1 and 2, respectively. As in [3], we consider the problem in the regime of complete redistribution and no physical scattering of photons.

The equations governing the atomic population fractions n₁ and n₂ are

$${\frac{{{dn_2 }}}{{{dt}}}} \left(\vphantom{ {K_{12} n_1 - K_{21} n_2 } }\right. \left.\vphantom{ {K_{12} n_1 - K_{21} n_2 } }\right)$$ =

$$\int_{{ - 1}}^{1} \mu \int_{0}^{\infty} \nu \phi \nu \mu \nu$$ (3)

and

n₁ + n₂ = 1, (4)

where C₁₂ and C₂₁ are rate constants for the collisional transitions 1 $\rightarrow$ 2 and 2 $\rightarrow$ 1, respectively. One must also add in appropriate boundary conditions and initial state to the above equations.

Using (4), equations (1) and (3) may be rewritten as

$${\frac{{{\partial f}}}{{{\partial t}}}} \mu {\frac{{{\partial f}}}{{{\partial x}}}} {\frac{{n}}{{2}}} \phi \left[\vphantom{ {K_{12} - \left( {K_{21} + K_{12} } \right)n} }\right. \left(\vphantom{ {K_{21} + K_{12} } }\right. \left.\vphantom{ {K_{21} + K_{12} } }\right) \left.\vphantom{ {K_{12} - \left( {K_{21} + K_{12} } \right)n} }\right] \phi$$ (5)

and

$${\frac{{{dn}}}{{{dt}}}} \left(\vphantom{ {C_{12} + C_{21} + A_{21} } }\right. \left.\vphantom{ {C_{12} + C_{21} + A_{21} } }\right) \left[\vphantom{ {K_{12} - \left( {K_{21} + K_{12} } \right)n} }\right. \left(\vphantom{ {K_{21} + K_{12} } }\right. \left.\vphantom{ {K_{21} + K_{12} } }\right) \left.\vphantom{ {K_{12} - \left( {K_{21} + K_{12} } \right)n} }\right]$$ =

$$\int_{{ - 1}}^{1} \mu \int_{0}^{\infty} \nu \phi \left(\vphantom{ \nu }\right. \nu \left.\vphantom{ \nu }\right) \left(\vphantom{ {\mu ,\nu } }\right. \mu \nu \left.\vphantom{ {\mu ,\nu } }\right)$$ (6)

respectively, where n is the upper level population fraction.

We can generate a finite differencing scheme in time for (6) by using the standard IMC technique [3]. We integrate (6) from t₀ to t₀ + $\Delta$t. In the spontaneous emission and collision terms, we approximate n(t) by n(t₀ + $\Delta$t). In the absorption term, we substitute n(t₀) for n(t) and obtain

$$\Delta \left[\vphantom{ {C_{12} - \left( {C_{12} + C_{21} + A_{21} } \right)n \left( {t_0 + \Delta t} \right)} }\right. \left(\vphantom{ {C_{12} + C_{21} + A_{21} } }\right. \left.\vphantom{ {C_{12} + C_{21} + A_{21} } }\right) \left(\vphantom{ {t_0 + \Delta t} }\right. \Delta \left.\vphantom{ {t_0 + \Delta t} }\right) \left.\vphantom{ {C_{12} - \left( {C_{12} + C_{21} + A_{21} } \right)n \left( {t_0 + \Delta t} \right)} }\right] \Delta$$ =

$$\left[\vphantom{ {K_{12} - \left( {K_{21} + K_{12} } \right)n\left( {t_0 } \right)} }\right. \left(\vphantom{ {K_{21} + K_{12} } }\right. \left.\vphantom{ {K_{21} + K_{12} } }\right) \left(\vphantom{ {t_0 } }\right. \left.\vphantom{ {t_0 } }\right) \left.\vphantom{ {K_{12} - \left( {K_{21} + K_{12} } \right)n\left( {t_0 } \right)} }\right]$$

$$\int_{{t_0 }}^{{t_0 + \Delta t}} \int_{{ - 1}}^{1} \mu \int_{0}^{\infty} \nu \phi \left(\vphantom{ \nu }\right. \nu \left.\vphantom{ \nu }\right) \left(\vphantom{ {\mu ,\nu ,t} }\right. \mu \nu \left.\vphantom{ {\mu ,\nu ,t} }\right)$$ (7)

In the standard IMC technique [3], we substitute n(t₀ + $\Delta$t) from (7) into the spontaneous emission term of (5) while using n(t₀) in the absorption term. After some algebra, we obtain

$${\frac{{{\partial f}}}{{{\partial t}}}} \mu {\frac{{{\partial f}}}{{{\partial x}}}} {\frac{{{\gamma A_{21} \phi }}}{{2}}} \left[\vphantom{ {n\left( {t_0 } \right) + \Delta tC_{12} } }\right. \left(\vphantom{ {t_0 } }\right. \left.\vphantom{ {t_0 } }\right) \Delta \left.\vphantom{ {n\left( {t_0 } \right) + \Delta tC_{12} } }\right] {\frac{{{\gamma cA_{21} \phi \Delta t}}}{{2}}}$$ =

$$\left[\vphantom{ {K_{12} - \left( {K_{21} + K_{12} } \right)n\left( {t_0 } \right)} }\right. \left(\vphantom{ {K_{21} + K_{12} } }\right. \left.\vphantom{ {K_{21} + K_{12} } }\right) \left(\vphantom{ {t_0 } }\right. \left.\vphantom{ {t_0 } }\right) \left.\vphantom{ {K_{12} - \left( {K_{21} + K_{12} } \right)n\left( {t_0 } \right)} }\right] \int_{{ - 1}}^{1} \mu \int_{0}^{\infty} \nu \phi \left(\vphantom{ \nu }\right. \nu \left.\vphantom{ \nu }\right) \left(\vphantom{ {\mu ,\nu } }\right. \mu \nu \left.\vphantom{ {\mu ,\nu } }\right)$$

$$\phi \left[\vphantom{ {K_{12} - \left( {K_{21} + K_{12} } \right)n\left( {t_0 } \right)} }\right. \left(\vphantom{ {K_{21} + K_{12} } }\right. \left.\vphantom{ {K_{21} + K_{12} } }\right) \left(\vphantom{ {t_0 } }\right. \left.\vphantom{ {t_0 } }\right) \left.\vphantom{ {K_{12} - \left( {K_{21} + K_{12} } \right)n\left( {t_0 } \right)} }\right]$$ (8)

where $\gamma$ is defined as

$$\gamma {\frac{{1}}{{{1 + \Delta t\left( {C_{12} + C_{21} + A_{21} } \right)}}}}$$ (9)

Equation (8) can be interpreted as a transport equation with a net absorption, $\sigma_{a}^{}$, and effective scattering, $\sigma_{s}^{}$, contribution as given by

$$\sigma_{a}^{} \nu \hat{f} \phi \nu \left[\vphantom{ {K_{12} - \left( {K_{21} + K_{12} } \right)n(t_0 )} }\right. \left(\vphantom{ {K_{21} + K_{12} } }\right. \left.\vphantom{ {K_{21} + K_{12} } }\right) \left.\vphantom{ {K_{12} - \left( {K_{21} + K_{12} } \right)n(t_0 )} }\right]$$ (10)

and

$$\sigma_{s}^{} \nu \left(\vphantom{ {1 - \hat f} }\right. \hat{f} \left.\vphantom{ {1 - \hat f} }\right) \phi \nu \left[\vphantom{ {K_{12} - \left( {K_{21} + K_{12} } \right)n(t_0 )} }\right. \left(\vphantom{ {K_{21} + K_{12} } }\right. \left.\vphantom{ {K_{21} + K_{12} } }\right) \left.\vphantom{ {K_{12} - \left( {K_{21} + K_{12} } \right)n(t_0 )} }\right]$$ (11)

where the fraction, $\hat{f}$, is given by

$$\hat{f} {\frac{{{1 + \Delta t\left( {C_{12} + C_{21} } \right)}}}{{{1 + \Delta t\left( {C_{12} + C_{21} + A_{21} } \right)}}}}$$ (12)

The fraction, $\hat{f}$, determines how much effective scattering the problem contains. As $\hat{f}$ approaches unity, effective scattering vanishes. The disadvantage of this method is that the effective scattering term dominates the execution time in optically thick problems, resulting in very long problem runs. Reducing the effective scattering by making $\hat{f}$ approach 1 requires a smaller time step which also increases execution time and may reduce accurracy (see Sec. 2.5).

Symbolic Implicit Monte Carlo [4] achieves implicit time integration using a different point of view. Photons produced by spontaneous emission are given a symbolic weight that remains undetermined until the end of the time step. The Monte Carlo procedure is the same as for IMC, except there is no effective scattering term. Particle scoring results in a linear system of equations that are solved for the upper atomic population fraction at the end of a time step. After spatial discretization (see Sec. 2.3), the upper atomic population is updated with

$$\left(\vphantom{ {t_0 + \Delta t} }\right. \Delta \left.\vphantom{ {t_0 + \Delta t} }\right)_{i}^{} \left[\vphantom{ {C_{12} - \left( {C_{12} + C_{21} + A_{21} } \right)n\left( {t_0 + \Delta t} \right)_i } }\right. \left(\vphantom{ {C_{12} + C_{21} + A_{21} } }\right. \left.\vphantom{ {C_{12} + C_{21} + A_{21} } }\right) \left(\vphantom{ {t_0 + \Delta t} }\right. \Delta \left.\vphantom{ {t_0 + \Delta t} }\right)_{i}^{} \left.\vphantom{ {C_{12} - \left( {C_{12} + C_{21} + A_{21} } \right)n\left( {t_0 + \Delta t} \right)_i } }\right] \Delta$$ =

$$\left[\vphantom{ {K_{12} - \left( {K_{21} + K_{12} } \right)n(t_0 )_i } }\right. \left(\vphantom{ {K_{21} + K_{12} } }\right. \left.\vphantom{ {K_{21} + K_{12} } }\right) \left.\vphantom{ {K_{12} - \left( {K_{21} + K_{12} } \right)n(t_0 )_i } }\right]$$

$$\left[\vphantom{ {FN_i + \sum\limits_j {FS_{ij} n\left( {t_0 + \Delta t} \right)_j } } }\right. \sum\limits_{j}^{} \left(\vphantom{ {t_0 + \Delta t} }\right. \Delta \left.\vphantom{ {t_0 + \Delta t} }\right)_{j}^{} \left.\vphantom{ {FN_i + \sum\limits_j {FS_{ij} n\left( {t_0 + \Delta t} \right)_j } } }\right]$$ (13)

The quantities in the square brackets are short hand for the following integral

$$\int \int_{{t_0 }}^{{t_0 + \Delta t}} \int_{{ - 1}}^{1} \mu \int_{0}^{\infty} \nu \phi \left(\vphantom{ \nu }\right. \nu \left.\vphantom{ \nu }\right) \left(\vphantom{ {\mu ,x,\nu ,t} }\right. \mu \nu \left.\vphantom{ {\mu ,x,\nu ,t} }\right)$$ (14)

In Eq. (13), the term containing FN_i is the contribution to the change in upper state population within zone i from photons that were present at the beginning of the time step. They have therefore explicit, numerical weights. The term containing FS_ijn(t₀ + $\Delta$t)_j is a similar contribution coming from bundles with symbolic weights that were born in zone j during the time step; n(t₀)_i is the upper level atomic population fraction in zone i at start of time step; n(t₀ + $\Delta$t)_i is the unknown upper atomic population fraction in zone i at the end of the time step; and V_i is the thickness of zone i. For further details, we refer to Ref. [4]. Although the problems we present will involve small numbers of zones, and as a result the solution of small linear systems using a direct solver, there is some concern that the size of the linear system to be solved will become intractable using direct solution methods for problems with larger zone counts. Jacobi iteration [8], or other more sophisticated techniques, do very well at solving the system of equations given in (13) by using the values in the previous time step for the initial starting point.