Chromato-axial memory effect through a forward-scattering slab

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. Chromato-axial memory effect through a forward-scattering slab Longdi Zhu, Jacques Boutet de Monvel, Pascal Berto, Sophie Brasselet, Sylvain Gigan, Marc Guillon


INTRODUCTION
Imaging through or inside biological media is especially difficult due to the scattering nature of tissues wherein coherent light propagation results in complex speckle patterns. Imaging and remote focusing thus require beam engineering by typically controlling the phase of the beam in a Fourier plane using spatial light modulators. Optimization can be achieved in an iterative way [1] or by prior measurement of the transmission matrix [2]. Furthermore, nonlinear imaging by two-photon fluorescence excitation [3], second-harmonic generation [4], sum frequency generation [5], or coherent anti-Stokes Raman scattering [6], especially efficient for deep imaging, requires using broadband ultrashort pulses. In this case, not only the spatial profile of the beam must be controlled but also its spectral phase, in order to control the pulse shape. In these situations, full characterization of the hyperspectral transmission matrix allows perfect spatiotemporal wavefield control behind the scattering medium [7]. However, measuring the transmission matrix requires having access to the other side or the scattering medium, which should remain unchanged between the characterization and imaging.
Alternative approaches were then investigated to achieve such remote wavefield control by exploiting correlation properties of the transmission matrix. So-called "memory effects" allow speckle interferometry [8] and even coherent imaging [9] as well as fluorescent speckle imaging [10] through scattering samples. Until now, only spatial correlations of the transmission matrix were studied, especially in the case of thin diffusers (for which the transmission matrix is diagonal in the spatial domain [11,12]) and for forwardscattering media [13,14]. The spatial memory effect (ME) can be further increased by time gating the detection of an ultrashort pulse passing through a thick scattering medium [15]: early arriving photons are those experiencing fewer scattering events and propagating mostly in the forward direction. The mean arrival time of the pulse was clearly shown to increase with the output scattering angle [15] due to longer optical paths.
In the reciprocal spectral domain, a chromato-axial (χ -axial) ME was recently observed [16], yielding an axial displacement of a focused spot in response to a ∼100 nm chromatic shift through 1-mm-thick forward-scattering biological tissues. Through thin diffusers, such correlations could be shown to be useful for 3D imaging [17,18] but can be mostly accounted for by invariance properties of the Fresnel propagation equation [16,17,19]. Conversely, thick diffusers involve optical path lengths much longer than the wavelength, and narrow spectral correlations would be expected at the diffuser output. Therefore, the χ -axial ME observed experimentally in Ref. [16] cannot be solely attributed to the Fresnel diffraction equation but results from still unexplored intrinsic correlation properties of forward-scattering media. The χ-axial ME of thick forward-scattering media thus calls for both specific experimental characterization and theoretical modeling.
A scattering medium is typically characterized by its phase function p(θ ) describing the angular spread distribution at each scattering event and by its first momentum g = cos θ , its anisotropy factor. When g is close to unity, the impinging beam experiences an adiabatic broadening of its angular spread, and light dominantly propagates in the forward direction until the transport mean free path * = s /(1 − g ) (with s the scattering mean free path). Both time-gated detection of scattered light [15] and the χ-axial ME [16], which are similar phenomena related by a Fourier transform, are due to the correlations between the output scattering angle and the optical path length in the diffuser. Optical rays trajectory is usually identified to the Wigner transform of the field [20]. In the short wavelength approximation k s 1, the Wigner transform satisfies the radiative transfer equation (RTE) in scattering media. The average optical path in scattering media could be experimentally demonstrated to depend only on the geometry of the medium [21]. In the large medium limit (L * ), the diffusion theory approximation becomes valid [22], allowing the derivation of the distribution function of the optical paths [23,24]. Alternatively to the RTE, short-pulse propagation can be studied by considering the two-wavelength mutual coherence function (TWMCF) (also known as the time-dependent C (1) correlation function [25]) and was derived inside forward-scattering diffusers [26] yielding qualitative agreement with experiments [27] in the intermediate scattering regime L < s . Although both the RTE and the TWMCF models can be derived in the paraxial approximation for forward-scattering media [26,28], the optical path distribution has been little studied in the specific regime s L * of special interest for bio-imaging applications [29]. Here, the experimental evidence of a strong χ-axial ME is demonstrated over a 200 nm spectral bandwidth through a 250-µm-thick forward-scattering medium, under collimated illumination. The field propagation in the slab volume is modeled as a two-dimensional random walk in the transverse k-space and written as a path integral. The TWMCF is calculated in two different axial planes based on the joint probability density function of optical path lengths as a function of the output scattering angle. Analytical results predict a homothetic dilation of speckles with spectral detuning, in perfect agreement with experimental results. The center of the homothety is a virtual plane located at L/(3n) before the output surface, inside the scattering medium (where n is the mean refractive index contrast). Very large spectral correlation widths of the medium are also theoretically predicted, scaling as k = 3 √ 2/(L 2 0 ) (with 0 the average output scattering angle), in qualitative agreement with experiemental results. For comparison, the angular ME is also derived from our model and is shown to scale as = √ 3λ/(2L 0 ).

