Multiple scattering theory in one dimensional space and time dependent disorder: average field [Invited]

Alexandre Selvestrel; Julia Rocha; Rémi Carminati; Rémi Carminati; Romain Pierrat

doi:10.1364/OME.517488

1. Introduction

Light (or more generally wave) propagation in spatially disordered media has been an active topic for many decades, stimulated by basic questions in fundamental physics and by a large number of applications. On the fundamental side, the existence of Anderson localization for different kinds of waves is an emblematic example, among many other questions in mesoscopic physics [1]. On the applied side, imaging and sensing [2] or light control in complex materials [3] are highly developed research themes. The basic concepts and theoretical tools for modeling light propagation in spatially disordered media are known to a large extent [4].

Beyond spatial modulation of the medium, there has recently been a surge in research on propagation of different kinds of waves in time-dependent media, including electromagnetic [5], optical [6–9], acoustic [10] or water waves [11,12]. This emerging field opens new perspectives in terms of applications. For example, periodic space-time metamaterials offer new degrees of freedom for wave control [13–15]. It also stimulates the development of appropriate theories, in an area that has been largely unexplored so far. For example, some of us have highlighted the atypical behavior of wave propagation in a time-varying disorder, showing that the average energy of the field grows exponentiall at long times [16], providing a theoretical support to observations based on numerical simulations [7] or experiments [17]. Another recent study has focused on the role of correlations in time disorder in providing innovative optical properties [18]. These bricks contribute to the development of theories of wave propagation in time-varying disordered media, which remains a widely open topic.

In this article, we address the question of light propagation in a medium exhibiting both space and time disorders. To start with a simple model, we consider a one-dimensional space disorder combined to a time modulation, resulting in a medium described by a fluctuating dielectric function $\epsilon (x,t)$ considered to be a random variable, with $x$ and $t$ the space and time coordinates, respectively. We assume that the medium exhibits short-range correlations in both space and time, without cross correlation between them. The main objective is to develop a theory for the average field (or intensity) proving the existence of a scattering mean-free path $\ell _s$ and a scattering mean-free time $\tau _s$, and to provide explicit expressions in the weak scattering regime. The paper is organized as follows: In Sec. 2, we develop the theory that extends the standard multiple scattering theory to a situation with both space and time disorders. We provide expressions for $\ell _s$ and $\tau _s$ using a perturbative approach. In Sec. 3, we consider the particular case of disorder with a gaussian correlation in space and time, and show that the expressions of the mean-free path and mean-free time are in full agreement with numerical simulations performed without approximations.

2. Multiple scattering theory for space-time disorder

In this section we build a theory to compute the average electric field, from which we will define $\ell _s$ and $\tau _s$, and derive their explicit expressions. To proceed, we generalize the standard multiple scattering theory to account for space-time disorder. The interested reader can find detailed presentations of multiple scattering theory in various textbooks [1,4,19,20]. In a medium with one-dimensional space-time disorder described by a random dielectric function $\epsilon (x,t)$, an electric field linearly polarized along the $y$-direction obeys the equation

(1)$$-\frac{\partial^2 E(x,t)}{\partial x^2}+\frac{1}{c^2}\frac{\partial^2}{\partial t^2}\left[\epsilon(x,t)E(x,t)\right]=S(x,t) \, ,$$

which is easily derived form Maxwell’s equations. Here $E(x,t)$ is the real amplitude of the field in the time domain, $c$ is the speed of light in vacuum, and $S(x,t)$ is a source term that we do not need to specify. It is interesting to note that in Eq. (1) the dielectric function remains within the time derivative operator, which has important consequences as will be seen later. This is a feature of scattering problems involving two types of disorder, for example with both permittivity and permeability disorders [21], or in acoustics with mass density and compressibility disorders [22]. We also note that working with the displacement field $D$, as in Ref. [16], does not simplify the equation when space and time disorders coexist.

2.1 Lippmann-Schwinger equation

The first step of the derivation consists in defining a homogeneous reference (or background) medium with permittivity $\epsilon _b$. The reference field $E_b$ in this medium satisfies

(2)$$-\frac{\partial^2 E_b(x,t)}{\partial x^2}+\frac{\epsilon_b}{c^2}\frac{\partial^2 E_b(x,t)}{\partial t^2}=S(x,t).$$

The choice of $\epsilon _b$ will be specified later, with the constraint that it should be close to the typical value of $\epsilon (x,t)$ to ensure the accuracy of the perturbative approach.

Subtracting Eq. (2) from Eq. (1) leads to

(3)$$-\frac{\partial^2 E_s(x,t)}{\partial x^2}+\frac{\epsilon_b}{c^2}\frac{\partial^2 E_s(x,t)}{\partial t^2} =\frac{\epsilon_b}{c^2}\frac{\partial^2 E(x,t)}{\partial t^2}-\frac{1}{c^2}\frac{\partial^2}{\partial t^2}\left[\epsilon(x,t)E(x,t)\right]$$

where $E_s=E-E_b$ is the scattered field. Equation (3) shows that the scattered field can be seen as a field propagating in the reference medium and caused by a complex source term given by the right-hand side. We now introduce the Green function $G_b$ defined as the solution to

(4)$$-\frac{\partial^2 G_b(x-x',t-t')}{\partial x^2}+\frac{\epsilon_b}{c^2}\frac{\partial^2 G_b(x-x',t-t')}{\partial t^2} =\delta(x-x')\delta(t-t')$$

where $\delta$ is the Dirac delta function, satisfying Sommerfeld’s radiation condition in space (or equivalently causality in time). The Green function can be understood as the electric field radiated in the reference medium by a point source emitting an infinitely short pulse. Its detailed calculation is given in App. A. In the Fourier domain, the Green function is

(5)$$G_b(k,\omega)=\operatorname{PV}\left[\frac{1}{k^2-\omega^2/v^2}\right] +\frac{i\pi}{2k}\delta\left(k-\frac{\omega}{v}\right) -\frac{i\pi}{2k}\delta\left(k+\frac{\omega}{v}\right)$$

where $v=c/\sqrt {\epsilon _b}$ is the phase velocity in the reference medium and $\operatorname {PV}$ stands for the Cauchy principal value operator. In the following, we will consider two different problems: (1) the evolution of the wave in space for a monochromatic incident beam, and (2) the evolution of the wave in time for an incident pulse with a fixed wave number. Expressions for the reference Green function in the $(x,\omega )$ and in the $(k,t)$ domains are thus required for problems (1) and (2), respectively. They are given by

(6)$$G_b(x,\omega)=\frac{i}{2k_b}\exp\left(ik_b|x|\right) \quad;\quad G_b(k,t)=\frac{\operatorname{H}(t)v^2}{\omega_b}\sin\left(\omega_b t\right)$$

