Controlling Light Transmission Through Highly Scattering Media Using Semi-Definite Programming as a Phase Retrieval Computation Method

Complex Semi-Definite Programming (SDP) is introduced as a novel approach to phase retrieval enabled control of monochromatic light transmission through highly scattering media. In a simple optical setup, a spatial light modulator is used to generate a random sequence of phase-modulated wavefronts, and the resulting intensity speckle patterns in the transmitted light are acquired on a camera. The SDP algorithm allows computation of the complex transmission matrix of the system from this sequence of intensity-only measurements, without need for a reference beam. Once the transmission matrix is determined, optimal wavefronts are computed that focus the incident beam to any position or sequence of positions on the far side of the scattering medium, without the need for any subsequent measurements or wavefront shaping iterations. The number of measurements required and the degree of enhancement of the intensity at focus is determined by the number of pixels controlled by the spatial light modulator.

The ability to focus light through highly scattering translucent or 'opaque' random media has been a long-standing challenge. Strong scattering impedes information transfer through random media and is a limiting factor in many optical imaging and characterization systems. The last decade however has witnessed a tremendous advance in techniques for focusing fields that pass through complex media through shaping of incident wavefronts. We refer the reader to Rotter and Gigan 1 which provides an excellent and timely review of the topic.
Vellekoop and colleagues 2 were the first to show that by shaping incident wavefronts with a spatial light modulator (SLM) one can compensate for scattering phenomena to produce light focused onto a point inside or beyond the scattering medium [3][4][5] . Their technique and variants thereof have been used to focus light through white paint layers, eggshells, clouds, dense fogs 6 , and biological tissue 7 . The ability to focus light at a desired location within or beyond highly scattering media has potential applications ranging from biomedical engineering and microscopy to endoscopy 8 (and references therein), optical trapping 9 , super resolution imaging 10,11 , and nano-positioning 12 .
Their enormous potential notwithstanding, widespread use of the above and related schemes in real-world applications has somewhat lagged expectations. Needless to say, adoption of any wavefront shaping method to focus light through highly scattering media hinges on the practicality of the scheme for determining the optimal wavefront. Methods developed to date by and large belong to one of two categories.
(i) Iterative methods that use measurements of transmitted field magnitudes (square roots of field intensities). The majority of wavefront shaping methods developed to date use measurements of transmitted field magnitudes to progressively improve field focus by sequentially changing the phase retardation imposed on each pixel of the input beam using an SLM. These methods maximize intensity of the output field at a defined location by changing one SLM pixel at a time 13 . While very powerful, these methods tend to be slow and prone to convergence to a local optimum. They also must be repeated every time a new target output intensity profile is specified. Oftentimes, speckle correlations can be exploited to refocus fields onto a neighboring pixel with only minor degradation in intensity, without restarting the algorithm. Recent advances have increased the depth of the medium to where the "memory effect" (i.e, correlations in the fields) can be reliably exploited 14 . In the deep medium regime, where the fields are decorrelated and the memory effect no longer holds, refocusing still requires the procedure to be repeated every time a new target output profile is specified.
(ii) Non-iterative methods that use measurements of transmitted field magnitudes and phases. Knowledge of the transmission matrix (TM) of the medium allows for the optimal wavefront to be computed non-iteratively for any desired output. If the desired output is a focal spot, then only a portion (e.g. a row or set of rows) or the TM are required. The first measurement of an optical TM was undertaken by Popoff et al. 15 . Using a common-path interferometer equipped with an SLM, the TM of an opaque layer of ZnO nanoparticles was measured and then used to generate optical foci in the output plane. This method was also used to demonstrate image transmission through the opaque layer 16 . Most approaches developed to date to measure TMs rely on holographic methods. Needless to say, methods for determining TMs that do not require a reference beam would be tremendously useful in this context. A first step in this direction recently was taken by Drémeau et al. 17 , who used a phase retrieval algorithm to measure the complex TM of a highly scattering medium using a digital micro-mirror.
This work introduces a new waveform-shaping technique for producing focused fields that uses semidefinite programming (SDP) to construct portions of the medium's TM from intensity measurements only; once the relevant parts of the TM are calculated the method proceeds like any of the above non-iterative schemes. The SDP approach leverages a rigorous yet flexible computational framework that utilizes the algorithm developed by Waldspurger et al. 18 to retrieve the phase of elements of the TM after recording intensities in the desired focal spot produced by several randomly structured illuminations of the scattering medium. We note that SDP has been successfully applied to many problems in other fields, ranging from X-ray and crystallography imaging and  A 633 nm CW laser is expanded then spatially filtered through a lens and a 100 μm pinhole (not shown). The polarization optics is placed before the SLM. The homogeneous beam is then reflected by the SLM that is 4f imaged onto the sample. The resulting speckle pattern from the scattering sample is collected by a lens and imaged on the CMOS camera. diffraction imaging to Fourier optics and microscopy 19 . We believe this work to be the first to apply this powerful algorithm to problems related to optical phase retrieval enabled wavefront control. We demonstrate that SDP allows the construction of the incident wavefronts that generate intense foci anywhere beyond the medium from a single set of measurements of intensities generated by random illuminations. Fig. 1. A highly scattering, random medium is sandwiched in between input and output apertures A and B. In what follows, we assume both apertures are rectangular and comprised of M = K × L pixels. Light in the input aperture impinging on the random medium therefore is characterized by M complex-valued electric field samples a (k,l) with 1 ≤ k ≤ K, 1 ≤ l ≤ L. Depending on the experimental setup used, we can control both the amplitude and phase of a (k,l) , or as in more common, only its phase. Likewise, light in the output aperture exiting the random medium is characterized by M complex-valued electric field samples b (k,l) with 1 ≤ k ≤ K, 1 ≤ l ≤ L. Our experimental setup only allows for measurements of the magnitude of each b (k,l) . Note: for expositional simplicity we assume that the number of pixels in the input and output apertures are equal. In an experimental setup, they often are different with the number of pixels in B exceeding those in A. The arguments below however are easily generalized to this setting. Likewise, we characterize electric fields in A and B by scalars. Again, the arguments below are easily modified to account for the vector nature of the fields. Alternatively, polarizers can be added to the input and output apertures to effectively scalarize the field. To simplify the notation, we replace all double indices (k, l) by a single index m = (k − 1)K + l and immediately "flatten" all inherently multidimensional tensors introduced below accordingly.

