Focusing inside Disordered Media with the Generalized Wigner-Smith Operator

TL;DR

This paper presents a wavefront shaping method using a generalized Wigner-Smith operator to focus inside disordered media, leveraging the scattering matrix and its derivative, enabling targeted focusing without optimization.

Contribution

It introduces a novel wavefront shaping protocol based on a generalized Wigner-Smith operator, applicable to disordered media, demonstrated experimentally in the microwave regime.

Findings

Eigenstates sorted by focusing strength

Effective focusing without optimization or phase-conjugation

Potential applications in optical imaging

Abstract

We introduce a wavefront shaping protocol for focusing inside disordered media based on a generalization of the established Wigner-Smith time-delay operator. The key ingredient for our approach is the scattering (or transmission) matrix of the medium and its derivative with respect to the position of the target one aims to focus on. A specifc experimental realization in the microwave regime is presented showing that the eigenstates of a corresponding operator are sorted by their focusing strength - ranging from strongly focusing on the designated target to completely bypassing it. Our protocol works without optimization or phase-conjugation and we expect it to be particularly attractive for optical imaging in disordered media.

Figures

Equations

$$Q{\vec{u}}_n=-i{S}^{-1}\frac{\mathrm{d}S}{\mathrm{d}\omega }{\vec{u}}_n={\theta }_n{\vec{u}}_n,$$

$${\frac{\mathrm{d}{\vec{v}}_n}{\mathrm{d}\omega }|}_{{\omega }_0}=i{\theta }_n{\vec{v}}_n\left({\omega }_0\right)\to {\vec{v}}_n\left({\omega }_0+\Delta \omega \right)\approx e^{i{\theta }_n\Delta \omega }{\vec{v}}_n\left({\omega }_0\right).$$

$${\vec{u}}^{†}{Q}_{\alpha }\vec{u}={\vec{u}}^{†}{C}_{\alpha }\vec{u}-{\vec{v}}^{†}{C}_{\alpha }\vec{v}={⟨C_{\alpha }⟩}_{\mathrm{in}}-{⟨C_{\alpha }⟩}_{\mathrm{out}},$$

$$q_{y_{\mathrm{b}}}:=-i{t}^{-1}\frac{\mathrm{d}t}{\mathrm{d}{y}_{\mathrm{b}}},$$

$${\vec{v}}_{n,R}\left(y_{\mathrm{b}}+\Delta y_{\mathrm{b}}\right)\approx e^{(i{\beta }_n^{y_{\mathrm{b}}}-{\kappa }_n^{y_{\mathrm{b}}})\Delta y_{\mathrm{b}}}{\vec{v}}_{n,R}\left(y_{\mathrm{b}}\right).$$

$$[C_{\alpha }{]}_{mn}={\left[-i{\int }_{\partial \Omega }d{\xi }^D{\psi }_m^{\star }(x,\vec{\xi })\frac{\mathrm{d}{\psi }_n(x,\vec{\xi })}{\mathrm{d}\alpha }\right]}_{x=x_s},$$

$$S(\alpha +\Delta \alpha )=e^{-i{C}_{\alpha }\Delta \alpha }S(\alpha )e^{i{C}_{\alpha }\Delta \alpha }.$$

$$\begin{array}{rl}S(\alpha +\Delta \alpha )& \approx (1-i{C}_{\alpha }\Delta \alpha )S(\alpha )(1+i{C}_{\alpha }\Delta \alpha )\\ & \approx S(\alpha )-i{C}_{\alpha }S(\alpha )\Delta \alpha +iS(\alpha )C_{\alpha }\Delta \alpha ,\end{array}$$

$$\frac{\mathrm{d}S}{\mathrm{d}\alpha }=\underset{\Delta \alpha \to 0}{\lim }\frac{S(\alpha +\Delta \alpha )-S(\alpha )}{\Delta \alpha }=iS(\alpha )C_{\alpha }-i{C}_{\alpha }S(\alpha ).$$

$$Q_{\alpha }=-i{S}^{†}\frac{\mathrm{d}S}{\mathrm{d}\alpha }=C_{\alpha }-S^{†}{C}_{\alpha }S.$$

$$⟨Q_x⟩={\vec{u}}^{†}{Q}_x\vec{u}={\vec{u}}^{†}{k}_{\mathrm{in}}\vec{u}-{\vec{v}}^{†}{k}_{\mathrm{out}}\vec{v},$$

$$Q_x=\sum _{i=1}^{i_{\mathrm{tot}}}-i{S}^{†}\frac{\mathrm{d}S}{\mathrm{d}{x}_i}:=\sum _{i=1}^{i_{\mathrm{tot}}}{Q}_{x_i}=C_x-S^{†}{C}_{x}S,$$

$$⟨Q_{x_i}⟩={\vec{u}}^{†}{Q}_{x_i}\vec{u}={\vec{u}}^{i†}{k}_{\mathrm{in}}{\vec{u}}^i-{\vec{v}}^{i†}{k}_{\mathrm{out}}{\vec{v}}^i=\Delta k^i,$$

$$k_{\mathrm{in}}=\left(\begin{array}{cc}k& 0\\ 0& -k\end{array}\right)=-k_{\mathrm{out}},$$

$$⟨Q_{x_i}⟩=$$

$$=$$

$$Q_{x_i}=\left(\begin{array}{cc}{Q}_{x_i}^{11}& Q_{x_i}^{12}\\ Q_{x_i}^{21}& Q_{x_i}^{22}\end{array}\right),$$

$${\vartheta }^{x_i}=\frac{1}{|t{|}^2}\left(\Delta k^i-⟨i{r}^{†}\frac{\mathrm{d}r}{\mathrm{d}{x}_i}⟩\right),$$

