Path sampling and integration method to calculate speckle patterns

: A stable speckle pattern is generated when a coherent beam illuminates a stationary scattering medium that contains numerous scatterers with fixed positions. To date, there has been no valid method to the best of our knowledge for calculating the speckle pattern of a macro medium with a large number of scatterers. Here, a new method based on possible path sampling with corresponding weights and coherent superposition is presented for the simulation of optical field propagation in a scattering medium and output speckle patterns. In this method, a photon is launched onto a medium with fixed scatterers. It propagates in one direction; upon collision with a scatterer, its direction is updated. The procedure is repeated until it exits the medium. A sampled path is obtained in this manner. By repeatedly launching photons, numerous independent optical paths can be sampled. A speckle pattern, corresponding to the probability density of the photon, is formed by the coherent superposition of sufficiently sampled path lengths ending on a receiving screen. This method can be used in sophisticated studies of the influences of medium parameters, motion of scatterers, sample distortions on speckle

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][24][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][27][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.

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. 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.
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 θ and ϕ. The probability distributions of θ and ϕ can be expressed as where g is the anisotropy factor. Accordingly, the scattered intensity is expressed as where I 0 represents the intensity of the incident beam on the spherical particle.I 0 = E 0 2 , thus, we have the following: Taking the square root yields Only one root was selected for this study. We used p 2 (θ, ϕ), which denotes the phase function of the scattered field.
Here, k is the normalized coefficient. The azimuthal angle ϕ is uniformly distributed over [0 2π). The cumulative distribution function can be expressed as follows: After the derivation of P(θ), a marginal probability density function was obtained In the simulation, cos θ instead of θ was sampled [31]. Hence, the following was obtained: The sampled cos θ for the scattered field is within [-1, 1], thus, the integration of p 2 (cos θ) should equal 1, as follows: Then, the coefficient k was calculated to be Eq. (8) was thus transformed as follows: 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 10 photons per core per hour. Typically, 80 CPUs were allocated to run the simulation.

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 mm was normally incident on a scattering layer. In this and the following simulations, the medium size was L x × L y × L z = 3 × 3 × 0.5 mm 3 , and the cell size wasl x × l y × l z = 10 × 10 × 10 µm 3 . In the medium, the refractive index of the background was n b = 1, and the refractive index of the particles was n p = 1.06. All the particles had the same radius (r p = 1 µm), and the calculated anisotropy coefficient was g = 0.98 [31]. The scattering mean free path is expressed as follows: where ρ represents the particle concentration. In this study, we used ρ = 1 × 10 6 mm −3 , and l s = 0.32 mm, and the corresponding optical thickness (OT) of the scattering layer was L z /l s = 1.57. The recording plane was 5 mm from the scattering medium. The area of interest on the plane had an array size of 2000×2000 and a pixel pitch of 0.2µm. The intensity distributions for different numbers of launched photons were obtained. Only scattered photons were collected, because we were interested in the scattered light. Eventually, 1.54 × 10 13 photons were launched and 1.46 × 10 11 photons were collected in this investigation. The computational time was approximately 10 h. When the internal distribution of the scattering layer was changed, the speckle pattern changed accordingly. We also calculated the speckle pattern of the medium under the same macroscopic parameters but with a different internal distribution.
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]: where D b = k B T/6πηr p is a diffusion coefficient, τ represents the time interval between two instants, k B is Boltzmann's constant, T represents the absolute temperature, and η is the viscosity coefficient. Over time, the particles shifted, and the speckle pattern changed. Assuming that D b = 2.5 µm 2 /s, we obtained speckle patterns at different instants and calculated their temporal cross-correlation coefficients (CCCs). Except for the positions of the particles, all the parameters remained the same. The temporal correlation of speckle patterns can be used to measure the internal changes and physical properties of the dynamic scattering media, such as the estimation of fluid diffusion [33] and blood flow [34] within biological tissues. 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 13 photons were launched, and approximately 1.46 × 10 11 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]: where Φ represents the diameter of the selected aperture on the rear surface. The premise of this equation is that a coherent beam impinges on a screen with a completely random phase distribution and illuminates an observation plane in the Fresnel region. When the coherent beam interacts with a thick volumetric scattering medium, the calculated speckle sizes agree with measurements. However, the speckle grain sizes obtained in our simulation were significantly larger than those predicted using Eq. (14). We assumed that parameters such as the anisotropy factor can affect the speckle size. Hence, the influence of the anisotropy factor on the average speckle size was investigated. We changed this factor to 0.97, 0.96, 0.95, and 0.94 while the other parameters remained the same, including the output aperture diameter of Φ = 1 mm. In every speckle pattern, 7.68 × 10 12 photons were launched, and 5.61 × 10 10 , 4.84 × 10 10 , 4.22 × 10 10 , and 3.71 × 10 10 photons were collected, respectively.

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 13 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 µm 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. 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 6.03 µm, although they had different heights. After downsampling, the FWHMs of the two curves shown in Fig. 2(f) were approximately 6.39µ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µ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).

Speckle patterns of a dynamic scattering medium
As the time interval τ 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 τ are shown in Fig. 3(b). Owing to shot noise, the CCC of the two independently calculated speckle patterns of the medium at τ = 0 was 0.94, which was smaller than 1. The CCC decreased exponentially with an increase in τ, 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. With regard to blood flow, the aforementioned diffusion model can be used for a complex capillary network. Generally, the measured relative change in D b is proportional to the change in the blood velocity in the capillaries [36]. Thus, we can obtain D b from the rate of speckle decay and estimate the blood-flow velocity accordingly. The particle size of spherical colloids, which is inversely proportional to D b , can also be estimated via speckle decorrelation [37]. 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 θ in = 0.18 • and 0.36 • were 16.8 and 32.8µm, which were consistent with the theoretical values of 16.5 and 33µm, respectively. Compared with the speckle pattern at θ in = 0 • , the speckle patterns at these two angles had distortions, although they appeared similar. The CCC curves of the two speckle patterns at θ in = 0 • 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 0.6 • . This value was an order of magnitude larger than the theoretical prediction of ∆θ in = λ/2πL z = 9.7 × 10 −3 • , 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. 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.

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.
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µm) and second (250 < z2 < 500µm) halves of the medium into S1 and S2, respectively. As the incident beam tilted to an angle θ in , the beam on the second half of the layer not only tilted to an angle θ in but also shifted to L z θ in /2, as shown in Fig. 5(d). When θ in was 0.55 • , the spatial shift in the second half was 2.3µ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 θ in = 0 • , the speckle patterns of S1 and S2 shifted to (Z + L z /2)θ in and Zθ 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. 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.
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][41][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.

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µ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µm.

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. 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.
Funding. National Natural Science Foundation of China (81930048).