Shaping the branched flow of light through disordered media.

Electronic matter waves traveling through the weak and smoothly varying disorder potential of a semiconductor show a characteristic branching behavior instead of a smooth spreading of flow. By transferring this phenomenon to optics, we demonstrate numerically how the branched flow of light can be controlled to propagate along a single branch rather than along many of them at the same time. Our method is based on shaping the incoming wavefront and only requires partial knowledge of the system's transmission matrix. We show that the light flowing along a single branch has a broadband frequency stability such that one can even steer pulses along selected branches-a prospect with many interesting possibilities for wave control in disordered environments.

When waves propagate through a disorder landscape that is sufficiently weak and spatially correlated, they form branched transport channels in which the waves’ intensity is strongly enhanced. This phenomenon of “branched flow” was first discovered for electrons gliding through semiconductor heterostructures (1). Instead of an isotropic spreading into all possible directions, the electron density injected through a quantum point contact was observed to form esthetically very appealing branch patterns. This intriguing behavior can be attributed to ripples in the background potential that are always present in such structures (1, 2), which act like an array of imperfect lenses, giving rise to caustics (3) and thereby, to distinct intensity enhancements along branches (4–6). Although first discovered as a nanoscale wave effect, branched flow was soon understood to occur on a wide range of length scales up to the formation of hot spots in tsunami waves as a result of the propagation through the rough ocean sea bed (7–13).

Whereas a number of previous studies have already focused on the statistics of this phenomenon (3, 10, 14) and on its origins (2, 3, 15–18), the question of how branched flow can be controlled and thereby, put to use for steering waves through a complex medium has not been addressed so far. This is probably due to the fact that the possibilities to shape and manipulate electrons or ocean waves are, indeed, very limited. In other words, for the experiments where branched flow was observed so far, the incoming wavefront as well as the potential that the wave explores were considered as predetermined and immutable. These limitations are currently about to be overcome in a new generation of experiments, where coherent laser light was observed to exhibit branching when propagating through very thin disordered materials, such as the surface layer of a soap bubble (19). Specifically, we expect that the transfer of branched flow to the optical domain will open up the whole arsenal of photonics to shape the wavefront of such branched light beams (20, 21).

A particularly exciting question that we will explore here from a theoretical point of view will be whether optical wavefront shaping tools, like spatial light modulators (SLMs), can be used to manipulate an incoming light beam in such a way that it follows only a single branch through a disorder landscape rather than many of them in parallel. The protocol that we will introduce based on our analysis will open up ways of sneaking a beam of light through a disordered medium while maintaining its focus throughout the entire propagation distance—like a highway for light through a scattering medium. More generally, we expect our approach to be useful in a variety of different contexts, where steering waves through a complex environment to a predetermined target is a key goal, like in wireless communication (22), adaptive optics (23), underwater acoustics (24), wave focusing (25–28), and biomedical imaging (29, 30) as well as for wave control in disordered systems at large (20, 21, 31, 32).

Results

System.

The system that we consider is shown in Fig. 1 and consists of a rectangular scattering region of length and width that is attached to two straight semiinfinite waveguides (leads) of the same width on the left and right (only the left lead is shown). In transverse direction, hard-wall boundary conditions are applied (i.e., the wavefunction is zero at these boundaries). In all of the calculations reported below, we choose the number of propagating open lead modes to be and a fixed wavenumber of the incoming light (the first three lead modes are indicated in Fig. 1). For simplicity, we set , resulting in the following simple expressions for the wavenumber and the wavelength . The length of the scattering region is chosen as .

Fig. 1.

(A) Illustration of the setup under study: a disordered waveguide (center region) of width is attached to an incoming lead on the left (through an aperture; blue) and an outgoing lead on the right (not shown). The smooth disorder potential in the waveguide of length is illustrated in blue/red colors. On top of the potential, the wave intensity is plotted, corresponding to a superposition of the wave intensities of modes 1–100 injected through the left lead and resulting in a pronounced branched structure. (The first three lead modes are depicted in orange on the left-hand side.) The main goal of our study is to separate these branches by suitably shaping the incoming wavefront in the left lead. (B) Intensity output profile as a function of transverse coordinate at the right end of the disordered region (). The seven maxima labeled from one to seven (highlighted in light blue/magenta) are produced by the arrival of different branches at the output.

