Light focusing is at the basis of a remarkable number of groundbreaking advances in numerous fields ranging from imaging [1], characterization [2] and fabrication [3]. To achieve the maximal localization, a coherent source is typically coupled to an objective lens: the well-defined propagation direction and shape of a laser beam is crucial to control light while any photon deviating from the ballistic path reduces the confinement efficiency. Indeed, experiments in presence of turbidity typically cannot rely on confined light structures. Starting from 2007, a technique known as wavefront shaping [4] (WS) has been borrowed from astrophysics to bench-top optics [5] promising strong focusing also for light experiencing strong scattering [6,7].

Indeed, each scattered light path $n$ is affected by a phase delay $\Delta \phi _n$ with respect to its non-scattered counterpart. The fundamental idea of WS is to introduce an adaptive optics device (such as a spatial light modulator, SLM), capable to tune phase independently for each mode. The ballistic propagation can be restored just by introducing a proper phase delay to cancel the scattering-induced delay [8].

Finding the phase delays $\Delta \phi _n$ may be a challenging and time-consuming task. Two main alternative approach have been proposed so far: the so called Sequential Optimization (SO) or the Transmission Matrix Measurement (TMM). SO is a “brute force” approach in which every segment of the SLM is randomly varied and the change is accepted if intensity at the focusing point increases. The main idea of SO is thus to explore a huge configurations space in a process which becomes extremely time consuming process when the number of controlled modes is large.

The alternative approach, TMM, is based on the Transmision Matrix ($\mathbb {T}$) representation of the scattering process [5,9–24]. $\mathbb {T}$ is composed by a set of complex parameters describing attenuation and dephasing of each output with respect to each input. The TM theory predicts a special role for eigenvectors of the disordered, complex $\mathbb{T}\mathbb{T}^\dagger$ matrix: for example if all the input degrees of freedom are controlled the high transmitting “open channels” [14] may be activated. In TMM, the elements of $\mathbb {T}$ are individually measured [12,25] to find eigenvectors of the $\mathbb{T}\mathbb{T}^\dagger$ matrix which enable to the maximize the intensity. TMM efficiency is however limited by the requirement of a reference beam, which reduces the focusing performance. To address this point, several alternative to TMM have been proposed: the use of sophisticated phase retrieval algorithms [17,26], self reference method [27], dual reference with Hamadard basis segments [28], pseudo-linearity approximation [20], and intensity transmission matrix [29,30] .

Here we realize experimentally a new reference-less transmission matrix measurement technique, the Complete Couplings Mapping (CCM) and demonstrate its advantages. Moreover we compare the CCM method with the SO, demonstrating the effects of roughness in the potential energy landscape

The idea stems from an a parallelism between the $\mathbb{T}\mathbb{T}^\dagger$ matrix and the Hopfield model coupling matrix (CM, [31] matrix)

In the Hopfield model, the coupling matrix(CM) represents the outer product of memory patterns stored into a neural network [32,33]. Given a stimulation pattern of neuronal activities the CM dives the subsequent dynamics of the firing activity. A memory pattern is retrieved when the neuronal configuration is aligned with that pattern, i. e., when the dynamics reaches a steady state, also termed an “attractor” in the space of neuronal configurations. This steady state is an eigenvector of the CM. In this sense memory retrieval process and light focusing are similar as both are the outcome of an eigenvector retrieval problems.

In this paper, by applying the CCM, we will measure the CM without the help of phase retrieval algorithms and without the requirements of a reference beam. Then we will demonstrate that our approach enables to reach the theoretically predicted limits in efficiency.

1. Theoretical background

Indeed, a critical feature for WS is its efficiency. If we term $I_{\rm focus}$ the intensity at a given location at the end of the shaping process and $I_{\rm speckle}$ the intensity distributed in a naturally speckle pattern generated with a randomly input pattern, the parameter defining the efficiency of the wavefront shaping technique is [34]

This parameter, $\eta$, measures the effectiveness of a modulation technique to deliver light at a desired location. Regardless of the specific WS method employed, the maximum achievable efficiency $\eta$ is connected to the number of total controlled segments $N$ in the input as [34]

