Blind focusing through strongly scattering media using wavefront shaping with nonlinear feedback

Scattering prevents light from being focused in turbid media. The effect of scattering can be negated through wavefront shaping techniques when a localized form of feedback is available. Even in the absence of such a localized reporter, wavefront shaping can blindly form a single diffraction-limited focus when the feedback response is nonlinear. We developed and experimentally validated a model that accurately describes the statistics of this blind focusing process. We show that maximizing the nonlinear feedback signal does not always result in the formation of a focus. Using our model, we can calculate the minimal requirements to blindly focus light through strongly scattering media.


Introduction
Refractive index inhomogeneities in a turbid medium scatter light in a complex manner. Consequently, a focus inside these types of media becomes more aberrated with increasing depth, until ultimately no ballistic light is left and the focus decays into a random speckle pattern. Even at this depth, light can still be focused through wavefront shaping techniques [1,2]. These techniques have been used to focus light inside scattering media for various applications, such as optical manipulation [3], optogenetics [4] and fluorescence microscopy through an intact skull of a mouse [5]. Using wavefront shaping techniques, light has been focused through several centimeters of biological tissue [6].
The main limitation of these wavefront shaping techniques is the requirement of a localized reporter (e.g. guide star), providing feedback of the light intensity at the focus location. Many different mechanisms can act as localized reporters, such as point detectors, fluorescent or nonlinear markers, acoustically tagged light and photoacoustic absorbers [2]. However, a localized reporter is not always available. For instance in multiphoton fluorescence excitation microscopy, generally all structures of interest are stained with fluorescent markers, which will all generate a signal when illuminated. When the excitation light is spread over a larger area, the total feedback signal will originate from multiple indistinguishable reporters. However, when a weak ballistic unscattered component is still preserved, this mixed form of feedback can be used with wavefront shaping to correct an aberrated focus [7,8,9].
Katz et al. [10] showed that even in the absence of ballistic light, light can be 'blindly' focused through a strongly scattering layer when feedback was generated by many fluorescent sources behind the opaque layer. Moreover, wavefront shaping with nonlinear feedback is not limited to continuous wave sources, but can also be employed to temporarily focus a scattered light pulse at one of the feedback sources [11,12]. Therefore, this blind focusing technique could potentially be a powerful technique for nonlinear deep-tissue microscopy. So far, a theoretical understanding of this iterative optimization process has been missing, and, therefore, it is not known under which conditions the blind focusing method will converge to a single diffraction-limited focus.
Here, we present and validate a model that accurately describes the statistics of blind focusing. Using this model, we show that focusing is only possible when certain minimal requirements are met. Additionally, we are able to predict the evolution of the optimized speckle pattern during the optimization process.
We will first derive the exact solution of the scattered field behind a strongly scattering layer for a single (known) realization of the scattering medium. Afterwards, we formulate a statistical model that describes the probability density function of the light intensity, averaged over the ensemble of possible samples. Our predictions are validated in a set of experiments with first, second and third order feedback.
2 Theory of blind focusing Figure 1 shows a simplified illustration of the blind focusing experiment. The goal of this experiment is to form a single diffraction limited focus through a scattering layer at the target plane using nonlinear feedback from M targets. A spatial light modulator (SLM) is used to control the field E a at the input plane. The feedback signals from the individual targets cannot be distinguished and only the total generated signal is recorded. We define this feedback signal as where n is the order of feedback, and E b represents the electric field at target b.
At the start of the experiment, an arbitrary starting wavefront is applied to the SLM. We then apply an optimization algorithm known as stepwise sequential wavefront shaping [1] to optimize the feedback signal. In contrast to the genetic algorithms used in other studies [10,13], our method is deterministic and can be analyzed analytically.
Each iteration of the algorithm results in an optimized incident wavefront E a , which is used as the starting wavefront of the next iteration of the algorithm. This process is repeated until the algorithm has reached a fixed point, where the optimized wavefront does not change anymore. In Appendix A, we prove that such a fixed point always corresponds to a (local) maximum of S.
In this section, we will derive how the target field E b changes as the nth-order feedback signal S is optimized by iteratively running the algorithm. Segment by segment, this algorithm optimizes S for a single incident field segments E a , while keeping the other field segments static. By introducing the transmission matrix T, we can describe the target field as a linear combination of input fields Here, t ba is an element of the transmission matrix, N is the number of independently controlled segments on the SLM and the superscript number k between brackets indicates the iteration number. We take the total incident intensity N a |E a | 2 = 1.
For the (k + 1)th wavefront shaping iteration, we use the optimized wavefront from the previous iteration as our new starting wavefront. Then we shift the phase of a single incident field segment a', such that E . We measure S as we vary φ from 0 to 2π in P steps. Using Eq. (1) and Eq. (2), the intermediary feedback signal can now be written as where * represents the complex conjugate. Since we only change a single segment of the incident wavefront, the effect on S is expected to be small. Therefore, we expandS (k+1) in terms of E with the nonlinearly weighted target field W (see appendix A for a rigorous mathematical treatment).
Following the processing steps as described in [14], we then find the optimized wavefront by isolating the contribution of φ using the following expression where c (k+1) is a normalization factor, which normalizes total incident intensity again to 1. The phase shift φ at segment a is set back to 0 and the process is then repeated for all other SLM segments. Finally, by inserting Eq. (5) into Eq. (2), we arrive at an expression for the resulting target field after the (k + 1)th wavefront shaping iteration We recognize that the optimized target field in Eq. (6) is the nonlinearly weighted sum over the phase conjugated fields propagated from the target locations. The nonlinear weighted field scales with the field strength to the power 2n − 1, and thus the brighter targets will have a stronger contribution to the optimized incident field. Performing multiple iterations of wavefront can therefore ultimately result in a focus being formed at one of the targets.
The expression in Eq. (6) can be connected to previous iterative wavefront shaping experiments. For linear feedback, W In this case, after performing the wavefront shaping algorithm multiple times, E (k) b will converge to the eigenvector of matrix T with the highest eigenvalue [15]. In other words, instead of forming a focus, the algorithm will instead optimize the total transmission through the scattering sample. Alternatively, for n > 1 in a weakly scattering medium, where T is close to unitary, the optimized target field becomes E b . Now the brightest targets are enhanced more than the dimmer targets until finally a focus is formed at the brightest target [16]. However, we will show that when T is not close to unitary, as is the case in strongly scattering media, maximizing S is not necessary equivalent to forming a focus.
3 Statistical model for blind focusing through strongly scattering media Next, we want calculate the optimized target intensity distribution through a strongly scattering medium with a nonunitary transmission matrix. Finding the exact value of the target field, using Eq. (6), requires full knowledge of the transmission matrix, which often cannot be obtained. Therefore, we will instead calculate the probability density function of E , averaged over the ensemble of possible samples. We assume that W is known and independent of T. The transmission matrix is assumed to be a random matrix with independent elements, such that the expectation value t ba t * b a = δ b b δ a a |t ba | 2 , where δ is the Kronecker delta. In Appendix B, we derive that for an imperfect wavefront shaping setup the probability density function of E (k+1) b is a complex normal distribution of the form: Here, the average optimized field µ (k+1) b and I 0 are given by where the quality of the wavefront modulation is described by the fidelity parameter |γ| 2 . This parameter's value ranges between 0 and 1, where a value of 1 corresponds with a perfect wavefront modulation. Unlike regular speckle fields, E will have a non-zero average value µ because of the performed wavefront shaping iteration.
follows a complex normal distribution with a non-zero mean, the corresponding optimized target intensity I | 2 follows a modified Rice distribution [17]. For µ (k) b = 0, this probability density function reduces to an exponential distribution as normally seen in speckle statistics. For now, we are mainly interested in the average optimized intensity and the corresponding standard deviation at target b, which are given by respectively. Note that I After the (k + 1)th wavefront shaping iteration, the average intensity at target b depends on the nonlinearly weighted fields at all targets. The targets which are generating a strong signal will contribute more to the feedback signal S, and will therefore, on average, be enhanced more than the weaker targets. We recognize that if nearly all of the feedback signal is coming from a single target, then the average intensity enhancement can be approximated by This expression is equivalent to the enhancement in a conventional single target wavefront shaping experiment [14]. However, when a large number of targets are contributing to the feedback signal, the expected enhancement is reduced to 1, and as a result, no focus is formed at all.
To summarize, in the blind focusing experiment, the statistics of the optimized intensities at the feedback targets are described by modified Rice distribution. Using Eq. (9), we can predict how much the light intensity will be optimized at the feedback targets after an iteration of wavefront shaping.

