## Abstract

Focusing light into highly disordered biological tissue is a major challenge in optical microscopy and biomedical imaging due to scattering. However, correlations in the scattering matrix, known as “memory effects”, can be used to improve imaging capabilities. Here we discuss theoretically and numerically the possibility to achieve three-dimensional ultrashort laser focusing and scanning inside forward scattering media, beyond the scattering mean free path, by simultaneously taking advantage of the angular and the chromato-axial memory effects. The numerical model is presented in details, is validated within the state of the art theoretical and experimental framework and is finally used to propose a scheme for focusing ultra-short laser pulses in depth through forward scattering media.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Exploiting the optical properties of biological tissues [1] is of crucial interest for many biomedical applications such as optical imaging. One critical challenge concerns focusing broadband and ultrashort laser pulses into such complex media, especially for non-linear microscopy [2–5]. However, light beam manipulation is hampered by inhomogeneities of the refractive index map. Coherent light propagation through tissues results in random speckle patterns [6], limiting the penetration depth of optical imaging techniques [7]. Light perturbation by refractive index mismatch has historically been categorized into geometrical aberrations and scattering. Geometrical aberrations typically includes phase-only perturbations described by low-order Zernike modes, such as introduced by the tissue surface roughness. Scattering refers to phase and intensity perturbations of higher spatial frequencies, such as introduced by micron-sized organelles inside the tissue volume. Several approaches have been proposed to increase the light penetration depth by improving the amount of focused ballistic photons by using longer wavelength lasers [8], higher energy pulses [9] and also by optimizing collection of fluorescence light [10]. Adaptive optics has proved to restore imaging capabilities, at least partially, by compensating geometrical aberrations, especially in the case of non-linear optical microscopy [11–13]. All these techniques are focusing on optimizing ballistic light. However, beyond the scattering mean free path ($l_s\simeq 100\mu m$ in tissues), most of the energy is scattered, thus limiting the efficiency of these methods, all the more in the presence of screw phase dislocations involving $0-2\pi$ phase-steps which cannot be properly addressed by deformable mirrors [14]. Recent advances in the field of wavefront shaping [15,16] have overcome this problem by using segmented phase actuators to precompensate the input wavefront. Such an optimization can be achieved by using an iterative process [16], phase conjugation [17], or by measuring the transmission matrix in advance [3,18,19]. These techniques have proven to be powerful tools for optical imaging but they either decrease imaging speed or require access to the plane of interest for characterizing the scattering medium.

Noteworthy, non-invasive [20–23] and high-speed imaging [24,25] could be achieved without any active wavefront control by making use of correlation properties of speckles behind diffusers. These correlations of the transmission matrix have been known as the “memory effect” and allow deterministic transformations of speckles behind diffusers such as transverse [26–29] and axial [30–32] translations as well as topological transforms of critical points [33,34]. The spatial correlations of speckle patterns have been widely studied especially in the case of thin diffusers [26,27] and forward scattering media [28,29,32] (See Fig. 1(a)). In this latter case, the angular range of the ME scales as $\lambda /(L\sin \Theta _0)$, where L is the diffuser thickness and $\Theta _0$ the mean output angular spread under plane wave illumination [28,32,35]. In tissues, scattering mostly occurs in the forward direction, which is quantified by the anisotropy factor $g=\left \langle\cos \theta \right \rangle $, where $\theta$ is the scattering angle at each scattering event. Consequently, energy propagates until the transport mean free path $\ell ^\ast =\ell _s/(1-g)$, which is typically one order of magnitude larger than $\ell _s$ [1]. For $L\ll \ell ^\ast$, the mean angular spread after a slab of such material is then $\Theta _0^2=L/\ell ^\ast$ [28] and the rotation center of the angular memory effect lies at axial coordinate $z=L/2$, inside the slab (Fig. 1(a)). Noteworthy, for a material made of non-dispersive and smooth refractive index map, $\ell ^\ast$ is achromatic and $\ell _s$ scales as $\lambda ^2$ [36], in contrast with the well-known $\lambda ^4$-scaling of the Rayleigh scattering regime. Numerical fitting of experimental data can be obtained by using a combination of these two regimes [1].

Nowadays, much efforts focus on improving imaging capabilities in the intermediate depth range where light energy is propagating forward but with negligible amount of ballistic light ($\ell _s\ll L \ll \ell ^\ast$). In this regime, solutions have been explored to increase the spatial memory effect, including combined angular and translational transforms [28] and time-gated detection [37]. In addition, a strong spatio-spectral coupling has recently been observed [31] wherein a spectral detuning over $~200~\textrm {nm}$-widths was experimentally demonstrated to result in a simple axial translation of the focused laser beam through a $1~\textrm {mm}$-thick brain slice. This effect has been shown to be of interest for 3D imaging with thin diffusers [30,38]. The chromato-axial memory effect ($\chi$-axial ME) was further theoretically investigated and was shown to be especially important in the case of forward scattering samples, with convincing confrontation with experiment [32]: Under plane wave illumination, a spectral shift results in an axial homothetic dilation of the output speckle pattern, with origin located at a virtual plane located at axial coordinate $z=2L/3$ inside the diffuser where the speckle is achromatic both in phase and amplitude, at second order approximation (Fig. 1(b)). The spectral correlation width of the medium is then $\Delta k=3\sqrt {2}/(L\Theta _0^2)$, half-width-at-half-maximum (HWHM). The model developed in Ref. [32] is based on a two dimensional random walk in the k-space and thus relies on the hypothesis of a large number of scattering events, which in practice, accounts for physical samples exhibiting typical phase-correlation dimensions much smaller than the slab thickness. Because studying optical properties of biological tissues is so important, a numerical model accounting for memory effects in thick scattering media is of high interest for achieving 3D imaging in the biomedical imaging domain.