$${\vartheta }^{x_i}=\frac{\mathrm{d}{\phi }_t}{\mathrm{d}{x}_i}-i\frac{\mathrm{d}\ln \left(|t{\vec{u}}_L|\right)}{\mathrm{d}{x}_i},$$

$$\mathrm{Im}({\vartheta }^{x_i})=-\frac{\mathrm{d}\ln \left(|t{\vec{u}}_L|\right)}{\mathrm{d}{x}_i}=-\mathrm{Im}\left(\frac{1}{|t{|}^2}⟨i{r}^{†}\frac{\mathrm{d}r}{\mathrm{d}{x}_i}⟩\right),$$

$$⟨Q_{x_i}⟩$$

$$=\sum _{n=1}^N{\left|(r{\vec{u}}_L{)}_n\right|}^2\frac{\mathrm{d}{\phi }_{r,n}}{\mathrm{d}{x}_i}+\sum _{n=1}^N{\left|(t{\vec{u}}_L{)}_n\right|}^2\frac{\mathrm{d}{\phi }_{t,n}}{\mathrm{d}{x}_i},$$

$$\mathrm{Re}({\vartheta }^{x_i})=\frac{\mathrm{d}{\phi }_t}{\mathrm{d}{x}_i}=\frac{1}{|t{|}^2}\left(\Delta k^i-\sum _{n=1}^N{\left|(r{\vec{u}}_L{)}_n\right|}^2\frac{\mathrm{d}{\phi }_{r,n}}{\mathrm{d}{x}_i}\right).$$

$${\tilde{q}}_{x_i}=-i\tilde{V}({\tilde{U}}^{†}t\tilde{V}{)}^{-1}{\tilde{U}}^{†}\tilde{U}{\tilde{U}}^{†}\frac{\mathrm{d}t}{\mathrm{d}{x}_i}\tilde{V}{\tilde{V}}^{†}$$

$$c_s(x,y)=s\ast K(x,y)=\sum _{x^{\prime },y^{\prime }}K(x^{\prime },y^{\prime })s(x-x^{\prime },y-y^{\prime }),$$

$$K=\frac{1}{4}\left(\begin{array}{ccc}0.25& 0.5& 0.25\\ 0.5& 1& 0.5\\ 0.25& 0.5& 0.25\end{array}\right),$$

Full text

One of the most formidable challenges for imaging in complex environments is to overcome the limitations imposed by the presence of disorder. In particular, if the wave scattering induced by a disordered medium is strong enough to suppress the ballistic contribution in the imaging process entirely, seeing through this medium or focusing on a target inside of it becomes a highly non-trivial exercise – a difficulty that is particularly evident in the field of biological and medical imaging. A promising new approach to image and focus also in the regime of multiple scattering is to exploit the information stored in a system’s scattering matrix Popoff et al. (2010a); Yu et al. (2013). In this emerging new field of “wavefront shaping” Mosk et al. (2012); Rotter and Gigan (2017) spectacular advances have recently been made, such as to focus light behind an opaque layer Vellekoop et al. (2010); Vellekoop and Mosk (2007); Katz et al. (2011); McCabe et al. (2011); Aulbach et al. (2012a), or to send and retrieve images across it Popoff et al. (2010b); Mounaix et al. (2016a, b); Paudel et al. (2013); Aulbach et al. (2011). Thanks to these advances also the focusing of light inside highly disordered media could recently be demonstrated using embedded fluorescent probes and nano-crystals Vellekoop et al. (2008); Aulbach et al. (2012b); Chaigne et al. (2014) or using digital optical phase conjugation to focus light onto a target moving inside an otherwise static environment Xu et al. (2011); Ma et al. (2014); Judkewitz et al. (2013); Zhou et al. (2014).

Here we present a new approach for focusing inside a disordered material that has the considerable advantage of working without the requirement to implant a fluorescent body at the focus or to phase-conjugate a wave scattered at the focus position. Our technique also allows us to tune the degree of focus on a designated position inside the disorder, including the case where the target is entirely avoided by the scattered wavefront. Our starting point for achieving this goal is the time-delay operator $Q$ introduced by Eugene Wigner and Felix Smith Wigner (1955); Smith (1960). Originally devised for nuclear scattering problems to deduce the time associated with a scattering event from stationary measurements of the asymptotic scattering amplitudes, this concept prominently resurfaced in mesoscopic physics Kottos (2005) and very recently in attosecond physics Pazourek et al. (2015) as well as in the newly emerging community of wavefront shaping Rotter and Gigan (2017); Fan and Kahn (2005); Rotter et al. (2011); Carpenter et al. (2015); Xiong et al. (2016); Gérardin et al. (2016).

The Wigner-Smith time-delay operator $Q$ is constructed based on a system’s scattering matrix $S$ by way of a frequency derivative, $Q=-iS^{-1}{\rm d}S/{\rm d}\omega$. The eigenvalues of $Q$, also called “proper delay times”, measure the time-delay associated with the scattering by a given potential Rotter and Gigan (2017); Fyodorov and Sommers (1997); Brouwer et al. (1997); Savin and Sommers (2003). The corresponding eigenvectors, also called “principal modes”, are states that can be associated with this well-defined time-delay – a property that makes them dispersion-free Fan and Kahn (2005) in the sense that a small variation of their input frequency does not change their spatial output profile. Moreover, principal modes have been shown to be “particle-like” in situations where ballistic scattering occurs Rotter et al. (2011). The potential for applications of these states as efficient communication channels in systems with multiple in- and output ports has recently taken center stage when first experiments reported on the successful implementation of principal modes in optical multi-mode fibers Carpenter et al. (2015); Xiong et al. (2016) as well as in resonant scattering media Gérardin et al. (2016); Böhm et al. (2016).