In analogy to the first observation of branched flow, where electrons were injected through a constriction (quantum point contact) into a high-mobility electron gas (1), we also include such a constriction in the form of an aperture of width between the left lead and the disordered scattering region (at ). Whereas in many previous studies, the width of the constriction was chosen such that it only allows for 1 or 2 modes to propagate, the 50 modes that we allow to pass provide us with many tunable degrees of freedom as required for shaping the incoming wavefront (for a narrow quantum point contact that transmits only a single mode, our method cannot be applied). The smooth and long-range disorder necessary to observe branched flow is modeled by a spatially dependent index of refraction throughout the whole scattering region indicated by the light red/blue color in Fig. 1. This correlated refractive index is characterized by a correlation length , minimum/maximum values and , and a mean value (more details on the correlated disorder are in ).

The scalar scattering problem in this two-dimensional setup is described by the two-dimensional Helmholtz equationwith representing the out-of-plane component of the electric field and being the incoming wavenumber. To solve this equation numerically, we discretize space on a Cartesian grid and use the modular recursive Green’s function technique (33, 34) ( has more details) to efficiently evaluate the scattering states and the unitary scattering matrixHere, the transmission (reflection) matrix () contains the complex amplitudes () for transmission (reflection) from mode from the left lead to mode in the right (left) lead. The primed quantities and (not used here) are the corresponding matrices for injection from the right. (More details on the used scattering formalism can be found in .)

To observe the branched flow of light, we inject the different lead modes from the left into the constriction and superimpose the corresponding wave intensities that they give rise to. In the superposition, we consider only the first 100 lead modes (of 200) to avoid high-angle scattering and to ensure a high visibility of the individual branches. The branched structure in the propagation of waves through our setup is clearly visible in Fig. 1.

The challenge that we rise to in the next step is to address these branches individually through a suitable coherent superposition of incoming modes in the left lead. The methods that we choose for this purpose involve only the transmission matrix from Eq. , which is available in optics through interferometric measurements involving an SLM (21, 35). As the branched flow in our system naturally leads to a concentration of intensity at certain spots at the output, we find here that the knowledge of the transmission matrix for modes concentrated around these spots is sufficient for a clean separation of branches. In other words, we may restrict ourselves to those regions in space at the output where the branches arrive. These regions are determined from the intensity profile at the output facet of our system at (Fig. 1). Seven intensity maxima corresponding to the arrival of different branches are clearly visible in Fig. 1 and highlighted in light blue/magenta. For each intensity maximum, we manually set lower and upper boundaries, which are indicated by vertical lines in Fig. 1, and define a reduced transmission matrix connecting the incoming lead with the corresponding region at the output (labeled from one to seven). The elements of this matrix hence describe the coherent transmission amplitudes from all points at the input (we choose 200 equidistant points in the input lead corresponding to 200 open lead modes) to all points around a specific intensity maximum at the output ( has more details on the transmission matrix in a spatial basis).

First Approach: Transmission Eigenstates.

Our first approach to achieve clean branch separation is to use a singular value decomposition of : that is,allowing us to access the transmission eigenvalues in this truncated spatial basis as contained in the real diagonal matrix . The matrix consists of eigenvectors of , and consists of eigenvectors of . The largest transmission eigenvalues correspond to those transmission eigenchannels (contained in the columns of ) that transmit the most intensity to the desired region at the output. At this point, one might be tempted to think that these highly transmitting channels will already constitute the branches that we are after. To test this hypothesis, we inject for each of the seven transmission matrices the corresponding transmission eigenchannels with the largest transmission eigenvalues. Checking the corresponding scattering wavefunctions (Fig. 2 ), we find that a clean branch separation is, indeed, possible for a number of cases. We also find, however, that several among the highly transmitting eigenchannels follow two different branches in parallel instead of only one. Fig. 2 shows an example of such a state, where one can clearly see a mixing of one branch propagating directly into the selected region (marked with blue bars at the output) with another branch bouncing off the lower boundary. Demanding high transmission into a desired region by choosing high-transmission eigenchannels is thus clearly not enough to guarantee clean branch separation, since high transmission can also be obtained by propagating along multiple branches at once. In a possible optical experiment, such a mixing can be expected to be even more prevalent than in our numerical example, simply because optical implementations can typically involve a large number of branches (19).

