# Optimal spatiotemporal focusing through complex scattering media

###### Abstract

We present a new approach for spatiotemporal focusing through complex scattering media by wave front shaping. Using a nonlinear feedback signal to shape the incident pulsed wave front, we show that the limit of a spatiotemporal matched filter can be achieved, i.e., the wave amplitude at the intended time and focus position is maximized for a given input energy. It is exactly what is also achieved with time-reversal. Demonstrated with ultrasound experiments, our method is generally applicable to all types of waves.

## I Introduction

Most imaging systems are based on the ability to focus a wave beam inside the area of interest. In the context of echographic imaging, focusing of ultrasound in the human body can be achieved with a transducer array and electronic delay lines: the same pulsed waveform is sent from each transducer with the appropriate delay making all the waveforms converge in synchrony at the desired focus location. However, that principle is not practicable any more as soon as the sample thickness becomes larger than the mean free path, i.e., the average distance between two random scattering events. In such a strongly scattering medium, it has been shown with ultrasonic waves that spatiotemporal focusing can be achieved using time-reversal: a training pulse is first sent from a source located at the intended focal point, travels through the scattering medium and is captured at a transducer array, the ’time reversal mirror’Fink (1997). The waveforms received on the time reversal mirror are flipped in time and sent back, resulting in a wave converging at the desired focus location. Time-reversal focusing is optimal is the sense that it achieves a spatiotemporal matched filter, a term well-known from signal processing; for a given input energy, the amplitude of the pulse at the focal spot is maximal Tanter et al. (2001). Time-reversal has also been implemented for microwaves Lerosey et al. (2007). For optical waves, time-reversal has not been demonstrated yet, since optical time-reversal mirrors are challenging to realize.

Therefore in optics a different approach, named ’wave front shaping’ (WFS), has been taken to steer light through strongly scattering media Vellekoop and Mosk (2007). The principle is based on spatial phase modulation by Spatial Light Modulators (SLM), which enables one to modulate the complex amplitude of the independent modes of the incident wave front. The intensity in the selected output mode, e.g. the focal spot, is used to match the wave front to the scattering medium by an adaptive algorithm. Owing to their enormous amount of pixels, state-of-the-art spatial light modulators allow achieving transmission of substantial amounts of the input energy Vellekoop and Mosk (2008).

A big advantage of adaptive wave front shaping approaches lies in the fact that a direct access to the amplitude at the focus is not required, such that any type of intensity probe can be used. For instance the fluorescence from dye molecules allows to focus on objects hidden inside a complex medium Vellekoop et al. (2008). The concept has been applied to other types of waves, such as surface waves Gjonaj et al. (2011) and ultrasonic waves in the single scattering regime for medical imaging Herbert et al. (2009) and MR-guided ultrasonic therapy Larrat et al. (2010).

Initially introduced in the monochromatic regime, WFS has recently been extended to broadband light for spatiotemporal focusing Aulbach et al. (2011). In this work it was shown that spatial control of the input waves is sufficient to achieve spatiotemporal control of the scattered waves. The key idea is to use coherence gating for optimizing a single specific predetermined point in time. Performing the optimization at a specific time in the speckle amounts to combine all the scattering paths that have about the same length. This principle also works efficiently without defining a focusing time in the case when the temporal delay spread induced by the scattering medium is comparable to the input pulse duration Katz et al. (2011). In both cases, the missing control of temporal degrees of freedom is compensated by manipulating spatial ones. In a complementary way, temporal control of the input light based on pulse shaping Weiner (2000) can be used to achieve spatiotemporal focusing by spatially localized phase compensation McCabe et al. (2011). Both approaches exploit the fact that temporal and spatial degrees of freedom are mutually convertible in a complex medium Lemoult et al. (2009).

In the present paper, we demonstrate how to fully control both spatial and temporal degrees of freedom to achieve optimal spatiotemporal focusing for a broadband ultrasonic wave propagating through a complex multiple scattering medium. Like in the context of time-reversal focusing, we use the term ’optimal’ in the sense of a spatiotemporal matched filter, i.e., for a given input energy, our method maximizes the wave amplitude at the intended time and focus position. The difference with time reversal is that no amplitude measurement is required at the focus location. Instead spatial and frequency resolved wave front shaping can be achieved with a nonlinear feedback intensity signal and is capable to reach optimal spatiotemporal focusing.Thus our concept should be generalizable to all types of waves.

## Ii The spatiotemporal matched filter approach

### ii.1 Optimal focusing by time reversal

The experimental system is shown in Fig. 1. Our goal is to maximize the signal amplitude at some focus point in the receiver plane at time for a given energy radiated from points in the emitter plane. In signal processing, this is known as the spatiotemporal matched filter approach.

