Finding the field transfer matrix of scattering media.

When illuminated by temporally coherent light, multiply scattering media produce speckle patterns that in many situations are unpolarized on spatial averaging. As a result, the underlying field statistics are assumed to be Gaussian and information about them can be extracted from intensity-intensity correlations. However, such an approach cannot be applied to any scattering medium where the interaction leads to partially developed speckle patterns. We present a general procedure to directly measure the field transfer matrix of a linear medium without regard to the scattering regime. Experimental results demonstrate the ability of our procedure to correctly measure field transfer matrices and use them to recover the polarization state of incident illumination.


Introduction
Multiple scattering is typically considered to degrade the information in a beam propagating through a random medium. For temporally coherent radiation propagating through an optically thick medium, the light scattered out of the medium will produce an interference pattern with alternating bright and dark regions known as speckles. A speckle pattern having a negative exponential intensity distribution is called fully developed and obeys Gaussian field statistics [1]. Since the field distribution is a Gaussian random variable, it was thought to contain neither information about the underlying scattering medium nor information about the beam incident on the medium.
While an individual speckle pattern does not contain a large amount of useful information, the correlation between speckle patterns can show properties of the underlying medium if the illumination is held constant or of the incident field if the scattering medium is stable in time [2,3]. The former idea was pursued first, and many techniques for speckle photography and speckle interferometry were developed to measure the stress deformations, vibrational modes, and other mechanical properties of materials under constant illumination.
Later, it was recognized that speckle patterns could track changes in the direction of an incident beam; the direction information was not completely lost. Instead, the direction change was encoded into the speckle pattern [4,5]. Techniques have been developed that recover the field information through correlating speckle patterns generated by a fixed material when illuminated by known and unknown fields [2]. It was suggested that performing speckle correlations over a fixed realization of a random medium could offer the possibility to do many manipulations that normally required precise optical instruments, and proposals were made to use random media as a lens, a spectrum analyzer, or many other optical devices [6].
It is interesting to note that, in many respects, the speckle correlations in a stationary pattern are complimentary to the measurements in the classical Hanbury Brown and Twiss (HBT) experiment [7]. In the HBT experiment, the correlation function is mapped out very well within a single speckle by comparing the temporal fluctuations in the intensities recorded at different points in space. The speckle correlation techniques compare the spatial intensity fluctuations across a speckle pattern at given moment in time. Generally, intensity correlation measurements, like HBT and speckle correlation, are used to characterize properties of the electric field, so they rely on the speckle pattern being fully developed in order to relate the field correlation to the measured intensity correlation; they are based on an assumption regarding the statistics of the field.
In this paper, we suggest a different approach to recovering information from a speckle pattern that does not rely on ensemble averaging or assumptions about the underlying field statistics. We propose to actually measure components of the field transfer matrix across a particular output plane. In a proof of concept demonstration, we show successful recovery of the polarization state of an unknown beam from the speckle pattern it produces. Our technique is a deterministic procedure performed simultaneously in many different spatial locations. Because the technique does not rely on any assumptions about the statistics of the field distribution and because it is carried out by point operations rather than image correlations, it can be used in regimes ranging from no scattering to high-order multiple scattering.

