## Abstract

Wavefront shaping allows for ultimate control of light propagation in multiple-scattering media by adaptive manipulation of incident waves. We shine two separate wavefront-shaped beams on a layer of dry white paint to create two enhanced output spots of equal intensity. We experimentally confirm by interference measurements that the output spots are almost correlated like the two outputs of an ideal balanced beam splitter. The observed deviations from the phase behavior of an ideal beam splitter are analyzed with a transmission matrix model. Our experiments demonstrate that wavefront shaping in multiple-scattering media can be used to approximate the functionality of linear optical devices with multiple inputs and outputs.

© 2014 Optical Society of America

## 1. Introduction

Linear optical components like lenses, mirrors, polarizers and wave plates are the essential building blocks of optical experiments and applications [1, 2]. One has often little freedom to modify the linear optical circuit, besides rearranging components or including adaptive elements. For computational applications and on-chip light processing [3–6], it would be fantastic to have a programmable linear optical circuit that can be controlled during operation. Spatial light modulation in combination with random light scattering provides an excellent platform to accomplish this [7]. The collective interference of scattered light propagating through a scattering material results into a speckle pattern. In wavefront shaping [8, 9] incident light is modulated by a spatial light modulator to obtain a desired speckle pattern for functionality. Although one has to tolerate some losses, this technique makes it possible to transform opaque scattering media in linear optical elements that are flexible in performance [10–15]. If one controls light propagation inside strongly scattering media, one can make the most exotic linear optical circuits within a fraction of a mm^{3}.

In earlier wavefront-shaping experiments, multiple target spots have been simultaneously optimized with a single incident wavefront [8, 15, 16]. These experiments essentially demonstrate 1 × *m* linear optical circuits with 1 incident mode projected to *m* output modes. If one is capable to manipulate *n* incident modes with wavefront shaping, it becomes possible to program *n* × *m* optical circuits with a desired transmission matrix **T**. To our knowledge, no experiment demonstrating this capability has been reported.

In this article we describe an experiment in which a layer of white paint is used as an optical beam splitter. We apply wavefront shaping to intensity enhance 2 output spots for 2 separate incident modes, forming a 2 × 2 optical circuit. The enhanced spots are almost correlated like the outputs of a 50:50 beam splitter. This behavior is verified by an optical interference experiment. Our measurements indicate that the system approaches the behavior of an ideal beam splitter when the intensity enhancement *η* increases, which is confirmed by simulations. This surprisingly suggests that one can program beam splitters without prior knowledge of the transmission matrix **T**. We explain this with random matrix calculations. Our experiment demonstrates that wavefront shaping can be extended to multiple incident modes that can interfere in a controlled manner, allowing for programmable linear optical circuits.

## 2. Interference on a beam splitter

The scattering matrix of a lossless beam splitter represents a unitary transformation that can be written in its most general form as the product of three matrices [17]:

*θ*applied by the beam splitter on the incident and outgoing modes respectively. In this article we use the term input mode for an incident wave that describes a single orthogonal input of the normal beam splitter or wavefront-shaped beam splitter. The phase angle Θ in the center matrix determines the splitting ratio, which has to be Θ = (2

*n*+ 1)

*π*/2 for a 50:50 beam splitter, with

*n*an integer value.

Now consider the experiment in Fig. 1(a) in which two modes (1, 2), carrying coherent light of equal power, frequency, polarization, and wavefront are incident on a 50:50 beam splitter giving two output modes (1′, 2′). The relative phase difference Δ*θ* between mode 1 and 2 can be controlled. The power *P*_{1′}, *P*_{2′} as a function of Δ*θ* is shown in Fig. 1(b). *P*_{1′} and *P*_{2′} will oscillate as a function of Δ*θ* with a phase difference of the power oscillation *δ* = *π* for the ideal lossless beam splitter, independent of phases Δ*θ* and Ψ in Eq. (1). For this graph we have set Ψ = 0. Any nonzero value for Ψ will provide a phase offset to both the output modes, essentially shifting both *P*_{1′} and *P*_{2′} along the horizontal axes by the same amount. If one of the incident modes is blocked, a constant power is detected in both output modes that is 4 times lower than the maximum power in one of the output modes when both inputs are present. Note that energy is conserved: *P*_{1′} + *P*_{2′} = *P*_{1} + *P*_{2}, and a fringe visibility of 100% is observed.