Experimental validation 4.1 Experimental setup
We will now experimentally validate the expected optimized intensity as predicted by Eq. (9). Our experimental setup is illustrated in Fig. 2. Light from a HeNe laser is expanded and modulated by a phase-only spatial light modulator (Hamamatsu X13138-07). With two lenses in a 4f-configuration, the SLM-modulated wavefront is conjugated to the back focal plane of a microscope objective (Zeiss A-Plan 100x/0.8), which focuses the light onto the surface of a scattering sample. As our source of feedback, we chose to use a CMOS camera (Basler acA640-750um), which measures the intensity distribution at the back surface of the scattering sample through an identical microscope objective. On the camera, we set groups of 3 × 3 pixels (corresponding to 0.16 × 0.16 µm 2 ) as independent targets for the wavefront shaping algorithm. Using this experimental setup, we can choose the order of feedback, and we can accurately set the number of targets contributing to the feedback signal. We would like to emphasize that in these experiments only the total n-th order intensities summed over all targets was used as the feedback signal, mimicking the nonlinear response of fluorescent markers in a multiphoton fluorescence excitation microscope.
As a scattering sample, we used a glass substrate dip-coated in a suspension of 5% zinc-oxide (Sigma Aldrich, average grain size 200 nm) and demineralized water. The thickness of the layer were measured to be 37.6 ± 9.8 µm. Based on previous work [18], the transport mean free path of the zinc-oxide sample is expected to be around 0.6 µm, ensuring that all light passing through the sample is multiple scattered. The sample is mounted on a translation stage (Zaber T-LSM050A), to allow the interrogation of different sample locations for statistical averaging. The average wavefront shaping fidelity of our setup was measured to be |γ| 2 = 0.39 with N = 208.

