## Abstract

The optical transmission matrix (TM) characterizes the transmission properties of a sample. We show a novel experimental procedure for measuring the TM of light waves in a slab geometry based on sampling the light field on a hexagonal lattice at the Rayleigh criterion. Our method enables the efficient measurement of a large fraction of the complete TM without oversampling while minimizing sampling crosstalk and the associated distortion of the statistics of the matrix elements. The procedure and analysis described here is demonstrated on a clear sample, which serves as an important reference for other systems and geometries, such as dense scattering media.

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

## 1. Introduction

In wave physics, the scattering operator of a sample links the incident waves to the ballistic and scattered outgoing waves [1,2]. By definition it considers only waves that propagate to, or from, the far field. If the sample has a slab geometry, the part of the scattering operator that relates the incident field to the transmitted one is called the transmission operator (TO). By choosing a basis of incident and outgoing modes one can construct a matrix representation of the transmission operator known as the transmission matrix (TM) [3–8]. A TM measurement that gives full information on the transmission operator requires a basis that spans all incident and outgoing modes of the wave field. The number of independent modes $N_s$ incident on an area $A$ on the surface of a slab is equal to the number of modes in a waveguide of the same area, $N_s = 2\pi A / \lambda ^2$ [5,9], where $\lambda$ is the wavelength. It is however not straightforward to find an orthogonal and complete basis on a finite sample in free space.

One of the reasons to measure the TM is to explore the distribution and properties of its singular values which, if the TM accurately represents the TO, correspond to the transmission eigenchannels of the sample. These channels are useful in imaging and communication [5,10] and are an important building block of quantum transport theory in scattering samples [11–15]. The predicted theoretical distribution of the singular values in a multiple scattering waveguide is a bimodal distribution known as the Dorokhov-Mello-Pereyra-Kumar curve [16,17] and the properties of transmission channels are a subject of recent intense research [7,10,18–32]. In an open optical system, the finite field of view and finite numerical aperture of the optics and the escape of energy parallel to the surface of the sample make it impossible to measure the complete matrix. The statistics of a partial transmission matrix has been shown to differ strongly from that of a full one [33,34]. This is less of a problem for guided waves such as microwaves [19,20], ultrasound [35,36] and elastic waves [21] where it is easier to access a large fraction of the scattered field. It is an ongoing challenge to accurately measure the largest possible fraction of the TM of a three-dimensional sample in an open optical system [8].

Popoff and coauthors [4] published the first measurements of the optical TM of a multiple scattering sample, using a spatial light modulator (SLM) to project light fields on a Hadamard basis and with an unmodulated part of the light as the phase reference. They demonstrated transmitting predetermined focused spots and images [37]. The low sampling density in their experiment precluded the observation of open channels or other mesoscopic correlations. Other methods employed a basis of spots in real or Fourier space, using a scan-mirror to enhance measurement speed as well as a separate reference arm to retrieve the transmitted field with a single shot method [27,34,38–40]. Interferometric methods are sensitive to phase drifts, and to this end it has been shown that the TM can also be measured without a reference beam [41]. However, this method requires a much larger number of measurements than the number of independent modes. Optical TM measurements close to the critical sampling density have been performed by scanning the beam on a square lattice [27,31,32,34,42]. When approaching the critical sampling density, where the number of sampled incident modes equals the total number of linearly independent modes $N_s$, bandwidth-limited input spots become correlated as the nearest neighbor spot wavefunctions (usually Airy disks) show increasing overlap, leading to distorted statistics of the TM. The amount of overlap depends on the lattice configuration.

In this paper, we demonstrate TM measurements using a hexagonal lattice of Airy disks spaced at the Rayleigh criterion, where the field overlap between neighboring spots is exactly zero and only small overlap exists between more distant spots. We note that measurements on an undersampled hexagonal grid far from the Rayleigh criterion have been demonstrated in Ref. [43]. We show that the hexagonal grid leads to a better representation of the statistics of the TM than square-grid sampling only close to the Rayleigh criterion. As in [27,34], we measure a polarization-complete TM. The spots are scanned in an outward hexagonal spiral which helps determine efficiently the point after which the beam escapes the microscope objective’s field of view or numerical aperture. We describe the experimental method and data analysis and perform a measurement of the TM of a zero-thickness medium, i.e., a clear medium where the focal planes of the input and output optics are the same. In this case the incident and transmitted fields are identical and the transmission operator is diagonal. Crucially, the orthogonality of the basis determines whether the transmission matrix is also diagonal, and we can check directly how the statistics of the TM are influenced by the sampling lattice.

