#  Correlation function of a random scalar field evolving with a rapidly
fluctuating Gaussian process

Jared C. Bronski<EMAIL_ADDRESS>Department of Mathematics, University
of Illinois Urbana-Champaign, Urbana, IL 61801, United States  Lingyun Ding
<EMAIL_ADDRESS>Richard M. McLaughlin<EMAIL_ADDRESS>Department of
Mathematics, University of North Carolina, Chapel Hill, NC, 27599, United
States

###### Abstract

We consider a scalar field governed by an advection-diffusion equation (or a
more general evolution equation) with rapidly fluctuating, Gaussian
distributed random coefficients. In the white noise limit, we derive the
closed evolution equation for the ensemble average of the random scalar field
by three different strategies, i.e., Feynman-Kac formula, the limit of
Ornstein-Uhlenbeck process, and evaluating the cluster expansion of the
propagator on an $n$-simplex. With the evolution equation of ensemble average,
we study the passive scalar transport problem with two different types of
flows, a random periodic flow, and a random strain flow. For periodic flows,
by utilizing the homogenization method, we show that the $N$-point correlation
function of the random scalar field satisfies an effective diffusion equation
at long times. For the strain flow, we explicit compute the mean of the random
scalar field and show that the statistics of the random scalar field have a
connection to the time integral of geometric Brownian motion. Interestingly,
all normalized moment (e.g., skewness, kurtosis) of this random scalar field
diverges at long times, meaning that the scalar becomes more and more
intermittent during its decay.

###### keywords:

Passive scalar , Scalar intermittency, Method of moment , Random shear flow ,
Turbulent transport

###### MSC:

[2010]37A25, 37H10, 37N10, 82C70, 76R50

††journal: XXX

##  1 Introduction

An important and difficult problem is determining the statistics of scalar
inherits from some given external random field. These statistics are often of
great physical interest. Some examples include the understanding of the
evolution of the probability distribution function for fluctuations in a
passive scalar diffusing in an incompressible turbulent random flow field [24,
36, 26, 25, 41, 22, 33], the prediction of escape times for a particle
diffusing in a random fluctuating potential [39, 38], diffusing diffusivity
model [11, 20, 43, 44], light propagating through random media [40] and random
water waves impinging on a step [5, 28] and numerous other examples.

The major question to be addressed can be formulated as follows: Is there a
universal probability distribution of a scalar in a random velocity field? If
this distribution exists, then what are its properties? Naturally, the
determination of these statistics requires understanding the evolution of
moments of the scalar field. However, it is in general difficult to derive
closed equations governing the evolution of moments since these problems
involve variable coefficients.

In certain situations, however, full moment closure is available. Such
situations typically involve problems whose coefficients are rapidly
fluctuating in time. For the case of a passive scalar, Kraichnan [24, 22]
considered a stationary isotropic turbulent velocity field admitting Gaussian
white noise time statistics and was able to derive a closed expression for the
time evolution of the scalar’s mean and covariance. Majda [26, 27], studying
the passive scalar problem in a rapidly fluctuating shear layer, obtained
closed equations for the evolution of the $N$-point correlation function by
manipulating a Feynman-Kac path integral representation for the scalar’s
Fourier transform in the shear direction. For the special case where the shear
layer is linear, these equations may be exactly solved to explicitly determine
the scalar statistics inherited from the random flow and to illustrate
situations leading to non-Gaussian distributions [26, 8, 9, 45].
Unfortunately, it is not at all obvious how to generalize the shear layer
example to more general flow fields.

The main purpose of this paper is to demonstrate that such moment closures
involving Gaussian white noise in shear layer example are special cases of a
more general result. We present three strategies to derive the evolution
equation of the $N$-point correlation function of the scalar field, i.e, the
moment closure equation. Different strategies are suitable for different types
of random flow and different governing equations. The first strategy is
suitable for the advection-diffusion equation with the flow multiplied by a
Gaussian white noise process. We derive the moment closure equation by
manipulating the Fokker-Planck equation or Feynman-Kac formula. The second
method uses the fact that the stationary Ornstein-Uhlenbeck (OU) process
converges to the Gaussian white noise process as the correlation time of the
OU process vanishes. With the moment closure equation for the OU case derived
by Resnick [34], the perturbation method yields the equation for the white
noise case. The methods may be generalized to study the advection-diffusion
equation with a more complicated flow that is a function of a Gaussian white
noise process. In the third method, we consider a very general time-dependent
equation with multiplicative white noise. We first derive the moment closure
equation by formally using the pseudo-differential operator and the
characteristic function of normal random variables. Next, we rigorously prove
this result by evaluating the complete Picard iteration expansion. The
assemble average of the multi-dimensional integral in the expansion is
difficult to compute in general. Thanks to the correlation structure of the
white noise process, we can find the rule for contracting of the multi-
dimensional integral and show the expansion coverges to the moment closure
equation we formally derived before.

The paper is organized as follows. In section 2, we formulate the evolution
equation of passive scalar and the evolution equation of its $N$-point
correlation function. By employing the strategies we described in the previous
paragraph, we derive the ensemble averaged equation of the advection-diffusion
equation with rapidly fluctuating random flows in section 3 and 4, and derive
the ensemble averaged equation of a general evolution equation in section 5.
In section 6, as applications of the averaged equations we derived, we explore
the statistics of the random scalar field which is advected by the random
strain flows and random periodic flows. In the strain flows case, we show that
the statistics have a close relationship with the time integral of geometric
Brownian motion which has an important application in Asian option pricing.
Based on the homogenization method, we derive the effective equation of
$N$-point correlation function and calculate the limiting probability
distribution at long times.

##  2 Preliminaries

We consider the following general evolution equation with initial condition
$T_{I}(\mathbf{x})$

$\displaystyle\partial_{t}T=$ $\displaystyle A(\partial_{\mathbf{x}},t)T,\quad
T(\mathbf{x},0)=$ $\displaystyle T_{I}(\mathbf{x}),$ (1)

where $A(\partial_{\mathbf{x}},t)$ is a random linear operator which depends
on the derivative or higher derivative with respect to the spatial variable
$\mathbf{x}$, for example, the Laplacian
$\Delta_{\mathbf{x}}=\sum\limits_{j=1}^{N}\partial_{x_{j}}^{2}$.

Usually, the random scalar field $T(\mathbf{x},t)$ can be characterized by the
$N$-point correlation function
$\Psi_{N}(\mathbf{x}_{1},..,\mathbf{x}_{N},t)=E_{A}\left(\prod\limits_{j=1}^{N}T(\mathbf{x}_{j},t)\right)$
where $T(\mathbf{x}_{j},t)$ is the solution of equation (1) with the same
realization of $A$ and different initial condition $T_{I}(\mathbf{x}_{j})$.
There are two reasons for that. First, the mean, variance and skewness which
provides important statistical information regarding the random scalar field
can be computed from the first three correlation functions. Second, with the
computed $N$-point correlation function, one can analyze the tail of PDF [8,
9], or even obtain the analytic expression for the full PDF [6, 14] via the
conclusion of moment problem [35].

A nice property of equation (1) which makes easier to obtain the the $N$-point
correlation function is that the product
$\Theta(\mathbf{x}_{1},...,\mathbf{x}_{N},t)=\prod\limits_{j=1}^{N}T(\mathbf{x}_{j},t)$
satisfies a same type of equation with enlarged coordinates space,

$\displaystyle\partial_{t}\Theta=$
$\displaystyle\sum\limits_{j=1}^{N}A(\partial_{\mathbf{x}_{j}},t)\Theta,\quad\left.\Theta\right|_{t=0}=$
$\displaystyle\prod\limits_{j=1}^{N}T_{I}(\mathbf{x}_{j}).$ (2)

Hence, once we could derive the evolution equation for the mean of $T$, we
could also obtain the evolution equation for the mean of $\Theta$ which is the
$N$-point correlation function of the random scalar field $T$.

##  3 Feynman-Kac formula based approach

In this section, we will study the following $d$ dimensional random advection-
diffusion equation with the initial condition $T_{I}\left(\mathbf{x}\right)$,

$\displaystyle\partial_{t}T+\xi(t)\mathbf{v}(\mathbf{x},t)\cdot\nabla
T=\kappa\Delta T,\quad T(\mathbf{x},0)=T_{I}(\mathbf{x}),$ (3)

where $\kappa$ is the diffusivity, $\mathbf{v}(\mathbf{x},t)$ is a
$d$-dimensional deterministic vector field, $\xi(t)$ is a zero-mean, Gaussian
random process with the correlation function given by
$\left\langle\xi(t)\xi(s)\right\rangle=R(t,s)$. A special case of flow is the
shear flow $\xi(t)\mathbf{v}=(u(y)\xi(t),0)$. This type of flow can originate
from either a time varying pressure field, or by randomly moving portions of
the boundary, in a high viscosity fluid [13]. Here, we assume that $\xi(t)$ is
a Gaussian white noise in time so that $R(t,s)=g^{2}\delta(t-s)$. In most
applications of this equation, the white-noise term emerges as limit of a
colored noise. Therefore, we interpret this equation in the sense of
Stratonovich, which will be indicated by the symbol $\circ$ in this paper.

For a fixed realization of velocity $\xi(t)\mathbf{v}(\mathbf{x},t)$, the
Feynman-Kac formula yields the path integral representation for the solution
of equation (3),

$\displaystyle T(\mathbf{x},t)=$ $\displaystyle
E_{\mathbf{B}}(T_{I}(\mathbf{X}(t))),$ (4)

where the random process $\mathbf{X}(s)=\mathbf{X}^{(\mathbf{x},t)}(s)$ is the
solution of stochastic differential equation (SDE),

$\displaystyle\mathrm{d}\mathbf{X}(s)=$
$\displaystyle-\xi(t-s)\mathbf{v}(\mathbf{X}(s),t-s)\mathrm{d}s+\sqrt{2\kappa}\mathrm{d}\mathbf{B}(s),\quad\mathbf{X}(0)=\mathbf{x}.$
(5)

Here $\mathbf{B}(t)$ is the $d$-dimensional standard Brownian motion. Notice
that both the Gaussian white noise process is stationary and temporally
homogeneous, so $\xi(t-s)=\xi(s)$ in the sense of distribution. With the fact
$\xi(s)\mathrm{d}s=g\mathrm{d}\tilde{B}(s)$, where $\tilde{B}(t)$ is
independent of $\mathbf{B}(t)$, we have

$\displaystyle\mathrm{d}\mathbf{X}(s)=$
$\displaystyle-g\mathbf{v}(\mathbf{X}(s),t-s)\circ\mathrm{d}\tilde{B}(s)+\sqrt{2\kappa}\mathrm{d}\mathbf{B}(s).$
(6)

It will be convenient to convert the Stratonovich integral into the Itô
integral,

$\displaystyle\mathrm{d}\mathbf{X}(s)=$ $\displaystyle
g^{2}(\nabla\mathbf{v}(\mathbf{X}(s),t-s))\mathbf{v}(\mathbf{X}(s),t-s)\mathrm{d}t-g\mathbf{v}(\mathbf{X}(s),t-s)\mathrm{d}\tilde{B}(s)+\sqrt{2\kappa}\mathrm{d}\mathbf{B}(s).$
(7)

By utilizing the Feynman-Kac formula, the ensemble average of
$T(\mathbf{x},t)$ with respect to the flow,
$\Psi=E_{\xi(t)}(T(x,t))=E_{\xi(t)}(E_{\mathbf{B}}(T_{I}(\mathbf{X}(t))))$,
satisfies the equation

$\displaystyle\partial_{t}\Psi=\left(\kappa\Delta+\frac{g^{2}}{2}(\mathbf{v}\cdot\nabla)^{2}\right)\Psi,\quad\Psi(\mathbf{x},0)=T_{I}(\mathbf{x}).$
(8)

When the velocity field in equation (3) is the shear flow,
$\xi(t)\mathbf{v}(\mathbf{x},t)=(\xi(t)y,0)$, the advection operator in the
equation of $\Theta=\prod\limits_{j=1}^{N}T(\mathbf{x}_{j},t)$ is
$\sum\limits_{j=1}^{N}y_{j}\partial_{x_{j}}$. With equation (8), we obtain the
conclusion in [26]

$\displaystyle\partial_{t}\Psi=\left(\kappa\Delta+\frac{g^{2}}{2}\left(\sum\limits_{j=1}^{N}y\partial_{x_{j}}\right)^{2}\right)\Psi,\quad\Psi(\mathbf{x},0)=T_{I}(\mathbf{x}).$
(9)

##  4 Limit of Ornstein-Uhlenbeck process

In this section, we present the second method for deriving the moment
evolution equation, which may allow us to study more general random flows.
This method utilizes the fact that the stationary Ornstein-Uhlenbeck (OU)
process converges to the Gaussian white noise process as the correlation time
of the OU process vanishes.

We consider the advection-diffusion equation with more general random flow,

$\displaystyle\partial_{t}T+\mathbf{v}(\mathbf{x},\xi(t),t)\cdot\nabla
T=\kappa\Delta T,\quad T(x,y,0)=T_{I}(x,y).$ (10)

where $\xi(t)$ is the stationary Ornstein-Uhlenbeck process with the damping
$\gamma$ and dispersion $\sigma$ is the solution of stochastic differential
equation (SDE)
$\mathrm{d}\xi(t)=-\gamma\xi(t)\mathrm{d}t+\sigma\mathrm{d}B(t)$ with initial
condition $\xi(0)\sim\mathcal{N}(0,{\sigma^{2}}/{2\gamma})$. Here $B(t)$ is
the standard Brownian motion and $\mathcal{N}(a,b)$ is the normal distribution
with mean $a$ and variance $b$. The correlation function of $\xi(t)$ is
$R(t,s)=\frac{\sigma^{2}}{2\gamma}e^{-\gamma\left|t-s\right|}$. $\gamma^{-1}$
is often referred to as the correlation time of the OU process. It is easy to
check that the stationary Ornstein-Uhlenbeck process converges to the Gaussian
white noise process as the correlation time vanishes with fixed
$g={\sigma}/{\gamma}$. Due to this property, we can obtain the equation for
the white noise case by applying the perturbation method on the equation for
the OU process case.