Modeling lossy beam splitters is an interesting subject on its own with a wide variety of approaches [18–22]. We describe the effective transmission matrix **T** as a non-unitary version of Eq. (1). The phase difference *δ* = *π* will still hold if there are losses in any of the input or output modes. In such a case, the damping can be modeled in the left or right matrix in Eq. (1) by making them non-unitary with a determinant smaller than 1. This is valid for the typical beam splitter one uses. However, when the scattering in the beam splitter does not conserve energy (*i.e.* the damping is in the center matrix in Eq. (1) dictating that not all energy goes to the considered output modes), *δ* could in principle take on any value as long as the total output power does not exceed the total input power. To model the most general form of a lossy beam splitter we make the following assumptions: 1. The system is described by an effective transmission matrix **T** consisting of 2×2 elements; 2. Transmission matrix **T** provides an equal power splitting ratio; 3. The power losses for both input modes are identical.

The system behaves as a lossy beam splitter with equal splitting ratio. For convenience to compare to the ideal beam splitter, we consider that a single input mode is reflected with power reflectance |*r*|^{2} and transmitted with transmission |*t*|^{2}. For a lossy balanced beam splitter this means |*r*| = |*t*| and |*r*|^{2} + |*t*|^{2} ≤ 1. We relate the complex input amplitudes *A*_{1} and *A*_{2} to the output amplitudes *A*_{1′} and *A*_{2′} as *A*_{1} → |*r*|(*A*_{1′} + *e ^{i}^{δ}*

*A*

_{2′}) and

*A*

_{2}→ |

*r*|(

*A*

_{1′}+

*A*

_{2′}), with relative phase term

*δ*. Now we set |

*r*|

^{2}= 1/

*N*, with splitting factor

*N*and

*N*≥ 2. This leads to the following transmission matrix:

**T**is only unitary when

*N*= 2, with for example

*δ*=

*π*, and would always result to interference as shown in Fig. 1(b). The eigenvalues

*λ*

_{1}and

*λ*

_{2}of the matrix in Eq. (2) are:

*δ*. Since

**T**is a square matrix, from the singular values of Eq. (2) ${\tau}_{1}^{2}={\left|{\lambda}_{1}\right|}^{2}=\left(2+2\text{cos}\left(\frac{\delta}{2}\right)\right)/N$ and ${\tau}_{2}^{2}={\left|{\lambda}_{2}\right|}^{2}=\left(2-2\text{cos}\left(\frac{\delta}{2}\right)\right)/N$ we obtain: If this relation is fulfilled, it is guaranteed that |

*r*|

^{2}= |

*t*|

^{2}. In addition

*τ*

_{1},

*τ*

_{2}≤ 1 to guarantee a transmission not exceeding 1. This restricts the possible

*δ*to the range 2cos

^{−1}(

*N*/2 − 1) ≤

*δ*≤ 2cos

^{−1}(1 −

*N*/2) for 2 ≤

*N*≤ 4, as marked by the gray area in Fig. 2.

With wavefront shaping it is possible to use a multiple-scattering material as a balanced beam splitter. In our experiment we are working with *N* ≫ 10^{2} and therefore any *δ* is allowed. The scattering statistics of the sample, such as the the singular value distribution, and the intensity enhancement defining *N* in Eq. (4), determine the combination of *τ*_{1} and *τ*_{2} that satisfy Eq. (4). Therefore one would not expect in general a constant probability distribution for *δ* in the gray marked area in Fig. 2. We would like to approximate the behaviour of a beam splitter where *δ* → *π* since this mimics the beam splitter one normally uses.

