1. Introduction

Owing to numerous randomly distributed scatterers, it is difficult to accurately describe light propagation in a scattering medium. Compared to the analytical radiative transfer equation, Monte Carlo (MC) simulation provides a simple alternative for predicting the photon flux in and out of a medium [1–5]—particularly when the medium is not uniform inside or irregular in shape. MC simulation is a standard tool that is widely used in many areas, such as biomedical engineering [6–8], fundamental physics [9–11], and atmospheric optics [12–15]. The scattered field can be coherently superposed via the electric field Monte Carlo (EMC) method [16–19] by treating each photon as a wavelet and adding a path length-determined phase to its optical field. However, without fixed particles, the EMC method deals with statistical implementations of scattering media, resulting in no difference between any two points within the medium [20]. In the temporal dimension, the polarized thermal light intensity of each pixel on a detection screen obeys negative exponential statistics, that is, it has a negative exponential probability density [21]. However, there is no stable intensity correlation between any two points on a detection plane; therefore, an unstable speckle pattern is obtained using the EMC method. In 2005, Xiong et al. proposed a particle-fixed Monte Carlo (PF-MC) [22] method for investigating the photon flux in a stationary non-uniform medium, where the step size is not determined by a statistical probability function but by the distribution of internal scatterers. In some studies, MC simulation was used to investigate the temporal—but not spatial—coherence of light after propagation in a scattering medium [23–25]. In general, scatterers are not fixed in MC methods, which involves ensemble averaging over many realizations of the medium. Speckles do not survive through ensemble averaging. Thus, a key step in achieving the simulation of coherent speckles is to choose an instant realization of the medium with fixed scatterers inside.

Compared with incoherent photon flux, coherent speckle patterns are more useful for understanding optical field propagation in random media. Unfortunately, no methods for calculation or simulation of coherent field propagation in a macro-scattering medium with a large number of scatterers have been developed, to the best of our knowledge. There are tools for computational electromagnetic analysis, such as finite difference time domain (FDTD) methods, to solve Maxwell or Helmholtz equations directly and calculate the field in a complex medium. However, large amounts of computer memory and computing power are required, which limits the size of the medium to nanometers or micrometers, generally with only a few particles [26–28]. Thus, they are not suitable for most applications in the aforementioned fields.

In this paper, we present a path sampling and integration (PASIN) method for calculating field propagation in scattering media and exit speckle patterns. The principles and workflow are introduced in Section 2. In Section 3, the results of the proposed method in different scenarios are presented. Subsection 3.1 discusses the evolution of the intensity distribution with an increased number of launched photons, which eventually converges to a stable speckle pattern. Different internal distributions of scatterers resulted in different speckle patterns. Together, the results indicated that the PASIN method can be used to calculate the internal field and exit speckle pattern of a medium. For a dynamic scattering medium, speckle patterns at different times and their temporal correlations are presented in Subsection 3.2. Using the PASIN method, the angular optical memory effect (ME) was examined in depth, and the effects of the anisotropy factor on the speckle size and appearance were investigated; the results are presented in the next two subsections. Finally, we discuss the potential impact of the proposed method and present conclusions.

2. Methods

2.1. Principle and workflow

In traditional MC methods, photons are considered energy packages [2]. After a photon impinges on a turbid medium, along its direction a random step size is sampled, and at the end of the step size, a random direction is generated according to a phase function. This procedure is repeated until the photon exits the medium or is absorbed, if the absorption effect is considered. A speckle pattern is the fingerprint of a scattering medium with a fixed internal structure. Therefore, an instant distribution of fixed scatterers should first be generated in the PASIN method. An incident photon propagates along its direction until it collides with a scatterer, after which its direction changes. Importantly, the step size here is not random; it is determined by the photon direction and the distribution of the scatterers. Here, the photon is considered a localized wave packet. The possible propagation paths are sampled with the corresponding weights. The final speckle pattern is determined by coherent superposition of all the sampled paths. There are two major distinctions from the conventional MC methods. First, the particles are fixed. The cohort of the scattering paths is determined, which differs from a conventional MC simulation. Second, the propagation paths in the medium are tracked in the conventional MC simulation; however, the path lengths have a trivial influence on the intensity pattern, because light is added incoherently. In our method, sampled photons with different phases determined by path lengths are coherently added to obtain an interferogram, that is, the speckle pattern. Figs. 1(a) and (b) present a comparison between the MC method and the proposed method.

 