Resnick [34] derived the PDE for the ensemble average $\Psi=E_{\xi(t)}(T)$ by
utilizing the Feynman-Kac formula. By introducing an extra variable $z$ to
represent the initial value of the stationary OU process, we have
$\Psi(\mathbf{x},t)=\frac{1}{\sqrt{\pi}}\int\limits_{-\infty}^{+\infty}\psi(\mathbf{x},z,t)e^{-z^{2}}\mathrm{d}z$,
where $\psi(\mathbf{x},z,t)$ satisfies the following partial differential
equation:

$\displaystyle\frac{\partial\psi}{\partial
t}+\mathbf{v}(\mathbf{x},g\sqrt{\gamma}z,t)\cdot\nabla\psi+\gamma
z\frac{\partial\psi}{\partial z}$
$\displaystyle=\kappa\Delta\psi+\frac{\gamma}{2}\frac{\partial^{2}\psi}{\partial
z^{2}},\quad\psi(\mathbf{x},z,0)=T_{I}(\mathbf{x}).$ (11)

We note that by using the same strategy, when the flow depends on $M$
independent OU processes, one can obtain the PDE for the $N$-point correlation
function with $M$ new variables $z_{i}$, $1\leq i\leq M$. This generalization
allows one to study the flow which is driven by a more general stochastic
process, for example, the general stationary Gaussian random process which can
be represented as a linear combination of OU processes [34, 12], the Cox-
Ingersoll-Ross process in mathematical finance which is a sum of squared OU
processes (chapter 6 of [21]).

For the linear shear case,
$\mathbf{v}(\mathbf{x},t,\xi(t))\cdot\nabla=\xi(t)y\partial_{x}$, Resnick
derived the exact expression for $\Psi$ and showed it converges to the white
case in the limit $\gamma\rightarrow\infty$ of the damping OU parameter. Here,
we will generalize this conclusion. We use the velocity filed
$\mathbf{v}(\mathbf{x},\xi(t),t)=\xi(t)\mathbf{v}(\mathbf{x},t)$ as an example
to show the derivation of the effective equation in large $\gamma$ limit. For
this type of flow, equation (11) becomes

$\displaystyle\frac{\partial\psi}{\partial
t}+g\sqrt{\gamma}z\mathbf{v}(\mathbf{x},t)\cdot\nabla\psi+\gamma
z\frac{\partial\psi}{\partial z}$
$\displaystyle=\kappa\Delta\psi+\frac{\gamma}{2}\frac{\partial^{2}\psi}{\partial
z^{2}},\quad\psi(\mathbf{x},z,0)=T_{I}(\mathbf{x}).$ (12)

Inspired by the coefficient $\sqrt{\gamma}$ in equation, we assume that the
solution has the series expansion
$\psi=\sum\limits_{n=0}^{\infty}\psi_{n}\gamma^{-\frac{n}{2}}$. Substituting
the expansion into the equation (12) and grouping all terms of the same order
of $\gamma$, we find that we have to solve the sequence of equations. The
first three equations are

$\displaystyle\left(z\partial_{z}-\frac{1}{2}\partial_{z}^{2}\right)\psi_{0}=$
$\displaystyle 0,$ (13)
$\displaystyle\left(z\partial_{z}-\frac{1}{2}\partial_{z}^{2}\right)\psi_{1}=$
$\displaystyle-gz\mathbf{v}(\mathbf{x},t)\cdot\nabla\psi_{0},$
$\displaystyle\left(z\partial_{z}-\frac{1}{2}\partial_{z}^{2}\right)\psi_{2}=$
$\displaystyle-
gz\mathbf{v}(\mathbf{x},t)\cdot\nabla\psi_{1}-\frac{\partial\psi_{0}}{\partial
t}+\kappa\Delta\psi_{0}.$

We have $\psi_{0}(\mathbf{x},t,z)=\psi_{0}(\mathbf{x},t)$,
$\psi_{1}(\mathbf{x},t,z)=-gz\mathbf{v}(\mathbf{x},t)\cdot\nabla\psi_{0}+F(\mathbf{x},t)$.
The equation for $\psi_{2}$ becomes

$\displaystyle\left(z\partial_{z}-\frac{1}{2}\partial_{z}^{2}\right)\psi_{2}=$
$\displaystyle
g^{2}z^{2}\left(\mathbf{v}(\mathbf{x},t)\cdot\nabla\right)^{2}\psi_{0}-gz\mathbf{v}(\mathbf{x},t)\cdot\nabla
F-\frac{\partial\psi_{0}}{\partial t}+\kappa\Delta\psi_{0}.$ (14)

Then $\psi_{2}$ takes this form

$\displaystyle\psi_{2}=(4z^{2}-2)F_{2}+F_{0}-gz\mathbf{v}(\mathbf{x},t)\cdot\nabla
F.$ (15)

Hence, we have

$\displaystyle\psi$
$\displaystyle=\psi_{0}+\frac{-gz\mathbf{v}(\mathbf{x},t)\cdot\nabla\psi_{0}+F(t,\mathbf{x})}{\sqrt{\gamma}}+\frac{(4z^{2}-2)F_{2}+F_{0}-gz\mathbf{v}(\mathbf{x},t)\cdot\nabla
F}{\gamma}+\mathcal{O}(\gamma^{-\frac{3}{2}}).$ (16)

Since the ensemble average $\Psi=E_{\xi(t)}(T)$ is the integral of $\psi$,
$\Psi(\mathbf{x},t)=\frac{1}{\sqrt{\pi}}\int\limits_{-\infty}^{+\infty}\psi(\mathbf{x},z,t)e^{-z^{2}}\mathrm{d}z$,
we take the integration at both side of equation (12)

$\displaystyle\partial_{t}\Psi+\frac{1}{\sqrt{\pi}}\int\limits_{-\infty}^{+\infty}g\sqrt{\gamma}z\mathbf{v}(\mathbf{x},t)\cdot\nabla\psi
e^{-z^{2}}\mathrm{d}z$
$\displaystyle=\kappa\Delta\Psi,\quad\Psi(\mathbf{x},0)=T_{I}(\mathbf{x}).$
(17)

Substituting the expansion (16) into equation (17), we have

$\displaystyle\partial_{t}\Psi-\frac{g^{2}}{2}\left(\mathbf{v}(\mathbf{x},t)\cdot\nabla\right)^{2}\Psi$
$\displaystyle=\kappa\Delta\Psi+\mathcal{O}(\gamma^{-\frac{1}{2}}),\quad\Psi(\mathbf{x},0)=T_{I}(\mathbf{x}).$
(18)

By taking the limit $\gamma\rightarrow\infty$, we obtain equation (8) again.

##  5 General Evolution Equation

In this section, we will show that the conclusions of advection-diffusion
equations we obtained in the previous sections are special cases of a more
general result. We consider
$A(\partial_{\mathbf{x}},t)=C(\partial_{\mathbf{x}},t)+\xi(t)V(\partial_{\mathbf{x}},t)$
in equation (1), $C(\partial_{\mathbf{x}},t),V(\partial_{\mathbf{x}},t)$ are
deterministic linear operators which depend on the derivative or higher
derivative with respect to the spatial variable $\mathbf{x}$, or random linear
operators which are independent of the Gaussian white noise process $\xi(t)$.
Hence, equation (1) becomes,

$\displaystyle\partial_{t}T=$
$\displaystyle\left(C(\partial_{\mathbf{x}},t)+\xi(t)V(\partial_{\mathbf{x}},t)\right)T,\quad
T(\mathbf{x},0)=$ $\displaystyle T_{I}(\mathbf{x}).$ (19)

Obviously, equation (19) reduces to the advection-diffusion equation (3) by
taking $C(\partial_{\mathbf{x}},t)=\Delta$ and
$V\left(\partial_{\mathbf{x}},t)\right)=\mathbf{v}(\mathbf{x},t)\cdot\nabla$.
Another interesting example which can be treated using this method is the
linearized Kuramoto-Sivashinsky equation [48]

$\displaystyle T_{t}=-T_{xx}-\epsilon T_{xxxx}+\xi(t)U(x,t)T_{x}.$ (20)

Such an equation can be analyzed in the constant coefficient case by Fourier
methods [2] but a new approach is needed to treat the non-constant coefficient
case. Also note that path integral based approaches fail, since there is no
known path integral for the above equation.

We can formally write down the solution with the functional of operators
(propagator) $U(t)$,

$\displaystyle T(\mathbf{x},t)=U(t)T_{I}(\mathbf{x})=$
$\displaystyle\exp\left(\int\limits_{0}^{t}C(\partial_{\mathbf{x}},s)+\xi(s)V(\partial_{\mathbf{x}},s)\mathrm{d}s\right)T_{I}$
(21) $\displaystyle=$
$\displaystyle\exp\left(\int\limits_{0}^{t}C(\partial_{\mathbf{x}},s)\mathrm{d}s+g\int\limits_{0}^{t}V(\partial_{\mathbf{x}},s)\mathrm{d}B(s)\right)T_{I}.$

The second step follows the fact
$\xi(s)V(\partial_{\mathbf{x}},s)\mathrm{d}s=gV(\partial_{\mathbf{x}},s)\circ\mathrm{d}B(s)=gV(\partial_{\mathbf{x}},s)\mathrm{d}B(s)$
in the sense of distribution. This functional calculus method was used in
studying the advection-diffusion equation in [42]. Here, we formally treat the
operators as ordinary variables. For a normal random variable
$X\sim\mathcal{N}(\mu,\sigma)$, the characteristic function of $X$ yields
$E_{X}(e^{X})=e^{\mu+\frac{1}{2}\sigma^{2}}$. With this conclusion and Itô
isometry
$E_{B}\left(\int\limits_{0}^{t}V(\partial_{\mathbf{x}},s)\mathrm{d}B(s)\right)^{2}=\int\limits_{0}^{t}V(\partial_{\mathbf{x}},s)^{2}\mathrm{d}s$,
we obtain

$\displaystyle\Psi(\mathbf{x},t)=E_{\xi(t)}(U(t))T_{I}=$
$\displaystyle\exp\left(\int\limits_{0}^{t}C(\partial_{\mathbf{x}},s)\mathrm{d}s+\frac{g^{2}}{2}\int\limits_{0}^{t}V(\partial_{\mathbf{x}},s)^{2}\mathrm{d}s\right)T_{I}.$
(22)

From this expression, we observe that $\Psi$ satisfies the equation

$\displaystyle\partial_{t}\Psi=\left(C(\partial_{\mathbf{x}},t)+\frac{g^{2}}{2}V(\partial_{\mathbf{x}},t)^{2}\right)\Psi,\quad\Psi(\mathbf{x},0)=T_{I}(\mathbf{x}).$
(23)

The above formal proof could be rigorously justified by the pseudo-
differential operator. Alternatively, here, we rigorously prove this result by
evaluating the Picard iteration expansion. The basis for our calculation is
the time ordered product expansion for the propagator $U(t)$, which is given
by

$\displaystyle U(t)=$ $\displaystyle
I+\int\limits_{0}^{t}A(\partial_{\mathbf{x}},s)\mathrm{d}s+\int\limits_{0}^{t}\int\limits_{0}^{s_{2}}A(\partial_{\mathbf{x}},s_{2})A(\partial_{\mathbf{x}},s_{1})\mathrm{d}s_{1}\mathrm{d}s_{2}$
(24)
$\displaystyle+\int\limits_{0}^{t}\int\limits_{0}^{s_{3}}\int\limits_{0}^{s_{2}}A(\partial_{\mathbf{x}},s_{3})A(\partial_{\mathbf{x}},s_{2})A(\partial_{\mathbf{x}},s_{1})\mathrm{d}s_{1}\mathrm{d}s_{2}\mathrm{d}s_{3}+\ldots$

The usual practice at this point is to introduce the time ordering operator,
change the integration to be over the cube of side $t$, rather than over the
$n$-simplex
$\Omega_{n}(t)=\left\\{(s_{1},...,s_{n})|0\leq s_{1}\leq s_{2}\leq
s_{3}\ldots\leq s_{n}\leq t\right\\}$, and introduce a factor of $1/n!$. For
reasons that will become clear shortly we choose to leave the series in the
form written above, with the $n$-fold integration over the $n$-simplex.

Next, our goal is to show that after ensemble average with respect to
$\xi(t)$, the expansion of random propagator (24) becomes

$\displaystyle
E_{\xi(t)}(U(t))=I+\int\limits_{0}^{t}C(\partial_{\mathbf{x}},s)+\frac{g^{2}}{2}V(\partial_{\mathbf{x}},s)^{2}\mathrm{d}s$
(25)
$\displaystyle+\int\limits_{0}^{t}\int\limits_{0}^{s_{2}}\left(C(\partial_{\mathbf{x}},s_{2})+\frac{g^{2}}{2}V(\partial_{\mathbf{x}},s_{2})^{2}\right)\left(C(\partial_{\mathbf{x}},s_{1})+\frac{g^{2}}{2}V(\partial_{\mathbf{x}},s_{1})^{2}\right)\mathrm{d}s_{1}\mathrm{d}s_{2}+...$

In doing so, we compute the ensemble average of the above propagator over the
statistics of $\xi$ utilizing a cluster expansion. We use the fact that $\xi$
has Gaussian white noise statistics to “break the bars” in the usual way.
Generally, it is difficult, if not impossible, to interpret this averaged
propagator as the solution of an evolution equation. However, in the white
noise limit, we will establish that this averaged propagator does in fact obey
a simple evolution equation.