## 3. Optimization algorithm

The optimization procedure for wavefront-shaped balanced beam splitters is illustrated in Fig. 3. We start with a single incident mode. All segments of a phase-only spatial light modulator (SLM) are subsequently addressed. The phase *ϕ _{n}* ⊂ [0, 2

*π*) of the

*n*

^{th}segment is randomly chosen and the output powers

*P*

_{1′}and

*P*

_{2′}are monitored.

*ϕ*is accepted and kept on the SLM if the summed output powers of both spots has increased and the difference power has decreased:

_{n}*P*_{1′,new}+*P*_{2′,new}>*P*_{1′,old}+*P*_{2′,old}+*ε*_{1}- |
*P*_{1′,new}−*P*_{2′,new}| < |*P*_{1′,old}−*P*_{2′,old}| +*ε*_{2}

*ε*

_{1},

*ε*

_{2}→ 0 to compensate for noise. Otherwise the previous

*ϕ*was restored. Next the (

_{n}*n*+ 1)

^{th}segment is addressed, etc. After the final segment has been addressed, the entire optimization is repeated until the desired convergence is reached.

The same procedure is repeated for the second incident mode for the same two target spots. Finally both modes are incident to perform an interference experiment as described in Fig. 1, where the relative phase Δ*θ* is controlled with the SLM.

We have decided to select this optimization algorithm because of ease of implementation and the guaranteed convergence to spots of equal power. There are algorithms that work faster and more efficiently resulting in spots of higher intensity enhancement [23]. At any rate, in the next sections we demonstrate that our algorithm is adequate for this experiment.

## 4. Experimental setup

The setup is illustrated in Fig. 4(a). The light source is a mode-locked Ti:Sapphire laser (Spectra-Physics, Tsunami) emitting transform-limited pulses at a repetition rate of 80 MHz with a pulse width of approximately 0.3 ps and a center wavelength of 790.0 nm. The pulses are spectrally filtered by a Fabry-Perot cavity with a linewidth of 1.5 nm. The beam is split and coupled into two separate single-mode fibers. The output modes have identical polarization and waist and form the input modes 1 and 2. The two modes are phase-modulated with a SLM (Hamamatsu, LCOS-SLM). The two modes are spatially overlapped with a half-wave plate (HWP) and polarizing beam splitter (PBS) cube, resulting in co-linear propagation of modes with orthogonal polarization. This allows us to completely fill the aperture of the objective (NA=0.95, Zeiss) that focuses the light on a layer of white paint. Both pulses arrive simultaneously at the sample to within 20 fs. We make sure that the power of both incident modes on the objective are identical (power of approximately 0.5 mW per mode). The layer of white paint consists of ZnO powder with a scattering mean free path of 0.7 ± 0.2 *μ*m. The layer is approximately 30 *μ*m thick and spray painted on a glass microscope slide of 1 mm thickness. The transmitted speckle pattern is collected with a second objective (NA=0.55, Nikon) and imaged on a CCD camera after reflection on a PBS, see for example Fig. 4(b). The intensity values for the CCD pixels that correspond to the target spots are integrated to obtain the output powers for the enhanced spots. The optimized spots can be transmitted through the PBS, towards a different part of the setup for applications, by rotating the HWP.

The SLM was divided into segments of 20×20 pixels. Each segment was sequentially addressed with a random phase as described in the optimization algorithm. The phase is applied by writing a gray value between 0 and 255. This corresponds to a phase modulation depth of (2.0±0.1)*π* rad. This algorithm was repeated approximately 15 times for all segments to obtain two enhanced spots of equal power at 1′ and 2′, see Fig. 4(b). The total optimization procedure for both incident modes takes about 3 hours. In our experiments we achieved *η* ∼ 5. At ≈ 10^{3} speckles this is not yet close to the theoretical maximum [8] and was limited by a noisy pulsed laser and the not completely optimized SLM segment size. It nevertheless served its purpose and we confirm interference between the output modes by adding a phase offset to the encircled area in Fig. 4(c).