Here we demonstrate that the concept underlying the principal modes is not at all restricted to the time-delay operator $Q$ from above; specifically, we show that in the same way as the conventional principal modes are invariant with respect to a frequency variation, we may also create wave states that are invariant with respect to changes in the system configuration, like a local shift of a designated scatterer inside a disordered medium. While the frequency-insensitive principal modes are the eigenstates of the time-delay operator $Q=-iS^{-1}{\rm d}S/{\rm d}\omega$ (involving a frequency derivative), the class of states we introduce will be the eigenstates of a corresponding operator $Q_{\alpha}=-iS^{-1}{\rm d}S/{\rm d}\alpha$, where the parameter $\alpha$ stands, e.g., for the position of a movable scatterer. What is special about our new approach is that the eigenstates of $Q_{\alpha}$ are not only insensitive to a small change of $\alpha$, but that the associated eigenvalues indicate how strongly the corresponding conjugate quantity to $\alpha$ is affected by the scattering process. This insight allows us to maximally focus on or maximally avoid a specific target inside multiply scattering media just based on the medium’s scattering matrix $S$ and its variation due to a shift of the target particle that one aims to focus on.

We start by reviewing the principal modes’ remarkable property of being insensitive with respect to a change in frequency of the incident wave. The corresponding eigenvalue equation for the principal modes $\vec{u}_{n}$ (given as a coefficient vector in a certain basis) and for the proper delay times $\theta_{n}$ reads as follows,

[TABLE]

with $\omega$ being the frequency. In a waveguide system, as shown in Fig. 1, the incident principal modes can be decomposed into waves injected from the left and right lead respectively, i.e., $\vec{u}_{n}=(\vec{u}_{n,L},\,\vec{u}_{n,R})^{T}$. For a unitary scattering matrix, $S^{\dagger}S=\mathds{1}$, the time-delay operator $Q$ is Hermitian with real eigenvalues $\theta_{n}$. Using the input-output relation $\vec{v}_{n}(\omega)=S(\omega)\vec{u}_{n}$, we can rewrite Eq. (1) for a (static) time-delay eigenstate $\vec{u}_{n}$ evaluated at a chosen frequency $\omega_{0}$

[TABLE]

The expression on the right-hand side of Eq. (2) states that the output vector of a principal mode, ${\vec{v}_{n}=(\vec{v}_{n,L},\,\vec{v}_{n,R})^{T}}$, retains its original orientation when shifting the input frequency in the vicinity of $\omega_{0}$. This, in turn, translates into output patterns that are invariant (to first order) with respect to a frequency change (apart from a complex global factor). An aspect of the above derivation that has so far been unexploited is the fact that the derivative in Eq. (1) does not necessarily have to be taken with respect to the frequency. In other words: the symbol $\omega$ can in principle represent any other parametric dependence of the scattering matrix and the stability property of the corresponding principal modes will still hold with respect to a variation of this new parameter. Accordingly, we shift our attention to a whole class of generalized Wigner-Smith (GWS) operators $Q_{\alpha}=-iS^{-1}{\rm d}S/{\rm d}\alpha$ in which the frequency $\omega$ is replaced by the arbitrary parameter $\alpha$. We also expect that these GWS operators may have interesting connections to earlier works where statistical properties of parametric variations of the scattering matrix have been studied Lee and Stone (1985); Gorin et al. (2006); Fyodorov and Savin (2012).

While the eigenstates of $Q_{\alpha}$ are invariant with respect to a small parametric shift of $\alpha$ already by construction [as in Eq. (2) with $\omega\to\alpha$], we still need to clarify how to interpret the corresponding eigenvalues $\theta^{\alpha}_{n}$. Already from the dimensions it is clear that $\theta^{\alpha}_{n}$ must be associated with the conjugate variable to $\alpha$, in the same way as the delay time $\theta$ is the conjugate quantity to the frequency $\omega$. To make this more evident we now define $C_{\alpha}:=-i{\rm d}/{\rm d}\alpha$ as the corresponding conjugate operator and assume that the parameter $\alpha$ is global. It follows after a short derivation (see supplemental material Amb ),

[TABLE]

where $\vec{v}(\alpha)=S(\alpha)\vec{u}$ is the output vector and $\langle\ \ldots\ \rangle_{\mathrm{in/out}}$ denotes the expectation value evaluated in the input/output scattering state. Following Eq. (3), $Q_{\alpha}$ is the appropriate operator to measure the shift in the conjugate variable to $\alpha$, which the wave experiences due to the scattering process – in perfect analogy to the time-delay operator measuring a shift in time, i.e., in the conjugate quantity of the frequency $\omega$. One specific example for such a global parameter $\alpha$ could be the displacement in $y$-direction (i.e., $\alpha=y$) of the entire scattering landscape from its initial position. The conjugate variable to the position is the momentum. The eigenvalues $\theta^{y}_{n}$ of the operator $Q_{y}$ then measure the momentum shift (in $y$-direction) that the corresponding eigenstates experience as a result of scattering at the entire potential landscape. Solving the eigenvalue problem for the operator $Q_{y}$ thus provides us with an orthogonal and complete basis of eigenstates that are sorted by this momentum shift. To be more precise, the operator $Q_{y}$ generally measures the shift in the wavenumber in displacement direction, but for the sake of simplicity we refer to this shift as momentum shift $\Delta k$.