We start defining the problem. The field $E^{(\nu )}$ at a given target location $\nu$ on the output screen is determined by the input fields $E_n$, $n=1,\ldots, N$ and by the structure of the transmission matrix $\mathbb {T}$ whose elements $\mathbb {T}_{n\nu }={t}^{(\nu )}_n$ connect input and output fields. The intensity transferred from all inputs to the $\nu$ output target is written as

In typical wavefront shaping experiments, the light field ${E}$ before the modulator is steady and, thus, it is convenient to define a vector $\mathbf {\Theta }$ merging field and transmission contributions

By employing $N$ segments, the intensity on the output location at the far side of the disordered material is fully defined by the coupling matrix $\mathfrak {J}=\mathbf {\Theta } \otimes \mathbf {\Theta }^{\dagger}$ whose elements $\mathfrak {J_{nm}}=\Theta _n \bar \Theta _m$ contain all the interference pairs.

It is well known (see S.I. and Ref. [12]) that the eigenvector of $\mathfrak {J}$ corresponds to the input configurations providing the highest amount of intensity at the target. As often the ${\mathbf {\Theta }}$ vector, and consequently the coupling matrix $\mathfrak {J}$, is inaccessible or extremely difficult to measure, light focusing behind disordered materials is usually achieved by a sequential optimization (SO) protocol in which many random iterations (amplitude or phase changes on the input modes) are performed in order to maximize the intensity at a target.

In a recent paper by the authors [31], a connection was established between the optimization of light intensity through disordered materials and the physics of disordered systems. Specifically, it was found that the mathematical equations governing light diffusion can be linked to the physics of spin glass and Hopfield networks.

Hopfield networks, introduced in [38], provide a mathematical framework for understanding biological memory.

In this model memories are represented by patterns of the neuronal activity. A stored memory is retrieved when the configuration of a set of neurons is very similar to its pattern. The set of memories stored in a network are defined by the strength of the interconnections between neurons, represented by a connectivity matrix (also known as the synaptic matrix).

The network dynamics evolves by iteratively multiplying the neural state vector (containing ones for active neurons and zeros otherwise) with the connectivity matrix and applying a thresholding operation. Memory elements are recognized as "fixed points" (stationary states) of the network dynamics.

The connection between optics and Hopfield networks arises from a formal similarity: the equation describing the energy of a Hopfield network is identical to the equation used to estimate optical intensity (see (9) and [39]).

While Hopfield networks can store multiple memories using the Hebb’s rule, a limitation on memory storage was found by Amit, Gutfreund, and Sompolinsky [40,41]. This limitation arises from the emergence of uncontrolled spurious states, leading the system into a complex regime.

Interestingly, when the Sequential Optimization (SO) algorithm is employed to maximize the intensity at multiple output targets simultaneously, its driving equation is identical to the energy expression of the Hopfield model for memory storage. This connection highlights the correspondence between the optimization processes in optics and the memory in neural networks.

Overall, the relationship between optics and Hopfield networks provides valuable insights into memory storage mechanisms and optimization processes, opening up new avenues for research in both fields.

When few targets are involved, in the statistical mechanical description, the transmission through the random optical system corresponds to the regime of memory recovery in the Hopfield neural network. We will, therefore, focus on this regime. Indeed, in the analogy, the transmission matrix elements linking each input to a single output correspond to a memory pattern stored in a neural network. And the input light mode fields $E_n$ play the role of a neuron activity.

The intensity optimization corresponds, from a statistical mechanics point of view, to a zero temperature dynamics in a corrugated random landscape of mode configurations whose lowest minima are the optimized solutions. Each random transmission matrix yields a different corrugation, i.e., a different energy landscape. This implies that any small energy/intensity barrier around a non-optimal configuration incidentally encountered in the optimization procedure would not allow to optimize any further. For instance, if the input fields are combined to provide $N$ variables only taking $\pm 1$ values, the system relaxing towards equilibrium can be stuck in a non-optimal solution because of a barrier that could be overcome just by a single sign switch of the field. We will refer to this sign switch as “spin-flip”, since the $\pm 1$ variable is termed Ising spin in the statistical mechanical analogue. The system is, thus, stuck in a “one-spin-flip-stable” minimum corresponding to a non-optimal intensity.

