Angular-spectrum modeling of focusing light inside scattering media by optical phase conjugation

Focusing light inside scattering media by optical phase conjugation has been intensively investigated due to its potential applications, such as in deep tissue imaging. However, no existing physical models explain the impact of the various factors on the focusing performance inside a dynamic scattering medium. Here, we establish an angular- spectrum model to trace the field propagation during the entire optical phase conjugation process in the presence of scattering media. By incorporating fast decorrelation components, the model enables us to investigate the com- petition between the guide star and fast tissue motions for photon tagging. Other factors affecting the focusing performance are also analyzed via the model. As a proof of concept, we experimentally verify our model in the case of focusing light through dynamic scattering media. This angular-spectrum model allows analysis of a series of scattering events in highly scattering media and benefits related applications.

The PBR of the DOPC focus is the most important factor that determines the performance of its applications. The transmission-matrix theory is commonly utilized to analyze the PBR of the DOPC focus in a scattering medium [9,30,31]. Nevertheless, the transmission-matrix theory is unable to account for the specific process of light field propagation inside the scattering media. The influence of some key factors, such as the signal-to-noise ratio (SNR) of the measurement, size of the guide star, and curvature of the spatial light modulator, on the PBR has not been well investigated. More significantly, both the guide star and the fast decorrelation components such as blood flow [32] could tag photons in the process of OPC. This phenomenon inevitably leads to the competition between the fast decorrelation components and the guide star in terms of light focusing. Thus, it is still a challenge to analyze the impact of all these factors on the PBR of the focus inside dynamic scattering media, which poses a barrier for improving the performance of DOPC in dynamic applications.
Here, we establish an angular-spectrum model to trace both the light field propagation in scattering media and the entire OPC process. As an approximation, we divide the volume of a scattering medium into multiple layers, and "squeeze" each layer into an infinitesimally thin plane with non-scattering free space between neighboring planes [33]. We develop a method to distribute the refractive index in each plane, so that a scattering medium with a desired optical property can be simulated. Then, the angular-spectrum method is adopted to propagate an optical field in the free space from one plane to the next plane in the frequency domain [34][35][36]. Using this model, we study the influence of the SNRs of the holograms, size of the guide star, curvature of the SLM, and presence of fast decorrelation components on the PBR of the DOPC focus. To demonstrate the feasibility of this model, a single-shot focusing-through DOPC system [9,37] is built and the PBRs of the DOPC focus, as a function of the system latency, are measured for a dynamic scattering medium. The experimental results agree well with the simulation results. Our findings demonstrate that this model works well for highly scattering media, and might benefit other analyses of scattering events. Figure 1 illustrates the angular-spectrum model established to trace the light field propagation for a DOPC system. To focus light inside scattering media by DOPC, two steps are involved:(1) guide-star modulation and hologram recording [ Fig. 1(a)], and (2) conjugated-light playback and focusing [ Fig. 1(b)]. Figures 1(c) and 1(d) show the specific processes in the scattering medium corresponding to the two steps in our model.

A. Angular-Spectrum Model for DOPC
Here, we adopt the angular-spectrum method (see Section 1 and Fig. S1 in Supplement 1) to study the entire process of focusing light in a scattering medium by OPC, including the light field propagation in free space and scattering medium. As shown in Fig. 1(a), we assume that the guide star is located inside a scattering medium at distances of l 1 and l 2 from the front and rear medium surfaces, respectively. The guide star has two different states [14]: high absorption and no absorption. When in the high-absorption state, we assume for simplicity that all light passing through the guide star is absorbed. Assume that the light beam impinging upon the scattering medium has the complex amplitude of E in (x, y), where x and y are the transverse Cartesian co-ordinates. After propagating through a distance of l 1 in the scattering medium, the optical field becomes E a (x, y) on the guide-star plane. The optical field becomes E b (x, y) after being modulated by the guide star, and E c (x, y) after propagating through a distance of l 2 . Then, the field propagates in free space with a distance of d 1 , passes through a collecting lens, and finally arrives at the surface of a camera with a distance of d 2 , resulting in E det (x, y).
In this process, free-space propagation can be calculated by the angular-spectrum method, and the lens is modeled by a transmission function exp[−iπ(x 2 + y 2 )/λf], where f is its focal length and λ is the wavelength of light. For the scattering medium, we divide it into M layers and compress each layer into a plane with non-scattering free space between neighboring layers, as shown in Fig. 1(c), and calculate the field propagation layer by layer. We assume that each layer can be divided into many speckles with an average speckle size of D speckle . In each speckle, the refractive index is modeled as a single constant. The refractive indices for different speckles are Gaussian randomly distributed with a mean value of n mean and a standard deviation of n dev . n dev is controlled to yield the desired transport mean free path l′ (see Section 2B below), and the absorption of the scattering medium is neglected. The interval between two adjacent layers is Δl = l∕M, where l is the thickness of the scattering medium. Thus, the field propagation from E m (x, y) on layer m to E m+1 (x, y) on layer m+1 can be modeled in two steps. First, a random phase map P m (x, y) = 2π · n m (x, y) · Δl∕λ is directly imposed on E a (x, y), where n m (x, y) is the refractive index at point (x, y) of layer m. Then, it propagates in free space with a distance of Δl, which is calculated by the angular-spectrum method. As for layers g and g + 1 which sandwich the guide star, the transmission function of the guide star should be multiplied in addition to the two preceding steps.
To measure the field E det (x, y), a collimated reference beam E ref (x, y) is introduced for interferometry and a hologram is recorded. By switching the state of the guide star, two holograms I holo1 (x, y) and I holo2 (x, y) can be obtained to produce a binary phase map [see Eq. (S5) in Supplement 1] [38,39]. In the DOPC system, the phase map is used to modulate the playback beam E mod (x, y) by an SLM, which is located at the conjugate position of the camera, as shown in Fig. 1(b). For the playback process, the same angular-spectrum procedure just described is adopted to trace the optical field propagation. Once the light propagates layer by layer back to the guide-star plane, an optical focus with the field of E focus (x, y) is formed, as shown in Fig. 1(d).

