Multiple scattering theory in one dimensional space and time dependent disorder: Average field

We theoretically study the propagation of light in one-dimensional space- and time-dependent disorder. The disorder is described by a fluctuating permittivity $\epsilon(x,t)$ exhibiting short-range correlations in space and time, without cross correlation between them. Depending on the illumination conditions, we show that the intensity of the average field decays exponentially in space or in time, with characteristic length or time defining the scattering mean-free path $\ell_s$ and the scattering mean-free time $\tau_s$. In the weak scattering regime, we provide explicit expressions for $\ell_s$ and $\tau_s$, that are checked against rigorous numerical simulations.


I. 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][7][8][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][14][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 onedimensional space disorder combined to a time modulation, resulting in a medium described by a fluctuating dielectric function ϵ(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 corre-lation 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 ℓ s and a scattering mean-free time τ s , and to provide explicit expressions in the weak scattering regime.The paper is organized as follows: In Sec.II, we develop the theory that extends the standard multiple scattering theory to a situation with both space and time disorders.We provide expressions for ℓ s and τ s using a perturbative approach.In Sec.III, 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.

II. 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 ℓ s and τ 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 ϵ(x, t), an electric field linearly polarized along the y-direction obeys the equation 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.

A. Lippmann-Schwinger equation
The first step of the derivation consists in defining a homogeneous reference (or background) medium with permittivity ϵ b .The reference field E b in this medium satisfies The choice of ϵ b will be specified later, with the constraint that it should be close to the typical value of ϵ(x, t) to ensure the accuracy of the perturbative approach.Subtracting Eq. ( 2) from Eq. ( 1) leads to 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 where δ 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 where v = c/ √ ϵ b is the phase velocity in the reference medium and 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, ω) and in the (k, t) domains are thus required for problems (1) and (2), respectively.They are given by where H is the Heaviside step function, k b = ω/v and ω 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 +∞ (respectively −∞) for x > 0 (respectively x < 0).In time, this condition is equivalent to causality which translates into the presence of the 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 The total field E = E b + E s obeys the integral equation 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 where the bullets are to be uderstood as the quantity on which the operator acts.In this formalism, the Lippmann-Schwinger equation ( 8) can be rewritten as 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 V.

B. 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 ϵ(x, t, )), we first expand Eq. (10) in the form (11) which is known as the Born series.Performing a statistical ensemble average, we find that where ⟨. ..⟩ denotes the average value.The problem now reduces to the computation of terms of the form ⟨V(G b V) n ⟩.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 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 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 and similarly for higher-order terms.
The key idea to obtain an equation for the average field consists in defining a new operator S containing all non-factorizable terms, i.e.
With this definition, Eq. ( 12) can be factorized in the form 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 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 We now need to make a choice for the reference medium.
Taking ϵ b = ⟨ϵ(x, t)⟩ 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 δϵ(x, t) = ϵ(x, t) − ⟨ϵ(x, t)⟩ = ϵ(x, t) − ϵ b we find that For the second order term, we obtain where we have used the relationship ⟨δϵ(x, t)δϵ(x ′ , t ′ )⟩ c = ⟨δϵ(x, t)δϵ(x ′ , t ′ )⟩, 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 The second derivation relates to t ′ and applies to δϵ(x ′ , t ′ )•.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 Similar transformations can be performed on the higher order terms in S, but are not written here since they will not be useful in practice.Finallly, S can be written as where Σ is the self-energy and is here a simple multiplicative function (not an operator).Using the self-energy, the average field can be written which is the integral form of the Dyson equation.