In this article, we first present our numerical tool consisting in modeling a thick forward scattering diffuser by a stack of thin diffusers [39,40], in the intermediate range $\ell _s \ll L\ll \ell ^\ast$. Then this model is validated quantitatively based on theoretical expectations published in the literature, both regarding the angular memory effect [28,32] and the $\chi$-axial ME [32]. In Sec. 2., the model parameters are expressed as a function of the three main parameters characterizing a slab of tissue: L, $\ell _s$, $\ell ^\ast$ (or equivalently: $g$). Then, in Sec. 3., this model is first validated in the case of a very large anisotropy factor ($g=0.999$) which provides a correspondingly very large thickness validity range from $\ell _s$ to $\ell ^*$. In Sec. 4., the model is also validated for an anisotropy factor comparable with typical biological tissues ($g=0.975$). Finally, making use of both $\chi$-axial and angular ME simultaneously, in Sec. 5 this study proposes a strategy for three-dimensional ultra-short beam focusing and scanning inside a scattering tissue beyond the scattering mean free path, from a single monochromatic characterization of the medium.

## 2. Description of the numerical model

Our numerical model is similar to those published in the literature [39,40] and the main programs can be found online in Code 1 [41]. It consists in slicing a forward scattering slab of thickness $L$ into a stack of $N_\textrm {diff}$ thin diffusers separated by distance $d$ as illustrated in Fig. 2(a). Backscattering is neglected because the index contrasts are assumed to be low enough and because the slab thickness $L$ is assumed to be much smaller than the transport mean free path $\ell ^*$. Light propagation through the stack is then achieved in the scalar approximation model and simply consists in the succession of thin phase mask and Fresnel propagations. Thin diffusers are generated by spatial filtering of a zero-meaned pixelated random thickness map thanks to a Gaussian filter with spatial autocorrelation width $w$ and yielding a mean squared optical path delay $\overline {\delta }=\sqrt {\left \langle\delta ^2\right \rangle -\left \langle\delta \right \rangle ^2}$ [40,42] (Fig. 2(b)). The spatial correlation function of the thin diffusers is then $C(\textbf {r}_{\textbf{1}},\textbf {r}_{\textbf{2}})=\exp \left [-(\textbf {r}_{\textbf{1}}-\textbf {r}_{\textbf{2}})^2/(4w^2)\right ]$. Noteworthy, the exact spatial spectrum profile of the diffusers does not matter because the validity of the model relies on the approximation of large number of scattering events, for which the central-limit theorem is valid. The spatial spectrum must only be bound to low scattering angles for Fresnel approximation to be valid. Here, we further assume that the scattering medium is made of non-dispersive materials. The phase delay introduced by the $n^{th}$ thin diffuser is thus assumed to result in a product by the following transmission coefficient:

where $\varphi _n(x,y)=k\delta (x,y)$ is the phase delay introduced by diffuser number $n$, and $\delta (x,y)$ the corresponding optical path delay.The characteristics of the scattering medium $L$, $\ell _s$ and $\ell ^\ast$ can then be tailored solely by tuning the geometrical parameters: the number of thin diffuser plates $N_\textrm {diff}$, their separation distance $d$, mean squared amplitude of optical path delays $\overline {\delta }$ and spatial correlation width $w$. We thus dispose of one more geometrical control parameter than the number of physical constraints, so enabling to set one parameter arbitrarily. The relations between all these parameters is detailed below. They can be obtained by considering first the thickness of the thick diffuser $L = (N_\textrm {diff}-1)d$, second the fraction of ballistic light after the diffuser and third the output scattering angle.

The fraction of ballistic light after a single thin phase mask $n$ may be calculated by considering that $\varphi _n$, for a single thin diffuser, is a zero-mean Gaussian random variable with spatial correlation $w$ and variance $\left \langle\varphi _n^2\right \rangle =\left (k\overline {\delta }\right )^2$. Considering the central moments of a normal distribution, a Taylor expansion of $e^{i\varphi _n}$ provides the average intensity of the ballistic light as:

A numerical illustration of this expression is shown in Fig. 3(a) as a function of $\overline {\delta }$ ($\lambda =0.8~\mu\textrm{m}$ and $w=3~\mu\textrm{m}$). After crossing $N_\textrm {diff}$ thin diffusers, the fraction of ballistic light is then $I_{bal}(L) =\exp \left [-N_\textrm {diff}\left (k\overline {\delta }\right )^2\right ]$, which can be identified to the definition of the scattering mean free path: $I_{bal}(L)=\exp (-L/\ell _s)$, so giving $\ell _s$ as a function of numerical and geometrical parameters.

Finally, the transport mean free path $\ell ^\ast$ is related to the output angular spread of an impinging plane wave. Considering the Gaussian spatial correlation function of thin diffusers, a second order Taylor expansion for small $k\overline {\delta }$ values allows getting an analytical expression of the angular variance of the intensity per transverse dimension:

This expression may then be identified to $\Theta _0^2 = L/\ell ^\ast$ [28]. The validity of this expression is illustrated in Fig. 3(b) for a single thin phase plate. A slight discrepancy is observed for $\overline {\delta }$ values approaching $\lambda$ due to the limited validity of the Taylor expansion: for too large $\overline {\delta }$-values, the power spectrum of the field deviates from a Gaussian-shape.

