Optimizing light storage in scattering media with the dwell-time operator

TL;DR

This paper introduces a dwell-time operator derived from the scattering matrix to optimize light energy storage in disordered media, revealing eigenstates that deposit more energy than open channels and enabling energy storage in complex targets.

Contribution

It develops a novel dwell-time operator from the scattering matrix and characterizes its eigenvalue distribution, advancing control of light storage in scattering media.

Findings

Eigenvalue distribution has finite support with a maximum dwell-time much larger than the typical value.

Highest dwell-time eigenstates deposit more energy than open channels.

Dwell-time operator can be used to store energy in resonant targets within complex media.

Abstract

We prove that optimal control of light energy storage in disordered media can be reached by wavefront shaping. For this purpose, we build an operator for dwell-times from the scattering matrix, and characterize its full eigenvalue distribution both numerically and analytically in the diffusive regime, where the thickness $L$ of the medium is much larger than the mean free path $\ell$. We show that the distribution has a finite support with a maximal dwell-time larger than the most likely value by a factor $(L/\ell)^2\gg 1 $. This reveals that the highest dwell-time eigenstates deposit more energy than the open channels of the medium. Finally, we show that the dwell-time operator can be used to store energy in resonant targets buried in complex media.

Equations

$$ϵ(\mathbf{r})|\psi (\mathbf{r}){|}^2=\frac{c^2}{2\omega }\nabla .\left({\partial }_{\omega }\psi \nabla {\psi }^{\ast }-{\psi }^{\ast }{\partial }_{\omega }\nabla \psi \right),$$

$$U=\frac{{ϵ}_0{c}^2}{4\omega }{\int }_{\mathcal{S}}\text{d}\mathbf{r}\mathbf{n}.\left({\partial }_{\omega }\psi \nabla {\psi }^{\ast }-{\psi }^{\ast }{\partial }_{\omega }\nabla \psi \right),$$

$$U={\varphi }^{\text{in}}⟨{\psi }^{\text{in}}|Q_d|{\psi }^{\text{in}}⟩,$$

$$Q_d=Q+Q_e+Q_i.$$

$$\mathcal{H}=2N{\tau }_s\sum _{n=1}^N{\gamma }_n-\sum _{n<m}\text{ln}|{\gamma }_n-{\gamma }_m|,$$

$$p(\tau )\simeq \frac{2{\tau }_s}{\pi {\tau }^2}\sqrt{\left(\alpha \frac{\tau }{{\tau }_s}-1\right)\left(1-\beta \frac{\tau }{{\tau }_s}\right)}\left(1+\gamma \frac{\tau }{{\tau }_s}\right),$$

$${\tau }^{\text{max}}\simeq \frac{4}{9}⟨\tau ⟩\left(\frac{⟨\tau ⟩}{{\tau }_s}-4\right).$$

$${\tau }^{\text{oc}}\simeq \frac{\pi }{12}⟨\tau ⟩\left(\frac{⟨\tau ⟩}{{\tau }_s}+\frac{12}{\pi }\right).$$

$${\tau }^{\text{max}}\simeq {\tau }_s\left[\zeta \left(\frac{⟨\tau {⟩}^2}{{\tau }_s^2}+4\frac{⟨\tau ⟩}{{\tau }_s}-4\right)-1\right]$$

Full text

Recent developments in wavefront shaping protocols have allowed spectacular demonstrations of light manipulation in complex media Mosk et al. (2012); Rotter and Gigan (2017), such as non-invasive imaging in biological tissues Horstmeyer et al. (2015); Yu et al. (2015), focusing Vellekoop and Mosk (2008) or enhanced power delivery Kim et al. (2012); Popoff et al. (2014); Hsu et al. (2017) behind opaque media, or focusing Choi et al. (2013); Ambichl et al. (2017); Katz et al. (2019) and enhanced absorption Liew et al. (2016) inside scattering materials. The large number of degrees of freedom supported by disordered systems has also been proposed as a ressource for imaging with high resolution Vellekoop et al. (2010); Park et al. (2014), controlling the strength of light-matter interaction Bachelard et al. (2014); Davy and Genack (2018); Pichler et al. (2019), or performing optically complex and reconfigurable operations Matthès et al. (2019).

