Controlling Light Through Optical Disordered Media : Transmission Matrix Approach

We experimentally measure the monochromatic transmission matrix (TM) of an optical multiple scattering medium using a spatial light modulator together with a phase-shifting interferometry measurement method. The TM contains all information needed to shape the scattered output field at will or to detect an image through the medium. We confront theory and experiment for these applications and we study the effect of noise on the reconstruction method. We also extracted from the TM informations about the statistical properties of the medium and the light transport whitin it. In particular, we are able to isolate the contributions of the Memory Effect (ME) and measure its attenuation length.

3 for optical systems, in particular for multiple scattering media, and detail the experimental acquisition method. We will show that the TM can be exploited for different applications. The most straightforward ones are focusing and image detection using well-established operators used for inverse problems in various domains (telecommunication [13], seismology [12], optical tomography [16,17], acoustics [18,19] and electromagnetics [20]). We will confront theory and experiment for applications and study the statistical properties of the TM. We will, in particular, study the effect of the ballistic contributions on the statistics of the TM and its consequences on image transmission. We will finally study how short-range correlation of the memory effect (ME) [21] affects the TM and how to quantify its correlation length with the knowledge of the TM.

Modeling transmission through a linear optical system
In free space, a plane wave propagates without being modified and its k-vector is conserved since plane waves are the eigenmodes of the free space propagation. In contrast, a plane wave illuminating a multiple scattering sample gives rise to a seemingly random output field. This output field is different for different input k-vectors. In order to use a scattering medium as an optical tool, one has to characterize the relation between input and output free modes. At a given wavelength, we will model such a system with a matrix linking the k-vector components of the output field to those of the input field. We will justify the validity and the interest of such a model.
For any linear propagation medium, the propagation of an optical wave is entirely characterized by its Green function [22]. If the sources are linearly polarized and the observation is made on the same polarization, one can use a scalar model. The scalar Green function G(r, r , t, t ) describes the influence of the optical field at position r at time t on the optical field at position r at time t. For a propagation in a medium stationary over time, the time dependence of the Green function is governed by t − t . For a surface S src containing all the sources, the optical field on r at t reads This expression can be written in the spectral domain: where E(r, ω) (resp. G(r, r , ω)) is the temporal Fourier transform of E(r , t) (resp. G(r, r , t)) at the pulsation ω. We want to characterize the relation between the optical field on a surface S src containing the sources and the optical field on the elements of an observation surface S obs . Experimentally, the sources and receptors have finite size. We denote by E out m = S out m E(r) dr the average optical field on the mth receptor of surface S out m ⊂ S obs . In a similar way, we denote by E in n = S in n E(r) dr the average optical field on the nth controlled area of surface S in n ⊂ S src . We denote by N and M, respectively, the number of sources and receptors. 4 Therefore, we can define the mesoscopic TM of an optical system for a given wavelength as the matrix K of the complex coefficients k mn connecting the optical field (in amplitude and phase) at the mth output element to the one at the nth input element element. Thus, we have In essence, the TM gives the relationship between input and output pixels at a given frequency, notwithstanding the complexity of the propagation for a stationary medium. In this paper, we will perform the singular value decomposition (SVD) of TMs to study their transmission properties. The SVD consists in writing where † denotes the transpose conjugate. The SVD decomposes the system into independent transmission channels characterized by their transmission coefficients and by their corresponding input and output modes. V is a unitary change of the basis matrix linking input freemodes with transmission channel input modes of the system. is a diagonal matrix containing real and positive values called singular values of K and denoted by λ m . These values are the square roots of the energy transmission values of the transmission channels. U is the unitary change of the basis matrix between transmission channels' output modes and output freemodes. For practical purposes, we will study the normalized singular values defined by The SVD is a powerful tool that gives access to the statistical distribution of the incident energy injected through the multiple transmission channels of the system. Within this matrix model, the number of channels of the observed system is min (N , M). Nevertheless, this number does not quantify the information that could simultaneously be conveyed in the medium.
Quantifying information is of high interest and this work has been done in acoustics [23,24] for time-reversal focusing through a multiple scattering medium. In those theories, the number of spatial and temporal degrees of freedom (N DOF ), also called the number of information grains, is introduced. These degrees of freedom represent as many independent information grains as can be conveyed through the system. For a monochromatic experiment there are only spatial degrees of freedom and N DOF is given by the number of uncorrelated elements used to recover information [2]. This number will be different depending on the type of experiment considered, focusing or image detection. A focusing experiment consists in finding the optimal input wavefront to send in order to maximize energy at a desired position (or a combination of positions) on the output of the sample. An image detection experiment consists in retrieving an image sent at the input of the system by measuring the output complex field. For focusing experiments N DOF is given by the number of independent input segments, which corresponds to the number of degrees of freedom the experimentalist can control. Conversely, for detection experiments, the experimentalist cannot access more than the number of independent output elements recorded. In the general case, we have N DOF N for focusing and N DOF M. We choose to fix the input (resp. output) pixel size at the input (resp. output) coherence length of the system. For this ideal sampling, N DOF is given by the input or output dimension of the TM according to the type of experiment: N DOF = N for focusing and N DOF = M for detection.