## 2. Sampling the transmission matrix

Theoretically, the scattering matrix is expressed using a basis set of orthogonal flux-normalized far-field modes. We have natural basis sets in the case of a sample in a waveguide (flux-normalized waveguide modes) as well as in the case of a small scatterer in free space (flux-normalized partial waves) [10]. For a slab-type sample in free space the S-matrix is conventionally written as $S_{\sigma \boldsymbol {k} \sigma ' \boldsymbol {k'}}$, where $\boldsymbol {k}$ and $\boldsymbol {k'}$ are the wavevectors and $\sigma$ and $\sigma '$ the polarizations corresponding to the incoming and scattered waves respectively [1,44,45]. The scattering matrix element $S_{\sigma \boldsymbol {k} \sigma ' \boldsymbol {k'}}$ represents the amplitude scattered into the flux-normalized channel $\boldsymbol {k'}$ with polarization $\sigma '$ due to an incident flux with normalized amplitude in channel $\boldsymbol {k}$ and polarization $\sigma$. If $k$ and $k'$ are on opposite sides of a slab-type sample we speak of a transmission matrix element. On a finite sample, infinite plane waves are not a useful basis set and one has to resort to modes that are confined both in real space and Fourier space, such as diffraction limited spots.

Here we choose a real-space basis of diffraction limited spots that are generated by our microscope objective. We note that the TM can equivalently be measured in a $k$-space basis, where similar considerations apply. An objective fulfilling the Abbe sine condition [46] filled with a uniform collimated beam produces an Airy disk in the focal plane, which is a band-limited field defined by the circular shape of the pupil. The optimal sampling configuration for band-limited signals by a periodic lattice on a 2D plane is a hexagonal lattice (also known as a triangular lattice) [47]. If the density of sampling points is at least the critical or Nyquist density, corresponding to one sampling point per mode, the entire signal can be reconstructed from the sampled points.

On a square grid with a lattice constant $a$, as shown on the left image in Fig. 1, the diagonal neighbors have a separation of $a\sqrt {2}$. Consequently, the transmitted fields from two spots separated by $a$ or by $a\sqrt {2}$ yield different amounts of overlap. However, in our hexagonal sampling grid as illustrated on the right of Fig. 1, the six nearest neighbors are all at a distance $a$, while the six second neighbors are at a distance $a\sqrt {3}$. Critical sampling (within the band limit imposed by the NA) is achieved for the square lattice at

and for the hexagonal lattice atRemarkably, near critical sampling the nearest neighbor Airy spots on the hexagonal grid are almost exactly at the Rayleigh criterion

Next, we compare the effect of sampling density for both lattices, for simplicity assuming $\mbox {NA}=1$ and a single polarization component. A good measure to gauge these effects is the normalized rank of the basis of sampling modes, which is defined as the number of linearly independent modes divided by the total number of modes in the basis. We numerically compute the normalized rank of a basis of 933 spots within a circular area for the hexagonal grid and 931 spots within the same area for the square grid, as a function of the area density of sampling points. Computationally we find the normalized rank by constructing the matrix of inner products between the modes of the basis and calculating the fraction of its singular values $s$ that are larger than $c \cdot s_{\textrm {rms}}$, where $s_{\textrm {rms}}$ is the RMS (root mean square) singular value and $c=0.05$ is the numerical tolerance factor below which we consider a singular value to be effectively zero. Our results are quite insensitive to the choice of $c$.

Figure 2 shows the normalized rank as a function of the sampling density $N/A$, where $N$ is the number of sampling points, normalized to the critical density of $\pi / \lambda ^2$. The vertical line indicates the critical sampling density and the dashed line indicates the Rayleigh criterion for a hexagonal lattice. Ideally, the spots should be orthogonal and therefore lead to a normalized rank of 1. This is indeed the case for either lattice if we sample far below the critical sampling density. Close to the critical sampling density, the hexagonal lattice exhibits less spot overlap as indicated by a normalized rank closer to 1. At the critical sampling density, a hexagonal grid yields a normalized rank of 0.981 while a square one gives 0.935. In other words, close to critical sampling the square lattice induces three times more spurious near-zero singular values than the hexagonal one. Oversampling beyond the critical density can increase signal to noise but only with a square root law in measurement time. Without further data analysis and regularization procedures, it leads to additional spurious zero singular values in the TM as the neighboring lattice spots increase in overlap. Beyond a sampling density of 1.3, both lattices exhibit a similar amount of spurious zeros. This is not necessarily a problem since any critically or oversampled basis may be resampled into a strictly orthogonal set of functions that is complete within the FOV and NA, such as Bessel function modes. Such a resampling procedure requires accurate knowledge of the incident fields and introduces extra parameters and complexity in the data analysis which may be undesirable in real-time applications. The hexagonal grid of Airy spots spaced at the Rayleigh criterion is almost a complete basis, close to orthogonal and straightforward to implement in an experiment. To maintain symmetry in the system configuration, we choose to use the same basis in the detection plane, yielding a square TM. In case $\mbox {NA}\;<\;1$, the maximum achievable sampling density unavoidably scales with $\mbox {NA}^{2}$ irrespective of the sampling method or lattice, and the results in Fig. 2 are still valid with the $x$-axis now representing the scaled sampling density.

The polarization-complete TM can be composed from four measurements of sub-matrices as [48]

## 3. Measurement and data processing

#### 3.1 Apparatus

Our experimental setup to perform polarization-complete TM measurements is depicted in Fig. 3. Light from a Helium-Neon continuous wave (CW) laser (JDS Uniphase, 5 mW) with a wavelength of 633 nm first passes through a polarizer to ensure a linear and well defined polarization after which it is split into a signal (transmitted) and reference (reflected) beam. The signal beam is expanded, and directed with a two-axis scan-mirror (Newport FSM-300) towards the sample. A 4-$f$ relay telescope images the scan mirror to the pupil of the microscope objective O1 (Zeiss Achroplan 63x, 0.95-NA). A $\lambda /2$-plate on a motorized rotation stage controls the incident beam polarization. The transmitted light is collected with an oil immersion microscope objective O2 (Olympus PlanApo N 60x, 1.42-NA). The vertical and horizontal polarization components of the light are split by a polarizing beamsplitter (PBS) and imaged through a tube lens on two separate cameras C1 and C2 (Basler daA1280-54um) respectively, allowing simultaneous measurements in both polarization channels. The complex fields are retrieved using digital off-axis holography [49–51] by interfering the signal beam with a slightly tilted reference beam. The reference beam is delivered through a single mode polarization-maintaining fiber which ensures a near-Gaussian spatial profile. The reference beam is expanded to overfill the camera chip and provide a nearly flat intensity profile. The signal and reference path lengths are equal to within a few cm, which is crucial to avoid drifts due to changes of the laser frequency and coherence length. The $\lambda /2$-plate in the reference arm controls the ratio of the reference beam power directed to the cameras C1 and C2. Two light-emitting diodes (LED1 and LED2, 633 nm), which are switched off during the measurement, enable wide field illumination to locate desired regions of the sample and to focus the microscope objectives on the sample.

To measure the polarization-complete TM (Eq. (4)), we first set the incident polarization to H, after which the scan mirror steers the beam across the front surface of the sample. We start the spot scan at the center and then spiral outwards, which enables straightforward cropping of the resulting matrix to a smaller spatial area. For each incident spot, the cameras C1 and C2 record the V and H polarized transmitted light respectively. The process is repeated for V incident polarization. Dark and reference beam images, averaged over four images to even out the background noise, are taken before and after the measurement.

#### 3.2 Experimental sampling density

The effective area covered by the TM measurement is $A= {156}$ µm^{2}, which corresponds to the hexagonal envelope containing all measured points including half a unit cell outside the outermost spot centers so that the Airy spots belonging to this layer are contained entirely inside the measurement area. This area corresponds to $2\pi A / \lambda ^2 = 2443$ optical modes. However, our experimental TM has a dimension of $1838 \times 1838$, resulting in a sampling density of $1838/2443 \approx 0.75$.

The specified NA of our microscope objective O1 is $0.95$, but the effective NA retrieved by examining the first zero of the Airy spot is somewhat lower and estimated to be $0.9$ possibly due to apertures and slight aberrations in the optical system. The sampling density is therefore estimated to be at $0.75/0.9^2 = 0.93$ of the maximum sampling density allowed by the NA and close to the Rayleigh criterion at $0.988$. Hence, according to Fig. 2, we expect to be in a situation of slight undersampling but preserving the statistics of the TM very well with only about 1% of spurious zero singular values.

#### 3.3 Digital filtering of the transmitted fields

The transmitted fields are measured through the optical system with an NA of 1.4 (objective O2), which exceeds the NA of the incident fields. For the case of a clear medium, and in general whenever we want to obtain a square TM, we need the transmitted fields to be filtered to the same NA as the incident fields. We accomplish this by digitally Fourier filtering the measured fields. We filter with a supergaussian disc $SG$ with radius $R$,

where $r$ is the radial coordinate and the supergaussian exponent is chosen to be $n=50$. The supergaussian filter reduces the Gibbs phenomenon [52]. We set $R= {1.59}$ µm^{−1}, which corresponds to a filter that drops off sharply beyond $\mbox {NA}=0.95$, so that it just passes the spatial frequencies present in the incident beam.

#### 3.4 Sampling the transmitted field

We convert the measured and filtered fields to a vector by sampling at the positions that correspond to the spot centers of our incident field, which we obtain from the centers of mass (COM) of the transmitted images of the input Airy spots focused on a plain glass slide. The lattice constant of the hexagonal grid of the output basis is the mean separation of the COM of the input spots. We sample the fields at the output lattice points by linearly interpolating the field values from neighboring pixels. In this manner, the output sampling vector has the same size as the incident sampling vector which is equal to the number $N$ of input spots. The measured TM is therefore a square $N \times N$ matrix.

Due to a rotation of the CCDs and possible conformal distortion induced by optical elements, the input and output bases do not overlap exactly. This problem is usually ignored when measuring TMs of strongly scattering samples as there are no theoretical values of the matrix elements to compare to, but it leads to incorrect results for clear samples. To correct for this mismatch with a minimum number of parameters, we use an affine transformation to map the input grid to the measured grid of COM points.

A least squares algorithm is used to minimize the difference between the COM coordinates and the transformed grid coordinates. The result of this transformation is shown in Fig. 4. In (a), the generated (output) and measured (input) grids do not initially overlap. After performing an affine transformation on the generated lattice, Fig. 4(b) demonstrates that the input and output grids are superimposed to within 60 nm.

#### 3.5 Phase drift correction

During the course of the measurement slow phase drifts occur because of temperature, mechanical and laser wavelength drifts. These phase drifts affect the statistics and reproducibility of the measured TM. A typical measurement of a complete TM of dimension $1838 \times 1838$ elements takes approximately $11$ minutes, and phase drifts of order 1 rad occur in this time. Similar to the method used in [53], we account for these phase drifts by taking a reference image after every 20 scan positions of the input laser beam. The reference position is chosen to be the central spot of the hexagonal spiral. The phase of these reference fields is then compared to the phase of the starting spot by calculating the angle of their inner product. The phase drifts are interpolated with a one-dimensional smoothing spline fit. Finally, the computed corrections are applied to the retrieved fields to compensate for the phase drifts, reducing the inter-frame phase error to 0.07 rad.

## 4. TM of a zero-thickness reference

We study the accuracy and precision of our method by measuring the TM of a zero-thickness sample, i.e., objectives O1 and O2 are focused in the same plane. The focal plane we choose is the interface between a crown glass microscope cover slide and air. We note that this zero-thickness TM does not contain any near field terms according to the definition of the scattering matrix, and can consequently be measured with far-field optics. We scan 919 incident spots, corresponding to 18 complete hexagon layers across the interface, which leads to a polarization-complete TM of $1838 \times 1838$ elements.

A simulated TM is obtained by scanning a numerically generated Airy field on a hexagonal grid, spaced by the Rayleigh criterion just as in the experiment. The number of spots is 919, corresponding to a single polarization component as cross-polarization terms are absent in the numerical case.

The magnitude of one component (H incident, H transmitted) of the measured and simulated TM elements is shown in Fig. 5. The magnitude is normalized to the RMS singular value of $T_{\mathrm {HH}}$. In Fig. 5(a), the main diagonal of the TM is close to 1 and dominates the off-diagonal elements, with 84% of the total power in the diagonal elements. This indicates that the fields associated to different input spots are nearly orthogonal. Other lines, which are clearly visible due to the logarithmic color scale, emerge from the crosstalk between neighbor and next-neighbor Airy disks. This crosstalk is also visible, albeit at a slightly lower amplitude, in the simulated TM elements in Fig. 5(b). In this case 95% of the total power is in the diagonal elements. In the histogram in Fig. 5(c) we plot the distribution of the magnitude of the diagonal elements, which is peaked close to 1, as expected for a situation with only small crosstalk. The broadened width of the peak compared to the calculated one in Fig. 5(d) is attributed to experimental noise and imperfect imaging of the scan mirror on the pupil of the objective O1. In Fig. 5(e) we show the magnitude of the diagonal elements as a function of their input position. The nonuniform distribution favoring the upper positions of the grid indicates a slightly angle-dependent transmission of our optical system. In the numerical case displayed in Fig. 5(f), the diagonal elements exhibit very small fluctuations of around $1\%$ which arise due to the conversion from a rectangular pixel-based sampling to a hexagonal sampling. We conclude that our experimental data agrees very well with numerical simulations even though this direct comparison detects imperfections in the experiment very sensitively.

A zero-thickness glass sample is the ideal test case to verify whether the sampling procedure leads to distortions of the TM statistics. An important quantitative probe of these statistics is the distribution of the singular values [10,15]. The singular value histograms of the individual polarization sub-matrices and of the complete TM are plotted in Fig. 6. All histograms are normalized to the RMS singular value of the respective co-polarized sub-matrix. In Fig. 6(a,e), associated to polarization conserving measurements, the singular value histograms exhibit a pronounced peak at $s=1$ corresponding to completely transmitting channels. The histogram exhibits a low pedestal extending from $s=0$ to almost $s=2$. In Fig. 6(b,d), the singular value distributions for the cross-polarization matrices follow the Marchenko-Pastur law [54], which predicts the singular value probability function for a random uncorrelated matrix [55,56]. Since the glass sample does not cause appreciable cross-polarization scattering, we attribute the random values to a combination of camera noise and imperfections in the polarizing beamsplitter. In Fig. 6(c) we show the singular value histogram for the full experimental TM, and it resembles the $T_{\mathrm {HH}}$ and $T_{\mathrm {VV}}$ sub-matrices. The width is a little broader because of the noise from the cross-polarization channels. The singular value histogram of a simulated TM including 0.5% RMS Gaussian random noise on all matrix elements is displayed in Fig. 6(f). We observe that the width of the peak and the pedestal of the experimental histogram are well reproduced by this model. Without noise, as shown in Fig. 6(g), the simulated histogram is much narrower, but the pedestal remains. We conclude that the peak width is caused by experimental noise while the pedestal is a sampling feature. The low density of effectively zero singular values ($s\;<\;0.05$) indeed matches the results shown in Fig. 2. For comparison we plot in Fig. 6(h,i) the corresponding histograms for the case of square-lattice sampling with the same lattice constant without and with noise, respectively. The square-lattice sampling introduces much stronger spurious side lobes in the histogram, which obscure the true statistics of the sample TM. This proves clearly the advantage of sampling the TM on a hexagonal grid at the Rayleigh criterion.

## 5. Conclusion

In summary, we have measured the optical transmission matrix (TM) of a zero-thickness medium using a novel sampling method, namely a hexagonal sampling grid of Airy disks at the Rayleigh criterion. We show in numerical simulations that this sampling method induces far less distortions in the singular value statistics than the usual procedure of sampling on a square grid. We provide experimental verification of the statistics in a reference case of a zero-thickness sample. Experimentally verifying that sampling does not lead to distortion of the statistics is an important prerequisite for probing more complex samples such as scattering media. The advantage of sampling on a Rayleigh-criterion spaced hexagonal grid applies to any kind of waves including ultrasound and microwaves, and it can be used to measure the transmission matrix of fibers [42,57–60], waveguides [19] and for full scattering matrix measurements of samples in other geometries.

## Funding

Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) (68047618).