In this context, great attention has been given to the statistical properties of the transmission matrix $t$ Popoff et al. (2010); Shi and Genack (2012); Goetschy and Stone (2013); Gérardin et al. (2014); Hsu et al. (2015). This matrix admits a significant fraction of singular states, called open channels, responsible for complete destructive (constructive) interference in reflection (transmission), even if the medium is opaque on average Dorokhov (1984); Nazarov (1994). The intensity map inside the medium resulting from the propagation of open channels has also been elucidated Davy et al. (2015a), revealing a bell-shape profile along the propagation direction very different from the familiar linear decay obtained with plane wave illumination van Rossum and Nieuwenhuizen (1999). This property makes open channels good candidates for enhancing significantly energy deep inside disordered media Sarma et al. (2016); Koirala et al. (2017); Hong et al. (2018). However, these states are, by construction, those maximizing the output flux (they are eigenstates of the operator $t^{\dagger}t$ associated with the largest eigenvalues), and not necessarily the stored energy. The transverse localization of open channels, discovered very recently Yilmaz et al. (2019), also supports the idea that they are not necessarily the ones optimizing energy storage.

In this Letter, we explicitly build an operator for the dwell time (or stored energy) in complex media illuminated with monochromatic light. Its expression can be obtained directly from the scattering matrix (including evanescent channels) of the disordered material and the dispersion properties of the surrounding medium. First, we show that the dwell-time (DT) operator is not strictly identical to the Wigner-Smith matrix, introduced historically to characterize the duration of a scattering process Smith (1960), by identifying a contribution resulting from the interference between the incident and scattered fields, similar to that predicted in 1D for electrons Winful (2003a). Second, we study its eigenvalue distribution $p(\tau)$ for wave propagating through a disordered slab of thickness $L\gg\ell$, where $\ell$ is the light mean free path. We find that for non-resonant scattering this distribution is parametrized by two time scales only: the scattering time $\tau_{s}\sim\ell/c$ ($c$ being the speed of light in vacuum), as well as the mean time $\langle\tau\rangle\sim L/c$, which is known to be remarkably independent of the disorder strength Pierrat et al. (2014); Savo et al. (2017). It also exhibits a dominant peak at $\tau\sim\tau_{s}$, and has a finite support with a maximal DT eigenvalue $\tau^{\textrm{max}}\sim\langle\tau\rangle^{2}/\tau_{s}$. This last result implies that the maximal energy that can be stored in a disordered medium by wavefront shaping with fixed input power $\phi^{\textrm{in}}$ scales as $U^{\textrm{max}}\sim\phi^{\textrm{in}}\tau_{\textrm{Th}}$, where $\tau_{\textrm{Th}}$ is the Thouless time. Finally, we also show that the DT operator is a powerful tool to selectively deposit energy on local resonant targets embedded in a given realization of a complex medium.

Let us start with the construction of the DT operator. For clarity, we restrict the present discussion to the propagation of scalar waves in non-resonant and non-absorbing materials, described by the equation $\left[\nabla^{2}+k^{2}\epsilon({\mathbf{r}})\right]\psi({\mathbf{r}})=0$. Here, $\epsilon({\mathbf{r}})$ is the (real) dielectric function, and $\psi$ is the complex amplitude of the monochromatic wave with frequency $\omega=ck$. The quantity to maximize is the electromagnetic energy $U=\epsilon_{0}\int_{\mathcal{V}}\textrm{d}{\mathbf{r}}\,\epsilon({\mathbf{r}})|\psi({\mathbf{r}})|^{2}/2$, where $\mathcal{V}$ is the volume occupied by the disordered slab. From the wave equation, we readily obtain the relation

[TABLE]

which allows us to express the energy as the surface integral

[TABLE]

