Optical scattering randomly redirects the flow of light. It is a ubiquitous phenomenon that has wide-ranging effects. Since imaging relies on light travelling in straight lines from a scene to a camera, scattering prevents image formation through fog and precludes high-resolution microscopy inside biological tissue^1,2. Scattering also impairs optical communications through air and water, and disrupts the transmission of microwave and radio signals³. Overcoming the adverse effects of light scattering is an extremely challenging problem⁴. Nonetheless, owing to its prominence, substantial progress has been made over the past decades⁵.

When light propagates through a strongly scattering environment (also known as a ‘complex’ medium¹), the wavefront of the incident optical field is distorted, corrupting the spatial information that it carries. Elastic scattering from a static medium is deterministic, meaning that the precise way in which light has been perturbed can be characterized and subsequently corrected. By sending a series of probe measurements through the medium, a digital model of its effect on light can be created: represented by a linear matrix operator known as a transmission matrix (TM)⁶. The TM reveals how to pre-distort an input optical field so that it evolves into a user-defined state at the output—a technique known as wavefront shaping⁷.

Despite these successes, control of light through time-varying complex media remains a largely open problem². Evidently, wavefront shaping can only be successfully applied if the medium in question stays predominantly stationary for the time taken to make probe measurements and apply a wavefront correction. Yet many application scenarios feature complex media that rapidly fluctuate on a timescale of milliseconds or faster⁸—rendering wavefront shaping approaches exceedingly difficult. So far, the main strategies to control light through moving complex media have focused on achieving the task of wavefront shaping as quickly as possible^{8,9,10,11,12,13,14,15,16}. Approaches rely on ultra-fast beam shaping^{17,18,19,20,21,22,23} or reducing the number of probe measurements by, for example, spectral multiplexing^21,24 or exploiting prior knowledge about the scattering medium^{25,26,27,28,29}.

Here we introduce an alternative way to control the propagation of light through dynamic scattering media. We begin by classifying complex media into three categories, based on the type of motion exhibited over the timescale required for wavefront shaping, τ_ws. Class 1 represents static complex media that remain fixed over time τ_ws—established TM-based methods can be applied in this case. Class 2 represents moving complex media, which undergo substantial motion everywhere over time τ_ws and elude current wavefront shaping approaches. There is a third class of medium that falls between classes 1 and 2. Class 3 comprises partially moving scattering media, which, over the timescale τ_ws, exhibit localized time-varying pockets embedded within a comparatively static medium. This situation describes, for example, slowly moving tissue through which the small capillaries conducting blood flow typify faster changing regions, and pockets of turbulent air above hot chimneys within calmer air over a city skyline. In this article, we focus on how to identify light fields that thread through networks of static material within such partially moving complex media.

When a light field u is incident on a dynamic medium, the time-dependent transmitted field v(t) is given by

where T(t) is the time-dependent TM of the medium, and u and v are complex column vectors. Our aim is to find an input u that minimizes the temporal fluctuations in the output field v(t) by analysing externally scattered light.

To begin experimentally investigating this scenario, we emulate a three-dimensional time-varying forward scattering medium using a cascade of three computer-controlled diffractive optical elements, each separated by free space, as shown in Fig. 1. Cascades of phase planes can emulate atmospheric turbulence^30,31 and mimic multiple scattering samples^32,33,34,35. In practice, our set-up is implemented using multiple reflections from a liquid crystal spatial light modulator (SLM), allowing the phase profiles to be arbitrarily digitally reconfigured. We choose this test bed as it is a versatile way to control the number and location of dynamic regions for proof-of-principle experiments.