Blind focusing experiment with two feedback targets
We performed a blind focusing experiment, where two targets were simultaneously optimized using the wavefront shaping algorithm as described before. The targets are separated by a distance 4.4 µm. On the SLM, a random pattern of N = 208 square segments was displayed, matching the size of the light beam on the SLM. The experiment was performed 100 times, changing the initial SLM pattern and sample lateral position in between every experiment. In Fig. 3(a) and (b), examples of the intensity distributions measured on the back side of the sample, before (a) and after (b) the optimization, are shown. Here, targets 1 and 2 are indicated by the red (right) and blue (left) circle, respectively. In these figures, the starting intensity of target 2 is much higher than the intensity target 1. As a result, the contribution of target 2 to the feedback signal is much larger than target 1, and therefore, only the intensity of target 2 is enhanced. Fig. 3(c)-(e) show the optimized intensities I (1) of target 1 (red circles) and target 2 (blue squares) as function of the ratio of the initial intensities I (0) of the 2 targets, for first-order (c), second-order (d), and third-order (e) feedback. The optimized intensity is normalized by I max = N |γ| 2 I 0 , which is the average intensity obtained by performing a single-target wavefront shaping experiment. The solid lines and the shaded areas indicate the average optimized intensities and the standard deviation as predicted by Eq. (9) and Eq. (10).
In Fig. 3(c)-(e), we see that even when both targets contribute to the feedback signal, the target intensities are not necessarily equally enhanced. Neither are the optimized intensities I (1) randomly distributed, but rather the optimized intensities are directly related to the starting intensities of target 1 and 2, as expected from Eq. (9). In the case of n = 1, the optimized intensities at the targets are linearly proportional to the initial intensities, since in that case |W For higher-order feedback measurements, the relation between the initial intensity ratio and the optimized intensity becomes nonlinear. As a result, a small difference in I (0) between the two targets can result in a big difference between optimized intensities. For instance for n = 3, only approximately 65% of the total initial intensity is required in one target to ensure, in almost all cases, that only the intensity of that target will be optimized. This thresholding effect for n > 1, observed in Fig. 3(d-e), is the mechanism that allows the blind focusing method to form a single diffraction-limited focus even when multiple targets contribute to the feedback signal. The data is in good agreement with our theoretical predictions and most of the data points fall within the statistical variation as described in Eq. (10). Other deviations might be explained by fluctuations in the fidelity during the experiment.