where $\mathcal{S}$ denotes the input and output surfaces of the slab and ${\mathbf{n}}$ the outward normal on them. For a Schrödinger wave $\psi$, a similar relation holds for the probability $\int_{\mathcal{V}}\textrm{d}{\mathbf{r}}\,|\psi({\mathbf{r}})|^{2}$, which is independent of the explored potential Smith (1960), contrary to $U$ that depends on $\epsilon({\mathbf{r}})$. Next, we express the field $\psi$ on each surface in terms of the incident field $\psi^{\textrm{in}}$, and the reflection and transmission matrices. Care must be taken since the scattered field at the sample surface can have contributions from both propagating and evanescent channels of the surrounding medium. Here, we consider a disordered slab embedded in a multimode waveguide supporting $N$ propagating channels, so that the reflection and transmission matrices restricted to this channel basis, noted $r$ and $t$ respectively, are of size $N\times N$. Relegating the technical derivation to the Supplementary Material (SM) SI , we find that the stored energy can be expressed as the expectation value

[TABLE]

where the DT operator $Q_{d}$ is a sum of three contributions, with clear and distinct meanings discussed below:

[TABLE]

Depending on the physical situation of interest (choice of $\psi^{\textrm{in}}$, size and scattering strength of the medium), each of these terms can produce the dominant contribution to $U$. We discuss their expressions below for an incident field without evanescent component, and injected from one side of the slab only, which corresponds to the most common experimental situation. More general expressions are given in the SM SI .

The first term in the right-hand side in Eq. (4) is the well-known Wigner-Smith matrix $Q=-\textrm{i}(t^{\dagger}\partial_{\omega}t+r^{\dagger}\partial_{\omega}r)$, which characterizes the duration of the scattering process for quasi-monochromatic signals measured in the far field of the medium Smith (1960); Texier (2016); Rotter and Gigan (2017). The utility of this operator for controlling wave propagation in multimode fibers and scattering media has been demonstrated in recent years Rotter et al. (2011); Carpenter et al. (2015); Davy et al. (2015b); Gérardin et al. (2016); Xiong et al. (2016); Ambichl et al. (2017). The second term captures scattering contributions into evanescent channels, and contributes even if $\psi^{\textrm{in}}$ has no evanescent component. It reads $Q_{e}=(t_{e}^{\dagger}D^{e}t_{e}+r_{e}^{\dagger}D^{e}r_{e})/2$, where $r_{e}$ and $t_{e}$ are the reflection and transmission matrices into evanescent channels of the waveguide Mello and Kumar (2004). The matrix $D^{e}$ is diagonal, with elements $D^{e}_{\alpha\beta}=\frac{\partial_{\omega}\kappa_{\alpha}}{\kappa_{\alpha}}\delta_{\alpha\beta}$, where $\kappa_{\alpha}=\sqrt{q_{\alpha}^{2}-k^{2}}$ is the inverse decay length of the evanescent channel $\alpha$ ($q_{\alpha}=\alpha\pi/W$ in a 2D waveguide of width $W$). The contribution of $Q_{e}$ cannot be neglected close to the onset of a new propagating mode of the waveguide Bagwell (1990); Gómez-Medina et al. (2001); SI . The important impact of evanescent channels on dwell-times has also been revealed recently in the case of scattering by subwavelength particles Shen and Dogariu (2019). However, the contribution of $Q_{e}$ to the distribution $p(\tau)$ studied below turns out to be negligible for all frequencies except a discret set (see SM), and will not be discussed further.

The third term in the decomposition (4) describes a qualitatively different contribution, due to the interference between the incident and reflected propagating fields. Since the total field is $\psi\sim\psi^{\textrm{in}}+r\psi^{\textrm{in}}$ on the front surface and $\psi\sim t\psi^{\textrm{in}}$ on the back surface, the field products in Eq. (2) involve cross-terms in reflection only. The associated matrix reads $Q_{i}=-\textrm{i}(Dr-r^{\dagger}D)/2$, where $D$ has elements $D_{\alpha\beta}=\frac{\partial_{\omega}k_{\alpha}}{k_{\alpha}}\delta_{\alpha\beta}$, and $k_{\alpha}=\sqrt{k^{2}-q_{\alpha}^{2}}$. In the SM SI , we provide an alternative proof of Eq. (4) based on the continuity equation, that highlights the close connection between $Q_{i}$ and the interference term between the incident and scattered field at the origin of the optical theorem Akkermans and Montambeaux (2007). The contribution of $Q_{i}$ becomes appreciable for states with large reflection. For this reason, it cannot be neglected in strongly scattering media that reflect most of the light. In particular, it contributes substantially to the lower part of the DT eigenvalue distribution in the regime $L\gg\ell$ (see SM for details).

