1. Introduction

In biophotonics, the interaction of light waves with biological specimens is the basis for optical diagnosis, phototherapy, and light stimulation. In in vivo applications, a typical prerequisite condition is that sufficient light energy should be delivered to the target region to initiate the light-specimen interaction. However, multiple light scattering in biological tissues leads to the spatial spread of the propagating waves, resulting in the attenuation of the light energy with an increase in the propagation depth. The total energy delivered to depth z is approximately given by $I(z )\approx {l_t}/z$ if z is sufficiently larger than the transport mean free path ${l_t}$, which is typically on the order of a millimeter or so [1–3]. Accordingly, the working depth of biophotonics technologies has been limited to the superficial layers.

While light propagation through complex biological tissues is considered a random diffusion process, the wave nature of light can be preserved when the sensitivity of the detection system exceeds the ratio between the broadening of the path length by multiple light scattering and the coherence length of the light source. In this context, the control of wave propagation through wavefront shaping of the illuminating beam has emerged as a promising solution in the past decade [4–7]. The main concept of wavefront shaping is that the incident wavefront is adjusted such that constructive interference of multiply scattered waves occurs in the plane of interest. Initial studies on wavefront shaping were aimed at improving light penetration through a scattering layer in the transmission geometry [8–13]. These studies made significant conceptual advances, but they were not suitable for in vivo applications because a detector should be placed on the other side of the scattering layer. To promote their applicability to in vivo applications, researchers introduced the reflection mode of operation, first in the steady state [14–16] and then in the time-gated regime [17–21]. In particular, time-gated detection proved to be effective in focusing light energy on embedded targets as it preferably optimizes the backscattered signals that have interacted with the target relative to the background multiple scattering noise.

A few types of operations have been introduced for optimizing light energy delivery, and each of them has its own pros and cons. The most ideal approach is to measure either the transmission [5,13] or the reflection matrix [14,15,20] and to couple light to the individual eigenchannel [11,20]. While this matrix approach ensures the maximum light energy delivery that the given experimental configuration can support, it requires the scattering medium to be stable during the recording of the matrix, which is often challenging in the case of living specimens. Optimization approaches based on iterative feedback control are advantageous in this respect as they can cope with mechanical instability. One of the most straightforward iterative feedback methods is the point optimization method, where the signal at a single speckle grain is optimized [8]. The optimization ensures the maximization of the light intensity at the target speckle grain, and therefore, the method has been particularly suitable for imaging [22]. Although the operation also enhances the light energy delivery [23], the solution cannot converge to the eigenchannel with the largest eigenvalue. In an ideal waveguide geometry filled with absorption-free scattering medium, the theoretical limit on the enhancement of transmittance is 2/3. Another approach is to directly optimize the total light energy through iterative feedback control [16,24]. Half of the illumination area is randomly selected in each iteration, and the total output intensity is maximized by a three-step phase-shifting method. Theoretically, the optimized solution can converge closer to the open eigenchannels, especially when the maximum eigenvalue is distinctly larger than the other eigenvalues. However, the convergence is rather slow. A more effective optimization method is iterative phase conjugation, which was initially used in acoustics and later in microwaves, optics, and acousto-optics [25–28]. The phase conjugation algorithm effectively purifies eigenchannels with large eigenvalues as it enhances the amplitude of each eigenchannel by powers of its eigenvalue. However, careful alignment is required to match the imaging system that measures the scattered waves and the wavefront shaping system performing the phase conjugation of the measured wavefront. So far, the use of iterative phase conjugation has been limited to steady-state measurements.

Here, we present a time-gated iterative phase conjugation method for focusing light energy on a target embedded in a scattering medium with a faster convergence rate than previous time-gated approaches. We measure the backscattered waves from the sample arriving at a certain flight time and send their phase conjugation back to the sample by using a spatial light modulator. By iterating this process, the illumination wave converges to the eigenchannel with the largest eigenvalue; notably, only a few iterations are required. The optimization time is shortened by hundreds of times compared with that of a random optimization approach. Since the projected wavefront is updated in each step, the process can adaptively cope with vibrations or displacement of the sample. Moreover, the convergence rate is almost unchanged even for an enlargement of the field of view, in contrast to other methods. All these advantages are expected to facilitate the use of wavefront shaping based on time-gated detection in in vivo and in situ optical imaging and stimulation.