In a nutshell, the equations relating the geometrical parameters of the model to the physical parameters of the diffuser are:

Since the number of parameters in the model is $4$ and the number of physical parameters to set is only $3$, one parameter can be set arbitrarily. Noteworthy, a thick scattering sample characterized only by the three parameters $\ell _s$, $\ell ^*$ (or $g$) and $L$ (and satisfying $\ell _s \ll L\ll \ell ^\ast$) can be modeled by a stack of any number of thin diffusers, so long as $N_\textrm {diff}$ is large enough for the central limit theorem to be considered as valid [32]. Regarding this constraint, we quantitatively discuss the accurateness of results as a function of $N_\textrm {diff}$ in Sec. 3. Taking $N_\textrm {diff}>5$ is typically sufficient but for more accurate matching with the model, we took $N_\textrm {diff}$ larger than, or equal to, $10$ in our simulations. Furthermore, it is also convenient to take $N_\textrm {diff}$ large enough so that the mean squared phase delay $\sqrt {\left \langle\varphi _n^2\right \rangle }=k\overline {\delta }$ remains smaller than $2\pi$. This second requirement makes former Taylor expansions valid and conveniently allows making use of analytical expression given in Eqs. (3) and (6) to determine geometrical parameters. Interestingly, we point out that that in the framework of this model, the scattering mean free path $\ell _s$ scales as $\lambda ^2$ (Eq. (5)) and the transport mean free path $\ell ^\ast$ is achromatic (Eq. (6)).

## 3. Angular and $\chi$-axial memory effects: numerical *vs.* theoretical results

#### 3.1 Angular memory effect

As a first step, we study the angular memory effect in a thick scattering slab. According to this phenomenon a tilt of the incident wavefront results in a tilt of the output wavefront within a limited angular range [28,32] (See Fig. 1(a)). Since the rotation coordinate of the angular memory effect is at $z=L/2$, the correlation products between output fields are thus considered in this plane for varying input tilt angles. The analytical expression of the cross-correlation product between speckle intensity patterns when tilting the angle of incidence by $\Theta$ is:

The angular correlation width – or angular ME range – full width at half maximum (FWHM), is then:

In Fig. 4(a) the results of simulations is compared to the analytical expression given in Eq. (7), for a diffuser with $\ell _s=100~\mu\textrm{m}$, $L=10\ell _s= 1~\textrm {m m}$, $\ell ^\ast =1000\ell _s= 100~\textrm {m m}$. The angular correlation width (Eq. (8)) is plotted as the function of $N_\textrm {diff}$ in Fig. 4(b) and exhibits quick convergence for a few thin diffuser plates. It is then plotted in Fig. 4(c) for $N_\textrm {diff}=10$ as the function of thickness which varies from $L=5\ell _s$ to $L=20\ell _s$ and is also compared to the analytical expression given in Eq. (8).

It is interesting to relate this model to results published in Ref. [39] wherein the authors compare their experimental results obtained for forward scattering samples, with the angular memory effect derived for media thicker than $\ell ^\ast$, in which case $\Delta \Theta \approx 3.0 \, (kL)^{-1}$ [27,43]. Fitting experimental data with this expression for forward scattering samples, the authors define a effective thickness $L_{eff}$ which turns out to be much thinner than the actual thickness of samples. This smaller effective thickness is due to the larger angular memory effect of forward scattering slabs than would be obtained for isotropically scattering samples of same thickness. In the frame of our model, the effective thickness can be identified to $L_{eff}\approx 0.52 L^{3/2}\ell ^{\ast -1/2}$, in qualitative agreement with experimental data published in Ref. [39].

#### 3.2 $\chi$-axial Me

Another important correlation property of forward scattering media is a $\chi$-axial ME, observed and studied in Refs. [31,32]. We now aim at validating this effect with our multi phase-plates model. In a nutshell, under plane wave illumination condition, a chromatic shift of the beam results in an axial dilation of the speckle pattern over large spectral bandwidths. The analytic expression quantifying this effect could be derived in Ref. [32] thanks to a model based on a two-dimensional random walk in the $k$-space. The two-wavelength mutual coherence function at two different axial positions $z$ and $z^\prime$ for two copropagating impinging plane waves with wavenumbers $k$ and $k'$, respectively, was calculated to be [32]:

The correlation product $C_{\chi }$ in Eq. (9) is then a double Lorentzian function both along spectral detuning parameter $k_0$ and along the axial z-coordinate. As pointed out in Ref. [32], along the propagation z-coordinate, the two-wavelength correlation function $C_{\chi }$ in Eq. (9) is maximum for:

A spectral shift of the incident beam thus results in an axial homothetic dilation of the speckle with origin at $z=2L/3$, as illustrated in Fig. 1(b). The Lorentzian evolution of $C_\chi$ (Eq. (9)) as a function of $z=z^\prime$ is plotted in Fig. 5(a) and compared with numerical data. In this example, the zero-mean cross-correlation product of two speckles ($\lambda =800~\textrm {nm}$ and $\lambda ^\prime =820~\textrm {nm}$) is computed for a scattering slab of thickness $L=1~\textrm {mm}$. At plane $z=z^\prime =2L/3$, the speckles are achromatic at second order approximation in $k_0$. The location of the achromatic plane was computed numerically from thicknesses ranging from $L=\ell _s$ to $L=20\ell _s$ and compared to the $2L/3$ value (Fig. 5(c)).

At the achromatic plane ($z=z^\prime =2L/3$), the zero-mean cross-correlation product is given by:

