# Universal statistics of waves in a random time-varying medium

R. Carminati<EMAIL_ADDRESS>Institut Langevin, ESPCI Paris, PSL
University, CNRS, 1 rue Jussieu, F-75005 Paris, France  H. Chen  Institut
Langevin, ESPCI Paris, PSL University, CNRS, 1 rue Jussieu, F-75005 Paris,
France  R. Pierrat  Institut Langevin, ESPCI Paris, PSL University, CNRS, 1
rue Jussieu, F-75005 Paris, France  B. Shapiro  Department of Physics,
Technion-Israel Institute of Technology, Haifa IL-32000, Israel

###### Abstract

We study the propagation of waves in a medium in which the wave velocity
fluctuates randomly in time. We prove that at long times, the statistical
distribution of the wave energy is log-normal, with the average energy growing
exponentially. For weak disorder, another regime preexists at shorter times,
in which the energy follows a negative exponential distribution, with an
average value growing linearly with time. The theory is in perfect agreement
with numerical simulations, and applies to different kinds of waves. The
existence of such universal statistics bridges the fields of wave propagation
in time-disordered and space-disordered media.

Introduction. In recent years there has been a growing interest in space-time
metamaterials for electromagnetic Caloz2020  or acoustic waves Fleury2019 .
These are materials whose properties are modulated in space and time.
Homogeneous materials modulated only in time, referred to as “temporal” or
“pure-time”, offer new approaches for the control of waves, e.g. through the
design of active metasurfaces. They also stimulate basic studies in wave
physics. For example, it is known for electromagnetic waves that when the
dielectric function is suddenly changed from one value to another, a backward
propagating (time-reversed) wave appears Morgenthaler1958 . This phenomenon
has been recently put into a general framework, and demonstrated
experimentally with water waves Bacot2016 ; Bacot2019 . A periodic modulation
of the dielectric function has also been investigated Holberg1966 ; Zurita2009
, leading to the appearance of bands and gaps in the wave propagation constant
$k$, as well as topological phases Lustig2018 . Building on the analogy
between space and time, new approaches have been proposed for the control of
waves Caloz2018 ; Pacheco2020 ; Pacheco2020b . There is a lot to expect in the
interaction between waves and more complex temporal materials, including
random time-varying media, a domain that has remained unexplored to a large
extent.

The purpose of this Letter is the study of wave propagation in a medium with a
dielectric function $\varepsilon(t)$ fluctuating randomly in time. Our focus
is on electromagnetic waves, but the developed theory and the results
encompass other kind of waves, such as acoustic or water waves. The question
of time evolution of a pulse subjected to random“kicks”, due to sudden changes
in $\varepsilon(t)$, has been posed and studied in Ref. Sharabi2021 . It has
been shown that, after a sufficiently long time, the energy $U(t)$ of the
pulse increases exponentially. A similar regime has been found in a recent
study, in which water waves propagate in a disordered time-periodic lattice
Apffel2021 . This behavior suggests a connection with Anderson localization of
waves in a spatially disordered medium Anderson . Here we consider a general
model of a disordered time-varying medium, with an emphasis on weak disorder
(the criterion for weak disorder will be stated later). In this case, one can
develop a detailed analytical theory. The theory shows that at times larger
than a crossover time $\tau_{c}$, $\langle\ln U(t)\rangle$ becomes
proportional to $t$, with the brackets denoting the average value, in
agreement with experimental and numerical observations Sharabi2021 ;
Apffel2021 . Interestingly, there is an intermediate regime $\tau_{m}\ll
t\ll\tau_{c}$, with $\tau_{m}$ the microscopic time (defined as the typical
time of the modulation of $\varepsilon(t)$), in which the average energy
$\langle U(t)\rangle$ grows linearly with $t$. The full statistical
distribution of the energy $U$ can be determined in both regimes. In the
intermediate regime, the energy follows a negative exponential distribution.
For long times $t\gg\tau_{c}$, the statistics becomes log-normal, in agreement
with known results in one-dimensional wave transport.

General framework. We consider the propagation of electromagnetic waves in a
homogeneous, isotropic and nonmagnetic medium, described by its time-dependent
dielectric function $\varepsilon(t)$ such that the displacement and electric
fields are related by ${\bf D}({\bf r},t)=\varepsilon_{0}\varepsilon(t){\bf
E}({\bf r},t)$. When the displacement field has a single component $D$, and
depends only on one space coordinate $x$, it satisfies

$\frac{\partial^{2}D}{\partial
x^{2}}(x,t)-\frac{\varepsilon(t)}{c^{2}}\frac{\partial^{2}D}{\partial
t^{2}}(x,t)=0\,,$ (1)

together with appropriate boundary conditions. It will prove useful to perform
the analysis in $k$ space. The space Fourier transform, defined as
$D(k,t)=\int_{-\infty}^{+\infty}D(x,t)\exp(-ikx)dx$, satisfies

$\frac{\partial^{2}D}{\partial t^{2}}(k,t)+\Omega^{2}(t)D(k,t)=0\,,$ (2)

where $\Omega^{2}(t)=c^{2}k^{2}/\varepsilon(t)$. Equation (2), supplemented
with two initial conditions for $D(k,t)$ and its time derivative, constitutes
a Cauchy problem. We note that since $D(x,t)$ is real, $D(-k,t)=D^{*}(k,t)$,
where the superscript $*$ stands for complex conjugate. Therefore the analysis
can be limited to $k\geq 0$. We also point out that this description is not
limited to a fully homogeneous space. Indeed, the only requirement is
homogeneity along the propagation direction $x$. For example, the analysis
could apply to a waveguide filled with a homogeneous medium having
$\varepsilon(t)$ depending on time, with the plane wave replaced by a guided
wave with a given transverse profile.

Over a time interval in which $\varepsilon$ is a constant, the general
solution to Eq. (2) is of the form

$D(k,t)=D^{+}(k,t)+D^{-}(k,t)\,,$ (3)

with $D^{+}(k,t)\sim\exp(-i\Omega t)$ and $D^{-}(k,t)\sim\exp(i\Omega t)$,
corresponding to plane waves propagating in the positive and negative $x$
direction, respectively. In this study, the observable of interest is the
electromagnetic energy $U(t)=[2\varepsilon_{0}\varepsilon(t)]^{-1}\int{\bf
D}^{2}({\bf r},t)d^{3}r+(2\mu_{0})^{-1}\int{\bf B}^{2}({\bf r},t)d^{3}r$
Landau_Book . For a one-dimensional and linearly polarized field, the energy
can be rewritten as

$U(t)=\frac{1}{2\pi\varepsilon_{0}\varepsilon}\int_{-\infty}^{+\infty}\left[|D^{+}(k,t)|^{2}+|D^{-}(k,t)|^{2}\right]dk\,.$
(4)

Note that the electric and magnetic contributions to the energy contain
interference terms that exactly compensate, resulting in the simple expression
above.

Transfer matrix. In $k$ space, the time evolution of the field can be
described using transfer matrices. To model a medium with a dielectric
function $\varepsilon(t)$ changing randomly in time, we can take
$\Omega^{2}(t)$ to be a series of instantaneous kicks ($\delta$-kicks) on top
of a background value $\Omega^{2}_{b}$, as represented in Fig. 1. In this
model the kick strength $v_{j}$ and times $t_{j}$ are independent random
variables.

Figure 1:  Random chain of $\delta$ kicks, with
$\Omega^{2}(t)=\Omega_{b}^{2}+\sum_{j}v_{j}\,\delta(t-t_{j})$. The kick
strengths $v_{j}$ and the times $t_{j}$ are random variables. After each kick
the medium recovers the background value $\Omega_{b}$. The statistical
distribution of $v_{j}$ is independent of $j$.

The transfer matrix ${\bf M}_{j}$ of a single block $j$ (see Fig. 1) connects
the field $D_{j}$ at time $t_{j}+0$ to the field $D_{j-1}$ at time
$t_{j-1}+0$, in such a way that

$\begin{bmatrix}D^{+}_{j}\\\ D^{-}_{j}\end{bmatrix}={\bf
M}_{j}\begin{bmatrix}D^{+}_{j-1}\\\ D^{-}_{j-1}\end{bmatrix}\,.$ (5)