2. Experiment

2.1 Experimental setup

A schematic of the experimental setup for the time-gated iterative phase conjugation is shown in Fig. 1. We constructed a reflection-mode interference microscope using an ultrafast Ti:sapphire laser (center wavelength: 780 nm, spectral bandwidth: 30 nm) as the light source. The wavefront of the sample beam was shaped by a spatial light modulator (SLM, X10468-02, Hamamatsu). The wave reflected from the sample was captured by an objective lens, OL_R (40×, 0.8 NA), and delivered to the camera (CAM_R, sCMOS, pco.edge 4.2). We adjusted the optical path length of the reference arm by using a scanning mirror to choose the gating time of the sample beam. The temporal resolution of the time gating is equivalent to the coherence length of 10 µm in terms of the reference mirror displacement, which was close to the theoretical expectation set by the bandwidth of the light source. A tilted reference beam generated by a diffraction grating (Ronchi 72/mm, Edmund Optics) was used to measure the time-gated and complex-field images through off-axis digital holography [29]. To monitor the intensity of light delivered to the target object, we measured the transmission image using a CCD camera (CAM_T, LM135M, Lumenera) and an objective lens, OL_T (20×, 0.4 NA). To implement optical phase conjugation, we should perform pixel-to-pixel matching of the image plane (CAM_R) and illumination (SLM). Translational and rotational misalignments were finely corrected through hardware alignment and numerical compensation [30].

 

Fig. 1. Experimental schematic of the time-gated iterative phase conjugation. Light source: Ti:sapphire laser, SLM: spatial light modulator, OL_R and OL_T: objective lenses for reflection and transmission measurements, respectively, CAM_R and CAM_T: cameras for reflection and transmission imaging, respectively, SM: scanning mirror for adjusting the path length of the reference arm, DG: diffraction grating. The detailed layout of the sample stage is shown in the red box. Sample wave, reference wave, reflected wave, and transmitted wave are colored in red, golden, yellow and green for visibility although their optical wavelengths are the same.

Download Full Size | PDF

2.2 Sample preparation

The detailed layout of the sample stage is shown in the red box of Fig. 1. We used a mouse skull as a scattering medium and measured the light energy delivery enhancement for a target underneath the skull. The skull was excised from a nine-week-old mouse and was then embedded in agar gel and placed on top of a target object. The thickness of the mouse skull was about 200 µm. We used a silver disk embedded in a 1-µm-thick polymethylmethacrylate (PMMA) layer as the target object. The disk diameter was 10 µm, and its thickness was 30 nm; the disk was sufficiently thin to have a transmittance of about 70%.

2.3 Optimization process

The detailed procedure for the time-gated iterative phase conjugation is shown in Fig. 2. The laser input was shaped by the SLM to generate an incident pulse with electric field ${\boldsymbol E}_j^{in}$ for the j^th iteration. Initially, we prepared a random phase pattern ${\boldsymbol E}_{j = 1}^{in}$ as the incident wave. We then measured the backscattered wave by making it interfere with a pulsed reference beam. We obtained the time-gated reflection field ${\boldsymbol E}_j^{out}({t = \tau } )$ at the camera CAM_R, where $\tau = {z_t}/c^{\prime}$. Here, ${z_t}$ is the depth of the target object and $c^{\prime}$ is the average speed of light in the scattering medium. The SLM generated the wavefront of the conjugation of the gated backscattered wave ${[{{\boldsymbol E}_j^{out}({t = \tau } )} ]^\ast }$ for the next input ${\boldsymbol E}_{j + 1}^{in}$. We continued this process until ${|{{\boldsymbol E}_j^{out}({t = \tau } )} |^2}$ was converged.

 