B. Model Implementation
To implement the model, all the parameters for the simulation were initialized. The incident beam was a collimated laser beam with a wavelength of 532 nm and diameter of 1 mm. The thickness of the scattering medium was 8 mm, and the interval Δl between the divided layers of the scattering medium was set to 20 μm, which was much smaller than the transport mean free path of the scattering medium. The mean value of the refractive index of the scattering medium was n mean =1.4, approximately the value of biological tissue [40]. A guide star was set at the middle of the scattering medium with a diameter of 50 μm. The distances d 1 and d 2 were 30 mm and 80 mm, respectively. The focal length of the collecting lens was 40 mm. Two-dimensional fast Fourier transformation [41] was used to realize the decomposition and superposition of the optical field. The optical fields on each plane were discretized into 2000 × 2000 points with an interval of 5 μm. Accordingly, the speckle size D speckle was set to the same value as that of the interval. Both the camera and the SLM were represented by 1000 × 1000 pixels, with a pixel size of 5 × 5 μm 2 .
In our model, the transport mean free path l′ of the scattering medium is determined by n dev . It is necessary to first express the relationship between n dev and l′ In photon diffusion theory, the transport mean free path is the mean distance after which a photon's direction becomes random compared with its original incident direction [42]. When a collimated beam with a beam size much larger than the optical wavelength illuminates normally on a scattering medium, the ballistic light corresponds to the direct-current (DC) component of the angular spectrum (k space), while the scattered light (which has a random propagation direction) corresponds to the whole region of the available k space. Thus, we can calculate the transport mean free path in k space by identifying the critical depth at which I DC ∕I AC falls to unity, where I DC is the intensity of the DC component, by reading the value at the original point in k space; I AC is the average intensity of the alternating-current (AC) components. Assuming that n dev = 1.5 × 10 −3 , the k-space intensity distribution I(k x , k y ) of the scattered light at different depths were simulated. To eliminate the influence of speckles, I(k x , k y )was averaged 100 times over different refractive index distributions. As shown in Fig. 2(a), the left three panels of the intensity distribution show the single simulation result of I(k x , k y )at depths of 0.5, 0.7, and 1.2 mm, while the right three panels show the corresponding results after averaging 100 times. Figure 2(b) shows I DC ∕I AC as a function of depth. When n dev =1.5 × 10 −3 , the transport mean free path l′ was found to be 1.16 mm at I DC ∕I AC =1. Similarly, we found the transport mean free path l′ to be 0.68 mm when n = 2.0 × 10 −3 . Accordingly, the relationship between n dev and l′ was calculated and shown in Fig.  2(c) when Δl = 20 μm. Based on these data, we set n dev as 1.6 × 10 −3 in our simulation to mimic biological tissues with l′ = 1 mm [42,43]. Note that the relationship between n dev and l′ should be recalculated when Δl is changed. However, the simulation outcome would remain the same. Figure 3 shows the simulation of light propagation during the entire DOPC process via the angular-spectrum method. In the recoding process [ Fig. 3(a)], the optical field associated with two different guide-star states, high absorption and no absorption, were simulated. Figures 3(a1)-3(a3) show the intensities of the scattered light before and after being modulated by the guide star in the high-absorption state, and on the rear plane of the scattering medium. Figure 3(a4) represents the hologram captured by the camera.