C. 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) = δ(x)δ(t) in an infinite medium.In this case, E b (x, t) = G b (x, t) and ⟨E(x, t)⟩ = ⟨G(x, t)⟩, with G the Green function of the medium in the presence of disorder.We assume statistical homogeneity in space and time, such that Σ(x ′ , x ′′ , t ′ , t ′′ ) only depends on x ′ − x ′′ and t ′ − t ′′ .In these conditions, Eq. ( 23) reduces to This equation can be solved by performing a space-time Fourier transform, which leads to (26) This expression of the average Green function will be used to derive expressions for the scattering mean-free path ℓ s and scattering mean-free time τ s .We start by considering a monochromatic source term, oscillating at a frequency ω, and we focus on the spatial behavior of the average field given by The computation of this inverse Fourier transform requires additional hypotheses.Considering the weakscattering regime defined by the condition |Σ(k, ω)| ≪ k 2 b , the self-energy has a significant contribution only when k ≃ ±k b .Assuming that the disorder is statistically isotropic, we also have Σ(k, ω) = Σ(−k, ω).As a result, the self-energy can be taken on-shell for k = k b in Eq. ( 27).Under this assumption, Eq. (27) becomes 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 Assuming that Im Σ(k b , ω) > 0, which will be justified below, we obtain For x < 0, we use the contour presented in Fig. 1 (b) and we obtain which finally leads to 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ℓ s ), with ℓ s the scattering meanfree path, and k r = n r ω/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 which will be the expression used later for comparison with numerical simulations.In the weak scattering regime, we have We also note that the approximation k r ≃ k b holds in the weak-scattering regime.A more refined expression would involve the real part of the self-energy.
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 The weak-scattering regime amounts to assuming that Under this assumption, the selfenergy takes significant values for ω ≃ ±ω b .For statistically isotropic disorder, such that Σ(k, ω) = Σ(−k, ω), and making use of the fact that Σ(x, t) is real valued, we find that As a result, the self-energy can be replaced by Σ(k, ω b ) * in Eq. (34) in the vicinity of −ω b , and by Σ(k, ω b ) in the vicinity of ω b .This is the counterpart of the onshell 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 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 Defining ω e = ω r − i/(2τ s ), with τ s the scattering meanfree time, we finally obtain 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π/ω r ≪ T ≪ τ s .This leads to which will be the expression used for comparison with numerical simulations.Under the weak-scattering approximation, the scattering mean-free time is We note that ω r ≃ ω b in the weak-scattering regime.We also stress that having Im Σ(k, ω 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 ℓ s and a scattering mean-free time τ s for a space and time dependent disorder.The reason for this is that self-energy Σ is a simple multiplicative function, even when the scattering potential is an operator.Moreover, we clearly see from the Dyson equation (25) that there is no change in frequency or wavevector during propagation of the average field.For a monochromatic source at frequency ω, this means that the average field propagates at ω and the scattering mean-free path can be defined for a fixed frequency ω.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].

III. DISORDER WITH A GAUSSIAN CORRELATION IN SPACE AND TIME
To get explicit expressions for the scattering mean-free path ℓ s and mean-free time τ 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.

A. Practical expressions for ℓs and τs
In the weak-scattering regime, we can derive expressions for ℓ s and τ s restricted to the leading term in the perturbative expansion of the self-energy.The selfenergy reads as which can be reorganized in the form in which the correlation function of disorder appears explicitly.We now assume that this correlation function factorizes into two components, in the form 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 where A(k, ω) is the Fourier transform of A(x, t) = G b (x, t)α(x).We note that for a pure static disorder, with β(t − t ′ ) = 1, we would recover the standard result involving the spatial correlation function of disorder and the Green function, namely β(ω ′ ) = 2πδ(ω ′ ) and Σ(k, ω) = ω 4 /c 4 A(k, ω) [24].
The gaussian-correlated disorder model amounts to considering that where ℓ and τ 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 where and z −s ds (47) is the generalized Meijer G-function, γ being an appropriate path in the complex plane, and Γ 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 where and Q is the regularized incomplete gamma function defined by Q(a, z) = Γ(a, z)/Γ(a), with Γ(a, z) the incomplete gamma function.Equation (48) together with Eq. (33) provide the expression of the scattering meanfree path ℓ 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., τ → ∞), which gives (50) Similarly, the expression in the limit of a vanishing space disorder (i.e., ℓ → ∞) is To get an expression for the scattering mean-free time τ 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 ω are replaced by kv.The limited cases are given by for a vanishing time disorder and for a vanishing space disorder.