yielding the spectral correlation width (FWHM) at first order in $\Delta k=k-k^\prime$: or equivalently for $\Delta \lambda$:The analytical Lorentzian profile as a function of $\lambda$ (given in Eq. (12)) is compared with the numerical computation of the spectral correlation width in Fig. 6(a). For $L=1~\textrm {mm}$, numerical data exhibit excellent agreement with theory. However, for $L=0.5~\textrm {mm}$, a significant discrepancy appears for longer wavelengths. This discrepancy can be understood by plotting the spectral dependence of the ballistic intensity $I_{bal}=exp(-z/l_s)$ in Fig. 6(b). Since the scattering mean free path $\ell _s$ scales as $\lambda ^2$ in the frame of this model (Eq. (5)), the approximation $L\gg \ell _s$ becomes wrong for longer wavelengths: the model prerequisites are thus not satisfied for longer wavelengths. In this case $L\sim l_s$, not only some ballistic light remains but also, the scattered light does not exhibit Gaussian statistics. The observed remaining correlations at long wavelengths is due to both of these consequences. Equations (13) and (14) are finally compared to numerical results in Fig. 7(b) and 7(c) showing again good agreement for thick enough scattering samples. The spectral correlation width is shown to scale as $1/L^2$, in agreement with Eqs. (13) and (14).

#### 3.3 Combined angular and $\chi$-axial MEs

In a final step, we study the simultaneous effect of both angular and $\chi$-axial ME in the perspective to perform three-dimensional intensity scanning inside or through forward scattering diffusers. We thus considered numerically a simultaneous angular tilt and chromatic shift of the impinging beam. Numerical simulations were achieved by computing speckle intensities in the achromatic plane ($z=2L/3$). For a $\Theta$ angular tilt, speckles were checked to be translated by an amount $\Theta L/6$ since the tilt rotation origin is at $z=L/2$. Making use of Eqs. (7) and (9), the combined angular-spectral memory effect is shown, in Fig. 8 to behave as:

## 4. $\chi$-Axial and Angular ME for a Realistic Anisotropy Factor

So far, we used a large anisotropy factor ($g=0.999$) in order to check our model over a large range of slab thickness from $\ell _s$ to $\ell ^*$. In particular, we could demonstrate that in this regime, spectral correlation widths are so large that the amount of ballistic light becomes significant within the spectral correlation range of the medium. However, such g-values are not really typical for biological tissues. We thus now illustrate the validity of our model for $g=0.975$ in Fig. 9.

## 5. Spectral cophasing plane & pulse broadening

In section 3.2, the speckles at different wavelengths were shown to exhibit a correlation maximum in the virtual plane $z=2L/3$, meaning that if sending a broadband light beam, a “white” speckle is obtained in this plane. However this result does not presume anything about the spectral phase of the beam or about ultra-short pulse broadening. Pulse broadening in scattering media reduces non-linear optical processes, especially those based on electronic virtual quantum states (two-photon absorption, second harmonic generation $\cdots$). Although usually considered as detrimental for imaging, Oron *et al.* considered turning the pulse-broadening effect to the experimentalist advantage for achieving spatio-temporal focusing in non-linear microscopy applications [44–47].

Here, we thus derive the analytical expression of the pulse beam broadening after the scattering sample, based on the analytical expression of the spectro-axial correlation properties of forward scattering diffusers (Eq. (9)). The temporal profile of a pulse can be obtained by considering the corresponding broadband electric field:

yielding, after the diffuser, the time-dependent intensity (averaged over transverse spatial coordinates):Simple integral manipulation leads to:

Eq. (21) is plotted in Fig. 10 together with the result of numerical simulations and demonstrate a very good agreement for a $L=500~\mu \;\textrm {m}$ scattering sample whose scattering properties are consistent with biological tissues’ ($l_s=100~\mu \;\textrm {m}$, $g=0.975$). Interestingly, the minimum pulse width is obtained in the virtual plane $z=2L/3$ wherein speckles are achromatic. Away from this plane, the difference in homothetic axial dilation for the different spectral components leads to a linear increase of the pulse width. As a final remark, Eq. (21) has been derived for a impinging plane wave but exhibit a similar form as the expression calculated for a focused beam [47].

## 6. Broadband spot imaging and focusing

Until now, we only considered collimated beam illumination of the diffuser at normal incidence. Noteworthy, the former conclusions are also trivially valid for tilted plane wave illumination. In this section, we study the implications of the $\chi$-axial ME in the case of spherical wave illumination. This problem is of interest for a wide range of applications like broadband (*e.g.* fluorescent) point source imaging and broadband beam focusing through or inside scattering media. Imaging broadband point source through a scattering medium demands special care for collecting light [10]. Furthermore, recent blind speckle imaging [20–22] and focusing [23,48] approaches typically require maximizing the speckle contrast because of the limited dynamic of detectors.

The different cases studied in this sections are illustrated in Fig. 11. First, in Sec. 6.1, we study the speckle-wavefield generated in the achromatic plane $z^\prime =-2L/3$ by a single broadband point source located behind a forward scattering medium (Fig. 11). Then, in Sec. 6.2, based on field conjugation of the wavefield generated by a monochromatic point source (located at axial coordinate $z$), we consider scanning this spot axially by two different means: by shifting the wavelength of the impinging beam (Fig. 11(b)), making use of the $\chi$-axial ME, or by imprinting a parabolic wavefront modulation to the impinging beam (Fig. 11(c)), so taking advantage of the angular ME of the medium. Finally, in Sec. 6.3, a strategy is proposed to achieve ultra-short, broadband, pulse focusing based on the wavefield generated by a monochromatic point source (Fig. 11(d)): dispersive parabolic wavefront modulation is applied to a broadband light beam to compensate the $\chi$-axial ME.