Fig. 1. Principle and workflow of the PASIN method. (a) Schematic of photon propagation in conventional MC methods. Initially, there are no scatterers in the medium. After a photon is launched, the step size is sampled, and at the end of the step, a scatterer appears, which is an ideal point with zero diameter. Immediately after the scattering event, the scatterer disappears and has no influence on other paths unless another sampled step ends at the same position. The procedure continues until the photon is absorbed or exits the medium. The end of the step can be anywhere; that is, scatterers can be at any position inside the medium. The scatterers are virtual points, not real particles. (b) Schematic of photon propagation in the PASIN method. All the photons propagate in a medium with the same fixed scatterers, which are real particles. (c) Configuration of the simulation. A laser beam with a diameter of 1 mm illuminates a scattering layer, and a speckle pattern is acquired on a screen behind the layer. (d) Workflow of the PASIN method.

Download Full Size | PDF

The workflow of the PASIN method is illustrated in Fig. 1(d). First, a scattering sample is created with spherical particles randomly distributed inside it. The internal distribution of the sample is fixed once it is generated. To increase the computation speed, the sample is segmented into numerous cells. Second, a photon is launched onto the sample. If a particle exists within a cell in the path of the photon, scattering occurs; the position (x, y, z) and direction (µ_x, µ_y, µ_z) of the photon are updated, and the passed path length is accumulated. Otherwise, the search continues within the cell. Third, after searching all particles in the cell, the photon determines whether it will enter a neighboring cell or exit. If it enters a neighboring cell, the cell index is updated, and Step 2 is repeated. If it exits and is received on the detection plane, its information is recorded; otherwise, the photon is dropped. It is then determined whether this photon is the last photon. If so, the path sampling is terminated; otherwise, the photon-launching procedure is repeated. On the receiving plane, the area of interest is meshed into grids with a size smaller than the wavelength of the photon. The initial field of each grid is set as 0. Once a photon reaches a grid element, its field value is updated by adding a localized wave packet, whose phase is determined by the accumulated path length. The speckle pattern is calculated from the final field distribution.

In traditional MC methods, the Henyey–Greenstein phase function [29,30], which corresponds to the angular distribution of the intensity of scattered light, is used to obtain the deflection angles ${\kern 1pt} \theta {\kern 1pt}$ and ${\kern 1pt} \phi$. The probability distributions of $\theta {\kern 1pt}$ and ${\kern 1pt} \phi {\kern 1pt}$ can be expressed as

By coherently summing all the received scattered photons, we obtained a speckle pattern. For simplicity, medium absorption and photon polarization were not considered in this study.

The calculations were performed using a commercial supercomputer center with multiple Intel Xeon Platinum 8153 central processing units (CPUs) at 2.00 GHz. The computation speed was 1.9 × 10¹⁰ photons per core per hour. Typically, 80 CPUs were allocated to run the simulation.

2.2. Simulation parameters

To validate the method, the evolution of the intensity distribution with an increasing number of launched photons was investigated. The setup is shown in Fig. 1(c). A 532-nm laser beam with a diameter of $1\textrm{ }mm{\kern 1pt}$ was normally incident on a scattering layer. In this and the following simulations, the medium size was ${\kern 1pt} {L_x} \times {L_y} \times {L_z} = 3 \times 3 \times 0.5\textrm{ }m{m^3}$, and the cell size was${\kern 1pt} {l_x} \times {l_y} \times {l_z}{\kern 1pt} = 10 \times 10 \times 10\textrm{ }\mu {m^3}$. In the medium, the refractive index of the background was ${\kern 1pt} {n_b} = 1$, and the refractive index of the particles was ${\kern 1pt} {n_p} = 1.06$. All the particles had the same radius (${r_p} = 1\textrm{ }\mu m$), and the calculated anisotropy coefficient was ${\kern 1pt} g = 0.98$[31]. The scattering mean free path is expressed as follows:

In a dynamic scattering medium, the speckle pattern changes with time. By driving the internal particles to move according to motion laws, we simulated the dynamic scattering media and calculated the speckle pattern at each instant using the PASIN method. The time dependence of speckles with Brownian motion was investigated. In Brownian motion, the displacement of particles is completely random and normally distributed, with an expected value of 0 for the averaged displacement and a second moment that increases linearly with time. The second moment can be expressed as follows [32]:

The angular ME is a phenomenon in which a speckle pattern shifts according to the tilting of the incident laser beam, as the pattern remembers where the laser originates from. Using the PASIN method, we observed the shift in the speckle pattern by tilting the incident beam and adding a phase ramp to the incident photons. We simulated speckle patterns at different incidence angles. In every speckle pattern, 1.54 × 10¹³ photons were launched, and approximately 1.46 × 10¹¹ photons were collected. The relationship between the shift distance of the speckle pattern and the tilt angle of the incident beam was examined. By recording the number of scattering events that each photon experiences, we separated the scattered photons, compared their ME ranges, and explained the counterintuitive narrower ME range of singly scattered photons relative to that of doubly scattered photons in our case. The results indicate that the optical ME is not yet fully understood.

Theoretically, the average speckle size on the observation screen is given as follows [35]:

3. Results

3.1. Convergence and characteristics of calculated speckle patterns

The calculated intensity distributions versus the numbers of incident photons are shown in Fig. 2(a). Initially, the intensity pattern had a shot noise appearance, and as the photon number increased, it gradually resembled speckles, eventually evolving into a stable speckle pattern. The color bar of each pattern represents the energy in units of photon number. The CCCs of the two independently calculated intensity distributions versus the numbers of photons are shown in Fig. 2(b). As the number of incident photons increased, the CCC increased, which is consistent with the observations. However, owing to the small pitch size and severe fluctuation in the shot noise, the rising speed was low. The CCC was smaller than 0.6 even at 1.54 × 10¹³ incident photons. Proper downsampling can suppress shot noise while preserving the features of speckles. The speckle patterns were downsampled via binning (4 × 4 pixels), that is to say, the new pixel pitch was larger than wavelength but smaller than the grain size. After the downsampling, the speckles were smoother, as shown on the left side of Fig. 2(d), and the CCC was increased. It approached 1 with an increase in the number of incident photons, as indicated by the orange curve in Fig. 2(b). Therefore, downsampling to a pitch of $0.8\textrm{ }\mu m{\kern 1pt}$ was applied in the subsequent simulations. Fig. 2(c) shows the probability density of the intensity of the generated speckle pattern, which agreed with the exponential probability density function of the polarized thermal light.

 

Fig. 2. (a) Evolution of the speckle pattern with an increasing number of photons and downsampled speckles. (a) Speckle patterns corresponding to different numbers of incident photons (N_inc = 3 × 10¹¹, 8.4 × 10¹¹, 5.52 × 10¹², and 1.54 × 10¹³). (b) CCCs of two independent speckle patterns versus the incident photon numbers. (c) Probability density of the intensity of the generated speckle pattern. (d) Autocorrelation curves of the different speckle patterns. The height of the curve—excluding the peak of noise—increased with N_inc. (e) Speckle patterns of two different scattering media after downsampling. (f) Autocorrelation curves of the downsampled speckle patterns. Their peaks, to which noise contributed, were significantly smaller than those of the autocorrelations in (d) and were almost negligible. The patterns were different, but the average speckle sizes were similar.

Download Full Size | PDF