B. 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 ϵ(x, t)] that will be used to perfom an ensemble average.The statistics of ϵ(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 ϵ(x, t) = 1 + δϵ(x)δϵ(t), with ϵ b = 1 for the sake of simplicity, and δϵ(x) and δϵ(t) statistically independent.In this case, we immediately find that Under this assumption, we only have to generate two independent one-dimensional disorders for δϵ(x) and δϵ(t), with gaussian correlation functions.Let us illustrate this process with space disorder with the correlation function α.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 ∆x = x 2 − x 1 .Next, we generate a white-noise gaussian disorder [standard normal distribution N (0, 1)] that is finally convolved with which gives one realization δϵ m .Restarting the process allows us to generate a set of disorder configurations.After averaging, the correlation function tends to the function α, 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 ∆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).
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 which corresponds to a point source oscillating at a given frequency ω 0 .To avoid numerical artifacts due to a discontinuous source term in time, we apply a C ∞ pseudo step function given by 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 where d is a C ∞ pseudo Dirac delta function given by 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 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., δ m,m0 /∆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., ∆t ≤ ∆x/c).The resolution is performed for each disorder configuration of the ensemble, allowing us to estimate the ensemble averaged electric field.Let us start with the spatial evolution of the field, with the source term S ω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 |⟨E(x, t)⟩| 2 at a fixed long time t is compared to the square modulus of the average Green function (i.e., |⟨G(x, ω 0 )⟩| 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 meanfree path ℓ 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.
Next, we study the temporal evolution of the field with the source term S k0 .In order to compare the numerical results to the square of the average Green function (i.e., ⟨G(k 0 , t)⟩ 2 , with ω r = ω b ), we first compute numerically the average field ⟨E(x, t)⟩ 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π/ω b ≪ T ≪ τ s .This eliminates rapid oscillations and keeps the decaying envelope that depends on the scattering meanfree time τ s .We obtain 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 τ 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.

IV. 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 ℓ s and mean-free time τ 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.

AKNOWLEDGMENTS
This work has received support under the program "Investissements d'Avenir" launched by the French Government.

DISCLOSURES
The authors declare no conflicts of interest.The 1D scalar Green function of the wave equation in the reference medium described by its relative permittivity ϵ b is given by Eq. ( 4) which reads in the Fourier domain where we recall that v = c/ √ ϵ b .k and ω 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 where λ and µ 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 (A4) The causality requires that G b (k, t < 0) = 0 which leads to −λ = µ = π/(2ik).For t > 0, we use the contour described in Fig. 5  The two previous results combined give It is interested to note that this last expression automatically fulfiled the outgoing wave condition thanks to causality.

FIG. 1 .
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 Im Σ(k b , ω) > 0 for (a) and (b) and Im Σ(k, ω b ) > 0 for (c) and (d), as explained in the main text.

FIG. 2 .
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: √ AB = 2 × 10 −2 and ϵ b = 1.1680 disorder configurations are used to perform the statistical average.

FIG. 3 .
FIG.3.Intensity of the average field versus the normalized space variable k0x, with k0 = ω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. τ → ∞, and green dotted line for the model taking into account the time disorder only, i.e. ℓ → ∞).The plot corresponds to the normalized time ω0t = 4000.The parameters are: k0L = 8000, k0ℓ = 4, ω0τ = 4, √ AB = 2 × 10 −2 , ϵ b = 1 and ω0tr = 100.10 4 disorder configurations are used to perform the statistical average.

FIG. 4 .
FIG.4.Intensity of the average field as a function of the normalized time variable ω0t with ω0 = k0c.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. τ → ∞, and green dotted line for the model taking into account the time disorder only, i.e. ℓ → ∞).The plot corresponds to a fixed normalized k0x = 0.The parameters are: k0L = 8000, k0ℓ = 4, ω0τ = 4, √ AB = 2 × 10 −2 , ϵ b = 1, ω0tr = 0.1 and ω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.
Appendix A: Calculation of the Green function G b

G
b (k, t) = PV +∞ that ω 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 toG b (k, t < 0) = − v 2 2ω b sin(ω b t) + λv 2π e −iω b t + µv 2π e iω b t .

FIG. 5 .
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.