#### 6.1 Broadband point source imaging

We study here the wavefield generated by a broadband point source located behind a diffuser (Fig. 11(a)). We assume a Gaussian point source with numerical aperture $\textrm {NA}=0.1$ to remain in the validity domain of small angles approximation. Although the angular memory effect is a factor five below the NA-value ($\Delta \Theta = 2.10^{-2}$), the speckle intensity is achromatic in the plane $z^\prime = -2L/3$ as illustrated in Fig. 12(a) where two speckles are shown for wavelengths exhibiting a spectral shift of $220~\textrm {nm}$ (smaller than the spectral correlation width of the diffuser $\Delta \lambda = 345~\textrm {nm}$, FWHM). The phase difference $\Delta \varphi$ between these two random waves is shown in Fig. 12(b) where it appears as mostly parabolic. This result may be understood by considering that pencil beams arising from the point source are mostly scattered in the forward direction. The resulting phase delay is thus dominated by a spherical component (as shown in Fig. 12(d)): $\varphi = k\textbf {r}^2/(2R)$, with $R = z+2L/3$. As a result, the different wavelengths of the outgoing beam share a same wavefront, resulting in parabolic relative phase difference between two different wavelengths: $\Delta \varphi = k\textbf {r}^2/(2R_{\delta \lambda })$, where $R_{\delta \lambda }$ may be expressed at first order approximation as:

The difference between this model (Fig. 12(c)) and numerical results (Fig. 12(b)) is shown in Fig. 12(e) where the phase difference appears as mostly flat.

In a nutshell, when considering plane wave illumination, the speckle through the diffuser is achromatic both in phase and amplitude in the virtual plane located at $2L/3$ inside the diffuser. In contrast, for point source illumination, the speckle is still achromatic in amplitude in this plane but the achromatic spherical wavefront results in a relative wavelength-dependent parabolic term.

#### 6.2 Axial spot scanning by chromatic shift *vs.* parabolic wavefront shaping

In the former section, we observed that the $\chi$-axial ME for a spherical wave results in a parabolic relative-phase-delay between two wavelengths, in the plane $z^\prime =-2L/3$. This result strongly recall the possibility to axially shift a spot through a thin diffuser by exploiting the angular memory effect by simply adding a parabolic phase delay.

Here, we thus compare the two approaches to axially scan a spot focused behind a forward scattering diffuser. The field originating from a point source with $\textrm {NA}=0.1$ was then computed in the virtual plane located at $z^\prime = -2L/3$, as represented in Fig. 11(a). Refocusing through the diffuser was then obtained by performing field conjugation. Prior to propagation through the diffuser, the field was modified in two different ways: propagation was achieved either for a shifted wavelength (scheme proposed in Fig. 11(b)), exploiting the $\chi$-axial ME, or after impinging a parabolic phase modulation (scheme illustrated in Fig. 11(c)), exploiting the angular ME. In both cases, these changes result in a axial shift of the focused spot. For a chromatic shift $\delta k$, the axial shift of the spot is

A same amount of defocus can also be obtained by imprinting a parabolic phase modulation:

The resulting focused spot intensities are plotted in Fig. 13 for three values of the initial focus location $z$: $0$, $L$, and $2L$; showing that the two approaches are almost perfectly equivalent. This result suggest that the ability to shift axially a beam arise from a single intrinsic invariant property of the scattering medium.

#### 6.3 Broadband beam focusing: parabolic wavefront compensation of the chromatic shift

Because of the $\chi$-axial ME, achromatic spatial phase-shaping results in an axially spread focus [31] (as illustrated in Fig. 11(b)). For a spot located at a distance $z$ behind the diffuser, the axial shift resulting from a spectral shift $\delta k$ is given by Eq. (23). In order to superimposes focii of all the spectral components, an achromatic wavefront – corresponding to a wavelength-dependent parabolic phase modulation – must then be applied (opposite of Eq. (24)) to compensate for this axial shift (Fig. 11(c)).

In practical implementations, achromatic parabolic wavefront modulation can be obtained thanks to an achromatic lens [31,49]. In contrast, phase only spatial light modulators introduce little dispersion so inducing almost achromatic phase modulation. Moreover dispersive spectral defocus can be simply achieved by using regular lenses made of dispersive materials. In our numerical simulations, we simulated the pulse refocused behind the diffuser when achieving (quasi-)phase conjugation by only considering the wavefield required for focusing at the central wavelength $\lambda _0$ of the illumination source. This required wavefield $E_{\lambda _0}$ is obtained by simulating a “guide star” located at $z=0$, at the central wavelength $\lambda _0$ only, in the virtual plane $z^\prime =-2L/3$. Like in previous sections, the considered “guide star” was actually a Gaussian spot with numerical aperture $NA=0.1$. Then the quasi-conjugate impinging beam defined as

was propagated back through the diffuser. Such a beam could be generated by combining a spatial light modulator achieving achromatic phase conjugation together with a focusing lens achieving achromatic spherical wavefront modulation. As a result, the spectrum at the focus is observed to be slightly filtered (blue curve in Fig. 14(a)) as compared to the initial spectrum (dashed black curve in Fig. 14(a)) due to imperfect phase conjugation. However, such*a priori*chromatic defocus yields a much broader spectrum than when performing achromatic phase conjugation only (pink curve in Fig. 14(a)) – as would be obtained using a spatial light modulator only – in which case the focii at different wavelengths are axially spread. In addition, it can be seen in Fig. 14(a) that not only the spectrum is wider in the case of