Theory
Now, we want to use our statistical model to find out under which circumstances the blind focusing method is able to form a single focus, when M targets are contributing to the feedback signal. Instead of analyzing convergence behaviour starting from a random speckle pattern, we analyze the case where light is already focused to one target before 3 µm 3 µm Figure 3: Results of the blind focusing experiment optimizing for two targets simultaneously. (a) Example image of the initial intensity distribution, and (b) intensity distribution after the first wavefront shaping iteration. The red (right) and blue (left) circles correspond to targets 1 and 2, respectively. The optimized intensities of target 1 (red circles) and target 2 (blue squares) are plotted as function of the initial ratio of the intensities in the two targets for (c) first-order, (d) second-order, and (e) third-order feedback. The predicted mean value for the optimized intensity and the corresponding standard deviation are represented by the colored solid lines and the shaded areas, respectively.
the optimization. The focus should be preserved after the wavefront shaping iteration when this focus corresponds to a fixed point of the algorithm. To quantify the intensity in the focus, we use the enhancement, which is defined as η ≡ I/I 0 . The starting enhancement in the focus is given by η (0) . Furthermore, we assume that all other (M − 1) targets are exponentially distributed with an average starting enhancement of 1. Inserting these parameters into Eq. (9) produces an expression for the expected optimized focus enhancement after a single optimization iteration We recognize that η (1) can be smaller than η (0) for a large M . In order for the blind focusing method to be able to form a focus, the average focus enhancement after the optimization should be equal to the enhancement before the optimization for some η (0) > 1. In Appendix C, we show that (given N , |γ| 2 and n) an upper limit for the number of feedback targets can derived, which is given by Whenever M > M max , a focus can on average not be formed. Rather, the wavefront shaping algorithm optimizes the intensity and contrast of the full speckle pattern to maximize the feedback signal. We see that the order of feedback n can be increased to guarantee blind focusing convergence for a larger number of contributing targets. Moreover, by increasing N by a factor of α, M max increases by a factor of α 2n−1 .

