Resampling the transmission matrix in an aberration-corrected Bessel mode basis

The study of the optical transmission matrix (TM) of a sample reveals important statistics of light transport through it. The accuracy of the statistics depends strongly on the orthogonality and completeness of the basis in which the TM is measured. While conventional experimental methods suffer from sampling effects and optical aberrations, we use a basis of Bessel modes of the first kind to faithfully recover the singular values, eigenvalues and eigenmodes of light propagation through a finite thickness of air. We further show that the latter can be mapped to the Fox-Li modes of a confocal resonator with finite-sized mirrors.


Introduction
The optical transmission response of a linear medium is described by the transmission operator (TO), which links all input modes to the transmitted output modes [1][2][3][4]. In a measurement, the TO is represented as a matrix in a specific basis of modes, be it on a square lattice as in the first reported measurement [5] or an optimized triangular lattice of free space modes [6], or a basis of waveguide modes [7,8]. In principle any system of free space modes [9] can be used to express the TM. In disordered systems, the TM is an important quantity because its statistics are related to the transmission eigenchannels of the sample, which are the building blocks of theories of transport and mesoscopic fluctuations [10][11][12][13][14]. This is still a hot topic of current research, with several recent experiments probing the transmission eigenchannels and their properties [15][16][17][18][19][20]. Other correlations within the TM, or between the transmission and reflection matrices [21], have been proposed for deep learning applications [22,23], parameter estimation [24,25] or used for imaging [26][27][28][29][30][31][32][33][34][35]. Moreover, a method to manipulate the TM of a fiber system directly was recently reported [36].
The statistics of the TM can be severely distorted if the sampling basis is non-orthogonal, resulting in the occurrence of spurious correlation that manifest as spurious high (as well as low) singular values. Since the basis sets that are accessible in optical experiments are typically non-orthogonal, resampling to an orthogonal basis set is necessary to obtain a singular value spectrum that is free of such artefacts. In an experiment that takes place in a non-lossy waveguide there is typically a complete orthogonal set of waveguide modes, which is the most logical choice for resampling [8]. On the other hand, in optical experiments taking place on slab-type samples, measurements are typically performed in a basis of spots in real space, approximate plane waves (spots in Fourier space), or random fields, neither of which are orthogonal.
Ideally, the TM should not only represent the singular value spectrum of the TO correctly, but also its eigenvectors and singular vectors. The correct reproduction of the eigenvalues can be taken as a second measure of an accurate retrieval of the TM, which concerns other relevant information than only the spectrum. For instance, while the singular value spectrum is quite insensitive against optical aberrations, the eigenvalues and eigenvectors are very sensitive since, as we show, they are closely related to the modes of an optical resonator that contains the medium. In this paper, we present a procedure that removes artefacts due to non-orthogonality and aberrations from a measured TM by resampling onto a set of orthogonal Bessel modes in an aberration-correcting way. We investigate the link between the eigenmodes of an optical cavity and the transmission eigenvectors of the TM, which result from a very similar eigenvalue equation, and hence, in a transparent system, have similar solutions and similar dependence on aberrations. We find that when the TM is measured straightforwardly in a basis of diffraction limited Airy spots [6], the propagation eigenvectors are strongly affected by system aberrations. When astigmatism is the leading-order aberration, the eigenvectors of the TM of a transparent system correspond to Ince-Gaussian modes (IGM), which are known as the third complete set of orthogonal solutions of the paraxial wave equation [9,[37][38][39][40]. They are also the resonant modes of a cavity with broken symmetry [9,37,38,41,42]. After the aberration-correcting resampling procedure the eigenvectors of the TM are no longer IGM but correspond to the Fox-Li modes of a perfect resonator with finite mirrors [43,44].

Resampling procedure
For any basis composed of a lattice of spots, the TM of a medium measured in such a basis will always be accompanied by some sampling crosstalk. For a hexagonal Airy spot basis sampled at the Rayleigh criterion, although nearest neighbors are chosen to have no overlap, the same cannot be said for second and further neighbors since the Airy spots have long tails. These nonzero overlap integrals between distant spots prevent the Airy basis from being completely orthogonal and therefore still mar the singular value statistics [6]. In this regard, an alternate choice for resampling the light fields is to use a basis of Bessel modes of the first kind, which is complete and orthogonal on an infinite 2D-plane and on a disk since the modes are solutions of the circular infinite square well [45] and of a cylindrical waveguide [46]. This is a problem normally studied in quantum mechanics [47] and in fiber optics, but we can use the same approach of constructing an orthogonal basis to express the TM of an open 3D system. The concept of such a resampling procedure is sketched in Fig. 1. Initially, the TM in an Airy disk basis is measured by scanning focused spots, in two orthogonal polarization states H and V, on a hexagonal lattice on the sample surface and recording simultaneously the transmitted fields in both polarizations with two cameras. The measured fields are sampled on the same hexagonal lattice to form the columns of the 2 × 2 TM in the Airy spot basis. The resulting matrix is subsequently resampled into a Bessel mode basis using the procedure detailed below.
In a basis of Airy spots, the spot-to-spot transmission matrix relates the input and output vectors ì and ì as ì = · ì .
To resample this matrix into a basis of Bessel modes of the first kind, we first have to define the circular region that contains the Bessel modes. We choose a circle that is inscribed inside the hexagonal scan area. In this way, we ensure that the entire circle is filled with spots that are scanned by the incident laser beam. In fact, we select a radius that is half a hexagon-layer-thickness larger than the outermost layer so that the Airy spots belonging to this layer are contained entirely inside the circle. We then set the appropriate boundary conditions that the Bessel beams fall to zero at the circle. We also truncate them to zero outside the circle since we do not have any control or information of this area where we do not scan input spots. This effectively determines the field of view of our system. Once the disk of interest is fixed, we find all the integer-order Bessel modes of the first kind that exist on the disk. A light field can be separated in two polarization components H and V , and each component can subsequently be decomposed in a basis of Bessel modes as where and are polar coordinates, is the Bessel mode of order , , is the th zero of the Bessel function, and , ∈ C are the contributing coefficients of each Bessel term. The expression above can equivalently be expressed using complex notation by combining the cosine and sine terms, but we choose to use a real basis because it is computationally more efficient to calculate and store the Bessel functions. We sort the Bessel modes in ascending order based on the eigenvalue of the Laplacian, which is equivalent to ordering them as a function of the radial zeros of the Bessel functions. A basis of dimension 2 is then formed by selecting the top Bessel modes for 2 orthogonal polarizations (defined in the pupil of the objective) H and V. The number is taken as the number of Airy spots contained entirely inside the circle. The low eigenvalue indices are more concentrated at the center of the circle, whereas the high ones live more on the edge. Also, as can be seen from indices (1,2) and (3,4), the Bessel modes of order > 0 come in pairs, where the spatial profile is the same but they are rotated with respect to one another. These pairs, degenerate in , are modes that possess the same radial number but opposite azimuthal numbers .
When the fields have been retrieved using digital holography [48][49][50], they can directly be decomposed into the Bessel modes. Since the Bessel modes are orthogonal, the decomposition is performed by taking inner products between the complex fields and the individual Bessel functions. This then yields a spot-to-Bessel transmission matrix . For each incident spot, the decomposition gives one column of . If Airy spots are scanned for each polarization and there are Bessel modes per polarization, the dimensions of are (2 × 2 ). This can mathematically be expressed as where ì is still an input field described in an Airy spot basis and ì is a transmitted field in a Bessel basis.
To obtain the TM in the Bessel-to-Bessel basis we resample the input Airy spot basis into Bessel modes. This is achieved using a resampling matrix such that ì = · ì [51].
with the TM in a Bessel-to-Bessel basis given by The difficulty arises in finding the resampling matrix . This matrix could be calculated as = ì | ì , but this would entail using theoretical Airy disks that do not contain the aberrations present in the optical setup. However aberrations typically exist in the experimental apparatus, especially on the input (wavefront generation) side. We perform resampling and correction of these aberrations on the input side in one step by studying the case of a TM where the incident and transmitted planes are the same, e.g. by focusing input spots on the surface of a microscope cover slide and imaging the same surface. In this situation, the incident field "is" the transmitted field, so ì = ì and consequently = 1 and = . The resampling matrix could contain low singular values which blow up when computing its inverse. To prevent this from occurring, we compute a pseudoinverse using a Tikhonov regularization [52], with a regularization parameter 2 = 0.1 as a reasonable threshold.
To recap, we find the TM in a Bessel-to-Bessel basis by first computing the spot-to-Bessel matrix and then multiplying it on the right by −1 , which is computed only once in the case of a zero-thickness TM where the incident and transmission planes coincide. The matrix is square with dimensions (2 × 2 ). By using the measured transmitted Airy spots in the resampling process, rather than the calculated spot functions, we not only resample the incident field but also correct for aberrations in the incident field optics.

Implementation for a zero-thickness reference
The spot-to-spot TM for a zero-thickness reference case is converted to a Bessel-to-Bessel TM by resampling the incident spots into a Bessel basis. When performing the analysis we select a circular region that is inscribed in the hexagon and the number of spots that lie completely in this circle is 778. Since we assume we are close to critical sampling, we choose 778 as the size of the Bessel basis per polarization component. Consequently the total number of Bessel modes is 2 × 778 = 1556. Theoretically, the number of modes of a slab-type sample should be approximately equal to 2 / 2 , where is the area of interest and is the wavelength of light [46,53]. In our optical setup where we use a 633 nm He-Ne laser, this translates to 2460 modes. However, we cannot measure the theoretically predicted full TM because of a finite field of view and a limited numerical aperture of the microscope objectives on either side of the sample [54,55]. The effective NA of our system is therefore √︁ 1556/2460 ≈ 0.80. The Tikhonov regularized pseudoinverse is computed on an independently measured zerothickness TM so that the noise in the singular vectors does not cancel out. The resampling yields a diagonal matrix containing some noise. The magnitudes of the elements of this Bessel-to-Bessel TM are depicted in Fig. 2(a). In Fig. 2(b), we show the magnitudes of the diagonal elements of only the HH sub-matrix because the behavior of the diagonal elements is the same for the HH and VV sub-matrices of a zero-thickness medium. It is clear that the Bessel modes with a low index are more prominent than the higher ones. The reason for this is that the high order Bessel modes live mostly on the edge of the measured area and thus contribute less. The diagonal elements of the simulated TM are also plotted, and we observe that as in the experimental case the magnitude of the diagonal elements drops off towards the end. However, the diagonal elements drop off less rapidly than in the experiment. In comparison, Fig. 2(c) and (d) depict the same results for the TM measured in an Airy spot basis at the Rayleigh criterion, discussed elaborately in Ref. [6]. Evidently, there is more crosstalk in the spot basis indicated by the presence of off-diagonal lines, and the experimental diagonal elements fluctuate due to an angle-dependent transmission in the optical setup [6]. The simulated result is free of these fluctuations since it is modeled aberration-free.
The histogram of the singular values is presented in Fig. 3(a). It resembles the one for the TM in the Airy spot basis, but here there is no pedestal on the right side of the peak. This demonstrates that using the Bessel mode basis eliminates the largest spurious singular values. The tail on the left side of the peak still persists because the low singular values result from the high order Bessel modes with a smaller contribution. Truncating the number of modes in the basis would cut this tail, but this would result in an even lower effective NA. Therefore, we do not reduce the number of modes since the high modes still have a non-zero contribution. If there were a singular value peak at zero, that would imply that the associated modes are redundant and we could remove them without lowering the effective NA. The other crucial statistics, viz. the eigenvalues of the TM, are displayed in Fig. 3(b) and (c). We observe that the width of the eigenvalue histogram is narrower than that of the singular values. With the Bessel basis, we ideally expect the singular values and eigenvalues to be identical. However, due to experimental noise is not exactly Hermitian and consequently the histograms in Fig. 3(a) and (b) are slightly different. Moreover, the major difference with the Airy basis is that the phases of the eigenvalues are not randomized anymore. The phase distribution shown in Fig. 3(c) is sharply peaked at 0, implying that the eigenvalues are almost completely real. Hence our resampling procedure, which uses the measured zero-thickness matrix, automatically corrects for the misalignment and aberrations of the optical setup. The corresponding statistics of the simulated TMs with 0.5% RMS random Gaussian noise on all elements are plotted in Fig. 3(d-f). The histograms confirm our experimental findings and reveal that the singular values and eigenvalue moduli are peaked at 1 and the eigenvalue phases at 0. Crucially, the difference between the two representation bases is clearly visible -the pedestal towards high singular values disappears after resampling. We checked that this pedestal is purely a sampling artefact and not an effect of noise. From the good agreement between the experimental and simulated results, we conclude that our resampling procedure is very effective at removing artefacts such as the spurious pedestal of high singular values.

TM eigenmodes and their parallel in the Fox-Li analysis of a confocal cavity
In cavity physics, a characteristic quantity of interest is the set of eigenmodes of the cavity. One way of finding the eigenmodes of an open cavity is with the computational method proposed by Fox and Li [43], where the field on the cavity mirrors is successively evaluated by solving diffraction integrals after every round trip. In fact, field propagation inside the cavity can be described by a Fox-Li operator whose eigenvectors are the cavity modes. In a confocal cavity with finite-sized mirrors, as sketched in Fig. 4, the field on the curved where 12 air and 21 air represents propagation through air from Mirror 1→2 and from Mirror 2→1 respectively.
To make a connection between a confocal cavity and the TM of a medium, it is reasonable to investigate the eigenvectors of the TM. The analogy with the TM setup is that there is a real space filter given by the finite field of view and a Fourier space filter given by the finite numerical aperture of the microscope objectives. For the TM of air, we can in fact find a direct relation between the transmission operator and the Fox-Li operator. The transmission operatorˆcan be expressed asˆ= 2 air 1 , where 1 and 2 are filtering operators by objectives 1 and 2 respectively. The eigenvectors and eigenvalues ofˆtherefore satisfy the equation If we introduce = 1 and multiply both sides of the equation by 1 from the left, we have Renaming 1 2 = and 2 = , we can simplify the physics of the transmission operator to ( air ) ( air ) = .
(transmission modes) For air we have 12 air = 21 air = air , and by comparing Eqs. (5) and (9), we notice that if 1 = 2 = the Fox-Li and transmission eigenmodes satisfy identical equations and consequently experience similar physics. In other words, predictions and experimental observations made from the TM of any medium that fulfils the reciprocity relation 12 = 21 can be mapped onto a confocal cavity with identical and possibly imperfect mirrors filled with the same medium; and vice versa. For example, given the principal modes of an optical fiber [57,58], it is possible to construct a confocal resonator that has the exact same mode structure.
The process of computing the eigenmodes between the mirrors in the resonator is similar to calculating the eigenstates of the TM. Additionally, the mode at the end of round trips is accompanied with an eigenvalue . Since we have 0 ≤ | | ≤ 1 without a gain medium in the cavity, the eigenvectors with smaller eigenvalues disappear more rapidly as the number of round trips increases. The loss encountered for every round trip in the cavity is attributed to higher order modes more than lower order ones. In this regard, the cavity is similar to the TM setup where high order modes fall outside the region of interest. With a greater air slab thickness, the high order modes rapidly obtain very small eigenvalues and only the low order modes survive. The choice of thickness of the air slab is thereby a suitable method to achieve mode filtering.

Propagation modes in air
To investigate the link between a TM and a free space cavity, we measure the TM of a finite thickness of air. We perform this in both an Airy spot basis and the aberration corrected Bessel mode basis and compare, in relation to the true Fox-Li modes of the cavity, the effect of the resampling and aberration correction procedure on the propagation modes.

Ince-Gaussian modes in a spot basis
We measure the complete 1838 × 1838 TM (2 incoming and 2 outgoing polarizations) of a 5 µm-thick air layer using the procedure described in Ref. [6]. Since propagation of light through air maintains polarization, it is sufficient to study just one sub-matrix of the complete experimental TM with the same incoming and outgoing polarization (as defined in the pupil plane of the objectives). The eigenvalue spectrum of the HH sub-matrix is plotted in Fig. 5(a). There are two outstanding features that can be observed in these spectra. The first distinct feature is that the eigenvalues spiral towards zero, and the second feature is that there exists a nonuniform separation between the eigenvalues along the spiral.
To understand these salient galaxy-like characteristics better, it is instructive to look at the corresponding eigenvectors. Their spatial profiles are obtained digitally by superposing the experimentally recorded transmitted field responses of the incident Airy spots weighted by the appropriate components of the column eigenvectors. The eigenvector corresponding to the highest eigenvalue is a single spot (Fig. 5(b)), which is the fundamental mode of the system. The next two largest eigenvalues are grouped together and the spatial profiles of their eigenvectors are two lobes with a central node (Fig. 5(c,d)). The only difference in their spatial profiles is that they have a nearly orthogonal orientation with respect to each other. We see the formation of a trend, with the next three eigenvalues also clubbed together. The eigenvectors are indeed the first few of a larger family of modes. Higher order modes, some of which are depicted in Fig. 5(e-g), have a larger spatial profile than the first few orders and therefore distortions in their symmetry are more evident. It is in fact observed that the higher the mode order, the more noisy the profile. For some of the highest orders (not plotted here), their spatial profiles are not distinguishable as they are drowned out by the experimental noise.
The numerically calculated eigenvectors of the same air slab are calculated and the ones corresponding to (b,c,e-g) are shown in Fig. 5(h-l). The fundamental mode (h) and first order mode (i) are similar in shape, but not in size, to the ones measured in our experiment. The higher order modes (j-l) also look similar to our measurements and exhibit elliptical symmetry.  However, they occupy the entire hexagon in the simulation while they are restricted to one side in the experiment. We established that this can be caused by an asymmetry in the optical system, such as a small intensity or phase gradient over the aperture. The family of modes found from the TM of an air slab with a thickness of 5 µm are neither Hermite-Gaussian nor Laguerre-Gaussian because all modes don't possess rectangular or circular symmetry respectively. The modes closely match a continuous transition between Laguerre-Gaussian and Hermite-Gaussian modes, known as Ince-Gaussian modes (IGM) [37,38]. These modes form a complete set of modes that fulfill the paraxial wave equation and exhibit inherent elliptic symmetry. Additionally, they are the eigenmodes of an astigmatic cavity [38].
The known properties of the Ince-Gaussian beams explain the spiraling eigenvalue plots. First, they possess a Gouy phase which is a function of the mode order. Ideally, modes of nonzero order come as degenerate pairs. However, although we do observe the bunching of the experimental eigenvalues in Fig. 5(a), they do not overlap in the complex plane. We conclude that aberrations or other imperfections lift the degeneracy of the experimental eigenvalues. The second relevant property of the Ince-Gaussian beams is that the spatial profile of the modes grows in size as the mode order increases. Consequently, for large modes, some part of their profile falls outside the region where the TM is measured. The extent that crosses the boundary is not taken into account when generating the TM, and its information is consequently lost. This loss in turn is translated as a reduced amplitude of the eigenvalue. Hence, the eigenvalue amplitude decreases as the mode order increases and thereby elucidates the inward spiraling of eigenvalues towards zero. Since a misaligned setup will endure more losses than an aligned one, the rate at which the eigenvalues spiral towards zero is an indicator of the alignment and aberrations of the optical setup.

Fox-Li modes in a Bessel basis
The illumination and imaging part of our experimental setup is ideally cylindrically symmetric because we have microscope objectives with circular apertures. We know that the eigenvectors of propagation in free space of a cylindrically symmetric system, e.g., an optical fiber, are the Bessel modes [46]. Therefore we study the eigenwaves of the TM of a finite thickness of air in a Bessel mode basis. The same experimental TM measured in the Airy spot basis is resampled using 778 Bessel functions per polarization, resulting in a total of 1556 modes. However, as in the study of the Airy basis, we focus on only the HH sub-matrix. The complex eigenvalue spectrum is plotted in Fig. 6(a), and we notice that unlike the eigenwaves in the Airy basis, the highest eigenvalues are less bunched. It is again insightful to look into the corresponding eigenstates. For the zero-thickness case, the TM in the Bessel basis is almost perfectly diagonal. Hence, the eigenvectors must be Bessel modes as a consequence of the TM construction procedure. However, with increasing thickness, the TMs are not solely diagonal anymore and this implies that the Bessel modes start mixing as a result of truncation at the outer radius .. The intensity of the experimental and simulated eigenvectors associated to high eigenvalues for an air slab of 5 µm are depicted in Fig. 6(b-l). It is noteworthy that the TM in the Bessel basis is constructed from the same data set as the TM in the Airy basis. The numerical eigenvectors in the Bessel basis are, unsurprisingly, Bessel modes up to a small amount of mode mixing due to truncation. It is clear that the experimental ones resemble the simulations but are noisy and distorted. The higher order modes (not shown) are affected by mode mixing and loss Fig. 7. (a-f) Experimental and (g-l) simulated complex eigenvalue spectra of the HH sub-matrix of the air TM in the Bessel mode basis. The air thicknesses range from 0-5 µm in steps of 1 µm. All phases are made to start at 0.
(low eigenvalues) and can't be recognized anymore. Nonetheless, the main contrast with the Airy basis is that we don't observe any IGM. This is crucial because the real principal modes of air are the (non-truncated) Bessel modes, and the TM in this basis does provide these states. The other notable feature is that the eigenmodes here fill the entire measurement area and are centered whereas those in the spot basis occupy only a fraction of the area and are located to one side. For completeness, we study the evolution of the eigenvalue spiral as a function of the air layer thickness, with the spectra for thicknesses up to 5 µm plotted in Fig. 7. We observe that the phases are initially zero (Fig. 3) and then gradually spiral inward as the modes propagate through air. This becomes progressively more apparent as the thickness of the air slab increases. Moreover, the center of the spiral is filled more uniformly for larger thicknesses. We compare these spirals with simulated ones without noise, plotted below in the same figure. As in the experiment, the eigenvalues are initially real and then spiral inward gradually with increasing thickness. The eigenvalues at the inner end of the spiral, corresponding to the higher order modes, are noisier than those at the start of the spiral. This confirms that the high order eigenvectors living primarily on the boundaries of the defined area suffer from edge effects and losses. We note that the simulated eigenvectors are similar in both bases, but not identical due to the non-orthogonality of the spot basis [6].
Nonetheless, it is still important to note that the TM in the Airy spot basis is useful. First, it is more straightforward to obtain and yet does not yield wildly deviating statistics. Second, it accounts for system imperfections and should consequently not be discarded. It is in fact meritorious to use such a TM and its Tikhonov inverse for imaging and projecting fields through the same system that has been used for the TM measurement because in this case aberrations cancel [59].

Conclusion
We demonstrated a procedure for simultaneous resampling and aberration correction of a measured transmission matrix, and tested the procedure on transmission measurements of the simplest possible sample, empty space. Resampling to an orthogonal set of Bessel modes eliminates the occurrence of spurious high singular values that could otherwise distort any measurements of transmission channel statistics. Correction of aberrations also allows one to retrieve the correct eigenvectors and eigenvalues of the transmission matrix.
The transmission eigenvectors are relevant because they correspond to the Fox-Li modes of an open confocal cavity containing the same medium. Without aberration corrections, we found that the transmission eigenvectors of a slightly aberrated TM are the Ince-Gaussian modes and they possess eigenvalues that spiral inwards on the complex plane. The spiral shape was explained by the different loss experienced by individual modes of the system. We observed that the spiral decays more rapidly towards zero for greater propagation distances in air, effectively making the air thickness a suitable parameter to filter out higher order modes.
In the aberration-corrected Bessel basis we show that not only the singular values are more faithfully represented, but also the eigenvalues and eigenvectors of the measured matrix which indeed correspond to the expected Fox-Li modes. Our method therefore is highly relevant for the interpretation of statistics such as singular value histograms and eigendecompositions of transmission matrices of simple and complex systems.

Funding
Netherlands Organization for Scientific Research NWO (Vici 68047618).