## 5. Experimental results

Figure 5 presents the main result of this article: optical interference with a wavefront-shaped nearly balanced beam splitter. Figure 5(a) shows a series of camera images when on incident mode 1 a phase offset is applied. The number in the left top corners represent the phase offset Δ*θ* in gray values. All pictures are subsequently taken with the same integration time. The intensity clearly oscillates between the two target spots.

Figure 5(b) shows the output power *P*_{1′} (red squares) and *P*_{2′} (blue diamonds) as a function of the applied phase difference Δ*θ*. Both curves show sinusoidal behavior and are approximately out-of-phase, mimicking the behavior of an ideal beam splitter as is shown in Fig. 1(b). We expect an error in the phase of about Δ(Δ*θ*) = 0.1 rad due to interferometric stability during data collection and an additional systematic error of 0.1 rad due to phase calibration (both not shown). We have fitted two functions of the form *A*sin(Δ*θ* + *b*) + *c* to the measured power, which is in good agreement with the data points. From *b* we determined the phase difference |*δ*| = 2.30 ± 0.14 rad, close to but significantly different from the value of *δ* = *π* of an ideal beam splitter. Note that since the ordering of the two spots is arbitrary, we only consider the reduced |*δ*| when discussing the experimental results. Since the Both *P*_{1′} and *P*_{2′} show a fringe visibility of approximately 100%, which indicates a near-perfect mode matching between the output modes for the two seperate incident modes. The maximum measured power in both spots is approximately the same to within 5%. When one of the incident modes is blocked in the interference experiment, the output power is approximately constant (white diamonds and squares). The small spatial separation between mode 1 and 2 on the SLM gives a small crosstalk, causing fluctuations within 10%. The output power is approximately 4 times lower than the maximum power in one output mode when both input modes are incident, in excellent agreement with Fig. 1(b).

We have repeated this interference experiment 5 times for different target spots and determined |*δ*|. The result is shown in Fig. 6. All measurements were performed under comparable circumstances. Although the number of measurements are not sufficient for any statistical relevant conclusion, our measurements suggest a tendency for |*δ*| to cluster close to *π*. This is somewhat unexpected and therefore interesting to explore. In the next section we therefore present a model that predicts this behavior.

## 6. Model for the phase difference

In this section we model the observed phase difference *δ* in the interference experiment based on random matrix theory. Light undergoes isotropic multiple scattering in a sample with a thickness much larger than the scattering mean free path *l* and *kl* ≫ 1. We therefore expect the transmission matrix **T** to follow the statistics of a random matrix, as was demonstrated experimentally for ZnO [16]. One could argue that **T** is a subset of the scattering matrix **S**, and therefore *δ* can take in principle any value between [0, 2*π*) with equal probability. This would naively result in a constant probability distribution for *δ*, which is not observed.

The scattering matrix **S** has to be unitary, which sets restrictions on the allowed values for each element *s _{a,b}*. Consider a random

**S**in a basis where one input mode is one element of the input vector and one target output spot is one element of the output vector. If

**S**contains a beam splitter of equal splitting ratio, there have to be 2 rows and 2 columns in

**S**with corner elements of approximately the same amplitude:

*i*,

*j*,

*m*,

*n*} ≤ Dim(

**S**) positive integers and

*c*

_{1}≈

*c*

_{2}≈

*c*

_{3}≈ 1. For a random scattering matrix with dimension Dim(

**S**) = 2, the only possibility for a balanced beam splitter is that

*δ*=

*π*. From matrix algebra it follows that for Dim(

**S**) = 3,