Fig. 2. Workflow of the time-gated iterative phase conjugation process. The green arrows indicate the flow of the process. ${\boldsymbol E}_j^{in}$ and ${\boldsymbol E}_N^{in}$: incident pulses shaped by the SLM after the j^th and N^th iterations, respectively, ${\boldsymbol E}_j^{out}(t )$: waves reflected from the sample as a function of the flight time t, ${\boldsymbol E}_j^{out}(\tau )$: backscattered wave detected at the gating time $t = \tau \; $.

Download Full Size | PDF

The experimental procedure shown in Fig. 2 can be described mathematically by successive operations on the time-gated reflection matrix $R(\tau )$, as shown in Fig. 3. The input-output response of the system can be described by $R(\tau )$, whose column index corresponds to the SLM pixel $({{x_{in}},{y_{in}}} )$ and row index to the camera pixel $({{x_{out}},{y_{out}}} )$. In the iterative phase conjugation process, $R(\tau )$ is not measured directly here, but it is used to describe the phase conjugation process. $R(\tau )$ can be decomposed as $R(\tau )= U{\Sigma }{V^\dagger }$, where ${\Sigma }$ is a diagonal matrix whose diagonal element ${\sigma _m}$ (called singular value) is a non-negative real number. Typically, ${\sigma _m}$ is sorted in the descending order with respect to the eigenchannel index m, namely, ${\sigma _m} \ge {\sigma _{m + 1}}$. V and U are unitary matrices whose column vectors ${{\boldsymbol v}_m}$ and ${{\boldsymbol u}_m}$ correspond to the eigenchannel in the input and output planes, respectively. Since the initial incident wave is a random phase pattern, ${\boldsymbol E}_{j = 1}^{in}$ can be written as a linear superposition of the input eigenchannels, ${\boldsymbol E}_{j = 1}^{in} = \mathop \sum \limits_{m = 1}^{{N_{ch}}} {c_m}{{\boldsymbol v}_m}$, where ${|{{c_m}} |^2} $ is $1/{N_{ch}}$ regardless of m. Here, ${N_{ch}}$ is the number of eigenchannels. The reflected wave is then given by ${\boldsymbol E}_{j = 1}^{out}(\tau )= R(\tau ){\boldsymbol E}_{j = 1}^{in} = \mathop \sum \limits_{m = 1}^{{N_{ch}}} {c_m}{\sigma _m}{{\boldsymbol u}_m}$. In the phase conjugation process, ${\boldsymbol E}_{j = 1}^{out}(\tau )$ is sent back to the medium as the next input. This can be described as ${\boldsymbol E}_{j = 2}^{in}(\tau )= {({{R^\textrm{T}}{\boldsymbol E}_{j = 1}^{out}{{(\tau )}^\ast }} )^\ast } = R{(\tau )^\dagger }{\boldsymbol E}_{j = 1}^{out}(\tau )= \mathop \sum \limits_{m = 1}^{{N_{ch}}} {c_m}\sigma _m^2{{\boldsymbol v}_m}$., and it is then used as the input in the next iteration, ${\boldsymbol E}_{j + 1}^{in}$. After the number of iterations is increased to j = N, the input is given by ${\boldsymbol E}_N^{in}(\tau )= \mathop \sum \limits_{m = 1}^{{N_{ch}}} {c_m}\sigma _m^{2({N - 1} )}{{\boldsymbol v}_m}$. When N is increased until ${\left( {\frac{{{\sigma_1}}}{{{\sigma_2}}}} \right)^N} \gg 1$, i.e., $\sigma _1^N \gg \sigma _2^N \gg \cdots \gg \sigma _{{N_{ch}}}^N$, ${\boldsymbol E}_j^{in}(\tau )$ converges to ${{\boldsymbol v}_1}$, the eigenchannel with the largest singular value. Therefore, the convergence is mainly dictated by the ratio $\frac{{{\sigma _1}}}{{{\sigma _2}}}$. Considering constant input energy, ${\boldsymbol E}_{j = N}^{in}$ is given by ${\boldsymbol E}_{j = N}^{in} = {A_N}\mathop \sum \limits_{m = 1}^{{N_{ch}}} {c_m}\sigma _m^{2({N - 1} )}{{\boldsymbol v}_m}$, where ${A_N}$ is a normalizing constant, and can be written as ${|{{A_N}} |^2} \propto \; \frac{{{N_{ch}}}}{{\mathop \sum \nolimits_{m = 1}^{{N_{ch}}} \sigma _m^{4({N - 1} )}}}$. Then, the output energy after N iterations becomes ${I_{j = N}} = \frac{{\mathop \sum \nolimits_{m = 1}^{{N_{ch}}} \sigma _m^{4({N - 1} )+ 2}\; }}{{\mathop \sum \nolimits_{l = 1}^{{N_{ch}}} \sigma _l^{4({N - 1} )}}}$, and the enhancement of iterative phase conjugation process can be described as ${\eta _{re}} = {I_{j = N}}/{I_{j = 1}}$.

 