In the present context, we will be specifically interested in the case where $\alpha$ is the position $y_{i}$ of the $i$-th scattering element inside a strongly disordered medium consisting of altogether $i_{\rm tot}\geq i$ such individual elements (see, Fig. 1). In this case, $\alpha$ does not represent a global variable since only the $i$-th scatterer is shifted rather than the whole system. It turns out that Eq. (3) still remains valid for this more general case, but that the corresponding expectation values are then evaluated based on the local field amplitudes in the vicinity of the $i$-th scatterer (see derivation in Amb ). The eigenvalues $\theta^{y_{i}}_{n}$ of $Q_{y_{i}}$ are thus equal to the momentum shift $\Delta k^{i}$ that the wave experiences locally when colliding with the $i$-th scatterer (in the presence of all other scatterers and boundaries of the medium). Note that a large momentum shift requires that the wavefront of the incident eigenstate is back-scattered significantly at the target position $x_{i}$, which is equivalent to focusing on the $i$-th scatterer. On the other hand, when the wave does not get scattered at the target at all (e.g., by not reaching the $i$-th scatterer), the momentum transfer will, correspondingly, be very small. In this way the GWS operator $Q_{y_{i}}$ provides us the means to focus on a designated scatterer inside a disordered medium or to omit this target simply by generating an eigenvector of the GWS operator corresponding to a large or small eigenvalue, respectively.

To demonstrate the efficiency of this approach also in the experiment, we implemented a microwave scattering setup as displayed in Fig. 1, consisting of a rectangular multi-mode waveguide made of aluminum with dimensions $L$ $\times$ $W$ $\times$ $H$ = 2.38 m $\times$ 10 cm $\times$ 8 mm (see caption in Fig. 1). Additional absorbers (types: LS-14 and LS-16 from EMERSON&CUMING) at the ends of the waveguide mimic semi-infinite leads. Ten antennas are placed equidistantly in the incident part of the waveguide controlling ten transmission channels individually. The antennas are addressed by IQ-modulators (GTM 1 M2L-68A-5 of GT Microwave Inc.), which allow us to adjust phase and amplitude of the signal emitted by each antenna. The shaping of the incident wave is explained in more detail in Böhm and Kuhl (2016). This array of emitting antennas is connected via a power splitter (Microot MPD16-060180) to a 4 port vector network analyzer (VNA, Agilent E5071C). The transmission through the scattering system is measured with a single movable antenna placed at the output side of the waveguide. In the middle we place 18 cylindrical Teflon scatterers of radius $r=2.55~{}$mm and one brass scatterer of radius $r_{\rm b}=8.825~{}$mm forming the disordered scattering system. The relatively high refractive index of the Teflon elements ($n=1.44$) lets all waves undergo multiple scattering events before being transmitted. We operate the waveguide at 15.5 GHz, at which frequency ten propagating transverse electric modes TE0i ($i=1,\dots,10$) are supported. In the following, the parameter $\alpha$ corresponds to the position $y_{\rm b}$ of the central brass scatterer and the operator $Q_{y_{\rm b}}$ can be computed from the shift $\pm\Delta y_{\rm b}$ from its initial position, which we choose in the experiment to be of the same size as the in-plane radius of the brass scatterer.

As in most experiments Popoff et al. (2010a); Yu et al. (2013), we also only have access to a subpart of the entire scattering matrix. Specifically, we can measure only the $10\times 10$ transmission matrix $t$, where the complex matrix elements $t_{ji}$ stand for the transmission of the $i$-th antenna of the input antenna array to the $j$-th position of the movable antenna at the output. Even for flux-conserving scattering (without gain or loss), the $t$-matrix is generally non-unitary, since the reflected part of the incident wave is not contained in $t$. The derivation of Eq. (2) can, however, be easily adapted by replacing $S$ with $t$ and $\omega$ with $y_{b}$ in Eq. (1) (see also Juarez et al. (2012)). Most importantly, the resulting non-Hermitian GWS operator,

[TABLE]

inherits the property from its Hermitian counterpart $Q_{y_{\rm b}}$ that its eigenstates are invariant with respect to a small change in the parameter $y_{\rm b}$. In contrast to eigenstates of $Q_{y_{\rm b}}$, the eigenstates of $q_{y_{\rm b}}$ feature injection only from one lead, i.e., $\vec{u}_{n}\rightarrow\vec{u}_{n,L}$. The transmitted state, i.e., the outgoing state to the right, can be calculated via $\vec{v}_{R}=t\vec{u}_{L}$. Specifically, when adapting Eq. (2) to feature the complex eigenvalues $\vartheta^{y_{\rm b}}_{n}=\beta^{y_{\rm b}}_{n}+i\kappa^{y_{\rm b}}_{n}$ of $q_{x_{\rm b}}$, we have for the corresponding transmitted states

[TABLE]

Since the construction of the operator $q_{y_{\rm b}}$ involves only the transmission matrix (the reflected part is omitted), its complex eigenvalues $\vartheta^{y_{\rm b}}_{n}$ do no longer correspond directly to the local momentum shift $\Delta k^{\rm b}$ at the brass scatterer (see details in the supplemental material Amb ). We do find, however, that a strong correlation between these two quantities persists (see Fig. S3 in Amb ), basically since the trace left in the transmission matrix by moving the target scatterer is in itself already quite indicative of the focus on this target. Since this trace in the transmitted signal appears both in the phase and amplitude of a $q_{x_{b}}$-eigenstate as measured, respectively, by the real ($\beta^{y_{\rm b}}_{n}$) and imaginary ($\kappa^{y_{\rm b}}_{n}$) parts of $\vartheta^{y_{\rm b}}_{n}$, we work with the absolute value of $\vartheta^{y_{\rm b}}_{n}$ to quantify this trace. A more detailed analysis Amb shows that the correlation between $|\vartheta^{y_{\rm b}}_{n}|$ and $|\Delta k^{\rm b}|$ can be further increased by normalizing the latter term with the transmission (i.e., by working with $|\Delta k^{\rm b}|/|t^{2}|$). We can thus apply the concept we derived for the Hermitian GWS operator $Q_{y_{\rm b}}$ also to its non-Hermitian counterpart $q_{y_{\rm b}}$ with the essential difference being that the eigenvalues are now complex and sorted by their absolute value.