*δ*can only lie on the boundary lines of the gray area of Fig. 2:

*δ*= 2cos

^{−1}(1/(2|

*r*|

^{2}) − 1) or, equivalently by reordering the two outputs,

*δ*= 2cos

^{−1}(1 − 1/(2|

*r*|

^{2})), with |

*r*|

^{2}the intensity reflectivity as defined in section 2. The amplitude coefficients of

**S**should satisfy |

*s*|

_{a,b}^{2}≥ 1/4. For Dim(

**S**) ≥ 4 any phase becomes accessible within the gray marked area of Fig. 2. However, the corresponding phase distribution is strongly dependent on Dim(

**S**), as illustrated in Fig. 7. There we have generated many random scattering matrices with different dimension that contain a balanced beam splitter [24]. The corresponding intensity enhancement is given by

*η*= |

*s*|

^{2}/〈|

*s*|

_{a,b}^{2}〉, with |

*s*|

^{2}= (1/4)(|

*s*|

_{i,j}^{2}+ |

*s*|

_{n,j}^{2}+ |

*s*|

_{i,m}^{2}+ |

*s*|

_{n,m}^{2}), and 〈|

*s*|

_{a,b}^{2}〉 = 1/Dim(

**S**). An increased probability for

*δ*=

*π*is observed with higher

*η*. The probability distribution becomes flat for small

*η*. How this scales, depends strongly on Dim(

**S**). It becomes extremely difficult to observe

*η*> 4 for 〈|

*s*|

_{a,b}^{2}〉 < 0.01, because the probability to get these realizations out of random unitary matrices becomes astronomically small. The SLM in our experiment, however, is programmed to realize exactly such situations.

Therefore we have simulated our experiment by applying our optimization algorithm on large random **S**. Figure 8 shows the observed distribution for |*δ*| as a function of intensity enhancement *η*. The main observation is that P(|*δ*|) has a global maximum at |*δ*| = *π* that increases with *η*. The simulations in Fig. 8 are performed for Dim(**S**) = 300. We control the phase of 40 input elements representing incident mode 1, and 40 input elements representing incident mode 2. Each controlled input element has a normalized input power of 1. We have set margins *ε*_{1} = 0 and *ε*_{2} = 0.001. We apply the optimization algorithm 2500 times per mode to guarantee convergence. We select output elements for which the total power of the optimization for mode 1 is within 10% of the optimization for mode 2, approximating our experiment. The intensity enhancement *η* is given by the observed power in a target spot, divided by 40 × 〈|*s _{a,b}*|

^{2}〉 = 40/Dim(

**S**), where the factor 40 comes from the number of channels that are controlled per incident mode.

It is beyond the scope of this model to match experimental conditions. In particular the large number of channels is difficult to implement. We have repeated our simulations for several Dim(**S**) and several amounts of controlled input channels, always demonstrating a global maximum at a |*δ*| = *π* that increases with *η*. This demonstrates that the two optimized spots approximate better the behavior of a balanced beam splitter with increasing enhancement, using our optimization algorithm. Based on the model with large unitary matrices, it is likely that this result is independent of the type of optimization algorithm used.

We can also explain our findings by considering an analytic model. Assume that the wavefront shaping process has been completed, where each enhanced spot has an identical intensity enhancement *η*. We are free to define a basis for the complete scattering matrix **S**, which includes the SLM and the scattering material. We write **S** as:

**S**for the beam splitter are the four top left elements. We have chosen a basis where the phase difference

*δ*of the beam splitter, as observed in our interference experiment, is included in

*s*

_{2,1}. The phase difference between the input modes

*α*

_{in}is included in

*s*

_{1,2}and

*s*

_{2,2}, the phase difference between the enhanced spots

*α*

_{out}is included in

*s*

_{2,1}and

*s*

_{2,2}. Since

**S**has to be unitary, each column of the matrix should be orthogonal to the other columns. Therefore the innerproduct between the first two columns becomes:

We define
$B={\sum}_{i=3}^{i=\text{Dim}(\mathbf{S})}{s}_{i,1}{s}_{i,2}^{*}$. For *η* ≪ Dim(**S**) we assume that all elements *s _{i,j}* still follows the statistics of the elements of a randomly generated unitary matrix, except for the elements describing the beam splitter. We are dealing with a system where Dim(

**S**) ≫ 2 and therefore we can approximate Dim(

**S**) − 2Dim(

**S**). We assume that for a random scattering matrix all elements

*s*are complex Gaussian distributed with mean 〈

_{i,j}*s*〉 = 0 and standard deviation

_{i,j}*σ*= 〈|

_{s}*s*〉 = 1. From the rules of multiplication and adding Gaussian distributions it follows that

_{i,j}*B*should be a complex Gaussian distributed value with 〈

*B*〉 = 0 and ${\sigma}_{B}=\u3008\left|B\right|\u3009=\sqrt{\frac{\text{Dim}(\mathbf{S})}{2}}$. Therefore if $\eta >\sqrt{\frac{\text{Dim}(\mathbf{S})}{2}}$,

*B*can be ignored in Eq. (7) and we expect

*δ*→

*π*to satisfy this equation. For the simulations in Fig. 8 this should be the case for

*η*= 12.2, which we indeed observe in the simulated distributions. In our experiments we have

*η*∼ 5 and Dim(

**S**) ∼ 10

^{3}, and therefore according to this model

*δ*→

*π*for

*η*∼ 10

^{1}. Therefore our experimental observations of Fig. 6 are not convincingly explained by this model. On the other hand, Eq. (7) sets restrictions on the allowed combinations of

*B*,

*α*

_{in},

*η*, and

*δ*, which we have ignored up till now. Therefore more advanced modeling is required that is outside the scope of the present paper.

## 7. Discussion

We have experimentally created two optimized spots that are correlated like the output of a lossy balanced beam splitter. The interference experiment suggests that |*δ*| → *π* with increasing *η*, which is indicated by a random matrix model. The model could not match the dimension of our experiment, nevertheless, the agreement is gratifying. Our model demonstrates that the probability distribution for *δ* depends on the number of modes in the system. It would be intriguing to perform these kind of experiments with systems of lower dimension, such as, multi-mode fibers with embedded disorder to confirm this scaling [25].

It would be fascinating to measure the transmission matrix of the sample prior to optimization. This allows for an experimental study on the influence of the optimization algorithm and *η* on the observed *δ*. In addition this would also allow to address speckle patterns with more complicated correlations.

The loss in our experiment can be reduced by several orders of magnitude by implementing more efficient wavefront-shaping procedures and using continuous-wave lasers. On the other hand, our pulsed experiment reveals the opportunity to apply this beam splitter on incident light produced in nonlinear processes, such as entangled photon pairs or higher order quantum states produced with spontaneous parametric down-conversion [26, 27]. This makes it possible to program linear optical circuits for incident quantum states to exploit quantum correlations in disordered media [28–35].

## 8. Conclusions and outlook

In summary, we have controlled light propagation in opaque scattering media with phase modulation of the two incident light modes. We have optimized two spots that show interference like a 50:50 beam splitter. Advanced optimization algorithms make it possible to create more complicated linear circuits out of multiple-scattering media [36].

## Acknowledgments

We thank J. P. Korterik, F. B. Segerink, M. P. van Exter, T. A. W. Wolterink, J. L. Herek, I. M. Vellekoop, and W. L. Vos for discussions and support. This work was supported by the Stichting Fundamenteel Onderzoek der Materie (FOM) that is financially supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO). A.P.M. acknowledges ERC grant 279248. P.W.H.P. acknowledges NWO Vici.

## References and links

**1. **E. Hecht, *Optics* (Addison Wesley, 4, 2002).