In the following we will study how the presence of such metastable (shallow, one-spin-flip-stable) minima affect the SO in the $[0$ $\pi ]$ phase modulation. We will demonstrate that the SO is not always capable to reach the maximal efficiency for light focusing and why this comes about. Then we will provide a new method (the CCM) to retrieve light behind a disordered medium without the use of a reference. We will demonstrate that CCM is capable to return an intensity enhancement very close to the theoretical predictions.

2. Experimental results

To investigate the effect of dynamic arrest in the SO, we perform repeated sequential optimizations. Starting from a random configuration, the SO algorithm flips the DMD segments (Spins) many times (a pseudocode for the SO is reported in Supplement 1). At each segment/spin flip, the intensity (energy) is measured and the change is accepted upon an intensity increase (energy decrease). The process is iterated up to when intensity saturates.

In Fig. 1(b)-e we report the (sorted) results for two illustrative target locations ($\nu _1$ and $\nu _2$ locations also reported in Fig. 3(d))).

 

Fig. 1. Experimental setup and experimentally retrieved Focusing States. Panels a) shows a picture of the experimental setup: LS Laser source; M Mirror; P Power control; H heaters; IA input angle; C camera; D diffusive sample. Lens Focals: L$_1$ 200 mm; L$_2$ 60 mm L$_3$ 25 mm ; L$_4$ 100 mm.Panels b) and d) report the configurations $\sigma _r$ obtained for different replicas of the experiments respectively in output targets $\nu _1$ and $\nu _2$. Note that the $r$ index has been sorted in order to have $\sigma _r$ pertaining to the same state clustered. b) shows two different clusters (the first from row 1 to row 28) the second 7 clusters (the first from row 1 to row 15). Panels c) and e) (respectively corresponding to $\nu _1$ and $\nu _2$) show the degree of similarity $Q$ calculated for all the replicas to all the replicas (with self-overlap reference $Q(\sigma _{r},\sigma _{r})=1$). The existence of multiple clusters is confirmed by this representation of $Q$.

Download Full Size | PDF

In particular, panels 1(b)) and 1(c)), refer to a location ($\nu _1$) displaying only two minima and 1(d)) and 1(e)), refer to a location ($\nu _2$) with more (meta)stable configurations. Here the final states $\sigma _r$ are represented with black and white pixels corresponding to DMD segments placed in the −1 or +1 segment configuration, respectively. In panel 1(b)) the first $28$ $\sigma _r$ rows, (indicated by a red brace on the left) resulting from as many sequential optimization processes, are very similar, pertaining to a first cluster of measurements. On the other hand, rows $29$-$60$ (indicated by a red brace on the right) are, instead, pertaining to a second cluster of $\sigma _r$. Panel 1(c)) reports the degree of similarity $Q$ between all the measurements acquired at target $\nu _1$. The degree of similarity is defined as the scalar product between configurations:

A value of $Q$ close to 1 indicates that $\sigma _{r1}$ and $\sigma _{r2}$ are very similar, while $Q = 0$ indicates orthogonal input vectors. $Q= -1$ indicates spin reversed vectors.

At target instance $\nu _1$ two different states can be clearly identified ($\sigma ^\uparrow$ and $\sigma ^\downarrow$) and they correspond to two “spin reversed” solutions ($\{\sigma _n\}\rightarrow \{-\sigma _n\}, ~\forall ~n$). In other words, at location $\nu _1$ a single optimal configuration is found, together with its additive opposite. The intensity landscape, akin to the potential energy landscape in statistical mechanics, only displays a couple of global minima (degenerate and symmetric under spin inversion) making the optimum solution easy to recover. Indeed location $\nu _1$ has a number of clusters $N_C$ equal to $2$.

The situation is quite different at location $\nu _2$ where multiple clusters are identified. The presence of many clusters is connected to the roughness of the potential energy landscape of the system and the presence of shallow metastable minima and more than a single couple of stable minima. The instance at location $\nu _2$ turns out to display $N_C=7$.

Notice that, even if $\nu _1$ and $\nu _2$ are distant in space and, thus, are identified by two completely different coupling matrices ${\mathfrak {T}}$, they have been selected for a fair comparison as they display very similar pre-optimization average intensity.