##### Spatial matched filter

At first we recapitulate the conditions for optimal spatial focusing in the sense of a spatial matched filterTanter et al. (2001). Let be the input signals in the emitter plane, and correspondingly the signals received in the image plane. The propagation through the medium from the emitter points to the control points at one frequency is described by:

(1) |

where is the propagation operator, or transfer matrix, between the points at a given frequency . The amplitude received at the focal point at a single frequency is given by

(2) |

where

(3) |

is the projection onto . The inequality of Cauchy-Schwartz sets the upper bound for this expression:

(4) |

Equality holds, if , or

(5) |

The maximum amplitude at the focal point is reached when all are proportional to the complex conjugate of the respective transmission coefficient. This situation corresponds to a phase conjugation experiment or a time-reversal experiment at a single frequency.

##### Temporal matched filter

Let us denote by the impulse response between emitter and receiver , which is the Fourier Transform of the respective transfer function in the frequency domain . When one aims at maximizing the amplitude received at time , it is well-known that has to be sent from . The signal at then writes as

(6) |

which is the autocorrelation of the impulse response, or the inverse Fourier transform of the power spectrum, respectively. It is an even function with its maximum value at when all frequency contributions add in phase. Hence the magnitude at is directly proportional to the transmitted energy. As is the Fourier Transform of , it also achieves the spatial matched filter (Eq. 5).

This signal exactly corresponds to the one transmitted in a broadband time-reversal experiment. In the first step, the control point in the image plane behaves like a source, such that the wave field recorded by the two-way transducers in the emitter plane is given by In the second step, the recorded fields are emitted in a time-reversed manner. Under assumption of reciprocity the emission writes as and thus the focus is given by Eq. 6. Hence, time-reversal exactly achieves a spatiotemporal matched filter since it implements a spatial filter at each frequency of the incident pulse and the correct phase relation between the frequency components is intrinsically recovered.

### ii.2 Optimal focusing by wave front shaping

In the following we lay out the wave front shaping algorithm and investigate its performance as a matched filter. In order to achieve optimal spatiotemporal focusing of a broadband pulse, we perform the wave front shaping algorithm in two stages. During the first phase, we achieve optimal spatial focusing. Based on this result, optimal temporal focusing is achieved after the second phase.

#### ii.2.1 Optimal spatial focusing

The direct implementation of wave front shaping as applied in optics to monochromatic acoustic waves would be described as follows: at first an arbitrary wave front is sent from the emitter plane. The intensity is then measured at the intended focal point. This can be done directly in measuring the pressure field with a hydrophone or, in the case of a soft elastic medium, in measuring the displacement induced at the focal point by the radiation force, which is proportional to the acoustic intensity Herbert et al. (2009). Then the phase of the first transducer is cycled between and while the intensity is recorded. The phase value for which intensity at focus is found maximum is then stored in the memory and the same operation is repeated for each transducer. Once the optimal phase has been stored for each transducer, the wave front can be synthesized by exciting each transducer with the determined optimal phase.

##### Hadamard method

Instead of using such an optimization scheme element by element, we perform a basis transformation to construct virtual transmitters , defined as linear combinations of real ones, which greatly improves the sensitivity of the method and reduces the number or iterations Herbert et al. (2009). The coefficients of the combination are taken as the columns of an by Hadamard matrix , with elements . This choice ensures that the amplitude transmitted from each transducer is maximal. The basis transformation is performed by

(7) |

and the inverse transformation by

(8) |

The field in the receiver plane is given by

(9) |

where is the transfer matrix between the emitters in the hadamard basis and the receivers in the canonical basis. For the optimization algorithm, we arbitraritly chose the first column as the reference. 4 coded beams are sent into the scattering medium as:

(10) |

(11) |

After each transmission the resulting intensity at the focal point is given by

(12) |

where is the corresponding transmission coefficient from the transfer matrix . From the set of measured intensities we are now able to calculate the optimal wave front. The complex amplitude to be transmitted from each virtual transducer to optimize the intensity is given by the following expressions:

(13) |

where the compensate for the phase differences between the contributions at the focal point. The choice of ensures that WFS achieves a spatial matched filter. Consequently, all are equivalent to their corresponding up to the global phase difference ,

(14) |

where selects the -th row of the transmission matrix. Finally, a change of basis leads back to the transducer basis in the emitter plane

(15) |