**2. **B. E. A. Saleh and M. C. Teich, *Fundamentals of Photonics* (Wiley-Interscience, 2, 2007).

**3. **M. A. Nielsen and I. L. Chuang, *Quantum Computation and Quantum Information* (Cambridge University, 1, 2000).

**4. **E. Knill, R. Laflamme, and G. J. Milburn, “A scheme for efficient quantum computation with linear optics,” Nature **409**, 46 (2001). [CrossRef]

**5. **P. Kok, W. J. Munro, K. Nemoto, T. C. Ralph, J. P. Dowling, and G. J. Millburn, “Linear optical quantum computing with photonic qubits,” Rev. Mod. Phys. **79**, 135–174 (2007). [CrossRef]

**6. **J. L. O’Brien, A. Furusawa, and J. Vučković, “Photonic quantum technologies,” Nat. Photonics **3**, 687–695 (2009). [CrossRef]

**7. **I. Freund, “Looking through walls and around corners,” Physica A **168**, 49–65 (1990). [CrossRef]

**8. **I. M. Vellekoop and A. P. Mosk, “Focusing coherent light through opaque strongly scattering media,” Optics Lett. **32**, 2309–2311 (2007). [CrossRef]

**9. **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**, 283–292 (2012). [CrossRef]

**10. **I. M. Vellekoop, A. Lagendijk, and A. P. Mosk, “Exploiting disorder for perfect focusing,” Nat. Photonics **4**, 320–322 (2010). [CrossRef]

**11. **E. G. van Putten, D. Akbulut, J. Bertolotti, W. L. Vos, A. Lagendijk, and A. P. Mosk, “Scattering lens resolves sub-100 nm structures with visible light,” Phys. Rev. Lett. **106**, 193905 (2011). [CrossRef]

**12. **J. Aulbach, B. Gjonaj, P. M. Johnson, A. P. Mosk, and A. Lagendijk, “Control of light transmission through opaque scattering media in space and time,” Phys. Rev. Lett. **106**, 103901 (2011). [CrossRef]

**13. **O. Katz, E. Small, Y. Bromberg, and Y. Silberberg, “Focusing and compression of ultrashort pulses through scattering media,” Nat. Photonics **5**, 372–377 (2011). [CrossRef]

**14. **D. J. McCabe, A. Tajalli, D. R. Austin, P. Bondareff, I. A. Walmsley, S. Gigan, and B. Chatel, “Spatio-temporal focussing of an ultrafast pulse through a multiply scattering medium,” Nat. Comm. **2**, 447 (2011). [CrossRef]

**15. **Y. Guan, O. Katz, E. Small, J. Zhou, and Y. Silberberg, “Polarization control of multiply scattered light through random media by wavefront shaping,” Opt. Lett. **37**, 4663–4665 (2012). [CrossRef]

**16. **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**, 100601 (2010). [CrossRef]

**17. **R. A. Campos, B. E. A. Saleh, and M. C. Teich, “Quantum-mechanic lossless beam splitter: SU(2) symmetry and photon statistics,” Phys. Rev. A **40**, 1371–1384 (1989). [CrossRef]

**18. **T. Kiss, U. Herzog, and U. Leonhardt, “Compensation of losses in photodetection and in quantum-state measurements,” Phys. Rev. A **52**, 2433–2435 (1995). [CrossRef]

**19. **U. Leonhardt, *Measuring the Quantum State of Light* (Cambridge University, 1, 1997).

**20. **S. M. Barnett, J. Jeffers, A. Gatti, and R. Loudon, “Quantum optics of lossy beam splitters,” Phys. Rev. A **57**, 2134–2145 (1998). [CrossRef]

**21. **L. Knöll, S. Scheel, E. Schmidt, D. -G. Welsch, and A.V. Chizhov, “Quantum-state transformation by dispersive and absorbing four-port devices,” Phys. Rev. A **59**, 4716–4726 (1999). [CrossRef]