Experimental setup
To measure the TM of an optical system, one needs a way of controlling the incident light and recording the complex optical field. Control is now possible thanks to the emergence of spatial light modulators (SLMs) and the output complex field can be measured using an interferometric method detailed in section 2.3. We want to study an optical system governed by multiple scattering in which the propagation function cannot be a priori calculated. The material chosen is zinc oxide, a white pigment commonly used for industrial paints. The sample is an opaque 80 ± 25 µm thick deposit of a ZnO powder (Sigma-Aldrich 96479) with a measured mean free path of 6 ± 2 µm on a standard microscope glass slide. The thickness of the sample being one order of magnitude higher than the mean free path, light transport inside the sample is in the multiple scattering regime. To generate the incident wavefront, the beam from a diodepumped solid-state single longitudinal mode laser source at 532 nm (Laser Quantum Torus) is expanded and spatially modulated by a liquid crystal SLM (Holoeye LC-R 2500). Using a polarizer and an analyzer, we choose an appropriate combination of polarizations before and after the SLM to achieve up to a 2π phase modulation depth with less than 10% residual amplitude modulation. Only a part of the SLM is modulated and the rest of the surface is used to generate a static reference (details in section 2.3). The surface of the SLM is then imaged on the pupil of a 20× objective with a numerical aperture (NA) of 0.5. Thus a pixel of the SLM matches a wave vector at the entrance of the scattering medium. The beam is focused at one side of the sample. We image the output pattern 0.3 mm from the opposite surface of the sample through a 40× objective (NA = 0.85) onto a 10-bit charge-coupled device (CCD) camera (AVT Dolphin F-145B). We do not directly image the surface of the material to mitigate the effects of the possible presence of ballistic light. This point is detailed in section 4.4.
The duration of a typical measurement of the TM is a few minutes, much less than the typical decorrelation time of the medium (more than an hour). From now on, we will ignore any time dependence of the medium. It is an important prerequisite for the measurement and exploitation of the TM.

Acquisition method
Measuring a TM in optics raises several difficulties coming from the impossibility of having direct access to the phase of the optical field. Generally, one can have access to the complex optical field using interferences, for instance with a known wavefront and a full field phaseshifting interferometry method [25]. For any input vector, if the reference phase is shifted by a constant α, the intensity in the mth output mode is given by where s m is the complex amplitude of the optical field used as the reference in the mth output mode.
Thus, if we inject the nth input mode and measure I 0 m , I π/2 m , I π m and I 3π/2 m , respectively, the intensities in the mth outgoing mode for α = 0, π/2, π and 3π/2 give Commercial SLMs cannot control independently the phase and amplitude of an incident optical field. Our setup uses a phase-only modulator since the phase is the most important parameter to control in wavefront shaping [26]. A recent method [27] allows phase and amplitude modulation with similar commercial modulators; nevertheless, we use in this work a phase-only configuration to keep the setup as simple as possible.
Controlling only the phase forbids us to switch off or modulate in amplitude an incident mode. Thus we cannot directly acquire the TM in the canonical basis. We want to use a basis that utilizes every pixel with identical modulus. The Hadamard basis in which all elements are either +1 or −1 in amplitude perfectly satisfies our constraints. Another advantage is that it maximizes the intensity of the received wavefront and consequently decreases the experimental sensibility to noise [28]. For all Hadamard basis vectors, the intensity is measured on the basis of the pixels on the CCD camera. We obtain the TM in the input canonical basis by a unitary transformation. The observed TM K obs is directly related to the physical one K of formula (3) by where S ref is a diagonal matrix of elements, where s ref mm = s m represents the static reference wavefront in amplitude and phase.
Ideally, the reference wavefront should be a plane wave to directly have access to the K matrix. In this case, all s m are identical and K obs is directly proportional to K . However, this requires the addition of a reference arm to the setup, and requires an interferometric stability. To have the simplest experimental setup and a higher stability, we modulate only 65% of the wavefront going into the scattering sample (this corresponds to the square inside the pupil of the microscope objective as seen in figure 1), the speckle generated by the light coming from the 35% static part being our reference. S ref results from the transmission of this static part and is now unknown and no longer constant along its diagonal. Nevertheless, since S ref is stationary over time, we can measure the response of all input vectors on the mth output pixel as long as the reference speckle is bright enough. We will show in the next subsection that we can go back to the amplitude elements of S ref and that the only missing information is the phase of the reference. We will quantify the effect of the reference speckle and show that it neither impairs our ability to focus or image using the TM, nor does it affect our ability to retrieve the statistical properties of the TM.