## Acknowledgments

The authors thank Sanli Faez and Dries van Oosten for fruitful discussions and Paul Jurrius, Dante Killian and Cees de Kok for technical support.

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **M. Nieto-Vesperinas and E. Wolf, “Generalized Stokes reciprocity relations for scattering from dielectric objects of arbitrary shape,” J. Opt. Soc. Am. A **3**(12), 2038–2046 (1986). [CrossRef]

**2. **A. Lagendijk and B. A. van Tiggelen, “Resonant multiple scattering of light,” Phys. Rep. **270**(3), 143–215 (1996). [CrossRef]

**3. **C. W. J. Beenakker, “Random-matrix theory of quantum transport,” Rev. Mod. Phys. **69**(3), 731–808 (1997). [CrossRef]

**4. **S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, “Measuring the transmission matrix in optics: An approach to the study and control of light propagation in disordered media,” Phys. Rev. Lett. **104**(10), 100601 (2010). [CrossRef]

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

**6. **M. C. W. van Rossum and T. M. Nieuwenhuizen, “Multiple scattering of classical waves: microscopy, mesoscopy, and diffusion,” Rev. Mod. Phys. **71**(1), 313–371 (1999). [CrossRef]

**7. **S. Rotter and S. Gigan, “Light fields in complex media: Mesoscopic scattering meets wave control,” Rev. Mod. Phys. **89**(1), 015005 (2017). [CrossRef]

