Spatiotemporal focusing in opaque scattering media by wave front shaping with nonlinear feedback

Jochen Aulbach; Bergin Gjonaj; Patrick Johnson; Ad Lagendijk

doi:10.1364/OE.20.029237

1. Introduction

Achieving optimal control of light propagation in any type of complex scattering medium is highly desirable for applications such as in imaging, micro- and nanomanipulation, coherent control and spectroscopy. However, in numerous systems of interest of current research, for example in thick biological specimens, random multiple scattering is a major obstruction for the performance of conventional optical techniques.

Recently, Vellekoop and Mosk demonstrated that random multiple light scattering can be exploited to increase rather than hamper the control of light propagation [1]. The approach, named ’wave front shaping’ (WFS) is based on spatial modulation of the complex amplitude of the light incident on the medium. With an adaptive algorithm the optimal wave front is determined, which matches the medium such that coupling of light to a desired output mode is maximized, e.g. a focal spot behind the medium. In this way the scattering medium can be turned into a lens to focus light for trapping nanoparticles or for high-resolution imaging as recent applications have demonstrated [2, 3]. Requiring only an intensity measurement as feedback signal makes the approach powerful and versatile, enabling focusing on fluorescent beads embedded inside a complex medium [4] without the need for phase-sensitive measurement techniques.

For acoustic waves and for microwaves, related experiments have demonstrated broadband focusing in complex scattering media using a time reversal mirror [5]. Based on this technique, Lerosey et al. [6] focused electromagnetic waves onto a small antenna inside a dense assembly of resonant scatterers with precision far below the wavelength of the radiation. However, for optics, the realization of an optical time reversal mirror remains challenging, requiring both an accurate broadband measurement and the synthesis of a complex electric field. The existing optical analog, phase conjugation, is effectively limited to monochromatic light [7–9].

While the initial optics experiments were performed with monochromatic light, WFS was recently extended to the broadband regime. Spatial wave front shaping in combination with a scattering medium allows both spatial and temporal control of the transmitted scattered light using optical gating as feedback [10] or pulse compression through the medium using two-photon fluorescence [11].

However, spatiotemporal focusing onto nanoparticles embedded inside a strongly scattering medium has not been demonstrated yet. For this, nanoparticles with a nonlinear optical response, such as two-photon fluorescence or second-harmonic generation [12], need to be embedded inside the medium. Candidates for these probe particles are fluorescent dyes such as are commonly used for two-photon microscopy or quantum dots [13]. Especially interesting is second harmonic generation (SHG) in nanoparticles [14], and in nanocrystals from wide band-gap materials with a high SHG efficiency in particular [15, 16]. Providing a stable and coherent signal, and flexibility for use in a wide spectral range, these class of particles recently gained considerable attention as markers for novel microscopy techniques [17] and makes them ideal probes for WFS experiments.

Here we experimentally demonstrate spatiotemporal focusing on single nanocrystals embedded inside a strongly scattering medium. Our approach is based on wave front shaping of short pulses, using second harmonic generation inside the target nanocrystals as feedback signal. We develop a model both for the achieved pulse duration at the particle position as well as the observed enhancement of the feedback signal which is in good agreement with the experiment. Our approach has implications for applications in which for control of light propagation in complex media is required and provides a powerful tool to study light propagation in the presence of strong multiple scattering.

This paper is structured as follows. At first we describe the experiment in section 2. It is followed by a detailed and extensive theory section 3 before we provide the experimental results in section 4.

2. Experiment

The general idea of the experiment can be described as follows. The pulse front impinges on the surface of a thick multiple scattering medium. The light reaching positions inside the medium at depths larger than a few mean free paths has been multiply scattered such that the spatial structure of the incident beam is lost and temporally the pulses are elongated to random speckle pulses. The goal is to partially reverse this effect and to focus the light spatiotemporally on a single point in the medium at a single moment in time. For this, we apply spatial wave front shaping, using second-harmonic active nanocrystals to deliver the feedback signal. The particles act as a local probe of the light intensity, converting the present light from the fundamental wavelength to the second-harmonic wavelength, which can be readily detected separately from the fundamental light. Since the particle is a local probe, we can expect that the algorithm will locally maximize the light intensity. The local second-harmonic generation not only acts as a spatial filter, but is also biased towards shorter pulse durations, leading to additional temporal focusing.

The experiment is illustrated in Fig. 1. The light from the Ti:sapphire laser (repetition rate 80MHz, output power 1.2W, center wavelength 841nm, bandwidth 19nm, pulse duration 67fs) is sent through a Michelson-type interferometer. One interferometer arm is of fixed length, the other has a variable delay (Newport XMS100). Subsequently the beam is expanded 6:1 (not shown) to fill the surface of a two-dimensional phase-only spatial light modulator (SLM, Hamamatsu X10468). The SLM pixels are grouped into N independent segments each of which induces a controllable phase shift. Two lenses project the SLM surface onto the back focal plane of an air-immersion objective (numerical aperture (NA) 0.9, x100). The focal plane of the objective is a conjugate plane of the SLM surface, with a total magnification of 1/300. The sample is placed in the focal plane of the objective. The transmitted light through the sample is collected by an oil-immersion objective (magnification 100x, NA 1.4). After filtering (short-pass filter 675nm and bandpass filter 420nm), only the generated second harmonic is imaged by an electron multiplying CCD camera (Andor DV885, frame rate 0.25 – 2Hz, total magnification 93.5/1). Optionally the laser light can be directly coupled into a Fourier transform spectrometer instead of entering the setup (not shown).

Fig. 1 Schematic illustration of the experiment. (Ti:Sa) Ti:sapphire laser, (MI) Michelson-type interferometer, (SLM) spatial light modulator, (NA) numerical aperture, (EMCCD) electron-multiplying CCD camera.

Download Full Size | PDF