Multiple scattering and random matrices
Under certain conditions, the rectangle N by M TM of a system dominated by multiple scattering [29] amounts to a random matrix of independent identically distributed entries of Gaussian statistics. We will study the properties of such matrices thanks to the random matrix theory (RMT) and investigate the hypothesis for which the TM of multiple-scattering samples are well described by this formalism. The modulated beam is focused on the multiple-scattering sample and the output intensity speckle pattern is imaged by a CCD camera. L, lens; P, polarizer; D, diaphragm.
We introduced previously that the SVD of a physical matrix is a powerful tool for studying the energy repartition in different channels. For matrices of independent elements of Gaussian distribution, RMT predicts that the statistical distribution ρ( λ) of the normalized singular values follows the Marcenko-Pastur law [30]. For γ = M/N , it reads with λ min = (1 − √ 1/γ ) being the smallest singular value and λ max = (1 + √ 1/γ ) the maximum singular value. In particular, for γ = 1, the Marcenko-Pastur law is known as the 'quarter circle law' [31] and reads Those predictions are valid for a TM of independent elements, i.e. without any correlation between its elements. This hypothesis can be broken by adding short-range correlation but also nontrivial correlations such as the ones appearing when introducing the flux conservation [32]. Since we only experimentally record and control a small part of the field on both sides of the sample, we will not be affected by this conservation condition.
Another straightforward parameter that could have an influence on the hypothesis of independent elements is the size of the segments used as sources (on the SLM) and receptors (on the CCD). If the segment size is less than the coherence length (i.e. the size of a speckle grain) two neighboring segments will be strongly correlated. The size of input (resp. output) segments 8 is set to fit the input (resp. output) coherence length. This way, the TM K should respect the hypothesis of independent elements and follow the Marcenko-Pastur law.
Nevertheless, the method of acquisition introduces correlation in the TM recorded. For M = N , K should respect the quarter-circle law but the methods detailed in section 2.3 give only access to the observed TM matrix K obs = S ref · K whose elements read To study the singular value distribution of our TM, we have to make sure that we can go back to the precise statistics of K using the observed matrix K obs . The effect of the reference S ref is static at the output. This means that its influence is the same for each input vector. Calculating the standard deviation of output elements over the input basis vectors, we have We made the hypothesis that each input k-vector gives an output speckle of the same mean intensity; it reads |k mn | 2 n = |k mn | 2 mn . This is valid in the diffusive regime as long as the illumination is homogeneous. We select for the illumination only the central part of an expanded Gaussian beam. We also neglected the correlation effects that are expected to be small regarding our system, such as the ME [21] (see section 5) and long-range correlations [33]. We can now define a normalized TM K fil whose elements k fil mn are normalized by |k obs mn | 2 m .
This matrix is not influenced by the amplitude of the reference and reads where S φ is a diagonal matrix of complex elements of norm 1 representing the phase contributions of the output reference speckle.
is a unitary matrix. Therefore, K fil and K have SVDs with the same singular values.
We show in figure 2 that the statistics of the observed TM K obs do not follow the quarter-circle (formula (10)) law, whereas the SVD of the filtered matrix K fil is close to it. The remaining correlation effects are due to inter-element correlations between nearby pixels. Taking a sub-matrix of the TM with only one element out of two, its statistics finally follows correctly the expected quarter-circle law. Those results are comparable to the ones obtained for transmission matrices in acoustics [29]. It is to be noted that we do not expect absorption to shows the matrix filtered to remove the reference amplitude contribution and • shows the matrix obtained by filtering and removing neighboring elements to eliminate interelement correlations.
affect the statistical properties of the matrix measured since absorption does not induce any new correlation and since the measure is not sensitive to the total energy conservation.
This clearly demonstrates that we can directly study the properties of K even when we have access only to the observed matrix K obs . We also want to point out that one has to be careful when interpreting deviation from a given prediction on SVD results. Once the matrix is measured and its singular values analyzed, one has to investigate the causes of the deviation from the Marcenko-Pastur law. Those causes could come from physical phenomena involved in light propagation in the medium that break the hypothesis of the absence of correlation in the TM (like the appearance of ballistic contributions; see section 4.4) but could also be inherent to the experimental method. However, K obs brings all information on singular values that describe the physics of the propagation in a given setup.
So far, we have compared the experimental SVDs for a symmetric case (M = N ). But what would be the effect of breaking the input/output asymmetric ratio M/N ? In the same setup, we experimentally acquired TMs with a fixed number of controlled pixels on the SLM (N = 1024) and changed the number of independent output pixels recorded by the CCD (M = γ N ) to achieve various ratio values (γ ∈ {1; 1.5; 2; 3; 4; 5}). Normalized singular value distributions of experimental TM filtered are shown in figure 3 and follow qualitatively the RMT prediction of equation (9).
We see that with increasing γ , the range of λ narrows. This effect has been measured in experiment in the acoustic field [34]. Physically this means that when γ increases, all recorded channels tend to converge toward the mean global transmission value λ = 1. When M > N , we record more independent information than the rank of the TM and averaging effects lead to a distribution that is peaked around a unique transmission value when M N . We notice that the experimental distributions deviate from the theoretical ones as γ increases. To have TMs with different values of γ , we record a large TM with γ = 11 and then create submatrices by taking lines selected randomly among the lines of the original TM. By increasing γ we increase the probability of having neighboring pixels in the TM and then are more sensitive to correlations between nearby pixels. This modifies the statistics and decreases the number of independent output segments. This explains the deviation from the theory. Other effects that could bring about correlations in certain experimental conditions are discussed is section 4.4.