**8. **E. van Putten and A. P. Mosk, “The information age in optics: Measuring the transmission matrix,” Physics **3**, 22 (2010). [CrossRef]

**9. **B. Saleh and M. Teich, * Fundamentals of Photonics*, Wiley Series in Pure and Applied Optics (Wiley, 2007).

**10. **D. A. B. Miller, “Waves, modes, communications, and optics: a tutorial,” Adv. Opt. Photonics **11**(3), 679–825 (2019). [CrossRef]

**11. **Y. V. Nazarov and Y. M. Blanter, * Quantum Transport: Introduction to Nanoscience* (Cambridge University Press, 2009).

**12. **E. N. Economou and C. M. Soukoulis, “Static conductance and scaling theory of localization in one dimension,” Phys. Rev. Lett. **46**(9), 618–621 (1981). [CrossRef]

**13. **H. Baranger, D. DiVincenzo, R. Jalabert, and A. Stone, “Classical and quantum ballistic-transport anomalies in microjunctions,” Phys. Rev. B **44**(19), 10637–10675 (1991). [CrossRef]

**14. **D. S. Fisher and P. A. Lee, “Relation between conductivity and transmission matrix,” Phys. Rev. B **23**(12), 6851–6854 (1981). [CrossRef]

**15. **J. Pendry, A. MacKinnon, and A. Pretre, “Maximal fluctuations–a new phenomenon in disordered systems,” Phys. A **168**(1), 400–407 (1990). [CrossRef]