Fig. 3. Mathematical representation of iterative phase conjugation. ${\boldsymbol E}_j^{in}$ is multiplied with $R(\tau )$, and its output field ${\boldsymbol E}_j^{out}(\tau )$ is multiplied with $R{(\tau )^\dagger }$ to obtain ${\boldsymbol E}_{j + 1}^{in}$, which is used as an input for the next iteration.

Download Full Size | PDF

3. Experimental results

3.1 Experimental results

The experimental demonstration of the time-gated iterative phase conjugation is shown in Fig. 4. The temporal response of the time-gated reflection intensity is shown in Fig. 4(a) before (blue curve) and after (red curve) the iterations. The flight time of the reflected wave from the upper surface of the skull was set to zero ($t = 0$). The large peak at the flight time of $\tau = 3.68\; ps$ corresponds to the reflection from the target object. The temporal width of the peak was considerably broader than the temporal resolution of the time-gated detection (∼60 fs) owing to multiple scattering by the skull. Signals in the flight time range between 0 and $\tau $ result from backscattering from the skull. We observed that the signal from the target object increased after the optimization. We defined the reflection enhancement factor ${\eta _{re}}$ as the ratio of the optimized intensity to the initial intensity at $t = \tau $. Figure 4(b) shows the reflection enhancement factor of the time-gated reflection intensity as a function of the number of iterations. ${\eta _{re}}$ increased rapidly and saturated in about 10 iterations. The reflection intensity at the target flight time was enhanced 22.3 times compared with the initial pattern by the optimization process. Furthermore, the signal transmitted through the target area was enhanced 10.4 times; in other words, the light energy delivered by the optimized pattern to the target was about 10 times greater than that delivered by the random input.

 

Fig. 4. Experimental demonstration of the enhanced light energy delivery by the time-gated iterative phase conjugation. (a) Temporal response of the reflection intensity on the logarithmic scale for the initial random input (blue curve) and after 50 iterations (N = 50) of phase conjugation (red curve). (b) The reflection enhancement factor ${\eta _{re}}$ (red curve) and transmission enhancement factor ${\eta _{tr}}$ (blue curve) at the target area, which is indicated by yellow circles in (e) and (f), as a function of the number of iterations. (c, d) Intensity images of the time-gated reflection for the initial random input and after 50 iterations, respectively. Scale bar, 10 µm. (e, f) Intensity images of transmission for the initial random input and optimized solution. Scale bar, 10 µm.

Download Full Size | PDF

Intensity images of the initial and optimized reflection signals at $t = \tau $ are shown in Figs. 4(c) and 4(d), respectively. After the iterations, the light wave became more concentrated and intensified near the center of the target object’s location. In addition, the speckle size of the reflected wave was slightly increased from 0.85 µm to 1.3 µm. When the wave is focused to the target object, the backscattered signal emanates from the small target area and thus develops much less speckles. Therefore, the enlargement of the speckle grain is an additional evidence of wave focusing to the target. The initial and optimized intensity images measured from the transmission side of the target object are shown in Figs. 4(e) and 4(f), respectively. As the silver disk target had a transmittance of about 70%, the transmission intensity at the target object indicated by the yellow circles is proportional to the energy delivered to the target. From the transmission images, we found that the average intensity at the target object increased by about 10 times after the optimization.