showing, that the wave front shaping algorithm fulfils the condition of the spatially matched filter (Eq. 5). The emission is equivalent to the time reversal experiment. With the definition of the Hadamard matrix, , the prefactor N arises since we omitted the equivalent transformation in the receiver plane. Since this proportionality factor does not impair any conclusions (see Eq. 5), we omit it in the following.

##### Enhancement factor

After optimization for focusing on point , the emitted signal from transducer is

(16) |

where the denominator normalizes the transmission. The received signal at focus point is

(17) |

The average energy received at the focus point is

(18) |

where the brackets denote the statistical average. We assume that the transmission coefficients are independent random variables and follow a circular Gaussian distribution. Without optimization, the normalized initial emission is

(19) |

The average received signal on point is

(20) |

The average enhancement, which also gives the signal to background ratio to other non-optimized modes, is therefore given by

(21) |

Initial works in optics Vellekoop and Mosk (2007) used spatial light modulators which were limited to phase-only control of the emitted amplitudes, such that the emission is given by

(22) |

The amplitude at the focus is

(23) |

In this case the received energy at focus point m is

(24) |

Hence the enhancement is lowered by the well-known prefactor

(25) |

##### Algorithm for optimal spatial focusing

During the first phase, we achieve spatial focusing for all frequency components contained in emission spectrum of the broadband pulse with its bandwidth . Since the transmission coefficients (Eq. 1) are random variables of frequency with a correlation length , the optimal emission (Eq. 15) has to be determined frequency-resolved for intervals . Consequently, the Hadamard algorithm is performed frequency by frequency in steps of , each time using a narrowband pulse with a Gaussian spectrum of a bandwidth for the emission. When afterwards the full broadband signal with the optimal coefficients is emitted, spatially matched focusing is achieved. Contrary to the time-reversal experiment, temporal focusing can not yet be expected. For this, all frequency components have to be in phase, which is intrinsically fulfilled by time reversal. However, after the wave front shaping optimization, the unknown phase factor remains.

#### ii.2.2 Optimal temporal focusing

During the second phase, the spatially focused signal is additionally focused in the time-domain. For that, the phase factor needs to be determined and its conjugate equally multiplied to all emissions . Such spatially invariant phase filtering leads only to a temporal redistribution of the signal at the focal point. The time-integrated intensity at the focus is unaltered by the filtering process, excluding it as an appropriate feedback signal. Instead of detecting the linear intensity at the focal spot, we place a nonlinear detector which allows us to measure the time-integrated second harmonic intensity. This detector is treated in more details in the following paragraph. The second harmonic energy is sensitive for a temporal redistribution of the pulse energy. Used as a feedback signal, it reaches its maximum when the transmitted signal is shortest in duration, corresponding to the so-called Fourier-limit, when the waves at all frequencies are in phases at the focus. This technique has similarly been used in optics for the compression of dispersed ultrashort laser pulses Meshulach et al. (1997).

The most straight-forward approach is an adaptive algorithm equivalent to the spatial optimization scheme described in the beginning of the section. The frequency range of the broadband emission spectrum is subdivided into intervals . Within each iteration of the optimization, we emit the full signal as determined for optimal spatial focusing, , but modified with a spatial invariant phase mask . Step by step, all frequency intervals are addressed successively. Before the first iteration, the phase mask is . In step , the phase mask is modified only at the frequency . We perform several emissions with the phase increasing from to , while the second harmonic energy is recorded. As a function of each , follows a cosine behaviour with a phase offset , which depends on the actual phases of all other frequency components. The optimum phase offset is found when reaches its maximum. The modified phase mask is subsequently used in the next step address, in which the next frequency interval is handled in the same way. This way, the independent frequency components of the signal are set in phase at the focus step by step, i.e. the phase mask converges towards the conjugate of .

#### ii.2.3 Nonlinear detector

For the first phase, to achieve spatial focusing, the detection of the linear intensity would be sufficient. However, for an experimentally convenient approach, we use the nonlinear detector for both spatial and temporal focusing. We consider a system with a second-order nonlinearity, which could be realized, e.g., by a nonlinear scatterer, such as microbubbles in acoustics Leighton (1994) or nonlinear nanocrystals Hsieh et al. (2009) in optics, combined with a time-integrated detection of the scattered second-order response. We calculate the detector response by

(26) |

where is an arbitrary prefactor and the amplitude at the focus point in the time domain. For the narrowband pulses emitted in the first phase, we can assume that the linear intensity can be readily obtained by for the calculation of the transmission coefficients (Eq. 13).

## Iii Experiment

### iii.1 Transfer matrix measurement