Since $\xi(t)$ is mean zero any terms with an odd number of factor of $V(t)$
vanish upon taking the expectation over $\xi$. The terms which appear at
$n$-th order are products of $2k$ factors of $V(s_{i})$ with $n-2k$ factors of
$C(s_{i})$, multiplied by $k$ delta functions. The arguments of the delta
functions are difference between the arguments of the $V$’s. The integral is
over the $n$-simplex $\Omega_{n}$. For instance the terms appearing at fourth
order are

$\displaystyle
C(s_{4})C(s_{3})C(s_{2})C(s_{1})+\delta(s_{4}-s_{3})V(s_{4})V(s_{3})C(s_{2})C(s_{1})$
(26) $\displaystyle+$
$\displaystyle\delta(s_{4}-s_{2})V(s_{4})C(s_{3})V(s_{2})C(s_{1})+\delta(s_{4}-s_{1})V(s_{4})C(s_{3})C(s_{2})V(s_{1})$
$\displaystyle+$
$\displaystyle\delta(s_{3}-s_{2})C(s_{4})V(s_{3})V(s_{2})C(s_{1})+\delta(s_{3}-s_{1})C(s_{4})V(s_{3})C(s_{2})V(s_{1})$
$\displaystyle+$
$\displaystyle\delta(s_{2}-s_{1})C(s_{4})C(s_{3})V(s_{2})V(s_{1})$
$\displaystyle+$
$\displaystyle\left(\delta(s_{4}-s_{3})\delta(s_{2}-s_{1})+\delta(s_{4}-s_{2})\delta(s_{3}-s_{1})+\delta(s_{4}-s_{1})\delta(s_{3}-s_{2})\right)V(s_{4})V(s_{3})V(s_{2})V(s_{1}).$

The utility of leaving the integration over the simplicity regions now becomes
apparent. The integral over the delta functions can be done easily. The result
of this integration is, of course, the characteristic function of the support
of the delta function. The important thing to realize for the purposes of this
calculation is that the intersection of the $n-1$ dimensional hyperplane
$s_{i}=s_{i+k}$ with the $n$-simplex $0\leq s_{1}\leq s_{2}\leq
s_{3}\ldots\leq s_{n}\leq t$ is the $n-k$ dimensional simplex $0\leq
s_{1}\leq...\leq s_{i}=s_{i+1}=...=s_{i+k}\leq\ldots\leq s_{n}\leq t$. Thus
the intersection of the support of the delta function with the region of
integration is a region of dimension $n-k$. Integrating the characteristic
function of an $n-k$ dimensional set over an $n-1$ dimensional set gives no
contribution unless $k=1$. The upshot is that the only contribution which do
not vanish are those from adjacent factors of $V$ in the product. It is
precisely the non-adjacent terms which make it difficult to interpret the
averaged propagator as the solution to an evolution equation. In the white
noise limit, these terms do not contribute. In addition, we should also notice
that since the support of delta function lies on the boundary of
$\Omega_{n}(t)$, the integration produces the prefactor $1/2$. We summarize
the above discussion as the following formula

$\displaystyle\int\limits_{\Omega_{n}(t)}\delta(s_{i}-s_{j})f(\mathbf{s}_{n})\mathrm{d}\mathbf{s}_{n}=$
$\displaystyle\begin{cases}\frac{1}{2}\int\limits_{\Omega_{n-1}(t)}f(s_{1},...,s_{j-1},s_{i},s_{j+1},...,s_{j})\mathrm{d}\mathbf{s}_{n-1},&i=j+1,\\\
0,&\text{otherwise,}\\\ \end{cases}$ (27)

where $i>j$, $\mathbf{s}_{n}=(s_{1},s_{2},...,s_{n})$. We also relabel $s_{n}$
as $s_{j}$ in the right hand side integral to make the formula more concise.

It is worthwhile to verify the above claim for the first several terms. After
taking ensemble average, the first time integral in (24) becomes

$\displaystyle
E_{\xi(t)}\left(\int\limits_{0}^{t}A(\partial_{\mathbf{x}},s)\mathrm{d}s\right)=$
$\displaystyle\int\limits_{0}^{t}C(\partial_{\mathbf{x}},s)\mathrm{d}s.$ (28)

The second term becomes

$\displaystyle
E_{\xi(t)}\left(\int\limits_{0}^{t}\int\limits_{0}^{s_{2}}A(\partial_{\mathbf{x}},s_{2})A(\partial_{\mathbf{x}},s_{1})\mathrm{d}s_{1}\mathrm{d}s_{2}\right)=\int\limits_{0}^{t}\int\limits_{0}^{s_{2}}C(\partial_{\mathbf{x}},s_{2})C(\partial_{\mathbf{x}},s_{1})\mathrm{d}s_{1}\mathrm{d}s_{2}+\int\limits_{0}^{t}\frac{g^{2}}{2}V(\partial_{\mathbf{x}},s)^{2}\mathrm{d}s.$
(29)

The third term becomes

$\displaystyle
E_{\xi(t)}\left(\int\limits_{\Omega_{3}(t)}\prod\limits_{i=1}^{3}A(\partial_{\mathbf{x}},s_{i})\mathrm{d}\mathbf{s}_{3}\right)=\int\limits_{\Omega_{3}(t)}\prod\limits_{i=1}^{3}C(\partial_{\mathbf{x}},s_{i})\mathrm{d}\mathbf{s}_{3}+\frac{g^{2}}{2}\int\limits_{0}^{t}\int\limits_{0}^{s_{2}}C(\partial_{\mathbf{x}},s_{2})V(\partial_{\mathbf{x}},s_{1})^{2}\mathrm{d}s_{1}\mathrm{d}s_{2}$
(30)
$\displaystyle+\frac{g^{2}}{2}\int\limits_{0}^{t}\int\limits_{0}^{s_{2}}V(\partial_{\mathbf{x}},s_{2})^{2}C(\partial_{\mathbf{x}},s_{1})\mathrm{d}s_{1}\mathrm{d}s_{2}.$

Now, with these examples and conclusion(27), we can show that equation (24)
becomes equation (25) after taking ensemble average. After taking ensemble
average, any term in the cluster expansion of equation (24) is a a integral of
the product of $C(s_{i})$ and $\frac{g^{2}}{2}V(s_{j})$. Clearly, that is a
term in (25). On the other hand, for any term in (25) we can find the unique
corresponding term in equation (24) by writing $\frac{g^{2}}{2}V(s_{i})^{2}$
as $\xi(s_{i})V(s_{i})\xi(s_{i+1})V(s_{i+1})$ and suitable relabeling the
subscripts.

Last, it should be stressed that the crucial element in this Calculation is
the fact that the bracket of non-adjacent $V(s_{i})$’s gives no contribution
after integration. This observation allows us to interpret the averaged series
as the propagator for the above problem. In a more general situation, for
instance one where $\xi(t)$ is Gaussian with some finite time correlation,
terms like

$\displaystyle\int\limits_{0}^{t}\int\limits_{0}^{s_{3}}\int\limits_{0}^{s_{2}}V(s_{3})C(s_{2})V(s_{1})\mathrm{d}s_{1}\mathrm{d}s_{2}\mathrm{d}s_{3},$
(31)

have a non-vanishing contribution. Such terrns do not appear in the series for
the propagator in equation (24), and prevent us from interpreting the
resulting series as a time-ordered exponential.

##  6 Application

In this section, we apply our general result to the case of a passive scalar
diffusing in a rapidly fluctuating fluid flow and study closed equations for
the scalar’s $N$-point correlation function. We are interested in two types of
the flow, the periodic flow and strain flow.

###  6.1 Random periodic flow

For the shear flow $\xi(t)\mathbf{v}(\mathbf{x},t)=\xi(t)(u(y),0)$ and the
parallel-plate channel domain
$\Omega=\left\\{(x,y)|x\in\mathbb{R},y\in[0,1]\right\\}$ with no-flux boundary
conditions, Camassa et al.[10] proved that the $N$-point correlation function
satisfies an effective heat equation at long times. Interestedly, in this
example, the limiting distribution of the scalar field is not Gaussian and
exhibit intermittency. In fact, it is negatively skewed at sufficiently low
Péclet number and could be positively skewed at large Péclet number.

For the parallel-plate channel domain, one can convert the no-flux boundary
condition into periodic boundary condition by evenly extending the scalar
field. Hence, we expect the similar result holds for a general periodic
velocity field. We assume $\mathbf{v}(\mathbf{x},t)$ is a periodic
incompressible velocity field with the spatial period $L_{x}$ and temporal
period $L_{t}$. We consider equation (3) with the slow varying initial
condition $T(\mathbf{x},0)=T_{I}(\frac{\mathbf{x}}{a})$, $a\gg 1$. Notice that
$g^{2}\sim\text{Time}^{-1}$. With the change of variables,

$\displaystyle
a\mathbf{x}^{\prime}=\mathbf{x},\quad\frac{a^{2}}{\kappa}t^{\prime}=t,\quad
U=Lg^{2},\quad\epsilon=\frac{L_{x}}{a},$ (32) $\displaystyle
U\xi^{\prime}(t^{\prime})\mathbf{v}^{\prime}(\mathbf{x}^{\prime},t^{\prime})=\xi(t)\mathbf{v}(\mathbf{x},t),\quad\frac{g\sqrt{\kappa}}{a}\xi^{\prime}(t^{\prime})=\xi(t),$

we can drop the primes without confusion and obtain the nondimensionalized
version of (3):

$\partial_{t}T+\text{Pe}\xi(t)\mathbf{v}\left(\frac{\mathbf{x}}{\epsilon},\frac{t}{\epsilon^{2}}\right)\cdot\nabla
T=\Delta T,\quad T(\mathbf{x},0)=T_{I}(\mathbf{x}),$ (33)

where $\mathbf{v}(\mathbf{x},t)$ has the spatial period $1$ and temporal
period $\kappa L_{t}/L_{x}^{2}$, and we have introduced the Péclet number
$\text{Pe}={UL}/{\kappa}={L^{2}g^{2}}/\kappa$. With equation (8), the ensemble
average of the scalar field $\Psi=E_{\xi(t)}(T)$ satisfies the equation

$\displaystyle\partial_{t}\Psi=\left(\Delta+\mathrm{Pe}^{2}\left(\mathbf{v}\left(\frac{\mathbf{x}}{\epsilon},\frac{t}{\epsilon^{2}}\right)\cdot\nabla\right)^{2}\right)\Psi,\quad\Psi(\mathbf{x},0)=T_{I}(\mathbf{x}).$
(34)

We follow the similar multiscale analysis procedure described in the appendix
7.1 of [13] or [30]. We seek the asymptotic approximation to
$\Psi(\mathbf{x},t)$ in the limit $\epsilon\rightarrow 0$ that has the
following multiscale expansion,

$\Psi(\mathbf{x},t)=\Psi_{0}(\mathbf{x},\mathbf{y},t,\tau)+\epsilon\Psi_{1}(\mathbf{x},\mathbf{y},t,\tau)+\epsilon^{2}\Psi_{2}(\mathbf{x},\mathbf{y},t,\tau)+\mathcal{O}(\epsilon^{3}),$
(35)

with two different scales in $\mathbf{x}$-direction: $\mathbf{x}$ (slow),
$\mathbf{y}={\mathbf{x}}/{\epsilon}$ (fast), and in $t$-direction: $t$ (slow),
$\tau={t}/{\epsilon^{2}}$ (fast), respectively. We use the angle bracket to
denote the average of a function,

$\displaystyle\left\langle\mathbf{v}(\mathbf{x},t)\right\rangle_{\mathbf{x},t}=\int\limits_{0}^{L_{t}}\int\limits_{[0,L_{x}]^{d}}\mathbf{v}(\mathbf{x},t)\mathrm{d}\mathbf{x}\mathrm{d}t.$
(36)

By taking the average on both side of equation (35), we have

$\displaystyle\left\langle\Psi(\mathbf{x},t)\right\rangle_{\mathbf{y},\tau}=\left\langle\Psi_{0}(\mathbf{x},\mathbf{y},t,\tau)\right\rangle_{\mathbf{y},\tau}+\epsilon\left\langle\Psi_{1}(\mathbf{x},\mathbf{y},t,\tau)\right\rangle_{\mathbf{y},\tau}+\epsilon^{2}\left\langle\Psi_{2}(\mathbf{x},\mathbf{y},t,\tau)\right\rangle_{\mathbf{y},\tau}+\mathcal{O}(\epsilon^{3}).$
(37)

Since equation (35) holds for arbitrarily small $\epsilon$, we have
$\left\langle\Psi(\mathbf{x},t)\right\rangle_{\mathbf{y},\tau}=\left\langle\Psi_{0}(\mathbf{x},\mathbf{y},t,\tau)\right\rangle_{\mathbf{y},\tau}$,
and
$\left\langle\Psi_{n}(\mathbf{x},\mathbf{y},t,\tau)\right\rangle_{\mathbf{y},\tau}=0$,
$n\geq 1$. In addition, the differential operators along the $\mathbf{x}$ and
$t$ direction will be replaced,

$\displaystyle\nabla_{\mathbf{x}}\rightarrow\nabla_{\mathbf{x}}+\frac{1}{\epsilon}\nabla_{\mathbf{y}},\quad\Delta_{\mathbf{x}}\rightarrow\Delta_{\mathbf{x}}+\frac{2}{\epsilon}\nabla_{\mathbf{x}}\cdot\nabla_{\mathbf{y}}+\frac{1}{\epsilon^{2}}\Delta_{\mathbf{y}},\quad\partial_{t}\rightarrow\partial_{t}+\frac{1}{\epsilon^{2}}\partial_{\tau}.$
(38)

We would have a hierarchy of equations such that the expansion of
$\Psi(\mathbf{x},t)$ holds for arbitrarily small $\epsilon$. For
$\mathcal{O}(\epsilon^{-2})$, we have:

$\mathcal{L}\Psi_{0}=0,\quad\Psi_{0}(\mathbf{x},\mathbf{y},t,\tau)|_{t=0,\tau=0}=T_{I}(\mathbf{x}),$
(39)

where
$\mathcal{L}\Psi=\left(\partial_{\tau}-\Delta_{\mathbf{y}}-\mathrm{Pe}^{2}\left(\mathbf{v}(\mathbf{y},\tau)\cdot\nabla_{\mathbf{y}}\right)^{2}\right)\Psi$.
Since the initial condition is a function of the variable $\mathbf{x}$ only,
we have $\Psi_{0}(\mathbf{x},\mathbf{y},t,\tau)=\Psi_{0}(\mathbf{x},t)$.

For $\mathcal{O}(\epsilon^{-1})$, we have

$\displaystyle\mathcal{L}\Psi_{1}=2\nabla_{\mathbf{x}}\cdot\nabla_{\mathbf{y}}\Psi_{0}+\mathrm{Pe}^{2}(\mathbf{v}(\mathbf{y},\tau)\cdot\nabla_{\mathbf{x}})(\mathbf{v}(\mathbf{y},\tau)\cdot\nabla_{\mathbf{y}})\Psi_{0}+\mathrm{Pe}^{2}(\mathbf{v}(\mathbf{y},\tau)\cdot\nabla_{\mathbf{y}})(\mathbf{v}(\mathbf{y},\tau)\cdot\nabla_{\mathbf{x}})\Psi_{0},$
(40) $\displaystyle\Psi_{1}(\mathbf{x},\mathbf{y},0,0)=0.$

The first and second term on the right hand side are zero. Fredholm
solvability states that the equation $\mathcal{L}\Psi=f$ has a solution if and
only if $\left\langle fg\right\rangle=0$ for any solution of equation
$\mathcal{L}^{*}g=0$, where $\mathcal{L}^{*}$ is the adjoint operator of
$\mathcal{L}$. Here, the solvability condition is guaranteed by the
periodicity of $\mathbf{v}$,
$\left\langle(\mathbf{v}(\mathbf{y},\tau)\cdot\nabla_{\mathbf{y}})\mathbf{v}(\mathbf{y},\tau)\right\rangle_{\mathbf{y},\tau}=\frac{1}{2}\left\langle\nabla_{\mathbf{y}}(\mathbf{v}(\mathbf{y},\tau)\cdot\mathbf{v}(\mathbf{y},\tau))\right\rangle_{\mathbf{y},\tau}=\mathbf{0}$.
Due to the linearity of the equation, the general form of the solution is
$\Psi_{1}=\bm{\theta}(\mathbf{y},\tau)\cdot\nabla_{\mathbf{x}}\Psi_{0}+C(\mathbf{x},t)$.
Therefore, we have

$\mathcal{L}\theta_{i}=\mathrm{Pe}^{2}(\mathbf{v}(\mathbf{y},\tau)\cdot\nabla_{\mathbf{y}})v_{i}(\mathbf{y},\tau),\quad\theta_{i}(\mathbf{y},0)=0,\quad
1\leq i\leq d,$ (41)

where $\theta_{i}$, $v_{i}$ are $i$th component of $\bm{\theta}$, $\mathbf{v}$
respectively. In addition, we have an extra condition
$\left\langle\bm{\theta}(\mathbf{y},\tau)\right\rangle_{\mathbf{y},\tau}=\mathbf{0}$
from
$\left\langle\Psi_{1}(\mathbf{x},\mathbf{y},t,\tau)\right\rangle_{\mathbf{y},\tau}=0$.
For $\mathcal{O}(\epsilon^{0})$, we have

$\displaystyle\mathcal{L}\Psi_{2}=-\partial_{t}\Psi_{0}+\left(\Delta_{\mathbf{x}}+\mathrm{Pe}^{2}\left(\mathbf{v}(\mathbf{y},\tau)\cdot\nabla_{\mathbf{x}}\right)^{2}\right)\Psi_{0}+2\nabla_{\mathbf{x}}\cdot\nabla_{\mathbf{y}}\Psi_{1}$
(42) $\displaystyle\hskip
34.14322pt+\mathrm{Pe}^{2}(\mathbf{v}(\mathbf{y},\tau)\cdot\nabla_{\mathbf{x}})(\mathbf{v}(\mathbf{y},\tau)\cdot\nabla_{\mathbf{y}})\Psi_{1}+\mathrm{Pe}^{2}(\mathbf{v}(\mathbf{y},\tau)\cdot\nabla_{\mathbf{y}})(\mathbf{v}(\mathbf{y},\tau)\cdot\nabla_{\mathbf{x}})\Psi_{1},$
$\displaystyle\Psi_{2}(\mathbf{x},\mathbf{y},0,0)=0.$

Due to the incompressibility and periodicity of the velocity field, the third
term and last term on the right hand side of the above equation is zero under
the average. The solvability condition yields the effective diffusion equation
with $\epsilon\rightarrow 0$

$\displaystyle\partial_{t}\Psi_{0}$
$\displaystyle=\nabla_{\mathbf{x}}\cdot\left(\Lambda_{d}\nabla_{\mathbf{x}}\Psi_{0}\right),\quad(\Lambda_{d})_{i,j}=\delta_{i,j}+\mathrm{Pe}\left\langle
v_{i}v_{j}+v_{i}\sum\limits_{k=1}^{d}v_{k}\partial_{y_{k}}\theta_{j}\right\rangle_{\mathbf{y},\tau},$
(43) $\displaystyle\Lambda_{d}$
$\displaystyle=\mathbf{I}_{d}+\mathrm{Pe}^{2}\left\langle\mathbf{v}(\mathbf{y},\tau)^{T}\mathbf{v}(\mathbf{y},\tau)+\mathbf{v}(\mathbf{y},\tau)^{T}\left(\mathbf{v}(\mathbf{y},\tau)\nabla_{\mathbf{y}}\right)\bm{\theta}(\mathbf{y},\tau)\right\rangle_{\mathbf{y},\tau},$

where $\mathbf{I}_{d}$ is the $d\times d$ identity matrix. $\delta_{i,j}$ is
Kronecker delta, $\delta_{i,j}=1$ if $i=j$, otherwise, $\delta_{i,j}=0$. Due
to molecular diffusion, the solution of equation (3) with a non-slowly varying
initial condition will be a slowing varying function eventually. Hence,
equation (43) is also the effective equation of equation (3) with a general
initial condition at the diffusion time scale.

Next, we will compute the $N$-point correlation function for a 2 dimensional
scalar field. In doing so, as we discussed in section 2, we need to consider
the problem in higher spatial dimensions. We assume the two dimensional flow
is $(u_{1}(x,y,t),v_{1}(x,y,t))$. The $N$-point correlator will live in the
space
$\mathbf{z}=(x_{1},y_{1},x_{2},y_{2},...,x_{N},y_{N})\in\mathbb{R}^{2N}$. The
flow is

$\displaystyle\mathbf{v}(\mathbf{z},t)=(u_{1}(x_{1},y_{1},t),v_{1}(x_{1},y_{1},t),u_{2}(x_{2},y_{2},t),v_{2}(x_{2},y_{2},t),...,u_{N}(x_{N},y_{N},t),v_{N}(x_{N},y_{N},t)).$
(44)

Then, we can apply the results we derived in previous section to this flow
directly. The cell problem (41) becomes

$\displaystyle\left(\partial_{t}-\partial_{x_{1}}^{2}-\partial_{y_{1}}^{2}-\mathrm{Pe}(u\partial_{x_{1}}+v\partial_{y_{1}})^{2}\right)\theta_{1}=\mathrm{Pe}^{2}(u\partial_{x_{1}}+v\partial_{y_{1}})u,\quad\theta_{1}(\mathbf{z},0)=0,$
(45)
$\displaystyle\left(\partial_{t}-\partial_{x_{1}}^{2}-\partial_{y_{1}}^{2}-\mathrm{Pe}(u\partial_{x_{1}}+v\partial_{y_{1}})^{2}\right)\theta_{2}=\mathrm{Pe}^{2}(u\partial_{x_{1}}+v\partial_{y_{1}})v,\quad\theta_{2}(\mathbf{z},0)=0.$

Notice that $\theta_{1},\theta_{3},...\theta_{2N-1}$ are identical under the
permutation of variables. More precisely, we have

$\displaystyle\theta_{3}(x_{1},y_{1},x_{2},y_{2},...,x_{N},y_{N})=$
$\displaystyle\theta_{1}(x_{2},y_{2},x_{1},y_{1},...,x_{N},y_{N}),$ (46)
$\displaystyle...$
$\displaystyle\theta_{2N-1}(x_{1},y_{1},x_{2},y_{2},...,x_{N},y_{N})=$
$\displaystyle\theta_{1}(x_{N},y_{N},x_{2},y_{2},...,x_{1},y_{1}).$

This properties also holds for $\theta_{2},\theta_{4},...,\theta_{2N}$. With
this fact, once we solve the cell problem (45) for $\theta_{1},\theta_{2}$, we
can determine the leading order approximation of $N$-point correlation
function for arbitrary $N$. The diffusion tensor in the effective equation
(43) now becomes a $N\times N$ block matrix with $N^{2}$ $2\times 2$ matrix as
blocks. More accurately, when $i\neq j$, we have

$\displaystyle\left(\Lambda_{2N}\right)_{2i-1,2j-1}=$
$\displaystyle\mathrm{Pe}^{2}\left(\left\langle
u\right\rangle_{\mathbf{y},\tau}^{2}+\left\langle
u\right\rangle_{\mathbf{y},\tau}\left\langle(u\partial_{x}+v\partial_{y})\theta_{1}\right\rangle_{\mathbf{y},\tau}\right),$
(47) $\displaystyle\left(\Lambda_{2N}\right)_{2i-1,2j}=$
$\displaystyle\mathrm{Pe}^{2}\left(\left\langle
u\right\rangle_{\mathbf{y},\tau}\left\langle
v\right\rangle_{\mathbf{y},\tau}+\left\langle
u\right\rangle_{\mathbf{y},\tau}\left\langle(u\partial_{x}+v\partial_{y})\theta_{2}\right\rangle_{\mathbf{y},\tau}\right),$
$\displaystyle\left(\Lambda_{2N}\right)_{2i,2j-1}=$
$\displaystyle\mathrm{Pe}^{2}\left(\left\langle
v\right\rangle_{\mathbf{y},\tau}\left\langle
u\right\rangle_{\mathbf{y},\tau}+\left\langle
v\right\rangle_{\mathbf{y},\tau}\left\langle(u\partial_{x}+v\partial_{y})\theta_{1}\right\rangle_{\mathbf{y},\tau}\right),$
$\displaystyle\left(\Lambda_{2N}\right)_{2i,2j}=$
$\displaystyle\mathrm{Pe}^{2}\left(\left\langle
v\right\rangle_{\mathbf{y},\tau}^{2}+\left\langle
v\right\rangle_{\mathbf{y},\tau}\left\langle(u\partial_{x}+v\partial_{y})\theta_{2}\right\rangle_{\mathbf{y},\tau}\right).$

When $i=j$, we have

$\displaystyle\left(\Lambda_{2N}\right)_{2i-1,2i-1}=$ $\displaystyle
1+\mathrm{Pe}^{2}\left(\left\langle
u^{2}\right\rangle_{\mathbf{y},\tau}+\left\langle
u(u\partial_{x}+v\partial_{y})\theta_{1}\right\rangle_{\mathbf{y},\tau}\right),$
(48) $\displaystyle\left(\Lambda_{2N}\right)_{2i-1,2i}=$
$\displaystyle\mathrm{Pe}^{2}\left(\left\langle
uv\right\rangle_{\mathbf{y},\tau}+\left\langle
u(u\partial_{x}+v\partial_{y})\theta_{2}\right\rangle_{\mathbf{y},\tau}\right),$
$\displaystyle\left(\Lambda_{2N}\right)_{2i,2i-1}=$
$\displaystyle\mathrm{Pe}^{2}\left(\left\langle
vu\right\rangle_{\mathbf{y},\tau}+\left\langle
v(u\partial_{x}+v\partial_{y})\theta_{1}\right\rangle_{\mathbf{y},\tau}\right),$
$\displaystyle\left(\Lambda_{2N}\right)_{2i,2i}=$ $\displaystyle
1+\mathrm{Pe}^{2}\left(\left\langle
v^{2}\right\rangle_{\mathbf{y},\tau}+\left\langle
v(u\partial_{x}+v\partial_{y})\theta_{2}\right\rangle_{\mathbf{y},\tau}\right).$

When the two dimensional flow is a shear flow
$(u_{1}(x,y,t),v_{1}(x,y,t))=(u(y,t),0)$, the cell problem (45) only has a
trivial solution $\theta_{1}=\theta_{2}=0$. The diffusion tensor would be
extremely simplified. In this case, we have
$\Lambda_{2N}=\mathbf{I}_{2N}+\mathrm{Pe}^{2}\left\langle\mathbf{v}(\mathbf{y},\tau)^{T}\mathbf{v}(\mathbf{y},\tau)\right\rangle_{\mathbf{y},\tau}$,
which reproduces the result in [10, 14, 7]. Then one can compute the PDF of
the scalar field by using the procedure described in [14, 6]. We remark that
[10, 14] only consider the case where $u(y)$ is independent of $t$. Our result
presented in this section generalizes it to the time dependent case, $u(y,t)$.

###  6.2 Random strain flow and time integral of geometric Brownian motion

It is well known that any velocity field can be locally decomposed into the
sum of a shear flow, a strain flow and a rotational flow. Many studies focus
on the random shear flow [26, 31], few studies have addressed the strain flow
or rotational flow. Here we apply our general results to the 2-dimensional
random strain flow which takes the form
$\xi(t)\mathbf{v}(\mathbf{x},t)=\xi(t)(x,-y)$.