The sample, depicted in the inlay of Fig. 1, consists of a disordered layer of silica beads with a rough surface which covers sparsely distributed second-harmonic active nanocrystals on a standard glass coverslip (thickness 0.17mm). The layer of silica spheres (diameter 400nm to 500nm) has been fabricated by spray-coating a colloidal dispersion in ethanol, which after drying forms a rough surface with a layer thickness between L = 25 μm and L = 50 μm. The mean free path of this medium (l = 3.5 μm) has been determined by measuring the enhanced backscattering cone. The second-harmonic active nanocrystals consist of barium titanate (BaTiO₃) with a tetragonal crystal structure and an average diameter of 200nm. Prepared using the method described in [18], about 90% of the nanocrystals are not clustered and are well-separated from other particles, which we confirmed by imaging a sample of nanocrystals without the scattering medium. During the application of the silica dispersion with the spray-coating technique most of the barium titanate nanocrystals are detached from the glass coverslip, but remain in the vicinity of the substrate interface within the first few layers of silica beads (rather than diffusing deep into the silica layer) such that they can be identified as individual sources of second harmonic generation.

For every wave front shaping experiment presented in this paper we follow the following procedure. At first, we take the average image of a large number of images of random wave fronts (with a segmentation of SLM into 768 segments), which serves as a reference measurement for the obtained enhancement of the feedback signal. Then a flat phase pattern is applied to the SLM, and we perform a spatially resolved second-harmonic auto-correlation measurement [19] by recording images while the movable arm of the Michelson interferometer is scanned in discrete steps. The step resolution of the delay stage (10nm) allows us to record the full interferometric auto-correlation trace. To speed up the data acquisition we measure only parts of this trace, applying a two-step algorithm. After a large step (2.5 μm or 5 μm), the stage samples one optical cycle (841nm) in eight steps, from which the cycle-averaged intensity auto-correlation is calculated. Then the wave front shaping algorithm is started. We use a sequential algorithm; it addresses segment by segment one by one, scanning the phase from φ = 0 to φ = 2π in N_φ = 8 steps. For each step we record the time-integrated second-harmonic radiation from a single nanocrystal retrieved from the camera image as feedback signal. The phase value which maximizes the feedback signal is determined by fitting the measured behavior of the feedback signal vs. phase with a cosine function. We perform several consecutive sequences over all segments of the SLM, starting with a low segmentation of the SLM (48 segments) which is increased twice after the first sequences (to 192 and 768 respectively), always starting the new sequence with the phase pattern obtained from the previous one. After the algorithm has finished, another auto-correlation measurement is performed.

3. Theory

3.1. Intensity-intensity auto-correlation of speckle pulses

Our experimental setup enables us to measure the second-harmonic autocorrelation AC(τ) [19] spatially resolved at the particle position as a function of delay τ introduced in the Michelson interferometer. For the discussion, we consider ensemble averaged intensity auto-correlation (AC) with background, which are

〈 A C (τ) 〉 = \int d t 〈 I^{2} (t) 〉 + 2 \int d t 〈 I (t) I (t + τ) 〉 .

where I(t) denotes the intensity at the fundamental frequency and the brackets denote the ensemble average over all realizations of the scattering medium. For short pulses the observed contrast ratio between the maximum and the background is AC(τ = 0) : AC(τ → ∞) =3:1, regardless of the exact temporal shape of the intensity I(t).

Characteristic for the time-dependent average diffuse transmission 〈I_d(t)〉 through a strongly scattering disordered layer is an exponential decay with the diffuse decay time τ_d (or Thouless time [20]) as characteristic time scale. It is approximately related to the thickness of the medium by τ_d = L²/6D, where D is the diffusion constant [21]. An exact expression for the temporal behaviour of 〈I_d(t)〉 for the diffuse transmission through a slab is given in [22]. For a single experimental realization, the transmitted intensity I_d(t) is dominated by temporal speckle. Characteristically for such a speckle pulse, it has a limited temporal coherence time τ_c, which is given by the Fourier transform of its intensity spectrum. Therefore, the speckle pulses can be described as a random sequence of short pulses, each approximately of the duration τ_c, distributed randomly in a time window of the diffuse decay time of the medium τ_d.

The particular shape of the AC of a speckle pulse, observed after transmission through a slab of turbid material, can be described as follows. The maximum of the AC for τ = 0 decays on the time scale τ_c, since for larger delay times the intensities I_d(t) and I_d(t + τ) within a speckle pulse are uncorrelated random variables given by the Gaussian statistics of the transmission process. However, for thick scattering media with τ_d ≫ τ_c, the average intensity 〈I_d(t)〉 is approximately constant on timescales of τ_c. E.g. for τ > τ_c the value of the integral on the right-hand side can be approximated by separating the product into 〈I_d(t)〉〈I_d(t + τ)〉. Therefore the correlation term does not directly decay to the background level, but shows an intermediate regime resulting from the incoherent part of the speckle pulse. Since for Gaussian statistics 〈I_d(t)〉² = 0.5〈I_d(t)〉², the coherent peak is separately visible on top of the diffuse peak with a contrast ratio 3:2:1 over the background. Altogether, we can model the normalized AC of the speckle pulses by

\frac{〈 A C (τ) 〉}{〈 A C (τ = 0) 〉} = 1 + \frac{\int d t 〈 I_{d} (t) 〉 〈 I_{d} (t + τ) 〉}{\int d t {〈 I_{d} (t) 〉}^{2}} + \frac{\int d t I_{0} (t) I_{0} (t + τ)}{\int d t {(I_{0} (t))}^{2}} .

The second term on the right-hand side of Eq. (2) is the normalized contribution of the average diffuse transmission 〈I_d(t)〉, the third term is the normalized contribution from the non-scattered bandwidth limited short pulse I₀(t) provided by the laser.

3.2. Enhancement of the time-integrated second harmonic intensity

3.2.1. Enhancement vs. number of segments

The wave front shaping algorithm maximizes the time-integrated second-harmonic intensity generated at the particle position. In the following, we derive an expression for the average enhancement of this feedback signal which is achieved after wave front shaping.

In the first step, we derive the enhancement expected from a hypothetical monochromatic experiment with the second-harmonic intensity as feedback. The polarization at the second harmonic frequency P_2ω at the particle positioned at point b can be calculated by