Blind focusing experiment with M targets
To verify the prediction in Eq. (11), we performed a second experiment using a large feedback area, containing M feedback targets. For this experiment, the same experimental setup was used as in the first experiment. We start our experiment with a pre-optimized focus at a single target, whereas the remaining targets in the feedback area are illuminated by a random speckle pattern. The pre-optimized focus is obtained by first optimizing for a single target in the center of the camera frame. Afterwards, a wavefront shaping iteration is performed using the total second-order feedback signal from all targets within a small and a large region of interest (ROI), which have a radius of 2.1 µm and 10.5 µm. Based on the average speckle size, the number of targets in the ROIs are estimated to be M = 96 and M = 2400 for the small and the large ROIs, respectively. The blind focusing experiments were performed for N = 80 and N = 208, and were performed 100 times for both ROIs. In between each experiment, the lateral sample position was changed and the intensity of pre-optimized focus was varied by adding a controlled amount of uniformly-distributed noise to the starting wavefront.
In Fig. 4(a), an example image of a starting speckle pattern with a pre-optimized focus is shown. The red and purple circles in the figure indicate the size of the small and large ROIs. In Fig. 4(b) and (c), the enhancement of the preoptimized focus intensity after the blind focusing experiment η (1) is plotted as function of the starting enhancement of the pre-optimized focus, for N = 80 (Fig. 4b) and N = 208 (Fig. 4c). The experiments with small and the large ROI are represented by the red circles and the purple squares, respectively. The expected value for η (1) , as predicted by Eq. (11), and the corresponding standard deviation are represented in Fig. 4(b) and (c) by the colored solid lines and the shaded areas, respectively. The black solid lines indicates the identity lines, where η (1) = η (0) . In Fig. 4(b) and (c), two general regions can be recognized. The data points above the identity line represent experiments where the focus enhancement increased, whereas the data points under the identity line represent a decrease in the focus enhancement. In Fig. 4(b), in the experiments with M = 96 and N = 80, nearly all measured η (1) values lie above the identity line, meaning that the focus enhancement increases after the iteration of wavefront shaping. However, in the experiments with the large ROI, η (1) is almost always lower than η (0) . These results suggest that the large ROI contains too many competing targets for the pre-optimized focus to be preserved. In other words, the wavefront shaping algorithm favors optimizing the intensity of multiple targets within the ROI over preserving the intensity in the pre-optimized focus.
In Fig. 4(c), the number of controlled SLM segments was increased to N = 208. In this figure, the η (1) curves cross the identity curves for both the small and large ROIs, indicating that the intensity in the pre-optimized focus was further enhanced. Increasing N thus allows blind focusing to preserve a focus for a higher number of targets that are contributing to the feedback signal. The experimental data is in excellent agreement with the prediction in Eq. (11).
These results demonstrate that for second-order feedback, the blind focusing method will not form a focus using only 80 SLM segments when an area with 2400 targets is illuminated. This is in accordance with the prediction in Eq. (12), since the number of feedback targets in the large ROI (M = 2400) exceeds the upper limit, M max = 751, whereas the experiment with N = 208 (in Fig. 4(c)) has a higher M max of 13181 targets.

Discussion
We studied wavefront shaping with nonlinear feedback through a strongly scattering sample using the stepwise sequential algorithm [14]. This deterministic algorithm allowed us to derive the exact solution for the optimized field at the feedback targets, which is given in Eq. (6). This expression shows that in weakly scattering samples, the optimization algorithm takes higher orders of the initial target field with every iteration. As a result, only the brightest speckles at the target plane will be enhanced [16]. Therefore, wavefront shaping and adaptive optics techniques can be used to improve the intensity of the focus in multiphoton fluorescence excitation microscopy without the need for a guide star [7,8].
Katz et al. [10] have shown that a focus can be formed through a strongly scattering layer without the need of a localized reporter. We showed that in strongly scattering samples, the statistics of the optimized target intensities are accurately described by the Rice distribution. However, due to the large standard deviation in the distribution, it remains hard to predict at which target the focus will be formed when the targets are illuminated with a random speckle pattern.
In Appendix A, we showed that a fixed point of our optimization algorithm always corresponds to a maximum of the feedback signal. Moreover, we experimentally demonstrated that when the number of targets exceeds M max from Eq. (12), our algorithm is, on average, unable to form a focus. As a result, maximizing S does not always guarantee the formation of a focus, but will instead optimize the intensity and contrast of the speckle pattern. This limitation applies to all optimization algorithms using the total nonlinear signal as feedback, including the commonly used genetic wavefront shaping algorithm [13].
In our experiments, we only considered targets distributed across a flat two-dimensional image plane. However, our model can also be employed to point targets on an arbitrarily shaped plane. Therefore, our model can easily be generalized for structure with targets distributed in three spatial dimensions. In such samples, the light intensity will spread out more the deeper the light propagates into the sample. As a result, the number of illuminated targets will also rapidly increase with imaging depth.
In our statistical model, we assumed that the nonlinearly weighted field W is independent of the transmission matrix. As a result, we ignore certain properties of the matrix T known from random matrix theory, such as the distribution of transmission eigenvalues [19]. Nevertheless, our statistical model is in a surprisingly good agreement with the experimental data. At the moment, our model is unable to predict the rate of convergence. Therefore, we believe that studying these properties of random matrices will be interesting to get a full understanding of this nonlinear optimization process.