As a first step of the experiment, we recorded the transfer matrices of a multiple scattering medium over a broad range of frequencies. The experimental setup consists of two identical transducers, one used as emitter and the other as receiver respectively, placed in water on opposite sides of the medium (Fig. 1). The latter is an arrangement (thickness L = 80 mm) of parallel steel rods with density and a rod diameter 0.8 mm. Both the emitting and the receiving transducer can be translated on a line parallel to the medium, on which we defined each 128 emitter points and 128 receiver points spaced by 0.4 mm. From each emitter point, we sent an ultrasonic pulse at a central frequency of 3.5 MHz with a relative bandwidth of 31%. In this frequency range the mean free path of the medium is as determined in Derode et al. (2001). The transmitted signals were recorded on all receiver points with a sampling frequency of 10 MHz. From the Fourier transform of the temporal signal, we obtain the transfer matrices in the considered frequency range (Eq. 1). We determined the correlation of the transmission coefficients both in the spatial and in the frequency domain by

(27) |

(28) |

where the average is performed over the denoted matrix coefficients in Eq. 27 and over the first matrix entry and frequencies in Eq. 28. From the full width at half maximum of the correlation functions, and , we can determine that we have independent frequencies per frequency unit and independent emitter points available van Beijnum et al. (2011).

With the transfer matrices at hand, we can calculate the linear response of the system for arbitrary signals emitted either from the ’emitter’ or ’transducer’ side. For all following experiments, we emitted pulses with a Gaussian spectral function of 10% bandwidth around the central frequency from all transducers. A typical signal in the receiver plane is plotted in Fig. 2(a) and Fig. 3(a). The energy is spread widely both temporally and spatially. The ballistic part of the wave has disappeared, which confirms that the medium is strongly scattering.

### iii.2 Spatial focusing

Starting from this emission we performed a frequency-resolved wavefront shaping optimization as described in section (II.2.1). The bandwidth of the narrowband pulses is chosen to to be smaller than the correlation . In two experiments we use both emission with optimal phase and amplitude (Eq. 16) and the emission with optimal phase only (Eq. 22). The deposed energy in the receiver plane for both experiments is plotted in Fig. 4. It is calculated by integrating over time the square of the wave amplitude at the receiver points (see Fig. 2(b) for the signal obtained in optimizing both amplitude and phase). The optimization leads to a spatial focusing in both cases. For phase and amplitude control, the increase in energy at the focal spot is , which corresponds perfectly to the number of estimated independent emitters (Eq. 21). The enhancement by the time-reversal focus is slightly higher, indicating remaining deviations between the wave front shaping emission and the time-reversal emission. In further calculations, we observed that a further reduction of eliminates these deviations, which shows that the single-frequency TR emission and the WFS emission are identical (see Eq. 16 and Eq. 5).

For phase-only optimization, we observe an enhancement lowered by a factor 0.71 (). The factor (Eq. 25) is not fully reached, which we attribute to the fact, that due to the geometry of the setup not all emitters contribute equally at the focal point, effectively reducing the number of emitter points Aulbach et al. (2011). The resulting time-resolved signals for the first case are shown in Fig. 2(b) and Fig. 3(b). The energy is still spread temporally, since temporal focusing is impeded by the remaining random phase relation between independent frequency components (Eq. 14). In the next step, we determine and correct for this phase factor.

### iii.3 Spatiotemporal focusing

Starting with the emission for spatial optimization, we perform frequency phase filtering as described in section (II.2.2). After two cycles over all frequencies, the received signals are not only spatially, but also temporally focused (Fig. 2(c) and Fig. 3(c)). The resulting peak intensity can be put in context with the degrees of freedom available for the optimization. We have independent emitters at hand. In the frequency domain, degrees of freedom are available calculated from the frequency correlation, and a bandwidth of . Hence, the total number of degrees of freedom is . We estimate the intensity enhancement to , comparing the peak intensity to the intensity before optimization (obtained from the transmission averaged over the transducers around the focus around the peak of the diffuse transmission, smoothed by a window to remove temporal speckle). This number corresponds well to the order of magnitude of the number of degrees of freedom.

As a reference, we performed an equivalent time-reversal experiment, which intrinsically reaches spatiotemporal focusing (Fig. 2(d) and Fig. 3(d)). The quality of the focusing from our adaptive wave front shaping and time-reversal are nearly identical (Fig. 5). The small deviations are a result of small remaining phase differences which can be minimized by a reduction of the frequency steps during the spatial optimization (see above) and further iterations of the adaptive phase filtering process. As an alternative to the step-by-step algorithm we used here, so-called stimulated annealing algorithms Meshulach et al. (1997) or genetic algorithms Brixner et al. (1999) should equivalently find the optimal phase relation. These algorithms are known to be robust in the presence of noise, but have the disadvantage to turn inefficient for an increasing number of degrees of freedom.