EXPERIMENTAL EVIDENCE
The χ-axial ME is experimentally evidenced in Fig. 1. In this experiment, a femtosecond titanium sapphire laser (Micra, Coherent) pumps a super-continuum generation photonic crystal fiber (femtowhite800, NKT photonics) at 800 nm. The supercontinuum beam is then collimated and sent to a homemade monochromator with a slit delivering spectral width of 10 nm for a spectrum ranging from 480 nm to 680 nm. The beam is finally spatially filtered, collimated, and sent onto a sample diffuser (paraffin layer or surface holographic). The speckle behind the diffuser is imaged onto a CCD camera using a Plan-Fluorite microscope objective (×4, 0.13 NA, Olympus) in combination with an = 740 nm, absorption α = 1 mm −1 ). However, in our collimated illumination conditions with beam diameters on the order of 1 mm, the contribution of surfaces to scattering leaves no observable ballistic light. This point is discussed further in Section 3.F. Scattering at surfaces thus makes the actual number of scattering events on the order of three for the L = 125 µm parafilm (one at each interface and one inside the bulk) and on the order of five or six when stacking two of them to get the 250-µm-thick sample. The scattering-spread along the propagation axis in the paraffin samples thus calls for an axially distributed scattering model.
In Fig. 1(a), images of speckle intensities I λ (r ⊥ , z) are shown to be similar for various wavelengths λ when appropriately tuning the distances z from the output surface of an L = 125 µm thick parafilm film (r ⊥ ∈ R 2 designates the twodimensional transverse vector). In Fig. 1 is plotted for the two considered paraffin-film thicknesses. C (λ, z) is defined here as the zero-mean correlation product. In Fig. 1(b), we arbitrarily considered reference speckles at λ ref = 571 ± 5 nm and z ref = 760 µm. The cross-correlation product of speckles was formerly noticed by I. Freund to allow interferometric measurements between two beams, altough only intensities are captured [31]. For zero-mean signals with Gaussian statistics [32], C can indeed be expressed as (1) C (λ, z) thus corresponds to a virtual interference involving the spatial phase only, which is all the more interesting when considering waves at different wavelengths: such interferences could not be achieved physically because of oscillation with time. In Fig. 1(b), a hyperbolic correlation trend is obtained between λ and the axial position z of correlation maximum. This observed behavior can be related to the Fresnel propagation expression [16,19]. In the frame of this approximation (for weak enough scattering angles), the intensity at coordinate z from the output surface of the diffuser can be expressed as where it appears that I k (r ⊥ , z) is unchanged when preserving the product λz constant. For a surface amplitude diffuser imprinting an achromatic output wavefield E k (r ⊥ , 0), a wavelength shift thus yields a perfectly homogeneous homothetic dilation along the propagation axis.
In Fig. 1(c), the spectral width of the correlation product along the hyperbolic line is observed to expand much beyond 200 nm for a paraffin film of thickness L = 125 µm, beyond the value obtained for a thin holographic diffuser and despite the larger thickness, the larger output scattering angle and the volumetric multiple scattering nature of paraffin.
Although the Fresnel propagation equation seems to account for the axial homothetic dilation behind the diffuser, the question of the center of the dilation is not clear for a volumetric diffuser: in Eq. (2), the output field E k (r ⊥ , 0) cannot be considered as achromatic. Furthermore, the large spectral correlation width measured experimentally is rather counterintuitive considering the volumetric multiple-scattering nature of the paraffin samples. The wavelength dependence of the field at the output of the diffuser is typically characterized by the spectral correlation width of the transmission matrix of the diffuser. Here, it appears that for such forward-scattering diffusers, the spectral correlation width cannot be estimated from a single 2D image, since the result would depend on the distance from the diffuser.
In the following, we propose a model to calculate the TWMCF after a thick forward-scattering diffuser, wherein the angular spread of the initial impinging collimated beam is adiabatically broadened by a volumetric multiple-scattering process.