Appendix A: Convergence to local maximum
We used the stepwise sequential optimization algorithm [1] in an attempt to maximize the total feedback signal S. For linear feedback from a single target, it is well known that this algorithm finds a global maximum for the intensity at that target in a single iteration [14]. In the case of nonlinear feedback originating from multiple targets, however, the optimization problem has multiple local maxima. In this situation, it is not directly trivial that the algorithm finds a local maximum of S; it is not even directly clear that the stepwise sequential algorithm increases S at all.
Here, we will prove that all local maxima of S correspond to attractive fixed points of the algorithm (Eq. (6)) and vice versa. Consequentially, when the algorithm converges, it converges to a local maximum of S as well. Furthermore, this proof implies that each local maximum of S has a finite region of attraction for which the algorithm will converge to that local maximum. Note that we cannot exclude the existence of initial conditions for which the algorithm does not converge. However, we have not observed these cases in our experiments. In order to simplify the derivations, we introduce a compact vector notation, replacing E a by a, E b by b etc. For example, Eq. (2) can now be written with T the transmission matrix. The incident field is normalized so that a = 1.

Wirtinger calculus
To analyze convergence to a local maximum, we will expand the system around a fixed point a * . Here a technical complication arises: due to the complex conjugate in S, it is not a holomorphic function of a, so the complex derivative ∂S/∂a does not exist.
In order to avoid this complication, we use Wirtinger calculus to calculate the derivatives [20]. In Wirtinger calculus, S(a) is replaced by a function S(a, a), where a and a are independent variables: technically, a is not equal to the complex conjugate of a. However, as long as we restrict ourselves to only evaluate S(a, a ) in the subspace a = a * , we can perform all derivatives in the usual manner and obtain correct results. Using the compact notation and Wirtinger calculus, we can now write Eq. (1) as where () b denotes the b-th element of the vector in the parentheses. As an example, we use Wirtinger calculus to calculate the first order Taylor expansion of S for small perturbations ∆ around a S(a + ∆, a + ∆ ) = S(a, a) We can now evaluate the derivative towards component a i of the incident field a where the elements of w and w are given by By evaluating the derivative for all elements a i , we get the vector derivatives We can use the Taylor expansion Eq. (A3) to calculate what happens during the optimization process, when the phase of a single element of the incident field is changed. The perturbation corresponding to changing the phase of segment a is given by Inserting ∆ into Eq. (A3) gives Eq. (4) from the main text.

Attractivity of fixed point
The algorithm in Eq. (5) can be thought of as a cyclical series of mappings a (k) → b (k) → w (k) → T † w (k) → a (k+1) , etc., where T † is the conjugate transpose of T. We rewrite Eq. (5) as f(a, a) = T † w(a, a ) where f now defines the mapping from a (k) to a (k+1) . Using Wirtinger calculus, we can make a first order approximation of the mapping for a small perturbation ∆ around a fixed point of the mapping, a * = f(a * ).
where the matrix is J f (a * ): the Jacobian of f, evaluated at the fixed point a * . When the fixed point is attractive in some finite region, each iteration brings a closer to a * , so we must have a (k+1) − a * ≤ q a (k) − a * , with 0 ≤ q < 1. From Eq. (A9), we see that this condition is equivalent to saying that the spectral radius of the Jacobian ρ(J f ) < 1. We will show below that this condition is always met at local maxima of S.
From the definition of f it is clear that any perturbation in the direction of a * will have no effect at all. Therefore, we restrict ourselves to perturbations perpendicular to a * , i.e. ∆ T a * = 0, hence ∆ T T † w * = 0. Under this condition, we can find a simple expression for terms of the form ∆ T ∂f/∂a which we will use in the next section. Here the product rule was used in the first step, and orthogonality of ∆ and a * was used in the second step.

Local maximum
To find local maxima of S under the constraint that a = 1, we apply the method of Lagrange multipliers and minimize the Lagrange function S L (a, a) = S(a, a) − λ(a T a − 1), (A11) where λ is a Lagrange multiplier and a T a − 1 represents the constraint that a = 1.
The first order conditions for a local maximum follow by equating the first derivatives of S L to zero, giving with the solution λ = T † w * and a * = T † w * /λ, proving that a fixed point of f is also a critical point of S L . In order to prove that this stationary point is a local maximum (and not a local minimum or a saddle point), we need show that any small perturbation that maintains the first order conditions decreases the value of S L . In order to do so, we evaluate the second derivative (the Hessian matrix H) of S L considering only perturbations perpendicular to a * , i.e. perturbations that maintain the constraint a = 1 to the first order. Using Eq. (A6), we arrive at Using Eq. (A10), we find with I the identity matrix. Since ρ(J f ) < 1, from Eq. (A16), J f − I is negative definite, proving the last step: every perturbation ∆ decreases the signal. Hence, every fixed point a * corresponds to a local maximum of the constrained optimization problem.