Equation (4) generalizes to arbitrary scattering media and arbitrary dimension (arbitrary $N$) the relation established for electrons in the case $N=1$ Winful (2003a), or for electromagnetic waves Winful (2003b), interacting with a simple barrier (for which $Q_{e}$ is zero). For electrons, the trace of the matrix $Q_{i}$ is known as a correction to the Friedel sum rule Texier (2016), that relates the density of states ($\sim\textrm{Tr}Q_{d}$) to the Wigner time delay ($\sim\textrm{Tr}Q$). In its operator form, the difference between $Q_{d}$ and $Q$ has also been pointed in Refs. Lagendijk and van Tiggelen (1996); Sokolov and Zelevinsky (1997); Rotter and Gigan (2017), but not expressed in the explicit and computationally useful expansion given by Eq. (4) amb .

To characterize the properties of the matrix $Q_{d}$, we performed extensive numerical simulations of scalar wave propagation through two-dimensional disordered slabs placed in a multimode waveguide, using the recursive Green’s function method Goetschy and Stone (2013), and computed $Q_{d}$ as defined in Eq. (4). The eigenvalue distribution $p(\tau)$ of $Q_{d}$ is represented in Fig. 1 for three values of the disorder strength $1/k\ell$, in the diffusive regime $kL\gg k\ell\gg 1$. We find that $p(\tau)$ has a pronounced peak that shifts towards small time as $k\ell$ decreases. This illustrates the fact that most of the light experiences a few scattering events before being reflected after a time $\sim\ell/c$. On the other hand, a close look at the largest eigenvalues (see inset) reveals that $p(\tau)$ is bounded, with an upper edge $\tau^{\textrm{max}}$ that increases with the disorder strength. This effect is triggered by light crossing the sample by diffusion, a process that takes more time when the mean free path is reduced since the number of scattering events is increased.

To support the previous observations, we developed an analytical model for the distribution $p(\tau)$. Since $Q_{d}$ and $Q$ have similar spectra, with differences observed at small times only (see SM), we work with $Q$ in the theoretical development. First, following Ref. Beenakker and Brouwer (2001), we use the simple relation that links the scattering matrix $S=\left(\begin{smallmatrix}r\\ t\end{smallmatrix}\right)$ in absence of absorption, to the scattering matrix $S^{a}(\omega)=S(\omega+\textrm{i}/2\tau_{a})$ with uniform absorption, where $\tau_{a}$ is the absorption time. A first-order expansion gives $Q(\omega)\simeq\tau_{a}[1-S^{a}(\omega)^{\dagger}S^{a}(\omega)]$ for $\omega\tau_{a}\gg 1$. The advantage of this relation lies in the fact that the joint probability distribution (JPD) of the eigenvalues of the operator $S^{a}(\omega)^{\dagger}S^{a}(\omega)$ is known, for disordered media excited from one side, in the limit $L\to\infty$ Bruce and Chalker (1996); Beenakker et al. (1996). Denoting by $\tau_{n}$ the eigenvalues of $Q$, the JPD of the decay rates $\gamma_{n}=1/\tau_{n}$ takes the form of the Gibbs distribution $p(\{\gamma_{n}\})\sim e^{-\mathcal{H}}$, with

[TABLE]

where $\tau_{s}$ is the scattering time. In 2D, for non-resonant scattering, it has to be defined as $\tau_{s}=(\pi/2)\ell/v$, where $v=c/n$ is the group velocity, $n$ being the effective refractive index of the disordered slab (see SM). The Laguerre distribution defined by Eq. (5) implies that the matrix $Q^{-1}$ behaves as a Wishart matrix in a disordered medium, a property which is also true in a chaotic cavity Brouwer et al. (1997); Grabsch and Texier (2016). In the first case, $Q^{-1}$ has the same JPD as the matrix $HH^{\dagger}/\tau_{s}$, where $H$ is a $N\times N$ complex Gaussian random matrix, while in the second case, $Q^{-1}$ has the JPD of $HH^{\dagger}/\left<\tau\right>$, where $H$ is of size $N\times N/2$.