As shown in Fig. 1e (top row), we show a static random phase pattern on each phase screen, spatially distorting optical signals flowing through the system. On each plane, we also define an area within which the phase profile is programmed to randomly fluctuate in time. A second SLM is used to shape the incident light, and a digital camera records the level of temporal intensity fluctuations in transmitted light. We quantify these fluctuations by calculating the temporal fluctuation level ({F}_{l}={bar{sigma }}_{I}/{bar{mu }}_{I}), where ({bar{sigma }}_{I}) denotes the standard deviation of the time-fluctuating intensity, averaged over all illuminated camera pixels, and ({bar{mu }}_{I}) is the average transmitted intensity. For example, F_l = 0 for an unchanging transmitted field, and F_l = 1 for a series of uncorrelated Rayleigh speckle patterns (Methods). Here we use the term ‘fluctuations’ to refer to temporal fluctuations.

Unguided optimization

We first pursue a straightforward optimization method: iterative adjustment of the input field u to suppress the level of temporal intensity fluctuations (F_l) in the transmitted field. Our approach is shown schematically in Fig. 1a (Methods and Supplementary Section 1). Figure 1b shows a typical optimization convergence curve, which plateaus at F_l ~0.1. Noiseless simulations of this experiment reproduce a similar convergence plateau (Supplementary Section 8). We speculate that this may be because we optimize only the phase of the input field in this experiment. Figure 1c,d shows examples of the output temporal fluctuations of a randomly chosen initial trial field (F_l = 0.25) and an optimized field (F_l = 0.1), respectively. See also Supplementary Video 1, which shows that intensity fluctuations have been heavily suppressed.

The artificial nature of our dynamic scattering medium renders it possible to directly observe the evolution of the optimized field inside by digitally ‘peeling back’ the outer scattering layers (Supplementary Section 2). We find the optimized field takes the form of a speckle pattern that evolves to exhibit near-zero intensity on fluctuating regions on each plane (Fig. 1e, bottom row)—thus almost entirely avoiding these dynamic areas.

Our unguided optimization strategy is analogous to the first methods used to shape light through static scattering media⁷, and as such may be improved using more advanced algorithms^36,37. Furthermore, the form of the objective function can be arbitrarily chosen. For example, intensity shaping terms could potentially be included to simultaneously reduce fluctuations and shape the output^38,39. However, undirected optimization is a relatively slow process requiring many iterations to converge (~2,500 in our experiment). Therefore, we next ask: is there a way to find optimized fields more rapidly?

Physical adjoint optimization

In our first strategy, on each iteration, we measure how one single spatial component of the input field should be adjusted to reduce F_l. We now describe how to simultaneously measure how all spatial components should be adjusted in parallel. Our approach can be understood as gradient descent optimization using a physical realization of fast adjoint methods—which enable the efficient determination of the gradient of an objective function with respect to the optimization variables.

Our scheme is shown in Fig. 2a (Methods and Supplementary Section 4). To render an adjoint optimization approach feasible, it is necessary to modify the optimization objective function and send light back and forth in both directions through the dynamic scattering medium (reminiscent of work placing a scattering medium inside a laser cavity⁹). We now aim to maximize the correlation C between all N measured output fields over time, given by

See Supplementary Section 3 for derivation of this method and proof of its convergence.

Figure 2b shows the increase in the objective function throughout the optimization process. For comparison with our first strategy, we also monitor the level of temporal fluctuations F_l. After only ~15 iterations the procedure converges. See also Supplementary Video 2. Here we attribute the higher level of residual fluctuations to additional noise induced by the asynchronous uncontrolled flickering of the three SLMs used in this experiment. Supplementary Section 8 provides supporting simulations and discussion of the noise floor. Once again looking inside the dynamic medium, we see that we have found a more smoothly varying optical field that avoids the moving regions, as shown in Fig. 2c—an effect also reproduced in simulations (Supplementary Section 8).

This adjoint optimization strategy is reminiscent of iterative time reversal⁴⁰ and recently proposed in situ methods to train photonic neural networks⁴¹. Indeed, our work may be considered one of the first real-world implementations of a photonic adjoint optimization routine—a challenging yet powerful method to realize experimentally⁴².

Temporal fluctuation eigenchannels of the time-averaged TM