**16. **P. Mello, P. Pereyra, and N. Kumar, “Macroscopic approach to multichannel disordered conductors,” Ann. Phys. **181**(2), 290–317 (1988). [CrossRef]

**17. **T. Martin and R. Landauer, “Wave-packet approach to noise in multichannel mesoscopic systems,” Phys. Rev. B **45**(4), 1742–1755 (1992). [CrossRef]

**18. **I. M. Vellekoop and A. P. Mosk, “Universal optimal transmission of light through disordered materials,” Phys. Rev. Lett. **101**(12), 120601 (2008). [CrossRef]

**19. **Z. Shi and A. Z. Genack, “Transmission eigenvalues and the bare conductance in the crossover to Anderson localization,” Phys. Rev. Lett. **108**(4), 043901 (2012). [CrossRef]

**20. **M. Davy, Z. Shi, J. Wang, and A. Z. Genack, “Transmission statistics and focusing in single disordered samples,” Opt. Express **21**(8), 10367–10375 (2013). [CrossRef]

**21. **B. Gérardin, J. Laurent, A. Derode, C. Prada, and A. Aubry, “Full transmission and reflection of waves propagating through a maze of disorder,” Phys. Rev. Lett. **113**(17), 173901 (2014). [CrossRef]

**22. **X. Hao, L. Martin-Rouault, and M. Cui, “A self-adaptive method for creating high efficiency communication channels through random scattering media,” Sci. Rep. **4**(1), 5874 (2015). [CrossRef]