*a priori*chromatic defocusing phase conjugation but also that the spectral phase is much flatter. As a consequence, the pulse width at focus is three times as short as compared to achromatic phase conjugation (Fig. 14(c)). The peak intensity was measured to be decreased by $70\%$ when achieving achromatic phase modulation and by only $11\%$ when combining achromatic phase modulation with dispersive refocusing. In Fig. 14(c) most of the energy is focused but interestingly, away from the focus, the pulse width of residual unfocused energy exhibits a temporal broadening similar to so called X-waves studied in Ref. [46].

## 7. Conclusion and perspectives

In this paper, we have demonstrated that the memory effects induced by angular tilt, defocusing and spectral shift can be used to focus an ultrashort laser beam with achromatic wavefront shaping in order to achieve a three-dimensional scanning of the beam through the scattering tissue. First, we proposed a numerical model for a forward-scattering anisotropic medium based on a stack of $N_\textrm {diff}$ thin diffusers as shown in Code 1 [41]. Our diffuser model assumes that the material is fully characterized by three parameters: $\ell _s$, $\ell ^\ast$ and $L$. Furthermore, our model is based on the approximation $\ell _s \ll L\ll \ell ^\ast$. It thus assumes a large number of scattering events, and neglects backscattering. Noteworthy, for non-linear microscopy applications, wherein a non-linear optical response is induced by an ultrashort focused spot, backscattering can be neglected even if the medium is semi-infinite. In this case, the length $L$ is to be interpreted as the considered depth of focus.

Broadband beam focusing through forward scattering media results in an axial dilation of spectral components, in agreement with the $\chi$-axial ME described in Ref. [32]. The origin of the homothety, where the speckle is "achromatic" lies in a virtual plane located at a distance $L/3$ before the considered focal plane. Here we further demonstrate numerically that the spectral phase of the beam is flat in this same plane, resulting in a minimum pulse duration. For broadband light beams with achromatic impinging wavefronts the achromatic (and spectrally flat) plane does not coincide with the focal plane, which is a problem for microscopy applications. The achromatic plane is indeed inaccessible since in a virtual plane.

Here, we thus suggest a strategy to move forward this achromatic (and spectrally flat) plane to the focal plane. The possibility to get the achromatic plane and the focal plane to coincide in the medium is provided by the ability to achieve an axial displacement of the beam either by spectral detuning or parabolic modulation of the impinging wavefront. As a consequence *a priori* quasi phase-conjugation can be achieved based on the field at a single wavelength only. For a slightly detuned beam, a simple additional wavefront curvature $R_{\delta \lambda }=(z+2L/3) \lambda _0/ \delta \lambda$ at the virtual plane located at $L/3$ from the diffuser input surface, allows almost perfect refocusing behind (or inside) the diffuser. The key role of the virtual plane $L/3$ from the diffuser input surface is in perfect agreement with the result derived in Ref. [28] wherein wavefront manipulation for optimized transverse-scanning-range is achieved in this same virtual plane. Here, we get that this plane is also the best to compensate the axial chromatic drift of the focused spot by parabolic wavefront modulation. This result is of high interest for non-linear microscopy applications requiring focusing multiple wavelengths, in which case using a single spatial light modulator is much easier. Comparison between this approach and using the broadband transmission matrix would also be of high interest [49].

Interestingly, for a non-dispersive and non-absorbing diffuser, our model results in an achromatic transport mean free path. The scattering mean free path is also shown to scale as $\lambda ^2$. This behavior results from the model we considered for intermediate thin diffusers which are all assumed to be non-dispersive, smooth and continuous, so yielding phase function mostly pointing towards forward directions. In practice, typical phase biological objects exhibit sharp boundaries mitigating the validity of this approximation: a small fraction of light energy is thus scattered at large angles. Even though the contribution of this scattered light is negligible in the forward direction, it results in a energy-loss and in reduced memory effects. One solution to complete tissue modeling may consist in adding a Rayleigh scattering case, in which case scattering scales as $\lambda ^4$ [1]. Our results showed that the spectral correlation width scales at $\ell ^\ast /L^2$ and that the angular memory effect range scales as $\ell ^{\ast \, 1/2}/(kL^{3/2})$.

Finally, we presented results for non-dispersive and non-absorbing media for the sake of simplicity of physical interpretation but both dispersion and absorption can be trivially introduced into the presented numerical model as we show in Code 1 [41]. Achromatic bulk absorption in the medium is mostly expected to just result in an energy decrease at the diffuser output but should also attenuate the contribution of the longest optical path lengths which are both responsible for the largest outgoing angles and for dispersion in the optical-path-lengths distribution. Achromatic bulk absorption is thus expected to increase the forward-scattering nature of the medium and should thus increase its spectral correlation width. The presence of localized achromatic absorbers should also increase the achromatic nature of the medium as compared to phase masks. Conversely, chromatic-dependent absorption will result in reduced correlation widths both in the bulk or for localized absorbers. Dispersion is not expected to change significantly the random-walk nature of light propagation. The analytical formula derived in Ref. [32] have been derived whatever the refractive index. Numerical results (not shown) confirm this expectation and introducing either bulk or localized dispersion mostly results in a mere wavelength rescaling of plots.

## Funding