3.2 Comparison of different optimization methods

We compare the time-gated iterative phase conjugation (PC) method with other optimization methods such as the time-gated eigenchannel (EC) method and iterative feedback optimization (FO) method. The three methods have similar performance in terms of the enhancement factor, but there are significant differences in their operation times. As will be shown in the following, the PC method is orders of magnitude faster than the other methods in reaching the optimal solution. Therefore, it is the most relevant for applications involving sample fluctuations and mechanical perturbations.

In the EC method, a time-gated reflection matrix is measured, and singular value decomposition (SVD) of the measured matrix is performed. Light is then coupled to the eigenchannel with the largest eigenvalue. Consequently, the light energy delivered to the target is enhanced. The EC method guarantees the optimal solution as long as the sample is stable. However, there is long temporal gap between the measurement of a matrix and wavefront shaping since the method requires matrix acquisition and computation of singular value decomposition. If the sample moves or vibrates, the obtained solution is gradually decorrelated from the optimal solution.

In the FO method, the illumination wavefront is modulated to optimize the intensity of the time-gated reflection. The optimization is achieved by the three-step phase shifting of the randomly chosen half of the illumination area. Owing to limitations of the algorithm, this method requires more iterations than the EC method for finding the optimized pattern.

In order to compare the operation time and susceptibility to sample perturbations among the three methods, we performed numerical simulations using the experimentally measured time-gated reflection matrix $\boldsymbol{R(\tau )}$. This matrix was obtained using the same setup that was used for the PC method. We used the SLM to scan the angle of illumination of the sample, and the time-gated backscattered wave was imaged by the camera (CAM_R) for each illumination angle over the recording field of view of 40 µm ${\times} $ 40 µm. We chose 1600 angular steps to uniformly cover this recording area up to the numerical aperture of 0.4. All these images were used to construct the reflection matrix. Since the reflection matrix describes the input-output response of the sample of interest, we can obtain the reflection image from $\boldsymbol{E^{out}} = \boldsymbol {R(\tau ){E^{in}}}$ for any given input field $\boldsymbol{E^{in}}$. Therefore, we can simulate the FO and PC methods by using the recorded reflection matrix.

Figure 5(a) shows the reflection enhancement depending on the number of images acquired when the target object and the scattering medium are static. In the EC method (blue curve), the enhancement reaches the maximum only after the full set of illumination angles is scanned for constructing the matrix. The enhancement in the FO method (green curve) increases slowly and converges after 10,000 iterations or so. The number of iterations required for the enhancement to reach the maximum in the FO method is about an order of magnitude more than that required in the EC method, but the maximum enhancement is still lower than that in the EC method. By contrast, the PC method (red curve) rapidly converged to the maximum enhancement after only 10 iterations, which is more than a hundred times faster than the other methods. Furthermore, the maximum enhancement factor in the PC method is almost identical to that in the EC method. Figures 5(b) and 5(c) show how the enhancement factor is disturbed by sample perturbations. We introduced artificial vibrations to the system by shifting the images constituting the reflection matrix at every iteration step. We introduced random image shifts that followed a Gaussian distribution with standard deviations of 0.2 µm for the weak vibrations and 1 µm for the strong vibrations. Simultaneously, we added a sample drift by shifting the images in the reflection matrix at a rate of 1 µm (5 µm) per 1000 steps for weak (strong) perturbations. The sample perturbation modeled in our study causes the decorrelation of the speckles similar to that observed in real biological tissues [31]. With the perturbations, the EC method showed a reduced enhancement factor that was about half and one-sixth of the maximum enhancement for weak and strong vibrations, respectively. Furthermore, the obtained solution was gradually decorrelated with the continuously moving sample. Similarly, the enhancement in the FO method approached about half of the maximum enhancement for weak vibrations and fell to the random input for strong vibrations. By contrast, the enhancement in the PC method quickly recovered because of the fast convergence rate, suggesting that this method can be much more robust to sample perturbations.

 