**23. **S. Liew, S. Popoff, A. Mosk, W. Vos, and H. Cao, “Transmission channels for light in absorbing random media: from diffusive to ballistic-like transport,” Phys. Rev. B **89**(22), 224202 (2014). [CrossRef]

**24. **M. Kim, W. Choi, C. Yoon, G. H. Kim, S. hyun Kim, G.-R. Yi, Q.-H. Park, and W. Choi, “Exploring anti-reflection modes in disordered media,” Opt. Express **23**(10), 12740 (2015). [CrossRef]

**25. **M. Davy, Z. Shi, J. Park, C. Tian, and A. Z. Genack, “Universal structure of transmission eigenchannels inside opaque media,” Nat. Commun. **6**(1), 6893 (2015). [CrossRef]

**26. **A. Daniel, L. Liberman, and Y. Silberberg, “Wavefront shaping for glare reduction,” Optica **3**(10), 1104–1106 (2016). [CrossRef]

**27. **D. Akbulut, T. Strudley, J. Bertolotti, E. P. A. M. Bakkers, A. Lagendijk, O. L. Muskens, W. L. Vos, and A. P. Mosk, “Optical transmission matrix as a probe of the photonic strength,” Phys. Rev. A **94**(4), 043817 (2016). [CrossRef]

**28. **C. W. Hsu, S. F. Liew, A. Goetschy, H. Cao, and A. D. Stone, “Correlation-enhanced control of wave focusing in disordered media,” Nat. Phys. **13**(5), 497–502 (2017). [CrossRef]