The mean of scalar field, $\Psi_{1}(x,y,t)=E_{\xi(t)}[(T(x,y,t)]$, satisfies
the equation

$\displaystyle\partial_{t}\Psi_{1}=\left(\kappa\Delta+\frac{g^{2}}{2}\left(x\partial_{x}-y\partial_{y}\right)^{2}\right)\Psi_{1},\quad\Psi_{1}(x,y,0)=T_{I}(x,y).$
(49)

We can derive the exact solution for a special initial condition,
$T_{I}(x,y)=\delta(x)$. In this case, the solution is independent of $y$, the
equation (49) becomes

$\displaystyle\partial_{t}\Psi_{1}=\left(\kappa\Delta+\frac{g^{2}}{2}(x\partial_{x})^{2}\right)\Psi_{1},\quad\Psi_{1}(x,0)=\delta(x).$
(50)

Interestingly, when $\kappa=0$, $f=e^{\frac{1}{2}g^{2}t}\Psi_{1}(x,-t)$ is the
solution of Black-Scholes equation,

$\displaystyle\partial_{t}f+\frac{1}{2}\sigma^{2}x^{2}\partial_{x}^{2}f+rx\partial_{x}f-rV=0,$
(51)

with the volatility $\sigma=g$ and the risk-free interest rate
$r=\frac{1}{2}g^{2}$. Shortly, we will see that the random strain flow has a
more close relationship to the mathematical finance.

With the change of variable
$z=\int\limits_{0}^{x}\frac{\mathrm{d}s}{\sqrt{\frac{g^{2}s^{2}}{2}+\kappa}}=\frac{\sqrt{2}}{g}\sinh^{-1}\left(\frac{gx}{\sqrt{2\kappa}}\right)$,
equation (50) becomes the heat equation

$\displaystyle\partial_{t}\Psi_{1}=\partial_{z}^{2}\Psi_{1},\quad\Psi_{1}(z,0)=\delta(z)/\sqrt{\kappa}.$
(52)

Hence the solution of equation (50) is

$\displaystyle\Psi_{1}(x,t)=\frac{1}{2\sqrt{\pi\kappa
t}}\exp\left(-\frac{\left(\frac{\sqrt{2}}{g}\sinh^{-1}\left(\frac{gx}{\sqrt{2\kappa}}\right)\right)^{2}}{4t}\right).$
(53)

We can study this random strain flow problem via an alternative approach. With
a fixed realization of white noise $\xi(t)$, we can solve equation (3) with
the strain flow $\mathbf{v}=(x,-y)$ and the line source initial data
$T_{I}(x)=\delta(x-x_{0})$ by utilizing the Fourier transform and method of
characteristics,

$\displaystyle T(x,t)=$
$\displaystyle\frac{1}{\sqrt{4\kappa\pi\int_{0}^{t}e^{-2gB(s)}\mathrm{d}s}}\exp\left(-\frac{e^{-2gB(t)}\left(x-x_{0}e^{gB(t)}\right)^{2}}{4\kappa\int_{0}^{t}e^{-2gB(s)}\mathrm{d}s}\right).$
(54)

First, we consider a special case of this random field. When $x=0,x_{0}=0$, we
have
$X_{t}=T(0,t)=\frac{1}{\sqrt{4\kappa\pi\int_{0}^{t}e^{-2gB(s)}\mathrm{d}s}}$.
This random variable is related to the time integral of geometric Brownian
motion. In mathematical finance, the stock price is usually modeled by the
geometric Brownian motion. For Asian options, the payoff is determined by the
average underlying price over some pre-set period of time. The time integral
of geometric Brownian motion is directly related to the time average of the
stock price which makes it important in Asian option pricing [29].

Some properties of the time integral of geometric Brownian motion could help
us understand the statistical properties of scalar field. Let’s note
$B_{t}^{(\mu)}=B_{t}+\mu t$,
$A_{t}^{(\mu)}=\int_{0}^{t}e^{2B_{s}^{\mu}}\mathrm{d}s$. We omit the
superscript when $\mu=0$. With this notation, the random variable we
interested in is $X_{t}=g\left(4\pi A_{g^{2}t}\right)^{-\frac{1}{2}}$. Using
the time reversal and Itô’s formula (or otherwise), we have the recurrence
relation between negative moments of $A_{t}^{\mu}$ [16, 18].

###### Lemma 6.1

If $2m-\mu\geq 0$, it holds that

$E[(A^{(\mu)}_{t})^{-(m+1)}]=2(m-\mu)E[(A_{t}^{(\mu)})^{-m}]-\frac{1}{m}\frac{\partial}{\partial
t}E[(A_{t}^{(\mu)})^{-m}].$ (55)

Dufresne [17] obtained the following formula of the moments of $A_{t}^{\mu}$,

###### Lemma 6.2

For any $t>0$ and $r\in\mathbb{C}$ with $\Re(r)>-\frac{3}{2}$, it holds that

$E[(2A^{(\mu)}_{t})^{r}]=\frac{e^{-\frac{\mu^{2}t}{2}}}{\sqrt{2\pi
t^{3}}}\int\limits_{0}^{\infty}ye^{-\frac{y^{2}}{2t}}\varphi_{\mu}(r,y)\mathrm{d}y,$
(56)

where, for the Gauss hypergeometric function ${}_{2}F_{1}$ and

$\displaystyle\varphi_{\mu}(r,y)=\frac{\Gamma(r+1)}{\Gamma(2r+2)}e^{-\mu
y}(1-e^{-2y})^{1+2r}{}_{2}F_{1}(\mu+2r+1,1+r;2+2r;1-e^{-2y}),$ (57)
$\displaystyle\varphi_{\mu}(-1,y)=\frac{\cosh((\mu-1)y)}{\sinh(y)}.$

Donati-Martin et al. [16] proved that

###### Lemma 6.3

For any $t>0$, it holds that

$\displaystyle E\left[A_{t}^{-1}\right]=$
$\displaystyle\int\limits_{0}^{\infty}ke^{-k^{2}t/2}\coth\left(\frac{\pi}{2}k\right)\mathrm{d}k.$
(58)

Using lemma 6.2 or Bougerol’s identity [3], we have
$E[A_{t}^{-\frac{1}{2}}]=t^{-\frac{1}{2}}$. Hence,
$E[X_{t}]=g\left(4\pi\right)^{-\frac{1}{2}}\left(g^{2}t\right)^{-\frac{1}{2}}=(4\pi
t)^{-\frac{1}{2}}$, which is consistent with equation (53). With lemma 6.3 and
Laplace’s method [23], we obtain the asymptotic expansion of $E[A_{t}^{-1}]$
as $t\rightarrow\infty$,

$\displaystyle E[A_{t}^{-1}]$ $\displaystyle=\sqrt{\frac{2}{\pi
t}}+\frac{\pi^{3/2}}{6\sqrt{2}t^{3/2}}-\frac{\pi^{7/2}}{120\sqrt{2}t^{5/2}}+\mathcal{O}(t^{-\frac{7}{2}}).$
(59)

Hence, we have

$\displaystyle E[X_{t}^{2}]=$
$\displaystyle\frac{g}{2\sqrt{2}\pi^{3/2}\sqrt{t}}+\frac{\sqrt{\frac{\pi}{2}}}{24gt^{3/2}}-\frac{\pi^{5/2}}{480\sqrt{2}g^{3}t^{5/2}}+\mathcal{O}(t^{-\frac{7}{2}}),\quad
t\rightarrow\infty.$ (60)

The asymptotic expansion of variance as $t\rightarrow\infty$ is

$\displaystyle\mathrm{Var}[X_{t}]$
$\displaystyle=E[X_{t}^{2}]-(E[X_{t}])^{2}=\frac{g}{2\sqrt{2}\pi^{3/2}\sqrt{t}}-\frac{1}{4\pi
t}+\frac{\sqrt{\frac{\pi}{2}}}{24gt^{3/2}}-\frac{\pi^{5/2}}{480\sqrt{2}g^{3}t^{5/2}}+\mathcal{O}(t^{-\frac{7}{2}}).$
(61)

.

The asymptotic expansion of higher moments of $E[A_{t}^{n}]$ can be derived by
utilizing the recurrence relation provided in lemma 6.1,

$\displaystyle
E[A_{t}^{n}]=\frac{\Gamma(n)}{\sqrt{\pi}2^{\frac{1}{2}-n}}t^{-\frac{1}{2}}+\mathcal{O}(t^{-1}),\quad\frac{n}{2}\in\mathbb{Z}^{+}.$
(62)

where $\Gamma(x)$ is the gamma function. This relation also holds for a
positive real number $n$ [29]. Then the leading order expansion of
$E[X_{t}^{n}]$ is

$\displaystyle
E[X_{t}^{n}]=(2\pi)^{-\frac{n}{2}-\frac{1}{2}}g^{n-1}\Gamma\left(\frac{n}{2}\right)t^{-\frac{1}{2}}+\mathcal{O}(t^{-1}).$
(63)

It is easy to see that all normalized moment (e.g., skewness, kurtosis)
diverges as $t\rightarrow\infty$,
$\frac{E\left[\left(X_{t}-E[X_{t}]\right)^{n}\right]}{\mathrm{var}(E_{t})^{n/2}}\sim(2\pi)^{\frac{n-2}{4}}g^{\frac{n}{2}-1}\Gamma\left(\frac{n}{2}\right)t^{\frac{n-2}{4}}$,
$n\geq 3$, meaning that the scalar becomes more and more intermittent during
its decay. Whereas, the normalized moment in the random shear flow problem
converges to a constant at long times [9, 7]. Another example of divergent
normalized moments is given in [37], where the flow is also white in time but
has a more complicated spatial correlation. In that model, the kurtosis
(flatness) grows as $t^{\frac{7}{4}}$.

##  7 Conclusion and discussion

In this paper, we considered an advection-diffusion equation (or a more
general evolution equation) which involves a multiplicative Gaussian white
noise. We have derived the evolution equation for the ensemble average of the
random scalar field by three different approaches, the Feynman-Kac formula,
the limit of the Ornstein-Uhlenbeck process, and the calculus of white noise.
All three different strategies yield the same evolution equation of ensemble-
averaged scalar field, which demonstrates the validity of our results. With
our established evolution equation of ensemble average, we studied the passive
scalar transport problem with two different types of flows, the periodic flow,
and strain flow. For the periodic flow, by using the homogenization method, we
show that the $N$-point correlation function of the random scalar flow is
governed by a diffusion equation at long times, which generalized the
conclusion in [10, 14] from time-independent case to time-dependent case. For
the strain flow, we have explicitly calculated the mean of the random scalar
field and shown the statistic of this random field is related to the time
integral of geometric Brownian motion.

There are several future areas of exploration. First, in this paper, we
focused on the system which is driven by the white noise process. However,
many applications concern the problem which involves a stochastic process with
nonzero correlation. For example, the diffusing diffusivity model models the
diffusivity as a function of OU process[11, 20, 43]. The study of dynamo [4,
49] and the passive-scalar decay problem[46, 15, 1] concerns the renewing
process. The modeling of Black-Scholes market [19] involves fractional
Brownian motion. We expect that the strategies we presented in this paper
could be generalized to problems that involve a more general stochastic
process. Second, the standard homogenization method always yields an effective
equation with constant coefficients. That could limit the valid range of
parameters. Center manifold theory is a powerful tool to handle the scalar
transport problem in bounded domain [32], and to derive the self-similarity
solution of stochastic nonlinear diffusion reaction equation [47]. A recent
study [15] has shown that, for the scalar transport problem in channel domain,
the center manifold theory yields an effective equation with variable
coefficients, which yields a long time limiting distribution of the random
scalar field in the white noise case. Numerical simulation showed that this
approximation performs better than the result from the standard homogenization
method at an intermediate time scale. In addition, the effective equation
suggests that a shear depending upon a Gaussian white noise process or a
renewing process will induce a deterministic effective diffusivity at long
time. We expect that the center manifold theory could yield a better
approximation for the non-shear flow problem provided in equation (33). Last,
we expect this averaged propagator to be useful in the analysis of
fluctuations arising in many physical systems such as the ones alluded to in
this paper.

##  8 Acknowledgements

We acknowledge funding received from the following National Science Foundation
Grant Nos.:DMS-1910824; and Office of Naval Research Grant No: ONR
N00014-18-1-2490. Partial support for Lingyun Ding is gratefully acknowledged
from the National Science Foundation, award NSF-DMS-1929298 from the
Statistical and Applied Mathematical Sciences Institute. We would also like to
thank Andrew J. Majda for his comments and encouragement.

###  8.1 Lists of abbreviations

See table 1.

Full Form | Abbreviation
---|---
Ornstein-Uhlenbeck | OU
Partial differential equation | PDE
Probability distribution function | PDF
Stochastic differential equation | SDE
Table 1: Lists of abbreviations.

## References

  * Aiyer et al. [2017] Aiyer, A.K., Subramanian, K., Bhat, P., 2017.  Passive scalar mixing and decay at finite correlation times in the batchelor regime.  Journal of Fluid Mechanics 824, 785–817.
  * Avellaneda and Majda [1994] Avellaneda, M., Majda, A.J., 1994\.  Simple examples with features of renormalization for turbulent transport.  Philosophical Transactions of the Royal Society of London. Series A: Physical and Engineering Sciences 346, 205–233.
  * Bertoin et al. [2013] Bertoin, J., Dufresne, D., Yor, M., 2013.  Some two-dimensional extensions of bougerol’s identity in law for the exponential functional of linear brownian motion.  Revista matemática iberoamericana 29, 1307–1324.
  * Bhat and Subramanian [2015] Bhat, P., Subramanian, K., 2015\.  Fluctuation dynamos at finite correlation times using renewing flows.  Journal of Plasma Physics 81.
  * Bolles et al. [2019] Bolles, C.T., Speer, K., Moore, M., 2019.  Anomalous wave statistics induced by abrupt depth change.  Physical Review Fluids 4, 011801\.
  * Bronski et al. [2007] Bronski, J.C., Camassa, R., Lin, Z., McLaughlin, R.M., Scotti, A., 2007.  An explicit family of probability measures for passive scalar diffusion in a random flow.  Journal of Statistical Physics 128, 927–968.
  * Bronski and McLaughlin [1997] Bronski, J.C., McLaughlin, R.M., 1997\.  Scalar intermittency and the ground state of periodic schrödinger equations.  Physics of fluids 9, 181–190.
  * Bronski and McLaughlin [2000a] Bronski, J.C., McLaughlin, R.M., 2000a.  The problem of moments and the majda model for scalar intermittency.  Physics Letters A 265, 257–263.
  * Bronski and McLaughlin [2000b] Bronski, J.C., McLaughlin, R.M., 2000b.  Rigorous estimates of the tails of the probability distribution function for the random linear shear model.  Journal of Statistical Physics 98, 897–915.
  * Camassa et al. [2021] Camassa, R., Ding, L., Kilic, Z., McLaughlin, R.M., 2021\.  Persisting asymmetry in the probability distribution function for a random advection-diffusion equation in impermeable channels.  Physica D: Nonlinear Phenomena , 132930doi:https://doi.org/10.1016/j.physd.2021.132930.
  * Chechkin et al. [2017] Chechkin, A.V., Seno, F., Metzler, R., Sokolov, I.M., 2017\.  Brownian yet non-gaussian diffusion: from superstatistics to subordination of diffusing diffusivities.  Physical Review X 7, 021002\.
  * Csáki et al. [1991] Csáki, E., Csörgő, M., Lin, Z., Révész, P., 1991\.  On infinite series of independent ornstein-uhlenbeck processes.  Stochastic processes and their applications 39, 25–44.
  * Ding et al. [2021] Ding, L., Hunt, R., McLaughlin, R.M., Woodie, H., 2021\.  Enhanced diffusivity and skewness of a diffusing tracer in the presence of an oscillating wall.  Research in the Mathematical Sciences 8, 1–29.  doi:https://doi.org/10.1007/s40687-021-00257-4.
  * Ding and McLaughlin [2020] Ding, L., McLaughlin, R.M., 2020\.  Ergodicity and invariant measures for a diffusing passive scalar advected by a random channel shear flow and the connection between the kraichnan-majda model and taylor-aris dispersion.  arXiv preprint arXiv:2012.06610 doi:https://arxiv.org/pdf/2012.06610.pdf.
  * Ding and McLaughlin [2021] Ding, L., McLaughlin, R.M., 2021\.  Determinism and invariant measures for diffusing passive scalars advected by unsteady random shear flows.  arXiv preprint arXiv:2109.05617 doi:https://arxiv.org/pdf/2109.05617.pdf.
  * Donati-Martin et al. [2000] Donati-Martin, C., Matsumoto, H., Yor, M., 2000.  On positive and negative moments of the integral of geometric brownian motions.  Statistics & probability letters 49, 45–52.
  * Dufresne [2000] Dufresne, D., 2000.  Laguerre series for asian and other options.  Mathematical Finance 10, 407–428.
  * Dufresne [2001] Dufresne, D., 2001.  The integral of geometric brownian motion.  Advances in Applied Probability , 223–241.
  * Elliott and Van Der Hoek [2003] Elliott, R.J., Van Der Hoek, J., 2003\.  A general fractional white noise theory and applications to finance.  Mathematical Finance 13, 301–330.
  * Jain and Sebastian [2017] Jain, R., Sebastian, K., 2017\.  Diffusing diffusivity: a new derivation and comparison with simulations.  Journal of Chemical Sciences 129, 929–937.
  * Jeanblanc et al. [2009] Jeanblanc, M., Yor, M., Chesney, M., 2009.  Mathematical methods for financial markets.  Springer Science & Business Media.
  * Kimura and Kraichnan [1993] Kimura, Y., Kraichnan, R.H., 1993\.  Statistics of an advected passive scalar.  Physics of Fluids A: Fluid Dynamics 5, 2264–2277.
  * Kirwin [2010] Kirwin, W.D., 2010.  Higher asymptotics of laplace’s approximation.  Asymptotic Analysis 70, 231–248.
  * Kraichnan [1968] Kraichnan, R.H., 1968.  Small-scale structure of a scalar field convected by turbulence.  The Physics of Fluids 11, 945–953.
  * Majda [1993a] Majda, A.J., 1993a.  Explicit inertial range renormalization theory in a model for turbulent diffusion.  Journal of statistical physics 73, 515–542.
  * Majda [1993b] Majda, A.J., 1993b.  The random uniform shear layer: an explicit example of turbulent diffusion with broad tail probability distributions.  Physics of Fluids A: Fluid Dynamics 5, 1963–1970.
  * Majda and Kramer [1999] Majda, A.J., Kramer, P.R., 1999\.  Simplified models for turbulent diffusion: theory, numerical modelling, and physical phenomena.  Physics reports 314, 237–574.
  * Majda et al. [2019] Majda, A.J., Moore, M., Qi, D., 2019.  Statistical dynamical model to predict extreme events and anomalous features in shallow water waves with abrupt depth change.  Proceedings of the National Academy of Sciences 116, 3982–3987.
  * Matsumoto et al. [2005] Matsumoto, H., Yor, M., et al., 2005\.  Exponential functionals of brownian motion, i: Probability laws at fixed time.  Probability surveys 2, 312–347.
  * McLaughlin [1994] McLaughlin, R.M., 1994.  Turbulent transport.  Ph.D. thesis. Princeton University.
  * McLaughlin and Majda [1996] McLaughlin, R.M., Majda, A.J., 1996\.  An explicit example with non-gaussian probability distribution for nontrivial scalar mean and fluctuation.  Physics of Fluids 8, 536–547.
  * Mercer and Roberts [1990] Mercer, G., Roberts, A., 1990\.  A centre manifold description of contaminant dispersion in channels with varying flow properties.  SIAM Journal on Applied Mathematics 50, 1547–1565.
  * Pumir et al. [1991] Pumir, A., Shraiman, B.I., Siggia, E.D., 1991.  Exponential tails and random advection.  Physical review letters 66, 2984\.
  * Resnick [1996] Resnick, S.G., 1996.  Dynamical problems in non-linear advective partial differential equations.  Ph.D. thesis. The University of Chicago.
  * Shohat and Tamarkin [1943] Shohat, J.A., Tamarkin, J.D., 1943\.  The problem of moments.  1, American Mathematical Soc.
  * Sinai and Yakhot [1989] Sinai, Y.G., Yakhot, V., 1989\.  Limiting probability distributions of a passive scalar in a random velocity field.  Physical review letters 63, 1962\.
  * Son [1999] Son, D., 1999.  Turbulent decay of a passive scalar in the batchelor limit: Exact results from a quantum-mechanical approach.  Physical Review E 59, R3811.
  * Stein et al. [1990] Stein, D., Doering, C.R., Plamer, R., Van Hemmen, J., McLaughlin, R., 1990.  Escape over a fluctuating barrier: the white noise limit.  Journal of Physics A: Mathematical and General 23, L203.
  * Stein et al. [1989] Stein, D., Palmer, R., Van Hemmen, J., Doering, C.R., 1989\.  Mean exit times over fluctuating barriers.  Physics Letters A 136, 353–357.
  * Stephen [1988] Stephen, M.J., 1988.  Temporal fluctuations in wave propagation in random media.  Physical Review B 37, 1\.
  * Sukhatme [2004] Sukhatme, J., 2004.  Probability density functions of decaying passive scalars in periodic domains: An application of sinai-yakhot theory.  Physical Review E 69, 056302\.
  * Taylor [2012] Taylor, M., 2012.  Random walks, random flows, and enhanced diffusivity in advection-diffusion equations.  Discrete & Continuous Dynamical Systems-B 17, 1261.
  * Tyagi and Cherayil [2017] Tyagi, N., Cherayil, B.J., 2017\.  Non-gaussian brownian diffusion in dynamically disordered thermal environments.  The Journal of Physical Chemistry B 121, 7204–7209.
  * Uneyama et al. [2019] Uneyama, T., Miyaguchi, T., Akimoto, T., 2019.  Relaxation functions of the ornstein-uhlenbeck process with fluctuating diffusivity.  Physical Review E 99, 032127\.
  * Vanden Eijnden [2001] Vanden Eijnden, E., 2001.  Non-gaussian invariant measures for the majda model of decaying turbulent transport.  Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences 54, 1146–1167.
  * Vanneste [2006] Vanneste, J., 2006.  Intermittency of passive-scalar decay: Strange eigenmodes in random shear flows.  Physics of fluids 18, 087108\.
  * Wang and Roberts [2013] Wang, W., Roberts, A.J., 2013\.  Self-similarity and attraction in stochastic nonlinear reaction-diffusion systems.  SIAM Journal on Applied Dynamical Systems 12, 450–486.
  * Yakhot [1981] Yakhot, V., 1981.  Large-scale properties of unstable systems governed by the kuramoto-sivashinksi equation.  Physical Review A 24, 642\.
  * Zel’Dovich et al. [1984] Zel’Dovich, Y.B., Ruzmaikin, A., Molchanov, S., Sokoloff, D., 1984\.  Kinematic dynamo problem in a linear velocity field.  Journal of Fluid Mechanics 144, 1–11.


# Superconducting Phase-Diagram of H3S under High Magnetic Fields

Shirin Mozaffari  National High Magnetic Field Laboratory, Florida State
University, Tallahassee, FL 32310, USA  Dan Sun  Los Alamos National
Laboratory, Los Alamos, NM 87545, USA  Vasily S. Minkov  Max-Planck-Institut
fuer Chemie, Hahn-Meitner Weg 1, 55128 Mainz, Germany  Dmitry Knyazev  Max-
Planck-Institut fuer Chemie, Hahn-Meitner Weg 1, 55128 Mainz, Germany
Jonathan B. Betts  Los Alamos National Laboratory, Los Alamos, NM 87545, USA
Mari Einaga  KYOKUGEN, Graduate School of Engineering Science, Osaka
University, Machikaneyamacho 1-3, Toyonaka, Osaka 560-8531, Japan  Katsuya
Shimizu  KYOKUGEN, Graduate School of Engineering Science, Osaka University,
Machikaneyamacho 1-3, Toyonaka, Osaka 560-8531, Japan  Mikhail I. Eremets
Max-Planck-Institut fuer Chemie, Hahn-Meitner Weg 1, 55128 Mainz, Germany
Luis Balicas  National High Magnetic Field Laboratory, Florida State
University, Tallahassee, FL 32310, USA  Fedor F. Balakirev<EMAIL_ADDRESS>Los
Alamos National Laboratory, Los Alamos, NM 87545, USA

###### Abstract

We report the temperature dependence of the upper critical fields
$\mu_{0}H_{c2}(T)$ of the high temperature superconductor H3S under applied
pressures of 155 and 160 GPa through the electrical resistance transition
observed under DC and pulsed magnetic fields up to 65 T, a record high
combination of fields and pressures. We find that $H_{c2}(T)$ generally
follows the Werthamer, Helfand and Hohenberg (WHH) formalism at low fields,
albeit with noticeable deviations upon approaching our experimental limit of
$\mu_{0}H=65$ T. In fact, $H_{c2}(T)$ displays a remarkably linear dependence
on temperature over an extended temperature range also found in multigap or in
strongly-coupled superconductors. The best fit of $H_{c2}(T)$ to the WHH
formula yields a negligible value for the Maki parameter $\alpha$ and for
spin-orbit scattering constant $\lambda_{\text{SO}}$. However, its behavior is
relatively well-described by a model based on strong coupling
superconductivity with a coupling constant $\lambda\sim 2$. Therefore, we
conclude that H3S behaves as a strong-coupled orbital-limited superconductor
over the entire range of temperatures and fields used for our measurements.

The ongoing scientific quest to stabilize superconductivity at room
temperature led to the discovery of superconductivity, with a very high
critical temperature $T_{c}$ = 203 K, in sulfur hydride H3S under high
pressures of $p$ = 155 GPa Drozdov  _et al._ (2015). H3S along with other
hydridesDrozdov  _et al._ (2018); Somayazulu  _et al._ (2019) seems to be the
closest compound, so far, to metallic hydrogen, which is predicted to be a
high temperature superconductorAshcroft (1968); Ginzburg (1969); Maksimov and
Savrasov (2001). H3S forms as the result of the chemical instability of H2S
under high pressures, where H2S decomposes into elemental sulfur S and H3SLi
_et al._ (2014); Duan  _et al._ (2014); Errea  _et al._ (2015). The structure
of H3S is believed to be body-centered cubic Im-$\bar{3}$m over the pressure
range of 92 - 173 GPa Einaga  _et al._ (2016), and characterized by H atoms
situated midway between two body-centered S atoms. The pressure as a function
of temperature phase diagram of H3S is asymmetric, meaning that $T_{c}$ shows
a sharp increase from 95 K to 203 K in the pressure range of 110 - 155 GPa,
but decreases with further increasing the pressureDrozdov  _et al._ (2015).

The superconductivity in H3S is believed to be conventional, in which pairing
is mediated by phonons, with high-frequency optical modes due to the motion of
hydrogenDrozdov  _et al._ (2015); Gorí Ûkov and Kresin (2016); Flores-Livas,
JoséA.  _et al._ (2016); Akashi  _et al._ (2015); Gor’kov and Kresin (2018).
The high $T_{c}$ arises from both the metallization of the strongly covalent
bonds in H3SBernstein  _et al._ (2015) and the high phonon frequencies
displayed by its light elementsQuan and Pickett (2016). Band structure
calculations indicate that H3S is a multiband metal Bianconi and Jarlborg
(2015) having a large Fermi surface (broad energy dispersive bands) as the
result of the hybridization between the H 1s and the 3p orbitals of sulfurDuan
_et al._ (2014); Errea  _et al._ (2015); Quan and Pickett (2016); Heil and
Boeri (2015); Flores-Livas  _et al._ (2016); Akashi  _et al._ (2015). Band
structure calculations also yield small Fermi surface pockets for the
high-$T_{c}$ phase of H3S Bianconi and Jarlborg (2015).

Despite the very high $T_{c}$s reported for H3S, the studies of $H_{c2}(T)$
are limited to a narrow range of temperatures close to $T_{c}$ due to the
inherent experimental difficulties in performing ultra-high pressure
measurements under very high magnetic fields. The behavior of $H_{c2}(T)$
provides valuable information such as an estimation of the Cooper pair
coherence length, the strength of the electron-phonon coupling, the role of
the spin-orbit coupling, and the dominant mechanism breaking the Copper pairs.

In type-II superconductors, the orbital and the spin-paramagnetic pair-
breaking effects are the two main mechanisms depairing electrons under high
magnetic fields. The orbital pair breaking effect explains the suppression of
superconductivity via the formation of Abrikosov vorticesAbrikosov (2004) in
the presence of a field. The superconductivity is suppressed when the kinetic
energy associated with vortex currents exceeds the condensation energy of the
paired electrons. When the magnetic field approaches a critical value
$H_{c2}^{\text{orb}}=\phi_{0}/2\pi\xi^{2}$, referred to as the orbital
limiting field, the vortex cores overlap and the system returns to the normal
state. Here, $\phi_{0}$ is the flux quantum and $\xi$ is the coherence length.
The spin-paramagnetic effect explains the effect of the magnetic field on the
spin of the electron based on the classical work by ClogstonClogston (1962)
and ChandrasekharChandrasekhar (1962). When the magnetic energy exceeds the
superconducting gap, i.e. $1/2\chi_{p}H_{p}^{2}=1/2N(E_{F})\Delta^{2}$,
superconductivity will be suppressed, where $\chi_{p}$, $N(E_{F})$, and
$\Delta$ are the normal state paramagnetic susceptibility, the density of
states at the Fermi level, and the superconducting gap, respectively. If we
only consider the spin-paramagnetic effect, the zero-temperature Pauli
limiting field (Chandrasekhar-Clogston limit) for a weakly coupled
superconductor is approximately $\mu_{0}H_{p}(0)=1.86~{}T_{c}$, where
$2\Delta(T=0)\sim 3.52~{}k_{B}T_{c}$ (for a conventional BCS superconductor)
and $\chi_{p}=g\mu_{B}^{2}N(E_{F})$ ($g$ is the Landé $g$ factor and $\mu_{B}$
is the Bohr magneton). $H_{p}$ can be renormalized by the strength of the
electron phonon-coupling or by electronic correlations. In many
superconductors $H_{c2}(T)$ is affected by both pair-breaking mechanisms. The
MakiMaki (1966) parameter $\alpha=\sqrt{2}H_{c2}^{\text{orb}}(0)/H_{p}(0)$ is
a measure of the relative strength between the orbital and the paramagnetic
pair breaking mechanisms for a given type-II superconductor. The WHH
theoryWerthamer  _et al._ (1966) includes both pair-breaking effects through
the Maki parameter and the spin-orbit constant.

Here, we report measurements of the upper critical field in H3S samples under
extremely high pressures in high magnetic fields. The ultrahigh pressure can
only be obtained by a diamond anvil cell (DAC), however, measurements on
samples contained by a DAC are quite challenging, particularly in pulsed
fields due to the narrow magnet bore and large induced currents which heat the
cell. The small diameter of the DACDrozdov  _et al._ (2015) employed in our
study allows us to perform transport measurements up to 65 T without heating
the sample significantly. Our pulsed field measurements under $p>150$ GPa is a
significant pressure increase over latest pulsed field achievement of 4
GPaBraithwaite  _et al._ (2016). We find that H3S is an orbital limited
superconductor over the entire temperature range and likely a multigap or
strongly-coupled superconductor.

Two samples of H3S were synthesized in-situ inside DACs using two different
techniques. The first sample was prepared from condensed liquid H2S via
disproportionation reaction as described in Ref. [Drozdov  _et al._ , 2015],
the final pressure inside the sample was as high as 160 GPa. We refer to this
sample as the 160 GPa sample. The second sample was synthesized directly from
elemental sulfur and hydrogen at high pressure. For this sample a small piece
of elemental sulfur (purity of 99.98%) with a lateral dimension of about 20
$\mu$m and a thickness of $\sim 3-5$ $\mu$m was placed in a DAC having
diamonds with culets of $60-70$ $\mu$m. Excess hydrogen was introduced in the
DAC at a gas pressure of $\sim 130-150$ MPa. Subsequently, the DAC was
pressurized up to 150 GPa and then heated up to 1000 K with a YAG laser at
room temperature to initiate the chemical reaction. Pressure increased
slightly to 155 GPa after the synthesis. The vibrational properties of the
pressurized initial reactants and of the synthesized products were probed
using a triple grating Raman spectrometer equipped by a HeNe laser with a
wavelength of 633 nm. Sputtered gold electrodes were thoroughly isolated from
the metal gasket by a layer made from magnesium oxide, calcium fluorite and an
epoxy glue mixture. This layer also prevented hydrogen penetration into the
rhenium gasket. The pressure was estimated using the Raman shift of stressed
diamondEremets (2003) and the H2 vibron’s wavenumber, previously calibrated in
a separate experimentEremets and Troyan (2011). Both scales indicated a
pressure of 155 GPa after the laser heating. Under continuous fields up to
$\mu_{0}H=35$ T the resistance of the sample was measured using a commercial
AC resistance bridge. A custom made Lock-In amplifier was used under pulsed
fields up $\mu_{0}H=65$ T.

Figure 1: (Color online) Temperature dependence of resistance of the H3S
samples under 155 GPa and 160 GPa pressures. Red line is a fit to
$R=R_{0}+AT^{2}$ indicating Fermi liquid behavior. Note that $T_{c}$ for the
155 GPa sample shifts slightly to larger values over time. Figure 2: (Color
online) (a) Hall resistance $R_{xy}$ as a function of the magnetic field
$\mu_{0}H$ at a temperatures $T=205$ snd 250 K. Red line is a linear fit. (b)
Hall coefficient $R_{H}$ as a function of $T$ for the H3S sample pressurized
up to 155 GPa.

Figure 1 displays the resistance as a function of the temperature $T$ for H3S
samples pressurized up to 155 GPa and 160 GPa. The onset, i.e. first deviation
from normal state resistivity, of the superconducting transitions are
estimated to be $T_{c}=201$ K for the sample under 155 GPa, and $T_{c}=174$ K
for the other sample under 160 GPa. The drop in $T_{c}$ for the higher
pressure sample is consistent with the results in Ref. [Drozdov  _et al._ ,
2015]. The width of these transitions are $\Delta T_{c}=5.5$ K and 26 K for
the samples under 155 GPa and 160 GPa, respectively indicating that the latter
sample is less homogeneous, or that it is subjected to stronger pressure
gradients. Both samples display Fermi-liquid behavior in the normal state as
indicated by the red line which corresponds to a fit to $R(T)=R_{0}+AT^{2}$.
The value of $T_{c}$ for the sample under 155 GPa increases slightly over
time, indicating an evolution towards a more homogeneous sample or weaker
pressure gradients. Notice that this sample displays a step around 195 K in
the superconducting transition indicating that its superconducting state
indeed is inhomogeneous. For samples synthesized within the confines of the
DAC a certain degree of inhomogeneity is inevitable.

Figure 2 (a) displays the Hall resistance $R_{xy}$ for the sample under
$p=155$ GPa as a function of magnetic field $\mu_{0}H$ at 205 K or near the
superconducting transition. The $R_{xy}(\mu_{0}H)$ does not reveal any non-
linearity as a function of $\mu_{0}H$ all the way up to 35 T. This indicates
that near and above $T_{c}$ the electrical transport in H3S is dominated by
single electron-like carrier pocket despite the fact that H3S was predicted to
be a multiband metalBianconi and Jarlborg (2015). Most likey, this indicates
that electrons display a considerably higher carrier mobility than holes at
these temperatures. The temperature dependence of the Hall coefficient $R_{H}$
is shown in Fig. 2 (b). We observed a noticeable decrease in $R_{H}(T)$ as the
temperature is decreased, corresponding to an increase in the effective Hall
density $n_{H}=1/eR_{H}$. The decrease in the Hall coefficient at low $T$s
could come from a decrease in electron mobility, a progressive increase in
hole mobility, or from an evolution in their relative densities. If one
assumed that only one band contributed to the transport of carriers, the Hall
coefficient would yield an electron density $n=1/eR_{H}=8.5\times 10^{22}$
cm-3 at room temperature, which is relatively close to known values for
transition metals such as copper. At 205 K, the carrier concentration would
increase to $n=8.9\times 10^{22}$ cm-3.

Figure 3: (Color online) Resistance as a function of field $\mu_{0}H$ for the
H3S samples under (a) $p=155$ GPa and (b) 160 GPa and for several
temperatures. This data was collected under continuous and pulsed field.

To determine the temperature dependence of $H_{c2}$, we measured the
isothermal resistance as a function of the magnetic field $\mu_{0}H$ at
selected temperatures ranging from $T=60$ K to 200 K. Figures 3 (a) and (b)
show the magnetic field dependence of the resistance for the 155 GPa and 160
GPa samples under fixed temperatures. At each temperature the resistance
changes from zero to a finite value as $\mu_{0}H$ increases due to the field-
induced suppression of superconductivity. $H_{c2}$ is defined as the
intersection between an extrapolation of the resistance of the normal state
and a line having the slope the resistive transition at its middle point Meier
_et al._ (2016). The same criterium is used across DC and pulsed field traces.

The resistive onset has been found to match thermodynamic $H_{c2}$ obtained by
other experimental probesMielke  _et al._ (2001); Altarawneh  _et al._ (2008);
Grissonnanche  _et al._ (2014). This criterion prevents the contribution of
possible superconducting vortex related phases. The resistive transition
shifts to higher fields as the temperature is lowered while broadening
slightly. The $H_{c2}(T)$ values obtained from pulsed field measurements agree
with the curvature and the values extracted under DC fields. At very high
pulsed magnetic fields the self-heating of the metallic DACs by the induced
currents becomes noticeable at temperatures below 60 K.

Figure 4: (Color online) Upper critical fields $H_{c2}$ as a function of
reduced temperature $t=T/T_{c}$ for the H3S samples under 155 GPa and 160 GPa
pressure. Solid lines are fits to theoretical models: Ginzburg-Landau (GL),
Werthamer, Helfand and Hohenberg (WHH), and strong strong-coupling model using
the coupling strength parameter $\lambda=2$ (reproduced from a prediction
based on strong-coupling theory Bulaevskii  _et al._ (1988); Thomas  _et al._
(1996)). Inset: H3S upper-critical fields for two different pressures plotted
in reduced variables. Red line is a linear fit.

Figures 4 shows the extracted $H_{c2}$ as a function of $T$ for both samples
under $p=155$ GPa and 160 GPa as a function of the reduced temperature
$t=T/T_{c}$. To avoid the influence of possible superconducting vortex phases,
we plot just the onset of their resistive transition. These phase boundaries
were obtained from both isothermal field scans and from temperature scans
under fixed magnetic fields. For both samples $H_{c2}(T)$ increases almost
linearly upon decreasing the temperature under fields all the way up to 60 T
yielding the slopes $|dH_{c2}/dT|_{T_{c}}=0.62$ and 0.45 T/K for the samples
under 155 GPa and 160 GPa, respectively. The linear dependence of $H_{c2}$
over an extended range of temperatures was also observed in the two-band
superconductor MgB2Lee  _et al._ (2001); Gurevich (2003) and in the multigap
Fe-pnictide superconductors for fields along certain, but not for all
crystallographic orientations, and claimed to result from the orbital limiting
effect Zhang  _et al._ (2016); Jaroszynski  _et al._ (2008); Yuan  _et al._
(2009). Multigap effects might explain our phase-diagram for H3S in analogy
with those systems. However, multiband scenario contrast with the linear in
field Hall data in Fig. 2 (a), which indicates conduction dominated by a
single type of carrier. Measurements of, for example, the penetration depth as
a function of $T$ should help clarify the nature of the superconducting state
in this system. The estimated weak coupling Pauli limiting fields for the
samples under 155 GPa and 165 GPa are $H_{p}(0)\simeq$ 370 T and 344 T,
respectively. In the dirty limit, the $H_{c2}(0)$ values derived from the
slope of the $H$-$T$ phase boundary at $T_{c}$ are,
$H_{c2}^{orb}(0)=0.69~{}|dH_{c2}/dT|_{T_{c}}\times T_{c}=$ 88 T for the sample
under 155 GPa and 74 T for the one under 160 GPa. These latter values are much
closer to the extrapolation of our experimental $H_{c2}(T)$ towards zero-
temperature, when compared to the Pauli limiting values which differ by a
factor of $\sim 5$, thus indicating that H3S is an orbital limited
superconductor.

Subsequently, we fit our experimentally determined phase-boundary to two
different expressions in order to estimate the value of $H_{c2}$ at $T=0$ K.
First, the conventional Ginzburg-Landau (GL) expression:

$H_{c2}(T)=\dfrac{\phi_{0}}{2\pi\bar{\xi}(0)^{2}}\left(1-t^{2}\right),$ (1)

where $\phi_{0}$ is the quantum of flux and $\bar{\xi}(0)$ is an average
Ginzburg-Landau coherence length at $T=0$.

And second WHH formalism where the temperature dependence of $H_{c2}$ defined
by orbital and spin-paramagnetic effects in the dirty limit is given by
Werthamer  _et al._ (1966):

$\begin{split}\ln\left(\dfrac{1}{t}\right)&=\sum_{\nu=-\infty}^{\infty}\Bigg{\\{}\dfrac{1}{|2\nu+1|}-\Bigg{[}|2\nu+1|\\\
&+\dfrac{\bar{h}}{t}+\dfrac{(\alpha\bar{h}/t)^{2}}{|2\nu+1|+(\bar{h}+\lambda_{so})/t}\Bigg{]}^{-1}\Bigg{\\}}\end{split}$
(2)

where, $\bar{h}=(4/\pi^{2})[H_{c2}(T)/(-dH_{c2}/dT)_{T_{c}}]$, $\alpha$ is the
Maki parameter, and $\lambda_{\text{SO}}$ is the spin-orbit constant.

At temperatures near $T_{c}$, our experimental $H_{c2}$ values seem to follow
the usual GL-expression, but at low temperatures they deviate from it
considerably. This is particularly apparent for the sample under 160 GPa, for
which we were able to reach smaller $t$ values. Comparing the upper critical
fields to the WHH expression yields very small values for $\alpha$ which
points again to rather small orbital limiting fields relative to the Pauli
limiting ones, with no apparent saturation in $H_{c2}(t)$ at low $t$s. We
obtained the best fits to the WHH formula for $\alpha$ values ranging between
0.0 and 0.3 for the sample under 155 GPa and $\alpha=0-0.2$ for the sample
under 160 GPa. For such small values of $\alpha$, the WHH formula is almost
insensitive to the strength of the spin-orbit interaction. As for the effect
of pressure on the H3S samples, our experimental upper critical fields do not
reveal any change in the relative strength of both pair breaking mechanisms
between both pressures. As seen in the inset of Fig. 4, the superconducting
phase diagrams for both samples fall into a single curve when plotted in
reduced units. What is remarkable is the linearity of the phase-boundary as a
function of $t$ all the way down to $t\sim 0.25$. Again, this resembles
behavior reported, for example, for the Fe based superconductors Yuan  _et
al._ (2009).

This deviation from conventional WHH could be also explained by fact that it
was largely formulated in the weak limit of electron-phonon coupling constant
$\lambda\ll 1$. Recent ab initio calculations report substantial enhancement
of $\lambda$ in H3S due to the proximity to a structural instability between
Im-$\bar{3}$m and $R3m$ crystal symmetriesErrea  _et al._ (2016). As phonon
frequencies $\omega_{p}$ soften near this instability, Cooper pairing is
weakened at elevated temperatures by thermal phonons, while it remains robust
in $T\rightarrow 0$ limit, leading to a relative enhancement of $H_{c2}$ at
low temperatures. The calculated $\lambda\sim 2$ Errea  _et al._ (2015, 2016)
would correspond to a nearly linear $H_{c2}(T)$ for $t>0.25$ within the
framework of the strong coupling modelBulaevskii  _et al._ (1988); Thomas  _et
al._ (1996), which was also found to be applicable in an unusually strong
electron-phonon coupling case in Bi-III Brown  _et al._ (2018).
Coincidentally, the strong coupling model also predicts a very high
$T_{c}\sim$300-400 K for a lattice containing light atoms like HBulaevskii
_et al._ (1988). Our experimental $H_{c2}(T)$ values display a much better
agreement with the temperature dependence predicted by the strong coupling
modelBulaevskii  _et al._ (1988) with an electron-phonon coupling constant
$\lambda=2$. We obtained the superconducting coherence length by extrapolating
the strong-coupling fit to zero-temperature. This extrapolation of
$H_{c2}(T)$s yields $H_{c2}=97$ T and 73 T for the H3S samples pressurized up
to 155 GPa and 160 GPa, respectively. Thus, we find a coherence length
$\xi=1.84$ nm and of 2.12 nm for the samples under 155 GPa and 160 GPa,
respectively.

In summary, we have investigated the temperature dependence of the upper
critical fields of H3S under magnetic fields up to 65 T. At lower fields the
phase-boundary separating normal and superconducting states is relatively well
described by the Werthamer, Helfand and Hohenberg formula. Pronounced
deviations from the WHH formula are observed at lower temperatures due to the
linearity of $H_{c2}(T)$. Overall, the phase boundary between superconducting
and metallic states indicates that the orbital-effect suppresses
superconductivity over the entire temperature range. The linearity of
$H_{c2}(T)$ observed over an extended range of reduced temperatures suggests
that this system might indeed be a multigap superconductor as predicted
theoretically and in analogy with similar results for the Fe based
superconductors. Alternatively, enhanced electron-phonon coupling and
softening of hydrogen vibrational modes in the vicinity of structural
instability could also explain the deviation from WHH predictions. We extract
values for the coherence length $\xi$ ranging between $\xi=18.4~{}\text{\AA}$
and 21.2 Å. Above the transition temperature the Hall-effect remains linear up
to high magnetic fields indicating that H3S is a very good metal and that at
high temperatures carrier conduction is dominated by electrons.

We acknowledge A. V. Gurevich and H. H. Wen for useful discussions. This work
was performed at the National High Magnetic Field laboratory, which is
supported by National Science Foundation Cooperative Agreement No.
DMR-1644779, and the State of Florida. S.M. acknowledges support from the FSU
Provost Postdoctoral Fellowship Program. L.B. is supported by DOE-BES through
award DE-SC0002613. M.I.E. is thankful to the Max Planck community for the
invaluable support and U. Pöschl for the constant encouragement.

## References

## References

  * Drozdov  _et al._ (2015) A. P. Drozdov, M. I. Eremets, I. A. Troyan, V. Ksenofontov, and S. I. Shylin, Nature 525, 73 (2015).
  * Drozdov  _et al._ (2018) A. P. Drozdov, P. P. Kong, V. S. Minkov, S. P. Besedin, M. A. Kuzovnikov, S. Mozaffari, L. Balicas, F. Balakirev, D. Graf, V. B. Prakapenka, E. Greenberg, D. A. Knyazev, M. Tkacz, and M. I. Eremets, arXiv:1812.01561 (2018).
  * Somayazulu  _et al._ (2019) M. Somayazulu, M. Ahart, A. K. Mishra, Z. M. Geballe, M. Baldini, Y. Meng, V. V. Struzhkin, and R. J. Hemley, Phys. Rev. Lett. 122, 027001 (2019).
  * Ashcroft (1968) N. W. Ashcroft, Phys. Rev. Lett. 21, 1748 (1968).
  * Ginzburg (1969) V. L. Ginzburg, Phys. Usp. 12, 241 (1969).
  * Maksimov and Savrasov (2001) E. Maksimov and D. Savrasov, Solid State Commun. 119, 569 (2001).
  * Li  _et al._ (2014) Y. Li, J. Hao, H. Liu, Y. Li, and Y. Ma, J. Chem. Phys. 140, 174712 (2014).
  * Duan  _et al._ (2014) D. Duan, Y. Liu, F. Tian, D. Li, X. Huang, Z. Zhao, H. Yu, B. Liu, W. Tian, and T. Cui, Sci. Rep. 4, 6968 (2014).
  * Errea  _et al._ (2015) I. Errea, M. Calandra, C. J. Pickard, J. Nelson, R. J. Needs, Y. Li, H. Liu, Y. Zhang, Y. Ma, and F. Mauri, Phys. Rev. Lett. 114, 157004 (2015).
  * Einaga  _et al._ (2016) M. Einaga, M. Sakata, T. Ishikawa, K. Shimizu, M. I. Eremets, A. P. Drozdov, I. A. Troyan, N. Hirao, and Y. Ohishi, Nat. Phys. 12, 835 (2016).
  * Gorí Ûkov and Kresin (2016) L. P. Gorí Ûkov and V. Z. Kresin, Sci. Rep. 6, 25608 (2016).
  * Flores-Livas, JoséA.  _et al._ (2016) Flores-Livas, JoséA., Sanna, Antonio, and Gross, E. K.U., Eur. Phys. J. B 89, 63 (2016).
  * Akashi  _et al._ (2015) R. Akashi, M. Kawamura, S. Tsuneyuki, Y. Nomura, and R. Arita, Phys. Rev. B 91, 224513 (2015).
  * Gor’kov and Kresin (2018) L. P. Gor’kov and V. Z. Kresin, Rev. Mod. Phys. 90, 011001 (2018).
  * Bernstein  _et al._ (2015) N. Bernstein, C. S. Hellberg, M. D. Johannes, I. I. Mazin, and M. J. Mehl, Phys. Rev. B 91, 060511 (2015).
  * Quan and Pickett (2016) Y. Quan and W. E. Pickett, Phys. Rev. B 93, 104526 (2016).
  * Bianconi and Jarlborg (2015) A. Bianconi and T. Jarlborg, EPL (Europhysics Letters) 112, 37001 (2015).
  * Heil and Boeri (2015) C. Heil and L. Boeri, Phys. Rev. B 92, 060508 (2015).
  * Flores-Livas  _et al._ (2016) J. A. Flores-Livas, M. Amsler, C. Heil, A. Sanna, L. Boeri, G. Profeta, C. Wolverton, S. Goedecker, and E. K. U. Gross, Phys. Rev. B 93, 020508 (2016).
  * Abrikosov (2004) A. A. Abrikosov, Rev. Mod. Phys. 76, 975 (2004).
  * Clogston (1962) A. M. Clogston, Phys. Rev. Lett. 9, 266 (1962).
  * Chandrasekhar (1962) B. S. Chandrasekhar, Appl. Phys. Lett. 1, 7 (1962).
  * Maki (1966) K. Maki, Phys. Rev. 148, 362 (1966).
  * Werthamer  _et al._ (1966) N. R. Werthamer, E. Helfand, and P. C. Hohenberg, Phys. Rev. 147, 295 (1966).
  * Braithwaite  _et al._ (2016) D. Braithwaite, W. Knafo, R. Settai, D. Aoki, S. Kurahashi, and J. Flouquet, Rev. Sci. Instrum. 87, 023907 (2016).
  * Eremets (2003) M. I. Eremets, J. Raman Spectrosc. 34, 515 (2003).
  * Eremets and Troyan (2011) M. I. Eremets and I. A. Troyan, Nat. Mater. 10, 927 (2011).
  * Meier  _et al._ (2016) W. R. Meier, T. Kong, U. S. Kaluarachchi, V. Taufour, N. H. Jo, G. Drachuck, A. E. Böhmer, S. M. Saunders, A. Sapkota, A. Kreyssig, M. A. Tanatar, R. Prozorov, A. I. Goldman, F. F. Balakirev, A. Gurevich, S. L. Bud’ko, and P. C. Canfield, Phys. Rev. B 94, 064501 (2016).
  * Mielke  _et al._ (2001) C. Mielke, J. Singleton, M.-S. Nam, N. Harrison, C. C. Agosta, B. Fravel, and L. K. Montgomery, J. Phys. Condens. Matter 13, 8325 (2001).
  * Altarawneh  _et al._ (2008) M. M. Altarawneh, K. Collar, C. H. Mielke, N. Ni, S. L. Bud’ko, and P. C. Canfield, Phys. Rev. B 78, 220505 (2008).
  * Grissonnanche  _et al._ (2014) G. Grissonnanche, O. Cyr-Choinière, F. Laliberté, S. René de Cotret, A. Juneau-Fecteau, S. Dufour-Beauséjour, M.-È. Delage, D. LeBoeuf, J. Chang, B. J. Ramshaw, D. A. Bonn, W. N. Hardy, R. Liang, S. Adachi, N. E. Hussey, B. Vignolle, C. Proust, M. Sutherland, S. Krämer, J.-H. Park, D. Graf, N. Doiron-Leyraud, and L. Taillefer, Nat. Commun. 5, 3280 EP (2014).
  * Bulaevskii  _et al._ (1988) L. N. Bulaevskii, O. V. Dolgov, and M. O. Ptitsyn, Phys. Rev. B 38, 11290 (1988).
  * Thomas  _et al._ (1996) F. Thomas, B. Wand, T. Lühmann, P. Gegenwart, G. R. Stewart, F. Steglich, J. P. Brison, A. Buzdin, L. Glémot, and J. Flouquet, J. Low Temp. Phys. 102, 117 (1996).
  * Lee  _et al._ (2001) S. Lee, H. Mori, T. Masui, Y. Eltsev, A. Yamamoto, and S. Tajima, J. Physical. Soc. Japan. 70, 2255 (2001).
  * Gurevich (2003) A. Gurevich, Phys. Rev. B 67, 184515 (2003).
  * Zhang  _et al._ (2016) Q. R. Zhang, D. Rhodes, B. Zeng, M. D. Johannes, and L. Balicas, Phys. Rev. B 94, 094511 (2016).
  * Jaroszynski  _et al._ (2008) J. Jaroszynski, S. C. Riggs, F. Hunte, A. Gurevich, D. C. Larbalestier, G. S. Boebinger, F. F. Balakirev, A. Migliori, Z. A. Ren, W. Lu, J. Yang, X. L. Shen, X. L. Dong, Z. X. Zhao, R. Jin, A. S. Sefat, M. A. McGuire, B. C. Sales, D. K. Christen, and D. Mandrus, Phys. Rev. B 78, 064511 (2008).
  * Yuan  _et al._ (2009) H. Q. Yuan, J. Singleton, F. F. Balakirev, S. A. Baily, G. F. Chen, J. L. Luo, and N. L. Wang, Nature 459, 565 (2009).
  * Errea  _et al._ (2016) I. Errea, M. Calandra, C. J. Pickard, J. R. Nelson, R. J. Needs, Y. Li, H. Liu, Y. Zhang, Y. Ma, and F. Mauri, Nature 532, 81 (2016).
  * Brown  _et al._ (2018) P. Brown, K. Semeniuk, D. Wang, B. Monserrat, C. J. Pickard, and F. M. Grosche, Sci. Adv. 4 (2018), 10.1126/sciadv.aao4793.