The result in Eq. (5) was obtained for infinite-size systems. In this limit, the marginal distribution $p(\tau)$ depends on $\tau_{s}$ only, with infinite $\langle\tau\rangle=\left<\sum_{n=1}^{N}\tau_{n}\right>/N$. In non-resonant systems of finite size $L$, it is known that $\langle\tau\rangle$ scales as $\sim L/c$, both in the quasi-ballistic and diffusive regimes Pierrat et al. (2014); Savo et al. (2017); in 2D, $\langle\tau\rangle=(\pi/2)L/v$. To restore a finite mean time, we make the ansatz that $\mathcal{H}$ is still well approximated by expression (5) in the regime $L/\ell\gg 1$ and $N\gg 1$, and search for the marginal distribution $p(\tau)$ that minimizes $\mathcal{H}$ under the constraint $\int\textrm{d}\tau p(\tau)=\langle\tau\rangle$. We find (see the SM for details)

[TABLE]

where $\alpha$, $\beta$, and $\gamma$ obey a set of three coupled equations depending on the single parameter $b\equiv\langle\tau\rangle/\tau_{s}=L/\ell$. At large optical thickness ($b\gg 1$), they reduce to $\alpha\simeq 1$, $\beta\simeq 9/[4b(b-4)]$, and $\gamma\simeq 2\beta$. The distribution (6) is maximum for $\tau\simeq 4\tau_{s}/3$, and has a finite support $[\tau^{\textrm{min}},\tau^{\textrm{max}}]$, with $\tau^{\textrm{min}}\simeq\tau_{s}$ and

[TABLE]

These theoretical results are in excellent agreement with the numerical simulations, without adjustable parameter, as shown in Fig. 1. The distribution also exhibits a power $p(\tau)\sim\tau^{-3/2}$, for $\langle\tau\rangle\lesssim\tau\lesssim\tau^{\textrm{max}}$. This is a hallmark of diffusion, also observed in numerical simulation of the 2D kicked rotor dynamics Ossipov et al. (2003).

Our analysis reveals the existence of a finite maximal eigenvalue $\tau^{\textrm{max}}$ of $Q^{d}$. According to Eq. (3), the corresponding eigenstate, $|\psi^{\textrm{max}}\rangle$, should yield the largest amount $U^{\textrm{max}}$ of stored energy. To check this prediction, we compared the intensity pattern produced inside the slab by $|\psi^{\textrm{max}}\rangle$ with that resulting from the propagation of other remarkable wavefronts, such as the most open channel $|\psi^{\textrm{oc}}\rangle$ (the eigenstate of $t^{\dagger}t$ associated to the largest transmission eigenvalue $T\simeq 1$). Representative results are shown in Fig. 2(a). Contrary to the intensity profile created by the first mode of the waveguide (which behaves as a plane wave), both $|\psi^{\textrm{oc}}\rangle$ and $|\psi^{\textrm{max}}\rangle$ give rise to a concentration of energy deep inside the medium. In addition, the intensity pattern due to $|\psi^{\textrm{max}}\rangle$ is significantly larger than that due to $|\psi^{\textrm{oc}}\rangle$, when integrated over the transverse dimension $y$. This clearly shows that $U^{\textrm{max}}>U^{\textrm{oc}}$.

Using the analytical result (6), the ratio $U^{\textrm{max}}/U^{\textrm{oc}}$ can be evaluated precisely. In Ref. Davy et al. (2015a), it was shown that the intensity profile created by $|\psi^{\textrm{oc}}\rangle$ in 2D is $I^{\textrm{oc}}(x)=\int\textrm{d}y|\psi^{\textrm{oc}}(x,y)|^{2}=(\pi/2)[1+(\pi/2)(L/\ell)x(L-x)/L^{2}]/k$, where $x$ is the direction perpendicular to the slab. After integration over $x$, we obtain $U^{\textrm{oc}}=\phi^{\textrm{in}}\tau^{\textrm{oc}}$, with

[TABLE]