Following this protocol also in the experiment, we first measure the transmission amplitudes $t_{ji}$ for two slightly different positions of the brass scatterer ($\pm\Delta y_{\rm b}$). We then determine the GWS operator $q_{y_{\rm b}}$ by replacing the derivative in Eq. (4) with a finite-difference approximation based on the difference between the two prior transmission matrix measurements. In a next step we evaluate the eigenstates of $q_{y_{\rm b}}$ and inject them directly through the antenna array at the input port. To test if our focusing protocol works successfully, we then measure the intensity distribution of the microwave field in the vicinity of the brass scatterer by an additional scanning antenna. This antenna is attached to a movable arm and enters through a grid of holes (5 mm $\times$ 5 mm) in the top plate of the waveguide (extending $3$ mm into the cavity). The obtained intensity distributions for the eigenvectors with the three largest and the three smallest eigenvalues (in absolute magnitude) are shown in Figs. 2 and 3, respectively. The displayed intensity profiles demonstrate very clearly that the largest eigenvalues correspond to states focusing on the target, see top row in Fig. 2. The eigenstates with the smallest eigenvalues, in turn, produce intensity patterns which are reduced almost to noise level in the vicinity of the target, see top row in Fig. 3. As an additional test of our protocol, we also recorded the wave intensity patterns after removing the brass scatterer altogether. In case of the focusing states such an intervention drastically changes the overall configuration of the intensity profiles, with intensity maxima near the position of the removed scatterer remaining clearly visible, see bottom row in Fig. 2. In the case of the states that avoid the brass scatterer, the wave intensity pattern remains almost unchanged as a whole when the scatterer is removed, see bottom row in Fig. 3. For the sake of clarity we emphasize here that the measurement of the intensity distribution in the interior of the investigated system serves only for demonstration purposes and is not necessary for implementing our protocol in the first place.

To summarize, we present here an extension of the Wigner-Smith time-delay operator to a whole class of operators with the exciting property of providing eigenstates that focus on or avoid a designated target inside a disordered medium. These generalized Wigner-Smith (GWS) operators require the information stored in a system’s scattering matrix, but, as we show here, a construction based on the transmission (or reflection) matrix alone still maintains its useful features. Measuring the transmission matrix of a system has meanwhile not only been demonstrated for the presented case of microwaves, but also for acoustics Gérardin et al. (2014), seismology Larose et al. (2004), and recently also for optics Popoff et al. (2010a); Yu et al. (2013). As a “guidestar” Horstmeyer et al. (2015) for focusing deep in the multiple scattering regime, we use a movable scatterer inside the medium whose spatial shift leaves conspicuous traces in the measured transmission matrix that we exploit for our protocol. In the practical applications that we envision for future implementations, the spatial shift of the target scatterer could, e.g., result from the self-propelled movement of an object inside an otherwise static medium or be excited from the outside through an ultrasound focal spot that can be conveniently scanned through the medium (see Xu et al. (2011); Judkewitz et al. (2013) for recent implementations). Since, in particular, many biological media are turbid for optical light but only weakly scattering for ultrasound, such optical focusing techniques enabled by ultrasound have recently already been successfully explored in biomedical imaging (see Horstmeyer et al. (2015) for a review). So far, however, these techniques had to rely on “optical phase conjugation” techniques to time-reverse a scattered wave to a scattering target. Our new approach has the advantage that it works without such a time-reversal mirror for light and that the degree of focusing can be tuned up to the point where a target inside a medium can be entirely avoided rather than focused on. The latter feature may be particularly attractive for imaging techniques for which it is imperative that certain parts of an imaged tissue do not get exposed to radiation. While the presented experiment using ten guided modes serves as a proof-of-principle demonstration, we expect that the full potential of the method can be exploited once many modes are accessible as in the optical domain using spatial light modulators Mosk et al. (2012); Rotter and Gigan (2017).

I SUPPLEMENTAL MATERIAL

I.1 Definition of the $C_{\alpha}$-Operator

In the main text, we defined the GWS-operators $Q_{\alpha}$ using the derivative of the scattering matrix $S(\alpha)$ with respect to the (unspecified) parameter $\alpha$. According to our definition used in Eq. (3), $C_{\alpha}=-i{\rm d}/{\rm d}\alpha$ is the conjugate operator to the observable $\alpha$. To write this operator as a matrix we expand it in the modes that carry energy from the asymptotic region to the scattering region through its surface and vice versa. We write for the $i$-th mode in the asymptotic region $\psi_{n}(x,\vec{\xi})$, with $\vec{\xi}=(y,z)^{T}$ being the transversal coordinate (within the surface) and $x$ the longitudinal coordinate (perpendicular to the surface), respectively. The matrix elements of the operator associated to $C_{\alpha}$ can then be expressed as follows,

[TABLE]

where $\partial\Omega$ denotes the surface with dimension $D$ and $x_{s}$ is the longitudinal coordinate at which the integral is evaluated.

I.2 Wave Expectation Values of $Q_{\alpha}$

In order to obtain a general expression for the operators $Q_{\alpha}$ defined in the main text, we start with the definition of $C_{\alpha}$ in Eq. (S1) right above. Using the “translation” operator $\exp(iC_{\alpha}\Delta\alpha)$, where $\Delta\alpha$ is the shift in the parameter $\alpha$, we can write for the scattering matrix,

[TABLE]

For small deviations $\Delta\alpha$ we can thus approximate,

[TABLE]

where we neglected terms of the order $(\Delta\alpha)^{n}$ with $n\geq 2$. The derivation of $S(\alpha)$ as needed for the construction of $Q_{\alpha}$ can then be expressed as follows,

[TABLE]

With $S^{\dagger}S=\mathds{1}$, it follows for $Q_{\alpha}$,

[TABLE]