P_{2 ω} = χ_{(2)} E_{b}^{2} = χ_{(2)} \sum_{a = 1}^{N} t_{a b} E_{a} \sum_{a^{'} = 1}^{N} t_{a^{'} b} E_{a^{'}} .

The complex transmission coefficients t_ab describe the propagation of the electric field E from SLM segment a to the particle position b. Here, we have made a number of simplifications; we neglect the polarization of the light, as well as the tensor character of the second-order susceptibility χ₍₂₎ of the nanocrystals. Furthermore we assume that the volume speckle inside the medium is larger than the particle size, and the particle radiates like a dipole at the second harmonic frequency. The radiated second harmonic power can be calculated by

W_{2 ω} = \frac{c k^{4} V^{2}}{12 π ε_{0}} {| P_{2 ω} |}^{2},

with the speed of light c, the particle volume V and the vacuum permittivity ε₀. In a wavefront shaping experiment, the electric field on each SLM segment a, E_a = A_ae^φ_a, is adapted with the aim that all contributions are in phase at the particle position, e.g. φ_a = −arg(t_ab) for all segments a. In our experiment, the phase-only SLM does not significantly modify the amplitudes A_a. As described above, we use a sequential algorithm for the WFS experiment; it addresses all segments one by one, scanning the phase term φ from 0 to 2π in discrete steps. Given that the segment of interest contributes the field E₁e^iφ ≡ t_abA_ae^iφ at the point of interest, while the sum of unmodified segments contributes

E_{2} \equiv \sum_{a^{'}, a^{'} \neq a}^{N} t_{a^{'} b} E_{a^{'}}

, the second harmonic power during the phase scan is given by

\begin{array}{l} W_{2 ω} (φ) & \propto & {| E_{1}^{2} |}^{2} + {| E_{2}^{2} |}^{2} + 4 {| E_{1} |}^{2} {| E_{2} |}^{2} + ℜ {4 E_{1} E_{2}^{*} {| E_{1} |}^{2} e^{i φ}} + \\ ℜ {2 E_{1}^{2} {(E_{2}^{*})}^{2} e^{2 i φ}} + ℜ {4 E_{1} E_{2}^{*} {| E_{2} |}^{2} e^{i φ}}) \end{array}

Typically |E₂| ≫ |E₁| holds, such that the last term on the right hand side dominates the varying terms, and the behavior of W_2ω(ϕ) is well described by a cosine function with a constant offset.

By evaluating the right hand side of Eq. (3), the ensemble averaged second harmonic power can be calculated,

〈 W_{2 ω} 〉 = 2 N (N - 1) {〈 {| t_{a b} |}^{2} 〉}^{2} + N 〈 {| t_{a b} |}^{4} 〉 .

For simplicity, we dropped the constants from Eq. (3) and Eq. (4) and set the field amplitudes to one.

The average resulting second harmonic intensity after WFS writes as

\begin{array}{l} {〈 W_{2 ω} 〉}_{wfs} & = & 〈 W_{2 ω} 〉 \\ + N (N - 1) {〈 {| t_{a b} |}^{2} 〉}^{2} \\ + 6 N (N - 1) (N - 2) 〈 {| t_{a b} |}^{2} 〉 {〈 | t_{a b} | 〉}^{2} \\ + 4 N (N - 1) 〈 {| t_{a b} |}^{3} 〉 〈 | t_{a b} | 〉 \\ + N (N - 1) (N - 2) (N - 3) {〈 | t_{a b} | 〉}^{4} . \end{array}

Taking only the highest order terms in N into account, the observed enhancement is approximately given by

η_{cw} \equiv \frac{{〈 W_{2 ω} 〉}_{wfs}}{〈 W_{2 ω} 〉} \approx 0.5 {(\frac{{〈 | A | 〉}^{2}}{〈 {| A |}^{2} 〉})}^{2} {(\frac{{〈 | t_{a b} | 〉}^{2}}{〈 {| t_{a b} |}^{2} 〉})}^{2} N^{2},

where we reintroduced the distribution of the contributing field amplitudes A from each of the N segments.

3.2.2. Enhancement for speckle pulses

In our experiment, the measured second-harmonic signal is not generated from a continuous wave, but from speckle pulses at the particle position. Effectively the wave front shaping optimization affects the original speckle pulse only within a time window of duration of the coherence time τ_c around the point in time at which the optimized pulse is formed. Within this temporal window, the enhancement is given by Eq. (8). At earlier or later moments, the speckle pulse will be modified in a random fashion, but on average it will be unaffected by the optimization.

Since this non-optimized part will always contribute to the time-integrated second-harmonic intensity (i.e. energy), the enhancement of the energy η_pulsed will be lower compared to the hypothetical continuous-wave case. We define the reduction factor c_τ by

c_{τ} \equiv \frac{η_{pulsed}}{η_{cw}} .

In a first approximation, we expect the enhancement to be lowered by the effective number of independent temporal speckle grains within a speckle pulse, which is approximately given by the ratio c_τ ≈ τ_c/τ_d. As the diffuse decay time is proportional to the thickness squared (see section 3.1), c_τ approximately scales with L⁻². However, the ratio τ_c/τ_d does not exactly reflect the temporal distribution of the second harmonic intensity. For a known average intensity of the fundamental light at the particle position 〈I(t)〉 we can calculate the correction factor c_τ more precisely. The average generated second harmonic intensity is consequently proportional to 〈I(t)〉². The reduction factor c_τ for the enhancement compared to the monochromatic case is calculated from the ratio of the energy of the generated second harmonic in a time-window τ_c around the time t_max of the maximum of 〈I(t)〉² and the total second harmonic energy,

c_{τ} = \frac{\int_{t_{max} - τ_{c} / 2}^{t_{max} + τ_{c} / 2} d t {〈 I (t) 〉}^{2}}{\int_{- \infty}^{+ \infty} d t {〈 I (t) 〉}^{2}} .

3.2.3. Susceptibility tensor

In our model we treat the nanocrystal as a single radiating dipole. For a detailed analysis, the light polarization and the second-order susceptibility tensor needs to be considered. The second-order polarization at the second harmonic frequency can be calculated by the matrix equation [15, 23]

P_{2 ω} = [\begin{matrix} 0 & 0 & 0 & 0 & d_{15} & 0 \\ 0 & 0 & 0 & d_{15} & 0 & 0 \\ d_{31} & d_{31} & d_{33} & 0 & 0 & 0 \end{matrix}] (\begin{matrix} E_{c x}^{2} \\ E_{c y}^{2} \\ E_{c z}^{2} \\ 2 E_{c y} E_{c z} \\ 2 E_{c x} E_{c z} \\ 2 E_{c x} E_{c y} \end{matrix}),

where the E_ci are the orthogonal components of the electric field along the three axis in the crystal frame and the d_ij are the second-order susceptibilities of the bulk BaTiO₃ crystal. The considered values are d₁₅ = −41 · 10⁻⁹esu, d₃₁ = −43 · 10⁻⁹esu and d₃₃ = −16 · 10⁻⁹esu [23]. Note that the second-harmonic response is independent of a rotation around the z-axis. The position of the latter in the lab frame is sufficient to describe the second-harmonic response of the nanocrystals, assuming that they are spherical. From Eq. (11) we can see that all components of the vector on the right hand side of Eq. (11) with a non-zero second-harmonic response (

E_{c x}^{2}

,

E_{c y}^{2}

,

E_{c z}^{2}

, 2E_cyE_cz, 2E_cxE_cz) compete for optimization during the wave front shaping optimization. We assume that the transmission coefficients connecting the SLM segments with each of the crystal axis are independent.

For an illustrative purpose we first analytically analyze the case in which the algorithm optimizes all contributions in the E_cx component. Assuming that the detection efficiency is equal for the second harmonic radiation from all crystal axis, and that before optimization the ensemble averaged fields on the three crystal axis 〈|E_ci|〉 are equal, the average generated second harmonic power is proportional to ( $2 d_{31}^{2} + d_{33}^{2} + 4 d_{15}^{2}$ ). Since only the E_cx components, which are generating second harmonic proportional to $α_{31}^{2}$ , are enhanced (with a factor given by the formulas above), the total enhancement is modified by the factor $c_{α} = d_{31}^{2} / (2 d_{31}^{2} + d_{33}^{2} + 4 d_{15}^{2}) \approx 0.17$ . We thereby assume that the optimized component is significantly larger than E_cz after optimization, such that the cross-terms E_cyE_cz and E_cxE_cz can be neglected.

In order to obtain a correction factor close to our experimental situation, we performed numerical simulations of the WFS experiment. For each run of the simulation, we generate a set of random transmission coefficients connecting each of the SLM segments with the three orthogonal field contributions at the crystal position, assuming that the average fields 〈|E_ci|〉, {i = x,y,z} are equal. To calculate the feedback signal, we first apply Eq. (11) to calculate the second-order polarization in the crystal frame. Secondly, the polarization vector is calculated in the lab frame, depending on the orientation of the nanocrystal. Finally, we calculate the second-harmonic intensity as it is collected by a high-NA (NA = 1.4) objective corresponding to our experimental parameters. We apply the sequential optimization algorithm, such as applied in our experiment. As a result, we observe that the algorithm in generally optimizes both the E_cz and the E_cx or E_cy component, with a ratio which varies slightly with crystal orientation. Averaged over all crystal positions, we find that the enhancement of the feedback is modified by the factor c_α = 0.28 ± 0.04 compared to the scalar model (Eq. (8)). Due to the large collection angle, the light radiated from all crystal axis is approximately collected with equal efficiency. The dependence of the factor c_α on the crystal orientation is superseded by variations caused by random variations of the transmission coefficients.

3.2.4. Correction for tight focus

In the introduction of our model we assumed nanocrystals which are much smaller than the focal volume after the optimization. However, the focus is created in a high index medium (n = 2.3) and the medium effectively acts as a lens with a large acceptance angle of about 90°. To a first approximation, the optimization will minimize the focal volume, since it leads to the highest peak intensity and consequently the highest efficiency for the SHG process. Therefore, for a correct description of our experiment we have to consider that the focal volume after WFS will be smaller than the particle volume, whereas before WFS the whole particle volume contributes equally to the average feedback signal. However, the exact shape and polarization state of the tight focus which is formed during a WFS experiment depends on the particle size and shape, the crystal orientation and the specific realization of the photonic environment. Modeling the exact shape of the optimal focus and consequently the generated second harmonic signal in a barium titanate nanocrystal with a size beyond the electrostatic limit is far from trivial, considering that ’simpler’ systems with centrosymmetric materials already require an extensive theoretical treatment [24].

In the following we assess the problem in a simplified treatment. Similar to calculations by van Putten et. al for a linear feedback [25], we calculate the correction factor of the enhancement for second-harmonic generation by

c_{R} \equiv \frac{1}{V} \int_{0}^{R} \int_{0}^{π} \int_{0}^{2 π} \frac{I_{2 ω} (ϕ, θ, r)}{I_{2 ω, peak}} r^{2} sin θ d r d ϕ d θ,

where V is the particle volume of the nanocrystal, and the fraction is the position-dependent intensity I_2ω(ϕ, θ, r) of the second harmonic radiation integrated over the three crystal axis, normalized to its peak intensity I_2ω,peak at the center of the focus. To assess the focal volume quantitatively, we assume that the focus on the nanocrystal has the same profile as a focus created with a high NA lens with and acceptance angle of 90°. From [15] we can conclude that the crystal orientation determines whether linearly or circularly polarized light is more efficient for the SHG process. We assume that the WFS process always converges to the optimal polarization state. We calculate the field distribution of the fundamental radiation at the focus according to [26]. Using the field distribution, we obtain I_2ω(ϕ, θ, r) using the susceptibility tensor given in Eq. (11). We perform this calculation for a sufficient number of polarization states between linearly and circularly polarized light. We calculate c_R, averaged over all orientations of the crystal c-axis. For each angle we thereby consider only the respective polarization state which maximizes the SHG process. We used the parameters of the particle radius R = 100nm, and the refractive index n = 2.3. We find c_R = 0.57 ± 0.05. We observe that c_R is approximately proportional to the radius in the considered size regime. For the crystals used in our experiment, we estimate a size polydispersity of 25%, which will result in equivalent variations of c_R.

Note that the simplification of our approach in Eq. (12) is twofold. Firstly, photonic effects from the spherical shape are neglected. Secondly, the second harmonic intensities within the volume are integrated, and not the electric fields, neglecting interference between different dipole radiating at the second harmonic frequency. We assume that both effects will affect both the numerator of Eq. (12) (optimized focus) and the denominator (reference signal from the average wave field) in the same fashion and therefore tend to level out. Given that the correction factor calculated with our simplified model is rather moderate with 0.57 (compared to the value 1 for a sphere smaller than the focal volume), we do not expect a drastic deviation of the correction factor if these effect would be taken into account.

3.2.5. Noise

The presence of noise on the measurement of the feedback will reduce the observed enhancement. For segments whose contribution to the feedback is on the order or below the noise level, the correct phase value will not be found. Taking this effect into account, we extend Eq. (8) to

η_{cw} \approx 0.5 {(\frac{{〈 | A_{a} t_{a b} | \cdot γ (| A_{a} t_{a b} |, σ_{a}, N_{φ}) 〉}^{2}}{〈 {| A_{a b} t_{a b} |}^{2} 〉})}^{2} N^{2} .

The function γ(|At_ab|, σ, N_φ) is the correlation between the optimal phase ϕ_a = −arg(t_ab) and the phase in the presence of noise ϕ_σ for a given noise level σ_a and the magnitude of the contribution given by |At_ab|.

In the following we describe how γ(|At_ab|, σ_a, N_φ) is calculated. The phase which is obtained from the fit of the feedback signal with a cosine function (see sections 2 and 3.2.1) is equivalent to the phase of the first non-zero frequency component of a discrete N_φ-point Fourier transform of the feedback scan. Component k of the Fourier transform is given by

{FT}_{k, a} \equiv \sum_{n = 1}^{N_{φ}} {I_{2 ω, a} (φ_{n})} e^{- i 2 π \frac{k}{N} n} .

If a single oscillation of a cosine signal with amplitude |At_ab| is sampled with N_φ steps, the amplitude in the first non-zero frequency component of a discrete Fourier transform is

s = \frac{1}{2} N_{φ} | A t_{a b} |

. Gaussian white noise with a standard deviation σ_a results in a mean amplitude in same component of

n = \sqrt{\frac{π}{4}} \sqrt{N_{φ}} σ_{a}

. However, the noise n has a random phase with respect to the signal s, leading to the mentioned deviation from the optimal phase. The probability density function of the phase deviation θ = ϕ_σ − ϕ_a is given by [27]

p (θ; k) = \frac{e^{- \frac{1}{2} k^{2}}}{2 π} + \frac{k cos (θ)}{\sqrt{2 π}} e^{- \frac{1}{2} k^{2} {sin}^{2} θ} Ψ (k cos (θ))

with k = s/n and

Ψ (x) = \frac{1}{\sqrt{2 π}} \int_{- \infty}^{x} d y e^{- \frac{1}{2} y^{2}} .

Consequently the phase correlation can be calculated by

γ (| A t_{a b} |, σ_{a}, N_{φ}) = \int_{0}^{2 π} d θ cos (θ) p (θ; k (| A t_{a b} |, σ_{a}, N_{φ})) .

When both contributions are present, the ensemble averaged squared amplitude of the first non-zero frequency Fourier component is written as

〈 {| {FT}_{1, a} |}^{2} 〉 = \frac{1}{4} N_{φ}^{2} 〈 {| A_{a} t_{a b} |}^{2} 〉 + N_{φ} 〈 σ_{a}^{2} 〉 .

Similarly, the second component squared, 〈|FT_2,a|²〉, will be the sum of the contribution of the noise,

N_{φ} 〈 σ_{a}^{2} 〉

, and, according to Eq. (5), a contribution proportional to 〈|A_at_ab|〉⁴. The latter contribution, since it is amplified less by the other segments (see Eq. (5)), is expected to be about a factor

N_{s}^{2}

lower compared to the first term on the right hand side of Eq. (18) and should therefore be negligible. Experimentally determined values 〈|FT_2,a|²〉 can therefore be used to calculate the noise level, with which 〈|A_at_ab|²〉 can be determined from Eq. (18) using the experimental values 〈|FT_1,a|²〉.

4. Results and discussion

4.1. Spatiotemporal focus on a single nanocrystal

Figure 2(a) depicts the experimentally observed average image of several nanocrystals at the sample backside, which is the average of 200 images, each with a different randomly generated illuminating phase pattern. A number of isolated sources of second harmonic radiation are visible in the field of view. The spot sizes are larger than the nanocrystals themselves, which we explain by scattering from neighboring silica particles.

Fig. 2 Sample backside imaged at the second harmonic wavelength λ = 420nm. (a) Average image obtained from 200 images, each with a different randomly generated illuminating phase pattern. Several nanocrystals are visible in the field of view. (b) Image after wave front shaping. The feedback signal for the optimization was the average count rate in a square of 1μmx1μm around the position of the selected particle, indicated by the dashed square. Here, the feedback signal was enhanced by a factor of η_exp = 3.0 · 10².

Download Full Size | PDF

We performed the WFS experiment for five of the particles in the field of view and a sixth particle at a different position on the sample (not shown). As feedback signal we use the average count rate from a square of 1μmx1μm around the position of the selected particle. The size of the feedback area was chosen to balance between the collection of the largest possible amount of the scattered SHG signal and the increasing influence of camera noise with an increasing feedback area. The feedback signal increases largely during the first three WFS sequences and typically converges to its maximum value after the second sequence with the highest chosen segmentation (N=768).

For all particles the integrated feedback signal was significantly increased (by a factor η_exp = 0.7 · 10² to 5.5 · 10², see Table 1). The image after optimization for one of the particles is shown in Fig. 2(b). We observed that for different nanocrystals the focus of the detection objective had to be moved towards the scattering medium (estimated adjustment in the range of up to a few μm by the manual adjuster) to obtain the smallest spot size on the camera after optimization. For different particles we observed spot sizes (FWHM of a two-dimensional Gaussian fit) from 0.2 μm (close to the diffraction limit) to about 0.6 μm, where bigger spots were typically more of an irregular shape rather than a homogeneous spot. These observations reflect the distribution of the nanocrystals in the silica layer and the resulting scattering of the SHG radiation by surrounding silica particles.

Table 1. Summary of the results of six wave front shaping experiments, each with a different particle visible in Fig. 2. The pulse duration τ_pulse is calculated from the Gaussian fit to the AC curve after optimization. The experimentally observed enhancement η_exp is calculated by the ratio of the average count rate in a 1μmx1μm square around the particle position after and before optimization. For each particle we calculated the enhancement η_model of the time-integrated SH according to our model, calculated by η_model = c_τc_αc_rη_cw. The theoretical enhancement of the SH in a hypothetical continuous wave experiment η_cw, is calculated based on the experimentally determined amplitude contribution. The factor c_τ, which corrects the enhancement for speckle pulses, was calculated using Eq. (10), where the sample thickness to model 〈I(t)〉 is determined from the fits of the AC in Fig. 3(a). The particles 1–5 were all located in close vicinity, where we assume a constant thickness and an average value c_τ = 0.12 was determined. Particle 6 was located on a different spot on the sample, where we obtained c_τ = 0.17. The factors c_α = 0.28 and c_R = 0.57 include the polarization and the susceptibility tensor and the focal volume, respectively (see sections 3.2.3 and 3.2.4).

View Table

The AC measurements from the particle in Fig. 2(b) before and after WFS are shown in Fig. 3. The behavior of the AC before WFS is fitted with the speckle AC model explained above (see Eq. (2)). The AC fit function is calculated based on average time-resolved transmission 〈I_d(t)〉 at the particle position according to [22] using the experimental parameters (mean free path l = 3.5 μm, extrapolation length ratio at silica-air interface z_e₁ = 1.38 and z_e₂ =0.71 at the silica-glass interface, effective refractive index n_eff = 1.25, beam waist of illumination w_I = 150nm and detector size w_D = 100nm). Since the thickness of the sample varies between about 25 – 50 μm, we use the thickness L as fit parameter. After WFS, the AC shows a sharp peak, demonstrating that the light is focused as a short pulse. A contribution of the non-optimized part of the pulse (see section 3.2.2), which could be expected as a small signature next to the correlation peak of the focused part, is not visible due to a present higher noise level. The measured AC can be fitted very well with the AC based on a sech² pulse shape. We calculate the pulse duration from the fit with the well-known deconvolution factor τ_pulse = 0.65τ_AC, τ_AC being the FWHM of the fit. We have summarized the results obtained from all six particles in Table 1. For all pulses, the pulse shape changed from speckle pulse before the optimization to a single short pulse after WFS with a duration τ_pulse ranging between 102fs to 111fs.

Fig. 3 Second harmonic auto-correlation (AC) measured on an individual nanocrystal at the backside of the sample. The graphs show the mean count rate in a 1μmx1μm square around particle position. In particular, the measurements for particle 5 are shown (see Table 1). (a) AC before wave front shaping. The data shows a 3:2:1 contrast ratio between the quickly decaying correlation of the coherent part of the speckle pulse followed by a slower decay of the incoherent part of the speckle pulse towards the background. The data is fitted according to that model of Eq. (2). The only free parameter of the fit is the thickness of the medium. (b) AC measured on nanocrystal at the backside of a slab of disordered silica after wave front shaping. The FWHM of the fit based on the AC of a sech² shaped pulse is 168fs, indicating a pulse duration of the fundamental pulse of 109fs at the particle position.

Download Full Size | PDF

4.2. Comparison of the measured and the modeled enhancement

In the following we compare the measured enhancement with the model introduced above. In order to evaluate Eq. (13) and Eq. (10), we determine the amplitude distribution |A_at_ab| and the average time-resolved transmission 〈I_d(t)〉 from the experimental data.

In order to calculate η_cw from Eq. (13), we determine the amplitude distribution |A_at_ab|, the noise σ_a and the phase correlation γ(|At_ab|, σ_a, N_φ) from the data following the considerations of section 3.2.5. In particular, we analyze the recorded feedback values during the WFS experiment, during a sequence after the algorithm already has converged. We calculate the Fourier transform of the phase scan for each segment by Eq. (14). The calculated components |FT_1,a|²and |FT_2,a|² on the SLM for particle 5 are shown in Fig. 4. The distribution of the components |FT_1,a|² reflects the Gaussian spatial profile of the illuminating beam on the SLM and the random transmission through the sample. The components |FT_2,a|² show a weak dependence on the segment position, with slightly higher average values for the center segments of the SLM. Furthermore, we find the phase values of FT_1,a and FT_2,a to be uncorrelated. Using Eq. (18), we obtain 〈|A_at_ab|²〉 by fitting the values of |FT_1,a|² in Fig. 4(a) with a two-dimensional Gaussian after subtracting the noise $〈 σ_{a}^{2} 〉$ , which itself was determined from a two-dimensional Gaussian fit of |FT_2,a|². In order to calculate γ(|At_ab|, σ_a, N_φ), we numerically evaluate Eq. (17) for each segment using the noise σ_a and the amplitude 〈|A_at_ab|〉 obtained from the fits. We successfully tested our analysis on a simulated WFS experiment (vector model, monochromatic light). In the simulation the enhancement is correctly calculated with an accuracy of about 20% for all noise levels as long as the convergence of the optimization is not spoiled by the noise.

Fig. 4 Fourier analysis of the feedback signal as a function of phase φ_a per segment during a sequence of the wave front shaping experiment. The graph on the left (a) shows the squared amplitude of the first non-zero frequency Fourier component |FT₁|² calculated by Eq. (14), the graph on the right (b) depicts the second component |FT₂|². The distributions are used to determine the amplitude distribution of the contributing segments and the noise level (see sections 3.2.5 and 4.2).

Download Full Size | PDF

It is noteworthy that the components |FT_2,a|² in Fig. 4(b) show a weak dependence on the segment position, which can be caused by several effects. Firstly, cross-talk between pixels of the SLM and diffraction effects for larger phase shifts between neighboring segments can cause a noise term which depends on the intensity present on the segment. Secondly, non-linear effects which are not considered in the above model might be present, such as two-photon absorption in the disordered medium. For all cases, we assume that the linear transmission will be affected by noise or a noise-like contribution of the same magnitude and therefore treat the contributions |FT_2,a|² as noise term to analyze Eq. (18).

In order to evaluate the integral in Eq. (10) the speckle correlation time τ_c is calculated as the half width at half maximum of the measured Fourier transformed laser spectrum. We model the average time-resolved transmission 〈I(t)〉 at the particle position according to [22] as described above, with the thickness L obtained from the fit to the AC measurements (Fig. 3(a)).

With the calculated amplitude enhancement η_cw, and the factors c_τ, c_α and c_R, we calculate the enhancement predicted by our model, η_model = c_τc_αc_rη_cw. For all particles, the calculated values are listed in Table 1. For the particles 1, 3, 5 and 6 the modeled enhancement predicts the experiment value within the accuracy of our model. For particles 2 and 4 it overestimates the enhancement by a factor 4.6 and 4.1 respectively. An overestimation of this magnitude occurs, if two or more particles are placed in direct vicinity, a case which is difficult to identify from the camera images. All particles would contribute to the average reference signal, which would decrease the calculated enhancement accordingly. Furthermore, an overestimation of the enhancement is not surprising considering the composition of the signal which is being optimized. For a given experimental realization of the sample and the nanocrystal in our experimental setup, the feedback signal is solely a function of the phase values set on the SLM. However, next to a global maximum, for which we perform our calculations, the feedback function has a vast number of local maxima into which the algorithm will converge for a given set of starting phase values. A specific example is the factor c_τ which is calculated for the ideal case, that the optimized pulse is formed exactly at the moment in time where the average time-resolved transmission has its maximum. However, in the experiment the point in time at which the optimized pulse is formed will most likely not coincide with the maximum, leading to a lower observed enhancement. Applying different algorithms for the wave front shaping optimization, e.g. a genetic algorithm [28, 29], could be a way to increase the experimentally achieved enhancement. We assume that a similar experiment with lower complexity, e.g. by use of nanoparticles with a less complex second-order susceptibility (see section 3.2.3) or by using quasi-monochromatic light, will produce enhancement values which are predicted even more accurately by our model.

4.3. Pulse duration after WFS

Here we discuss the observed duration of the focused pulses after WFS. The pulse duration, derived from the AC after optimization varies only slightly for different particles, between 102fs to 111fs. Altogether the values are longer than the correlation time τ_c = 52fs calculated from the measured spectrum of the laser. This observation resembles that made by the first study of (far-field) spatiotemporal focusing through a turbid medium [10]. As shown in [10], the optimization process leads to spectral narrowing. The spatial optimization with the SLM, which for each segment induces a phase-shift independent of frequency, cannot optimize all frequencies equally well. The optimization is biased towards those frequencies, which contribute higher to the feedback signal, e.g. in average the frequencies in the center of the spectrum. In the approach of [10], linear time-resolved optical gating is applied to generated the feedback signal, and the resulting spectrum of the focused pulse can be quantitatively predicted exactly from the known laser spectrum.

Following this consideration, we can give an upper limit for the final pulse duration for the present experiment. Here the feedback signal is based on second-harmonic generation, of which the intensity is proportional to the field amplitude at the fundamental frequency to the power of four. Using this dependence as the weighting factor for the transmission of the original spectral amplitude results in a calculated pulse duration of 122fs, assuming that the focused pulses are bandwidth limited. However, this calculation simplifies the fact that spectral amplitude of the second harmonic radiation depends both on the spectral amplitude and spectral phase function of the fundamental radiation. The exact duration of the part of the speckle pulse which is optimized and respectively its spectral composition at the beginning and their dynamics during the optimization procedure depends in a complex way on the starting conditions and the exact transmission coefficients. From the experimental observations we conclude that the spectral function will in general not converge to the calculated limit but remain slightly broader, resulting in shorter pulse duration than the calculated value. Any pulse duration in between the original bandwidth limit and the upper limit can be reached, with the most likely final pulse durations in the observed range.

4.4. Peak-to-background ratio

For future applications of the presented experiment, it is useful to estimate the temporal peak-to-background ratio of the focused pulses. The focused pulses are concentrated in a narrow time window as discussed above, but will still be preceded or followed by a small diffuse speckle contribution. For the generated second-harmonic pulse, the theoretical value for the ratio is given by the calculated enhancement at a single point in time (Eq. (8)). From the experimental enhancement values divided by the effective number of independent speckle grains c_τ, we calculate a peak-to-background ratio between 6 · 10² and 3.2 · 10³ for the investigated particles. The analog peak-to-background ratios for the fundamental intensity are approximately given by the square root of these values.

5. Conclusions

In conclusion, we have demonstrated spatiotemporal focusing inside scattering media using second-harmonic generation in nanocrystals as a feedback signal for wave front shaping. We developed a model which predicts the observed increase of the feedback signal and the pulse duration at the focus well. Our method provides a means to locally generate a short coherent light pulse with a high signal to background ratio inside a complex scattering medium. This provides an exciting tool for various applications of the control of light propagation in complex and especially nonlinear media, and to study light propagation in the presence of strong multiple scattering.

Acknowledgments

We thank Mohamed Tachikirt for the fabrication of the samples, Timmo v.d. Beek and Elbert van Putten for helpful discussions, and Daryl Beggs for careful reading of the manuscript. This work is part of the Industrial Partnership Program (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 and links

1. I. M. Vellekoop and A. P. Mosk, “Focusing coherent light through opaque strongly scattering media,” Opt. Lett. 32(16), 2309–2311 (2007). [CrossRef] [PubMed]

2. T. Cizmar, M. Mazilu, and K. Dholakia, “In situ wavefront correction and its application to micromanipulation,” Nat. Photonics 4, 388–394 (2010). [CrossRef]

3. E. G. van Putten, D. Akbulut, J. Bertolotti, W. L. Vos, A. Lagendijk, and A. P. Mosk, “Scattering lens resolves sub-100 nm structures with visible light,” Phys. Rev. Lett. 106, 193905 (2011). [CrossRef] [PubMed]

4. I. M. Vellekoop, E. G. van Putten, A. Lagendijk, and A. P. Mosk, “Demixing light paths inside disordered metamaterials, ” Opt. Express 16(1), 67–80 (2008). [CrossRef] [PubMed]

5. M. Fink, “Time reversed acoustics,” Phys. Today 50(3), 34–40 (1997). [CrossRef]

6. G. Lerosey, J. de Rosny, A. Tourin, and M. Fink, “Focusing beyond the diffraction limit with far-field time reversal,” Science 315, 1120–1122 (2007). [CrossRef] [PubMed]

7. M. Cui and C. Yang, “Implementation of a digital optical phase conjugation system and its application to study the robustness of turbidity suppression by phase conjugation,” Opt. Express 18(4), 3444–3455 (2010). [CrossRef] [PubMed]

8. C. L. Hsieh, Y. Pu, R. Grange, and D. Psaltis, “Digital phase conjugation of second harmonic radiation emitted by nanoparticles in turbid media,” Opt. Express 18(12), 12283–12290 (2010). [CrossRef] [PubMed]

9. I. M. Vellekoop, M. Cui, and C. Yang, “Digital optical phase conjugation of fluorescence in turbid tissue,” Appl. Phys. Lett. 101, 081108 (2012). [CrossRef]

10. J. Aulbach, B. Gjonaj, P. M. Johnson, A. P. Mosk, and A. Lagendijk, “Control of light transmission through opaque scattering media in space and time,” Phys. Rev. Lett. 106, 103901 (2011). [CrossRef] [PubMed]

11. O. Katz, E. Small, Y. Bromberg, and Y. Silberberg, “Focusing and compression of ultrashort pulses through scattering media,” Nat. Photonics 5, 372–377 (2011). [CrossRef]

12. J. Aulbach, A. Bretagne, M. Fink, M. Tanter, and A. Tourin, “Optimal spatiotemporal focusing through complex scattering media,” Phys. Rev. E 85, 016605 (2012). [CrossRef]

13. F. Helmchen and W. Denk, “Deep tissue two-photon microscopy,” Nat. Meth. 2, 932–940 (2005). [CrossRef]

14. P. C. Ray, “Size and shape dependent second order nonlinear optical properties of nanomaterials and its application in biological and chemical sensing,” Chem. Rev. 110(9), 5332–5365 (2010). [CrossRef] [PubMed]

15. C. L. Hsieh, Y. Pu, R. Grange, and D. Psaltis, “Second harmonic generation from nanocrystals under linearly and circularly polarized excitations,” Opt. Express 18, 11917–11932 (2010). [CrossRef] [PubMed]

16. L. L. Xuan, S. Brasselet, F. Treussart, J.-F. Roch, F. Marquier, D. Chauvat, S. Perruchas, C. Tard, and T. Gacoin, “Balanced homodyne detection of second-harmonic generation from isolated subwavelength emitters,” Appl. Phys. Lett. 89(12), 121118 (2006). [CrossRef]

17. R. Grange, T. Lanvin, C. L. Hsieh, Y. Pu, and D. Psaltis, “Imaging with second-harmonic radiation probes in living tissue,” Biomed. Opt. Express 2(9), 2532–2539 (2011). [CrossRef] [PubMed]

18. C. L. Hsieh, R. Grange, Y. Pu, and D. Psaltis, “Three-dimensional harmonic holographic microcopy using nanoparticles as probes for cell imaging,” Opt. Express 17(4), 2880–2891 (2009). [CrossRef] [PubMed]

19. J. C. Diels, J. J. Fontaine, I. C. McMichael, and F. Simoni, “Control and measurement of ultrashort pulse shapes (in amplitude and phase) with femtosecond accuracy,” Appl. Opt. 24(9), 1270–1282 (1985). [CrossRef] [PubMed]

20. D. J. Thouless, “Maximum metallic resistance in thin wires,” Phys. Rev. Lett. 39, 1167–1169 (1977). [CrossRef]

21. R. Landauer and M. Buttiker, “Diffusive traversal time: effective area in magnetically induced interference,” Phys. Rev. B 36, 6255–6260 (1987). [CrossRef]

22. I. M. Vellekoop, P. Lodahl, and A. Lagendijk, “Determination of the diffusion constant using phase-sensitive measurements,” Phys. Rev. E 71, 056604 (2005). [CrossRef]

23. R. W. Boyd, Nonlinear Optics (Academic Press, 2008).

24. S. Roke and G. Gonella, “Nonlinear light scattering and spectroscopy of particles and droplets in liquids,” Annu. Rev. Phys. Chem. 63, 353–378 (2012) [CrossRef] [PubMed]

25. E. G. van Putten, A. Lagendijk, and A. P. Mosk, “Optimal concentration of light in turbid materials,” J. Opt. Soc. Am. B 28(5), 1200–1203 (2011). [CrossRef]

26. B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems. II. structure of the image field in an aplanatic system,” Proc. R. Soc. Lond. A 253, 358–379 (1959). [CrossRef]

27. J. W. Goodman, Statistical Optics (Wiley, 2000).

28. I. M. Vellekoop and A.P. Mosk, “Phase control algorithms for focusing light through turbid media,” Opt. Commun. 281(11), 3071–3080 (2008). [CrossRef]

29. D. B. Conkey, A. N. Brown, A. M. Caravaca-Aguirre, and R. Piestun, “Genetic algorithm optimization for focusing through turbid media in noisy environments,” Opt. Express 20(5), 4840–4849 (2012). [CrossRef] [PubMed]

Part. #	τ_pulse(fs)	η_exp	η_model	η_cw
	experiment
1	110	2.5 · 10²	3.1 · 10²	1.6 · 10⁴
2	102	0.7 · 10²	3.2 · 10²	1.7 · 10⁴
3	111	2.7 · 10²	2.7 · 10²	1.4 · 10⁴
4	104	0.7 · 10²	2.9 · 10²	1.5 · 10⁴
5	109	3.0 · 10²	3.8 · 10²	2.0 · 10⁴
6	109	5.5 · 10²	6.5 · 10²	2.4 · 10⁴

Spatiotemporal focusing in opaque scattering media by wave front shaping with nonlinear feedback

Abstract

1. Introduction

2. Experiment

3. Theory

3.1. Intensity-intensity auto-correlation of speckle pulses

3.2. Enhancement of the time-integrated second harmonic intensity

3.2.1. Enhancement vs. number of segments

3.2.2. Enhancement for speckle pulses

3.2.3. Susceptibility tensor

3.2.4. Correction for tight focus

3.2.5. Noise

4. Results and discussion

4.1. Spatiotemporal focus on a single nanocrystal

4.2. Comparison of the measured and the modeled enhancement

4.3. Pulse duration after WFS

4.4. Peak-to-background ratio

5. Conclusions

Acknowledgments

References and links

Cited By

Figures (4)

Tables (1)

Equations (18)

Optics Express