## Iv Conclusions

In conclusion, we have presented a new approach for spatiotemporal focusing through complex scattering media by feedback-based wavefront shaping. In contrast to previous works, our approach is capable of reaching the limit of a spatiotemporal matched filter; the maximum possible amount of the emitted energy is deposed at the target. We showed that to achieve this, phase and amplitude of the emission need to be controlled at each frequency. In contrast to the classical time reversal experiment that achieves the ideal focus, here the direct access to the field amplitude at the target is not required as only the intensity associated to a nonlinear response needs to be detected. This point of is particular interest for many applications in wave physics where the accessible information is restricted to field intensity. For example, in optics, our method can be realized by fluorescent dye molecules as they are used for two-photon microscopy. In ultrasound therapy, it could optimize both spatial and temporal focusing of therapeutic beam through the skull bone thanks to MR radiation force imaging. Being generally applicable to all types of waves, we believe that our method is promising for a wide range of applications in imaging and sensing, and for the control of wave propagation in combination with complex media such as new metamaterials.

This work is part of the Industrial Partnership Programme (IPP) Innovatie Physics for Oil and Gas (iPOG) of the Stichting voor Fundamenteel Onderzoek der Materie (FOM), which is supported financially by Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO). The IPP MFCL is co-financed by Stichting Shell Research.

## References

- Fink (1997) M. Fink, Phys. Today 50, 34 (1997).
- Tanter et al. (2001) M. Tanter, J. Aubry, J. Gerber, J. Thomas, and M. Fink, J. Acoust. Soc. Am. 110, 37 (2001).
- Lerosey et al. (2007) G. Lerosey, J. de Rosny, A. Tourin, and M. Fink, Science 315, 1120 (2007).
- Vellekoop and Mosk (2007) I. M. Vellekoop and A. P. Mosk, Opt. Lett. 32, 2309 (2007).
- Vellekoop and Mosk (2008) I. M. Vellekoop and A. P. Mosk, Phys. Rev. Lett. 101, 120601 (2008).
- Vellekoop et al. (2008) I. M. Vellekoop, E. G. van Putten, A. Lagendijk, and A. P. Mosk, Opt. Express 16, 67 (2008).
- Gjonaj et al. (2011) B. Gjonaj, J. Aulbach, P. M. Johnson, A. P. Mosk, L. Kuipers, and A. Lagendijk, Nat. Photonics 5, 360 (2011).
- Herbert et al. (2009) E. Herbert, M. Pernot, G. Montaldo, M. Fink, and M. Tanter, IEEE Ultr. Ferro. Freq. Cont. 56, 2388 (2009).
- Larrat et al. (2010) B. Larrat, M. Pernot, G. Montaldo, M. Fink, and M. Tanter, IEEE Ultr. Ferro. Freq. Cont. 57, 1734 (2010).
- Aulbach et al. (2011) J. Aulbach, B. Gjonaj, P. M. Johnson, A. P. Mosk, and A. Lagendijk, Phys. Rev. Lett. 106, 103901 (2011).
- Katz et al. (2011) O. Katz, E. Small, Y. Bromberg, and Y. Silberberg, Nat. Photonics 5, 372 (2011).
- Weiner (2000) A. M. Weiner, Rev. Sci. Instrum. 71, 1929 (2000).
- McCabe et al. (2011) D. J. McCabe, A. Tajalli, D. R. Austin, P. Bondareff, I. A. Walmsley, S. Gigan, and B. Chatel, Nat. Commum. 2, 447 (2011).
- Lemoult et al. (2009) F. Lemoult, G. Lerosey, J. de Rosny, and M. Fink, Phys. Rev. Lett. 103, 173902 (2009).
- Meshulach et al. (1997) D. Meshulach, D. Yelin, and Y. Silberberg, Opt. Commun. 138, 345 (1997).
- Leighton (1994) T. G. Leighton, The acoustic bubble (Academic Press, London, 1994).
- Hsieh et al. (2009) C. Hsieh, R. Grange, Y. Pu, and D. Psaltis, Opt. Express 17, 2880 (2009).
- Derode et al. (2001) A. Derode, A. Tourin, and M. Fink, Phys. Rev. E 64, 036606 (2001).
- van Beijnum et al. (2011) F. van Beijnum, E. G. van Putten, A. Lagendijk, and A. P. Mosk, Opt. Lett. 36, 373 (2011).
- Brixner et al. (1999) T. Brixner, M. Strehle, and G. Gerber, Appl. Phys. B 68, 281 (1999).