Number of Scatterings in Random Walks

In an optically thin sphere, photons will be scattered at most only once. Consequently, the mean number of scatterings eventually becomes equivalent to the probability that a photon undergoes scattering in the medium. The mean number of scatterings is then given by the following expression:

For a medium of large optical thickness, photons will be scattered multiple times. The displacement of a photon between the $(i-1)$th and $i$th scatterings is denoted by $\mathbf{r}_{i}$, as illustrated in Figure 1. The net displacement of the photon after $N$ scatterings is then $\mathbf{R}=\mathbf{r}_{1}+\mathbf{r}_{2}+\mathbf{r}_{3}+\cdots+\mathbf{r}_{N}$. The average net displacement is obtained by squaring $\mathbf{R}$ and then averaging it over all photons:

where the angle brackets $\left<\right>$ denote the average over all photons. The second term in the equation involves averaging the cosines of angles between the directions before and after scatterings. In the summation, $i$ and $j$ do not necessarily denote sequentially occurring scattering events. In the optically thick limit, this will vanish for any scattering type with front-back symmetry, such as isotropic, Thomson, or Rayleigh scattering. Then, only the first term will contribute to the mean net displacement, which will consequently become $N_{\rm scatt}\left<\mathbf{r}^{2}\right>$, where $N_{\rm scatt}$ and $\left<\mathbf{r}^{2}\right>$ are the mean number of scatterings and the mean square displacement for a single scattering event, respectively. The mean square displacement of a single scattering event can be calculated as follows:

The photon will escape after the symetem acquiring the net displacement of $N_{\rm scatt}\left<\mathbf{r}^{2}\right>=L^{2}$. Therefore, the mean number of scatterings undergone by photons is given by

Combining the above equations for optically thin and thick cases, we obtain an approximation that might be applicable in both optically thin and optically thick media:

Our result differs from that of Rybicki & Lightman (1986) by a factor of 2 in the optically thick limit. This difference arises from their assumption that the path lengths between scatterings are always constant, rather than following the probability distribution function (PDF) $P(\tau)=e^{-\tau}$. Their result is not self-consistent because they assumed the PDF when deriving the number of scatterings in the optically thin limit but not in the optically thick limit.

In an optically thin slab, photons will escape the system when the $z$ component of the displacement vector expected at the first scattering event is beyond the half height of the slab, i.e., when $\left|\mathbf{r}_{1}\cdot\hat{\mathbf{z}}\right|>L$. Therefore, unlike in the spherical geometry case, the condition for scattering depends on the polar angle $\theta$ of the photon direction. The probability of being scattered when a photon is emitted into a direction of $\mu=\cos\theta$ is $p_{\rm scatt}(\mu)=1-e^{-\tau_{0}/\left|\mu\right|}$ because the optical thickness from the center to the slab boundary along that direction is $\tau_{0}/\left|\mu\right|$. Therefore, the mean number of scatterings can be obtained by integrating the probability $p_{\rm scatt}(\mu)$ over $\mu$, which is uniformly distributed in the range of $-1\leq\mu\leq 1$. The resulting number of scatterings is then given as follows:

Here, $\Gamma(0,\tau_{0})\approx-\gamma-\ln\tau_{0}+\tau_{0}+\mathcal{O}(\tau_{0}^{2})$ is the upper incomplete gamma function, where $\gamma\simeq 0.57722$ is the Euler-Mascheroni constant. It is noteworthy that the number of scatterings for a slab is found to be higher than $N_{\rm scatt}\approx\tau_{0}$ obtained for a sphere. This is because the optical depth along a direction with $\theta>0$ is always higher than that at $\theta=0$.

In the optically thick limit, we need to consider the mean square of the $z$ component of the net displacement, which is given by:

Here, $r_{i}$ is the magnitude of the $i$th displacement $\mathbf{r}_{i}$ (i.e., $r_{i}=\left|\mathbf{r}_{i}\right|$), and $\mu_{i}=\cos\theta_{i}$ is the cosine of the polar angle of $\mathbf{r}_{i}$, measured from the positive $z$-axis direction, as illustrated in Figure 1. The angles $\theta_{i}$ and $\theta_{j}$ are independent, and thus the second term will vanish if the scattering is front-back symmetric. The first term becomes $N_{\rm scatt}\left<r^{2}\mu^{2}\right>$ if photons undergo scattering $N_{\rm scatt}$ times before escaping. Here, $\left<r^{2}\mu^{2}\right>$ represents the mean squared $z$ coordinate of the displacement vector between scattering events. The radiation field would be isotropic in the optically thick limit, even when the scattering is not isotropic, if it is front-back symmetric. In that case, the polar angle $\mu$ will be uniformly distributed, and $\left<r^{2}\mu^{2}\right>$ can be estimated as follows:

As a result, the mean number of scatterings can be obtained from the condition $N_{\rm scatt}\left<r^{2}\mu^{2}\right>=L^{2}$, as follows:

As for the spherical geometry, it is tempting to add the equations for the optically thin and thick case. However, it was found that the equation for the optically thin slab is applicable only to much lower optical thickness ($\tau_{0}\lesssim 5\times 10^{-2}$) compared to that for a sphere. This arise because the optical thickness of a slab along a direction parallel to the $xy$ plane is, in fact, infinite. The formula was also found to significantly deviate from the simulation results at intermediate optical depths. Therefore, instead of simply combining the formulas for optically thin and thick limits, we made slight modifications to Equation (6) for the optically thin case and obtained an approximate formula. The final approximate formula is given as follows:

In the next section, we demonstrate that this adjusted formula reproduces the simulation results fairly well.

The Monte Carlo simulation algorithms adopted in this study are similar to those described in Seon & Kim (2020) and Choe & Lee (2023). Photons are traced from the center, where they are generated, until they escape from the scattering region. The emitting source is assumed to be monochromatic and isotropic. The scattering cross section is assumed to be independent of the wavelength, akin to the Thomson scattering process. The optical depth $\tau$ traveled by a photon between scattering events is directly proportional to the physical path length $\left|\mathbf{r}_{i}\right|$.

Simulations were performed using two different approaches: with and without adopting the first forced scattering algorithm. In the first approach, the optical depth traveled by a photon before experiencing the next scattering event is randomly chosen as follows:

where $\xi$ is a uniformly distributed random number between 0 and 1. In this approach, the photon weight is always $w=1$. In the second approach, to minimize the inefficiency inherent in simple Monte Carlo simulations in the optically thin limit, the forced scattering algorithm is adopted for the first scattering event for every photon (e.g., Cashwell & Everett, 1959; Witt, 1977; Baes et al., 2011; Seon, 2023). The optical depth traveled between scattering events is randomly chosen as follows:

where $\xi$ is a uniform random number, and $w$ is the photon weight. The photon weight is $w=1-e^{-\tau_{\rm max}}$ for the first scattering event, and $w=1$ for subsequent scatterings. Here, $\tau_{\rm max}=\tau_{0}$ for a sphere, while for a slab, $\tau_{\rm max}=\tau_{0}/\left|\mu\right|$ when the photon is emitted along a direction corresponding to $\mu$. The results obtained using both approaches were found to agree very well. In this paper, the results obtained using the first approach are presented.

The next scattering location from the current position $\mathbf{r}$ is determined by $\mathbf{r}^{\prime}=\mathbf{r}+\tau\ell\mathbf{k}$, where $\mathbf{k}$ is the current propagation direction vector of the photon before the scattering. Because the scattering medium is assumed to have a constant density in this study, distances are measured in units of the mean free path (i.e., $\ell=1$). After determining the next scattering location, the radial distance of the photon is compared with the radius in the case of spherical geometry or the $z$ coordinate is compared with the height of the medium in the case of slab geometry. The comparison determines the fate of the photon – whether it escapes the system or continues scattering.