where $\operatorname {H}$ is the Heaviside step function, $k_b=\omega /v$ and $\omega _b=kv$. It is important to note that the observed asymmetry between space and time arises from the different boundary conditions in both cases. In space, we have considered the Sommerfeld’s radiation condition which states that the field should propagate towards $+\infty$ (respectively $-\infty$) for $x>0$ (respectively $x<0$). In time, this condition is equivalent to causality which translates into the presence of the $\operatorname {H}$ operator.

Equation (3) together with the definition of the Green function $G_b$ allows us to write the scattered field $E_s$ in the integral form

(7)$$E_s(x,t)={-}\iint G_b(x-x',t-t')\frac{\partial^2\left[\epsilon(x',t')-\epsilon_b\right]E(x',t')}{c^2\partial t'^2}\mathrm{d} x'\mathrm{d} t'.$$

The total field $E=E_b + E_s$ obeys the integral equation

(8)$$E(x,t)=E_b(x,t)-\iint G_b(x-x',t-t')\frac{\partial^2\left[\epsilon(x',t')-\epsilon_b\right]E(x',t')}{c^2\partial t'^2}\mathrm{d} x'\mathrm{d} t' \, ,$$

known as the Lippmann-Schwinger equation. This equation is the elementary building block of multiple scattering theory. For further developments, it will prove useful to manipulate formal operator expressions. To that end, we define

(9)$$\mathcal{V}={-}\frac{\partial^2\left[\epsilon(x',t')-\epsilon_b\right]\bullet}{c^2\partial t^2} \quad;\quad \mathcal{G}_b=\iint G_b(x-x',t-t')\bullet\mathrm{d} x'\mathrm{d} t',$$

where the bullets are to be understood as the quantity on which the operator acts. In this formalism, the Lippmann-Schwinger Eq. (8) can be rewritten as

(10)$$E=E_b+\mathcal{G}_b\mathcal{V} E.$$

We emphasize that the main difference with the usual Lippmann-Schwinger equation appearing in standard multiple scattering theory is the operator character of the scattering potential $\mathcal {V}$.

2.2 Born series and Dyson equation

In order to estimate the field averaged over an ensemble of realizations of disorder (i.e., of the random variable $\epsilon (x,t,)$), we first expand Eq. (10) in the form

(11)$$E=E_b+\mathcal{G}_b\mathcal{V} E_b+\mathcal{G}_b\mathcal{V}\mathcal{G}_b\mathcal{V} E_b+\mathcal{G}_b\mathcal{V}\mathcal{G}_b\mathcal{V}\mathcal{G}_b\mathcal{V} E_b+\ldots$$

which is known as the Born series. Performing a statistical ensemble average, we find that

(12)$$\left\langle E\right\rangle=E_b+\mathcal{G}_b\left\langle \mathcal{V}\right\rangle E_b+\mathcal{G}_b\left\langle \mathcal{V}\mathcal{G}_b\mathcal{V}\right\rangle E_b+\mathcal{G}_b\left\langle \mathcal{V}\mathcal{G}_b\mathcal{V}\mathcal{G}_b\mathcal{V}\right\rangle E_b+\ldots$$

where $\left \langle \ldots \right \rangle$ denotes the average value. The problem now reduces to the computation of terms of the form $\left \langle \mathcal {V}(\mathcal {G}_b\mathcal {V})^n\right \rangle$. Let us first consider the second order term (i.e., $n=1$). We define the connected part of the correlation function of the potentiel by

(13)$$\left\langle \mathcal{V}\mathcal{G}_b\mathcal{V}\right\rangle=\left\langle \mathcal{V}\right\rangle\mathcal{G}_b\left\langle\mathcal{V}\right\rangle+\left\langle \mathcal{V}\mathcal{G}_b\mathcal{V}\right\rangle_c.$$

This corresponds to a splitting of the correlation function into a factorizable and a non-factorizable (connected) part. Similar splittings for more complicated terms would require relatively heavy writing. A convenient way to manipulate such expressions is to use diagrams. For Eq. (13), we write

(14)

where circles, solid lines and dashed lines represent scattering events (interactions with the scattering potential), Green functions of the reference medium, and connections (non-factorizable part of the correlation function), respectively. Using diagrams, the third-order case ($n=2$) becomes

(15)

and similarly for higher-order terms.

The key idea to obtain an equation for the average field consists in defining a new operator $\mathcal {S}$ containing all non-factorizable terms, i.e.

(16)

With this definition, Eq. (12) can be factorized in the form

(17)$$\left\langle E\right\rangle=E_b+\mathcal{G}_b\mathcal{S}\left\langle E\right\rangle$$

which is known as the Dyson equation. Equation (17) is exact and all the complexity of the multiple scattering problem lies in the closed form of the equation and in the operator $\mathcal {S}$. In order to define the scattering mean-free path and time, and to derive explicit expressions, we need to simplify this operator. To that end, let us consider the first term corresponding to a single scattering event. Applying the operator to the average field leads to

(18)$$\mathcal{S}^{(1)}\left\langle E\right\rangle ={-}\frac{1}{c^2}\frac{\partial^2}{\partial t^2}\left[\left\langle\epsilon(x,t)-\epsilon_b\right\rangle\left\langle E(x,t)\right\rangle\right].$$

We now need to make a choice for the reference medium. Taking $\epsilon _b=\left \langle \epsilon (x,t)\right \rangle$ ensures the accuracy of the pertubation method that we will use, by implying a vanishing first order in the perturbative expansion. Indeed, by defining the fluctuating part of the permittivity by $\delta \epsilon (x,t)=\epsilon (x,t)-\left \langle \epsilon (x,t)\right \rangle =\epsilon (x,t)-\epsilon _b$ we find that

(19)$$\mathcal{S}^{(1)}\left\langle E\right\rangle ={-}\frac{1}{c^2}\frac{\partial^2}{\partial t^2}\left[\left\langle\delta\epsilon(x,t)\right\rangle\left\langle E(x,t)\right\rangle\right]=0.$$

For the second order term, we obtain

(20)$$\mathcal{S}^{(2)}\left\langle E\right\rangle =\frac{1}{c^4}\iint \frac{\partial^2}{\partial t^2}\left\{\left\langle\delta\epsilon(x,t)G_b(x-x',t-t') \frac{\partial^2}{\partial t'^2}\left[\delta\epsilon(x',t')\left\langle E(x',t')\right\rangle\right]\right\rangle\right\}\mathrm{d} x'\mathrm{d} t'$$

where we have used the relationship $\left \langle \delta \epsilon (x,t)\delta \epsilon (x',t')\right \rangle _c=\left \langle \delta \epsilon (x,t)\delta \epsilon (x',t')\right \rangle$, generally referred to as the correlation function of disorder. The first second-order time derivative relates to the variable $t$ and thus applies only to $\delta \epsilon (x,t)G_b(x-x',t-t')$. The second derivation relates to $t'$ and applies to $\delta \epsilon (x',t')\bullet$. It will prove useful to perform a double integration by parts. Assuming that the correlation function of disorder vanishes at long times, Eq. (20) reduces to

(21)$$\mathcal{S}^{(2)}\left\langle E\right\rangle =\frac{1}{c^4}\iint \left\langle\frac{\partial^2}{\partial t^2}\left[\delta\epsilon(x,t) \frac{\partial^2}{\partial t'^2}G_b(x-x',t-t')\right]\delta\epsilon(x',t')\right\rangle\left\langle E(x',t')\right\rangle\mathrm{d} x'\mathrm{d} t'.$$

Similar transformations can be performed on the higher order terms in $\mathcal {S}$, but are not written here since they will not be useful in practice. Finally, $\mathcal {S}$ can be written as

(22)$$\mathcal{S}=\iint \Sigma(x,x',t,t')\bullet\mathrm{d} x'\mathrm{d} t'$$

where $\Sigma$ is the self-energy and is here a simple multiplicative function (not an operator). Using the self-energy, the average field can be written

(23)$$\left\langle E(x,t)\right\rangle=E_b(x,t)+\idotsint G_b(x-x',t-t')\Sigma(x',x^{\prime\prime},t',t^{\prime\prime})\left\langle E(x^{\prime\prime},t^{\prime\prime})\right\rangle\mathrm{d} x'\mathrm{d} x^{\prime\prime}\mathrm{d} t'\mathrm{d} t^{\prime\prime} \, ,$$

which is the integral form of the Dyson equation.

2.3 Weak-scattering regime

To derive expressions for the scattering mean-free path and time, we now consider the particular case of a source term of the form $S(x,t)=\delta (x)\delta (t)$ in an infinite medium. In this case, $E_b(x,t)=G_b(x,t)$ and $\left \langle E(x,t)\right \rangle =\left \langle G(x,t)\right \rangle$, with $G$ the Green function of the medium in the presence of disorder. We assume statistical homogeneity in space and time, such that $\Sigma (x',x'',t',t'')$ only depends on $x'-x''$ and $t'-t''$. In these conditions, Eq. (23) reduces to

(24)$$\left\langle G(x,t)\right\rangle=G_b(x,t)+\idotsint G_b(x-x',t-t')\Sigma(x'-x^{\prime\prime},t'-t^{\prime\prime})\left\langle G(x^{\prime\prime},t^{\prime\prime})\right\rangle\mathrm{d} x'\mathrm{d} x^{\prime\prime}\mathrm{d} t'\mathrm{d} t^{\prime\prime}.$$

This equation can be solved by performing a space-time Fourier transform, which leads to

(25)$$\left\langle G(k,\omega)\right\rangle=G_b(k,\omega)+G_b(k,\omega)\Sigma(k,\omega)\left\langle G(k,\omega)\right\rangle \, ,$$

or equivalently

(26)$$\frac{1}{\left\langle G(k,\omega)\right\rangle}=\frac{1}{G_b(k,\omega)}-\Sigma(k,\omega) =k^2-\frac{\omega^2}{v^2}-\Sigma(k,\omega).$$

This expression of the average Green function will be used to derive expressions for the scattering mean-free path $\ell _s$ and scattering mean-free time $\tau _s$. We start by considering a monochromatic source term, oscillating at a frequency $\omega$, and we focus on the spatial behavior of the average field given by

(27)$$\left\langle G(x,\omega)\right\rangle=\int_{-\infty}^{+\infty}\frac{e^{ikx}}{k^2-k_b^2-\Sigma(k,\omega)}\frac{\mathrm{d} k}{2\pi}.$$

The computation of this inverse Fourier transform requires additional hypotheses. Considering the weak-scattering regime defined by the condition $|\Sigma (k,\omega )|\ll k_b^2$, the self-energy has a significant contribution only when $k \simeq \pm k_b$. Assuming that the disorder is statistically isotropic, we also have $\Sigma (k,\omega )=\Sigma (-k,\omega )$. As a result, the self-energy can be taken on-shell for $k=k_b$ in Eq. (27). Under this assumption, Eq. (27) becomes

(28)$$\left\langle G(x,\omega)\right\rangle\approx\int_{-\infty}^{+\infty}\frac{e^{ikx}}{k^2-k_b^2-\Sigma(k_b,\omega)}\frac{\mathrm{d} k}{2\pi}.$$

In order to compute the integral, we apply the residue theorem. For $x>0$, we use the contour plotted in Fig. 1 (a). The semicircle in the upper plane is chosen in order to apply Jordan’s lemma. The poles are $k^{\pm }=\pm \sqrt {k_b^2+\Sigma (k_b,\omega )}$. Assuming that $\operatorname {Im}\Sigma (k_b,\omega )>0$, which will be justified below, we obtain

(29)$$\left\langle G(x>0,\omega)\right\rangle=\frac{ie^{ik^+x}}{k^+{-}k^-}=\frac{ie^{ik^+x}}{2k^+}.$$

For $x<0$, we use the contour presented in Fig. 1 (b) and we obtain

(30)$$\left\langle G(x<0,\omega)\right\rangle=\frac{ie^{{-}ik^+x}}{k^+{-}k^-}=\frac{ie^{{-}ik^+x}}{2k^+}$$

which finally leads to

(31)$$\left\langle G(x,\omega)\right\rangle=\frac{ie^{ik_e|x|}}{2k_e} \, ,$$

where $k_e=k^+$. We clearly see from Eqs. (6) and (31) that the average field propagates in a homogeneous effective medium with an effective wavevector $k_e$. It is convenient to split $k_e$ into its real and imaginary parts. We write $k_e=k_r+i/(2\ell _s)$, with $\ell _s$ the scattering mean-free path, and $k_r=n_r\omega /c$, which defines the real part $n_r$ of the effective refractive index of the medium. The intensity of the average field is then given by

(32)$$\left|\left\langle G(x,\omega)\right\rangle\right|^2=\frac{e^{-|x|/\ell_s(\omega)}}{4\left|k_e\right|^2}$$

which will be the expression used later for comparison with numerical simulations. In the weak scattering regime, we have

(33)$$\ell_s(\omega)=\frac{k_b}{\operatorname{Im}\Sigma(k_b,\omega)}.$$

We also note that the approximation $k_r\simeq k_b$ holds in the weak-scattering regime. A more refined expression would involve the real part of the self-energy.

Fig. 1. Integration contours used to compute the average Green function. (a) and (b) are used for the inverse Fourier transform for positive and negative positions, respectively. (c) and (d) are used for the inverse Fourier transform for negative and positive times, respectively. In these representations, we have assumed $\operatorname {Im}\Sigma (k_b,\omega )>0$ for (a) and (b) and $\operatorname {Im}\Sigma (k,\omega _b)>0$ for (c) and (d), as explained in the main text.

Download Full Size | PDF

We now turn to the illumination by a pulse source term with a fixed $k$-vector, and we focus on the temporal evolution of the average Green function, which is given by

(34)$$\left\langle G(k,t)\right\rangle=\int_{-\infty}^{+\infty}\frac{e^{{-}i\omega t}}{\omega_b^2/v^2-\omega^2/v^2-\Sigma(k,\omega)}\frac{\mathrm{d}\omega}{2\pi}.$$

The weak-scattering regime amounts to assuming that $|\Sigma (k,\omega )|\ll \omega _b^2/v^2$. Under this assumption, the self-energy takes significant values for $\omega \simeq \pm \omega _b$. For statistically isotropic disorder, such that $\Sigma (k,\omega )=\Sigma (-k,\omega )$, and making use of the fact that $\Sigma (x,t)$ is real valued, we find that

(35)$$\begin{aligned}\Sigma(k,-\omega)=\int_{-\infty}^{+\infty}\Sigma(x,t)e^{{-}ikx-i\omega t}\mathrm{d} x\mathrm{d} t &=\left[\int_{-\infty}^{+\infty}\Sigma(x,t)e^{ikx+i\omega t}\mathrm{d} x\mathrm{d} t\right]^*\\ &=\left[\int_{-\infty}^{+\infty}\Sigma(x,t)e^{{-}ikx+i\omega t}\mathrm{d} x\mathrm{d} t\right]^* =\Sigma(k,\omega)^*. \end{aligned}$$

As a result, the self-energy can be replaced by $\Sigma (k,\omega _b)^*$ in Eq. (34) in the vicinity of $-\omega _b$, and by $\Sigma (k,\omega _b)$ in the vicinity of $\omega _b$. This is the counterpart of the on-shell approximation in the frequency domain. In order to compute the integral, we now make use of the residue theorem. For $t<0$, we use the contour in Fig. 1 (c). The poles are $\omega ^-=-v\sqrt {k^2-\Sigma ^*(k,\omega _b)}$ and $\omega ^+=v\sqrt {k^2-\Sigma (k,\omega _b)}$. If $\operatorname {Im}\Sigma (k,\omega _b)>0$, we get $\left \langle G(k,t<0)\right \rangle =0$. This is the signature of causality in the time-domain Green function. For $t>0$, considering the contour in Fig. 1 (d), we find that

(36)$$\left\langle G(k,t>0)\right\rangle=iv^2\left[ \frac{e^{{-}i\omega^+t}}{\omega^+{-}\omega^-} +\frac{e^{{-}i\omega^-t}}{\omega^-{-}\omega^+} \right] .$$

Defining $\omega _e=\omega _r-i/(2\tau _s)$, with $\tau _s$ the scattering mean-free time, we finally obtain

(37)$$\left\langle G(k,t)\right\rangle=\frac{\operatorname{H}(t)v^2}{\omega_r}\sin(\omega_r t)e^{{-}t/[2\tau_s(k)]}$$

for the average field. To wash out the rapid oscillations in the intensity, we square this expression and perform a time average over a window with width $T$ such that $2\pi /\omega _r \ll T\ll \tau _s$. This leads to

(38)$$\overline{\left\langle G(k,t)\right\rangle^2}=\frac{\operatorname{H}(t)v^4}{2\omega_r^2}e^{{-}t/\tau_s(k)} \, ,$$

which will be the expression used for comparison with numerical simulations. Under the weak-scattering approximation, the scattering mean-free time is

(39)$$\tau_s(k)=\frac{k^2}{\omega_b\operatorname{Im}\Sigma(k,\omega_b)}.$$

We note that $\omega _r \simeq \omega _b$ in the weak-scattering regime. We also stress that having $\operatorname {Im}\Sigma (k,\omega _b)<0$ is not possible since this would lead to a non vanishing average Green function for $t<0$, thus violating causality.

In summary, Eqs. (33) and (39) show that it is possible to define a scattering mean-free path $\ell _s$ and a scattering mean-free time $\tau _s$ for a space and time dependent disorder. The reason for this is that self-energy $\Sigma$ is a simple multiplicative function, even when the scattering potential is an operator. Moreover, we clearly see from the Dyson Eq. (25) that there is no change in frequency or wavevector during propagation of the average field. For a monochromatic source at frequency $\omega$, this means that the average field propagates at $\omega$ and the scattering mean-free path can be defined for a fixed frequency $\omega$. Similarly, for a source at a fixed wavevector $k$, the average field evolves at the same $k$ and the scattering mean-free time can be defined for this fixed wavevector $k$. This behavior is typical of an average (or ballistic) field, and is observed for example in dynamic multiple scattering (or diffusing-wave spectroscopy) where the Doppler shift vanishes for the average field [23].

3. Disorder with a Gaussian correlation in space and time

To get explicit expressions for the scattering mean-free path $\ell _s$ and mean-free time $\tau _s$, we need to define a specific model of disorder. A canonical choice is that of a disorder with Gaussian correlation in both space and time, which allows to derive analytical expressions that can be easily compared to numerical simulations. This comparison is a relevant test of validity of the pertubation theory developed above.

3.1 Practical expressions for $\ell _s$ and $\tau _s$

In the weak-scattering regime, we can derive expressions for $\ell _s$ and $\tau _s$ restricted to the leading term in the perturbative expansion of the self-energy. The self-energy reads as

(40)$$\Sigma(x-x',t-t')=\frac{1}{c^4} \frac{\partial^2}{\partial t^2}\left\{\left\langle\delta\epsilon(x,t) \frac{\partial^2}{\partial t'^2}\left[G_b(x-x',t-t')\right] \delta\epsilon(x',t')\right\rangle\right\} \, ,$$

which can be reorganized in the form

(41)$$\Sigma(x-x',t-t')=\frac{1}{c^4} \frac{\partial^2}{\partial t^2}\left\{ \frac{\partial^2}{\partial t'^2}\left[G_b(x-x',t-t')\right] \left\langle\delta\epsilon(x,t)\delta\epsilon(x',t')\right\rangle \right\} \, ,$$

in which the correlation function of disorder appears explicitly. We now assume that this correlation function factorizes into two components, in the form

(42)$$\left\langle\delta\epsilon(x,t)\delta\epsilon(x',t')\right\rangle=\alpha(x-x')\beta(t-t') \, ,$$

meaning that short-range correlations may exist in the space or time dependence, with cross space-time correlations excluded. Plugging this expression into Eq. (41), and taking the Fourier transform, leads to

(43)$$\Sigma(k,\omega)=\frac{\omega^2}{c^4}\int (\omega-\omega')^2\mathcal{A}(k,\omega-\omega')\beta(\omega')\frac{\mathrm{d}\omega'}{2\pi}$$

where $\mathcal {A}(k,\omega )$ is the Fourier transform of $\mathcal {A}(x,t)=G_b(x,t)\alpha (x)$. We note that for a pure static disorder, with $\beta (t-t')=1$, we would recover the standard result involving the spatial correlation function of disorder and the Green function, namely $\beta (\omega ')=2\pi \delta (\omega ')$ and $\Sigma (k,\omega )=\omega ^4/c^4\mathcal {A}(k,\omega )$ [24].

The gaussian-correlated disorder model amounts to considering that

(44)$$\alpha(x-x')=A\exp\left[-\frac{(x-x')^2}{2\ell^2}\right] \quad;\quad \beta(t-t')=B\exp\left[-\frac{(t-t')^2}{2\tau^2}\right]$$

where $\ell$ and $\tau$ are the correlation length and time of disorder, respectively, and $A$ and $B$ are amplitudes of the correlation functions. With this model, we find that

(45)$$\begin{array}{r} \operatorname{Im}\Sigma(k,\omega)= \frac{AB\ell\tau\omega^4 v^3}{8c^4\left(\ell^2+\tau^2v^2\right)^{3/2}}e^{-\eta^2} \left\{ \sqrt{2\pi}\left(\omega\tau^2v-k\ell^2\right)e^{\xi^2} +2\sqrt{2\pi}\left(\omega\tau^2v+k\ell^2\right)e^{2k\omega\ell^2/v+\xi^2} \right.\\\left. +2\sqrt{\ell^2+\tau^2v^2} +\sqrt{2\pi}\left(k\ell^2-\omega\tau^2v\right)\operatorname{erf}(\xi)e^{\xi^2} \right\} -\frac{AB\ell\omega^4\tau v^3}{4c^4\ell^2\sqrt{\pi}\left(\omega+kv\right)}e^{-\eta^2} \left\{ \omega G_{0,0}^{2,1}\left(-\xi,\frac{1}{2}\left|\begin{matrix}1\\1/2,1\end{matrix}\right.\right) \right.\\\left. +\frac{2v^2}{\ell^2\left(\omega+kv\right)} G_{0,0}^{2,1}\left(-\xi,\frac{1}{2}\left|\begin{matrix}1\\1,3/2\end{matrix}\right.\right) \right\} \end{array}$$

where

(46)$$\eta=\frac{\ell(\omega+kv)}{\sqrt{2}v} \quad;\quad \xi=\frac{\ell^2(\omega+kv)}{\sqrt{2}v\sqrt{\ell^2+\tau^2v^2}}$$

and

(47)$$\begin{aligned}&G_{p,q}^{m,n}\left(z,r\left|\begin{matrix}a_1,\ldots,a_n,a_{n+1},\ldots,a_p\\b_1,\ldots,b_m,b_{m+1},\ldots,b_q\end{matrix}\right.\right)\\ &\qquad=\frac{r}{2i\pi}\int_{\gamma}\frac {\Gamma(1-a_1-rs)\ldots\Gamma(1-a_n-rs)\Gamma(b_1+rs)\ldots\Gamma(b_m+rs)} {\Gamma(a_{n+1}+rs)\ldots\Gamma(a_p+rs)\Gamma(1-b_{m+1}-rs)\ldots\Gamma(1-b_q-rs)} z^{{-}s}\mathrm{d} s \end{aligned}$$

is the generalized Meijer G-function, $\gamma$ being an appropriate path in the complex plane, and $\Gamma$ is the Gamma function. To have a practical expression of the scattering mean-free path, we use the on-shell approximation in Eq. (45), which leads to

(48)$$\begin{aligned} \operatorname{Im}\Sigma(k_b,\omega)&=\frac{AB\ell\tau\omega^2 v^2}{4c^4}\left\{ \omega\sqrt{\frac{2\pi}{\mathcal{W}}}+\frac{v}{\mathcal{W}}e^{-\mathcal{Z}} +\frac{\omega}{\mathcal{W}}\sqrt{\frac{\pi}{2\mathcal{W}}}e^{-\mathcal{X}} \left(\ell^2-v^2\tau^2\right)\left[\operatorname{erf}\left(\sqrt{\mathcal{Y}}\right)-1\right] \right.\\ &\left. +\frac{\omega}{\mathcal{W}}\sqrt{\frac{2\pi}{\mathcal{W}}}e^{-\mathcal{X}} \ell^2\left[{-}2+Q\left(-\frac12,\mathcal{Y}\right)\right] -\omega\sqrt{\frac{\pi}{2\mathcal{W}}}e^{-\mathcal{X}} \left[{-}2+Q\left(\frac12,\mathcal{Y}\right)\right] \right\} \end{aligned}$$

where

(49)$$\mathcal{W}=\epsilon_r\ell^2+v^2\tau^2 \quad;\quad \mathcal{X}=\frac{2\ell^2\omega^2\tau^2}{\mathcal{W}} \quad;\quad \mathcal{Y}=\frac{2\ell^4\omega^2}{v^2\mathcal{W}} \quad;\quad \mathcal{Z}=\frac{2\ell^2\omega^2}{v^2}$$

and $Q$ is the regularized incomplete gamma function defined by $Q(a,z)=\Gamma (a,z)/\Gamma (a)$, with $\Gamma (a,z)$ the incomplete gamma function. Equation (48) together with Eq. (33) provide the expression of the scattering mean-free path $\ell _s$ for a spatio-temporal gaussian-correlated disorder. It is also interesting to extract the expression of the scattering mean-free path in the limit where the time disorder vanishes (i.e., $\tau \to \infty$), which gives

(50)$$\frac{1}{\ell_{s,\tau\to\infty}(\omega)}=\frac{AB\ell\omega^2v^2}{2c^4}\sqrt{\frac{\pi}{2}}\left[ 1+\exp\left(-\frac{2\ell^2\omega^2}{v^2}\right) \right].$$

Similarly, the expression in the limit of a vanishing space disorder (i.e., $\ell \to \infty$) is

(51)$$\frac{1}{\ell_{s,\ell\to\infty}(\omega)}=\frac{AB\tau\omega^2 v^3}{2c^4}\sqrt{\frac{\pi}{2}}\left[ 1-\exp\left({-}2\tau^2\omega^2\right) \right].$$

To get an expression for the scattering mean-free time $\tau _s$ in the presence of space and time disorder, we have to make use of Eq. (39) together with Eq. (48) where all occurrences of the variable $\omega$ are replaced by $kv$. The limited cases are given by

(52)$$\frac{1}{\tau_{s,\tau\to\infty}(k)}=\frac{AB\ell k^2v^5}{2c^4}\sqrt{\frac{\pi}{2}}\left[ 1+\exp\left({-}2\ell^2k^2\right) \right]$$

for a vanishing time disorder and

(53)$$\frac{1}{\tau_{s,\ell\to\infty}(k)}=\frac{AB\tau k^2 v^6}{2c^4}\sqrt{\frac{\pi}{2}}\left[ 1-\exp\left({-}2\tau^2k^2v^2\right) \right]$$

for a vanishing space disorder.

3.2 Numerical simulations

In this section we compare the predictions of the theoretical model with numerical simulations performed without approximations. The first step consists in generating numerically an ensemble of configurations of disorder [i.e., of $\epsilon (x,t)$] that will be used to perform an ensemble average. The statistics of $\epsilon (x,t)$ has to satisfy Eq. (42), which is the only assumption on the model of disorder in the theory. One way of achieving this is to consider the particular case of a permittivity in which the space and time dependences factorize. We choose a permittivity in the form $\epsilon (x,t)=1+\delta \epsilon (x)\delta \epsilon (t)$, with $\epsilon _b=1$ for the sake of simplicity, and $\delta \epsilon (x)$ and $\delta \epsilon (t)$ statistically independent. In this case, we immediately find that

(54)$$\begin{aligned}\left\langle\delta\epsilon(x,t)\delta\epsilon(x',t')\right\rangle=\left\langle\delta\epsilon(x)\delta\epsilon(t)\delta\epsilon(x')\delta\epsilon(t')\right\rangle &=\left\langle\delta\epsilon(x)\delta\epsilon(x')\right\rangle\left\langle\delta\epsilon(t)\delta\epsilon(t')\right\rangle\\ & \qquad =\alpha(x-x')\beta(t-t'). \end{aligned}$$

Under this assumption, we only have to generate two independent one-dimensional disorders for $\delta \epsilon (x)$ and $\delta \epsilon (t)$, with gaussian correlation functions. Let us illustrate this process with space disorder with the correlation function $\alpha$. We consider a finite-size medium with size $L$, and we discretize the space into $N_x$ points $x_m$ in the interval $[-L/2,L/2]$, with a step $\Delta x=x_2-x_1$. Next, we generate a white-noise gaussian disorder [standard normal distribution $\mathcal {N}(0,1)$] that is finally convolved with

(55)$$f_m=\left(\frac{2}{\pi}\right)^{1/4}\sqrt{A\ell\Delta x}\exp\left[-\frac{x_m^2}{\ell^2}\right] \, ,$$

which gives one realization $\delta \epsilon _m$. Restarting the process allows us to generate a set of disorder configurations. After averaging, the correlation function tends to the function $\alpha$, as expected. The same process can be followed to generate configurations of the time disorder, the time interval $[0,T]$ being discretized into $N_t$ points $t_n$ with a step size $\Delta t$. An example of disorder is plotted in Fig. 2 (a) together with a comparison between the numerical and theoretical correlation function in Fig. 2 (b).

Fig. 2. (a) Example of spatial disorder at a fixed time $t$. (b) Comparison between the disorder correlation function in space computed numerically and that given by Eq. (44). The parameters are: $\sqrt {AB}=2{\times}10^{-2}$ and $\epsilon _b=1$. $1680$ disorder configurations are used to perform the statistical average.

Download Full Size | PDF

We now briefly describe the numerical resolution of the wave equation for a given configuration of disorder. We need to solve Eq. (1) with the boundary conditions $E(-L/2,t)=E(L/2,t)=0$, and the initial condition $E(x,0)=0$. The source term $S$ depends on the type of situation to be addressed. To compute the spatial evolution of the field, in order to estimate the scattering mean-free path, we choose

(56)$$S(x,t)=S_{\omega_0}(x,t)=s(t)\delta\left(x\right)e^{{-}i\omega_0 t}$$

which corresponds to a point source oscillating at a given frequency $\omega _0$. To avoid numerical artifacts due to a discontinuous source term in time, we apply a $C^{\infty }$ pseudo step function given by

(57)$$s(t)=\begin{cases} 0 &\text{if}\; t < 0, \\ 1 & \text{if}\; t > t_r, \\ \exp\left\{-\left[1-(t-t_r)^2/t_r^2\right]^{{-}1}+1\right\} &\text{otherwise,} \end{cases}$$

where $t_r$ is the rising time. To estimate the temporal evolution of the field, in order to compute the scattering mean-free time, we use a source term of the form

(58)$$S(x,t)=S_{k_0}(x,t)=d(t)e^{ik_0 x}$$

where $d$ is a $C^{\infty }$ pseudo Dirac delta function given by

(59)$$d(t)=\frac{s(t)-s(t-t_r)}{t_r},$$

also chosen to avoid numerical artifacts. This source term corresponds to a temporal pulse oscillating in space with a fixed wavevector $k_0$.

To solve the wave equation, we simply discretize it in space and time, with the numerical scheme

(60)$$E_{m,n+1}=\frac{2\epsilon_{m,n}E_{m,n}-\epsilon_{m,n-1}E_{m,n-1}}{\epsilon_{m,n+1}} +\frac{c^2\Delta t^2}{\epsilon_{m,n+1}}\left[ S_{m,n} +\frac{E_{m+1,n}+E_{m-1,n}-2E_{m,n}}{\Delta x^2} \right]$$

where the first indices ($m-1$, $m$, $m+1$) correspond to space discretization, and the second indices ($n-1$, $n$, $n+1$) to time discretization. The Dirac delta function in the source term is discretized using a Kronecker delta (i.e., $\delta _{m,m_0}/\Delta x$ where $m_0$ is the index corresponding to $x=0$). As for any finite-difference scheme, the Courant–Friedrichs–Lewy condition must be fulfilled to ensure numerical convergence and stability (i.e., $\Delta t\le \Delta x/c$). The resolution is performed for each disorder configuration of the ensemble, allowing us to estimate the ensemble averaged electric field.

Fig. 3. Intensity of the average field versus the normalized space variable $k_0x$, with $k_0=\omega _0/c$. This intensity is computed numerically (red solid line) and analytically (blue solid line for the full model, black dotted line for the model taking into account the space disorder only, i.e. $\tau \to \infty$, and green dotted line for the model taking into account the time disorder only, i.e. $\ell \to \infty$). The plot corresponds to the normalized time $\omega _0t=4000$. The parameters are: $k_0L=8000$, $k_0\ell =4$, $\omega _0\tau =4$, $\sqrt {AB}=2{\times}10^{-2}$, $\epsilon _b=1$ and $\omega _0t_r=100$. $10^4$ disorder configurations are used to perform the statistical average.

Download Full Size | PDF

Let us start with the spatial evolution of the field, with the source term $S_{\omega _0}$. We plot in Fig. 3 the intensity of the average field obtained from the full numerical simulation and from the analytical expressions, for the parameters given in the figure caption. The numerical result of $\left |\left \langle E(x,t)\right \rangle \right |^2$ at a fixed long time $t$ is compared to the square modulus of the average Green function (i.e., $\left |\left \langle G(x,\omega _0)\right \rangle \right |^2$ given by Eq. (32), with $k_r=k_b$). Excellent quantitative agreement is observed, which supports the validity of the theoretical model for the scattering mean-free path $\ell _s$. We also see that taking into account the spatial disorder only does not lead to an accurate result. The full model given by Eqs. (33) and (48) is needed to provide a relevant prediction, showing that the time dependence of the disorder clearly affects the spatial attenuation of the field. We also note that the scattering mean-free path is larger for the full disorder model than for spatial disorder model only, meaning that adding time disorder reduces the effect of scattering from space disorder. This result may look counter-intuitive, but it is because energy is not conserved in the presence of time disorder.

Fig. 4. Intensity of the average field as a function of the normalized time variable $\omega _0t$ with $\omega _0=k_0c$. This intensity is computed numerically (red solid line) by applying Eq. (61) and analytically (blue solid line for the full model, black dotted line for the model taking into account the space disorder only, i.e. $\tau \to \infty$, and green dotted line for the model taking into account the time disorder only, i.e. $\ell \to \infty$). The plot corresponds to a fixed normalized $k_0x=0$. The parameters are: $k_0L=8000$, $k_0\ell =4$, $\omega _0\tau =4$, $\sqrt {AB}=2{\times}10^{-2}$, $\epsilon _b=1$, $\omega _0t_r=0.1$ and $\omega _0T=30$. $10^4$ disorder configurations are used to perform the statistical average. Short times are not represented since for $t<T$, the averaging procedure given by Eq. (61) leads to an oscillating signal because of the Heaviside step function at $t=0$.

Download Full Size | PDF

Next, we study the temporal evolution of the field with the source term $S_{k_0}$. In order to compare the numerical results to the square of the average Green function (i.e., $\overline {\left \langle G(k_0,t)\right \rangle ^2}$, with $\omega _r=\omega _b$), we first compute numerically the average field $\left \langle E(x,t)\right \rangle$ for a fixed $x=0$. Then, we take the square modulus and perform a rolling average over a time window with width $T$ satisfying $2\pi /\omega _b\ll T\ll \tau _s$. This eliminates rapid oscillations and keeps the decaying envelope that depends on the scattering mean-free time $\tau _s$. We obtain

(61)$$I(x,t)=\int_{-\infty}^{+\infty} w(t-t')\left|\left\langle E(x,t')\right\rangle\right|^2\mathrm{d} t'$$

where $w$ is a rectangular function of width $T$ and amplitude $1/T$. The comparison is plotted in Fig. 4. Again, we obtain excellent quantitative agreement between the numerical simulation and the analytic expressions, supporting the theoretical model for the scattering mean-free time $\tau _s$. We also observe that the full theoretical model taking into account both space and time disorder is required to correctly predict the time decay of the intensity.

4. Conclusion

In conclusion, we have studied the behavior in space and time of the averaged field propagating in a medium with both space and time disorders. We have developed a multiple scattering theory that predicts the space and time decay of the average field, and allows to derive practical expressions of the scattering mean-free path $\ell _s$ and mean-free time $\tau _s$ in the weak-scattering regime. The model has been compared to exact numerical simulations, showing quantitative agreement in the particular case of a spatio-temporal gaussian-correlated disorder with no space-time cross correlation. Counter-intuitively, in this regime the introduction of a time disorder on top of a space disorder tends to reduce the scattering strength, even in the absence of cross correlation between the two types of disorders. However, this theory does not seem to predict the existence of a situation where the attenuation by scattering is totally cancelled out, at least in the case of Gaussian-correlated disorder. The theory developed in this work and the results bring a brick in the widely open field of waves in complex space and time varying media. In particular, the next step will be to study the case of the average intensity, whose behavior is known to be very different from that of the average electric field in the presence of disorder.

A. Calculation of the Green function $G_b$

The 1D scalar Green function of the wave equation in the reference medium described by its relative permittivity $\epsilon _b$ is given by Eq. (4) which reads in the Fourier domain

(62)$$\left(k^2-\frac{\omega^2}{v^2}\right)G_b(k,\omega)=1$$

where we recall that $v=c/\sqrt {\epsilon _b}$. $k$ and $\omega$ are the dual variables for $x$ and $t$ respectively. In terms of distributions (the Green function is rigorously a distribution), the inversion of this equation leads to

(63)$$G_b(k,\omega)=\operatorname{PV}\left[\frac{1}{k^2-\omega^2/v^2}\right] +\lambda\delta\left(k-\frac{\omega}{v}\right) +\mu\delta\left(k+\frac{\omega}{v}\right)$$

where $\lambda$ and $\mu$ are constants that should be determined in order to fulfil the boundary conditions in space and time. For that purpose, we consider first the inverse Fourier transform in time which gives

(64)$$G_b(k,t)=\operatorname{PV}\int_{-\infty}^{+\infty}\frac{e^{{-}i\omega t}}{\omega_b^2/v^2-\omega^2/v^2}\frac{\mathrm{d}\omega}{2\pi} +\frac{\lambda v}{2\pi}e^{{-}ikvt} +\frac{\mu v}{2\pi}e^{ikvt}$$

where we recall that $\omega _b=kv$. To compute the first term, we apply the residue theorem and consider two cases. If $t<0$, we use the contour described in Fig. 5 (a). The semicircle in the upper plane is chosen in order to apply Jordan’s lemma. This leads to

(65)$$G_b(k,t<0)={-}\frac{v^2}{2\omega_b}\sin(\omega_b t) +\frac{\lambda v}{2\pi}e^{{-}i\omega_b t} +\frac{\mu v}{2\pi}e^{i\omega_b t}.$$

The causality requires that $G_b(k,t<0)=0$ which leads to $-\lambda =\mu =\pi /(2ik)$. For $t>0$, we use the contour described in Fig. 5 (b). This finally gives

(66)$$G_b(k,t)=\frac{\operatorname{H}(t)v^2}{\omega_b}\sin\left(\omega_b t\right).$$

We consider now the inverse Fourier transform in space with the values of $\lambda$ and $\mu$ determined above. This gives

(67)$$G_b(x,\omega)=\operatorname{PV}\int_{-\infty}^{+\infty}\frac{e^{ikx}}{k^2-k_b^2}\frac{\mathrm{d} k}{2\pi} +\frac{i}{2k_b}\cos\left(k_bx\right)$$

where we recall that $k_b=\omega /v$. Again, we consider two cases to compute the first term. If $x>0$, we use the contour described in Fig. 5 (c). This leads to

(68)$$G_b(x>0,\omega)=\frac{i}{2k_b}\left[\cos\left(k_bx\right)+i\sin\left(k_bx\right)\right] =\frac{i}{2k_b}\exp\left(ik_bx\right).$$

If $x<0$, we use the contour described in Fig. 5 (d) which gives

(69)$$G_b(x<0,\omega)=\frac{i}{2k_b}\left[\cos\left(k_bx\right)-i\sin\left(k_bx\right)\right] =\frac{i}{2k_b}\exp\left({-}ik_bx\right).$$

The two previous results combined give

(70)$$G_b(x,\omega)=\frac{i}{2k_b}\exp\left(ik_b|x|\right).$$

It is interested to note that this last expression automatically fulfilled the outgoing wave condition thanks to causality.

Fig. 5. Various integration contours used to compute the Green function of the reference medium. (a) and (b) are used for the inverse Fourier transform for negative and positive times respectively. (c) and (d) are used for the inverse Fourier transform for positive and negative positions respectively.

Download Full Size | PDF

Funding

French Government.

Acknowledgements

This work has received support under the program “Investissements d’Avenir” launched by the French Government.

Disclosures

The authors declare no conflicts of interest.

Data availability

The data that supports the results of this study are available within the article.

References

1. P. Sheng, Introduction to Wave Scattering, Localization and Mesoscopic Phenomena (Springer, 2006).

2. P. Sebbah, ed.,Waves and Imaging through Complex Media (Springer, 2001).

3. S. Gigan, O. Katz, H. B. de Aguiar, et al., “Roadmap on wavefront shaping and deep imaging in complex media,” J. Phys. Photonics 4(4), 042501 (2022). [CrossRef]

4. R. Carminati and J. Schotland, Principles of Scattering and Transport of Light (Cambridge University Press, 2021).

5. C. Caloz and Z.-L. Deck Léger, “Spacetime Metamaterials—Part II: Theory and Applications,” IEEE Trans. Antennas Propag. 68(3), 1583–1598 (2020). [CrossRef]

6. E. Lustig, Y. Sharabi, and M. Segev, “Topological aspects of photonic time crystals,” Optica 5(11), 1390 (2018). [CrossRef]

7. Y. Sharabi, E. Lustig, and M. Segev, “Disordered photonic time crystals,” Phys. Rev. Lett. 126(16), 163902 (2021). [CrossRef]

8. S. Saha, O. Segal, C. Fruhling, et al., “Photonic time crystals: a materials perspective [Invited],” Opt. Express 31(5), 8267 (2023). [CrossRef]

9. R. Tirole, S. Vezzoli, E. Galiffi, et al., “Double-slit time diffraction at optical frequencies,” Nat. Phys. 19(7), 999–1002 (2023). [CrossRef]

10. F. Zangeneh-Nejad and R. Fleury, “Active times for acoustic metamaterials,” Rev. Phys. 4, 100031 (2019). [CrossRef]

11. V. Bacot, M. Labousse, A. Eddi, et al., “Time reversal and holography with spacetime transformations,” Nat. Phys. 12(10), 972–977 (2016). [CrossRef]

12. V. Bacot, G. Durey, A. Eddi, et al., “Phase-conjugate mirror for water waves driven by the Faraday instability,” Proc. Natl. Acad. Sci. U.S.A. 116(18), 8809–8814 (2019). [CrossRef]

13. A. Akbarzadeh, N. Chamanara, and C. Caloz, “Inverse prism based on temporal discontinuity and spatial dispersion,” Opt. Lett. 43(14), 3297 (2018). [CrossRef]

14. V. Pacheco Peña and N. Engheta, “Antireflection temporal coatings,” Optica 7(4), 323 (2020). [CrossRef]

15. Y. Sharabi, A. Dikopoltsev, E. Lustig, et al., “Spatiotemporal photonic crystals,” Optica 9(6), 585 (2022). [CrossRef]

16. R. Carminati, H. Chen, R. Pierrat, et al., “Universal statistics of waves in a random time-varying medium,” Phys. Rev. Lett. 127(9), 094101 (2021). [CrossRef]

17. B. Apffel, S. Wildeman, A. Eddi, et al., “Experimental implementation of wave propagation in disordered time-varying media,” Phys. Rev. Lett. 128(9), 094503 (2022). [CrossRef]

18. J. Kim, D. Lee, S. Yu, et al., “Unidirectional scattering with spatial homogeneity using correlated photonic time disorder,” Nat. Phys. 19(5), 726–732 (2023). [CrossRef]

19. S. M. Rytov, Y. A. Kravtsov, and V. I. Tatarskii, Principles of Statistical Radiophysics, vol. 4 (Springer-Verlag, 1989).

20. E. Akkermans and G. Montambaux, Mesoscopic Physics of Electrons and Photons (Cambridge University Press, 2007).

21. M. Born and E. Wolf, Principles of Optics (Cambridge University Press, 1999).

22. I. Baydoun, D. Baresch, R. Pierrat, et al., “Scattering mean free path in continuous complex media: Beyond the Helmholtz equation,” Phys. Rev. E 92(3), 033201 (2015). [CrossRef]

23. R. Pierrat, “Transport equation for the time correlation function of scattered field in dynamic turbid media,” J. Opt. Soc. Am. A 25(11), 2840 (2008). [CrossRef]

24. K. Vynck, R. Pierrat, R. Carminati, et al., “Light in correlated disordered media,” Rev. Mod. Phys. 95(4), 045003 (2023). [CrossRef]

Multiple scattering theory in one dimensional space and time dependent disorder: average field [Invited]

Abstract

1. Introduction

2. Multiple scattering theory for space-time disorder

2.1 Lippmann-Schwinger equation

2.2 Born series and Dyson equation

2.3 Weak-scattering regime

3. Disorder with a Gaussian correlation in space and time

3.1 Practical expressions for $\ell _s$ and $\tau _s$

3.2 Numerical simulations

4. Conclusion

A. Calculation of the Green function $G_b$

Funding

Acknowledgements

Disclosures

Data availability

References

Data availability

Cited By

Figures (5)

Equations (67)

Optical Materials Express