**29. **P. Hong, O. S. Ojambati, A. Lagendijk, A. P. Mosk, and W. L. Vos, “Three-dimensional spatially resolved optical energy density enhanced by wavefront shaping,” Optica **5**(7), 844–849 (2018). [CrossRef]

**30. **P. Fang, C. Tian, L. Zhao, Y. P. Bliokh, V. Freilikher, and F. Nori, “Universality of eigenchannel structures in dimensional crossover,” Phys. Rev. B **99**(9), 094202 (2019). [CrossRef]

**31. **H. Yilmaz, C. W. Hsu, A. Goetschy, S. Bittner, S. Rotter, A. Yamilov, and H. Cao, “Angular memory effect of transmission eigenchannels,” Phys. Rev. Lett. **123**(20), 203901 (2019). [CrossRef]

**32. **H. Yilmaz, C. W. Hsu, A. Yamilov, and H. Cao, “Transverse localization of transmission eigenchannels,” Nat. Photonics **13**(5), 352–358 (2019). [CrossRef]

**33. **A. Goetschy and A. D. Stone, “Filtering random matrices: The effect of incomplete channel control in multiple scattering,” Phys. Rev. Lett. **111**(6), 063901 (2013). [CrossRef]

**34. **H. Yu, T. R. Hillman, W. Choi, J. O. Lee, M. S. Feld, R. R. Dasari, and Y. Park, “Measuring large optical transmission matrices of disordered media,” Phys. Rev. Lett. **111**(15), 153902 (2013). [CrossRef]

**35. **R. Sprik, A. Tourin, J. de Rosny, and M. Fink, “Eigenvalue distributions of correlated multichannel transfer matrices in strongly scattering systems,” Phys. Rev. B **78**(1), 012202 (2008). [CrossRef]

**36. **A. Aubry and A. Derode, “Random matrix theory applied to acoustic backscattering and imaging in complex media,” Phys. Rev. Lett. **102**(8), 084301 (2009). [CrossRef]

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

**38. **Y. Choi, T. D. Yang, C. Fang-Yen, P. Kang, K. J. Lee, R. R. Dasari, M. S. Feld, and W. Choi, “Overcoming the diffraction limit using multiple light scattering in a highly disordered medium,” Phys. Rev. Lett. **107**(2), 023902 (2011). [CrossRef]

**39. **M. Kim, Y. Choi, C. Yoon, W. Choi, J. Kim, Q.-H. Park, and W. Choi, “Maximal energy transport through disordered media with the implementation of transmission eigenchannels,” Nat. Photonics **6**(9), 581–585 (2012). [CrossRef]