So far, we have focused on strategies to find a single stable channel. We now consider how to determine a set of stable channels that all navigate around dynamic regions inside the medium. Our approach makes use of the information stored in the time-averaged TM of a fluctuating optical system, 〈T〉_t. Figure 3a shows how 〈T〉_t is measured (Methods and Supplementary Section 5).

〈T〉_t reveals the input fields that deliver high levels of coherently time-averaged energy to the output. Finding such input fields can be represented as an eigenvalue problem by noting that the total intensity of the time-averaged field arriving at the output, 〈P〉_t, can be expressed as

Therefore, the eigenvectors of matrix ({langle {bf{T}}rangle }_{t}^{dagger }{langle {bf{T}}rangle }_{t}) with the largest eigenvalues represent input fields that deliver the highest intensity of the time-averaged output fields. These eigenvectors also correspond to input fields that interact least with the dynamic regions inside the medium, under the assumption that the internal fluctuations of the medium are large enough to randomize the phase of dynamically scattered light, meaning that fluctuating fields can be effectively ‘time-averaged away’ (Supplementary Section 12). We term this basis of eigenvectors the ‘temporal fluctuation eigenchannels’ of the dynamic medium. This concept links to our previous approach: we show in Supplementary Section 3 that when using the objective function given in equation (2), physical adjoint optimization finds the most stable temporal fluctuation eigenchannel of the time-averaged TM.

Figure 3b shows the distribution of absolute eigenvalues of ({langle {bf{T}}rangle }_{t}^{dagger }{langle {bf{T}}rangle }_{t}), arranged in ascending order. We compare two 3-plane dynamic samples with different numbers of moving regions: (i) has a single dynamic patch on each plane similar to before; (ii) has randomly placed fluctuating patches covering approximately half of the area of each plane. We show in Supplementary Section 11 that medium (i) has at least one stable temporal fluctuation eigenchannel that avoids all dynamic patches, while medium (ii) has no completely stable channels.

We first demonstrate excitation of the fluctuation eigenchannels through the more weakly fluctuating sample, shown in Fig. 3d. We see that the transmitted fields corresponding to high eigenvalues remain stable. Conversely, the transmitted fields corresponding to low eigenvalues vary with time—as these modes interact strongly with the moving parts of the medium. See also Supplementary Video 3.

We now turn our attention to the more strongly fluctuating medium. Figure 3e shows the stability of transmitted fields when exciting channels with high and low eigenvalues. As expected, even light propagating through the most stable eigenchannel exhibits non-negligible output fluctuations over time, indicating that we have not found any fields that thread perfectly around all dynamic regions. Despite this, we find that a marked improvement in focusing at the output is possible using the information stored in the time-averaged TM. Figure 3f shows a focus created using conventional wavefront shaping, where the medium freely fluctuates throughout TM measurement (left column), compared with focusing using a sub-basis formed from the top 100 most stable temporal fluctuation eigenchannels (right column)—see Supplementary Section 6 and Supplementary Video 4 for details. Both the contrast and stability of the focus are enhanced. Beyond focusing, more elaborate beam shaping may also be possible using this stable sub-basis⁴³.

The time-averaged TM is related to several previously introduced matrix operators connected to physical quantities in scattering media^{6,44,45,46,47}. In particular, it has similar properties to the TM of a static medium with inhomogeneous absorption^39,48,49. As fluctuations in the medium go to zero, the time-averaged TM tends to the conventional TM, and the temporal fluctuation eigenchannels tend to the transmission eigenchannels of a static scattering medium^50,51. The ‘deposition matrix’⁴⁷ and the ‘generalized Wigner-Smith operator’^44,52 are both also capable of revealing light fields that circumnavigate predetermined regions within a complex medium. However, only the time-averaged TM does so without requiring access to internal fields within the medium⁴⁷ or the measurement of an entire TM while the medium is held static⁴⁴.

Stable light transmission through flexing optical fibre