A. Light Propagation Simulation
By switching the state of the guide star between high absorption and no absorption, two holograms were obtained to calculate the binary phase map shown in Fig. 3(b1). In the playback process, the binary phase map was displayed on the SLM to modulate the reference beam. After being reflected by the SLM, the reference beam propagated in the opposite direction. Figure 3(b2) shows the intensity distribution of the playback beam incident on the scattering medium. After propagating this beam layer by layer backward inside the scattering medium, we finally obtained a bright focus [ Fig. 3(b3)] at the position of the guide star. We calculated the averaged intensity within the bright focus as the peak value for the PBR calculation. The boundary of the focus was determined by the area of the guide star. The background value was computed by displaying a uniform phase map on the SLM and tracing the optical field along the same path, resulting in a fully scattered optical field on the guide-star plane. Averaging the intensity near the guide star yielded the background value. Dividing the peak value by the background value returned a PBR of 253 in this simulation. The bright focus observed at the position of the guide star proves that the established model works well for DOPC.

B. Impact of the Camera's SNR, Guide-Star Size, and SLM Curvature
Noise in the recording process from the camera inevitably decreases the PBR of the focus inside scattering media, which to our knowledge has never been quantitatively studied before. In our model, the holograms corresponding to the two states of the guide star were obtained. Thus, it is straightforward to add noise to the holograms and analyze the influence of the SNR on the PBR of the focus [ Fig. 3(c)]. Under the conditions set in this particular simulation, the SNR of the hologram should be at least 11.7 to have a PBR greater than 2. For an SNR of 200, which is considered high for most cameras, the PBR of the focus becomes 176, which is much lower than that without noise (PBR = 253). Thus, the impact of the SNR on the PBR of the focus is substantial, which cannot be ignored in the experiments.
As predicted by the transmission-matrix theory, the PBR is inversely proportional to the square of the guide star size [9,30,31]. However, experimentally, the wavefront measurement of the tagged photons would be extremely sensitive to noise when the guide star is small. Thus, the inverse-proportional relationship breaks down when the SNR of the camera is insufficient. As shown in Fig. 3(d), we simulated the PBRs as a function of the guide-star size with different SNRs of 100, 200, and 300. When the guide star is small, the influence of the SNR on the PBR of the focus is significant. In contrast, the PBR becomes insensitive to the SNR when the guide star is sufficiently large. Thus, our model enables us to precisely describe the relationship between the guide-star size and the PBR in the presence of noise.
Due to manufacturing and alignment errors, the curvature and tilt of the SLM are two main factors adding aberrations to the wavefront of the playback beam, which deteriorate the PBR of the focus inside scattering media. The introduced aberrations can be expressed by Zernike modes. Here, five serious aberrations are considered, including the second and third Zernike modes that represent tilts along two orthogonal directions, and the fourth, fifth, and sixth Zernike modes that respectively represent oblique astigmatism, defocus, and vertical astigmatism. In our model, the aberrations are added to the playback wavefront and their influences are analyzed. Figure 3(e) shows how they affect the PBR of the focus inside scattering media. As we can see, the tilt modes have the most serious effects, while the defocus mode has the least serious impact. Fortunately, these aberrations can be measured in advance and then compensated for by adding a conjugated phase map to the SLM.

C. Competition Between the Guide Star and Fast Tissue Motion for Photon Tagging
For in vivo applications, fast motion in tissues such as blood flow would tag photons simultaneously with the guide star in the process of OPC. Thus, the two tagging mechanisms compete for time-reversal focusing, which can greatly decrease the PBR of the focus at the guide star. Moreover, with the increase of the system latency, the ratio of the number of photons tagged by the fast motion components to that tagged by the guide star would increase accordingly. However, the controlled mode number determined by the pixel count of the SLM is fixed. Thus, the PBR of the focus decreases with the increase of the system latency. To study the competition and the relationship between the system latency and PBR, we embedded fast movement between layers p and p + 1 in the scattering medium, as shown in Fig. 4(a). The field transmittance through the movement at each pixel was randomly changed with an additional phase within [−απ, απ], and an additional amplitude within [−β, β] for each step. By changing the values of α and β, the speckle decorrelation time due to the motion can be controlled. Figure 4(b) shows the speckle correlation coefficient as a function of time t, when the size of the fast movement was 0.5 × 10.0 mm 2 , p = 150, α = 0.16, and β = 0.03, and the time interval of each step was 1 ms. Accordingly, the decorrelation time of the dynamic scattering medium was found to be τ c = 31.7 ms by fitting the speckle correlation coefficient α cor with the function α cor = A exp(−2t∕τ c ) + B, where A and B are constants. The R-squared value of the fitting is 0.998. Assume that the system latency t lat was divided equally into two parts to capture the two holograms without considering the time for calculating and displaying the binary phase map. Thus, we simulated hologram 1 at t = 0 ms with the guide star in the high-absorption state, and hologram 2 at t = t lat /2 with the guide star in the no-absorption state. Finally, we obtained the DOPC focus at t = t lat right after the two holograms were captured. Figure 4(c) shows the DOPC foci with t lat of 0, 2, 6, and 16 ms. Figure 4(d) shows the PBR of the DOPC focus as a function of system latency time t lat . The PBR of the focus decreases exponentially as the system latency increases. When t lat = 1.7 ms, i.e., 5% of the decorrelation time, the PBR decreases to half of the original value. When the system latency time is longer than 16.1 ms, i.e., 50% of the decorrelation time, the PBR reduces to only 10%. Thus, to obtain a focus with a relatively high PBR inside a dynamic scattering medium, the system latency should be much shorter than the decorrelation time.