Phase conjugation and degrees of freedom
It was demonstrated in [20,35] that one can take advantage of multiple scattering to focus on tight spots, thanks to the reversibility of the wave equation.
In such an experiment, the responses to short temporal signals emitted by an array of sources at a receiver at the output of a disordered medium are recorded. Those signals, linked to the Green's functions associated with the couples (source/receiver), are sent reversed in time. The waves generated converge naturally toward the targeted spot. The waves are focused in space and time. Phase conjugation is the monochromatic equivalent of the previous concept named time reversal. Those techniques are robust methods for achieving focusing and imaging. In this section, we will experimentally and theoretically study the efficiency of phase conjugation on focusing.
When using phase conjugation for focusing, it corresponds to shaping the input wavefront to put in phase the contributions of each input pixel at a desired output target. Denoting by E target out the output target vector, the input vector that estimates the desired output pattern E target out using phase conjugation is given by The effective output vector E out will be given by We will now demonstrate the dependence of the phase conjugation accuracy on N DOF for the simple case of a single focusing spot. We will compare the efficiency of perfect phase conjugation (with amplitude modulation) and of phase only phase conjugation. The intensity output field resulting from the phase conjugation focusing on the jth output pixel reads For m = j, the average intensity of the sum of the incoherent terms reads For m = j, the average intensity of the sum of the coherent terms reads We define the energy signal-to-noise ratio (SNR) for a single spot focusing as the intensity at the target point over the mean intensity elsewhere after phase conjugation. We directly have Similar results are obtained considering detection instead of focusing. If E out is the output pattern corresponding to an unknown input target mask E target in to be detected, similarly to equation (17), the reconstruction E image of the target is given by The reconstructed pattern for detecting a single target on the jth input pixel reads The SNR for the reconstruction of a single pixel reads

Experimental phase conjugation
To perform such an ideal phase conjugation focusing, one has to control the amplitude and phase of the incident light. Our experimental setup uses a phase-only SLM to modulate the light. Thus the input lth component of the input field to display in order to focus on the jth output pixel reads k * jl /|k * jl |. The resulting output intensity pattern for a phase-only experiment reads Calculating the intensity of incoherent and coherent sums in the same way as equations (19) and (20), we have We define SNR exp as the expected SNR that is given by Experimentally, for a setup and a random medium satisfying the hypothesis introduced in section 2.4 (independent elements of the TM with a Gaussian statistic), we have |k| 2 |k| 2 = π 4 ≈ 0.78 [36]. Finally, To achieve such an SNR exp using the experimental setup presented in section 2.2, one has to remove the reference pattern, after measuring the TM. If the reference part is still present, the intensity coming from the non-controlled segments on the SLM will increase the background level. This effect (detailed in the appendix) modifies SNR exp . In our experimental conditions, We will now use equations (16) and (22) to experimentally achieve focusing and detection. We tested phase conjugation for simple targets focusing and for point objects detection through the medium. The results are shown in figure 4.
To experimentally study the effect of N DOF on the quality of focusing, we measured the energy SNR for a single spot focusing. We present in figure 5 this experimental ratio as a function of the number of controlled segments N on the SLM and compare it to the theoretical prediction SNR ≈ 0.5N DOF . Those experimental results correctly fit the theoretical predictions.