MODELING OF A THICK FORWARD-SCATTERING DIFFUSER A. General Description of Our Modeling Approach
We describe the general approach for the model developed in the following sections. The model of the forward-scattering medium is first presented in Section 3.B, with corresponding hypotheses. The χ-axial ME and the angular ME are analytically derived in the frame of this model in Sections 3.E and 4, respectively.
First, the χ-axial ME is studied by expressing the TWMCF in the plane wave basis (Section 3.C). Relying on the ladder approximation [25], the TWMCF is studied by investigating correlations between the output angle of a plane-wave component of the output beam and the optical path length of its trajectory in the diffuser (i.e., in the sense of a "Feynman path"). The joint probability density function of and is then calculated (Section 3.D). The analytical expression of the zero-mean cross-correlation product defined in Eq. (1) is then given in Eq. (18), thanks to a Taylor expansion for small output scattering angles (Section 3.E.1). The χ -axial ME is then demonstrated, and the spectral correlation width is analytically derived in Section 3.E.2. Theoretical expectations are finally compared to experimental measurements in Section 3.F.
Second, in a similar manner, the angular ME is calculated in the frame of this model by using the same TWMCF in the specific case of two beams having identical wavelengths but different impinging angles inc (Section 4). Again, a joint probability density function is calculated (Section 5 in Supplement 1) in order to provide the average value and the variance of the phase delay φ( ) between the trajectories of these two impinging beams. The average value of φ( ) yields the angular ME and its variance the maximum tilt angle.

B. Model and Hypotheses
Here, the volumetric forward-scattering diffuser is modeled by a succession of weakly scattering events such as obtained by a stack of uncorrelated random phase plates made of a non-dispersive material [33]. Along the propagation direction, the correlation width of the refractive index map of the medium is then assumed to be much smaller than the slab thickness. Furthermore, backscattering is neglected under two assumptions: (i) the index contrasts are low enough, making backward scattering negligible at each scattering event, and (ii) the diffuser thickness is much smaller than the transport mean free path (L * ). We also assume that the output scattering angle is small enough to a allow second-order Taylor expansion of the phase. Together with assuming that the diffuser is made of an optically isotropic material, depolarization of the impinging beam at the diffuser output is neglected.
The scattering process is then discretized in multiple scattering events with step d [ Fig. 2(a)]. At each step, an impinging plane wave is partially scattered with a σ 2 0 -variance per transverse dimension [ Fig. 2(b)]. According to the central-limit theorem, for a large number of scattering events (i.e., much beyond the scattering mean free path (L s ), in the continuum approximation [22]), the scattering process can be treated as a Gaussian twodimensional random walk in the k-space [ Fig. 2(c)]. Noteworthy, no assumption is made about the angular distribution profile at each scattering plate, in particular about the fraction of ballistic light. At every step " p," the plane wave decomposition of the field between planes p and p + 1 may be written as where it is assumed that two-dimensional scattering angle vectors p ∈ R 2 all remain small enough to justify a second-order Taylor expansion of the phase ( p 2 stands here for the norm squared of p ). In this expression, z p is chosen to be zero at the  k (r ⊥ , z p = 0). At each step " p," the diffuser is modeled by a thin phase plate introducing an achromatic path delay δ p (r ⊥ ), so multiplying the impinging field by a transmission coefficient T In the Fourier domain, crossing a diffuser results in a convolution byT ( p) (k p ) = F(T p )(k p ). In addition, propagation between each scattering step results in a simple multiplication by the propagator function h k ( p ) = exp[ikd (1 − 2 p /2)]. The output field, after N scattering steps, can then be written as a multiple convolution product: where 0 = inc ∈ R 2 is the angle vector of the beam impinging onto the diffuser. In analogy with Brownian motion modeling, we will call a trajectory: t( ) = {( p ) p∈ [[1;N]] | 0 = 0, N = } a series of angle vectors p leaving the diffuser with the angle vector , after passing through the N scattering planes. For E (0) k ( 0 ) = δ( 0 ), Eq. (4) can then be rewritten as a path integral: where p = p − p−1 is the change in the propagation direction induced by the p th -diffuser according to the path t( ). For such a given trajectory t( ), the accumulated phase delay at the diffuser output is then and the field amplitude at the output of the diffuser is With these notations, at an axial distance z(= z N ) from the output surface of the diffuser, the fieldẼ k in Eq. (5) can be written as

