Transmission-matrix-based point-spread-function engineering through a complex medium

We report a method to design at will the spatial profile of transmitted coherent light after propagation through a scattering sample. We compute an operator based on the experimentally measured transmission matrix, obtained by numerically adding an arbitrary mask in the Fourier domain prior to focusing. We demonstrate the strength of the technique through several examples: propagating Bessel beams, thus generating foci smaller than the diffraction limited speckle grain, donut beams, and helical beams. We characterize the 3D profile of the achieved foci and analyze the fundamental limitations of the technique. Our approach generalizes Fourier optics concepts for random media, and opens in particular interesting perspectives for super-resolution imaging through turbid media.


INTRODUCTION
Generating a specific optical point-spread-function (PSF) has been one of the cornerstones of modern microscopy. This is conventionally done by inserting a phase or amplitude mask in the Fourier plane of the imaging system. For instance, Durnin et al. generated Bessel beams using a spatial filter for beam shaping [1]. Nowadays, holographic methods using a spatial light modulator (SLM) are the most versatile [2,3]. These techniques allow flexible beam shaping in a wide range of applications such as super-resolution microscopy [4,5], 3D microscopy [6,7], optical tweezers [8][9][10] and particle trapping [11]. However, all these studies typically require high quality optics and demand little or no sample aberrations.
Light propagation in materials with optical index heterogeneities is affected by scattering. In scattering materials, such as white paint or biological tissue, multiple scattering is at the origin of an intricate interference light field at the output of the medium, also known as speckle pattern [13]. Although the size of a speckle grain is diffraction-limited, this complex interference figure is detrimental for all conventional imaging systems. Nonetheless, this scattering process is linear and deterministic and thus even strongly scattering materials can be described by a transmission matrix (TM) [12]. Wavefront shaping techniques have recently emerged as a powerful technique to control the output field, using spatial light modulators or phase-conjugate mirrors [14]. SLMs provide up to several million degrees of freedom to design the input field at will, in order to control the propagation of light after the medium. Over the last decade, pioneer works have proven the capability to drastically increase light intensity at one or several output positions of a disordered system such as white paint [12,15], multimode fibers [16] or biological samples [17]. However, the resulting focal spots have sizes comparable to a speckle grain, thus diffraction-limited [18]. Recently, Di Battista et al [19] overcame this size limit by mechanically inserting an annular mask just after the scattering medium, prior to iteratively optimizing the focus intensity at a distance via wavefront shaping. Due to the filtering of the low spatial frequencies, the resulting speckle and optimized spot was narrower than the initial speckle. After removing the filter, the narrow spot -effectively a Bessel-like beam -remained intense, over a background speckle pattern wider than the focus.
Herein we report the first formulation of a TM-based operator that enables deterministic focusing after propagation through the medium, with a controllable PSF. We build this new operator by numerically applying a well-chosen mask (that can be arbitrarily designed in amplitude and/or phase) in a virtual Fourier plane of the output modes of the experimentally measured transmission matrix. We then demonstrate experimentally that a focus with the corresponding PSF can be obtained after the medium, by performing digital phase conjugation on this operator. To demonstrate the robustness and wide applicability of our technique, we generate and characterize a variety of useful PSF. Firstly, we generate a Bessel beam focus using an amplitude annulus mask, and show that its central FWHM is narrower than the size of a speckle grain, therefore demonstrating deterministic sub-speckle focusing without mechanical masking as in [19]. We also demonstrate donut-mode generation with various topological charges, and helical foci. Characterization of the of H enables focusing in any position of the imaging plane, with a size limited by diffraction [12]. 2' To go beyond the previous approach, a numerical step can be added to control the PSF. With a prior numerical computation of the optical field in a virtual pupil with a Fourier transform operation, the corresponding new operator H filt is obtained by numerically applying an arbitrary mask onto that pupil. 3' The DOPC of H filt allows focusing in the output plane after the medium, the shape of the focus being given by PSF (in this example a donut mode), defined by the Fourier transform of the arbitrary mask. As in [12], the complex focus stands over a speckle background.
axial properties of Bessel and helical PSFs shows he potential of the technique for 3D wave control.