D. Experimental Verification
To verify the present angular-spectrum model and its ability to analyze the influence of speckle decorrelation on the focusing contrast, we established a single-shot focusing-through experimental system [9] (see Section 2 and Fig. S2 in Supplement 1 for details). The dynamic scattering medium consisted of two pieces of 1 mm thick 1% intralipid gelatin phantoms and a tube filled with 1% intralipid solution. The cross-section of the tube was a 1 × 2 mm 2 rectangle, and the tube was sandwiched between the phantoms with a gap of 2 mm. A syringe pump was used to drive the intralipid solution in the tube. The decorrelation time can be controlled by tuning the driving speed of the syringe pump (see Section 3 and Fig. S3 in Supplement 1). To perform the DOPC focusing, we employed an SLM with a pixel count of 1920 × 1080 and a pixel size of 8 × 8 μm 2 . We used a camera lens to relay the SLM plane to the camera plane and match their pixels with a 1:1 ratio.
Accordingly, we applied the model to the focusing-through DOPC system, as illustrated in Figs. 5(a) and 5(b). By adjusting the parameters α and β, and the time interval for each step, the decorrelation time was set to 78 ms. In the experiments, the decorrelation time of the dynamic scattering medium was also adjusted to 78 ms by changing the driving speed of the syringe pump. The speckle correlation coefficient as a function of time t for the simulation and experiment are shown in Fig. 5(c). Figure 5(d) compares the PBR of the DOPC focus as a function of system latency between the simulation and experiment. Even though the experimental PBR values were less than the simulated ones, the trends of the two curves agree well, which verifies that our model works well for dynamic scattering media.
The reasons that the experimental PBRs were lower than the simulation ones likely include (1) imperfect alignment in the setup, (2) residual aberration in the playback beam, and (3) greater noise.

CONCLUSION
In conclusion, we have established an angular-spectrum model to trace the light field propagation during the entire OPC process in the presence of scattering media. We have also developed a method to simulate a scattering medium with a desired optical property. Utilizing this model, the influence of the key factors on the PBR of DOPC focus, including the SNR of the camera, size of the guide star, and curvature of the SLM, have been studied. Most importantly, we have applied this model to investigate the effect of the system latency on the PBR of the focus inside dynamic scattering media. To demonstrate the feasibility of this model, a single-shot focusing-through DOPC experimental system was built, and the relationship between the system latency and the PBR of the DOPC focus was measured. The experimental results showed good agreement with the results predicted by the model, which verifies that the present model works well for dynamic scattering media.
Our angular-spectrum model can be further optimized. In the current implementation, the sampling interval of the optical field was 5 μm, larger than the speckle size (~λ∕2) inside a thick scattering medium. This sub-Nyquist sampling of the light field may lead to deviations from the true PBR values [44]. We can improve it by using finer sampling in the simulation at the expense of computing resources. Additionally, in our current model, only the forward propagation of light inside the scattering medium was considered and the absorption of the medium was neglected. The model can be further improved by taking the backward propagation and absorption into account. Nevertheless, our current implementation is effective for studying the trend of the PBR impacted by different factors.
Apart from being used to study the focusing performance inside scattering media, this angular-spectrum model can also be exploited for analysis of other effects related to light propagation through scattering media, such as memory effects [45][46][47], time delay eigenstates [48], and depolarization. Therefore, our methodology is expected to benefit optical imaging [26], sensing [49], and communication [50] in scattering media.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.

1.
AU: A mismatch has been discovered between the e-mail address in the manuscript and the e-mail address in OSA's system. Please let us know if the e-mail address in the manuscript should be changed.

2.
AU: The funding information for this article has been generated using the information you provided to OSA at the time of article submission. Please check it carefully. If any information needs to be corrected or added, please provide the full name of the funding organization/institution as provided in the CrossRef Open Funder Registry (https://search.crossref.org/funding).       Optica. Author manuscript; available in PMC 2020 February 05.