Fig. 2.

Transmission eigenstates (B, E–I, O, P, T) and time-delay eigenstates (A, C, D, J–N, Q–S, U) of the system shown in Fig. 1 calculated from the transmission matrices connecting the incoming lead with the seven different output regions marked in Fig. 1 (see blue horizontal bars here). The branched structure (Fig. 1) is drawn here as a gray background. We show all eigenstates that are selected by our procedure (steps i–iii) with black text labels and two examples of eigenstates that are filtered out by our procedure, since they address two different branches at once (red labels). The states are ordered according to their output region on the right (numbers ①–⑦). The average transmittance of all states with black labels is around into the region marked by the blue bars.

Second Approach: Time-Delay Eigenstates with Large Transmission.

To also be able to address such mixed branches individually, we now introduce a more efficient method. Specifically, our aim is to set up an approach in terms of the scattering matrix and the Wigner–Smith time-delay operator (36–38) derived from it:Eigenstates of , also known as principal modes, are associated with scattering states that have a well-defined time delay and the remarkable property that their output profile is very robust to frequency changes (39–44). Some of these eigenstates have the additional feature of having a particle-like wavefunction (i.e., the scattering states follow classical particle trajectories) (40, 42, 44). Modifying the Wigner–Smith time-delay operator for our purpose now allows us to separate those eigenstates of with the largest transmission eigenvalues by their time delay (27, 44). The key idea here is that two branches that may both be highly transmitting (like those in Fig. 2) can be distinguished by their different time delays (as determined by the different branch lengths). To be specific, we only work with those transmission eigenvalues that are larger than some value and derive the matrices , , and from , , and by truncating all rows and columns corresponding to . With these truncated matrices, we can now replace the terms in Eq. ,to arrive at the reduced time-delay operator ,that operates in the subspace of highly transmitting states only. Note that Eq. involves a quasi-inverse “” of the rectangular matrix , the regularity of which is guaranteed by the restriction to only those transmission eigenvalues that are larger than the cutoff value . In practice, a value of proved suitable for all of our calculations. Note that, due to the nonunitarity of , the eigenvalues of the reduced time-delay operator in Eq. are complex (in contrast to the real eigenvalues of the Wigner–Smith time-delay operator in Eq. ). The imaginary parts of the complex eigenvalues are, however, very small, and the real parts can still be used as a good measure for the physical delay times (44).

To put this method directly to the test, we turn our attention to the state shown in Fig. 2 featuring a mixture of two branches with different path lengths and correspondingly, different time delays. A singular value decomposition of reveals that it contains nine singular values larger than . With this knowledge, we can now construct according to Eq. and indeed, find among its eigenstates the desired wave fields that follow the two involved branches individually (Fig. 2 ).

Restricting the construction of time-delay eigenstates to the subspace of high transmission thus yields already very good results. Using this method, we, however, also observed a few time-delay eigenstates that mix two different branches as, for example, shown in Fig. 2. These two branches, however, turn out to be individually addressable through those transmission eigenstates with the smallest time delays (Fig. 2 ) [the time delay of a transmission eigenstate can simply be calculated by taking the expectation value with the time-delay operator (i.e., )].

Combined Method.

One may thus also decide to turn the above strategy on its head and look for transmission eigenstates in the subspace of short time delays. Since neither one of these opposite strategies seems to have an a priori advantage, we now combine them with each other in a synergistic way to improve our results even farther. (i) We evaluate all eigenstates of and of for all of the seven transmission matrices corresponding to regions of maximum intensity shown in Fig. 1. (ii) We select those states that are identical in both eigenstate sets, since they turn out to be individual branch excitations in all of the observed cases. To do this, we project each eigenvector onto each eigenvector such that we end up with the matrix elements . (The matrix is not unitary, since the eigenstate sets are not complete.) For the case that two eigenvectors are the same, the matrix has only one significant nonzero element in the corresponding row/column. Practically, we consider two eigenvectors to be the same when . (iii) In a last step, we deal with those eigenstates that consist of more than one contribution from the respective other eigenstate set (i.e., that have more than one nonzero element in the corresponding row/column of ). Our task here is to select those states that consist of only single branches and to discard those states that propagate along more than one branch at once. Since, however, the coefficients do not indicate per se which states consist of single branches only, we first need to translate the eigenvector coefficients to a corresponding angular profile at the input aperture. Checking, in a next step, if this angular input pattern is sufficiently collimated provides us finally with the desired indicator for the excitation of a single branch ( has details).