C. Two-Wavelength Mutual Coherence Function
At the diffuser output, the correlation product defined in Eq. (1) can be calculated in the Fourier space, since Relying on the transverse-translation invariance of the system, we may make the following ergodic hypothesis that ensemble averaging is equivalent to space averaging over the r ⊥ -coordinates: In the limit N → ∞, the calculated average corresponds to the Wiener probability measure, over all possible Brownian paths. Making use of Eq. (8) yields The phase difference φ = φ k (t) − φ k (t ) is obtained from Eq. (6) and involves the difference between phases introduced by individual diffusers: Arg[T ( p) (k p )]. Since intermediate diffusers are random and independent, these phase terms are independent random variables uniformly distributed in ] − π, π ]. Therefore, the contribution of all pairs of trajectories involving different diffuser phase delays then vanishes when averaging over realizations. This assumption is known as the ladder approximation [25]. For two different wavenumbers k and k , only homothetic trajectories (t, t ) satisfying k p = k p for all scattering planes " p" will then contribute. Two such trajectories are shown in Fig. 2(a), differing only by a transverse scaling factor k /k. Since the input angle is 0, the former equality also implies k p = k p for all p, so yielding for φ: The field correlations at the volumetric diffuser output are expected to be chiefly driven by the phase term in the Fourier domain [34]. Assuming, in Eq. (10), that A k (t)A k (t ) does not depend on the considered path t but only on the output scattering angle , we now focus on the probability density function of φ.

D. Joint Distribution of Optical Path Length Delays
The phase difference φ in Eq. (11) involves the effective path The Nd term results just in a global constant phase that does not contribute to the correlation amplitude of fields. In order to reveal the correlations between the output scattering angle and the optical path length in the slab, we now calculate the joint probability density function of and : Here, it is useful to consider the characteristic function of ρ N according to : which can be interpreted as the average monochromatic output field at infinity. The random walk being considered as a Wiener process in the k-space, the function W N can be calculated relying on the equivalence principle between stochastic random variables and differential equations [35], as established by the Feynman-Kac formula [36]. In the continuous limit for large N values, the statistical average function W z can be shown to satisfy the following differential equation (see Section 2 in Supplement 1): where L = Nd is the diffuser thickness, and 2 0 = 2 x = 2 y = Nσ 2 0 is the angular spread of the beam at the diffuser output. It can be shown that in the frame of our model made of thin refractive diffusers (made of non-dispersive material), σ 0 is achromatic (see Section 1 in Supplement 1). This equation is a variant of the Schrödinger equation for the quantum harmonic oscillator whose solution is the (time-)propagator provided by the Mehler kernel [37]. In z = L, it yields (see Section 2 of Supplement 1) In the limit N → ∞, the probability density function ρ L ( , ) can be obtained from Eq. (15) by performing an inverse Fourier transform of W L relatively to k 0 [Eq. (13)]. Although ρ L ( , ) does not have an explicit expression, a color-coded illustration is shown in Fig. 3(a) from a numerical computation. The value of W L in the limit k 0 → 0 correponds to the integration of ρ L ( , ) along the coordinate, and yields the well-known Gaussian angular distribution at the diffuser output [ Fig. 3(b)]. E. χ-axial ME and Spectral Correlation Width