PRINCIPLE OF THE EXPERIMENT
A particular class of wavefront shaping methods relies on the measurement of the optical transmission matrix (TM), denoted H in this paper, which contains the relationship between the input field and the output field [20]. Its complex coefficients h X'X connect the optical complex field at the output [X' = (x , y ) CCD pixel coordinates] to the input field [X = (x, y) SLM pixel coordinates] by: Experimentally the TM is measured by displaying a set of input fields on the SLM and recording the corresponding output fields on the CCD, and requires the medium to be stable during the whole measurement process (a few minutes in our case). The digital phase conjugation of the TM, H † where † stands for the conjugate transpose, enables focusing at any output position [12], or scanning the focus [21]. The resulting spot has a size given approximately by the spatial correlation of the output speckle, i.e. diffraction-limited [13,18]. It effectively sets to a common phase all contributions arriving at this position. However, it is possible to completely tune the phase and amplitude distribution of the k-vectors forming this focus, and therefore control at will the PSF. In Fig. 1 we detail how we build a new operator based on the experimentally measured TM to generate an arbitrary mask, and generate the corresponding PSF at the focus. We first numerically perform a two-dimensional spatial Fourier transform on the TM, noted F 2D , of every output field. We is the wave vector associated to X'. This numerical operation is equivalent to compute the TM in a Fourier plane of the output imaging plane. In order to generate a given PSF, we then multiply in the k-space the field in the pupil plane by mask M (amplitude and/or phase), corresponding to this PSF. We thus obtain a new numerically filtered coefficient in the Fourier domain: We then return to the spatial domain X'X by taking the inverse Fourier transform ofĥ filt K'X : The resulting operator H filt can now be used instead of H to perform focusing and scanning of a focus, but with the chosen  Fig. 2. Experimental setup. A cw laser (λ = 800nm) illuminates an SLM that modulates the wavefront of light before propagation through a multiple scattering medium (ZnO nanoparticles). Transmitted light is collected with another microscope objective, placed onto a motorized stage to scan the z-axis. The output speckle pattern is recombined with a reference plane wave from the same laser, with a beam splitter (BS). One polarization state of the output field is selected with a polarizer (P) and imaged on a CCD camera. (L): lens (focal distance 200mm) PSF. Fig. 2 sketches the experimental setup to measure the transmission matrix H, compute H filt , and generate and characterize the formed PSF in three dimensions. A cw laser (λ = 800 nm, MaiTai, Spectra Physics) is split between a reference path and a sample path. In the sample path, a phase-only SLM (LCOS-SLM, Hamamatsu X10468) subdivided in 64 × 64 macropixels is conjugated with the back focal plane of a microscope objective, which illuminates a scattering medium made of ZnO nanoparticles (thickness 100µm). another microscope objective is used to image the transmitted speckle . The objective is placed onto a motorized stage (Thorlabs, Z825B) to scan the imaging plane axially. For the TM measurement, the output beam is recombined with the reference on a beam splitter and the hologram is recorded on a charged couple device (CCD) camera (Allied Vision, Manta G-046).The reference beam is blocked during focusing and PSF characterization.