Comparison with sequential optimization
The first breakthrough using phase-only SLM for focusing through a disordered medium was reported in [8]. The technique was a feedback optimization that leads to substantial enhancement of the output intensity at a desired target. We will present the principle of this approach and compare it to techniques requiring the measurement of the TM.
The experimental setup is roughly the same as that presented in figure 1. The wavefront phase is controlled onto the SLM on N groups of pixels, and the full output pattern is recorded by the CCD camera. A one-point target is chosen and the optical field E out m on this pixel is the linear combination of the effects of the N segments on the modulator: where k mn are the elements of the TM as defined in equation (3) and φ n is the phase imposed on the nth pixel of the SLM. The input amplitude term E 0 is constant for a phase-only SLM with homogeneous illumination. The idea is, for each pixel, to test a set of phase values and keep the one that gives a higher intensity on the target output point. The goal of this sequential algorithm is to converge toward a phase conjugation with an efficiency depending on the finite number of phases tested [8]. Different phase shifting algorithms were studied and their relative efficiencies in the presence of time decorrelation were analyzed in [37]; we will not discuss this problem in this paper.
In a single-spot focusing experiment for a sufficient number of phase shifts tested, the sequential algorithm and the phase conjugation computed with the knowledge of the TM are equivalent. Nevertheless, both techniques have different advantages and drawbacks. To be able to change the position or shape of the focal pattern, sequential optimization needs a new optimization for each new pattern. This forbids real-time applications. In contrast, once the TM is recorded, all information is known and this permits us to calculate the input wavefront for any output desired pattern, allowing, for example, to change the focus location at a video rate. Nevertheless, sequential optimization presents two strong advantages. With a setup comparable to the one shown in figure 1, sequential algorithm is more robust to calibration errors than the TM measurement since, even if the exact values of the phase shifts tested are unknown, the algorithm automatically selects the best phase mask. Another merit of this technique is that no interferometry measurement is needed. Since only the intensity is recorded, it allows experiments using incoherent phenomena such as focusing optical intensity on a fluorescence probe inside a disordered material [36].

Noise and reconstruction operator
Experimentally, one cannot carry out a measurement without noise. In our setup, the main noise sources are the laser fluctuation, the CCD readout noise and the residual amplitude modulation. Because of the error made in the TM measurement, one has to find, for focusing or image detection, reconstruction operators very resilient to experimental noise. These operators only give an estimation of the solution. From now on, we will make the distinction between the two sources of perturbation on the reconstructed signal: the experimental noise inherent to the measurement errors and the reconstruction noise inherent to the reconstruction operator. We will discuss the influence of those two types of noise on the performance of the technique used. Due to spatial reciprocity, focusing and image detection are equivalent. We will thus concentrate in this section on the case of image detection.
Inverse problems are widely common in various fields of physics. Optical image transmission through a random system is the exact analogue of multiple-input multiple-output (MIMO) information transmission through an unknown complex environment. This has been broadly studied in the last few decades since the emergence of wireless communications [13]. Similar problems have also been studied in optical tomography in incoherent [16] and coherent regimes [17]. In a more general way, inverse problems have been widely study theoretically for more than 50 years. Tikhonov [38] proposed a solution for an ill-posed system. The Tikhonov regularization method was then broadly used in various domains to approximate the solutions of inverse problems in the presence of noise.

Experimental noise
Prospecting an optimal operator for image reconstruction, the first straightforward option is the use of the inverse matrix K −1 obs for a symmetric problem (N = M) since K −1 obs K obs = I where I is the identity matrix. Inversion can be extended when N = M with the pseudoinverse matrix [K † obs · K obs ] −1 K † obs , denoted by K −1 obs . This operator provides perfect focusing but is unfortunately very unstable in the presence of noise. Indeed, the singular values of K −1 obs are the inverse ones of K obs . In other words, inverse operators tend to increase the energy sent through otherwise low-throughput channels. This way, the energy coming from all channels is normalized. Singular values of K obs below the noise level have the strongest contributions in K −1 obs but are aberrant. The reconstructed image is hence dominated by noise and bears no similarity to the input one. Another standard operator for image reconstruction or communication is time reversal. This technique has been studied in acoustics [18] and electromagnetic waves [39]. Being a matched filter [40], time reversal is known to be stable with regard to the noise level. Unlike inversion, it takes advantage of the strong singular values by forcing energy in the channels of maximum transmission. This way, singular values below the noise level carry little energy and do not degrade reconstruction. The monochromatic counterpart, phase conjugation, uses the operator K † obs (the operator used in section 3.2 for focusing and target detection). K † obs K obs has a strong diagonal but the rest of it is not null. We will see later that this implies that the fidelity of the reconstruction is not perfect and depends on the complexity of the image to be transmitted.
An intermediate approach is to use a mean-square optimized (MSO) operator. This operator, which we will call MSO and denote by W , is the Tikhonov regularization [38] for the linear perturbated systems mentioned earlier. It takes the experimental noise into account minimizing transmission errors, estimated by the expected value For an experimental noise of variance No σ on the output pixels, W reads This operator stabilizes inversion through the addition of a constraint depending on the noise level. This operator is intermediate between the two previous ones. If No σ = 0, W reduces to the inverse matrix K −1 obs , which is optimal in this configuration, while for a very high noise level it becomes proportional to the transpose conjugate matrix K † obs , the phase conjugation operator. For an intermediate noise level, the channels of transmission much greater than the noise level are used 'as in inversion' operator and channels much below the noise level are used 'as in phase conjugation'.