The analysis is performed at a fixed $k$, and the matrix elements depend on
$k$, but we drop the argument $k$ in the notations for simplicity. An explicit
calculation of the transfer matrix, detailed in the Supplemental Material SM ,
shows that

${\bf M}_{j}=\begin{bmatrix}a_{j}&b_{j}\\\ b_{j}^{*}&a_{j}^{*}\end{bmatrix}$
(6)

with

$\displaystyle a_{j}$ $\displaystyle=$
$\displaystyle(1-iu_{j})\exp(-i\theta_{j})$ $\displaystyle b_{j}$
$\displaystyle=$ $\displaystyle-iu_{j}\exp(i\theta_{j})\,,$ (7)

and $u_{j}=v_{j}/(2\Omega_{b})$. The matrix ${\bf M}_{j}$ has the following
properties:

$\displaystyle|a_{j}|^{2}\geq 1$ (8) $\displaystyle\mathrm{det}{\bf
M}_{j}=|a_{j}|^{2}-|b_{j}|^{2}=1\,.$ (9)

Moreover, due to property (9), it is easy to see that
$|D_{j}^{+}|^{2}-|D_{j}^{-}|^{2}=|D_{j-1}^{+}|^{2}-|D_{j-1}^{-}|^{2}$, showing
that the quantity $|D_{j}^{+}|^{2}-|D_{j}^{-}|^{2}$ is conserved. At time
$t<t_{1}$ (before the first kick), we can assume that the incident wave
propagates in the positive $x$ direction, with $D_{0}^{+}$ normalized such
that $|D_{0}^{+}|^{2}=1$. As a result, the following relation is satisfied for
all $j$:

$|D_{j}^{+}|^{2}-|D_{j}^{-}|^{2}=1\,.$ (10)

It is important to note that the subsequent analytical treatment is not
limited to the specific shape of modulation, like the $\delta$ kicks in Figs.
1, but is applicable to any modulation provided that $\Omega(t)$ recovers the
same background value $\Omega_{b}$ between the random kicks. For example we
show in the Supplemental Material SM  that a model based on rectangular time
barriers leads to a transfer matrix satisfying the same properties. In fact
the transfer matrix (6), with the properties (8-9) for its matrix elements, is
the most general $2\times 2$ transfer matrix MelloBook .

Statistical theory. We now develop a theoretical analysis of the statistical
properties of the quantity

$U_{j}=|D_{j}^{+}|^{2}+|D_{j}^{-}|^{2}\,,$ (11)

which, up to a prefactor that we omit, corresponds to the energy in the field
after $j$ kicks. It will prove useful to introduce new variables
$z_{j}=|D_{j}^{-}|^{2}$ and $\beta_{j}=|b_{j}|^{2}$. We note that due to
relations (9) and (10), we have $U_{j}=1+2z_{j}$ and
$|a_{j}|^{2}+|b_{j}|^{2}=1+2\beta_{j}$. From Eqs. (5) and (6), we find that

$\displaystyle\ln(1+2z_{j})=\ln(1+2z_{j-1})+\ln(1+2\beta_{j})$
$\displaystyle+\ln\left(1+2\gamma_{j}\,\frac{\sqrt{z_{j-1}(1+z_{j-1})}}{1+2z_{j-1}}\,\cos\Theta_{j}\right)\,,$
(12)

where $\Theta_{j}$ is the cumulative phase such that
$a_{j}b_{j}^{*}D_{j-1}^{+}D_{j-1}^{-*}=|a_{j}||b_{j}||D_{j-1}^{+}||D_{j-1}^{-}|\exp(i\Theta_{j})$
and $\gamma_{j}=2\sqrt{\beta_{j}(1+\beta_{j})}/(1+2\beta_{j})$. The recursion
relation (12) is our starting point in the analysis of the statistical
properties of the field energy.

We first analyze the behavior of the energy $U_{N}$ after $N$ kicks, in the
limit $N\to\infty$. In this limit, which is the same as $t\to\infty$, the
field amplitude $D(k,t)$ is expected to increase exponentially with time. A
more rigorous statement is that the Lyapunov exponent

$\lambda(k)=\lim_{t\to\infty}\frac{\ln|D(k,t)|}{t}$ (13)

takes a finite positive value. The Lyapunov exponent is a self-averaged
quantity, independent of the particular realization of disorder. It is also
related to the larger eigenvalue $\nu(k,N)$ of the transfer matrix
$M(N)=M_{N}*M_{N-1}...M_{1}$ corresponding to a chain of $N$ random kicks.
Actually $\lambda(k)=\lim_{N\to\infty}\ln\nu(k,N)/N$ Comtet2013 , which in
terms of the energy $U_{N}$ reads as $\lambda(k)=\lim_{N\to\infty}\ln
U_{N}/2N$. Thus, for large $j$, the variable $z_{j}$ is exponentially large
and we can assume

$2\sqrt{z_{j-1}(1+z_{j-1})}/(1+2z_{j-1})\simeq 1\,,$ (14)

which allows us to simplify Eq. (12) into

$\ln U_{j}=\ln
U_{j-1}+\ln(1+2\beta_{j})+\ln\left(1+\gamma_{j}\,\cos\Theta_{j}\right)\,.$
(15)

We now introduce the important assumption that the cumulative phase
$\Theta_{j}$ is completely random, with a uniform distribution over
$[0,2\pi]$. This assumption defines the so-called “random phase” model. It
cannot be generally valid for modulations like those in Fig. 1, with an
arbitrary degree of disorder. However, in many cases, and for a sufficiently
long chain of kicks, the cumulative phase does get randomized and its
distribution becomes close to uniform (this can be checked numerically, as
discussed in the last section). From Eq. (15) we can write

$\ln U_{N}=\sum_{j=1}^{N}X_{j}\,,$ (16)

and consider that the terms
$X_{j}=\ln(1+2\beta_{j})+\ln\left(1+\gamma_{j}\,\cos\Theta_{j}\right)$ are
independent and identically distributed random variables note_sum . In the
large $N$ limit, according to the central-limit theorem, this implies that
that $\ln U_{N}$ has a Gaussian distribution. It only remains to calculate the
average and the variance of that distribution. Averaging $X_{j}$ over the
random phase $\Theta_{j}$ note_average  followed by averaging over $\beta_{j}$
yields

$\langle\ln U_{N}\rangle=N\langle\ln(1+\beta)\rangle\,.$ (17)

Here $\langle...\rangle$ denotes the full statistical average, over the
statistical distributions of $\Theta_{j}$ and $\beta$ (we have dropped the
subscript $j$ since the statistical distribution of $\beta_{j}$ are taken to
be independent of $j$). We have found that in the large $N$ limit $\langle\ln
U_{N}\rangle$ grows linearly with $N$, with a slope
$\langle\ln(1+\beta)\rangle$. The linear increase of $\langle\ln U_{N}\rangle$
for large $N$, or equivalently of $\langle\ln U(t)\rangle$ at long times, is
reminiscent of known features of Anderson localization in a random spatially-
modulated medium. An analogy can be drawn from two standpoints. The Lyapunov
exponent defined in Eq. (13) also appears in the one-dimensional Anderson
localization problem Comtet2013 . In that problem one is interested in the
solution to the stationary Schrödinger equation, which is different from the
(time) Cauchy problem stated in Eq. (2). It can be useful, however, to
consider a Cauchy problem for Anderson localization, by fixing the wave
function and its spatial derivative at some point in space and then
calculating the Lyapunov exponent, that turns out to equal the inverse
localization length. We stress that in localization theory, this approach is
used as a trick to calculate the localization length, while in the present
work the Cauchy problem appears naturally due to physical initial conditions.
Another analogy, both physical and mathematical, exists between the resistance
of a one-dimensional spatially-disordered conductor, and the fraction of
backscattered energy $z_{N}=|D_{N}^{-}|^{2}$ in our problem. More precisely,
the resistance $\rho(L)$ of chain with length $L$ is known to grow
exponentially with $L$, which is another manifestation of Anderson
localization Anderson1980 . The backscattered wave energy $z_{N}$ grows with
$N$ in a similar way.