Up to this point, we have considered samples with well-defined moving regions embedded within static complex media. We now explore a more general scenario in which an entire medium deforms, leaving no clearly identifiable locations that remain static. To investigate this case, we search for stable channels within a step-index multimode fibre (MMF) as it is gradually bent—here all transmitted spatial modes will interact with the core-cladding boundary of the fibre as it flexes.

Figure 4a shows a schematic of our experiment (Methods and Supplementary Section 7). We measure the time-averaged TM of an MMF supporting ~750 spatial modes as it is smoothly moved through 9 different bend configurations. Figure 4b shows how the TM decorrelates as the fibre is bent. Transmitting a fixed random input field through the fibre yields different output speckle patterns for each bend state, corresponding to a temporal fluctuation level of F_l ~0.8.

Figure 4c shows the eigenvalues of matrix ({langle {bf{T}}rangle }_{t}^{dagger }{langle {bf{T}}rangle }_{t}). We experimentally measure the temporal fluctuations F_l in the output fields of 754 temporal fluctuation eigenchannels, which are shown in Fig. 4d. We find ~150 channels that exhibit lower fluctuation levels than that of a random input field, across this range of fibre movement. Figure 4e shows examples of the most stable channels—all of which show reduced intensity on the axis (see Supplementary Section 7 for discussion of this phenomenon). See also Supplementary Video 5. This approach may be useful for the transmission of data through MMFs.

We have shown that it is possible to identify and guide light through stable channels within dynamic scattering media. Our methods do not rely on prior knowledge of the location of dynamic regions or type of motion, and only require measurements of externally scattered light made on the same timescale as temporal fluctuations. To understand the potential applications of these techniques, a key question arises: to what extent do such stable channels exist within dynamic scattering media? While there is a large parameter space of possible scattering scenarios, we numerically examine some general trends for random disordered media.

We first use a layered forward scattering model to study how the proportion of stable transmission channels scales as a function of the area fraction of dynamic regions on each layer, and the number of randomly connected layers (that is, medium depth)—see Supplementary Section 11 for details. For the first few layers, we find that the number of stable channels decreases linearly with the dynamic area per layer and the number of layers. However, at greater depths, this trend becomes sub-linear, and a proportion of stable channels are capable of propagating deeply. For example, with a 5% dynamic area ratio per layer, we find that ~10% of the modes remain highly stable to a depth of 20 layers. As with all wavefront shaping approaches, our ability to excite stable channels will depend upon the number of SLM control elements relative to the mode capacity of the medium.

We next consider media in which backscattering is not negligible. Such media pose additional challenges owing to competing requirements: optical fields must both circumnavigate moving regions and also penetrate deeply enough into the medium to transmit substantial levels of power to the other side. In Supplementary Section 13, we explore how stable channels may be identified in this regime. We first show that if we are able to measure the full scattering matrix of a diffusive sample, coherently time-averaging the scattering matrix allows stable channels to be identified. Since the full scattering matrix is difficult to measure in most realistic scenarios, we next show how our physical adjoint optimization approach can be adapted to successfully guide light around dynamic regions while measuring transmitted light only. This modified optimization procedure requires light to be passed back and forth once per iteration on a timescale over which the entire medium is static (the medium is still free to move between iterations). This places more stringent constraints on the update rates of SLMs, yet offers an advantage over conventional wavefront shaping methods in which the medium must remain static throughout the entire optimization process. More broadly, we note that each new choice of adjoint optimization objective function will lead to a new physical optimization procedure, and it will be interesting to further explore this space in the future.

We now evaluate our adjoint methods in the context of entirely dynamic media containing no permanently static regions, but exhibiting multiple decorrelation rates ranging from τ_fast to τ_slow. If coherent scattered field averaging is carried out over a timescale longer than τ_slow, this average tends to zero and our methods fail. However, we find that as few as two medium realizations will suffice to create useful coherent averages to guide optimizations towards a more stable channel (Supplementary Sections 9, 10 and 12). Thus, if optimization is conducted on a timescale τ_opt that is τ_fast < τ