**22. **J. Jeffers, “Interference and the lossless lossy beam splitter,” J. Mod. Optic. **47**, 1819–1824 (2000). [CrossRef]

**23. **I. M. Vellekoop and A. P. Mosk, “Phase control algorithms for focusing light through turbid media,” Opt. Commun. **281**, 3071–3080 (2008). [CrossRef]

**24. **We use Matlab 2013 for creating random unitary matrices **S** with dimension Dim(**S**) using the following syntax: Random_Matrix=rand(N,N)+1i*rand(N,N); %generates uniformly distributed random complex numbers [**S**, r]=qr(Random Matrix); % **S** is the output of an orthogonal-triangular decomposition.

**25. **N. P. Puente, E. I. Chaikina, S. Herath, and A. Yamilov, “Fabrication, characterization and theoretical analysis of controlled disorder in the core of the optical fibers,” Appl. Optics **50**, 802–810 (2011). [CrossRef]

**26. **S. R. Huisman, N. Jain, S. A. Babichev, F. Vewinger, A. N. Zhang, S. H. Youn, and A. I. Lvovsky, “Instant single-photon Fock state tomography,” Opt. Lett. **34**, 2739–2741 (2009). [CrossRef]

**27. **T. J. Huisman, S. R. Huisman, A. P. Mosk, and P. W. H. Pinkse, “Controlling single-photon Fock-state propagation through opaque scattering materials,” Appl. Phys. B. (2013). [CrossRef]

**28. **P. Lodahl, A. P. Mosk, and A. Lagendijk, “Spatial quantum correlations in multiple scattered light,” Phys. Rev. Lett. **95**, 173901 (2005). [CrossRef]

**29. **S. Smolka, A. Huck, U. L. Andersen, A. Lagendijk, and P. Lodahl, “Observation of spatial quantum correlations induced by multiple scattering of nonclassical light,” Phys. Rev. Lett. **102**, 193901 (2009). [CrossRef]

**30. **Y. Bromberg, Y. Lahini, R. Morandotti, and Y. Silberberg, “Quantum and classical correlations in waveguide lattices,” Phys. Rev. Lett. **102**, 253904 (2009). [CrossRef]

**31. **A. Peruzzo, M. Lobino, J. C. F. Matthews, N. Matsuda, A. Politi, K. Poulios, X. -Q. Zhou, Y. Lahini, N. Ismail, K. Wörhoff, Y. Bromberg, Y. Silberberg, M. G. Thompson, and J. L. O’Brien, “Quantum walks of correlated photons,” Science **329**, 1500–1503 (2010). [CrossRef]

**32. **Y. Lahini, Y. Bromberg, D. N. Christodoulides, and Y. Silberberg, “Quantum correlations in Anderson localization of indistinguishable particles,” Phys. Rev. Lett. **105**, 163905 (2010). [CrossRef]

**33. **J. R. Ott, N. A. Mortensen, and P. Lodahl, “Quantum interference and entanglement induced by multiple scattering of light,” Phys. Rev. Lett. **105**, 090501 (2010). [CrossRef]

**34. **W. H. Peeters, J. J. D. Moerman, and M. P. van Exter, “Observation of two-photon speckle patterns,” Phys. Rev. Lett. **104**, 173601 (2010). [CrossRef]

**35. **D. Bonneau, M. Lobino, P. Jiang, C. M. Natarajan, M. G. Tanner, R. H. Hadfield, S. N. Dorenbos, V. Zwiller, M. G. Thompson, and J. L. O’Brien, “Fast path and polarization manipulation of telecom wavelength single photons in lithium niobate waveguide devices,” Phys. Rev. Lett. **108**, 053601 (2012). [CrossRef]

**36. **S. R. Huisman, Light Control with Ordered and Disordered Nanophotonic Media (PhD. thesis, University of Twente, 2013).