Reconstruction noise
The reconstruction noise is the perturbation brought by the reconstruction operator, even without measurement noise. This noise is intrinsic to the techniques used. The only operator to provide a perfect reconstruction, i.e. without any reconstruction noise, is pseudo-inversion (or inversion for M = N ). Nevertheless, this operator is not usable in a general case since it is very sensitive to experimental noise. Phase conjugation is stable in the presence of experimental noise but creates perturbations of the reconstructed signal. As seen in section 3.1, the SNR for a single focus is proportional to N DOF . The same goes for imaging a single spot; the target to detect generates a noise proportional to 1/N DOF in the rest of the image. For multiple targets, each target detection adds a perturbation for the reconstruction of the others; thus the intensity SNR is approximately given by N f /N DOF , where N f is the number of targets to detect. In other words, the fidelity of the reconstruction with phase conjugation is not perfect and depends on the complexity of the transmitted image. This limitation forbids the reconstruction of a complex image (with N f ≈ N ) by phase conjugation through a symmetric system (with N = M) since the reconstruction noise is of the same order of magnitude as the useful signal.
For any reconstruction operator Op, if K is known exactly (not affected by the experimental noise), the mean normalized reconstruction noise N rec for a single target detection can be estimated by N rec is defined here by the ratio between the mean intensity |Op × K | 2 mn m =n brought by the reconstruction of a target outside its position and the mean intensity |Op × K | 2 mn m=n at the target position. The average is over all input target positions. Graphically, N rec is the inverse of the SNR defined by the mean of the intensity of the diagonal of Op × K over the mean intensity elsewhere.
To study the reconstruction noise due to the MSO operator for different values of No σ , we show in figure 6 the simulation results for the reconstruction noise as a function of the No σ . It is important to understand that, for a given No σ , W is the ideal operator for an experimental energy noise level of No σ . In that study, the simulated TM K is noise-free. We see that for No σ = 0 (i.e. for the inversion operator), the reconstruction noise No rec is null and for No σ → ∞ (i.e. for the phase conjugation operator), No rec tends to N DOF (here N DOF = M = N ). This behavior is obvious when looking at the diagonal of |W × K | 2 compared to the other values.

Toward optimal reconstruction
In order to detect an input image with the best accuracy, one has to find a good balance between sensitivity to experimental noise and the influence of reconstruction noise on the recovered image. The mathematical optimum operator for those constraints is the MSO operator introduced in equation (32). To lower the reconstruction noise, one can only increase the parameter N DOF . We investigated two possibilities: averaging, and increasing the number of output modes M.
One possible way of virtually increasing N DOF is to perform the transmission of the same image through different realizations of disorder. With several realizations, the reconstruction noise, depending on the TM, will be different and averaging will lower perturbations on the final image. It is the monochromatic equivalent of using broadband signals, which takes advantage of temporal degree of freedom [24] where uncorrelated frequencies are used to increase N DOF . Iterating with different physical realizations of disorder is restrictive since it requires us to actively move or change the scattering sample. If we are only interested in detecting an amplitude object, another way to change the effective TM seen by the object is to illuminate it with different phase wavefronts. It is formally equivalent to transmitting the same image through different channels as if the image propagated through different realizations of disorder. For any amplitude object E obj with e obj m ∈ [0, 1], the effective incident field on the sample is S φ · E obj , where S φ is a diagonal matrix containing only phase terms. The output complex field pattern for a given TM K reads: with K = K · S φ . Physically, it is as if the same object E obj was transmitted through a different part of the medium. Another simple way to see it is that, by changing S φ (i.e. changing the 'virtual realization'), we change the projection of E obj on the input singular modes. Thus the same image propagates through different channels. For a given phase illumination, the amplitude object is estimated by E img : with W being the optimal MSO operator corresponding to the TM K for a given noise level. One now has to change the illumination wavefront S φ and average the different E img obtained to lower the reconstruction noise. To illustrate this effect we show in figure 7 simulation results for the reconstruction of a simple amplitude object using two different phase illuminations and averaging over ten different illuminations. The matrix is known exactly in the simulation, so all perturbations come from reconstruction noise. We see that the first two realizations reproduce the amplitude object with different noises with amplitude of the same order of magnitude. Averaging over multiple illuminations significantly reduces the noise background.  Since our SLM is a phase-only modulator, in order to simulate amplitude objects, we generate what we call 'virtual objects'. We use different combinations of random phase masks to generate the same virtual amplitude object E obj by subtracting two phase objects. From any phase mask E (1) phase one could generate a second mask E (2) phase where the phase of the mth pixel is shifted by s obj m · π. The optical complex amplitude of elements of the second mask are linked to those of the first mask by e (2) m = e (1) m · e is obj m π with e ( j) m being the mth element of E ( j) phase . We have the relation between the two masks and the amplitude of the 'virtual object': It is important to note that the norm of E (2) phase − E (1) phase is fixed and controlled whereas its phase is random, uncontrolled and differs from a 'virtual realization' to another. The term E amp = |E (2) phase − E (1) phase | can be estimated by where E (1) out (resp. E (2) out ) is the complex amplitude of the output speckle resulting from the display of the mask E (1) phase (resp. E (2) phase ). From a given set of random phase masks, we can define in the same way as in equation (34) a diagonal matrix S φ containing only phase terms. This term differs from one 'virtual realization' to another. This way, the reconstructed image reads This expression is formally equivalent to equation (35); thus generating 'virtual objects' using random phase masks is mathematically equivalent to illuminating a real object by a random phase-only wavefront.