To understand the origin of such unexpected variability in the $N_C$, we performed the same experiment on many different output locations $\nu$, namely $174$ targets, computing for each location the number of clusters $N_C$, that is a proxy for the number of arrested, one-spin-flip-stable, sub-optimal configurations in the system. In all cases where $N_C>2$, indeed, taken away the possible optimal state (and its spin-reversed one), there will be configurations of sub-optimal intensity.

We also compute the focusing efficiency $\alpha$ defined in Eq. (2). This is shown in Fig. 2) (gray squares) plotted versus the total number of states $N_C$ where the error bars indicate the standard error of the $\alpha$ distribution. It is possible to note how the optimization efficiency $\alpha$ decreases when $N_C$ increases and how $\alpha$ is always smaller than the theoretically predicted value $1/\pi$ (indicated by a dashed line in the graph, see also Supplement 1) of the maximal focusing efficiency.

 

Fig. 2. The Experimentally retrieved efficiency. Efficiency parameter $\alpha$ versus the number of clusters $NC$. Gray squares indicate results obtained with the SO approach while orange disks are obtained with the CCM. Error Bars represent standard error. The dotted line represent the theoretical value for the averaged-over-disorder CCM approach. See Supplement 1

Download Full Size | PDF

2.1 Complete coupling mapping method

To further investigate the role of the number of clusters, we resorted to a new approach to measure the coupling matrix Eq. (4), the Complete Couplings Mapping (CCM), which is depicted in Fig. 3. CCM is particularly fitted to exploit the DMD characteristics: the very high speed and the capability to “turn off” segments.

 

Fig. 3. Full Experimental CM measurements by CCM: panels a) and b) report the light intensity captured on the camera when only the mode n is $+1$ (a) or when only the segments n and m simultaneously are in the $+1$ state (b), enabling to retrieve both $I_{n,n}$ and $I_{n,m}$. The sketch c) represent the configuration of the DMD corresponding to b). The small squares are the DMD micromirrors which are organized in 2x2 segments. Two of them ($n$ and $m$) are “activated” (each represented as two larger azure squares arranged in diagonal fashion). In panel d) we report a typical light focus obtained from Eigenvectors obtained by the CCM method. Note the locations employed in 1 are indicated as $\nu _1$ and $\nu _2$. Panel e) is the same as d) with higher magnification: scale bars are 500 $\mu m$ (in a), b) and d)) and 150 in e) $\mu m$ respectively Inset of panel e) reports the intensity profiles obtained at different times during the fast scanning through the disordered medium. The color indicates the scanning row as indicated in the panel e) Note that the time axis in panel f has been constructed disregarding sensor related delays (data transfer and camera readout).

Download Full Size | PDF

Thus, by turning off all the segments except a pair (indices $n, m$, element is set in the $+1$ configuration), and experimentally measuring intensity at the target location $I_{nm}^{(\nu )}$, we are able to extract the individual elements of $\mathfrak {J}$ without the need of any reference beam. Indeed, if we express the elements of the transmission matrix eigenvector (5) in their real and imaginary parts, $\Theta _n=\xi _n +\imath \chi _n$, we obtain (see Supplement 1, below we neglected the positing index ($\nu$) to simplify notation):

To measure the full coupling matrix, we need to individually measure all the matrix elements, turning on sequentially all the $N$ segments and successively all the $N(N-1)/2$ pairs of segments. With $N=256$, e.g., we need to perform $32896$ measurements: a prohibitive task with a (slow) liquid crystal based spatial light modulator, but perfectly feasible experimentally in few minutes with a fast DMD.

