REVERSE MONTE CARLO SIMULATIONS OF LIGHT PULSE PROPAGATION IN NONHOMOGENEOUS MEDIA - PART I: THEORETICAL DEVELOPMENT

Abstract

This paper series presents a follow-up study of our previous work on the reverse Monte Carlo method solution of transient radiation transport in the homogeneous media. In this study, the method is extended to consider nonhomogeneous media, which exist in many practical problems. The transport process of ultra-short light pulse propagation inside the non-emitting, absorbing, and anisotropically scattering multi-layer media is studied. Although only onedimensional geometry is analyzed here, the method is applicable and easy to extend to multidimensional geometries. In Part I, the detailed algorithm of treating the source integration over a path length that contains non-uniform radiative properties is given. Validation and numerical simulations are given in Part II. NOMENCLATURE a absorption coefficient, m a1 linear anisotropic scattering phase function coefficient c propagation speed of radiation transport in the medium, m/s g pulse spatial and/or temporal shape function I radiation intensity, W/msr I radiation intensity at boundary z = 0 o Li the i-th layer thickness, m L total slab thickness, m o N number of samplings or energy bundles r radial coordinate, m r position vector of a location (x, y, z) in space s geometric path length, m ŝ unit vector along a given direction t time, s x, y, z rectangular coordinates, m Greek symbols κ extinction coefficient, κ = a + σ, m σ scattering coefficient, m μz direction cosine in z-direction τi optical thickness of the i-th layer, κiLi Φ scattering phase function Ω solid angle, sr ω scattering albedo Superscripts ' dummy variables '' dummy variables ^ vector ⊥ direction normal to z-axis Subscripts abs absorption d direct attenuation component; detector location ext extinction k path length index L lower limit of path length segment lk integration n index of N samplings o collimated pulse direction p pulse s multiply scattered or diffuse component seg path length segment, i.e., lk U upper limit of path length segment lk integration w wall or surface position INTRODUCTION The transient radiative process has diverse applications in astrophysics, thermal systems, biomedical imaging, remote sensing, etc. Various analytical models have been developed over the years, from the earlier diffusion approximation and it's variations, to the discrete ordinates method and the recently developed, rigorous integral equation models. Kumar, et al. [1] described several commonly used deterministic models. However, all the deterministic models have various computational, model fidelity, and scaling issues. For example, although diffusion approximations are computationally efficient, it failed to predict the photon behavior in the initial transient, incl. the propagation speed. The commonly used discrete ordinates method is relatively efficient in computational time. Nevertheless, the large memory requirement and the communication overhead make it difficulty to deploy an effective parallel scheme for large-scale problems. The problem also appears in other differential treatments of the radiative transport equation. The integral formulation is a bit awkward in the case of reflective boundary, although it has very good parallel efficiency. On the other hand, the stochastic method of Monte Carlo (MC) has long been considered a benchmark tool that is amicable to various model fidelity requirements. At the same time, it is not considered computationally efficient. An alternative is to use the reverse Monte Carlo (RMC) method, as demonstrated in our previous work [2] and described below. This study is a continuation of the prior effort to provide simulation tools that are efficient, accurate, and scalable to parallel computing system for large, complex problems. The conventional (or forward) Monte Carlo is efficient if radiation source is confined to a small volume and/or solid angle while the reverse Monte Carlo is efficient if radiative flux or incident radiation at a small surface and/or over a small solid angle is needed. Both methods become inefficient if radiation emits from a small source and the radiative flux onto a small detector is of interest, as shown in [3]. However, in transient processes the RMC method is still preferred over MC method. Since in some cases only a portion of the temporal radiative signal is needed and RMC can provide solution at a given time or over a specific time duration when needed in the socalled time-gating measurement [4]. On the other hand, the MC will require simulation from the very beginning to the time of interest [1]. Another important distinction between MC and RMC methods is the amount of samplings needed for a statistically converged solution of a general radiation transport problem, if the source or detector size is not an issue. For example, at a given standard deviation requirement of the solution the RMC can provide the necessary amount of samplings at different temporal or spatial points (see the sampling number versus time given in [1]). A MC solution with limited number of samplings will always have some temporal or spatial regions meet the standard deviation requirement but others do not. To ensure every temporal or spatial space meet the requirement, then more energy bundles will have to emit from the source. In MC simulations, the incident radiation or emission at any given position or time is proportional to the number of bundles arrive and received. The region or time that has larger statistically fluctuation is usually of lower magnitude or fewer arriving bundles. Therefore, only a smaller fraction of those newly emitted bundles will reach the needed region or time. The majority will simply go to the temporal or spatial spaces that have already met the standard deviation requirement. Overall, the total number of samplings needed for RMC is smaller. In some other cases, the need for statistically converged solutions with MC method requires huge number of samplings (greater than 10 is not unusual, for example, see [5]). The problem is further complicated if a good pseudo-random number generator was not used [6]. On the contrary, the reverse method is more suitable in general problems as accurate solutions can be obtained with the right amount of samplings. Since the total number of samplings needed for RMC is smaller, the requirement of good quality pseudo random number generators is less stringent. Although there are limited numbers of RMC simulations of the pulse propagation inside the scattering media, to the authors' knowledge, the detailed algorithm discussion is not available in the literature. A general algorithm to deal with nonhomogeneous media is especially needed. It is hoped this paper will provide an efficient, accurate, and scalable simulation tool for further study. RADIATION TRANSPORT MODEL OF LIGHT PULSE PROPAGATION The problem under consideration is a collimated, pulsed irradiation on the top surface of the multilayer slab medium. The light pulse is ultra-short, in the picoor femto-second order. The geometry is shown in Fig. 1.