Calculating the variance of $\ln U_{N}$ is also possible but one ends up with
an integral that cannot be evaluated analytically, unless the disorder is
weak. The case of weak disorder is of special interest because it is relevant
to experiments (indeed, the modulation of $\varepsilon$ is expected to be very
small), and it is amenable to complete analytical treatment.

Weak disorder. We define the weak-disorder regime by the condition
$\beta_{j}\ll 1$. In this case $\langle\ln U_{N}\rangle=N\langle\beta\rangle$,
which follows from (17). This relation shows that $\langle\beta\rangle/2$ is
the Lyapunov exponent for weak disorder. To determine the variance, we note
that $X_{j}\simeq
2\beta_{j}+2\sqrt{\beta_{j}}\cos\Theta_{j}-2\beta_{j}\cos^{2}\Theta_{j}$. To
leading order in $\beta_{j}$, averaging over the random phase leads to
$\langle X_{j}^{2}\rangle_{\Theta}\simeq
4\beta_{j}\langle\cos^{2}\Theta_{j}\rangle_{\Theta}=2\beta_{j}$ from which we
deduce

$\mathrm{Var}(\ln U_{N})=2N\langle\beta\rangle=2\langle\ln U_{N}\rangle\,.$
(18)

We conclude that for large $N$ the wave energy $U_{N}$ has a log-normal
distribution, with mean value and variance satisfying (18). An identical
result is known in Anderson localization along weakly disordered chains
Shapiro1988 , with the resistance being the analog of $U_{N}$. Our result
implies a high degree of universality in one-dimensional wave transport. Not
only a universal log-normal distribution is approached for large $N$, but the
variance and mean are related by a factor of two, which is a signature of
single parameter scaling.

The meaning of the large $N$ limit can be clarified. The above treatment is
based on Eq. (15) which, in turn, is based on the assumption $z_{j}\gg 1$.
Initially $z_{j}$ is very small and gradually grows to reach a value
$z_{j}\sim 1$ after a large number of kicks on the order of
$N_{c}=1/\langle\beta\rangle$. Thus, the condition for the log-normal
distribution and relation (18) is $N\gg N_{c}$.

It is also interesting to characterize the intermediate regime $1\ll N\ll
N_{c}$, in which we can also expect some universal - albeit different -
statistical distribution for the wave energy $U_{N}$. Going back to the
general Eq.(12), which has no restriction on the value of $z_{j}$, we set
there $z_{j},z_{j-1}\ll 1$ as well as $\beta_{j}\ll 1$. This leads to the
recursion relation

$z_{j}=z_{j-1}+\beta_{j}+2\sqrt{\beta_{j}z_{j-1}}\cos\Theta_{j}\ \ (j\ll
N_{c})\,.$ (19)

As before, we assume that the phase $\Theta_{j}$ is completely random or, at
least, that it gets randomized after some number of kicks $j_{0}\ll N_{c}$.
Next, we raise Eq. (19) to power $n$ and average first over $\Theta_{j}$, and
then over some arbitrary distribution of $\beta_{j}$, keeping only the leading
(linear) terms in $\beta_{j}$. This enables us to express the $n^{th}$ moment
of $z_{j}$ in the form (see Supplemental Material SM ) $\langle
z_{j}^{n}\rangle=(n!)\langle z_{j}\rangle^{n}$ for $j_{0}\ll j\ll N_{c}$, with
$\langle z_{j}\rangle=j\langle\beta\rangle$. This implies that after a
sequence of $N$ kicks, $z_{N}$ follows a negative exponential (or Rayleigh)
distribution:

$P(z_{N})=(N\langle\beta\rangle)^{-1}\exp[-z_{N}/(N\langle\beta\rangle)]\ \
(1\ll N\ll N_{c})\,.$ (20)

We conclude that the energy $U_{N}=1+2z_{N}$ in this regime has negative
exponential distribution, and that the average energy $\langle
U_{N}\rangle=1+2N\langle\beta\rangle$ grows linearly with the number of kicks.

It is possible to treat both regimes of short and long chains, or equivalently
short and long times, using a more formal approach based on a variant of a
Fokker-Planck equation, sometimes referred to as Melnikov’s equation
Melnikov_etal . To proceed, we start with the basic recursion relation (12)
for the variable $z$ and transform it into a recursion relation for the
distribution $P_{j}(z)$ for that variable after $j$ kicks (the derivation,
given in Ref. Melnikov_etal , is summarized in the Supplemental Materiel SM ).
In the weak-disorder regime $\langle\beta\rangle\ll 1$, this equation is

$P_{j}(z)=P_{j-1}(z)+\langle\beta\rangle\,\frac{\partial}{\partial
z}\left[(z+z^{2})\frac{\partial P_{j-1}}{\partial z}(z)\right]\,.$ (21)

Next, we transform the discrete time steps $t_{j}$ into the continuous time
$t$, by using the average time interval $\Delta t=\langle
t_{j}-t_{j-1}\rangle$. This leads to

$\frac{\partial P_{t}}{\partial t}(z)=\alpha\,\frac{\partial}{\partial
z}\left[(z+z^{2})\frac{\partial P_{t}}{\partial z}(z)\right]\,,$ (22)

with $\alpha=\langle\beta\rangle/\Delta t$. In principle, Eq. (22) should be
solved with an initial condition $P_{t=0}(z)$. Actually, the precise shape of
the initial distribution is rapidly forgotten and a universal function of $z$
(with $\alpha$ as a single parameter) is approached as time elapses. There are
two distinct regimes, which can be separated using the critical time
$\tau_{c}=1/\alpha$, which is the counterpart of $N_{c}$ in the continuous
time picture. At short times $t\ll\tau_{c}$, $z$ remains small and we can
neglect the $z^{2}$ term in (22) to obtain

$P_{t}(z)=(\alpha t)^{-1}\exp[-z/(\alpha t)]\ \ (t\ll\tau_{c})\,.$ (23)

We find that $z$ follows a negative exponential distribution, identical to Eq.
(20), but for continuous time. In the opposite limit $t\gg\tau_{c}$, the $z$
term in (22) can be neglected, and a log-normal distribution for $z$ is
obtained:

$P_{t}(z)=(z\sqrt{4\pi\alpha t})^{-1}\exp[-(\ln z-\alpha t)^{2}/(4\alpha t)]\
\ (t\gg\tau_{c})\,.$ (24)

This long-time statistics is in agreement with that obtained previously for a
discrete chain of kicks in the limit $N\gg N_{c}$.

Numerical results. In order to support and illustrate the theoretical
analyses, we have carried out numerical simulations, using the $\delta$-kicks
model defined in Fig. 1. The transfer matrix in this case takes the form (6),
with coefficients given by Eq. (7). In the simulations $u_{j}$ and
$\theta_{j}$ are taken to be uniformly distributed random variables, with
$u_{j}\in[0,0.05]$ (corresponding to weak disorder) and
$\theta_{j}\in[0,2\pi]$. By performing products of transfer matrices, we
simulate a random chain of $N$ kicks, and calculate numerically the energy
$U_{N}$. Doing this for many realizations of disorder (i.e. of $u_{j}$ and
$\theta_{j}$), we can compute the statistical distributions of $U_{N}$ or $\ln
U_{N}$, and compare the numerical results with the theoretical predictions.

Focusing first on the average energy, we show in Fig. 2 a plot of $\langle\ln
U_{N}\rangle$ versus the number of kicks $N$. For $N\gg N_{c}$, with
$N_{c}=1/\langle\beta\rangle\simeq 1200$ here, we find that $\langle\ln
U_{N}\rangle$ grows linearly with $N$, with a slope $\langle\beta\rangle\simeq
8.3\times 10^{-4}$ coinciding with that predicted theoretically, as indicated
by the straight line. In the region $1\ll N\ll N_{c}$, we observe a regime in
which $\langle U_{N}\rangle\sim 1+2N\langle\beta\rangle$, also predicted
theoretically. An enlargement of this intermediate regime is shown in the
inset. Although not shown for brevity, we have checked numerically that the
condition of the random phase model is satisfied as soon as $N\gg 1$, and that
the condition (14) is satisfied for $N\gg N_{c}$.