EXPERIMENTAL DEMONSTRATION
We first demonstrate the generation of a Bessel-like beam using an annular mask. The irradiance profile of the ideal beam is described by a zero-order Bessel function of the first kind, which propagates with an associated complex field ∝ J 0 (k r r) exp (−ik z z). The experimental result is represented in Fig. 3. The applied mask during the numerical filtering step is an annular amplitude mask, with an inner ring size 58% of the pupil size. This value has been found to be a good compromise between the loss in intensity at the output, the ability to detect the foci over the background speckle, and a significant narrowing of the FWHM of the central spot, as with a mechanical mask [19]. Using H filt † , we obtain in the imaging plane a Bessel-like focus, standing over a background speckle. We compare its FWHM to a standard focus obtain with the standard TM H † . In the same experimental conditions, the Bessel-like central lobe FWHM is 23% narrower than a standard speckle focus grain. Both intensity profiles are fitted to an ideal Bessel and Gaussian profile, respectively. The central lobe of a Bessel beam is surrounded by a decaying set of side-lobe rings. Each lobe carries approximately the same amount of energy as the central spot, whereas Gaussian beams contain 50% of their total energy within their FWHM [22]. Achieving a smaller FWHM than a diffraction-limited focus entails penalties such as loss of energy in the central peak (40% lower than a standard beam). Inherent to wavefront shaping in complex media techniques, the PSF is not perfect, but rather stands over a speckle pattern that remains in the background, with the peak to background ratio proportional to the number of input modes [15]. In the generation of engineered PSFs, an extra penalty is added due to the low transmission of the virtual annular mask, which steeply decreases the energy in the focus while the background remains the same. In order to detect the focus with sufficient SNR, a high number of degrees of freedom, i.e. SLM pixels, is required, here we used N = 4096 input pixels. We therefore demonstrate subspeckle focusing after propagation through a scattering medium, using a Bessel-like beam.
Since any arbitrary phase and/or amplitude mask can be computed and placed onto the virtual pupil field, we also demonstrate the generation of donut-shaped beams, which are closely related to single-ringed Laguerre-Gaussian beams (LG m 0 ). The corresponding mask to be applied in the numerical filtering step is a spiral phase plate distribution presenting a continuous and gradual phase change from 0 to 2mπ around the optical axis, where m is to the integer number of 2π cycles in the pupil plane [23] and corresponds to a topological charge. The ring pattern dimensions increase with m. Experimental results of focusing (LG m 0 ) beams, with m from 1 to 4, are presented in Fig. 4, using N = 1024 SLM pixels. The transmittance M of the numerical phase mask in the Fourier domain applied during the numerical filtering step, as defined in Eq. 2 and illustrated in Fig. 4, reads : where k 0 is the pupil size in k-space. Just as for with Bessellike beam of Fig. 3, the focus stands over a background speckle pattern. When m increases, the intensity is distributed over a wider area in the imaging plane. Therefore, for a fixed number of SLM pixels used, the total energy in the targeted area is about the Normalized intensity Fig. 4. Generation of donut modes of various diameters after propagation through the scattering medium. Diameter increases with topological charge m from 1 to 4 (related to the corresponding Laguerre-Gauss beam LG m 0 ). Applied masks in the pupil plane (left), intensity distribution retrieved with a CCD camera (center), and related intensity profile for each order (right).
same but distributed over a bigger area, while the background speckle remains at the same average intensity.
Several beam profiles, such as Bessel beams have also very interesting propagation properties along the Z-axis. Our setup enables the scan of the Z-axis thanks to a motorized stage under the collecting microscope objective. Experimentally, we can scan over a few Rayleigh ranges (z R = 10µm) on both sides of the focal plane. In Fig. 5a , yOz cross sections of Bessel-like beam obtained in the same conditions as in Fig. 3a are reported. As expected, we prove without ambiguity that the generated Bessel-like beam has a longer depth-of-focus during propagation along the z-axis relative to standard Gaussian beams, for which the depth of focus is the Rayleigh length. Here, we observe experimentally that our Bessel-like beam has a depth of focus 1.7 times longer than the Gaussian focus.
As a last example, we demonstrate the generation of doublehelix point spread functions (DH-PSF). This 3D-design has 2 dominant lobes in the imaging plane, whose angular orientation rotates with the axial (z) position [5,7]. This profile is realized by applying a particular phase mask during the numerical filtering step leading to H filt , illustrated in Fig. 5c, which reads [5]: is the number of vortices, (k j , θ j ) the position of the jth vortex and k 0 the pupil size. A stack of images taken on both sides of the focal plane illustrates the rotation of the DH-PSF in Fig. 5d. The pattern rotates approximately linearly with z, with a negative angle for negative z. The reference z = 0 corresponds to the focal plane position. The DH rotates from -30 • to +30 • between -2z R and +2z R . This type of profiled-PSF has been used for 3D localization with a single shot, where the rotation of the imaged spot is related to its depth [24].