If the photon does not escape, the next scattering propagation vector $\mathbf{k}^{\prime}$ is chosen by sampling the scattering polar angle $\vartheta$ and azimuthal angle $\varphi$ relative to $\mathbf{k}$. The scattering angles are defined in a coordinate system formed by three unit vectors: $\mathbf{e}_{z}=\mathbf{k}$, $\mathbf{e}_{y}=\hat{\mathbf{z}}\times\mathbf{k}/\left|\hat{\mathbf{z}}\times\mathbf{k}\right|$, and $\mathbf{e}_{x}=\mathbf{e}_{y}\times\mathbf{e}_{z}$, provided $\mathbf{k}$ is not parallel to $\hat{\mathbf{z}}$ (i.e., $\hat{\mathbf{z}}\times\mathbf{k}\neq 0$). If $\mathbf{k}$ is parallel to $\hat{\mathbf{z}}$, the three basis vectors are chosen to be $\mathbf{e}_{x}=\hat{\mathbf{x}}$, $\mathbf{e}_{y}=\hat{\mathbf{y}}$, and $\mathbf{e}_{z}=\hat{\mathbf{z}}$. Once the scattering angles ($\vartheta$, $\varphi$) were sampled according to an appropriate PDF, the next propagation vector is obtained as follows:

where ${\text{\textmu}}=\cos\vartheta$. The formulas given in Pozdnyakov et al. (1983) and Seon (2009) are equivalent to the above equation, except for the definition of $\varphi$, which is modified as $\varphi\rightarrow\varphi-\pi/2$. In the case of isotropic scattering, ${\text{\textmu}}=2\xi_{1}-1$ and $\varphi=2\pi\xi_{2}$ for two independent, uniformly distributed random numbers $\xi_{1}$ and $\xi_{2}$. In Thomson scattering, the scattering angle $\vartheta$ is sampled according to the distribution function proportional to $1+\cos^{2}\vartheta$. For Thomson scattering, µ is chosen as follow for a uniform random number $\xi$:

This equation is derived in Seon (2006) and further discussed in Seon & Kim (2020).

The number of scatterings for each photon is calculated by adding the photon weight $w$ every time it undergoes scattering. The process continues until the photon escapes. Finally, the mean number of scatterings is calculated by averaging the results obtained for a large number of photons. The number of photons used in this study ranges from $N_{\rm photons}=10^{4}$ to $10^{8}$, depending on the medium’s optical depth $\tau_{0}$. For optically thin cases ($\tau_{0}\leq 10$), $10^{7}$ or $10^{8}$ photons were used, while a relatively small number of photons were adopted for optically thick cases.

Figure 2 presents the results of the Monte Carlo RT simulations conducted in a sphere (left panel) and in a slab (right panel). In the upper panels, the results obtained by assuming the isotropic scattering are compared with the analytical approximations, Equations (5) and (10). Additionally, the formula proposed by Rybicki & Lightman (1986) is also compared. In the figure, the black lines with circles represent the simulation results, while the red lines indicate the analytic approximations. The red solid line represents our formula, while the red dashed line represents the formula from Rybicki & Lightman (1986). As can be seen in the figure, the simulation results are well reproduced using our formulas in both a sphere and a slab. On the contrary, in optically thick cases, the formula of Rybicki & Lightman (1986) overpredicts the simulation for a sphere by a factor of 2, while it underpredicts those for a slab by a factor of 1.5. In optically thin cases, it significantly underpredicts the results for a slab. It is also evident that, for a given optical thickness $\tau_{0}$, the number of scatterings in a slab is always higher than that in a sphere, as the optical thickness in a slab depends on the photon propagation direction and becomes very high in the $xy$ plane direction.

The red lines in the lower panels present the relative discrepancies between the approximate formulas and the simulations where the isotropic scattering is adopted. The figures illustrate that the approximate formulas reasonably account for the simulation results within approximately 9% for a sphere and 13% for a slab. It is evident that the accuracy of the approximate formula for a sphere is better than that for a slab. This arises from the fact the optical thickness toward the direction near the $xy$ plane in a slab can be very large even in the case of $\tau_{0}<1$.

The lower panels in Figure 2 also display the relative differences between the simulation results for isotropic scattering and Thomson scattering in the blue lines with crosses. Although the differences are minor, one might recognize that the Thompson scattering results yield slightly, but systematically, larger numbers of scatterings, particularly in the case of a slab geometry with optical thicknesses of $0.1\lesssim\tau_{0}\lesssim 1$. This minor difference is attributed to the fact that, in relatively low optical thicknesses, the radiation field in an asymmetric slab is less isotropic compared to the spherical symmetric case. However, at high optical thicknesses ($\tau_{0}\gtrsim 1$), the radiation field rapidly become isotropic even in a slab geometry.