As discussed in the main text, the GWS-operator $Q_{\alpha}$ measures the shift in $\langle C_{\alpha}\rangle$ the wave experiences due to the scattering process. Note that, strictly speaking, we define here the variable $\alpha$ as a tunable property of the incident wave rather than that of the system. This means that, when, e.g., spatially moving the scattering system, as in the case when $\alpha=x$, the parameter $x$ refers to the position of the incident wave relative to the system. Moving the scattering region in the positive $x$-direction while keeping the incoming wavefronts unmoved thus corresponds to a variation of $\Delta x<0$ (in analogy to active/passive transformations). For this case, the expectation values of $Q_{x}$ read as follows

[TABLE]

where we used $C_{x}=k_{\rm in/out}$ in order to emphasize that the operator of the conjugate variable is now the momentum.

For the case where the system consists of $i_{\rm tot}$ scatterers at positions $x_{i}$ (also the boundaries of the scattering region are considered as individual scatterers), the GWS-operator $Q_{x}$ reads as follows

[TABLE]

where $Q_{x_{i}}$ stands for the partial GWS-operator corresponding to the positions of the individual scatterers. The sum of all these operators $Q_{x_{i}}$ is identical to the operator $Q_{x}$, and thus measures the total $\langle C_{x}\rangle$-shift as discussed above [see Eq. (S5)]. Hence, each of the summands $Q_{x_{i}}$ measures the contribution to the total shift in $\langle C_{x}\rangle$, that can be attributed to the respective scatterer at position $x_{i}$.

By analogy to Eq. (S6) we now make the conjecture (to be verified numerically below) that the expectation value of the partial GWS-operator $Q_{x_{i}}$ can be written as

[TABLE]

where $\vec{u}$ is the incoming wave from the asymptotic regions. The vector $\vec{u}^{\,i}$, on the other hand, represents the local incoming wave directly impinging on the $i$-th scatterer and $\vec{v}^{\,i}=S_{i}\vec{u}^{\,i}$ is the local wave emanating directly from the considered scatterer. $S_{i}$ is the scattering matrix describing only the scattering process at the $i$-th scatterer. $k_{\rm in}=-k_{\rm out}$ are the momentum operators in the basis of incoming and outgoing waves. Eq. (S8) thus means, that the operator $Q_{x_{i}}$ measures the momentum difference $\Delta k^{i}$ between the local incoming and outgoing waves directly at the shifted scatterer by means of the scattering matrix of the whole system. For $\vec{u}$ being an eigenstate of the partial GWS-operator, the expectation value in Eq. (S8) is just the corresponding eigenvalue, i.e., $\langle Q_{x_{i}}\rangle=\theta^{x_{i}}$. This operator would thus provide us with a unique tool to extract information about the local scattering process by making use of information that is available in the asymptotic regions.

In order to give an intuitive picture for the above derivation, we consider now a one-dimensional (1D) system in which the momentum operators can be expanded in right- and left-propagating plane waves as

[TABLE]

where $k$ is the scalar wavenumber. Inserting Eq. (S9) into Eq. (S8) and using $\vec{u}^{\,i}=(u^{\,i}_{L},u^{\,i}_{R})^{T}$ and $\vec{v}^{\,i}=(v^{\,i}_{L},v^{\,i}_{R})^{T}$, we get

[TABLE]

