{'fraction_non_alphanumeric': 0.02002172517651706, 'fraction_numerical': 0.00461660000987508, 'mean_word_length': 5.221317769927815, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Radiation trapping in 1D using the Markov chain formalism: A computational physics project

A computational model study for complete frequency redistribution linear incoherent two-level atomic radiation trapping in optically dense media using the multiple scattering representation is presented. This model study discuss at length the influence of the spectral distributions, overall opacity and emission quantum yield to trapping distorted ensemble quantities stressing physical insight and with a non-specialist audience in mind. Macroscopic reemission yield, lifetime, steady state spectra and spatial distributions are calculated as a function of intrinsic emission yield, opacity and external excitation mode for Doppler, Lorentz and Voigt lineshapes. The work could constitute the basis for a final undergraduate or beginning graduate project in computational physics instruction and implements the analytical developments of the previous instalment of this contribution.


Introduction
In optically thick media, electronic excitation energy can undergo several reabsorption and reemission events before either escaping to the exterior or being converted in thermal energy by means of collisional deactivation. This resonant radiation trapping is important in areas as diverse as stellar atmospheres [1], plasmas and atomic vapours luminescence [2], terrestrial atmosphere and ocean optics [3], molecular luminescence [4], infrared radiative transfer [5] and cold atoms [6].
The starting point of the majority of the incoherent radiation trapping models is the Holstein-Biberman equation which is a Boltzmann-type integro-differential equation describing the spatial and temporal evolution of the excited state number density. The previous instalment [7] of this contribution outlined the two alternative ansatze commonly used to obtain solutions for the classical trapping problem, Holstein's original exponential mode expansion and the so called Multiple Scattering Representation (MSR). The MSR solution was given a simple stochastic formulation and trapping dependent quantities (overall relaxation parameters such as ensemble emission yield and lifetime, time-resolved and steady-state spatial distributions as well as spectra) were calculated with the Holstein fundamental mode singled out. This instalment will now use a simple Markov chain algorithm to quantify incoherent trapping in a computational model study for two-level atomic models. It could easily be adapted into a computational physics project valuable in the context of atomic or computational physics instruction for final undergraduate and beginning graduate students. With this objective in mind, trapping dependent quantities are estimated in a unidimensional geometry for a single line with Doppler, Lorentz and Voigt spectral distributions. The excitation relaxation dynamics is considered for conditions mimicking electron impact as well as photoexcitation. A thorough discussion of the conditions for which the use of Holstein's fundamental mode alone is a tolerable approximation is included. Section 2 discusses the way the dynamics of incoherent trapping is taken into account within the framework of the MSR. It gives explicit expressions for the overall relaxation parameters as well as steady-state spectra and spatial distribution summarizing the previous instalment of this work. Section 3 explains the rationale as well as critical implementation details of the Markov stochastic simulation algorithm, the one chosen in this model study. Section 4 presents the results and their discussion at length. Particular emphasis is placed in the discussion of the physical implications of the spectral distribution as well as quantifying the relative contribution of the fundamental mode. Section 5 summarizes the main conclusions. Finally, the Appendices show implementation fine details and add some possible routes for complementing the work presented here.

Dynamics of incoherent trapping
Radiation trapping studies should be cast in dimensionless coordinates since this increases computational efficiency and, more important, defines characteristic scales or universal conditions. The quantities most directly amenable to define characteristic scales in trapping are time, distance and optical frequency. The scaled time is t = Γt , where Γ is the global deactivation rate constant. The dimensionless distance is sometimes called the opacity or optical density and can be defined, along a given pathlength l and for homogeneously distributed species, as r = k (x) /Φ (0) = k 0 Φ (x) /Φ (0). x is used to represent the optical frequency (see below). k (x) is the single line monochromatic opacity (k 0 is the corresponding center-of-line value) and the absorption lineshape is given by the normalized spectral distribution Φ (x) (so that +∞ −∞ Φ (x) dx = 1). For two-level atomic models the intrinsic, trapping undistorted, spectra can be written as a function of a dimensionless optical frequency defined as x = ν−ν 0 ∆ν D . This is a normalized difference to the center of line frequency, where ∆ν D stands for the FWHM of the Doppler distribution at each given temperature. For two-level atomic models, it is common to describe the absorption lineshape by The Doppler distribution allow us to single out the pure Doppler broadening from the other broadening mechanisms while the Lorentz and Voigt's distributions are the ones to be used in pure radiation damping and combined radiation and collision broadening conditions, respectively. In this last case, a = ln (2) ∆ν L ∆ν D is the Voigt characteristic width, the relative Lorentz over Doppler spectral width and implicitly dependent upon both temperature and vapour pressure via the dependence on collisional cross-section values.
In the MSR ansatz for linear incoherent trapping with negligible time-of-flight for in-transit radiation, the spatial and temporal relaxation dynamics for excitation is given by n (r, t) = n a n p n (r) g n (t) , where n stands for the generation number of excited species (primordial excitation creates first generation, the trapping of this generation's reemission creates second generation and so forth; one can envisage each generation as the result of n − 1 previous scattering events of resonant radiation), a n is the population efficiency for each generation, and p n (r) and g n (t) are the (normalized) spatial and temporal excited species distributions.
From a practical point of view, the most important quantities are the macroscopic ensemble relaxation parameters, overall reemission yield φ and mean scaled lifetime τ , and the steady state emission spectrum I SS (x). These were derived in a previous contribution and are given by [7]: and In these, m stands for the first generation number that can be considered nonchanging (subscript nc is a remainder for nonchanging; see [7]) and α n ≡ a n+1 an is the mean trapping or reabsorption probability. Finally, is the mean escape probability in the direction defined by the solid angle Ω which, for a left escape from a unidimensional geometry (see next section), is just if one decides to start the opacity scale on the left side of the cell. It is also informative to know the steady-state spatial distribution, which can be cast as [7]: n SS (r) = m−1 n=1 a n p n (r) + anc 1−αnc p nc (r) m−1 n=1 a n + anc In the above expressions it is important to recognize that the trapping dynamics can be factored out into a generation varying part and another corresponding to generations that have the same spatial distribution (the fundamental mode) and the nonchanging part can be expressed explicitly as an analytical sum. The dynamics for this nonchanging part corresponds only to an attenuation of overall excitation in going from one generation to the next by a fixed α nc factor and, as a result, the contribution of this part corresponds to a monoexponential relaxation with a trapping dependent effective decay constant.

Markov stochastic algorithm
The stochastic formulation of the Multiple Scattering Representation (MSR) allows to quantify the trapping influence on observables essentially trough the estimation of mean reabsorption and escape probabilities. In this respect, it is simpler and more amenable to discussion at an elementary level than the alternative Holstein multiexponential expansion. The MSR numerical implementation is straightforward and can be easily discussed at a level that stresses physical intuition without actually getting bogged down in all the technicalities of the troublesome estimation of Holstein modes other than the fundamental. We have furthermore decided to use an unidimensional geometry, driven by the motivation to do a simple mimic of a cylindrical tube of a discharge fluorescence lamp and by the desire to keep computational detail and power adequate to the proposed audience of final undergraduate and beginning graduate students. In this 1D case the opacity scale comes naturally as the opacity along the cylindrical axis.
The model used is a two-level single line atomic model with Doppler, Lorentz or Voigt spectral distributions and particular attention is paid to the influence of the lineshape on the trapping efficiency. The Voigt case will illustrate the fact that a continuous variation of the characteristic width parameter will map the Doppler into the Lorentz distributions by changing the relative importance of the Lorentz-like wings over the Doppler-like core of the distribution (figure 1). Accordingly, the ability to compute the Voigt line to machine precision is mandatory and Appendix A gives some implementation details in this respect.
Complete frequency redistribution conditions are used [2], which means that the number of collisions during the lifetime of excited atoms is sufficiently high to render the reemitted photon's frequency completely uncorrelated with the frequency of the previously absorbed photon. In these conditions, the absorption and emission lineshapes coincide and the jump length distribution of the excitation random trajectory is independent of past history. This makes the formalism of Markov processes [8] especially adequate. Its rationale for the MSR implementation can be cast in the following way. The 1D cell is divided into several bins, each corresponding to a pure state of the system and characterized by a mean probability that the excitation resides in that state. The system dynamics corresponds to the evolution of the probability of excitation being inside each state. The stochastic process is completely specified by (i) a column vector with the (normalized) spatial probability distribution of the first generation species, p 1 = [p i 1 ], and (ii) a transition matrix P = [p ij ], whose entries are the one-step transition probabilities between states i and j. For complete frequency redistribution, there is an absence of memory effects (homogeneous chain) meaning that the transition probability between individual states depends only upon their relative opacity distance (it is independent of generation number and thus computed only once).  The spatial distribution functions for all the generations are calculated from the previous generation by: The sample cell is divided into h-length bins and the transition matrix elements are therefore given for a 1D geometry by [9]: corresponding to the Beer-Lambert law weighted according with the emission lineshape for an individual reemission-reabsorption (scattering) event. The integration takes into account all the possible emission frequencies, 1/2 is the left or right emission direction probability for a 1D geometry and it was assumed that the bin width is controlled in order to attain a satisfactory precision.
The complete specification of the Markov process is achieved once one specifies the initial spatial distribution p 1 . Two different cases were considered, one for homogeneous initial excitation (trivial) and another mimicking photoexcitation with the reabsorption undistorted line. For photoexcitation from the left side of the 1D cell, which afterwards must be properly normalized.
The binning of the spatial excitation distributions corresponds to the substitution of a continuous distribution for its discretized version, effectively transforming a continuous process into a discrete one realized in a lattice. It is therefore the Markov equivalent of a random walk defined over a regular spaced lattice [10]. The bin width (or the number of cells) is the critical parameter for the Markov algorithm. We have conducted several tests and found advisable to have a maximum bin size of 0.05 in an opacity scale. Otherwise, numerical artifacts associated with substituting the actual excitation migration for the jump between the mean coordinates of each bin of sample cell could exist. These were found to be the more important the higher the overall opacity.
We have mentioned in the beginning of this section that the Markov implementation of the MSR model aims at estimating mean reabsorption and escape probabilities. The reabsorption probabilities are estimated with the following procedure. Each time (8) is used, the fraction of the excitation remaining inside sample cell gives the n th mean reabsorption probability α T n . The excitation column vector is then (re)normalized and the process repeated. From the values of this parameter for each generation, the trapping population efficiency is a T n = n−1 n=1 α T n . This procedure is the equivalent of an importance sampling method (see Appendix B), in which it is assumed a unit intrinsic reemission yield in (9). The influence of the actual value of this yield (φ 0 = Γr Γr+Γnr , the ratio of radiative over radiative plus nonradiative relaxation rate constants) is introduced analytically. Using the notation developed in the previous instalment, the Markov algorithm directly estimates a T n , each generation population efficiency due to trapping (and geometry) alone, and the actual population efficiencies were then given by a n = a T n φ n−1 0 [7]. As for the escape probabilities, we have chosen the following implementation of (6). The monochromatic left escape probability is obtained from where p n is the spatial distribution and Q (x) is the escape matrix, whose entries are finally which comes directly from Beer-Lambert law.
Equation (8) is equivalent to linear response theory and therefore the evolution of the spatial excitation is given as a convolution integral between the excitation profile p n with the delta response function given by (9). This renders the Markov approach more efficient since this convolution can be easily made using FFT algorithms [11] paying attention to zero pad to double size the column vector containing the excitation distribution in order to avoid wrap-around effects due to the cyclic convolution [12]. The speed-up factors could rise up to several orders of magnitude (roughly 50 to 400 times for a number of Markov states of 2 000 to 200 000) making the FFT convolution the recommended implementation of the Markov algorithm. Figure 2 shows the overall relaxation parameters for an initial homogeneous excitation and an intrinsic quantum yield of φ 0 = 0.9 as a function of opacity for Doppler, Voigt and Lorentz lineshapes. The higher the opacity the more important trapping is, with the following implications: (i) an increase of the mean relaxation time (equivalent to a mean number of scattering events before escape) and (ii) a decrease of the ensemble reemission yield (the fraction of original excitation that eventually comes out; an increased importance of trapping translates into additional possibilities of nonradiative relaxation). The Voigt continuous transition from Doppler into Lorentz is evident as well as the relative importance of core and wings of the distributions. The higher the Lorentz character of the spectra the higher the weight of the wings and the higher the escape probabilities (lower mean lifetime and higher macroscopic emission yield). Figures 3 and 4 show the mean scaled lifetime and reemission yield for homogeneous initial excitation for the Doppler and Lorentz limiting spectral distributions and for several values of the intrinsic quantum yield. Three main conclusions can be drawn from these results. First of all, the two most important parameters controlling trapping efficiency are the spectral distribution and the value of φ 0 . The higher the overall opacity the more difficult is the escape of radiation for Doppler-like distributions and the more important the escape from the Lorentz-like wings of the distribution. Second, trapping for Doppler-like distributions is much more efficient since, especially in high opacity cases, the escape of excitation at optical frequencies far enough from the line center frequency is reduced due to the extremely small probability of reemission at those frequencies (for unit reemission yield, the lifetime for the higher opacity is about 200 for Doppler and only about 15 for Lorentz). Finally, under conditions rendering trapping efficient, the φ 0 value is of paramount importance; for unit intrinsic reemission probability all the excitation will eventually come out (thus giving a simple check for  consistency of computation) but, as soon as φ 0 is smaller than one, each new scattering event gives the excitation another chance of thermal degradation (in (2) each generation contribution is intrinsically dependent upon φ n−1 0 ). Note in these figures that there is a strong dependence of the relaxation parameters with the φ 0 value, especially for the more trapping influenced Doppler case.

Ensemble relaxation
All these conclusions are important in the discussion of atomic vapours ensembles for lighting applications, either electric discharge lamps or plasma display panel (PDP) devices. Better performance is achieved with higher macroscopic reemission yields. On top of that, the increase of the overall opacity is in principle desirable since this is related with the increase of the number of excited species. In a crude first order approximation, one can assume the overall lamp efficiency to be directly proportional to the product of φ times the overall opacity (directly proportional to initial excitation density): The actual behaviour of a lamp or a PDP can be quite involved since an increase in opacity means an increase of partial vapour pressure (the external dimensions of the device are fixed) and this induces several changes whose influence on the overall performance can be contradictory. The higher opacity means higher light throughput as long as the increase in the trapping efficiency does not substantially increase thermal degradation. The higher the opacity the more efficient the trapping (higher recapture probabilities) but, on the other hand, the resulting higher collisional effective rate constants could give rise to a smaller intrinsic φ 0 and render the Lorentz distribution only approximately valid. A reduced φ 0 value and a shift of the Voigt distribution towards a more Lorentz-like profile tend by themselves to decrease trappping efficiency and increase light throughput. To have a simple idea of the effect, and due to the paramount importance of the reemission quantum yield φ 0 , a series of results were made for both limiting Doppler and Lorentz distributions with the radiative quantum yield given by φ 0 = Γr Γr+Γq , where the quenching rate constant was assumed in a first order approximation to be linear with the cell opacity (Γ q ≡ k r, with the numerical values Γ r = 10 7 s −1 and k = 2.5 × 10 4 s −1 ). This corresponds to the well known Stern-Volmer equation for dynamical quenching by binary collisions for unitary intrinsic radiative yield in the absence of collisions. The results are shown in figure 5 which is judged more representative of the actual lamp behaviour than the results in figures 3 and 4. Figure 5 shows that, ultimately, a delicate balance will dictate the best operation conditions which manifest themselves in the peaks of the Ψ values. Of course, the results are very approximate since the assumed functional dependence of φ 0 is only approximate and, even within the CFR two-level Lorentzian assumption, one should use an opacity dependent characteristic width of the Voigt distribution. Nevertheless, figure 5 allows  the discussion of the qualitative behaviour emphasizing physical insight without the additional burden of fine grained details. It shows how critical the spectral distribution shape and the quantum reemission yield are. For φ 0 values sufficiently close to one, an increase in opacity corresponds to an increase in lighting efficiency due to the increased initial excitation number density. Due to the strong dependence of trapping on the φ 0 value, as soon as this value starts to be significantly smaller than one, trapping means an much increased importance of thermal degradation (compare the difference between the quantum and the ensemble yield or the reduced overall relaxation lifetime in the upper part of figure 5). From some point onwards this will be more important than the increase in initial excitation due to a higher vapour pressure, which originates optimal operation conditions, giving the best possible lighting efficiency. Figure 5 also shows that the Doppler distribution is associated with smaller throughput in lighting applications when compared with the Lorentz case due to the step reduction of the ensemble reemission yield with the increase of the overall opacity. This is of course related with the use of an inert gas filling to render collisions more important (increasing the Lorentz character of the spectral distribution and reducing trapping efficiency) in fluorescence lamps.
Finally, figure 6 shows the predicted ensemble relaxation parameters for both homogeneous and photoexcitation as a function of overall opacity. Up to opacities of the order of 10 no significant difference exists (photoexcitation is able to penetrate well deep into sample cell). But, for higher opacities, the importance of trapping continues to increase indefinitely for homogeneus excitation while it levels off for photoexcitation, a point to be revisited in section 4.3 when discussing spatial distribution functions. Figure 7 shows the estimated normalized spectral distribution in steady-state conditions for Doppler, Lorentz and Voigt lineshapes for both primary homogeneous and photoexcitation. The motivation for the homogeneous case is the excitation along the axis of a fluorescence lamp for lighting applications. It shows the well known selfreversal of spectral lines due to the higher attenuation of core optical frequencies. For photoexcitation there is a balance between reduced penetration of external excitation and higher attenuation at core frequencies which dictates a flattening of the spectra near the line center (of course, for left wall photoexcitation and right wall detection there is a self-reversal higher than the one for homogeneous excitation; not shown). In both cases, there is a considerable broadening of the detected spectra and the Voigt distribution has an intermediate character between core Doppler-like and wings Lorentz-like.     Figure 8 shows the spatial distributions for both limiting cases of Doppler and Lorentz distributions. This figure shows that the fundamental mode spatial distribution (limiting case for a relaxed, non-changing spatial distribution and thus independent upon the original excitation) could only give a reasonable approximation of the steady-state distribution for the homogeneous excitation case; for photoexcitation it is more convenient to choose the spatial distribution of the first generation species as a first approximation to the overall distribution. This illustrates the well known procedure of, whenever approximating the actual trapping dependent behaviour by the monoexponenial fundamental mode (easier to obtain by a variational procedure or given by Holstein's asymptotic approximations), design the experimental setup to mimic as much as possible the fundamental mode spatial distribution (symmetrical and well spread into the bulk of sample cell) with the external excitation. This can be accomplished with photoexcitation of high opacity samples using strongly detuned external radiation. Figure 8 shows some of the difficulties of quantifying trapping simply by using the fundamental mode, a point further illustrated in table 1. Several conclusions can be draw from its data: (i) Doppler distributions render trapping much more efficient and thus its fundamental mode contribution is always much higher than for the Lorentz case, (ii) the spatial spreading for Doppler is smaller, giving rise to higher generation number for the fundamental mode, (iii) the use of the fundamental mode alone for Lorentz distributions (and therefore, albeit with a lesser degree, for Voigt) is never justified and (iv) to approximate the actual behaviour for photoexcitation to the fundamental mode is never justifiable.

Spatial distribution
Two additional points related with the common practice of using only the fundamental mode to take into account trapping distortions should be stressed out. First, the fundamental mode is the slowest decaying possible and is located well (and symmetrical) into sample cell. To substitute the whole of the ensemble dynamics for the fundamental mode alone will always overestimate the lifetime and underestimate the reemission yield (spatial distribution giving the highest possible trapping efficiency) thus introducing a systematic error. Secondly, the use of the fundamental mode alone is too often misunderstood with the use of the asymptotic approximations proposed by Holstein [13], only valid in the high opacity limit and for ideal geometries [2]. The MSR has a clear cut advantage in this respect since it allows an easy estimation of the fundamental mode as the one corresponding to a nonchanging spatial distribution, irrespective of opacity and geometry.
The spatial distribution functions presented in this section draw some further insight into the previous results of figure 6. The leveling of the lifetime and reemission  yield at higher opacities for the photoabsorption case corresponds to absorption of external excitation complete within a layer smaller than the overall opacity making the ensemble relaxation effectively insensitive to the overall opacity. This approaches well the conditions of semi-infinite geometry, under the time scale of ensemble complete relaxation.

Conclusions
This instalment presented a unidimensional computational model study for complete frequency redistribution linear incoherent atomic radiation trapping illustrating the numerical implementation of the stochastic model developed previously. It illustrates the advantages of the multiple scattering representation (MSR) over the Holstein expansion, based on physical insight and computation feasibility at an elementary level. Holstein's ansatz has significant shortcomings when compared to the equivalent alternative of MSR: (i) Holstein spatial modes are unphysical except for the fundamental, (ii) their estimation is computationally much more troublesome than the simple algorithms used in this work, (iii) the wide spread use of original Holstein expressions for the fundamental mode are only valid in the asymptotic limit of high opacities while MSR allow an easy estimation of the fundamental at any opacity value and (iv) the higher Holstein modes are difficult to obtain while for Lorentz-like spectral distributions we found that their contribution must be always taken into account (the fundamental mode contribution to ensemble relaxation being always small; higher Holstein modes correspond in the MSR language to small number generations and are easy to obtain with the stochastic formulation presented). The dependence of the ensemble reemission yield and lifetime, relative efficiency for lighting applications, steady-state spectral and spatial distributions on quantum yield, opacity and homogeneous or external photoexcitation are discussed at length for the Doppler, Voigt and Lorentz lineshapes. We quantify the contribution of the nonchanging fundamental mode and found troublesome using uniquely this mode for Voigt and Lorentz like spectra. The results should appeal to a broad audience and provide insight in a wide range of more realistic situations in an atomic as well as in an astrophysical context.
Several possible developments of the general framework presented in this work are shown in Appendices B to D. In the first, we give an outlook of the implementation details that we have found critical for a Monte Carlo simulation alternative of the Markov algorithm. The Monte Carlo and Markov approaches constitute in fact two general purpose algorithms for incoherent radiation propagation problems. Each has in own advantages and shortcomings. The Monte Carlo constitute a simulation of particle like trajectories while the Markov model quantifies the evolution of mean probabilities. We advocate the second alternative in all but the more demanding cases (detailed 3D geometries, realistic multi-level atomic models, partial frequency redistribution and polarization dependent radiation transport, issues not addressed here). In Appendix C, we give a brief description of a more realistic but nevertheless still unidimensional Markov implementation of a radiation transport model, especially appealing for planeparallel stellar atmosphere theories. Finally, in Appendix D we give an analytical modification of the Markov algorithm for multiple specular reflections in boundaries. The appendices could be useful for more advanced projects.
vapours. However, given the current desktop computing capabilities, the numerical (careful) direct integration of the defining equation is perfectly adequate and the use of approximations to reduce computation time is not justifiable any more. The difficulty in the direct numerical integration cames from the behaviour of the integrand function: it differs from zero over two width scales, a broad scale centered in zero (corresponding to the exponential term in the numerator of the integrand) and a much narrower one centered in a frequency corresponding to the Voigt frequency (associated with the difference term in the denominator). The integration domain should thus be broken into smaller domains and an automatic adaptative integration algorithm should be used in each subdomain always starting at the integrand function maximum and with a initial stepsize adapted to the local scale of variation of the integrand [12]. We have used the 400 central frequencies for the central broad feature and a 0.4 frequency width for the floating narrow peak. Integrations further away from the central core were analytically mapped from an infinite to a finite integration range. In order to decrease the time for repetitive Voigt functions evaluations, the Voigt distribution was previously computed in a given table of frequency values and, whenever necessary, cubic spline interpolated [12]. We used a linear scale in the core (frequency range up to 100 with a 5 × 10 −2 spacing) and a log scale in the wings (frequency range from 100 to 10 8 with 1 × 10 −3 log 10 spacing). Natural cubic spline was not necessary since the derivatives of the Voigt distribution in the end points are analytical.

Appendix B. Monte Carlo simulation
The mean reabsorption and escape probabilities can be estimated either with a Markov chain algorithm or with Monte Carlo simulation mimic of the experiment. We will outline the computational details we found critical for the Monte Carlo alternative and give our best advice on each method strengths and limitations and under which physical conditions the Monte Carlo is especially adequate.
The Monte Carlo (MC) method makes a direct simulation of the trapping process using particle-like trajectories for radiation in cell [14]. The initial excitation coordinate is randomly chosen corresponding to either homogeneous or external photoexcitation. The (re)emission coordinate is the same as the absorption one and, after emission, a random direction and optical frequency x must be chosen from the appropriate distributions. The photon path in cell is followed and a reabsorption coordinate drawn from the Beer-Lambert exponential distribution with absorption coefficient given by Φ (x). This should be then tested for escape from cell; if the photon escaped, another excitation trajectory should be initiated from the first generation, otherwise the simulation must continue by increasing generation number and repeating. Appropriate counters keep track of the actual number of trajectories giving rise to n th generation species (the ratio to the total number of trajectories initiated for first generation gives the population efficiencies), and the mean escape probabilities (for each generation number, excitation coordinate r and optical frequency x, the escape probabilities counters are incremented with e −Φ(x)r for left escape; compare (6)). To decrease the computation burden, an importance sampling method should be used always assuming a unit intrinsic reemission yield and after the influence of the actual value of the yield introduced analytically, as described in the main text. The spatial distribution functions for each generation are most easily obtained for a discretized cell just by keeping track of the number of excited species in each cell bin but it is important to acknowledge the fact that this binning is used only for graphical representation purposes since the mean escape probabilities are computed from the actual spatial coordinates. This is an an important difference relative to the Markov algorithm since this last case effectively makes a simulation on a lattice model from the very onset.
The complete Monte Carlo description lacks only the fine details of the implementation of the transformation method to obtain non-uniform deviates [12]. These are used in the simulation to obtain excitation coordinates and optical frequencies.
The excitation coordinates are either homogeneous for initial excitation (trivial) or drawn form the Beer-Lambert law, giving the reabsorption coordinate relative to emission point in an opacity scale for a given absorption coefficient at each optical frequency Φ (x). Since Beer law corresponds to an exponential distribution, the random deviates are analytically given by −ln(y)/Φ (x), where y is a uniform deviate. For the optical frequencies x, we must distinguish between Doppler, Lorentz and the Voigt lineshapes. The Doppler and Lorentz distributions have analytical inverses and are therefore simpler. Lorentz deviates are given by x = tan [π (y − 1/2)] while Doppler correspond to zero mean 1/2 variance normal deviates. The reference [12] only gives directly a routine for zero mean unit variance deviates but it is easy to derive a general purpose formula for zero mean, arbitrary σ variance along the same lines of the used Box-Mueller algorithm. The Voigt case is more troublesome since the transformation method must be implemented numerically. The most straightforward implementation of this is to use a cubic spline of the cumulative distribution function, truncated to a very high frequency (se have obtained good results with an upper limit of 10 8 , on a logarithmic scale in the wings).
For the simple model study of this work, the Markov algorithm has distinctive advantages over the Monte Carlo alternative since it considers directly the evolution of mean excitation probabilities. These can be obtained with a fast algorithm and this is very important for trapping due to the critical need to proper quantify the fundamental mode. The estimation of the fundamental mode in MC is much more difficult due to the characteristic slow convergence of MC estimates and to the fact that MC simulations must be made with a maximum generation number specified at start. This is especially important for high opacity cases since in this case the MC computation time can be prohibitly large and we have found that a small error in the fundamental mode can be greatly amplified in overall relaxation parameters. MC is nevertheless more versatile with respect to generalizations to either more realistic geometries and atomic models as well as to include polarization effects. The influence of geometry deserves a small note. For realistic tridimensional geometries the demands introduced by a sufficiently fine spatial discretization in the Markov procedure can compromise the practical feasibility with current desktop technology. For MC there is no significant added overhead computational complexity since the level of discretization only changes the visual resolution of the spatial distribution functions. Balancing the shortcomings of both approaches, our advice would be to use Markov whenever possible. For 3D geometries, a possible procedure would be: (i) to identify the smallest dimension in cell, (ii) to do a 1D Markov estimation of the fundamental mode generation number for that smallest opacity and (iii) finalize by MC simulating the actual geometry until that generation number (assumed to describe well the fundamental mode).

Appendix D. Specularly reflecting walls
Another case of practical interest is a vapour cell with partially reflecting walls [2]. For the unidimensional geometry and for specularly reflecting walls, the Markov algorithm is amenable to a simple modification since the absorption, transition and escape probabilities can be cast as series which have analytical representations [16].
Consider a unidimensional geometry with zero at the left and r max maximum overall opacity in the right wall. The corresponding reflectances are assumed constant and will be represented as R L and R R . The stochastic interpretation of the Beer-Lambert law gives for the absorption and escape mean probabilities, after an r optical pathlength, Φ (x) e −Φ(x)r and e −Φ(x)r , respectively. Φ (x) is just the mean probability of reemission of a photon with frequency x. With all of this in mind, the Markov transition and escape probabilities can be reformulated taking into account the possibility of multiple specular reflections of the geometry boundaries. For the mean transition between Markov states i and j, instead of (9) one can write with the integrand function given by R n L R n R e −Φ(x)rmax 2n + The first line gives the transition probability due to direct absorption prior to any reflection. The other lines give the additional absorption probability in state j after at least one reflection in the cell walls. The second and third lines correspond to emission from state i to the left while the fourth and fifth correspond to emission to the right. Finally, the second and fifth lines consider photon absorption in state j with the radiation coming from the left of this state while the third and fourth lines correspond to absorption from the right.
Rearranging, one eventually obtains Using the same procedure, the left escape probability is given by The first line is the left escape, given that the emission was is state i to the left while the second is the left escape for initial right emission.
Finally, multiple reflections can also change the initial spatial distribution in the case of external photoexcitation. For the case of photoexcitation from the left side with the undistorted resonance line, (10) should be modified into The first line is the absorption in state i of radiation propagating from left while the second is the absorption of radiation coming from the right.
For specular reflection in this unidimensional model the advantages of the Markov algorithm over the alternative Monte Carlo are even more important than for the case in which reflection in boundaries is not considered. The multiple reflections are taken into account analytically for the Markov case while in the Monte Carlo simulation the increase in computation time with the increasing importance of reflections can be very high. Nevertheless, the increase in trapping efficiency with the added possibility of specular reflection is only significant if both reflection coefficients are important. In fact, a 1D geometry in which one of the boundaries is perfectly reflective while the other has a zero reflection coefficient is equivalent to a 1D sample of twice the size of the original one.