Programming balanced optical beam splitters in white paint

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 speckle spots of equal intensity. We experimentally confirm by interference measurements that the output speckle 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.


I. 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][4][5][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][11][12][13][14][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 speckle 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 speckle spots for 2 separate incident modes, forming a 2 × 2 optical circuit. The enhanced speckle 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 measurement 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. shift ∆θ applied to incident mode 1. All powers are normalized to the sum of the two input powers

II. 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]: The nonzero terms in the first and last matrix are phase differences Ψ and Φ 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 Θ = (2n + 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; this is the same as controlling Φ = ∆θ in Eq. (1). 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 δ = π 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][19][20][21][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: Assumption 1 The system is described by an effective transmission matrix T consisting of 2 × 2 elements.
Assumption 2 Transmission matrix T provides an equal power splitting ratio.

Assumption 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 input amplitudes A 1 and A 2 to the output amplitudes A 1 and A 2 as A 1 → |r|(A 1 + e iφ 1 A 2 ) and A 2 → |r|(e iφ 2 A 1 + A 2 ), with phase terms φ 1 and φ 2 . Now we set |r| 2 = 1/N , with splitting factor N and N ≥ 2. This leads to the following transmission matrix: For the ideal lossless 50:50 beam splitter, T is only unitary when N = 2, with for example φ 1 = φ 2 = π/2, and would always result to interference as shown in Fig. 1(b). The eigenvalues λ 1 and λ 2 of the matrix in Eq. (2) are: Next we make the substitution φ = (φ 1 + φ 2 )/2. The observed phase difference in an interference experiment is given by δ = φ 1 + φ 2 = 2φ. Since T is a square matrix, from the singular values of Eq. (2) τ 2 1 = |λ 1 | 2 = (2 + 2 cos(φ))/N and τ 2 2 = |λ 2 | 2 = (2 − 2 cos(φ))/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 be in the range between With wavefront shaping it is possible to use a multiple-scattering material as a balanced beam splitter, which is inherently lossy. In our experiment we are working with N 10 2 and therefore any φ and δ are 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. on the SLM, using both optimized phase patterns. A relative phase differenced ∆θ between mode 1 and 2 is applied with the SLM to confirm optical interference like in Fig. 1(b).

III. 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. φ n is accepted and kept on the SLM if the summed output powers of both spots has increased and the difference power has decreased: with positive tolerances 1 , 2 → 0 to compensate for noise. Otherwise the previous φ n was restored. Next the (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 speckles of higher intensity enhancement [23]. At any rate, in the next sections we demonstrate that our algorithm is adequate for this experiment.

IV. EXPERIMENTAL SETUP
The setup is illustrated in Fig. 4(a). The light source is a mode-locked Ti:Sapphire laser The SLM was divided into segments of 20x20 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 speckles of equal power at 1 and 2 , see Fig. 4(c). The total optimization procedure for both incident modes takes about 3 hours. We confirm interference between the output modes by adding a phase offset to the encircled area in Fig   4(b).  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. 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 speckle 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 π. In the next section we present a model that predicts this behavior.

VI. 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 φ 1 and φ 2 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 speckle 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: |s i,j | = c 1 |s n,j | = c 2 |s i,m | = c 3 |s n,m | with {i, j, m, n} ≤ N positive integers and c 1 ≈ c 2 ≈ c 3 . 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 and 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. [36] 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. Our experiment, however, is setup to achieve exactly such realizations by optimization.
Therefore we have simulated our experiment by applying our optimization algorithm on large random S. 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 speckle 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 speckle 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: where the first two elements of the input vector correspond with the field in the input modes and the first two elements of the output vector correspond to the field in the output modes. Therefore the relevant elements in 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 speckle 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 = i=N i=3 s i,1 s * i,2 . For η N 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 N 2 and therefore we can approximate N − 2 ≈ N . We assume that for a random scattering matrix all elements s i,j are complex Gaussian distributed with mean s i,j = 0 and standard deviation σ s = |s i,j | = 1. From the rules of multiplication and adding Gaussian distributions it follows that B should be a complex Gaussian distributed value with B = 0 and σ B = |B| = N 2 . Therefore if η > N 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 N ∼ 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.

VII. DISCUSSION
We have experimentally created two optimized speckle 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 [24].
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 [25,26]. This makes it possible to program linear optical circuits for incident quantum states to exploit quantum correlations in disordered media [27][28][29][30][31][32][33][34].

VIII. 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 speckles 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 [35].