**40. **W. Choi, C. Fang-Yen, K. Badizadegan, S. Oh, N. Lue, R. R. Dasari, and M. S. Feld, “Tomographic phase microscopy,” Nat. Methods **4**(9), 717–719 (2007). [CrossRef]

**41. **A. Drémeau, A. Liutkus, D. Martina, O. Katz, C. Schülke, F. Krzakala, S. Gigan, and L. Daudet, “Reference-less measurement of the transmission matrix of a highly scattering material using a DMD and phase retrieval techniques,” Opt. Express **23**(9), 11898–11911 (2015). [CrossRef]

**42. **M. Plöschner, T. Tyc, and T. Čižmár, “Seeing through chaos in multimode fibres,” Nat. Photonics **9**(8), 529–535 (2015). [CrossRef]

**43. **A. Boniface, I. Gusachenko, K. Dholakia, and S. Gigan, “Rapid broadband characterization of scattering medium using hyperspectral imaging,” Optica **6**(3), 274–279 (2019). [CrossRef]

**44. **M. Nieto-Vesperinas, * Scattering and Diffraction in Physical Optics*, 2nd Edition (World Scientific Publishing Company, 2006).

**45. **E. Wolf, “A scalar representation of electromagnetic fields: II,” Proc. Phys. Soc., London **74**(3), 269–280 (1959). [CrossRef]

**46. **E. Abbe Hon., “VII.–On the estimation of aperture in the microscope,” J. R. Microsc. Soc. **1**(3), 388–423 (1881). [CrossRef]

**47. **D. P. Petersen and D. Middleton, “Reconstruction of multidimensional stochastic fields from discrete measurements of amplitude and gradient,” Inf. Control. **7**(4), 445–476 (1964). [CrossRef]

**48. **S. Tripathi, R. Paxman, T. Bifano, and K. C. Toussaint, “Vector transmission matrix for the polarization behavior of light propagation in highly scattering media,” Opt. Express **20**(14), 16067–16076 (2012). [CrossRef]

**49. **E. N. Leith and J. Upatnieks, “Reconstructed wavefronts and communication theory,” J. Opt. Soc. Am. **52**(10), 1123–1130 (1962). [CrossRef]

**50. **M. Takeda, H. Ina, and S. Kobayashi, “Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry,” J. Opt. Soc. Am. **72**(1), 156–160 (1982). [CrossRef]

**51. **E. Cuche, P. Marquet, and C. Depeursinge, “Spatial filtering for zero-order and twin-image elimination in digital off-axis holography,” Appl. Opt. **39**(23), 4070–4075 (2000). [CrossRef]

**52. **Y. C. Eldar, * Sampling Theory: Beyond Bandlimited Systems* (Cambridge University Press, 2015).

**53. **M. Plöschner, B. Straka, K. Dholakia, and T. Čižmár, “GPU accelerated toolbox for real-time beam-shaping in multimode fibres,” Opt. Express **22**(3), 2933–2947 (2014). [CrossRef]

**54. **V. A. Marčenko and L. A. Pastur, “Distribution of eigenvalues for some sets of random matrices,” Math. USSR Sb. **1**(4), 457–483 (1967). [CrossRef]

**55. **M. Mehta, * Random Matrices*, ISSN (Elsevier Science, 2004).

**56. **J. Shen, “On the singular values of Gaussian random matrices,” Linear Algebra Appl. **326**(1-3), 1–14 (2001). [CrossRef]

**57. **S. Rosen, D. Gilboa, O. Katz, and Y. R. Silberberg, “Focusing and scanning through flexible multimode fibers without access to the distal end,” (2015).

**58. **A. Descloux, L. V. Amitonova, and P. W. H. Pinkse, “Aberrations of the point spread function of a multimode fiber due to partial mode excitation,” Opt. Express **24**(16), 18501–18512 (2016). [CrossRef]

**59. **R. N. Mahalati, R. Y. Gu, and J. M. Kahn, “Resolution limits for imaging through multi-mode fiber,” Opt. Express **21**(2), 1656–1668 (2013). [CrossRef]

**60. **N. Borhani, E. Kakkava, C. Moser, and D. Psaltis, “Learning to see through multimode fibers,” Optica **5**(8), 960–966 (2018). [CrossRef]