From the fitted curves of the autocorrelation coefficient calculated from the normalized covariance of the speckle pattern (Figs. 2(d) and (f)), the spatial coherence length, that is, the average speckle size, was estimated as the full width at half maximum (FWHM). Owing to the significant shot noise in the original speckles, the autocorrelation coefficient exhibited a steep decay at one pixel shift, as shown in Fig. 2(d). The width of the peak was neglected because the shot noise was white noise; i.e., the shot noise had a trivial effect on the average speckle size. The FWHMs of all four curves shown in Fig. 2(d) were ${\kern 1pt} 6.03\textrm{ }\mu m$, although they had different heights. After downsampling, the FWHMs of the two curves shown in Fig. 2(f) were approximately 6.39$\mu m$, indicating that the effect of the downsampling on the grain size was negligible. However, the spatial coherence length calculated using Eq. (14) was 3.25$\mu m$. This discrepancy can be explained by the anisotropy factor.

When the internal particles shifted, that is, when a different scattering medium with the same macroscopic parameters (denoted as scattering medium #2) was generated, the corresponding speckle pattern exhibited a different distribution, as shown on the right side of Fig. 2(e). However, this new speckle pattern had the same statistical properties; for example, the average speckle size was equal to that of the original speckle pattern, as confirmed by Fig. 2(f).

3.2. Speckle patterns of a dynamic scattering medium

As the time interval ${\kern 1pt} \tau {\kern 1pt}$ increased, the second moment of particle displacement increased, and the difference between the speckle pattern and the original pattern became more significant, as shown in Fig. 3(a). The CCCs of the two speckle patterns versus the time interval ${\kern 1pt} \tau {\kern 1pt}$ are shown in Fig. 3(b). Owing to shot noise, the CCC of the two independently calculated speckle patterns of the medium at ${\kern 1pt} \tau = 0{\kern 1pt}$ was 0.94, which was smaller than 1. The CCC decreased exponentially with an increase in ${\kern 1pt} \tau$, which is consistent with the theoretical prediction [34]. The rate at which the CCC decayed reflected the internal thermal motion of the dynamic scattering medium. Any changes in the internal scatterers led to variations in the speckle pattern. Hence, the PASIN method provides a unique quantitative tool for correlating the internal parameters of a scattering medium with the speckle patterns.

 

Fig. 3. Influence of Brownian motion on speckle patterns: (a) speckle patterns of the medium at different times;(b) CCCs of two speckle patterns (one at${\kern 1pt} t{\kern 1pt}$ and the other at 0 s) versus the time interval ${\kern 1pt} \tau = t - 0.$

Download Full Size | PDF

With regard to blood flow, the aforementioned diffusion model can be used for a complex capillary network. Generally, the measured relative change in ${\kern 1pt} {D_b}{\kern 1pt}$ is proportional to the change in the blood velocity in the capillaries [36]. Thus, we can obtain ${\kern 1pt} {D_b}{\kern 1pt}$ from the rate of speckle decay and estimate the blood-flow velocity accordingly. The particle size of spherical colloids, which is inversely proportional to ${\kern 1pt} {D_b}$, can also be estimated via speckle decorrelation [37].

3.3 Simulation of ME

3.3.1. Observation of angular ME

Fig. 4 shows the results for the angular ME. As the incident beam tilted, the speckle pattern shifted accordingly, but its shape was preserved within the angular ME range. The shift distances at ${\kern 1pt} {\theta _{in}} = 0.18^\circ {\kern 1pt} $ and ${\kern 1pt} 0.36^\circ {\kern 1pt} $ were 16.8 and 32.8$\mu m$, which were consistent with the theoretical values of 16.5 and 33$\mu m$, respectively. Compared with the speckle pattern at ${\theta _{in}} = 0^\circ $, the speckle patterns at these two angles had distortions, although they appeared similar. The CCC curves of the two speckle patterns at ${\theta _{in}} = 0^\circ $ and a certain angle are shown in Fig. 4(b). As the tilt angle increased, the peak of the CCC curve decreased because the distortion in the speckles increased. Nevertheless, the FWHMs of all the CCC curves were equal, indicating that the average grain sizes of the different speckle patterns were the same. From the dashed lines in Fig. 4(b), we estimated the angular ME range, that is, the range when the CCC decreased to half of the maximum normal incidence [38]. Using the simulation parameters, the angular ME range was estimated to be ${\kern 1pt} 0.6^\circ $. This value was an order of magnitude larger than the theoretical prediction of $\Delta {\theta _{in}} = \lambda /2\pi {L_z} = 9.7 \times {10^{ - 3}}^\circ $, owing to the large anisotropy factor of $g = 0.98$ and the small OT of 1.57 [39]. As shown in Fig. 4(b), the FWHMs of different CCC curves and the angular ME range were not as sensitive to noise as the speckle pattern.

 

Fig. 4. Speckle patterns at different incidence angles and their CCC curves at different incident photon numbers. (a) Speckle patterns at different incidence angles. The yellow circles indicate similar structures and their shifts. (b) CCC curves of speckle patterns obtained at different incidence angles $({\theta _{in}} = 0^\circ ,0.18^\circ ,0.36^\circ ,$ and $0.55^\circ )$ with the one at ${\theta _{in}} = 0^\circ $, with launched photons of 7.7 × 10¹² (left) and 1.54 × 10¹³ (right). The center of each CCC shifted accordingly, and the shift was proportional to the corresponding tilt angle. The dashed curves indicate the CCCs of the two speckle patterns versus the tilt angle causing the shift $Z{\theta _{in}}$. The FWHMs represent the angular ME range. Although the heights of the CCC curves for less incident photons were smaller, the FWHMs of all the curves were similar, and the estimated angular ME ranges for the two cases were equal.

Download Full Size | PDF

3.3.2. Angular ME ranges of singly and doubly scattered photons

When transmitted through a scattering medium, different photons experience different scattering events. The scattered photons are separable by the number of scattering events in the PASIN method, which allows detailed investigation of the ME ranges of different scattering components. Fig. 5 presents the simulated ME results for both singly and doubly scattered photons. As shown in Fig. 5(a), the speckle pattern of doubly scattered light had finer grains than that of singly scattered light, which was caused by more severe spatial fluctuations on the exit plane, similar to the effect of the anisotropy factor discussed in the next subsection. Surprisingly, doubly scattered light had a wider angular ME range than singly scattered light. To explain this counterintuitive outcome, we analyzed the effects of the geometric relationship and distribution of scattering events of singly and doubly scattered photons.

 

Fig. 5. Angular ME of singly and doubly scattered photons. (a) Speckle patterns of singly scattered photons (first row) and doubly scattered photons (second row) at different incidence angles (${\theta _{in}} = 0^\circ ,0.18^\circ ,$ and $0.36^\circ $). (b) CCC curves of the singly (S) and doubly (D) scattered light at different incidence angles. According to the two dashed curves, the doubly scattered light had a larger ME range than the singly scattered light. (c) Histogram of positions where scattering events occurred for singly scattered light. (d) Geometric schematic of the singly scattered light with the scattering event in the first (blue line, S1) and second (red line, S2) halves of the scattering layer. The pivot of tilting was on the input surface of the layer. (e) CCC curves of singly scattered light S1 and S2. The peak of the CCC curve of S2 was slightly shifted to the left of the peak of S1, indicating a smaller spatial shift of the speckle pattern with the same tilt angle. (f) Position histograms of the first and second scattering events for the doubly scattered light. The weight of the first scattering event decreased monotonically with an increase in z, and the weight of the second scattering event exhibited the opposite trend.

Download Full Size | PDF

First, a singly scattered photon continued to propagate when it did not exit the medium immediately after a scattering event. There was a probability that it was scattered a second time to become a doubly scattered photon. The field of the residual singly scattered photons contributed to the final wavefront of the singly scattered photons. Under normal incidence, the wavefront of the residual singly scattered photons had a specific distribution. When the incident beam tilted, the field of each singly scattered photons also tilted. However, the scatterers were fixed, and there was a relative shift between the scatterers and fields of singly scattered photons. Consequently, the wavefront of the residual singly scattered photons was altered. Thus, double scattering also influenced the wavefront of the singly scattered light. The effects of single and double scattering on the optical field were coupled. This suggests that the number of scattering events, which is widely adopted in conventional ray-tracing MC methods, may not be an appropriate parameter to label a scattered photon as a wave package.

Second, the position of the scattering event influenced the angular ME range. As indicated by the uniform histogram shown in Fig. 5(c), the single scattering events were equally distributed along the z-direction. We separated the singly scattered photons generated in the first (0 < z1 < 250$\mu m$) and second (250 < z2 < 500$\mu m$) halves of the medium into S1 and S2, respectively. As the incident beam tilted to an angle ${\kern 1pt} {\theta _{in}}{\kern 1pt} $, the beam on the second half of the layer not only tilted to an angle ${\kern 1pt} {\theta _{in}}{\kern 1pt} $ but also shifted to ${\kern 1pt} {L_z}{\theta _{in}}/2$, as shown in Fig. 5(d). When ${\kern 1pt} {\theta _{in}}{\kern 1pt} $ was $0.55^\circ $, the spatial shift in the second half was 2.3$\mu m$. However, it appeared that the shift only degraded the similarity between the speckle pattern and the original pattern at normal incidence. Specifically, the scatterer density of the second half was low, as observed from the incident beam. When there was a small shift of the normally incident beam, the distribution of the scattered field, that is, the speckle pattern formed by scattered photons, on the receiving plane hardly changed (no shift either). Compared with the speckle patterns at ${\kern 1pt} {\theta _{in}} = 0^\circ $, the speckle patterns of S1 and S2 shifted to ${\kern 1pt} (Z + {L_z}/2){\theta _{in}}{\kern 1pt} $ and ${\kern 1pt} Z{\theta _{in}}$, respectively, which agreed well with the offsets shown in Fig. 5(e). The probability of the singly scattered wavefront in S1, which was influenced by another scattering event, was higher than that in S2. Hence, the CCC of S1decayed faster than that of S2, which resulted in a smaller ME range for S1. Owing to the small thickness used in the simulation, the distortion of the speckle pattern caused by the position shift in the second half was negligible. As the thickness or OT increased, the speckle pattern became more sensitive to the shift, and the decay rates of S1 and S2 differed. These results indicated that the influence of the scattering position was related to the position shift of the incident beam, suggesting that the angular ME and shift ME were mixed. This confirmed that the different scattering components were coupled.

Third, for doubly scattered photons, the first and second scattering events had different position distributions, as shown in Fig. 5(f). Most of the first scattering events occurred in the first half of the layer, whereas the second scattering events were primarily concentrated in the second half of the layer. The different distributions of the scattering events also affected the angular ME ranges of the various scattering components.

The wider angular ME range of the doubly scattered light was an integrative result of these factors. The results may be different for samples with different parameters; however, there were still many unknown properties of the ME effect. For example, it is widely accepted that by selecting fewer scattered photons, the ME range can be expanded [40–42]. However, the weight of ballistic photons increased with an increase in the weight of less scattered photons. Whether the wider ME range is related to the less scattered photons remains unclear. To gain a comprehensive understanding, we investigated the field propagation inside the scattering media.

3.4. Influence of the anisotropy factor g on speckles

The results of the PASIN simulation indicated that the anisotropy factor affected the speckle size, as shown in Fig. 6. For media with different anisotropy factors, the corresponding speckle patterns had similar distributions but different contrasts and grain sizes. The reduced contrast at a smaller anisotropy factor was caused by the reduced signal-to-noise ratio on the detection screen. For a smaller anisotropy factor, the average scattering angle of photons was larger, and the number of received photons on the detection plane was smaller; thus, the noise was more significant. The speckle autocorrelation curves—excluding the peaks—had different heights, as shown in Fig. 6(b), which had a trivial influence on their FWHMs, similar to the results shown in Fig. 2(c), except that the speckle patterns were downsampled. The FWHMs at different anisotropic factors of $g =$0.94, 0.95, 0.96, 0.97, and 0.98 were 4.03, 4.26, 4.67, 5.30, and 6.03$\mu m$, respectively. The speckle size decreased with a reduction in the anisotropy factor. This is because a smaller $g$ resulted in a more randomly distributed angle and phase of emitted photons on the existing surface, reducing the spatial coherence length of the speckle patterns. This is supported by Fig. 6(c), which shows that the path-length distribution (and thus the phase distribution) become more random as $g$ decreased. Because the phase of the emitted light was more randomly distributed, it was closer to the precondition of incoherence given by Eq. (14); the spatial coherence length was closer to the calculated value of 3.25$\mu m$.

 

Fig. 6. Speckle patterns and their properties for different anisotropy factors. (a) Speckle patterns for anisotropy factors of $g$= 0.94, 0.95, 0.96, and 0.97. (b) Autocorrelation curves of the speckle patterns at different anisotropy factors. (c) Path-length (after subtracting the straight-line distance between the illuminating source and the detection surface) distributions for different anisotropy factors. The distribution was more random for larger $g$ values.

Download Full Size | PDF

4. Discussions and conclusions

In the PASIN method, the cohort of sampled paths is a subset of the cohort of paths transmitted by photons in a medium. Shot noise is unavoidable but can be suppressed by sampling more paths at the cost of time. However, shot noise does not affect the statistical characteristics of speckles, owing to its white noise feature. If weaker noise is preferred, smoothing filters [43,44] and algorithms (such as deep learning [45]) can be used to reduce noise.

In contrast to FDTD methods, the PASIN method can achieve a tradeoff between computing time and computing memory. Theoretically, there are no limits on the volume of the medium. Large scattering media with high orders of magnitude can be easily simulated and calculated using the PASIN method. From the ME ranges of singly and doubly scattered photons, we found that sophisticated field propagation in scattering media is influenced by many factors, and some factors are coupled. There are unresolved issues related to phenomena that are already thought to be understood. However, the PASIN method can be used for obtaining an in-depth understanding of the field propagation in scattering media.

To give the reader a comprehensive understanding of the PASIN method, major properties of the MC, PF-MC, EMC, and PASIN methods are presented in Table 1. In the conventional MC and PF-MC methods, photons are considered as particles, and the intensity profiles are obtained on the receiving plane. Increasing the number of photons renders the intensity profile smoother, without changing its shape. In the EMC and PASIN methods, photons are treated as waves. However, in the EMC method, scatterers are not fixed, and the fields are coherently superposed; thus, the pattern is not a speckle pattern corresponding to any specific distribution of scatterers. As the number of incident photons increases, it is difficult to predict what the pattern converges to. Because there is no difference between any two points inside the medium, the average of infinite realizations of different distributions of scatterers is a transparent medium. A reasonable assumption is that the coherent superposition of the EMC converges to a laser spot on the detection plane.

Table 1. Properties of the MC, PF-MC, EMC, and PASIN methods

View Table

In summary, we propose the PASIN method to simulate and calculate coherent light propagation in scattering media, which samples possible paths according to their corresponding probabilities. In contrast to traditional MC methods and EMC simulations, this method can coherently superpose all the collected scattered photons to produce a speckle pattern of a scattered medium. The speckle properties and speckle-related ME were re-examined using this method. The PASIN method is a powerful tool for speckle-based research, such as wavefront shaping, speckle autocorrelation imaging, and speckle temporal correlation. It can be used to scrutinize field propagation in volume scattering media under various types of illumination, geometries, and polarization conditions, which can provide insight into the propagation properties and utilization of scattered light.