χ-axial ME
In the frame of the model's approximations, a closed form of Eq. (10) can be written as represents the spectral detuning. A Taylor expansion of W L in 0 then yields (see Section 3 in Supplement 1) We may further assume that for small enough spectral shifts, A k A k t equals just the average output Gaussian intensity profile I = ρ L ( ) = exp[− 2 /(2 2 0 )]/(2π 2 0 ). The TWMCF [Eq. (9)], obtained by integrating Eq. (16) over ( , ), can then be calculated (see Section 3 in Supplement 1). Its square modulus yields an analytical expression of the zero-mean cross-correlation product [Eq. (1)]: For a zero spectral shift (k 0 = 0), we retrieve the Lorentzian axial correlation profile. For k 0 = 0, C reaches its maximum for An axial homothetic dilation is then obtained, with the center at a virtual plane located at z = −L/3 (inside the diffuser), where the speckle is achromatic. Here, we considered a mean refractive index contrast equal to one between the diffuser and the external medium. For an average refractive index contrast n, the Snell's law must be applied to correct for the index mismatch at the output surface, and the achromatic correlation plane appears at −L/(3n).

Spectral Correlation Width of the χ -axial ME
We now consider the spectral correlation amplitude in the plane where speckle correlations are maximum. When Eq. (19) is satisfied, Eq. (18) yields At first order for small spectral detuning, we may approximate k 0 by k = k − k , so giving a spectral correlation width scaling as The spectral correlation width is limited by the statistical spread in that introduces a random phase. A typical random phase of amplitude φ rand = (k − k ) k k σ then degrades the correlation amplitude. From Eq. (17), the standard deviation σ of can be obtained as a function of (see Sections 3 and 4 in Supplement 1): resulting in the -dependent spectral correlation width. The maximum spectral correlation is thus obtained in the impinging beam direction and decreases as the observation angle increases.

F. Experimental Validation
To compare our theoretical model with experimental data, we plotted the zero-mean cross-correlation product as a function of z and κ = kz Fig. 4, we did so considering three reference speckles lying at plane distances z ref = 253 µm, 506 µm, and 760 µm, respectively (λ ref = 571 ± 5 nm). As a result, a spectral shift of the beam results in an axial shift of the pattern with a slope equal to 1 ± 0.03 between κ and the axial shift z. The coordinate z 0 at origin, corresponding to the achromatic plane location, was measured to be z 0 = 110 ± 10 µm. According to Eq. (19), the achromatic correlation plane should be at L/(3n) = 57 µm before the output surface of a L = 250 µm thick slab (based on the tabulated average refractive index of paraffin n = 1.45). A qualitative agreement is thus obtained with experimental observation. Regarding the spectral width, Eq. (21) predicts a spectral width λ λ 2 0 k/2π on the order of 400 nm for the . Again a qualitative agreement is obtained with the 200 nm width measured experimentally. Noteworthy, in our experimental case, the term k/k cannot be considered as close to unity, which may contribute to the discrepancy between the two values. Moreover, scattering at the diffuser surface, as well as the dispersive nature of parafilm over such a large bandwidth, is likely to decrease the spectral correlation width. Here, we attribute discrepancies between theoretical and experimental results mostly to the scattering events occurring at the slab surfaces. Within parafilm, scattering is thus not so well axially distributed as in our volumetric-scattering model. Nevertheless, although the slab thicknesses are on the same order of magnitude as the scattering mean free path [30] s = 170 µm, no remaining ballistic light is observed, even for the thinner 125-µm-thick sample. The reason for this seemingly paradoxical observation appears when re-melting the paraffin film between glass coverslips: when stuck to the coverslips, paraffin is much more transparent than when unstuck. In addition to volumetric scattering, two more or four more scattering events occur at the slab interfaces for the L = 125 µm and L = 250 µm samples, respectively.
In contrast with scattering media in the diffusive regime [38], no extrapolation length is to be considered, since output angles are assumed to remain small and much below the critical angle.