Figure 2:  Numerical simuation of $\langle\ln U_{N}\rangle$ versus the number
of kicks $N$ (black line). For $N\gg N_{c}$ ($N_{c}\simeq 1200$ here) we
observe a linear growth with a slope $\langle\beta\rangle$ (indicated by the
red line). For $1\ll N\ll N_{c}$, we observe a regime where $\langle
U_{N}\rangle\sim 1+2N\langle\beta\rangle$, as shown in the inset. The
numerical simulations confirm the theoretical predictions. A $\delta$-kicks
model as in Fig. 1 is used in the simulations, with $u_{j}$ and $\theta_{j}$
uniformly distributed in $[0,0.05]$ and $[0,2\pi]$.

The numerical simulation also permits a computation of the full statistical
distribution of $U_{N}$. The distributions in the intermediate regime $1\ll
N\ll N_{c}$ and in the large $N$ limit $N\gg N_{c}$ are shown in Figs. 3 (a)
and (b). In the intermediate regime, we find that the energy $U_{N}$ follows a
negative exponential law, with average value $\langle
U_{N}\rangle=1+2N\langle\beta\rangle$. For $N\gg N_{c}$, we find that the
distribution of $U_{N}$ is log-normal ($\ln U_{N}$ is Gaussian), with
$\langle\ln U_{N}\rangle=N\langle\beta\rangle$ and $\mathrm{Var}(\ln
U_{N})=2\langle\ln U_{N}\rangle$. The calculated statistical distributions
perfectly match the theoretical predictions, and confirm the universal
character and single-parameter scaling of wave transport in randomly time-
varying homogeneous and isotropic media, in the regime of weak disorder.

Figure 3:  (a): Statistical distribution of $U_{N}$ in the intermediate regime
$1\ll N\ll N_{C}$ ($N=50$). A negative exponential distribution is observed,
in full agreement with theory. (b): Statistical distribution of $\ln U_{N}$ in
the asymptotic regime $N\gg N_{c}$ ($N=2\times 10^{4}$). A Gaussian
distribution is observed, also in agreement with theory. Same model as in Fig.
2.

Conclusion. In summary, we have presented a general model for wave propagation
in a random time-dependent medium, and demonstrated the existence of universal
statistical distributions of the wave energy $U$. We proved that, after a
sufficiently long time, $U$ approaches a log-normal distribution with
$\langle\ln U\rangle\sim t$, in agreement with well-established results in
one-dimensional transport. In the weak-modulation regime, a full analytic
theory was developed, which reveals two distinct regimes: For time smaller
than some crossover time $\tau_{c}$, the energy distribution follows a
negative exponential (Rayleigh) distribution, while for times $t\gg\tau_{c}$
the distribution crosses over to a log-normal law. The intermediate regime for
$t\ll\tau_{c}$ might be relevant to experiments in which the long-time regime
could be difficult to reach. The theory, in perfect agreement with numerical
simulations, lays some foundation in the emerging topic of waves in disordered
temporal media, with expected outcomes in the control of various kinds of
waves.

###### Acknowledgements.

We thank A. Piette for his help in the early stage of the project, and E. Fort
and M. Fink for stimulating discussions. Work supported by LABEX WIFI
(Laboratory of Excellence within the French Program Investments for the
Future) under references ANR-10- LABX-24 and ANR-10-IDEX-0001-02 PSL*.

## References

  * (1) C. Caloz and Z.L. Deck-Léger, IEEE Trans. Antennas Propag. 68, 1569 (2020); 1583 (2020).
  * (2) F. Zangeneh-Nejad and R. Fleury, Rev. Phys. 4, 100031 (2019).
  * (3) F. R. Morgenthaler, IRE Trans. Microwave Theory and Techniques 6, 167 (1958).
  * (4) V. Bacot, M. Labousse, A. Eddi, M. Fink, and E. Fort, Nat. Phys. 12, 972 (2016).
  * (5) V. Bacot, G. Durey, A. Eddi, M. Fink, and E. Fort, Proc. Natl. Acad. Sci. USA 116, 8809 (2019).
  * (6) D. Holberg and K. Kunz, IEEE Trans. Antennas Propag. 14, 183 (1966).
  * (7) J. R. Zurita-Sanchez, P. Halevi, and J. C. Cervantes-González, Phys. Rev. A 79, 053821 (2009).
  * (8) E. Lustig, Y. Sharabi, and M. Segev, Optica 5, 1390 (2018).
  * (9) A. Akbarzadeh, N. Chamanara, and C. Caloz, Opt. Lett. 43, 3297 (2018).
  * (10) V. Pacheco-Pena and N. Engheta, Optica 7, 323 (2020).
  * (11) V. Pacheco-Pena and N. Engheta, Light: Sci. Appl. 9, 129 (2020).
  * (12) Y. Sharabi, E. Lustig, and M. Segev, Phys. Rev. Lett. 126, 163902 (2021).
  * (13) B. Apffel, S. Wildeman, A. Eddi, and E. Fort, arXiv:2102.07583.
  * (14) P. W. Anderson, Phys. Rev. 109, 1492 (1958).
  * (15) L.D. Landau, E.M. Lifshitz and L.P. Pitaevskii, Electrodynamics of Continuous Media, 2nd ed. (Pergamon Press, Oxford, 1984), chap. IX, section 80.
  * (16) See Supplemental Material which includes (i) the transfer matrix for a rectangular barriers model, (ii) the derivation of the transfer matrix of the $\delta$-kicks model, (iii) the derivation of the moment relation leading to Eq. (20), and (iv) the derivation of Eq. (22).
  * (17) P.A. Mello and N. Kumar, Quantum Transport in Mesoscopic Systems: Complexity and Statistical Fluctuations. A Maximum Entropy Viewpoint (Oxford University Press, 2004), section 2.1.5.
  * (18) A. Comtet, C. Texier, and Y. Tourigny, J. Phys. A: Math. Theor. 46, 254003 (2013).
  * (19) Since phase randomization is established after a certain number of kicks, the sum in Eq. (16) should start from some model-dependent value $j_{0}>1$. This, however, makes no difference in the large $N$ limit.
  * (20) This amounts to calculating $\int_{0}^{2\pi}\ln(1+\gamma_{j}\cos\Theta)\,d\Theta=2\pi\ln[(1+\sqrt{1-\gamma_{j}^{2}})/2]$ for $\gamma_{j}<1$.
  * (21) P.W. Anderson, D.J. Thouless, E. Abrahams and D.S. Fisher, Phys. Rev. B 22, 3519 (1980).
  * (22) A. Cohen, Y. Roth, and B. Shapiro, Phys. Rev. B 38, 12125 (1988).
  * (23) V. I. Melnikov, Sov. Phys. Solid State 23, 444 (1981); N. Kumar, Phys. Rev. B 31, 5513 (1985); B. Shapiro, Phys. Rev. B 34, 4394 (1986); P. A. Mello, J. Math. Phys. (N.Y.) 27, 2876 (1986).

## Supplemental Material

##  I Random barriers model

In this section we describe a model for a time-disordered homogeneous medium
that could be used as an alternative to the $\delta$-kicks model presented in
the main text. $\Omega^{2}(t)$ is considered to be a chain of rectangular
barriers, as represented in Fig. 4. Each barrier is denoted as a kick, with
the times $t^{\prime}_{j}$ and $t_{j}$ defining the onset and end of kick
number $j$. The times $t^{\prime}_{j}$ and $t_{j}$, as well as the kick
strengths (barrier heights) $\Omega^{2}_{j}$ are random variables.

Figure 4:  Random chain of rectangular barriers, each of them referred to as a
“kick”. The barrier heights $\Omega^{2}_{j}$ and the times $t^{\prime}_{j}$
and $t_{j}$ are random variables. After each kick the medium recovers the
background value $\Omega^{2}_{b}$.

The transfer matrix ${\bf M}_{j}$ of a single block $j$ takes the form