Optimal reconstruction in the presence of ballistic light
We experimentally used this technique to virtually increase N DOF , and we average the results to lower the reconstruction noise. The complex amplitude image to detect is a grayscale 32 by 32 pixels flower image. The system is symmetric with N = M = 1024. The first step to reconstruct the input image with the best accuracy is to find the variance of the experimental noise No σ in order to obtain the optimal MSO operator. To that end, we numerically iterate the reconstruction operation for different values of the constraint in the MSO operator (equation (32)) and keep the one maximizing correlation between the reconstructed image and the image sent, hence obtaining an estimation of the experimental noise level. This step is made a priori knowing the image sent, but could be performed by any technique leading to an accurate estimation of the noise level. We have detailed in [15] the results in a purely diffuse system where ballistic contributions have no measurable effect. We study here the results in the presence of quantitative ballistic contributions.
In previous experiments, the focal planes of the two objectives are not the same. If the two focal planes coincide, if we remove the sample, an incident k-vector matches an output k-vector. Thus a pixel of the SLM corresponds to a pixel of the CCD. In the same configuration with a scattering medium, if the part of the energy ballistically transmitted is important enough, this ballistic part can be isolated in the TM.
To observe those contributions and study the accuracy of the reconstruction methods in that particular case, we experimentally use a thinner and less homogeneous area of the medium with a symmetric setup (M = N ) where input and output segments are of approximately the same size. We test the different reconstruction methods presented in section 4.3 with one realization and averaging over 40 'virtual realizations' with random phase masks. The results are shown in figure 8. As expected, we see that the inversion operation does not allow image reconstruction, even with averaging, whereas optimal MSO allows a 93% correlation for 40 averages (and a 60% correlation for one realization). Phase conjugation that is supposed to be stable in the presence of noise gives a weak correlation of 31% with averaging.
Those results are similar to the ones obtained when we measure the TM of a system governed by multiple scattering [15] for inversion (11% correlation with averaging) and for MSO operator (94% correlation with averaging) but are very different for phase conjugation (76% correlation with averaging). We see in figures 8(e) and (f) that most of the energy of the reconstructed field is localized on a very few points. In the presence of strong enough ballistic components, the direct paths associated give rise to channels of high transmission linking one input k-vector to the same output k-vector. Figure 9(a) shows normalized singular values for both experiments and for a simulated random matrix following the quarter-circle law. For the 'ballistic' experiment, two high singular values cannot be explained with a TM governed by multiple scattering only. We show for comparison in figure 9(b) the spatial aspect of the input singular vectors corresponding to the first and the 100th highest singular values, which should not be affected by balistic light. The 100th singular vector is spatially randomly distributed in energy (as expected for random multiple scattering), whereas the first singular vector mostly contains energy in one k-vector. This confirms that, in this particular case, the channel of maximum transmission is mostly dominated by the ballistic contribution.
Phase conjugation maximizes the energy transmission in channels of maximum transmission [40]. Therefore, ballistic high singular value contributions will be predominant in phase conjugation, regardless of the image E obj . Since the ballistic singular vectors are not 20 a.  homogeneously spatially distributed in energy, they will not efficiently contribute to an arbitrary image reconstruction. MSO is not affected by the contributions of ballistic waves, since it lowers the weight of the corresponding channels that do not efficiently contribute to the image reconstruction.