Following the above three steps (i–iii), which notably rely only on the experimentally accessible transmission matrices , we obtain well-separated branch states (Fig. 2 ) that stay collimated throughout the entire scattering region and that feature an average transmittance of over into one of the designated seven output regions. These results show that our method leads to a channeling of waves through the disordered region and to a well-controlled branched flow. An interesting detail that we emphasize here is that our approach not only yields a single state for each individual branch but in fact, also states that propagate along the same branch but with a higher transverse quantization (Fig. 2 ) (40).

Discussion and Summary

Injection into Empty Cavity.

To underscore the nontrivial nature of these collimated branch states that we identify here, we inject several of the states shown in Fig. 2 into a clean waveguide without any disorder. The results are displayed in Fig. 3, showing that these states feature a considerably reduced collimation compared with the case including the disorder (Fig. 3, Insets). This observation demonstrates that the states that we identify here do not just rely on a trivial injection with a narrow angular distribution at the input and that the disorder plays a crucial role for the states’ collimation.

Fig. 3.

(A–D) Eigenstates from Fig. 2 injected into an empty waveguide without the disordered refractive index. The fact that removing the disorder potential leads to a defocusing demonstrates that the formation of collimated branch states crucially relies on the presence of the underlying disorder landscape. The blue bars indicate the region where the branch exits the scattering region in the presence of the disorder (Insets).

Pulse Propagation.

In the last part of this study, we also demonstrate explicitly that our collimated single-branch states can be sufficiently stable in frequency to allow for the transmission of pulses along a branch. Consider here, as an example, the time-delay eigenstate shown in Fig. 4 that propagates along a certain branch. Taking a superposition of this branch state at different frequencies to form a Gaussian wave packet, we obtain a pulse propagating along the selected branch as shown in Fig. 4 at three different time steps (). We observe that the pulse transits the system while staying on this curved branch throughout the entire transmission process.

Fig. 4.

(A) Time-delay eigenstate propagating along one single branch in a disorder landscape (to prove the general applicability of our approach, a different disorder realization was used as in Fig. 1). (B–D) Pulse propagating along the branch shown in A at three different time steps (). The pulse remains spatially confined while traversing the disorder along the branch shown in A. The Fourier spectrum of the pulse is Gaussian shaped with an SD of , with being the wavenumber.

Summary.

In summary, this work demonstrates how to control the flow of waves through a correlated and weak disorder potential landscape. Such systems give rise to branches along which incoming waves travel through the disorder. We introduce a method that allows us to inject waves in such a way that almost all of the flow travels along a single branch alone. This nontrivial finding can even be extended to the temporal domain as we show by creating pulses that remain on a single branch throughout the entire transmission process. Implementing such concepts in optics requires only a small subpart of the transmission matrix and is thus within reach of present-day technology. We expect our work to be generalizable from scalar to vector waves and from two to three dimensions, where it may give rise to interesting applications in communication and imaging technology.

Methods

Numerical Method.

To solve the Helmholtz Eq. numerically, we discretize the central scattering region on a Cartesian grid with a grid spacing that is about a factor 20 smaller than the considered wavelength . This discretization allows us to reformulate the scattering problem as a standard matrix equation involving very large matrices. The central quantity in our approach is the so-called Green’s function that contains the information on how any incoming wave produces a certain wave pattern inside the entire waveguide. We calculated this Green’s function through an efficient “modular” approach (33, 34) that involves the inversion of large matrices on a computer cluster using efficient linear algebra packages.

Scattering Formalism.