ANGULAR MEMORY EFFECT
In order to compare how the χ-axial ME and the angular ME scale with the physical parameters of the diffuser, we now derive the mutual coherence function for a given wavenumber k but for two waves, one of which impinges with an angle inc = 0 onto the diffuser. The mutual coherence function expressed in Eq. (10) (with k = k and z = z = 0) now involves the phase difference: where again, the ladder approximation is applied. The joint probability density function ρ( , X), with X = N p=1 p , can be easily calculated and yields a Gaussian-like distribution function (see Section 5 in Supplement 1). As a result, the conditional average of the phase difference in the direction can be shown to be The first term is just a constant phase, and the second one represents a tilt in the angular domain, so corresponding to a translation by an amount of L/2 at the exit face of the diffuser. An angular ME is then obtained with an actual rotation point located in the middle plane of the diffuser.
It is now of interest to calculate the amplitude of the angular ME. Similar to what has been done before, the correlation between the two beams is lost when the standard deviation of the phase difference φ is on the order of π . The variance of φ can be calculated (see Section 5 in Supplement 1) as The standard deviation of φ remains below π for inc ≤ √ 3 2 This result can be easily understood since L 0 is the size enlargement of a pencil beam impinging onto the diffuser. The maximum amount of path difference between two such pencil beams separated by L 0 is indeed λ.

DISCUSSION
Relying on the experimental observation that a χ -axial ME operates over a large spectral width in the case of forward scattering media, we derived an analytical model. The forward-scattered wavefield is written as a path integral since it follows a twodimensional random walk in the k-space. The Feynman-Kac formula, typically used to describe Brownian motion, then allowed us to derive the joint probability density function of optical paths. The TWMCF was then obtained for two wavelengths in two different transverse planes. Theoretical conclusions are the following: (i) the χ-axial ME is validated in the frame of the Fresnel diffraction approximation, even for thick forward-multiple-scattering media; (ii) the achromatic correlation plane lies at 1/(3n) of the sample thickness L, before the output surface of the scattering slab; and (iii) the average spectral correlation width decreases for larger scattering angles, and its average value is 3 √ 2/(L 2 0 ), where 0 is the output scattering angle. Noteworthy, although the scattering mean free path does not explicitely appear in the spectral width, it is actually included in the 0 value. Furthermore, if considering two diffusers yielding the same output scattering angle 0 , the χ-axial ME will perform over a larger spectral width for the thinner diffuser. This result does not contradict our experimental observation about the surface diffuser (exhibiting a smaller spectral range for a smaller scattering angle), since a surface diffuser does not satisfy the adiabatic angular broadening hypothesis. Finally, the angular ME is also derived in the frame of this same model, and the angular amplitude is shown to scale as √ 3λ/(2L 0 ). Here, we assumed that the diffusers were made of a nondispersive and non-absorbing material, which is definitely abusive over a spectrum ranging from 480 nm to 680 nm. Considering chromatic dispersion by the material could be introduced as an additional degree of refinement to quantitatively match experimental results, but both experimental and theoretical results yield very large spectral widths. Assuming a non-dispersive and nonabsorbing material also results in an achromatic output scattering angle and transport mean free path (see Section 1 in Supplement 1). This would not be true if considering amplitude masks as intermediate diffusers for instance. In this latter case, σ 0 would depend on k, and so would the output scattering angle 0 .
As a remark, the homothethy center located at 1/(3n) of the sample thickness for the χ-axial ME strongly recalls the result obtained in the case of the generalized geometric ME studied by Osnabrugge et al. [14], although the reason for any connection with these two results, if any, would definitely deserve specific attention.
Finally, the χ-axial ME is expected to be of high interest for imaging applications [16], since it even allows changing the focus of a beam inside the diffuser, the center of the homothety being at L/(3n) before the plane of interest.