where the subscripts $L$ and $R$ denote the direction from which the waves enter or exit a scattering region. Fig. S1 shows a generic 1D scattering system consisting of Gaussian scatterers, one of which is the designated target that is shifted to take the derivative ${\rm d}/{\rm d}x_{i}$ in $Q_{x_{i}}$ (see black dashed double-arrow). The global incoming and outgoing waves $\vec{u}$ and $\vec{v}$ as well as the local incoming and outgoing waves directly before and after the scattering process with the shifted scatterer, $\vec{u}^{\,i}$ and $\vec{v}^{\,i}$, are labeled and indicated by blue arrows. Since for the eigenstates the expectation value $\langle Q_{x_{i}}\rangle$ is just the corresponding eigenvalue $\theta^{x_{i}}$, Eq. (S10) means that the eigenvalues of such a partial GWS-operator provide us with the information of the momentum shift the local wave acquires directly at the shifted scatterer (weighted with the corresponding intensities [see Eq. (S8) and Fig. S1].

In order to verify Eq. (S10) [which is the 1D version of Eq. (S8)], we now provide an explicit numerical verification based on random matrix theory. For this, we consider the system illustrated in Fig. S1 consisting of three scattering regions as described by three scattering matrices $S_{L}$, $S_{i}$ and $S_{R}$. $S_{L/R}$ is the scattering matrix of the system to the left/right of the shifted scatterer and $S_{i}$ is the scattering matrix of the $i$-th scatterer itself (see Fig. S1). For the numerical verification we test Eq. (S10) for many random configurations by assuming that $S_{L}$, $S_{i}$ and $S_{R}$ are symmetric and unitary random matrices drawn from the circular unitary ensemble (CUE). For the wavenumber $k$ and the position of the shifted scatterer we use random values. Note, that shifting a scatterer only changes the free propagation length before and after the actual scattering region – a quantity which we tune explicitly in our numerics. We calculate the eigenvalues $\theta^{x_{i}}$ of the GWS-operator for each configuration and compare them to the expression on the right-hand side of Eq. (S10), for which the local amplitudes $\vec{u}^{\,i}$ and $\vec{v}^{\,i}$ are accessible numerically. We find excellent agreement between these two calculated quantities, i.e., the maximal squared absolute error never exceeds $10^{-10}$ for all 1000 system configurations studied – a deviation which is extremely small compared to the value of $k$ that is of the order of $10^{1}$.

To show that Eq. (S8) is not only valid in 1D, but also in 2D (as in the experiment), we now consider a waveguide structure as shown in Fig. S2 in which several transverse lead modes can propagate. Also here we test the numerical agreement between the left-hand side of Eq. (S8) with its right-hand side for many random configurations, in analogy to the 1D procedure described above. The difference in 2D is, however, that all scattering matrices have dimension $2N\times 2N$, where $N$ is the number of propagating modes. In our simulations we take $S_{L}$, $S_{R}$ from th CUE, but the scattering matrix for the middle part, $S_{i}$, is obtained by an explicit numerical calculation for a single scatterer randomly placed in a 2D waveguide Rotter et al. (2000); Libisch et al. (2012). We calculate $Q_{x_{i}}$ through shifting this scatterer (in propagation direction) for each random configuration. In this way we obtain excellent agreement between the $Q_{x_{i}}$-eigenvalues [left-hand side of Eq. (S8)] and the local momentum shift [right-hand side of Eq. (S8)], which are independently calculated with our numerics. As in 1D, we find that the maximal squared absolute error stays below $10^{-10}$ for all system configurations, which is again small compared to the mode wavevectors that are of the order of $10^{1}$. This verifies our assumption that the GWS-operator connects asymptotic quantities stored in the scattering matrix of the system with local wave amplitudes at the shifted scatterer. For reasons of simplicity, we shifted the scatteres in the 1D and 2D verification in propagation direction. However, in the derivation of $Q_{\alpha}$ in Eq. (S5), no direction is distinguished from another one, such that Eq. (S8) also holds for shifting the scatterer in any other direction (such as the transverse direction considered in the experiment).

I.3 Wave Expectation Values of $q_{x_{i}}$

As in our experiment, where we only have access to a subpart of the full scattering matrix $S$, we use the non-Hermitian GWS-operator (4) based on the system’s transmission matrix $t$ in order to find states that focus on or omit a designated scatterer. Eigenstates of both, the Hermitian GWS-operator $Q_{x_{i}}$ as well as the non-Hermitian operator $q_{x_{i}}$ used in the experiment are to first order invariant with respect to a change of the position of the chosen scatterer (up to a global factor). Note that eigenstates of $Q_{x_{i}}$ generally require injection from both leads, whereas eigenstates of $q_{x_{i}}$ are injected only from one lead. The eigenvalues of $Q_{x_{i}}$ are real and correspond to the local momentum shift as shown above. The eigenvalues of $q_{x_{i}}$, however, are complex but can still be related to the local momentum shift as we will show in the following: Writing the Hermitian GWS-operator $Q_{x_{i}}$ in its block structure

[TABLE]

we consider now the upper left $Q_{x_{i}}^{11}$-block, which takes the explicit form $Q_{x_{i}}^{11}=-i(t^{\dagger}\frac{{\rm d}t}{{\rm d}x_{i}}+r^{\dagger}\frac{{\rm d}r}{{\rm d}x_{i}})$, with $t$ and $r$ being the transmission and reflection matrix, respectively. Calculating the expectation value of $Q_{x_{i}}$ for a wave entering only from the left lead, i.e., $\vec{u}=(\vec{u}_{L},0)^{T}$, we end up with $\langle Q_{x_{i}}\rangle=\vec{u}^{\dagger}Q_{x_{i}}\vec{u}=\vec{u}_{L}^{\dagger}Q_{x_{i}}^{11}\vec{u}_{L}=\langle Q_{x_{i}}^{11}\rangle=\langle k_{\rm in}\rangle^{i}-\langle k_{\rm out}\rangle^{i}=\Delta k^{i}$, where $\Delta k^{i}$ is the local momentum difference. Using the expression for $q_{x_{i}}=-it^{-1}\frac{{\rm d}t}{{\rm d}x_{i}}$, we can rewrite $Q_{x_{i}}^{11}=-i(t^{\dagger}\frac{{\rm d}t}{{\rm d}x_{i}}+r^{\dagger}\frac{{\rm d}r}{{\rm d}x_{i}})=t^{\dagger}tq_{x_{i}}-ir^{\dagger}\frac{{\rm d}r}{{\rm d}x_{i}}$. Calculating now the expectation value of $Q_{x_{i}}^{11}$ for an eigenstate of $q_{x_{i}}$ (in the following we will omit the eigenstate index for the sake of readability) and resolving for the complex eigenvalue $\vartheta^{x_{i}}$ of $q_{x_{i}}$ yields

[TABLE]

where $|t|^{2}=\langle t^{\dagger}t\rangle$ is the global transmittance of the eigenstate. As in the case of the Wigner-Smith time-delay operator, one can show that the eigenvalue can also be written in the following form Ambichl (2016)

[TABLE]

where the real part is the derivative of the global transmitted scattering phase $\varphi_{t}$ and the imaginary part can be interpreted as the change in the global output intensity with respect to a change of the local position of the $i$-th scatterer. It is worth noting that this relation still holds for a variation of a local parameter rather than a global one as in the case of time-delay. The reason for this is that in both the global and the local case the output vector is invariant to first order with respect to a change in a certain parameter (e.g., $\omega$ or $x_{i}$). Comparing Eq. (S12) with (S13) now shows that the imaginary part of $\vartheta^{x_{i}}$ can be related to (global) reflections, i.e.,

[TABLE]

which have been omitted in the construction of $q_{x_{i}}$. However, the reflection term in Eq. (S12) also features a real part which means that the real part of $\vartheta^{x_{i}}$ is not just given by $\Delta k^{i}/|t|^{2}$.

To assign further meaning to the real part of $\vartheta^{x_{i}}$, we write, in analogy to the time-delay operator, the expectation value of $Q_{x_{i}}$ as a sum of derivatives of the scattering phases of each mode weighted with its corresponding output intensity Ambichl (2016) and split the sum into two parts. One part covers the reflected waves and the other one covers the transmitted waves, respectively:

[TABLE]

where $N$ is the number of modes. For the evaluation of this expectation value with a $q_{x_{i}}$-eigenstate, we use the fact that the transmitted phase derivatives are the same for all modes, i.e., $\frac{{\rm d}\varphi_{t,n}}{{\rm d}x_{i}}=\frac{{\rm d}\varphi_{t}}{{\rm d}x_{i}}=\mathrm{Re}(\vartheta^{x_{i}})$. Thus, bringing this factor in front of the sum, using $\langle Q_{x_{i}}\rangle=\Delta k^{i}$ and resolving for the real part of $\vartheta^{x_{i}}$ gives

[TABLE]

Eq. (S16) shows that the real part of the $q_{x_{i}}$-eigenvalue, $\mathrm{Re}(\vartheta^{x_{i}})$, contains both the momentum shift $\Delta k^{i}$ as well as the global reflectivity of each mode multiplied by the corresponding phase derivatives. Thus, the desired correlation between a small/large $|\mathrm{Re}(\vartheta^{x_{i}})|$ (or $|\vartheta^{x_{i}}|$) and a small/large momentum shift $|\Delta k^{i}|$ is not as obvious as for the eigenvalues of $Q_{x_{i}}$. We therefore investigate with our numerical approach (see above), which quantity based on $\vartheta^{x_{i}}$ shows the strongest possible correlation with $\Delta k^{i}$. Testing several different expressions, we find the best results for $|\vartheta^{x_{i}}|$ and $|\Delta k^{i}|/|t|^{2}$ (see Fig. S3). Specifically, the probability for finding the two states (out of ten) with the highest/smallest values of $|\Delta k^{i}|/|t|^{2}$ among the states with the three highest/smallest $|\vartheta^{x_{i}}|$ is 99.9%/91.7% in our 2D random matrix model featuring 10 modes and 1000 random configurations. This observation allows us to relate the eigenvalues of $q_{x_{i}}$-eigenstates with its focusing strength even without knowledge of the reflection matrix. As we show in the next section, it is also possible to go beyond such a semi-empirical analysis by applying a special projection procedure to our $q_{x_{i}}$-operator, which restores a strong correlation between $|\mathrm{Re}(\vartheta^{x_{i}})|$ and $|\Delta k^{i}|$ itself.

I.4 Wave Expectation Values of $q_{x_{i}}$ for Singular Transmission Matrices $t$

Since the construction of the non-Hermitian operator $q_{x_{i}}$ involves the inverse of the transmission matrix $t$, we have to deal with the problem that $t$ can become singular. This occurs, e.g., when transmission channels are closed. As a result, a straightforward inversion of the transmission matrix fails, requiring instead the evaluation of an effective inverse by projecting the transmission matrix onto a subspace of highly transmitting channels. For the projection procedure we make use of a singular value decomposition (SVD) of the transmission matrix $t=U\Sigma V^{\dagger}$, where the matrices $U$ and $V$ contain column-wise the left and right singular vectors and the matrix $\Sigma=\mathrm{diag}(\{\sigma_{n}\})$ contains the singular values $\sigma_{n}$ on its diagonal. In a next step, we pick a certain subset of large singular values $\tilde{\Sigma}=\mathrm{diag}(\{\tilde{\sigma}_{n}\})$ and corresponding left and right singular vectors $\tilde{U}$ and $\tilde{V}$. With these matrices we can construct an effective inverse $t^{-1}=\tilde{V}(\tilde{U}^{\dagger}t\tilde{V})^{-1}\tilde{U}^{\dagger}=\tilde{V}\tilde{\Sigma}^{-1}\tilde{U}^{\dagger}$, where $\tilde{\Sigma}^{-1}=\mathrm{diag}(\{\tilde{\sigma}_{n}^{-1}\})$. Projecting also the derivative of $t$ onto this subspace, using the proper projection operators $P_{\tilde{U}}=\tilde{U}\tilde{U}^{\dagger}$ and $P_{\tilde{V}}=\tilde{V}\tilde{V}^{\dagger}$, we arrive at the following generalized construction rule for our operator

[TABLE]

which turns into our original expression for $q_{x_{i}}$ if all transmission channels, i.e., all singular values, are taken into account.

Apart from applying this procedure to non-invertible transmission matrices, it can also be used to restore a correlation between $|\mathrm{Re}(\tilde{\vartheta}^{x_{i}})|$ and $|\Delta k^{i}|$ since the projection on highly transmitting channels suppresses the reflection term in Eq. (S16) and increases the global transmittance of $\tilde{q}_{x_{i}}$-eigenstates. In order to numerically verify this, we again consider many random 2D configurations as above and look at the statistical correlation between these two quantities with and without cutting away reflecting channels. To quantify the strength of this correlation, we look at the probability for finding the state with the highest/smallest $|\Delta k^{i}|$ among the states with the two highest/smallest $|\mathrm{Re}(\tilde{\vartheta}^{x_{i}})|$. Without SVD these probabilities are 63.7%/79.3%, whereas the application of our procedure using the five highest transmitting states (out of ten) yields 95.5%/88.3%. This verifies that projecting onto highly transmitting channels provides the possibility to restore a correlation between $|\mathrm{Re}(\tilde{\vartheta}^{x_{i}})|$ and $|\Delta k^{i}|$.

I.5 Data Processing

For the processing of the data measured in the experiment we Fourier filter the transmission signal between the antenna array and the scanning antenna. The resulting intensity data are treated with a discrete two dimensional convolution of the form,

[TABLE]

with

[TABLE]

to generate Fig. 2 and 3. The parameter $s$ denotes the intensity of the measured transmission signal and $K$ is called the convolution kernel, which in our case takes only closely neighboring measurements (pixels) into account. Note that $x,~{}y,~{}x^{\prime}$ and $y^{\prime}$ are discrete quantities corresponding to the coordinates of the grid-holes of the scanning antenna.