We report the measured coupling matrix in the inset of Fig. 4(a). A couple of interference patterns and a sketch of a typical DMD configuration are reported in Fig. 3 In our experimental setup, intensity is recorded by a camera with a Region Of Interest (ROI) of $128 \times 128$ pixels. Thus $P=16384$ different locations $\nu$ can be recorded simultaneously in a single measurement run. Then, for each location $\nu$, the elements of the real part of the coupling matrix $\mathfrak {J}$ are extracted from intensities, employing Eqs. (6)–(7) and the eigenvalues $\lambda ^{(\nu )}$ and the eigenvectors $\mathbf {e}^{(\nu )}$ are retrieved. In this setting the measurement procedure is completed in about $5$ minutes. We realized a custom intensity analysis software which handles this extended data-set, retrieving all the elements of the coupling matrix. With the proposed formalism the focusing configuration can be excited presenting as an input pattern for $J$ the eigenvector $\mathbf {e}$, see Supplement 1. The CCM focus obtained experimentally, are reported in Fig. 3(d)-e. The focus size (15 $\mu$m) is defined by the numerical aperture of the last lens (100 mm focal Thorlabs LB1676-A-ML lens). The inset 3d plot reports about the systems speed when employed as a laser scan. Different plots report the experimentally retrieved intensity profile at the foci during the scanning. Note that the focus is displaced in 100 $\mu$s, thus CCM with DMD can reach laser scanning speed with a 10kHz rate. This result demonstrates CCM working at multiple location with a single measurement run. Light focusing at a diffraction limited spot, is fundamentally different from the optimization of the total (wide field) transmission though disordered media [14]. Control of total transmission is difficult in the visible, [42], but may be achieved in the microwave regime due to complete control of the channels with a modulator or by engineering a complementary scattering medium [43]. In both cases a full measure of the $\mathbb{T}\mathbb{T}^\dagger$ matrix (as in CCM) is needed.

 

Fig. 4. Eigenvalues retrieved from experimentally measured optical CM. a) Sorted eigenvalues of the bi-dyadic single target coupling matrix associated to the TM. The inset shows the real coupling matrix $J$ matrix: the values of $J_{nm}$ are reported in image format versus the $n$ and $m$ modulator segment index. The presence of highly valued rows or columns is resulting from the highly correlated nature of $J$ as it results from the external product $\mathbf {t} \mathbf {t}^{\dagger}$. In panel b) we report the $I_{\rm focus}$ obtained with SO. Panel c) reports the $I_{\rm focus}$ estimated employing Eq. (9) while panel d) provides the $I_{\rm focus}$ obtained experimentally. All b), c) and d) matrices are retrieved in the same ROI and constructed so that each pixel is reporting the maximal intensity retrieved at the end of the optimization at that location.

Download Full Size | PDF

2.2 Optical coupling matrix

Now we proceed to analyze the experimentally measured coupling matrix $J$. At difference with single memory Hebbian matrices (which are dyadic), CM is a bi-dyadic real-valued matrix constructed with the Hebb’s rule [44] starting from the real and imaginary coefficients of the transmission matrix $\boldsymbol {\Theta }=\boldsymbol {\xi }+i\boldsymbol {\chi }$. In particular,

The average structure of the eigenvalues found from experimental data, is reported in Fig. 4(a)). We measure two positive eigenvalues (the 224 $th$ and the 225 $th$) and a set of $N-2$ eigenvalues very close to 0, as expected for a bi-dyadic matrix. A small deviation from the ideal picture is due to noise in the interference measurement.

We can employ the eigenvector $\mathbf {e}^{\rm max}$, correspondent to the largest eigenvalue, to transmit with the maximal focusig efficiency. Eigenvectors are composed by $N$ real numbers, while the DMD is only capable to deliver $N$ $\pm 1$ valued inputs. Thus, to get the DMD input array which is best approximating the eigenvector $\mathbf {e}^{\rm max}$ elements, we define the vector of signs $\boldsymbol {\sigma ^{\rm max}}=\text {sign}(\mathbf {e}^{\rm max})$. The resulting modulation factor $\alpha$ obtained within this approach is reported in Fig. 2) as orange circles versus the number of clusters $N_C$.

3. Discussion