Preliminaries. The setup under consideration is shown in
The fields in A and B are related by the TM of the random medium. Define the M × 1 vectors For future reference, we denote where "m" denotes the row index of T while m′ denotes the column index of T. With this notation, the M × 1 vectors m m m contain the complex conjugates of the elements of row m of the transmission matrix T. Note: use of the terminology "transmission matrix" implies that the vectors a and b resolve all (propagating) modes in A and B. In practice, we often violate this condition, especially in A. Undersampling of fields in A however does not invalidate the method described below.
The objective is to sculpt the incident field so that it produces a maximally focused spot at a predefined location in the output plane subject to an energy constraint. Mathematically, our goal is to determine the incident field vector a that maximizes b m for a predefined m subject to = a 1 2 2 . We will consider two cases: (i) full amplitude and phase control of all a m and (ii) phase-only control of a m , i.e., |a m | 2 = 1/M for all m. It follows from (1) . Hence, the above objectives would be easily realized if † z m , i.e. the m th row of the transmission matrix T, was known. Unfortunately † z m is not known. Below, we first describe a procedure for determining † z m from measurements of magnitudes of transmitted fields at pixel m (Section 2.2). Expressions for the optimal a follow easily (Section 2.3).

Phase retrieval algorithm.
To determine z m , we illuminate the random medium with N randomly selected or "trial" fields a t t ( ) Note that we have abused notation in thus defining b m -we do so to keep our notation light. It follows from Eq. (1) that Or equivalently, that We therefore proceed to determine z m and u m by solving the following optimization problem: We solved the above problem using the PhaseCut algorithm described in ref. 18, restated below for expositional purposes in terms of the above notation. We start out to note that if u m were known, then the vector z m that minimizes Eq. (10) can be explicitly computed and is given by where A † is the Moore-Penrose pseudoinverse of A. Thus, following PhaseCut 18 , our strategy is to solve Eq. (10) for u m and then determine z m via Eq. (11). To that end, we note that plugging in z m in Eq. (11) into the objective function on the right hand side of Eq. (10) yields where we have utilized the fact that ( Consequently, for z m given by Eq. (11), the optimization problem involving z m and u m in Eq. (10) can be expressed as an optimization problem involving only u m given as: implies that the diagonal elements of U m will necessarily equal one. Equations (12) and (13) constitute different formulations of equivalent optimization problems that are still difficult to solve because of the non-convex rank constraint in Eq. (13). Dropping the non-convex rank constraint in Eq. (13)  Equation (14) is a convex optimization problem and can be solved efficiently using numerical solvers. We solve Eq. (14) using the cvx package 20 . If the N × N matrix U m thus obtained via Eq. (14), has rank 1 then we have exactly solved the original optimization problem in Eq. (12) because it is equivalent to Eq. (13). We can compute the N × 1 vector u opt from U m by computing the eigen-decomposition of U m and setting is the eigenvector of U m associated with its largest eigenvalue. We apply the same procedure as well when U m is not rank 1. The theoretical guaranties accompanying PhaseCut 18 establish that when N > O(MlogM), then this procedure will perfectly recover u m with extremely high probability in the noise-free setting when the columns of the matrix A are drawn at random as we have. We then obtain the desired M × 1 vector z m from u m via Eq. (13).
Having estimated z m , from N > O(MlogM) measurements, we now describe how to form a maximally focused spot at pixel m.
Focusing wavefront. We next determine the incident field vector a that maximizes z a m m for a predefined pair m subject to the energy constraint = a 1 2 2 . We first consider the case where we have full amplitude and phase control of all a m . This optimization problem can be expressed mathematically as: Equation (18) has a closed-form solution that is given by exp arg  Figure 2 shows the experimental setup. The light source is a single longitudinal mode CrystalLaser diode laser with wavelength λ ≈ 633 nm and output power 50 mW. The beam is expanded to a diameter of 20 mm by a beam collimator. A set of polarization optics is used to select the suitable polarization state of the incident beam to achieve phase only modulation. The SLM is the a phase only Holoeye PLUTO. It is a LCOS (Liquid Crystal on Silicon) micro-display with full HD resolution (1920 × 1080 pixel) and 8 μm pixel pitch. The surface of the SLM is 4f imaged and focused on the scattering sample by a lens with a 50 mm focal length. The Fourier plane of the backside of the sample is imaged onto a CMOS camera. The detector is a PhotonFocus camera series MV1-D2048 with 2048 × 2048 resolution and pixel size 8 μm.

Experimental Results
The samples used for this experiment are ground or frosted glass and yogurt. The glass diffuser is 2 mm thick and 120 grit. The yogurt sample is prepared using plain white yogurt spread in between two thin microscope slide that are pressed together to form a thin white translucent sample similar to white paint on glass. The yogurt sample thickness is L~150 μm. The transport mean free path length ( ⁎ l s ) of this sample at λ ≈ 633 is extracted from the angular width at half maximum of its coherent backscattering peak: The SLM display is subdivided into M equally sized squares that are dubbed superpixels. For each incident wavefront the superpixels' phases are randomly set between [0, 2π]. The corresponding transmitted intensities are measured. Here, ground glass is used as the scattering sample shown in Fig. 3(b). The 461 field intensities are measured at pixel k l ( , ) = (940, 540) (note: in this experiment we used a different number of pixels in A and B as alluded to in Section 2.1; in what follows, M refers to the number of pixels in A and the number of pixels in B is constant and much larger). The algorithm generates the incident field a opt that produces the maximum field intensity at the desired focus using Eq. (19). The 462-th intensity measurement displayed corresponds to that produced by the optimal wavefront a opt .
In what follows, we define the enhancement factor η m as 5 : where b m opt is the intensity at m ( ) when a opt and b m is the average intensity at m over the N training realizations. In Fig. 4(a), we show b m for values of m centered around (940, 540). Figure 4(b) displays the enhanced intensity at pixel (940, 540); this yields an enhancement factor of η = 48 m which in par the results of the ratio of the intensity at focus to the mean intensity of the speckle outside generated in ref. 15.
The SDP algorithm is stable and fast enough to be relevant for samples that are quasi-static such as a fresh yogurt sample. This sample has a speckle persistence time in the minutes. This is a compromise between the ground glass sample speckle pattern which is static for hours, and a live biological tissue which has speckle persistence time in the order of milliseconds. The measurements in Figs 5 and 6 show that for a plain yogurt sample the intensity enhancement at the desired point increases as the number of superpixels increases. Choi et al. 21 showed  that once the TM is generated, the medium can be used as an unconventional lens to achieve imaging beyond the diffraction limit of the optical system. To this end, we show in Fig. 5 that the generated focus can be confined to within a single detector pixel.
The SDP algorithm also allows for the generation of as many foci as desired. This is accomplished with the same initial set of N = MlogM random excitations used to generate a single focus. The number of rows of the TM required now exceeds one and equals the number of foci. Figure 7 shows a set of foci generated beyond the scattering medium that traces out MICHIGAN comprised of 157 foci. This method becomes a relevant alternative in the deep medium regime where fields de-correlation occur and the memory effect may no longer hold. Clearly, SDP-based wavefront shaping has the potential to be used in tandem with fast light modulators to perform fast measurements in dynamic media.

Conclusions
We demonstrate a new method for controlling light transmission through highly scattering random media by shaping wavefronts that optimally focus the incident beam to any position, or sequence of positions, on the far side of a scattering medium. The optimal wavefront for a particular focus position is determined in  closed-form from knowledge of the complex-valued transmission matrix, which we determine using a sequence of intensity-only measurements, without the need for a reference beam, using a semidefinite programming based phase retrieval method. Once the transmission matrix is thus determined, we demonstrate that increasing the number of modes increases the intensity of the focus and that we can steer the focus.