Our scattering system is a waveguide as described. The semiinfinite asymptotic regions (also called leads) feature a constant refractive index of , whereas the rectangular scattering region in the middle contains the randomly distributed but correlated refractive index . Solving the Helmholtz Eq. in the asymptotic regions yields straightforward solutions—the so-called lead mode basis functionswhere is the wavenumber in the direction of mode and is a flux-normalization factor. The transverse profile of each mode is given bywith and being the lead width. The transverse and longitudinal wavenumbers follow the relation . The first three lead modes () are depicted in Fig. 1 in orange. The lead modes in Eq. are orthogonal and complete, and they can thus be used as a basis to decompose any arbitrary wave such that the wave can be described by a corresponding coefficient vector , where the upper components represent a wave injected from the left lead and the lower components represent a wave injected from the right lead. For our simulations, we exclusively work with waves that approach the scattering region from the left [i.e., ]. Due to scattering, the wave injected through the aperture will be different after passing the scattering region [i.e., the outgoing wave is described by a different vector , with being the reflected wave and being the transmitted wave]. Incoming and outgoing waves are related by , with being the scattering matrixwhere are the reflection and transmission matrices containing the scattering amplitudes for injection from the left-hand side and consist of the corresponding amplitudes for injection from the right-hand side. As an example, the matrix element is the complex amplitude for the transmission from incoming mode to outgoing mode . If is the number of propagating modes in one of the two leads, the matrices are dimensional, and correspondingly, the scattering matrix is dimensional. When applying this formalism to the reduced time-delay operator in Eq. , the coefficient vector that corresponds to an eigenstate of this operator reads . Here, is the th eigenvector obtained by solving the eigenvalue equation for the reduced time-delay operator (i.e., , where is the corresponding eigenvalue).

Spatial and Angular Profile of Eigenstates.

To find individual branch excitations among all eigenstates and , it is essential to determine if either a time-delay eigenstate or a transmission eigenstate addresses only one single branch rather than many at the same time. As we show here, the spatial and/or angular distribution of an eigenstate at the input aperture provides us with sufficient information to perform this task, since an eigenstate exciting only one branch will be spatially more confined and will radiate into a smaller angular region than a state addressing multiple branches. Assuming that the transmission matrix is measured in the spatial pixel basis, the eigenvectors and are naturally given in this spatial basis as well. By plotting the absolute value of the coefficients , where is the th component of the vector or , as a function of the transverse coordinate at the aperture (), we can easily generate the spatial distribution of an eigenstate.

To estimate the angular distribution of an eigenstate at the aperture, we work with the Hermitian operator measuring the transverse component of the wavevector. The eigenvalue equation of the th eigenvector of this operator readswhere is the th eigenvalue. Since a well-defined transverse wavevector component corresponds to a well-defined angle of incidence, we can now decompose the eigenvectors and into the momentum basis spanned by the vectors and analyze the different angular components .

Fig. 5 display the spatial and angular components of the transmission eigenstate shown in Fig. 2 (red) and the time-delay eigenstate in Fig. 2 (blue). We see that the spatial profile of the transmission state is broader and that it features more angular components than the time-delay state. We can, therefore, conclude that the transmission state is more likely to address multiple branches, whereas the time-delay state addresses only one single branch, which is confirmed by the wave plots shown Fig. 2 . In Fig. 5 , we plot the same distributions for the time-delay state shown in Fig. 2 and the transmission state shown in Fig. 2. From Fig. 5, we deduce that the time-delay state consists of more than one branch due to the larger spatial distribution, which is confirmed by the wave plots. We successfully applied this procedure to all eigenstates and , from which we can conclude that the spatial and angular distributions of the time-delay and transmission states can be used to find those states out of both eigenstate sets (time-delay and transmission eigenstates) that excite only one single branch.

Fig. 5.

(A and B) Spatial (left) and angular distribution (right) at the input aperture (located at between and ) of the transmission (red) and time-delay states (blue) shown in Fig. 2 , respectively. The different widths of the distributions indicate that the transmission state is likely to excite more than one branch at once. (C and D) Spatial (left) and angular distribution (right) of the time-delay (blue) and transmission (red) states shown in Fig. 2 , respectively. From the widths of the spatial distributions shown in C, we can conclude that the transmission state is more likely to excite only one single branch, which is confirmed by the wave plots. The widths of the normalized distributions are quantified by the interval around the maximum value of the distribution (indicated by the vertical dashed lines) in which 60% of the distribution lies.