DISCUSSION
We have demonstrated focusing with various PSFs, defined by arbitrary controlling a mask in phase and amplitude in the virtual Fourier domain. It is to be noted that, as in all methods relying on focusing through a complex medium, the focus stands over a background speckle resulting from the incomplete phase conjugation [15], the signal to background of the focus increasing linearly with the number of segments controlled on the SLM. However, our technique goes well-beyond spatial-shaping techniques relying on intensity optimization in the imaging plane; namely, our approach allows for fine control of the focus shape, in amplitude and phase, that cannot simply be achieved using optimization approaches on the intensity. In the latter, the maximal intensity scales inversely with the number of target points [15] or the size of the focusing area, as in acousto-optic and photoacoustic techniques [25,26]. Similarly, in our approach, the maximal intensity in the PSF decreases with its complexity, but the amount of energy in the PSF remain the same.
This PSF engineering method through complex medium is also not limited to amplitude and phase modulation in the Fourier domain, and could be for instance extended to arbitrary polarization mask if measuring a polarization-resolved transmission matrix [27], and even to more complex spectrally-dependent PSFs thanks to the spectral dependence of the speckle [28,29]. One can therefore generate PSFs that could not easily be realized with physical masks or SLMs. Moreover, the Fourier domain of the plane of interest may not be accessible in practice, as in [19] where the amplitude mask was placed after the medium and the speckle observed at a distance. In contrast, our technique would work in any plane, even at the output surface of the medium. Furthermore, a single transmission matrix measurement can be used to focus at different positions, and allows rapid switching between various PSFs, as the numerical filtering step and phase conjugation are realized a posteriori. Using a simple quadratic phase mask, the focusing plane can also be translated longitudinally at will.
On the downside, in order to perform an accurate spatial Fourier-transform, the reference beam during the transmission matrix measurement needs to be a well-defined plane wave, which requires a reference arm with interferometric stability during the measurement process. Additionally, the resolution of the generated mask is related to the sampling in the Fourier domain, and therefore depends mostly on the sampling in the imaging plane (on the CCD): generating an accurate mask requires sampling at least a few CCD pixels per speckle grain over an extended output spatial region of at least a few times the PSF size. Note that, however, the resolution of the mask is advantageously independent of the number and resolution of the SLM segments, that only affect the signal to background of the focus after phase-conjugation. It is also worth pointing out that the signal to background is crucially affected by the total transmission of the virtual mask (as pointed out by [19] for physical masks).

CONCLUSION
In conclusion, we have reported the first formulation of an operator, built upon the experimental transmission matrix, that enables deterministic focusing of any arbitrary PSF after propagation through a scattering sample. We have illustrated the strength of this technique by generating Bessel, "donut" and Double Helix beams through a scattering sample with a simple use of this new operator, and characterizing their transverse and longitudinal properties. The method can readily be extended to other complex media, from biological tissues to multimode fibers. The possibility of arbitrarily generate complex PSF through or in scattering media opens up new perspectives for many fields, in particular for 3D and super-resolution microscopy and for optical manipulation and trapping [30].

Funding Information
This research has been funded by the European Research Council ERC COMEDIA (278025). S.G. is a member of the Institut Universitaire de France. R.P. acknowledges NSF support under awards 1611513 and 1310487.