Hence, $U^{\textrm{oc}}$ grows quadratically with the sample thickness $L$. It is much larger than the energy $U^{\textrm{pw}}$ deposited by a plane wave, that grows linearly with $L$ [$I^{\textrm{pw}}(x)\sim(L-x)/(kL)$ gives $U^{\textrm{pw}}\sim\langle\tau\rangle\phi^{\textrm{in}}$], but it is smaller than $U^{\textrm{max}}$. Indeed, Eqs. (7) and (8) give $U^{\textrm{max}}/U^{\textrm{oc}}\simeq 16/3\pi>1$. These predictions are confirmed by the results of numerical simulations presented in Fig. 2(b), where $\tau^{\textrm{max}}=U^{\textrm{max}}/\phi^{\textrm{in}}$ and $\tau^{\textrm{oc}}=U^{\textrm{oc}}/\phi^{\textrm{in}}$ are shown to be both larger than the Thouless time, $\tau_{\textrm{Th}}=L^{2}/(\pi^{2}D_{B})$, which is the longest mode lifetime of the diffusion equation. In 2D, the light diffusion constant is $D_{B}=\ell v/2$, so that $\tau^{\textrm{max}}\simeq(\pi^{3}/9)\tau^{\textrm{Th}}$ and $\tau^{\textrm{oc}}\simeq(\pi^{4}/48)\tau^{\textrm{Th}}$ for $L\gg\ell$. In Fig. 2(b), we also show that the ratio $\tau^{\textrm{max}}/\tau^{\textrm{Th}}$ can be increased by injecting light from both sides of the sample. In this case, the average intensity profile $I^{\textrm{max}}(x)$ presents a mirror symmetry with respect to the middle of the slab $x=L/2$, as imposed by statistical invariance [see Fig. 2(a)]. Inspired by the microscopic approach developed in Ref. Ossipov (2018), we could establish (see SM for details) an expression for $\tau^{\textrm{max}}$ in this situation, which reads

[TABLE]

at large optical thickness, where $\zeta\simeq 0.57$ is the solution of a transcendental equation. This prediction also agrees with numerical simulations [see Fig. 2(b)].

We have discussed the properties of the DT operator $Q_{d}$ in statistically homogeneous non-resonant disordered materials, and demonstrated quantitatively superior (yet qualitatively similar) performances of the largest DT eigenstates for energy storage, compared to open channels. However, $Q_{d}$ is an operator specifically constructed to optimize the quality factor of an arbitrary complex structure. This concept is radically different from the monochromatic scattering properties usually captured by matrices such as $t^{\dagger}t$ or $r^{\dagger}r$. To illustrate this last point, let us consider a set of small absorbers buried in an otherwise non-absorbing disordered medium. One of them, resonant at frequency $\omega_{0}$ with quality factor $\omega_{0}/\Gamma\gg 1$, is described by the Lorentzian dielectric function $\epsilon=1+\omega_{0}^{2}\mathcal{F}/(\omega_{0}^{2}-\omega^{2}-i\omega\Gamma)$. Our goal is to compare the performances of the absorption matrix $A=1-t^{\dagger}t-r^{\dagger}r$ and the DT matrix. In the presence of absorption, we can show that the DT matrix is related to the Hermitian part of $Q_{d}$ defined as $Q^{H}_{d}=(Q_{d}+Q_{d}^{\dagger})/2$ (see SM). We show in Fig. 3 the intensity patterns calculated inside the medium at the resonance frequency $\omega_{0}$, and resulting from the propagation of the states $|\psi^{a}\rangle$, associated to the largest absorption eigenvalue, and $|\psi^{H}\rangle$, associated to the extremal eigenvalue of $Q^{H}_{d}$. We clearly see that $|\psi^{a}\rangle$ deposits energy indistinctly on all absorbers, whereas $|\psi^{H}\rangle$ focuses specifically on the resonant scatterer that will induce the largest dwell time.

In summary, we have presented a general setting for tuning light storage in arbitrary complex media based on the dwell-time operator (4). We showed that the eigenvalue distribution of this operator takes the universal form (6) in the diffusive regime, and demonstrated the possibility to reach more than one hundred percent energy increase with respect to what can be achieved with open channels. Finally, we established that $Q_{d}$ can be used for addressing hidden resonant targets without need for guide stars.