Pseudo-inversion
We showed in section 2.1 that when we set the output pixel size to fit the output coherence length, the number of degrees of freedom for detection is equal to the number M of segments  read on the CCD. This way, a second way to increase N DOF and thus the quality of the reconstructed image is to enlarge the output observation window. An important advantage is that the limiting time in our experiment is proportional to the number of input basis vectors N ; thus, in contrast with a focusing experiment where N DOF is equal [8] to N , increasing the degrees of freedom for image reconstruction does not increase the measurement time.
Another consequence of increasing the size of the image recorded, that is, of increasing the asymmetric ratio γ = M/N 1, is that the statistical properties of the matrix are deeply modified. We saw in section 2.4 that by increasing γ , the range of normalized singular values decreases [30,34]. A direct consequence is that the smallest nonzero singular value increases and reads λ 0 γ = (1 − √ 1/γ ). Controlling less information on the input than measured on the output results in an averaging effect. An interesting consequence is that for a reasonable experimental noise level, one can find a value of γ for which all singular values (and thus all channels' energy transmission) are above this noise level. In such a situation, the TM recorded a.
b. is barely sensitive to the experimental noise and since no singular values are drowned in the noise, the pseudo-inverse operator can be efficiently used.
To highlight this effect, we experimentally measure the TM of our sample for different values of γ 1. For each TM, we tested a single realization of the reconstruction of a complex image using the optimal MSO operator and the pseudo-inverse operator. The results are shown in figure 10(a). For γ = 1, we are in agreement with previous results since optimal MSO gives a partially accurate image when the results of inversion are totally drowned in noise. Increasing γ and thus N DOF strongly improves the quality of the reconstructed image for both techniques and reaches a >85% fidelity for the largest value of γ = 11, without any averaging. As expected, the curves of fidelity reconstruction as a function of γ for both techniques become identical when the value of λ 0 γ reaches noise level ( figure 10(b)). Above this value, pseudo-inversion and optimal MSO are totally equivalent. One notes that experimental λ 0 γ are always smaller than their theoretical predictions. We explain this deviation by the correlation induced in the matrix by the reference pattern |S ref | (see section 2.4). 23 We saw in this part that increasing the degrees of freedom and using an appropriate operator that takes noise into account permits us to drastically enhance the reconstruction accuracy of an image detection. Similar results can be obtained for focusing or for output wavefront shaping with the same techniques. Such experiments are more difficult to carry out since they require amplitude and phase modulation. Increasing γ = N /M and thus N in that case also presents the drawback of significantly increasing the measurement time.

Memory effect
To define the statistical properties of the TM of multiple scattering media in section 2.4, we did not take into account weak correlation effects that could break the hypothesis of independent elements of the TM. The strongest (first order) of those effects that can be recorded for a system dominated by multiple scattering (no ballistic contribution) in the diffusive regime is the socalled memory effect (ME). It brings correlation with characteristic length L me , which depends, in transmission, on the geometry of the observation system [21]. If the speckle grain size L speckle is of the same order of magnitude or smaller than the ME range L me , ME brings additional correlations in the TM.
We define L speckle as the full-width at half-maximum of the autocorrelation function of a speckle intensity pattern and L me as the attenuation length of the ME. Both lengths are defined for a given observation plane.
The speckle size L speckle depends on the distance between the back of the sample and the observation plane D, the width of the illumination area W and the wavenumber k. It reads: In contrast, L me depends on k, D and the thickness e of the sample and reads Therefore, for W e we have L me L speckle . This kind of geometry can be convenient since it allows us to move a focal spot along a distance greater that L speckle by just tilting the SLM [41].
We experimentally fixed D = 1.5 ± 0.2 mm, W = 0.5 ± 0.1 mm and the sample has a thickness e = 80 ± 25 µm. For those conditions, we expect L speckle ≈ L me ≈ 1.5 µm. The TM is measured. Each column of the TM is the complex amplitude response of a plane wave with a given incident angle. One can easily reconstruct the intensity patterns for different incident angles and calculate the correlation function between two patterns: where ⊗ m denotes the spatial convolution over the variable m and j is the variable of the correlation function corresponding to a displacement in pixel in the focal plane. We average results over all combinations of n, n with n − n = constant. The results are shown in figure 11. It is clear that the output speckle is coherently translated for distances greater than L speckle . L speckle is deduced from the mean autocorrelation function of the speckle pattern f (n−n )=0, j as a response to an incident wave plane. We found L speckle = 1.3 ± 0.48 µm. We define C(x) as the correlation coefficient (maximum of the correlation function) between two output speckles shifted by x in the observation plane. Exponential decay of the correlation C(x) ∝ e −x/L me is generally observed for an open geometry [21]. The measured C(x) (figures 11 and 12) fit an exponential decay with L me = 1.72 ± 0.13 µm.
We found the expected correlation lengths with a good order of magnitude for autocorrelation (speckle grain size) and ME and we were able to qualitatively separate both effects.

Conclusion
We reinterpreted in this paper previous works on optical phase conjugation in terms of a linear system defined by a monochromatic TM. Within this formalism, we show how to take advantage of multiple scattering to achieve image transmission and wavefront shaping. We also pointed out how to improve the quality of measurements in the presence of noise, taking advantage of RMT and information theory. We finally showed how to take advantage of the statistical properties of the TM to go back to physical phenomena.