In Fig. 2, it is possible to note that $\alpha$ obtained with the CMM based focusing (red disks) is larger than the one obtained with SO (gray squares), and very close to the theoretically expected average value of $\alpha _{theor}=1/\pi \simeq 0.31831$ See Supplement 1 for all values of $N_C$. The CCM, indeed, provides an experimentally retrieved focusing efficiency $\alpha _{exp}=0.3125\pm 0.0021$, which is $98 {\% }$ of the theoretical value (estimate obtained averaging the values from in Fig. 2 for all $N_C$). Note that the value of $\alpha _{theor}$ is obtained with an average-over-disorder operation, and as such, individual realizations of scattering can reach values higher than this theoretical limit. Indeed error bars in Fig. 2, are obtained extracting the statistical error on different disorder realization, and can surpass the theoretical, disorder averaged, value. We stress that the $\alpha$ obtained with CCM focusing, is independent of the number of clusters $N_C$, since a deterministic approach is expected to be unaffected by the roughness of the intensity landscape. Note that , at the best of our knowledge, previous methods to experimentally determine $\mathbf {t}$, related to the eigenvector of $J$, are sub-optimal (see Supplement 1), reaching from 40${\% }$ to 60${\% }$ of the theoretically expected $\alpha$, typically because they are strongly influenced by the reference patterns [12]. Lower panels of Fig. 4, report the intensity obtained at the end of the optimization process with SO (b) and with the eigenvector found by CCM (d). The two pictures are similar but the eigenvalue approach based on CCM provides optimized intensity always higher than SO. Figure 4(c) reports the predicted intensity obtained by extracting the eigenvalue and the eigenvector $\boldsymbol {\sigma }^{\rm max}$ from the measured coupling matrix and then using the formula

To further investigate the relation between roughness of the intensity landscape and focusing efficiency $\alpha$ obtained with the SO we plot in Fig. 5 the overlap (degree of similarity) calculated with the normalized scalar product $\boldsymbol \sigma ^{(r)} \cdot \boldsymbol {\sigma ^{\rm max}}$ between the optimized input vector resulting from a sequential optimization replica $\boldsymbol {\sigma }^{(r)}$ and the input profile obtained by eigenvalue retrieval $\boldsymbol {\sigma ^{\rm max}}$.

 

Fig. 5. Overlap between experimentally retrieved Focusing states with CM and SO versus number of Clusters NC. a) reports the overlap $\boldsymbol \sigma ^{(r)} \cdot \boldsymbol \sigma ^{\rm max}$ versus $N_C$, $\boldsymbol \sigma ^{(r)}$ being the SO configuration, $\boldsymbol \sigma ^{\rm max}$ being the eigenvector of the retrieved transmission channel by means of CCM. The orange circles are the average values (error bars are the statistical error), and smaller blue dots are individual measurements. Panel b) reports the probability density function of the intensity ($N$=225 at a single location) for SO (azure) and CCM (orange). The hand drawn pictures in panels c) and d) report a qualitative view of how could be imaged the Intensify landscape for the low NC and for the high NC case.

Download Full Size | PDF

It is possible to note how higher $N_C$ are correlated with sequential optimizations far from the optimal solution $\boldsymbol \sigma ^{\rm max}$. A sketch of the potential energy landscape representing the two configuration is shown in the bottom sketches. Figure 5(b)) reports the probability density function of the final intensity, obtained for the CCM (red) and for the SO. The graph demonstrates that the CCM retrieves a much higher intensity ($+28 {\% }$ over the field shown in Fig. 4(b-c)) and provides a smaller intensity variance. The hand-drawn pictures of the landscape generating these results is reported in Fig. 5(c)-d).

4. Summary

We demonstrated an unexpected parallelism between the memory retrieval process for Hopfield Networks and the Open channel activation in disordered media light transmission. In Neural Networks a memory state is defined by the coupling matrix, which is typically constructed as a dyadic matrix. In the case of diffusing light focusing, the role of coupling matrix is played by the $\mathbb{T}\mathbb{T}^\dagger$ matrix which is bi-dyadic due to the complex nature of the transmission matrix elements.

Thus both memory and open channel retrieval problems are solved by finding the correspondent eigenvectors, however, in the case of optics, two separate eigenvalues exist, producing a rough intensity landscape resulting in a non-convex problem. Our analysis showed that this roughness influences dramatically the SO, which being a dynamic process of exploration, can bring the segments configuration in to an inefficient, one-spin-stable “trap” with intensity lower than that obtainable employing eigenvector. On the other hand CCM retrieves the entire CM while being not subject to local traps.

CCM may be only delivered by employing the DMD technology, because of its speed and flexibility. This novel approach can be used to retrieve CM at each location of the transmitted plane, on the far side of an opaque medium with a single measurement, making fast laser scanning through walls a realizable goal.