${\bf M}_{j}=\begin{bmatrix}a_{j}&b_{j}\\\ b_{j}^{*}&a_{j}^{*}\end{bmatrix}$
(25)

with

$\displaystyle a_{j}$ $\displaystyle=$ $\displaystyle
T_{j}T^{\prime}_{j}\exp[-i(\theta_{j}+\phi_{j})]+R_{j}R^{\prime}_{j}\exp[-i(\theta_{j}-\phi_{j})]$
(26) $\displaystyle b_{j}$ $\displaystyle=$ $\displaystyle
T_{j}R^{\prime}_{j}\exp[i(\theta_{j}-\phi_{j})]+R_{j}T^{\prime}_{j}\exp[i(\theta_{j}+\phi_{j})]\,.$
(27)

Here, $T_{j}=(\Omega_{b}+\Omega_{j})/(2\Omega_{b})$,
$T^{\prime}_{j}=(\Omega_{b}+\Omega_{j})/(2\Omega_{j})$,
$R_{j}=(\Omega_{b}-\Omega_{j})/(2\Omega_{b})$ and
$R^{\prime}_{j}=(\Omega_{j}-\Omega_{b})/(2\Omega_{j})$ are time-domain
transmission and reflection coefficients, that can be deduced from the
continuity of the field and its time derivative at the times
$t=t^{\prime}_{j}$ and $t=t_{j}$. The phases
$\phi_{j}=\Omega_{j}(t_{j}-t^{\prime}_{j})$ and
$\theta_{j}=\Omega_{b}(t^{\prime}_{j}-t_{j-1})$ correspond to free propagation
during the kick, and free propagation between two succesive kicks,
respectively.

We note that the transfer matrix takes the same form as Eq. (6) in the main
text, but with different matrix elements. It can be verified that the matrix
elements also satisfy Eqs. (8) and (9) of the main text, which are very
general and independent of the model of disorder (see Ref. MelloBook  for a
general derivation of these properties).

##  II Transfer matrix for $\delta$-kicks

In this section we deduce the expression of the transfer matrix for the
$\delta$-kicks model described in Fig. 1 of the main text, starting from the
random barriers model decribed in the previous section. For kick number $j$ we
take a barrier with duration $\delta t$ and amplitude such that

$t_{j}-t^{\prime}_{j}=\delta t,\quad\Omega_{j}^{2}-\Omega_{b}^{2}=v_{j}/\delta
t.$ (28)

Taking the limit $\delta t\to 0$ leads to
$\Omega^{2}(t)=\Omega_{b}^{2}+\sum_{j}v_{j}\delta(t-t_{j})$ which corresponds
to the $\delta$-kicks model with $v_{j}$ the amplitude of kick number $j$.
Taking the same limit for the coefficients of the transfer matrix ${\bf
M}_{j}$ in Eqs. (26) and (27) leads to

$\displaystyle a_{j}$ $\displaystyle=$
$\displaystyle[1-iv_{j}/(2\Omega_{b})]\exp(-i\theta_{j}),$ (29) $\displaystyle
b_{j}$ $\displaystyle=$ $\displaystyle-iv_{j}/(2\Omega_{b})\exp(i\theta_{j}).$
(30)

These expressions correspond to Eq. (7) of the main text.

##  III Derivation of the moment relation leading to Eq. (20)

In this section we derive the moment relation $\langle
z_{j}^{n}\rangle=(n!)\langle z_{j}\rangle^{n}$ that is used as a step in the
derivation of Eq. (20) in the main text. We start by raising Eq. (19) of the
main text to the power $n$. Keeping terms of order $\sqrt{\beta_{j}}$ and
$\beta_{j}$, we obtain

$z_{j}^{n}=z_{j-1}^{n}+n\beta_{j}z_{j-1}^{n-1}+2nz_{j-1}^{n-1}\sqrt{\beta_{j}z_{j-1}}\cos(\Theta_{j})\\\
+2n(n-1)z_{j-1}^{n-1}\beta_{j}\cos(\Theta_{j})^{2}+\operatorname{O}(\beta_{j}^{3/2}).$
(31)

We perform an average over $\Theta$, assuming that it is fully randomized
after a sufficiently large number of kicks (this hypothesis is discussed in
the main text) and independent of $\beta_{j}$, and a subsequent average over
an arbitrary distribution for $\beta_{j}$. This leads to the following
recursion relation for the moments of $z$

$\langle z_{j}^{n}\rangle=\langle
z_{j-1}^{n}\rangle+n^{2}\langle\beta\rangle\langle z_{j-1}^{n-1}\rangle\,,$
(32)

with the first moment given by $\langle z_{j}\rangle=j\langle\beta\rangle$. We
now assume that the moment of order $n$ is

$\langle z_{j}^{n}\rangle=n!\langle
z_{j}\rangle^{n}=n!j^{n}\langle\beta\rangle^{n}.$ (33)

Inserting the above relation in Eq. (32), it is easy to see that it is
satisfied up to terms of order $1/j^{2}$. This concludes the derivation of the
relation $\langle z_{j}^{n}\rangle=n!\langle z_{j}\rangle^{n}$, valid for
large enough $j$.

##  IV Derivation of Melnikov’s equation

In this section we derive Eq. (21) of the main text. The derivation can be
found in Ref. [S1], and we summarize the main steps here. We starts with the
basic recursion relation for the variable $z$ [Eq. (12) in the main text],
that can be rewritten as

$z_{j}=z_{j-1}+\beta_{j}(1+2z_{j-1})\\\
+2\sqrt{\beta_{j}(1+\beta_{j})z_{j-1}(1+z_{j-1})}\cos\Theta_{j}.$ (34)

Considering as above that $\Theta_{j}$ is fully randomized, the above relation
can be transformed into a recursion relation for the probability density
$P_{j}(z)$, which takes the form

$P_{j}(z)=\int_{0}^{\infty}f(\beta)d\beta\int_{0}^{2\pi}\frac{d\Theta}{2\pi}\\\
\times
P_{j-1}\left(z+\beta(1+2z)-2\sqrt{\beta(1+\beta)z(1+z)}\cos\Theta\right)\,,$
(35)

where $f(\beta)$ is the probability density asociated to the random variable
$\beta$. For weak disorder such that $\beta\ll 1$, we can perform a second
order Taylor expansion for $P_{j-1}$, which leads to

$P_{j-1}(z+\epsilon)=P_{j-1}(z)+\epsilon\frac{\partial P_{j-1}(z)}{\partial
z}+\frac{\epsilon^{2}}{2}\frac{\partial^{2}P_{j-1}(z)}{\partial z^{2}}\\\
+\operatorname{O}(\epsilon^{3})\,,$ (36)

where $\epsilon=\beta(1+2z)-2\sqrt{\beta(1+\beta)z(1+z)}\cos\Theta$.
Substituting in Eq. (35), we obtain

$P_{j}(z)=P_{j-1}(z)+\langle\beta\rangle(1+2z)\frac{\partial
P_{j-1}(z)}{\partial z}\\\ +\langle\beta\rangle
z(1+z)\frac{\partial^{2}P_{j-1}(z)}{\partial
z^{2}}+\operatorname{O}(\beta^{2})\,.$ (37)

This can be factorized in the form

$P_{j}(z)=P_{j-1}(z)+\langle\beta\rangle\frac{\partial}{\partial
z}\left[(z+z^{2})\frac{\partial P_{j-1}}{\partial z}(z)\right]\,,$ (38)

where terms beyond first order have been neglected. This is Eq. (21) in the
main text, which in the continuous limit leads to Melnikov’s equation (22).

### References

[S1] P.A. Mello and N. Kumar, Quantum Transport in Mesoscopic Systems:
Complexity and Statistical Fluctuations. A Maximum Entropy Viewpoint (Oxford
University Press, 2004).


$\rho=\rho(n,\lambda,\|A\|_{C^{1}(\overline{B}_{1})},\|b\|_{C^{0}(\overline{B}_{1})})$.
The functions $u(y_{j}+\rho\cdot)$ are stable solutions of an elliptic
equation with coefficients $A^{y_{j},\rho}=A(y_{j}+\rho\cdot)$ and
$b^{y_{j},\rho}=\rho b(y_{j}+\rho\cdot)$. Since $B_{\rho}(y_{j})\subset B_{1}$
and by (7.23), the coefficients satisfy the bounds