Conclusion
This final result in Eq. (A16) shows that any small perturbation around the fixed point decreases the signal. In conclusion, we demonstrated that, when a * is an attractive fixed point of mapping f, the resulting signal S(a * , a * ) must be a local maximum. The converse is also true, since both statements are equivalent to the condition that ρ(J f ) < 1.
Finally, as detailed in the main text, even though the algorithm maximizes S, such a maximum does not necessarily correspond to a focus.

Appendix B: Derivation of the blind focusing statistical model
In this appendix, we derive Eq. (7), the complex normal distribution of the optimized target field when blind focusing through a strongly scattering sample. We start by constructing the optimized incident field, which is given in Eq. (5).
In an imperfect wavefront shaping setup, the quality of the wavefront modulation is described by the fidelity, |γ| 2 . The constructed input field is given byÊ is the desired input field and ∆E a is a normalized field, which is by definition orthogonal to E (k+1) a [21]. When this imperfect input field is inserted into Eq. (2), we find that the constructed optimized target field becomeŝ Here, ζ b ≡ N a t ba ∆E a which is assumed to be an uncorrelated scattered field with ζ b = 0 and var(ζ b ) = |t ba | 2 . Note that for a perfect setup with γ = 1,Ê (6)). We recognize that Eq. (B1) can be written as a sum over N independent random variables χ â When N is large, by the central limit theoremÊ has a complex normal distribution.
To find the average and variance ofÊ , we start by calculating the first and second raw moments of the terms χ a . As stated in the main text, we assume that W is known and independent of T, and that t ba t * b a = δ b b δ a a |t ba | 2 . Under these assumptions, the first moment is given by and the second moment Realizing that for a Gaussian distribution, |t ba | 4 = 2 |t ba | 2 2 , we find Next, we calculate the value of c (k+1) , which was introduced in Eq. (5) to normalize the total incident intensity after the optimization. We assume that, for large N , the normalization factor is self-averaging, such that To obtain the mean ofÊ , we can simply add the means of χ a and ζ b , since the two variables are uncorrelated. We calculate the optimized field average by inserting Eq. (B3) and ζ Similarly, we can calculate the variance ofÊ When we substitute µ Appendix C: Derivation of the blind focusing requirements In this section, we calculate the minimal requirements for the formation of a focus using the blind focus method. We assume that, when a focus is formed, only the intensity at the focus location is enhanced and that all other targets have an exponentially distributed enhancement with an average of 1. A focus can be formed when the expected focus enhancement of the next iteration η (1) (as described by Eq. (11)) is larger than or equal to the current focus enhancement, for some η (0) > 1. For simplification, we instead consider η (1) ≥ η (0) + 1, which is a slightly more restrictive condition. We can write this condition as N |γ| 2 (η (0) ) 2n−1 (η (0) ) 2n−1 + (2n − 1)!(M − 1) Next, we derive the minimum requirements for this condition to be satisfied. Therefore, we proceed by finding to find out in which cases this maximum value satisfies the condition in Eq. (C2). We maximize this function by equating the derivative of this function towards η (0) to zero, giving (2n − 1)(η (0) ) 2n−2 = (2n − 2)N |γ| 2 (η (0) ) 2n−3 (C4) We insert this maximized value for η (0) into Eq. (C2) − 2n − 2 2n − 1 N |γ| 2 2n−1 + N |γ| 2 2n − 2 2n − 1 N |γ| 2 2n−2 ≥ (2n − 1)!(M − 1).

Funding
This work was funded by the European Research Council under the European Union's Horizon 2020 (ERC-2016-StG-678919).