Theory of field transfer
Consider a linear, multiply scattering medium illuminated by polarized light. The field will undergo many series of scattering events before emerging from the medium. In a particular series, the field after the first scattering event at r 1 that illuminates the second scatterer at position r 2 may be expressed as where α is a complex tensor that expresses the magnitude and phase of the coupling between incident and scattered field components. It also includes transferring the scattered field to the point r 2 . Now, replacing E inc with E scat,1 and ) r , (r , we obtain the field at position r 3 on the third scatterer. This process can be continued until the field reaches the detector at position r to obtain the contribution of a particular scattering path to the field at the detector: In this particular decomposition the scattering matrix that relates the input and output fields contains information about the propagation from r i to r i+1 as well as information about the scattering event at r i . Since we are interested in the transfer of the field through the medium, not in the particular details of how the transfer occurs, the mixing of scattering and propagation information is not important. Furthermore, there are many scattering paths connecting the illuminated points on the front surface of a multiply scattering medium to each detection point after the medium as shown in Fig. 1. If the scattering path lengths are shorter than the coherence length of the illuminating radiation, the scattering process will not depolarize the output light; it will only change the radiation's state of polarization. The resulting field at any detection point due to a particular input point is simply the coherent summation of the outputs of each of the different paths starting at the illuminated point on the front surface and ending at the detection point. Also shown in Fig. 1 are three random walks on the Poincare sphere that begin in the same polarization state (chosen here to be circular). However, since they interact with different configurations of the scattering medium, they arrive at the detector in different states of polarization, denoted by the large dots of the appropriate color. The total field at the detector is the coherent sum of the different paths and is in the state of polarization marked by a white dot labeled "detected" on the Poincare sphere. Thus, the resulting field, neglecting time dependence, at a point r due to all illuminated input points can be written as The normalized input field, inc ê , can be factored out of the sum so long as its polarization is spatially constant, and we can express the resulting output field in terms of a single transfer matrix, (r) α ~, which subsumes the intensity profile of the illumination. For a given experimental geometry and illumination source, (r) α ~ is only a function of detector location; however, if the spatial intensity profile of the illumination is modified, (r) α ~ will change as well because the intensity profile acts as a weighting function for the contribution of each scattering path. We emphasize that the tilded quantities represent the only measurable parameters of the scattering material because it is not practically possible to separate the contributions of individual paths from the detected intensity. It is also important to note that the illuminated points can have an arbitrary spatial extent and spatial intensity profile on the random scattering medium as long as the scattering paths remain coherent with one another. The problem of determining the elements of the transfer matrix to a given point in the detection plane does not depend on the precise nature of the process that produced the output field and thus applies to all scattering regimes. In general, α ~ has nine elements with eighteen unknowns: nine coupling magnitudes and nine phases. Because the scattering is not isotropic, α ~d epends on the direction of propagation of the incident light, and it is difficult to measure the full transfer matrix for an arbitrary geometry. However, if the scattering medium is surrounded by an isotropic medium, the electric field of the illumination is confined to a plane and can be decomposed into two orthogonal polarization states with a phase between them. If the scattered fields are allowed to propagate away from the scattering medium before detection, they can also be decomposed into two orthogonal polarization states and a phase term. In this situation, the transfer matrix has only four elements consisting of eight unknowns, which can be determined by illuminating the scattering medium with appropriately polarized light. Moreover, we can choose one of the elements of the transfer matrix to be real since we cannot measure absolute phase at a point and only compare the intensities between points. Additional simplifications can be introduced by realizing that it is not necessary to characterize all seven of the remaining unknowns simultaneously. If a polarizer oriented along the x-axis is placed between the scattering medium and the detector, 21 α = 22 α =0, and there are only three unknowns that need to be characterized. The detected intensity at point r is then given by whereθ is the phase between the x and y components of the incident field, and E x and E y are their respective magnitudes. In Eq. (4), φ is the phase introduced by the coupling of E y into a scattered x polarized field, and 11 α and 12 α are the magnitudes of the coupling of the incident x and y polarized fields, respectively, into scattered fields polarized along x. The polarization of the scattered field and the elements of α ~ measured are determined by the orientation of the final polarizer.

Calibration and field recovery procedure
The procedure for recovering an unknown incident state of polarization comprises three main steps: measuring the transfer matrix for many points in the detection plane, selecting unique field combinations or transfer matrices, and solving Eq. (4) for the incident field components using the transfer matrices. First, the needed components of the transfer matrix corresponding to each detector point are determined via a calibration with known fields. To determine the magnitudes of the elements of the transfer matrix and eliminate the sign ambiguity in the argument of the cosine in Eq. (4), we use four calibration states. The relative spatial intensity profiles of the unknown source and the calibration source should be the same since the intensity profile of the illumination weights the contribution of each scattering path to the detected intensity. A fully polarized field is characterized by three different parameters, and, as with the transfer matrix elements, at least three independent combinations of the incident field are needed to completely determine its polarization state. For a static system, each point in the detection plane sees a particular transformation of the input field resulting from the combination of the scattering paths that end at that point. Thus, determining the parameters of a fully polarized field requires that the intensity be measured at three or more points with independent transfer matrices. As a result, the detection system must resolve at least three speckles.
As an example, the use of a particular combination of detectors is illustrated in Fig. 2. The axes of this representation of transfer matrices are defined in a manner analogous to the Poincare sphere with the elements of the transfer matrices at a particular point taking the place of the input field components that they couple [8].  Since there is no reason to select any particular group of independent detection points, it is possible to form many different groups and then perform a statistical analysis on the recovered fields rather than relying on the result of a single combination. Note that a medium with transfer matrices covering only a small portion of the sphere may not be used to fully analyze any unknown field from its speckle pattern simply because the medium does not produce a sufficient number of independent combinations of the field.
Finally, for the selected groups of transfer matrices as illustrated in Fig. 2, the unknown field parameters, E x , E y , and θ , are determined by solving the following system of equations for each group of detectors where all of the transfer matrix elements are known from the calibration process and where r i denotes the location of the point in the detection plane rather than the location of a scatterer.

Experimental demonstration
In order to demonstrate that our procedure is valid across all scattering regimes, we measured transfer matrices for the extreme cases of deterministic, single scatterers and a heavy multiple scatterer. For the single scatterers, the measured transfer matrices can be qualitatively compared to the expected values using the spherical representation from Fig. 2. On the other hand, since the multiply scattering medium produces many diverse mixings of the incident field, it is difficult to asses the capability of our procedure by direct examination of the measured transfer matrices. In this case, we will use the transfer matrices to recover the states of polarization of plane waves in order to show our procedure works even in a regime of heavy multiple scattering.
In our experiment, a laser beam with a controlled state of polarization was incident on the scattering medium and the speckle pattern resulting from its transmission through the medium was recorded by a CCD camera. A polarizer with a fixed orientation was placed in front of the detector to simplify the analysis as explained before. The speckles produced by the multiply scattering medium were approximately 5 pixels across on the CCD and each measurement sampled a few thousand speckles. The illumination was provided by a 532 nm laser beam that was passed through a polarizer to ensure a pure linear polarization and then through a half-wave plate and quarter-wave plate to generate the desired states of polarization. Figure 3 illustrates the measured transfer matrices for a polarizer, a quarter-wave plate, and a multiply scattering medium. The multiply scattering medium is a composite dielectric material with a thickness of 100μm and characterized by a transport mean free path of 10μm. Figures 3(a) and (b) show measurements from a polarizer oriented at approximately ±45° and a quarter-wave plate rotated in 15° increments from 0° to 90°. From linear optics theory, we expect a polarizer oriented at 45° and 135 o to have equal coupling strengths through the final polarizer for both incident x and y field components. Also, when the polarizer is oriented at 45°, the transmitted field components will be in phase, and when it is oriented at 135°, the transmitted field components should have a π phase difference between them. In Fig. 3(a) we see comparable coupling of the orthogonal field components, although there is a slight misalignment of the polarizer, and rotating the polarizer from 45 o to 135 o introduces a π phase shift between the elements of the transfer matrix as expected. In Fig. 3(b), the green line denotes the path that the transfer matrix of a quarter-wave plate should follow on the sphere as the wave plate is rotated through 90 o , and the labeled groups of points are the measured transfer matrices for a quarter-wave plate in the indicated orientations. In both cases the measured transfer matrix elements show behavior consistent with what was expected and demonstrate the ability of our process to analyze deterministic single scatterers. In Fig. 3(c) we present a subset of the measured transfer matrices for the multiply scattering sample. For a truly random medium, we would expect the measured transfer matrices to uniformly cover the sphere of possible matrices; however, because of numerical instabilities involving calculations with small numbers, the points in a ring around Δ axis of the sphere are rejected by our processing algorithm. Near the axis, either 12 11α α >> or 12 11α α << , and the smaller of the two can be approximated as 0 so that transfer matrices very close to the Δ axis are moved onto the axis. Also, there seems to be some clustering of the points near φ =-1. Even though the transfer matrices are not uniformly distributed on the sphere, our measurements show that the sample produces a sufficiently large number of substantially different mixings of the incident field. This is the only requirement for recovering the state of polarization of the incident field.
Since the medium exemplified in Fig. 3(c) is in a regime of multiple scattering, it is difficult to assess the accuracy of the transfer matrix measurement by viewing of the matrices using this spherical representation. However, one can still examine their accuracy by using the determined matrices to infer the polarization states of different beams illuminating the medium. Fig. 4 shows typical experimental results for a +45 degree linear polarization and an elliptical input state, represented by the blue dots. The white dots denote the polarization states recovered by different pixel groups, and the red dots represent the geometric centers of the white data points. The experimental data shown above subtend a solid angle of approximately 0.158 steradians. The spread in the states recovered by the different detector combinations can be mitigated to a large extent by averaging the results of many combinations. The averaging can be done well because transfer matrices are measured for a large number of detectors. In order to quantify the error in the recovered polarization state, we can compare the normalized Stokes vector [9] of the recovered field to that of the incident field. In Fig. 4(a) we input a Stokes vector of (0, 1, 0) and measured (0.018, 0.9997, 0.014), and in (b) we input (0, 0.643, 0.766) and measured (0.067, 0.638, 0.764). As can be seen, the recovered Stokes vector components do not deviate by more than 1% from their expected value.
We have also simulated our experiment to study the effects of the detector selection criteria and detector noise on the measured field transfer matrices. We choose 11 α and 12 α randomly distributed uniformly between 0 and 1 and 12 φ to be randomly distributed uniformly between 0 and 2π. The magnitudes of the coupling matrices were then scaled so that the resulting calculated intensity distribution was similar to the actual data for unit strength electric field inputs. Gaussian white noise with a signal to noise ratio of 34 was then added to the intensity image to simulate the detector noise in a real measurement. Speckle images were generated for both the calibration and test states and processed using the same code as the experimental data. Our simulations indicate that the most significant source of error in the data collection and processing is the noise in the detector itself. For simulated data with no noise, the solid angle covered by the measurement data on the Poincare sphere (white dots in Fig. 4) is approximately 0.013 steradians. When noise comparable to the noise of the detector used in the experiment is added to the simulated data before processing, the spread of the recovered states increases to 0.048 steradians. The remaining error is likely due to mechanical instabilities in the experiment.

Conclusions
In scattering from inhomogeneous media, complex transfer functions relate the input field to each point of the emerging, "random", electromagnetic field. We have demonstrated that the field transfer matrices of a system can be practically measured for specific geometries. In doing so, no assumptions about the statistics of the random field and no specific description of the scattering process need to be made. Previous attempts to utilize multiple scattering have relied on averaging cross-correlations of speckle patterns. Those methods rely on assumed relationships between field and intensity cross-correlations. Further, they require ensemble averaging which is usually performed over many independent realizations of the random medium. When the discussed applications concern a single realization of the randomness, the medium is assumed to be ergodic so that averaging over the image from a single realization can be done instead of averaging images from many realizations. Of course, errors are introduced depending on how accurate the ergodic assumption is and how well the averaging is done. Our method, on the other hand, makes use of the deterministic transfer of the field to each point in the detection plane, and no averaging or further assumptions are needed.
We demonstrated our technique by calculating transfer matrices for known, deterministic, scattering media and for highly scattering media. For the known samples, the measured transfer matrices behave in an expected manner. We used the random medium to solve the inverse problem of determining the polarization state of a beam illuminating the far side of the multiply scattering medium. Due to the large number of measured transfer matrices, it is possible to determine the polarization state with many different combinations of detectors and use statistical techniques to improve the quality of the measurement. The remarkable precision in recovering the incident state is also supported by numerical simulations and suggests the possibility of using calibrated random media as efficient polarimeters. Since our technique does not rely on a Gaussian distribution of the scattered field, the direct measurement of the transfer matrix of a system opens up interesting possibilities for extracting information about the underlying scattering process in regimes ranging from single to heavy multiple scattering. In addition, potential applications may be found for forward problems in areas such as interferometry and communication.