$\begin{split}\|DA^{y_{j},\rho}\|_{C^{0}(\overline{B}_{1})}+\|b^{y_{j},\rho}\|_{C^{0}(\overline{B}_{1})}&\leq\left(\|DA\|_{C^{0}(\overline{B}_{\rho}(y_{j}))}+\|b\|_{C^{0}(\overline{B}_{\rho}(y_{j}))}\right)\rho\leq\varepsilon_{0},\end{split}$

therefore, we can apply Step 2 to deduce

$\|u\|_{C^{\alpha}(\overline{B}_{1/2})}\leq\sum_{j}\|u\|_{C^{\alpha}(\overline{B}_{\theta\rho}(y_{j}))}\leq
C\sum_{j}\|u\|_{L^{1}(B_{\rho}(y_{j}))}\leq C\|u\|_{L^{1}(B_{1})},$

where
$C=C(n,\lambda,\|A\|_{C^{1}(\overline{B}_{1})},\|b\|_{C^{0}(\overline{B}_{1})})$.
This concludes the proof of the theorem. ∎

##  Appendix A Stability is not equivalent to the integral inequality

Let $u\in C^{2}(\overline{\Omega})$ be a solution to $-Lu=f(u)$ in $\Omega$
with $u=0$ on $\partial\Omega$. Recall that $u$ is a stable solution if

$J_{u}\varphi=L\varphi+f^{\prime}(u)\varphi\leq 0\quad\text{ in }\Omega,$
(A.1)

for some function $\varphi\in C^{2}(\overline{\Omega})$ with $\varphi>0$ in
$\Omega$ and $\varphi=0$ on $\partial\Omega$. This is the stability condition
(1.3) presented in the Introduction and is equivalent to the nonnegativity of
the first Dirichlet eigenvalue of $J_{u}$ (with the sign convention
$J_{u}\varphi=-\mu\varphi$). There, we also showed that stable solutions
satisfy the integral inequality (1.10), which reads

$\int_{\Omega}f^{\prime}(u)\xi^{2}\mathop{}\\!\mathrm{d}x\leq\int_{\Omega}|\nabla\xi-\textstyle\frac{1}{2}A^{-1}(x)\widehat{b}(x)\xi|^{2}_{A(x)}\quad\text{
for all }\xi\in C^{\infty}_{c}(\Omega).$ (A.2)

Our goal in this appendix is to show that the integral inequality (A.2) does
not imply the stability condition (A.1) in general. The main reason is that
the problem is not variational, due to the drift in $L$. We also give
conditions under which the equivalence holds. Namely, writing the operator in
divergence form $Lu={\rm div}(A(x)\nabla u)+\widehat{b}(x)\cdot\nabla u$, we
show that if $A^{-1}(x)\widehat{b}(x)$ is the gradient of a scalar function,
then the problem is variational and the two conditions are equivalent.

First we write the integrals in (A.2) as the quadratic form associated to a
linear self-adjoint operator. Integrating by parts, we have

$\int_{\Omega}\left(|\nabla\xi-\textstyle\frac{1}{2}\xi
A^{-1}(x)\widehat{b}(x)|_{A(x)}^{2}-f^{\prime}(u)\xi^{2}\right)\mathop{}\\!\mathrm{d}x=-\displaystyle\int_{\Omega}\xi\widetilde{J_{u}}\xi\mathop{}\\!\mathrm{d}x,$
(A.3)

where $\widetilde{J_{u}}$ is the operator

$\widetilde{J_{u}}\xi:={\rm
div}(A(x)\nabla\xi)-\Big{\\{}\textstyle\frac{1}{2}{\rm
div}(\widehat{b}(x))+\textstyle\frac{1}{4}|\widehat{b}(x)|_{A^{-1}(x)}^{2}\Big{\\}}\xi+f^{\prime}(u)\xi.$
(A.4)

Hence, by the variational characterization of eigenvalues, (A.2) amounts to
the nonnegativity of the principal eigenvalue of $\widetilde{J_{u}}$.

We can now state our example of a solution satisfying (A.2) but not (A.1):

######  Example 1.

Consider the operator $Lv=\Delta v+\widehat{b}(x)\cdot\nabla v$ with vector
field

$\widehat{b}(x)=\frac{-x_{2}e_{1}+x_{1}e_{2}}{\sqrt{x_{1}^{2}+x_{2}^{2}}}.$

For each constant $c>0$, we let $f(u)=(\lambda_{1}+c)u+1$, where $\lambda_{1}$
denotes the least Dirichlet eigenvalue of the Laplacian in the unit ball
$B_{1}$.

If $c>0$ is sufficiently small, then the unique solution $u$ to the boundary
value problem