Fig. 5. Comparison of the eigenchannel method (EC), iterative feedback optimization (FO) method, and iterative phase conjugation (PC) method. (a)–(c) Enhancement factor η_re as a function of the number of image acquisitions in ideal conditions, in the presence of weak vibrations, and in the presence of strong vibrations, respectively. (d)–(f) Enhancement factors of the EC, FO and PC methods, respectively, depending on the size of the illumination area. The dashed lines indicate the achievable maximum enhancement factor for the EC method.

Download Full Size | PDF

For achieving efficient light energy delivery in practice, it is desirable that the working area set by the illumination area be large. In the case of the EC and FO methods, the operation time tends to be longer as the illumination area is increased. By contrast, the PC method is rather independent of the illuminating beam size as its convergence is mainly determined by the eigenvalue distribution. To confirm this, we numerically constructed a reflection matrix $\boldsymbol{R}_{\boldsymbol{sim}}$ for various illumination areas and compared the optimization time of the three methods. The reflection matrix was constructed using the relation $\boldsymbol{R}_{\boldsymbol{sim}} = \boldsymbol{T}_{\boldsymbol{s}}^{\textbf{T}}{\boldsymbol{R}_{\boldsymbol{target}}}{\boldsymbol {T}_{\boldsymbol{s}}}$ under the assumption that forward scattering is dominant. Here, ${\boldsymbol{T}_{\boldsymbol{s}}}$ is the transmission matrix of the scattering layer, and it is constructed from filtered random matrices [32]. ${\boldsymbol{R}_{\boldsymbol{target}}}$ denotes the reflection matrix of a thin target object in the space domain, and it is a diagonal matrix whose diagonal elements represent the amplitude reflectance of the target. $\boldsymbol{T}_{\boldsymbol s}^\textbf{T}$ denotes for the return trip [33–35]. The size of the target object was assumed to be identical, and it was set by maintaining the number of nonzero diagonal elements of ${\boldsymbol{R}_{\boldsymbol{target}}}$. The number of iterations required for reaching the maximum enhancement was approximately proportional to the illumination area for the EC (Fig. 5(d)) and FO (Fig. 5(e)) methods since the required number of iterations depended on the number of independent input channels, which was proportional to the illumination area. By contrast, the required number of iterations in the PC method (Fig. 5(f)) changed little with an increase in the illumination area. In the case of the PC method, the convergence depended on the power of the ratio among the eigenvalues, not the number of input channels. Therefore, the PC method is more useful than the other methods for a large illumination area.

4. Discussion

We have presented a method to focus multiple scattered waves on a target behind a scattering medium; the method is a combination of the time-gated detection and iterative phase conjugation of backscattered waves. The method enhanced the light energy delivered to a target hidden behind a 200-µm-thick mouse skull by more than ten times compared with the initial random input. Only 10 iterations were required to achieve the maximum enhancement, which is more than two orders of magnitude smaller than the number of iterations required by the time-gated reflection matrix method or by iterative FO of the time-gated reflection intensity. In other words, the proposed method is faster than the existing methods by about a hundred times or more. Through additional modeling and analysis, we verified that the proposed method is less sensitive to sample perturbations. We also confirmed that the number of images required for optimization remains almost unchanged with an increase in the illumination area in the proposed method. These advantages render the proposed method suitable for practical applications where efficient light energy delivery is critical.

In the present experiment, the slow refresh rate (∼5 Hz) of the liquid-crystal based SLM is a major limiting factor for the optimization of the speed (slightly over 0.2 s/step). The speed can be increased substantially by using a high-speed wavefront shaping device such as a digital micromirror device (∼20 kHz) [36], a micro-electro-mechanical system (∼4 kHz) [37], or an acousto-optic deflector (∼40 kHz) [38]. With these devices and a high-speed camera, the iteration rate of the proposed method can be increased to 1 kHz or more. The light energy enhancement factor can be optimized in 10 ms for 10 iteration steps, which is sufficiently fast for most of the in vivo applications in biosensing and light stimulation such as optogenetics [39] and optical therapy [40].