Agence Nationale de la Recherche (ANR-18-CE42-0008-01); The French-Israeli Laboratory ImagiNano; Centre National de la Recherche Scientifique; European Research Council (H2020, SMARTIES-724473); H2020 European Research Council Research and Innovation Program (677909).

## Acknowledgement

The authors thank Sophie Brasselet for stimulating discussions.

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. **58**(11), R37–R61 (2013). [CrossRef]

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

**3. **H. B. de Aguiar, S. Gigan, and S. Brasselet, “Enhanced nonlinear imaging through scattering media using transmission-matrix-based wave-front shaping,” Phys. Rev. A **94**(4), 043830 (2016). [CrossRef]

**4. **A. Hanninen, M. W. Shu, and E. O. Potma, “Hyperspectral imaging with laser-scanning sum-frequency generation microscopy,” Biomed. Opt. Express **8**(9), 4230–4242 (2017). [CrossRef]

**5. **J.-X. Cheng and X. S. Xie, “Vibrational spectroscopic imaging of living systems: An emerging platform for biology and medicine,” Science **350**(6264), aaa8870 (2015). [CrossRef]

**6. **J. W. Goodman, * Speckle phenomena in optics: theory and applications* (Roberts and Company Publishers, 2007).

**7. **V. Ntziachristos, “Going deeper than microscopy: the optical imaging frontier in biology,” Nat. Methods **7**(8), 603–614 (2010). [CrossRef]

**8. **N. G. Horton, K. Wang, D. Kobat, C. G. Clark, F. W. Wise, C. B. Schaffer, and C. Xu, “In vivo three-photon microscopy of subcortical structures within an intact mouse brain,” Nat. Photonics **7**(3), 205–209 (2013). [CrossRef]

**9. **E. Beaurepaire, M. Oheim, and J. Mertz, “Ultra-deep two-photon fluorescence excitation in turbid media,” Opt. Commun. **188**(1-4), 25–29 (2001). [CrossRef]

**10. **M. Oheim, E. Beaurepaire, E. Chaigneau, J. Mertz, and S. Charpak, “Two-photon microscopy in brain tissue: parameters influencing the imaging depth,” J. Neurosci. Methods **111**(1), 29–37 (2001). [CrossRef]

**11. **N. Ji, D. E. Milkie, and E. Betzig, “Adaptive optics via pupil segmentation for high-resolution imaging in biological tissues,” Nat. Methods **7**(2), 141–147 (2010). [CrossRef]

**12. **O. Albert, L. Sherman, G. Mourou, T. Norris, and G. Vdovin, “Smart microscope: an adaptive optics learning system for aberration correction in multiphoton confocal microscopy,” Opt. Lett. **25**(1), 52–54 (2000). [CrossRef]

**13. **D. Débarre, E. J. Botcherby, M. J. Booth, and T. Wilson, “Adaptive optics for structured illumination microscopy,” Opt. Express **16**(13), 9290–9305 (2008). [CrossRef]

**14. **G. A. Tyler, “Reconstruction and assessment of the least-squares and slope discrepancy components of the phase,” J. Opt. Soc. Am. A **17**(10), 1828–1839 (2000). [CrossRef]

**15. **A. P. Mosk, A. Lagendijk, G. Lerosey, and M. Fink, “Controlling waves in space and time for imaging and focusing in complex media,” Nat. Photonics **6**(5), 283–292 (2012). [CrossRef]

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

**17. **Z. Yaqoob, D. Psaltis, M. S. Feld, and C. Yang, “Optical phase conjugation for turbidity suppression in biological samples,” Nat. Photonics **2**(2), 110–115 (2008). [CrossRef]

**18. **S. Popoff, G. Lerosey, M. Fink, A. C. Boccara, and S. Gigan, “Image transmission through an opaque material,” Nat. Commun. **1**(1), 81 (2010). [CrossRef]

**19. **M. Mounaix, D. Andreoli, H. Defienne, G. Volpe, O. Katz, S. Grésillon, and S. Gigan, “Spatiotemporal coherent control of light through a multiple scattering medium with the multispectral transmission matrix,” Phys. Rev. Lett. **116**(25), 253901 (2016). [CrossRef]

**20. **J. Bertolotti, E. G. Van Putten, C. Blum, A. Lagendijk, W. L. Vos, and A. P. Mosk, “Non-invasive imaging through opaque scattering layers,” Nature **491**(7423), 232–234 (2012). [CrossRef]

**21. **O. Katz, P. Heidmann, M. Fink, and S. Gigan, “Non-invasive single-shot imaging through scattering layers and around corners via speckle correlations,” Nat. Photonics **8**(10), 784–790 (2014). [CrossRef]

**22. **X. Yang, Y. Pu, and D. Psaltis, “Imaging blood cells through scattering biological tissue using speckle scanning microscopy,” Opt. Express **22**(3), 3405–3413 (2014). [CrossRef]

**23. **G. Stern and O. Katz, “Noninvasive focusing through scattering layers using speckle correlations,” Opt. Lett. **44**(1), 143–146 (2019). [CrossRef]

**24. **N. Antipa, P. Oare, E. Bostan, R. Ng, and L. Waller, “Video from stills: Lensless imaging with rolling shutter,” in 2019 IEEE International Conference on Computational Photography (ICCP), (IEEE, 2019), pp. 1–8.

**25. **G. Weinberg and O. Katz, “100,000 frames-per-second compressive imaging with a conventional rolling-shutter camera by random point-spread-function engineering,” arXiv preprint arXiv:2004.09614 (2020).