$\left\\{\begin{array}[]{ll}-Lu=f(u)&\text{ in }B_{1}\\\ u=0&\text{ on
}\partial B_{1}\end{array}\right.$

satisfies the integral stability condition (A.2) but is not a stable solution,
i.e., the stability condition (A.1) does not hold.555The function $u$ can be
given explicitly in terms of Bessel functions of the first kind
$\mathcal{J}_{\alpha}$ as
$u(x)=\frac{1}{(\lambda_{1}+c)\mathcal{J}_{\frac{n-2}{2}}\left(\sqrt{\lambda_{1}+c}\right)}|x|^{\frac{2-n}{2}}\mathcal{J}_{\frac{n-2}{2}}\left(\sqrt{\lambda_{1}+c}\,|x|\right)-\frac{1}{\lambda_{1}+c}.$

###### Proof.

The problem for $u$ is equivalent to solving

$\left\\{\begin{array}[]{ll}-\Delta u-\widehat{b}(x)\cdot\nabla
u-(\lambda_{1}+c)u=1&\text{ in }B_{1}\\\ u=0&\text{ on }\partial
B_{1}.\end{array}\right.$ (A.5)

Notice that the drift $\widehat{b}\in L^{\infty}(B_{1})$ has a weak derivative
$D\widehat{b}\in L^{p}(B_{1})$ for $1\leq p<2$, and satisfies the identities
$|\widehat{b}(x)|=1$ and ${\rm div}\,\widehat{b}(x)=0$ for a.e. $x\in B_{1}$.
Moreover, since $\widehat{b}$ is tangent to spheres, the derivative
$\widehat{b}(x)\cdot\nabla$ vanishes on radial functions. In particular, the
principal eigenfunction of the Laplacian is also an eigenfunction of the
adjoint operator $L^{T}=\Delta-{\rm
div}\big{(}\widehat{b}(x)\cdot\big{)}=\Delta-\widehat{b}(x)\cdot\nabla$, with
eigenvalue $\lambda_{1}$. Since the point spectrum of $L^{T}$ is discrete, for
$c>0$ small, we deduce that $\lambda_{1}+c$ is not an eigenvalue of the
adjoint operator. The Fredholm alternative now gives that (A.5) has a unique
solution.

Let $\varphi_{1}$ and $\xi_{1}$ be positive principal eigenfunctions of
$J_{u}$ and $\widetilde{J_{u}}$, respectively. Since $\varphi_{1}$ and
$\xi_{1}$ are positive in $B_{1}$, they must be radial. It follows that

$J_{u}\varphi_{1}=\Delta\varphi_{1}+(\lambda_{1}+c)\varphi_{1}=-\mu_{1}\varphi_{1}\quad\text{
and
}\quad\widetilde{J_{u}}\xi_{1}=\Delta\xi_{1}+(\lambda_{1}+c-1/4)\xi_{1}=-\widetilde{\mu_{1}}\xi_{1},$

where $\mu_{1}$ and $\widetilde{\mu_{1}}$ are the least eigenvalues of each
operator. By uniqueness, the functions are multiples of the principal
eigenfunction of the Laplacian. Therefore, we have $\mu_{1}=-c<0$ and
$\widetilde{\mu_{1}}=1/4-c>0$ for $c$ sufficiently small. This means that $u$
is not stable but (A.2) holds, which was the claim. ∎

Next we investigate the relation between the failure of the equivalence and
the form of the drift $\widehat{b}$. Let $\varphi_{1}\in
C^{2}(\overline{\Omega})$ be the unique positive principal eigenfunction of
$J_{u}$ with $\int\varphi_{1}^{2}\mathop{}\\!\mathrm{d}x=1$. In particular,
the function satisfies $\varphi_{1}>0$ in $\Omega$, $\varphi_{1}=0$ on
$\partial\Omega$, and $J_{u}\varphi_{1}=-\mu_{1}\varphi_{1}$, where
$\mu_{1}\in\mathbb{R}$ is the least eigenvalue of $J_{u}$. Consider a test
function $\xi\in C^{\infty}_{c}(\Omega)$. Multiplying $J_{u}\varphi_{1}$ by
$\xi^{2}/\varphi_{1}$ and integrating by parts in $\Omega$, we have

$\begin{split}-\mu_{1}&=\int_{\Omega}(J_{u}\varphi_{1})\frac{\xi^{2}}{\varphi_{1}}\mathop{}\\!\mathrm{d}x=\int_{\Omega}\left(-A(x)\nabla\varphi_{1}\cdot\nabla\left(\frac{\xi^{2}}{\varphi_{1}}\right)+\widehat{b}(x)\cdot\frac{\xi^{2}}{\varphi_{1}}\nabla\varphi_{1}+f^{\prime}(u)\xi^{2}\right)\mathop{}\\!\mathrm{d}x\\\
&=\int_{\Omega}\left(|\xi\nabla\log{\varphi_{1}}|_{A(x)}^{2}-2A(x)\xi\,\nabla\log{\varphi_{1}}\cdot\nabla\xi+\xi\,\widehat{b}(x)\cdot\xi\,\nabla\log{\varphi_{1}}+f^{\prime}(u)\xi^{2}\right)\mathop{}\\!\mathrm{d}x.\end{split}$

Using that

$\begin{split}&|\xi\nabla\log{\varphi_{1}}|_{A(x)}^{2}-2A(x)\xi\,\nabla\log{\varphi_{1}}\cdot\left(\nabla\xi-{\textstyle\frac{1}{2}}\xi
A^{-1}(x)\widehat{b}(x)\right)\\\
&\quad\quad=\left|\nabla\xi-\textstyle\frac{1}{2}\xi
A^{-1}(x)\widehat{b}(x)-\xi\nabla\log{\varphi_{1}}\right|^{2}_{A(x)}-\left|\nabla\xi-{\textstyle\frac{1}{2}}\xi
A^{-1}(x)\widehat{b}(x)\right|^{2}_{A(x)},\end{split}$

in the integral above, by (A.3) we obtain the identity

$\begin{split}-\mu_{1}=\int_{\Omega}\xi\widetilde{J_{u}}\xi\mathop{}\\!\mathrm{d}x+\int_{\Omega}\left|\nabla\xi-\textstyle\frac{1}{2}\xi
A^{-1}(x)\widehat{b}(x)-\xi\nabla\log{\varphi_{1}}\right|^{2}_{A(x)}\mathop{}\\!\mathrm{d}x.\end{split}$
(A.6)

Now, assuming the integral stability inequality (A.2), we can minimize (A.3)
among smooth functions $\xi$ with $\xi=0$ on $\partial\Omega$ and
$\int_{\Omega}\xi^{2}\mathop{}\\!\mathrm{d}x=1$. The unique positive minimizer
$\xi_{1}$ satisfies $\widetilde{J_{u}}\xi_{1}=-\widetilde{\mu_{1}}\xi_{1}$,
where $\widetilde{\mu_{1}}\geq 0$ is the least eigenvalue of
$\widetilde{J_{u}}$. Letting $\xi=\xi_{1}$ in (A.6) yields

$\begin{split}-\mu_{1}=-\widetilde{\mu_{1}}+\int_{\Omega}\big{|}\nabla\big{(}\log{\xi_{1}}-\log{\varphi_{1}}\big{)}-\textstyle\frac{1}{2}A^{-1}(x)\widehat{b}(x)\big{|}^{2}_{A(x)}\xi_{1}^{2}\mathop{}\\!\mathrm{d}x,\end{split}$
(A.7)

and from (A.7) we see that we always have $\mu_{1}\leq\widetilde{\mu_{1}}$,
with equality if and only if

$\nabla\log\left(\frac{\xi_{1}}{\varphi_{1}}\right)=\frac{1}{2}A^{-1}(x)\widehat{b}(x).$
(A.8)

This can only happen when the drift $\widehat{b}$ is of a special form. Notice
that the vector field from Example 1 is the curl of
$\sqrt{x_{1}^{2}+x_{2}^{2}}\,e_{3}$ and so, by the Helmholtz decomposition,
cannot be written as the gradient of a function.

Conversely, assume that $\widehat{b}(x)=A(x)\nabla w(x)$ for some function
$w\in C^{2}(\overline{\Omega})$. In this case, the problem can be cast in
variational form and conditions (A.1) and (A.2) are equivalent. Indeed, the
solutions of $-Lu=f(u)$ in $\Omega$ are critical points of the functional
$\mathcal{E}(u)=\int_{\Omega}e^{w(x)}\left(\textstyle\frac{1}{2}|\nabla
u|_{A(x)}^{2}-F(u)\right)\mathop{}\\!\mathrm{d}x$, where
$F(u)=\int_{0}^{u}f(t)\mathop{}\\!\mathrm{d}t$. The integral stability
inequality (A.2) amounts to the nonnegativity of the second variation

$\frac{\mathop{}\\!\mathrm{d}^{2}}{\mathop{}\\!\mathrm{d}^{2}t}\Big{|}_{t=0}\mathcal{E}\big{(}u+t\varphi\big{)}=\int_{\Omega}e^{w(x)}\Big{(}|\nabla\varphi|^{2}_{A(x)}-f^{\prime}(u)\varphi^{2}\Big{)}\mathop{}\\!\mathrm{d}x=-\int_{\Omega}e^{w(x)}\varphi
J_{u}\varphi\mathop{}\\!\mathrm{d}x$

since, letting $\varphi=e^{-w/2}\xi$ in this expression, we have

$\frac{\mathop{}\\!\mathrm{d}^{2}}{\mathop{}\\!\mathrm{d}^{2}t}\Big{|}_{t=0}\mathcal{E}\big{(}u+te^{-w/2}\xi\big{)}=\int_{\Omega}\Big{(}|\nabla\xi-\textstyle\frac{1}{2}\xi\nabla
w(x)|^{2}_{A(x)}-f^{\prime}(u)\xi^{2}\Big{)}\mathop{}\\!\mathrm{d}x=\displaystyle-\int_{\Omega}\xi\widetilde{J_{u}}\xi\mathop{}\\!\mathrm{d}x.$

In particular, since
$-\int_{\Omega}\xi\widetilde{J_{u}}\xi\mathop{}\\!\mathrm{d}x\geq\widetilde{\mu_{1}}\|\xi\|_{L^{2}(\Omega)}^{2}$
and taking $\varphi=\varphi_{1}$ to be the principal eigenfunction of $J_{u}$
above, we have

$\mu_{1}\int_{\Omega}e^{w(x)}\varphi_{1}^{2}\mathop{}\\!\mathrm{d}x=-\int_{\Omega}(e^{w/2}\varphi_{1})\widetilde{J_{u}}(e^{w/2}\varphi_{1})\mathop{}\\!\mathrm{d}x\geq\widetilde{\mu_{1}}\int_{\Omega}e^{w(x)}\varphi_{1}^{2}\mathop{}\\!\mathrm{d}x$

and we obtain the reverse inequality $\mu_{1}\geq\widetilde{\mu_{1}}$.

##  Appendix B A trace inequality

First we prove a simple lemma to control the $L^{p}$ norm in the ball by the
$L^{p}$ norms of the trace and the gradient:

######  Lemma B.1.

For $p\geq 1$ and $u\in W^{1,p}(B_{1})$, we have

$\|u\|_{L^{p}(B_{1})}^{p}\leq 2^{p-1}\left(\|u\|_{L^{p}(\partial
B_{1})}^{p}+\|\nabla u\|_{L^{p}(B_{1})}^{p}\right).$

###### Proof.

By approximation, we may assume that $u\in C^{\infty}(\overline{B}_{1})$. For
$r\in(0,1)$ and $\sigma\in\partial B_{1}$, we have
$u(r\sigma)=u(\sigma)-\int_{r}^{1}\sigma\cdot\nabla
u(t\sigma)\mathop{}\\!\mathrm{d}t$ and hence

$\begin{split}r^{n-1}|u(r\sigma)|^{p}&\leq
2^{p-1}r^{n-1}|u(\sigma)|^{p}+2^{p-1}r^{n-1}\int_{r}^{1}|\nabla
u(t\sigma)|^{p}\mathop{}\\!\mathrm{d}t\\\ &\leq
2^{p-1}|u(\sigma)|^{p}+2^{p-1}\int_{0}^{1}t^{n-1}|\nabla
u(t\sigma)|^{p}\mathop{}\\!\mathrm{d}t.\end{split}$ (B.1)

Integrating (B.1) in $\int_{0}^{1}\mathop{}\\!\mathrm{d}r\int_{\partial
B_{1}}\mathop{}\\!\mathrm{d}\mathcal{H}^{n-1}(\sigma)$ now yields the claim. ∎

We prove a Sobolev trace inequality with best exponent:

######  Proposition B.2.

For $1<p<n$, let $p^{\star}:=\frac{n-1}{n-p}p$. Then

$\|u\|_{L^{p^{\star}}(\partial B_{1})}^{p}\leq C\left(\|u\|_{L^{p}(\partial
B_{1})}^{p}+\|\nabla u\|_{L^{p}(B_{1})}^{p}\right)$

for all $u\in W^{1,p}(B_{1})$, where $C$ is a constant depending only on $n$
and $p$.

###### Proof.

By approximation, we may assume that $u\in C^{\infty}(\overline{B}_{1})$.
Recall the standard Sobolev inequality

$\|u\|_{L^{p_{S}}(B_{1})}^{p}\leq C\big{(}\|u\|_{L^{p}(B_{1})}^{p}+\|\nabla
u\|_{L^{p}(B_{1})}^{p}\big{)},$ (B.2)

where $p_{S}:=\frac{n}{n-p}p$ is the Sobolev exponent and $C$ depends only on
$n$ and $p$.

By the divergence theorem we have

$\int_{\partial
B_{1}}|u|^{p^{\star}}\mathop{}\\!\mathrm{d}\mathcal{H}^{n-1}=\int_{B_{1}}{\rm
div}(x|u|^{p^{\star}})\mathop{}\\!\mathrm{d}x=n\int_{B_{1}}|u|^{p^{\star}}\mathop{}\\!\mathrm{d}x+p^{\star}\int_{B_{1}}|u|^{p^{\star}-2}u(x\cdot\nabla
u)\mathop{}\\!\mathrm{d}x,$

whence

$\int_{\partial
B_{1}}|u|^{p^{\star}}\mathop{}\\!\mathrm{d}\mathcal{H}^{n-1}\leq
n\int_{B_{1}}|u|^{p^{\star}}\mathop{}\\!\mathrm{d}x+p^{\star}\int_{B_{1}}|u|^{p^{\star}-1}|\nabla
u|\mathop{}\\!\mathrm{d}x.$ (B.3)

The last term in (B.3) can be bounded by the Hölder inequality as

$\int_{B_{1}}|u|^{p^{\star}-1}|\nabla
u|\mathop{}\\!\mathrm{d}x\leq\left(\int_{B_{1}}|u|^{(p^{\star}-1)\frac{p}{p-1}}\mathop{}\\!\mathrm{d}x\right)^{\frac{p-1}{p}}\|\nabla
u\|_{L^{p}(B_{1})},$

and noticing that $(p^{\star}-1)\frac{p}{p-1}=p_{S}$ we deduce

$\|u\|_{L^{p^{\star}}(\partial B_{1})}^{p^{\star}}\leq
n\|u\|_{L^{p^{\star}}(B_{1})}^{p^{\star}}+p^{\star}\|\nabla
u\|_{L^{p}(B_{1})}\|u\|_{L^{p_{S}}(B_{1})}^{p^{\star}-1}.$ (B.4)

Since $p^{\star}<p_{S}$, by Hölder we have $\|u\|_{L^{p^{\star}}(B_{1})}\leq
C\|u\|_{L^{p_{S}}(B_{1})}$, and applying the Sobolev inequality (B.2) in
(B.4), we obtain the trace Sobolev inequality

$\|u\|_{L^{p^{\star}}(\partial B_{1})}^{p}\leq
C\big{(}\|u\|_{L^{p}(B_{1})}^{p}+\|\nabla u\|_{L^{p}(B_{1})}^{p}\big{)},$
(B.5)

where $C$ depends only on $n$ and $p$. Applying Lemma B.1 in (B.5) now yields
the claim. ∎

##  Appendix C Two interpolation inequalities

We recall two interpolation inequalities in cubes by Cabré [CabreQuant]. The
first one states that the $L^{2}$ norm of the gradient can be bounded in terms
of the $L^{1}$ norm of the “Hessian times the gradient” and the $L^{2}$ norm
of the function. In the second one, the $L^{2}$ norm of the function is
controlled by the $L^{2}$ norm of the gradient and the $L^{1}$ norm of the
function.

######  Proposition C.1 (X. Cabré [CabreQuant]).

Let $Q=(0,1)^{n}\subset\mathbb{R}^{n}$ and $u\in C^{2}(\overline{Q})$.

Then, for every $\delta\in(0,1)$,

$\|\nabla u\|_{L^{2}(Q)}^{2}\leq 2n\delta\|\,|\nabla
u|\,D^{2}u\,\|_{L^{1}(Q)}+n\left(\frac{18}{\delta}\right)^{2}\|u\|_{L^{2}(Q)}^{2}.$

######  Proposition C.2 (X. Cabré [CabreQuant]).

Let $Q=(0,1)^{n}\subset\mathbb{R}^{n}$ and $u\in C^{2}(\overline{Q})$.

Then, for every $\widetilde{\delta}\in(0,1)$,

$\|u\|_{L^{2}(Q)}^{2}\leq 2n^{2}\widetilde{\delta}^{2}\|\nabla
u\|_{L^{2}(Q)}^{2}+2^{n+1}\widetilde{\delta}^{-n}\|u\|_{L^{1}(Q)}^{2}.$

##  Appendix D Absorbing errors in larger balls

We recall the following well-known lemma of Simon [Simon], which allows to
absorb errors in larger balls when controlling quantities in smaller balls:

######  Lemma D.1 (L. Simon [Simon]).

Let $\beta\geq 0$ and $C_{0}>0$. Let $\mathcal{B}$ be the class of all open
balls $B$ contained in the unit ball $B_{1}$ of $\mathbb{R}^{n}$ and let
$\sigma\colon\mathcal{B}\rightarrow[0,+\infty)$ satisfy the following
subadditivity property:

$\sigma(B)\leq\sum_{j=1}^{N}\sigma(B^{j})\quad\mbox{ whenever
}N\in\mathbb{Z}^{+},\\{B^{j}\\}_{j=1}^{N}\subset\mathcal{B},\text{ and
}B\subset\bigcup_{j=1}^{N}B^{j}.$

It follows that there exists a constant $\delta>0$, which depends only on $n$
and $\beta$, such that if

$\rho^{\beta}\sigma\left(B_{\rho/2}(y)\right)\leq\delta\rho^{\beta}\sigma\left(B_{\rho}(y)\right)+C_{0}\quad\mbox{whenever
}B_{\rho}(y)\subset B_{1},$

then

$\sigma(B_{1/2})\leq CC_{0}$

for some constant $C$ which depends only on $n$ and $\beta$.

## Acknowledgements

The author thanks Xavier Cabré for his patient guidance and useful discussions
on the topic of this article.

## References