**26. **I. Freund, M. Rosenbluh, and S. Feng, “Memory effects in propagation of optical waves through disordered media,” Phys. Rev. Lett. **61**(20), 2328–2331 (1988). [CrossRef]

**27. **S. Feng, C. Kane, P. A. Lee, and A. D. Stone, “Correlations and fluctuations of coherent wave transmission through disordered media,” Phys. Rev. Lett. **61**(7), 834–837 (1988). [CrossRef]

**28. **G. Osnabrugge, R. Horstmeyer, I. N. Papadopoulos, B. Judkewitz, and I. M. Vellekoop, “Generalized optical memory effect,” Optica **4**(8), 886–892 (2017). [CrossRef]

**29. **B. Judkewitz, R. Horstmeyer, I. M. Vellekoop, I. N. Papadopoulos, and C. Yang, “Translation correlations in anisotropically scattering media,” Nat. Phys. **11**(8), 684–689 (2015). [CrossRef]

**30. **X. Xu, X. Xie, A. Thendiyammal, H. Zhuang, J. Xie, Y. Liu, J. Zhou, and A. P. Mosk, “Imaging of objects through a thin scattering layer using a spectrally and spatially separated reference,” Opt. Express **26**(12), 15073–15083 (2018). [CrossRef]

**31. **A. G. Vesga, M. Hofer, N. K. Balla, H. B. De Aguiar, M. Guillon, and S. Brasselet, “Focusing large spectral bandwidths through scattering media,” Opt. Express **27**(20), 28384–28394 (2019). [CrossRef]

**32. **L. Zhu, J. B. de Monvel, P. Berto, S. Brasselet, S. Gigan, and M. Guillon, “Chromato-axial memory effect through a forward-scattering slab,” Optica **7**(4), 338–345 (2020). [CrossRef]

**33. **J. Gateau, H. Rigneault, and M. Guillon, “Complementary speckle patterns: Deterministic interchange of intrinsic vortices and maxima through scattering media,” Phys. Rev. Lett. **118**(4), 043903 (2017). [CrossRef]

**34. **J. Gateau, F. Claude, G. Tessier, and M. Guillon, “Topological transformations of speckles,” Optica **6**(7), 914–920 (2019). [CrossRef]

**35. **I. Freund, “Looking through walls and around corners,” Phys. A **168**(1), 49–65 (1990). [CrossRef]

**36. **L. Zhu, J. B. de Monvel, P. Berto, S. Brasselet, S. Gigan, and M. Guillon, “Supplementary document for chromato-axial memory effect through a forward-scattering slab,” Optica **7**(4), 338–345 (2020). [CrossRef]

**37. **M. Kadobianskyi, I. N. Papadopoulos, T. Chaigne, R. Horstmeyer, and B. Judkewitz, “Scattering correlations of time-gated light,” Optica **5**(4), 389–394 (2018). [CrossRef]

**38. **J. Liang, J. Cai, J. Xie, X. Xie, J. Zhou, and X. Yu, “Depth-resolved and auto-focus imaging through scattering layer with wavelength compensation,” J. Opt. Soc. Am. A **36**(6), 944–949 (2019). [CrossRef]

**39. **S. Schott, J. Bertolotti, J.-F. Léger, L. Bourdieu, and S. Gigan, “Characterization of the angular memory effect of scattered light in biological tissues,” Opt. Express **23**(10), 13505–13516 (2015). [CrossRef]

**40. **X. Cheng, Y. Li, J. Mertz, S. Sakadžić, A. Devor, D. A. Boas, and L. Tian, “Development of a beam propagation method to simulate the point spread function degradation in scattering media,” Opt. Lett. **44**(20), 4989–4992 (2019). [CrossRef]

**41. **P. Arjmand and M. Guillon, *X Correlation in Forward Scattering Media* (2020). https://github.com/laboGigan/XCorrFSM.

**42. **M. Guillon, B. C. Forget, A. J. Foust, V. D. Sars, M. Ritsch-Marte, and V. Emiliani, “Vortex-free phase profiles for uniform patterning with computer-generated holography,” Opt. Express **25**(11), 12640–12652 (2017). [CrossRef]

**43. **R. Berkovits, M. Kaveh, and S. Feng, “Memory effect of waves in disordered systems: A real-space approach,” Phys. Rev. B **40**(1), 737–740 (1989). [CrossRef]

**44. **D. Oron, E. Tal, and Y. Silberberg, “Scanningless depth-resolved microscopy,” Opt. Express **13**(5), 1468–1476 (2005). [CrossRef]

**45. **D. Oron and Y. Silberberg, “Spatiotemporal coherent control using shaped, temporally focused pulses,” Opt. Express **13**(24), 9903–9908 (2005). [CrossRef]

**46. **E. Small, O. Katz, Y. Eshel, Y. Silberberg, and D. Oron, “Spatio-temporal x-wave,” Opt. Express **17**(21), 18659–18668 (2009). [CrossRef]

**47. **E. Small, O. Katz, and Y. Silberberg, “Spatiotemporal focusing through a thin scattering layer,” Opt. Express **20**(5), 5189–5195 (2012). [CrossRef]

**48. **A. Boniface, B. Blochet, J. Dong, and S. Gigan, “Noninvasive light focusing in scattering media using speckle variance optimization,” Optica **6**(11), 1381–1385 (2019). [CrossRef]

**49. **M. Hofer, S. Shivkumar, B. El Waly, and S. Brasselet, “Coherent anti-stokes raman scattering through thick biological tissues by single wavefront shaping,” Phys. Rev. Appl. **14**(2), 024019 (2020). [CrossRef]