# On the semiclassical regularity of thermal equilibria

Jacky J. Chong  ,  Laurent Lafleche  Department of Mathematics, The University
of Texas at Austin, Austin, TX 78712, USA<EMAIL_ADDRESS><EMAIL_ADDRESS>and  Chiara Saffirio  Department of Mathematics and
Computer Science, University of Basel, 4051 Basel, Switzerland,
<EMAIL_ADDRESS>
###### Abstract.

We study the regularity properties of fermionic equilibrium states at finite
positive temperature and show that they satisfy certain semiclassical bounds.
As a corollary, we identify explicitly a class of positive temperature states
satisfying the regularity assumptions of [J.J. Chong, L. Lafleche, C.
Saffirio: arXiv:2103.10946 (2021)].

###### Key words and phrases:

semiclassical limit, fermionic thermal states, quasi-free Gibbs states,
equilibrium

###### 2010 Mathematics Subject Classification:

82C10, 35Q41, 35Q55 (82C05,35Q83).

######  Table of Contents

  1. 1 Introduction
    1. 1.1 Second quantization formalism
    2. 1.2 Harmonic oscillators
    3. 1.3 Thermal states
    4. 1.4 Quasi-free states
  2. 2 Semiclassical thermal regularity
    1. 2.1 Preliminaries
    2. 2.2 Bounds on the inverse fugacity
    3. 2.3 $\mathcal{L}^{\infty}$ bounds
    4. 2.4 $\mathcal{L}^{p}$ bounds for $2\leq p<\infty$

##  1\. Introduction

We consider a system of $N$ non-interacting fermions in a harmonic trap and
study the semiclassical structure of its one-particle density matrix
associated with the corresponding Gibbs (equilibrium) state of the trapped gas
at temperature $T>0$. More precisely, given a family of one-particle
normalized nonnegative trace class operators $\boldsymbol{\rho}_{\beta}$
acting on $L^{2}(\mathbb{R}^{d})$ and indexed by the parameter $\beta$,
proportional to $T^{-1}$, we consider the following quantities

(1)
$\boldsymbol{\nabla}_{\\!x}\boldsymbol{\rho}_{\beta}:=\left[{\nabla,\boldsymbol{\rho}_{\beta}}\right]\quad\text{
and
}\quad\boldsymbol{\nabla}_{\\!\xi}\boldsymbol{\rho}_{\beta}:=\left[{\tfrac{x}{i\hbar},\boldsymbol{\rho}_{\beta}}\right],$

which we refer to as the quantum gradients of $\boldsymbol{\rho}_{\beta}$, and
show that for any $p\in[2,\infty]$,

$\boldsymbol{\rho}_{\beta},\,\sqrt{\boldsymbol{\rho}_{\beta}}\in\mathcal{W}^{1,p}(\mathsf{m})$

uniformly in $\hbar$, where $\mathcal{W}^{1,p}(\mathsf{m})$ are the
semiclassical analogue of weighted Sobolev spaces equipped with the norm

$\left\lVert\boldsymbol{\rho}\right\rVert_{\mathcal{W}^{1,p}(\mathsf{m})}^{p}:=\left\lVert\boldsymbol{\rho}\,\mathsf{m}\right\rVert_{\mathcal{L}^{p}}^{p}+\left\lVert\boldsymbol{\nabla}_{\\!x}\boldsymbol{\rho}\,\mathsf{m}\right\rVert_{\mathcal{L}^{p}}^{p}+\left\lVert\boldsymbol{\nabla}_{\\!\xi}\boldsymbol{\rho}\,\mathsf{m}\right\rVert_{\mathcal{L}^{p}}^{p}$

and $\mathsf{m}:=1+\left\lvert\boldsymbol{p}\right\rvert^{n}$. Here
$\boldsymbol{p}=-i\hbar\nabla$ is the momentum operator and $\mathcal{L}^{p}$
are the semiclassical versions of the Lebesgue spaces endowed with the
rescaled Schatten norms

$\left\lVert\boldsymbol{\rho}\right\rVert_{\mathcal{L}^{p}}:=h^{\frac{3}{p}}\operatorname{Tr}\\!\left(\left\lvert\boldsymbol{\rho}\right\rvert^{p}\right)^{\frac{1}{p}}$

where $h=2\pi\hbar$.

The quantities (1) play key roles in obtaining semiclassical estimates, such
as in the problems of the derivation of the Vlasov equation from quantum
mechanics and of the Hartree and the Hartree–Fock equations from many-body
quantum mechanics. To the author’s knowledge, they first appeared in [3],
where the Hartree equation was derived in the mean-field regime for arbitrary
long times from a system of $N$ fermions interacting through a smooth
potential, under the assumption that the initial state
$\boldsymbol{\rho}^{\mathrm{in}}$ of the Hartree dynamics is a pure state,
i.e. $(\boldsymbol{\rho}^{\mathrm{in}})^{2}=\boldsymbol{\rho}^{\mathrm{in}}$,
and satisfies the following semiclassical structures

(2)
$\left\lVert\boldsymbol{\nabla}_{\\!x}\boldsymbol{\rho}^{\mathrm{in}}\right\rVert_{\mathcal{L}^{1}}\leq
C\quad\text{ and
}\quad\left\lVert\boldsymbol{\nabla}_{\\!\xi}\boldsymbol{\rho}^{\mathrm{in}}\right\rVert_{\mathcal{L}^{1}}\leq
C.$

Similarly, bounds on semiclassical Schatten norms of commutators are used in
[10] to obtain uniform in $\hbar$ mean-field estimates, in [13] in the context
of the semiclassical limit, in [9] in the context of the convergence of
numerical schemes, in [11] to get estimates on the quantum Wasserstein
distance.

When the interaction is the Coulomb potential only partial results are
available on a time scale of order one (see [14]) and they rely on a
semiclassical structure expressed in terms of $L^{p}$ norms of the diagonal of
the operator $|\boldsymbol{\nabla}_{\\!\xi}\boldsymbol{\rho}|$, with $p>3$
(see also [13]).

In [8] and [1], the authors obtained semiclassical bounds of the form (2) for
pure states in the non-interacting case, which in turn provided explicit
examples of initial pure states that satisfy the assumptions in [3] on the
semiclassical structures. In this paper, we extend this analysis to the finite
positive temperature case, where the one-particle reduced density matrix
$\boldsymbol{\rho}$ is no longer a projection. When dealing with mixed states
the analogue relevant quantities considered in [2, 6] are

(3)
$\left\lVert\boldsymbol{\nabla}_{\\!x}\sqrt{\boldsymbol{\rho}^{\mathrm{in}}}\right\rVert_{\mathcal{L}^{p}}\leq
C\quad\text{ and
}\quad\left\lVert\boldsymbol{\nabla}_{\\!\xi}\sqrt{\boldsymbol{\rho}^{\mathrm{in}}}\right\rVert_{\mathcal{L}^{p}}\leq
C.$

In [2], the interaction is smooth and Equation (3) is required to hold true
for the the initial state for $p=1$. In [6], singular interactions of the form
$\left\lvert x\right\rvert^{-a}$ with $a\leq 1$ are considered and Equation
(3) is required to hold for $p\in[2,\infty]$.

In this work, we focus on quasi-free Gibbs states, which is an important class
of examples as they model the equilibria of an ideal gas at finite positive
temperature. It is well-known that the one-particle reduced density matrix
associated to a Gibbs state can be computed explicitly (see Proposition 1.2).
Using the explicit expression (4), we prove bounds on the weighted
$\mathcal{L}^{p}$ norms of the quantum gradients (1) with
$\boldsymbol{\rho}_{\beta}$ given by the one-particle reduced density matrix
associated to the Gibbs state of a system of non-interacting fermions confined
by a harmonic potential, and on its partition function. Our main result is the
following.

######  Theorem 1.1.

Let $H=\frac{\left\lvert\boldsymbol{p}\right\rvert^{2}+\left\lvert
x\right\rvert^{2}}{2}$ be the $d$-dimensional harmonic oscillator Hamiltonian
and define the family of density matrix operators given by

(4)
$\boldsymbol{\rho}_{\beta}:=\dfrac{1}{N\,h^{d}}\left(\mathbf{1}+e^{\beta(H-\mu)}\right)^{-1}$

where $\mu$ depends on $\beta$, $N$ and $h$ and is chosen so that
$h^{d}\operatorname{Tr}(\boldsymbol{\rho}_{\beta})=1$ holds. Then for any
$p\geq 2$, $\boldsymbol{\rho}_{\beta}\in\mathcal{W}^{1,p}(\mathsf{m})$ and
$\sqrt{\boldsymbol{\rho}_{\beta}}\in\mathcal{W}^{1,p}(\mathsf{m})$ uniformly
with respect to $\hbar\in(0,1)$. More explicitly, let

(5) $Z_{\beta}:=h^{d}\operatorname{Tr}(e^{-\beta H}),\quad\text{ and }\quad
Z_{\mu}:=\lambda\,e^{-\beta\mu}$

with $\lambda:=Nh^{d}$, then the following holds. For any fixed
$\beta\in\mathbb{R}_{+}$, $\hbar\in(0,1)$ and for $p\in[1,\infty]$, we have
the bound

(6)
$\left\lVert\boldsymbol{\nabla}\boldsymbol{\rho}_{\beta}\right\rVert_{\mathcal{L}^{p}}\leq
C_{d,p}\,\frac{\beta^{\frac{1}{2}-\frac{d}{p}}}{Z_{\mu}}\frac{\max\left(2\sqrt{2},\beta\hbar\right)^{\frac{1}{2}-\frac{1}{p}}}{\left(\theta(\beta\hbar)\right)^{\frac{1}{p}}}$

where $\theta(x)=\operatorname{th}(x)/x$ with
$\operatorname{th}(x)=\tfrac{e^{x}-e^{-x}}{e^{x}+e^{-x}}$, and

(7) $\displaystyle C_{\lambda,\beta}^{-1}\,Z_{\beta}\leq Z_{\mu}\leq
Z_{\beta}$

with $C_{\lambda,\beta}=2$ if $\mu\leq d\hbar/2$ and
$C_{\lambda,\beta}=1+e^{\beta\lambda^{1/d}/\pi}$ if $\mu\geq d\hbar/2$.

######  Remark 1.1.

As a corollary of our result, we exhibit a class of initial states satisfying
the regularity assumptions of [6, Theorem 3.1].

######  Remark 1.2.

Our proof in the case $p=\infty$ is rather general and in particular can be
adapted without difficulty to other external trapping potentials $V(x)$, that
is to the case of an Hamiltonian of the form
$H=\frac{\left\lvert\boldsymbol{p}\right\rvert^{2}}{2}+V(x)$.

######  Remark 1.3.

As indicated in [6, Figure 1], there are three scaling regimes to keep in mind
when studying mean-field dynamics of quantum systems (both bosonic and
fermionic), namely, $Nh^{d}\rightarrow 0,Nh^{d}\rightarrow\infty$, and
$Nh^{d}$ converges to some constant (in fact, $Nh^{d}=1$ which we called the
critical scaling).

Let $\beta>0$ be fixed. In the lower density scaling regime $Nh^{d}\rightarrow
0$, we have that $e^{\beta\mu}\rightarrow 0$ or equivalently
$\mu\rightarrow-\infty$, which is an immediate consequence of (7) in the above
theorem and the fact that $Z_{\beta}$ is bounded (cf. Section 2.2). It is not
difficult to see that in this case
$\boldsymbol{\rho}_{\beta}=\left(Nh^{d}+Z_{\mu}\,e^{\beta H}\right)^{-1}$ will
converge to $Z_{\beta}^{-1}\,e^{-\beta H}$. In other words, in the
semiclassical limit $\hbar\rightarrow 0$, $\boldsymbol{\rho}_{\beta}$ will
converge to the Maxwell–Boltzmann distribution. This result should not be
surprising to the reader since in the lower density scaling regime $Nh^{d}\to
0$, the Fermi gas behaves like a classical ideal gas.

In the case of the higher density scaling regime $Nh^{d}\rightarrow\infty$, we
recall that as in [6, Equation (16)],
$\lambda=Nh^{d}\leq\frac{1}{\left\lVert\boldsymbol{\rho}\right\rVert_{\mathcal{L}^{\infty}}}=:\mathcal{C}_{\infty}^{-1}$.
This shows that
$\left\lVert\boldsymbol{\rho}\right\rVert_{\mathcal{L}^{\infty}}\rightarrow 0$
as $\lambda\rightarrow\infty$. In short, the higher density scaling regime is
not accessible to Fermi gases. From a more physical perspective, a high
density Fermi gas exhibits quantum “degeneracy” behavior, consequence of the
Pauli exclusion principle, which limits the gas from getting arbitrarily
dense.

Finally, in the critical scaling case, we see that $\boldsymbol{\rho}_{\beta}$
does not converge to the Maxwell–Boltzmann distribution but instead to the
Fermi–Dirac distribution. The divergence from classical behavior is attributed
to the fact that at the critical scaling we are in fact modeling the
degeneracy behavior of the Fermi gas. See e.g. [5, Chapter 17.5].

######  Remark 1.4.

In the semiclassical limit $\hbar\to 0$, then our main Inequality (6) implies

$\left\lVert\boldsymbol{\nabla}\boldsymbol{\rho}_{\beta}\right\rVert_{\mathcal{L}^{p}}\leq
C_{d,p,\lambda}\,\beta^{\frac{1}{2}+\frac{d}{p^{\prime}}}+o(\hbar)$

which is optimal in the sense that it is consistent with the behavior of the
classical norm $\left\lVert\nabla\left(Z_{\beta}^{-1}e^{-\beta\left\lvert
z\right\rvert^{2}}\right)\right\rVert_{L^{p}(\mathbb{R}^{2d})}$ in terms of
$\beta$. It is however interesting to see that this scaling seems to change at
low temperature and when $\hbar$ is not negligible (i.e. when $\beta\hbar$ is
large).

###  1.1. Second quantization formalism

To introduce quasi-free Gibbs states, we will need the second quantized
description of a many-body system of non-interacting fermions in an external
trap. Here we give a brief review of the Fock space formalism. For a more
comprehensive treatment, the reader could consult [4, 15].

The fermionic (antisymmetric) Fock space over the one-particle Hilbert space
$\mathfrak{h}=L^{2}(\mathbb{R}^{d})$ is defined to be the closure of

$\operatorname{\mathcal{F}}:=\bigoplus^{\infty}_{n=0}\mathfrak{h}^{\wedge
n}=\mathbb{C}\oplus\mathfrak{h}\oplus(\mathfrak{h}\wedge\mathfrak{h})\oplus\cdots$

with respect to the norm induced by the standard associated inner product.
Here the $n$-th sector of $\operatorname{\mathcal{F}}$ given by
$\mathfrak{h}^{\wedge n}$ denotes the subspace of
$\mathfrak{h}^{n}=L^{2}(\mathbb{R}^{dn})$ containing functions that satisfy
the antisymmetry property

$\psi(x_{1},\ldots,x_{n})=\operatorname{sgn}(\sigma)\,\psi(x_{\sigma(1)},\ldots,x_{\sigma(n)})$

for any permutation $\sigma$. An important vector to note is the vacuum vector
$\Omega=(1,0,0,\ldots)\in\operatorname{\mathcal{F}}$ which describes a pure
state with no particles.

The operator-valued distribution $a_{x}$ and $a^{\ast}_{x}$ on
$\operatorname{\mathcal{F}}$ are defined by their actions on the
$n^{\text{th}}$ sector of $\operatorname{\mathcal{F}}$ as follows. Let
$\psi\in\mathfrak{h}^{\wedge n}$, then

$\displaystyle(a_{x}\psi)(x_{1},\ldots,x_{n-1})$
$\displaystyle:=\sqrt{n}\,\psi(x,x_{1},\ldots,x_{n-1}),$
$\displaystyle(a_{x}^{\ast}\psi)(x_{1},\ldots,x_{n+1})$
$\displaystyle:=\sum^{n+1}_{j=1}\frac{(-1)^{j-1}}{\sqrt{n+1}}\,\delta_{x}(x_{j})\,\psi(x_{1},\ldots,\cancel{x_{j}},\ldots,x_{n+1}).$

Furthermore, for each $\varphi\in\mathfrak{h}$, we associate to it the
annihilation operator and its adjoint the creation operator, defined by

$a(\varphi)=\int_{\mathbb{R}^{d}}\overline{\varphi(x)}\,a_{x}\,\mathrm{d}x\quad\text{
and }\quad
a^{\ast}(\varphi)=\int_{\mathbb{R}^{d}}\varphi(x)\,a^{\ast}_{x}\,\mathrm{d}x.$

It can be shown that creation and annihilation operators are bounded on
$\operatorname{\mathcal{F}}$ with operator norm $\left\lVert
a(\varphi)\right\rVert_{\infty}=\left\lVert
a^{\ast}(\varphi)\right\rVert_{\infty}=\left\lVert\varphi\right\rVert_{\mathfrak{h}}$
and that they satisfy the canonical anticommutation relations, i.e.

$\\{a(\varphi),a(\phi)\\}=\\{a^{\ast}(\varphi),a^{\ast}(\phi)\\}=0\quad\text{
and
}\quad\\{a^{\ast}(\varphi),a(\phi)\\}=\left\langle\phi,\varphi\right\rangle_{\mathfrak{h}}\mathbf{1}$

hold for all $(\varphi,\phi)\in\mathfrak{h}^{2}$ where
$\\{\mathsf{A},\mathsf{B}\\}:=\mathsf{A}\mathsf{B}+\mathsf{B}\mathsf{A}$ is
the standard anticommutator of the operators $\mathsf{A}$ and $\mathsf{B}$.

In this work, we are interested in trapped non-interacting Fermionic systems
modeled by the one-particle Hamilton operator

(8) $H=\tfrac{1}{2}\left\lvert\boldsymbol{p}\right\rvert^{2}+V(x),$

where $\boldsymbol{p}=-i\hbar\nabla_{x}$ and $V$ is an external trapping
potential. The second quantization of the one-particle operator $H$ on
$\operatorname{\mathcal{F}}$ is defined to be the operator
$\mathsf{H}=\operatorname{\mathrm{d}\Gamma}(H):=\int_{\mathbb{R}^{d}}a^{\ast}_{x}Ha_{x}\,\mathrm{d}x$
whose action on the $n^{\text{th}}$ sector is given by

$\left(\operatorname{\mathrm{d}\Gamma}(H)\Psi\right)(x_{1},\ldots,x_{n})=\sum^{n}_{j=1}(-\tfrac{\hbar^{2}}{2}\Delta_{x_{j}}+V(x_{j}))\Psi(x_{1},\ldots,x_{n}).$

Another useful operator on $\operatorname{\mathcal{F}}$ is the number operator

$\mathcal{N}=\bigoplus^{\infty}_{n=0}n\,\mathbf{1}_{\mathfrak{h}^{\wedge n}}.$

###  1.2. Harmonic oscillators

To simplify the presentation, we choose the external potential $V$ to be the
harmonic trapping potential, that is,

(9)
$H=\tfrac{1}{2}\left\lvert\boldsymbol{p}\right\rvert^{2}+\tfrac{1}{2}\left\lvert
x\right\rvert^{2}=\sum^{d}_{i=1}\tfrac{1}{2}\,\boldsymbol{p}_{i}^{2}+\tfrac{1}{2}\,x_{i}^{2}=\sum^{d}_{i=1}H_{i}.$

Since $H$ is the $d$-dimensional quantum harmonic oscillator, it is well-known
that it emits a discrete spectrum where the eigenvalues and eigenvectors are
given by

$\displaystyle E_{n}$ $\displaystyle=\left(\left\lvert
n\right\rvert_{1}+\tfrac{d}{2}\right)\hbar$ $\displaystyle\text{for
}n=(n_{1},\ldots,n_{d})\in\mathbb{N}_{0}^{d},$ $\displaystyle u_{E_{n}}$
$\displaystyle=\phi_{n_{1}}\otimes\cdots\otimes\phi_{n_{d}}$
$\displaystyle\text{with
}H_{i}\phi_{n_{i}}=\left(n_{i}+\tfrac{1}{2}\right)\hbar\,\phi_{n_{i}}.$

where $\left\lvert n\right\rvert_{1}=n_{1}+\ldots+n_{d}$ and the functions
$\phi_{n_{i}}$ are the standard one-variable Hermite functions

$\phi_{n_{i}}(x_{i})=c_{n_{i}}\,e^{-x^{2}/2}P_{n_{i}}(x_{i})\quad\text{ with
}\quad
c_{n_{i}}^{-1}=2^{\frac{n_{i}}{2}}(n_{i}!)^{\frac{1}{2}}\left(\pi\hbar\right)^{\frac{1}{4}},$

with $P_{n_{i}}(x)$ the Hermite polynomials. Moreover, for each eigenvalue
$E_{n}$ the corresponding multiplicity is given by $g_{\left\lvert
n\right\rvert_{1},d}=\binom{\left\lvert n\right\rvert_{1}+d-1}{d-1}$. In
particular, by the spectral decomposition of $H$, we could now rewrite
$\mathsf{H}$ as

(10)
$\mathsf{H}=\sum_{n\in\mathbb{Z}_{0}^{d}}E_{n}\,a^{\ast}(u_{E_{n}})\,a(u_{E_{n}}).$

###  1.3. Thermal states

We consider quantum states that are given by density operators
$\boldsymbol{\omega}$ on $\operatorname{\mathcal{F}}$, that is, the
expectation of an observable $\mathsf{A}$ is given by
$\left\langle\mathsf{A}\right\rangle_{\boldsymbol{\omega}}=\operatorname{Tr}_{\operatorname{\mathcal{F}}}\left(\mathsf{A}\,\boldsymbol{\omega}\right)$
for every $\mathsf{A}\in\mathcal{B}(\operatorname{\mathcal{F}})$, where
$\mathcal{B}(\operatorname{\mathcal{F}})$ denotes the space of bounded
operators on the Fock space. The Gibbs equilibrium of a trapped non-
interacting Fermi gas associated to some positive temperature $T>0$ is defined
to be the unique minimizer of the Gibbs free energy functional

$\mathscr{F}(\boldsymbol{\omega})=\operatorname{Tr}_{\operatorname{\mathcal{F}}}\left(\mathsf{H}\,\boldsymbol{\omega}\right)-k_{\mathrm{B}}\,T\,\mathscr{S}(\boldsymbol{\omega})\quad\text{
with
}\quad\mathscr{S}(\boldsymbol{\omega})=-\operatorname{Tr}_{\operatorname{\mathcal{F}}}\left(\boldsymbol{\omega}\,\log\boldsymbol{\omega}\right)$

in
$\mathcal{X}_{N}=\\{\boldsymbol{\omega}\in\mathcal{B}(\operatorname{\mathcal{F}})\mid\boldsymbol{\omega}\geq
0,\
\operatorname{Tr}_{\operatorname{\mathcal{F}}}\left(\boldsymbol{\omega}\right)=1,\
\operatorname{Tr}_{\operatorname{\mathcal{F}}}\left(\mathcal{N}\,\boldsymbol{\omega}\right)=N\\}$
where $\mathscr{S}$ is the von Neumann entropy. It can be checked that the
Gibbs state associated to the temperature $T$ is given by the normalized
positive trace class operator

(11)
$\boldsymbol{\omega}_{N}=\frac{1}{\mathcal{Z}_{N}}\,e^{-\beta\left(\mathsf{H}-\mu\mathcal{N}\right)}$

where $\beta=1/(k_{\mathrm{B}}\,T)$ and the chemical potential $\mu$ is chosen
so that
$\operatorname{Tr}\\!\left(\mathcal{N}\boldsymbol{\omega}_{N}\right)=N$. Here
$\mathcal{Z}_{N}$ is the grand canonical partition function

$\mathcal{Z}_{N}=1+\sum_{n=1}^{\infty}e^{n\beta\mu}\operatorname{Tr}_{\mathfrak{h}^{\wedge
n}}\exp\bigg{(}-\beta\sum^{dn}_{j=1}H_{j}\bigg{)}=\prod_{n\in\mathbb{Z}_{0}^{d}}(1+e^{-\beta(E_{n}-\mu)}).$

###  1.4. Quasi-free states

Here we will only give a rudimentary introduction to the tools necessary for
the subsequent sections. A more comprehensive exposition of the following
content can be found in [4, 15]. The state $\boldsymbol{\omega}$ is said to be
quasi-free if for all $n\in\mathbb{N}$, we have that

$\displaystyle\small\operatorname{Tr}_{\operatorname{\mathcal{F}}}\left(a^{\sharp_{1}}(f_{1})\cdots
a^{\sharp_{2n+1}}(f_{2n+1})\boldsymbol{\omega}\right)=0,$
$\displaystyle\small\operatorname{Tr}_{\operatorname{\mathcal{F}}}\left(a^{\sharp_{1}}(f_{1})\cdots
a^{\sharp_{2n}}(f_{2n})\boldsymbol{\omega}\right)=\sum_{\sigma}(-1)^{\sigma}\prod^{n}_{j=1}\operatorname{Tr}\\!\left(a^{\sharp_{\sigma(j)}}(f_{\sigma(j)})a^{\sharp_{\sigma(j+n)}}(f_{\sigma(j+n)})\boldsymbol{\omega}\right),$

where $a^{\sharp_{j}}$ stands for either $a$ or $a^{*}$, and the sum is taken
over all permutations $\sigma$ satisfying
$\sigma(1)<\sigma(2)<\ldots<\sigma(n)$, and $\sigma(j)<\sigma(j+n)$, for all
$j\in\set{1,\ldots,n}$. The definition indicates that quasi-free states are
determined by the one-particle reduced density matrix operator defined via its
integral kernel

(12)
$\displaystyle\boldsymbol{\rho}(x,y)=\frac{1}{Nh^{d}}\operatorname{Tr}_{\operatorname{\mathcal{F}}}\left(a^{\ast}_{y}\,a_{x}\,\boldsymbol{\omega}\right),$

and the antisymmetric pairing function defined by

(13)
$\displaystyle\boldsymbol{\alpha}(x,y)=\frac{1}{Nh^{d}}\operatorname{Tr}_{\operatorname{\mathcal{F}}}\left(a_{y}\,a_{x}\,\boldsymbol{\omega}\right).$

More compactly, we introduce the generalized one-particle density matrix
$\boldsymbol{\Gamma}:\mathfrak{h}\oplus\mathfrak{h}^{\ast}\rightarrow\mathfrak{h}\oplus\mathfrak{h}^{\ast}$
defined by

(14)
$\displaystyle\boldsymbol{\Gamma}=\begin{bmatrix}\boldsymbol{\rho}&\boldsymbol{\alpha}\\\
-\mathds{J}\boldsymbol{\alpha}\mathds{J}&\mathbf{1}-\mathds{J}\boldsymbol{\rho}\mathds{J}^{\ast}\end{bmatrix}$

where $\mathds{J}:\mathfrak{h}\rightarrow\mathfrak{h}^{\ast}$ is the map
$\mathds{J}(\phi)=\left\langle\phi,\cdot\right\rangle$ and $\mathds{J}^{\ast}$
is its adjoint operator. This observation is summarized by the following
proposition (see e.g. [15, Appendix G]).

######  Proposition 1.1.

Let
$\boldsymbol{\Gamma}:\mathfrak{h}\oplus\mathfrak{h}^{\ast}\rightarrow\mathfrak{h}\oplus\mathfrak{h}^{\ast}$
be an operator of the form (14). Then $\boldsymbol{\Gamma}$ is the generalized
one-particle density matrix of some quasi-free state with finite particle
number if and only if $\boldsymbol{\Gamma}\geq 0$ and
$\operatorname{Tr}\boldsymbol{\rho}<\infty$.

In the case of $\boldsymbol{\omega}_{N}$, it can be readily shown that its
generalized one-particle density matrix is given by

$\boldsymbol{\Gamma}_{\beta}=\begin{bmatrix}\boldsymbol{\rho}_{\beta}&0\\\
0&\mathbf{1}-\mathds{J}\boldsymbol{\rho}_{\beta}\mathds{J}^{\ast}\end{bmatrix}$

where $\boldsymbol{\alpha}=0$, which we refer to as the gauge-invariant
condition. The following proposition, whose proof can be found for example in
[4, Proposition 5.2.23], makes the link between the Fock space Gibbs state
(11) and the one body Gibbs state considered in our main theorem.

######  Proposition 1.2.

Let $\boldsymbol{\omega}_{N}$ denotes the Gibbs state (11) with
$\operatorname{Tr}\\!\left(e^{-\beta\mathsf{H}}\right)<\infty$. Then
$\boldsymbol{\omega}_{N}$ is quasi-free and its one-particle reduced density
matrix operator is given by Equation (4).

##  2\. Semiclassical thermal regularity

In this section, we give the proof of the main theorem.

###  2.1. Preliminaries

For a trapped non-interacting Fermi gas with Hamiltonian (8), let us write
$G=e^{-\beta H}$ and $\lambda=Nh^{d}$. Consider the quantities

(15a) $\displaystyle\boldsymbol{g}_{\beta}$
$\displaystyle=Z_{\beta}^{-1}\,e^{-\beta H},$ (15b)
$\displaystyle\boldsymbol{\rho}_{\beta}$
$\displaystyle=\lambda^{-1}\left(\mathbf{1}+e^{\beta(H-\mu)}\right)^{-1}=G_{\mu}\left(\mathbf{1}+\lambda\,G_{\mu}\right)^{-1},$

where $G_{\mu}=\lambda^{-1}\,e^{-\beta\left(H-\mu\right)}=Z_{\mu}^{-1}\,G$ and
$Z_{\beta}$ and $Z_{\mu}$ are given in Equation (5). The quantum gradients of
these operators are given by the following lemma.

######  Lemma 2.1.

Let $\boldsymbol{\nabla}$ be a quantum gradient. Then we have

(16a) $\displaystyle\boldsymbol{\nabla}\boldsymbol{g}_{\beta}$
$\displaystyle=-\frac{\beta}{Z_{\beta}}\int_{0}^{1}G^{1-s}\left(\boldsymbol{\nabla}H\right)G^{s}\,\mathrm{d}s,$
(16b) $\displaystyle\boldsymbol{\nabla}\boldsymbol{\rho}_{\beta}$
$\displaystyle=-\beta\int_{0}^{1}G_{\mu}^{1-s}(\mathbf{1}+\lambda
G_{\mu})^{-1}\left(\boldsymbol{\nabla}H\right)(\mathbf{1}+\lambda
G_{\mu})^{-1}G_{\mu}^{s}\,\mathrm{d}s.$

###### Proof.

Notice that for any operators $A$ and $B$, we have that

$\partial_{t}\left[{A,e^{tB}}\right]=\left[{A,B\,e^{tB}}\right]=\left[{A,B}\right]e^{tB}+B\,\left[{A,e^{tB}}\right].$

Therefore,
$\partial_{t}\left(e^{-tB}\left[{A,e^{tB}}\right]\right)=e^{-tB}\left[{A,B}\right]e^{tB}$
and we obtain the Duhamel-like formula

(17)
$\left[{A,e^{tB}}\right]=\int_{0}^{t}e^{\left(t-s\right)B}\left[{A,B}\right]e^{sB}\,\mathrm{d}s.$

In particular, taking $t=1$ and $B=-\beta H$, we obtain the identity

$\boldsymbol{\nabla}G=-\beta\int_{0}^{1}G^{1-s}\left(\boldsymbol{\nabla}H\right)G^{s}\,\mathrm{d}s$

from which we deduce Identity (16a). Then observe that if $A$ is an invertible
operator, by the Leibniz rule for commutators,
$0=\boldsymbol{\nabla}\\!\left(AA^{-1}\right)=\left(\boldsymbol{\nabla}A\right)A^{-1}+A\boldsymbol{\nabla}\\!\left(A^{-1}\right)$
and so

(18)
$\boldsymbol{\nabla}\\!\left(A^{-1}\right)=-A^{-1}\left(\boldsymbol{\nabla}A\right)A^{-1}.$

In particular, since $\boldsymbol{\nabla}\left(\mathbf{1}+\lambda
G_{\mu}\right)=\lambda\boldsymbol{\nabla}G_{\mu}$, we deduce that

$\displaystyle\boldsymbol{\nabla}\boldsymbol{\rho}_{\beta}$
$\displaystyle=\left(\boldsymbol{\nabla}G_{\mu}\right)\left(\mathbf{1}+\lambda
G_{\mu}\right)^{-1}-G_{\mu}\left(\mathbf{1}+\lambda
G_{\mu}\right)^{-1}\lambda\left(\boldsymbol{\nabla}G_{\mu}\right)\left(\mathbf{1}+\lambda
G_{\mu}\right)^{-1}$ $\displaystyle=\left(\mathbf{1}+\lambda
G_{\mu}\right)^{-1}\left(\boldsymbol{\nabla}G_{\mu}\right)\left(\mathbf{1}+\lambda
G_{\mu}\right)^{-1}$

which leads to Formula (16b). ∎

Taking the square root of $G$ changes $\beta$ by $\beta/2$. Hence it follows
that

$\boldsymbol{\nabla}\sqrt{\boldsymbol{g}_{\beta}}=-\frac{\beta}{2\sqrt{Z_{\beta}}}\int_{0}^{1}G^{(1-s)/2}\left(\boldsymbol{\nabla}H\right)G^{s/2}\,\mathrm{d}s.$

On the other hand, it is more difficult to compute explicitly
$\boldsymbol{\nabla}\sqrt{\boldsymbol{\rho}_{\beta}}$. The following lemma
allows us to bound $\sqrt{\boldsymbol{\rho}_{\beta}}$ by reducing the problem
to estimating $G$.

######  Lemma 2.2.

Let $\mathsf{m}$ be a self-adjoint operator. Then

$\left\lVert\boldsymbol{\nabla}\sqrt{\boldsymbol{\rho}_{\beta}}\right\rVert_{\mathcal{L}^{p}(\mathsf{m})}\leq\left(\left\lVert\boldsymbol{\nabla}\sqrt{G_{\mu}}\right\rVert_{\mathcal{L}^{p}(\mathsf{m})}+\frac{\lambda}{2}\left\lVert\sqrt{\boldsymbol{\rho}_{\beta}}\right\rVert_{\mathcal{L}^{q}(\mathsf{m})}\left\lVert\boldsymbol{\nabla}G_{\mu}\right\rVert_{\mathcal{L}^{r}}\right)$

for any $1\leq p,q,r\leq\infty$ such that
$\frac{1}{p}=\frac{1}{q}+\frac{1}{r}$.

###### Proof.

By Identity (18), we have that

(19)
$\displaystyle\boldsymbol{\nabla}\sqrt{\boldsymbol{\rho}_{\beta}}=\left(\boldsymbol{\nabla}\sqrt{G_{\mu}}-\sqrt{\boldsymbol{\rho}_{\beta}}\left(\boldsymbol{\nabla}\sqrt{\mathbf{1}+\lambda
G_{\mu}}\right)\right)\left(\mathbf{1}+\lambda G_{\mu}\right)^{-\frac{1}{2}}.$

To bound $\boldsymbol{\nabla}\sqrt{\mathbf{1}+\lambda G_{\mu}}$, we proceed as
in [6, Lemma 7.1]. Let $A:=\lambda G_{\mu}$. Since $0\leq A\leq 2\,c:=\lambda
Z_{\mu}^{-1}$, we deduce that $\left\lVert
A-c\right\rVert_{\mathcal{L}^{\infty}}\leq c$ and so the following series is
absolutely convergent

$\sqrt{\mathbf{1}+A}=\sqrt{1+c}\,\sqrt{\mathbf{1}+\tfrac{1}{1+c}\left(A-c\right)}=\sum_{n=0}^{\infty}\binom{1/2}{n}\left(\tfrac{1}{1+c}\right)^{n-\frac{1}{2}}\left(A-c\right)^{n}.$

Since it follows from the Jacobi identity that

$\boldsymbol{\nabla}\left(A-c\right)^{n}=\sum_{k=1}^{n}\left(A-c\right)^{k-1}\left(\boldsymbol{\nabla}A\right)\left(A-c\right)^{n-k},$

then we deduce the estimate

$\displaystyle\left\lVert\boldsymbol{\nabla}\sqrt{\mathbf{1}+A}\right\rVert_{\mathcal{L}^{p}}$
$\displaystyle\leq\sum_{n=1}^{\infty}\left\lvert\binom{1/2}{n}\right\rvert\left(\tfrac{1}{c+1}\right)^{n-\frac{1}{2}}n\left\lVert\boldsymbol{\nabla}A\right\rVert_{\mathcal{L}^{p}}\left\lVert
A-c\right\rVert_{\mathcal{L}^{\infty}}^{n-1}$
$\displaystyle\leq\frac{1}{2\sqrt{c+1}}\,\sum_{n=1}^{\infty}\binom{-1/2}{n-1}\left(\tfrac{-c}{c+1}\right)^{n-1}\left\lVert\boldsymbol{\nabla}A\right\rVert_{\mathcal{L}^{p}}=\frac{1}{2}\left\lVert\boldsymbol{\nabla}A\right\rVert_{\mathcal{L}^{p}}.$

We conclude by applying the fact that
$\left\lVert\boldsymbol{\nabla}\sqrt{\boldsymbol{\rho}_{\beta}}\,\mathsf{m}\right\rVert_{\mathcal{L}^{p}}=\left\lVert\mathsf{m}\,\boldsymbol{\nabla}\sqrt{\boldsymbol{\rho}_{\beta}}\right\rVert_{\mathcal{L}^{p}}$,
which follows from taking the adjoint, and the triangle inequality to Formula
(19). This yields the desired result. ∎

###  2.2. Bounds on the inverse fugacity

In the rest of this paper, we assume $H$ is the harmonic oscillator
Hamiltonian $H=\frac{\left\lvert\boldsymbol{p}\right\rvert^{2}+\left\lvert
x\right\rvert^{2}}{2}$. In this case of the partition function in
$\boldsymbol{g}_{\beta}$ has the closed form

(20)
$Z_{\beta}=\left(\frac{2\pi}{\beta}\right)^{d}\frac{1}{\operatorname{shc}(\frac{\beta\hbar}{2})^{d}}$

where $\operatorname{shc}(x)=\operatorname{sh}(x)/x$ denotes the hyperbolic
sinc function. In particular $Z_{\beta}\leq(2\pi/\beta)^{d}$ and
$Z_{\beta}\sim(2\pi/\beta)^{d}$ as $\beta\hbar\rightarrow 0$. The function
$Z_{\mu}$, that one could call the inverse fugacity, can be compared to the
partition function $Z_{\beta}$ as proved in the following proposition.

######  Proposition 2.1.

Let $Z_{\mu}=\lambda\,e^{-\beta\mu}$ with $\lambda=Nh^{d}$. Then the following
inequality holds

(21) $C_{\lambda,\beta}^{-1}\,Z_{\beta}\leq Z_{\mu}\leq Z_{\beta}$

with $C_{\lambda,\beta}=2$ if $\mu\leq d\hbar/2$ and
$C_{\lambda,\beta}=1+e^{\beta\lambda^{1/d}/2\pi}$ if $\mu\geq d\hbar/2$.

###### Proof of Proposition 2.1.

Let $c_{\mu}=e^{\beta\mu}=Nh^{d}/Z_{\mu}$. Notice that

$1=h^{d}\operatorname{Tr}\\!\left(\boldsymbol{\rho}_{\beta}\right)=\frac{1}{N}\operatorname{Tr}\\!\left(\big{(}\mathbf{1}+c_{\mu}^{-1}e^{\beta
H}\big{)}^{-1}\right)\leq\frac{c_{\mu}}{N}\operatorname{Tr}\\!\left(e^{-\beta
H}\right),$

which implies the first inequality in Formula (21).

Next, define the function $g:\mathbb{R}_{+}\to\mathbb{R}_{+}$ by
$g(r)=\big{(}1+c_{\mu}^{-1}\,e^{\beta r}\big{)}^{-1}$, then for any $R>0$, it
holds $g(r)\geq\big{(}1+c_{\mu}^{-1}\,e^{\beta
r}\big{)}^{-1}\,\mathds{1}_{r\leq R}$, which leads to

$1=\frac{1}{N}\operatorname{Tr}\left(g(H)\right)\geq\frac{1}{N\left(1+c_{\mu}^{-1}\,e^{\beta
R}\right)}\,\operatorname{Tr}\left(\mathds{1}_{H\leq R}\right).$

Recalling the property of the harmonic oscillator given in Section 1.2, one
sees that the trace of the characteristic function is nothing but

$\operatorname{Tr}\\!\left(\mathds{1}_{H\leq
R}\right)=\left\lvert\Set{n\in\mathbb{N}^{d}_{0}:\left(\left\lvert
n\right\rvert_{1}+\tfrac{d}{2}\right)\hbar\leq R}\right\rvert.$

Since $\left\lvert n\right\rvert_{1}\geq\sup_{j}n_{j}=:\left\lvert
n\right\rvert_{\infty}$, then this can be crudely estimated by

(22) $\operatorname{Tr}\\!\left(\mathds{1}_{H\leq
R}\right)\geq\left\lvert\Set{n\in\mathbb{N}^{d}_{0}:\left\lvert
n\right\rvert_{\infty}\leq\tfrac{R}{\hbar}-\tfrac{d}{2}}\right\rvert={\left(\left\lfloor{\tfrac{R}{\hbar}-\tfrac{d}{2}}\right\rfloor+1\right)^{d}}.$

In particular, since $\mu\geq d\hbar/2>0$, then taking $R=\mu=\frac{\ln
c_{\mu}}{\beta}$ yields
$\tfrac{\mu}{\hbar}-\frac{d}{2}\leq\big{(}N^{1/d}+1\big{)}$ and so since
$N\geq 1$,

(23) $\mu\leq 2\,N^{1/d}\hbar+\frac{d\hbar}{2}.$

Now, let us obtain upper bounds for $c_{\mu}$. First, observe that we have the
lower bounds

(24)
$1=\frac{c_{\mu}}{N}\operatorname{Tr}\\!\left(\big{(}\mathbf{1}+c_{\mu}e^{-\beta
H}\big{)}^{-1}e^{-\beta
H}\right)\geq\frac{Z_{\beta}\,c_{\mu}}{Nh^{d}\left(1+c_{\mu}e^{-d\beta\hbar/2}\right)}.$

If $c_{\mu}\leq e^{d\beta\hbar/2}$ (i.e. $\mu\leq d\hbar/2$) then it follows
from Inequality (24) that

$\displaystyle c_{\mu}\leq 2\,Nh^{d}Z_{\beta}^{-1}.$

On the other hand, if $c_{\mu}\geq e^{d\beta\hbar/2}>1$ (i.e. $\mu\geq
d\hbar/2$), then it follows from Inequality (23) and Inequality (24) that we
have the bound

(25) $\displaystyle
c_{\mu}\leq\left(1+c_{\mu}\,e^{-d\beta\hbar/2}\right)Nh^{d}Z_{\beta}^{-1}\leq\left(1+e^{2\beta\hbar
N^{1/d}}\right)Nh^{d}Z_{\beta}^{-1}$

This completes our proof of the proposition. ∎

###  2.3. $\mathcal{L}^{\infty}$ bounds

This section is devoted to the proof of the following proposition.

######  Proposition 2.2.

Let $\beta>0$ and
$H=\frac{\left\lvert\boldsymbol{p}\right\rvert^{2}+\left\lvert
x\right\rvert^{2}}{2}$, then we have the estimates

$\displaystyle\left\lVert\boldsymbol{\nabla}\boldsymbol{g}_{\beta}\right\rVert_{\mathcal{L}^{\infty}}$
$\displaystyle\leq\frac{2}{Z_{\beta}}\max\left(\sqrt{\beta},\beta\sqrt{\hbar}\right),$
$\displaystyle\left\lVert\boldsymbol{\nabla}\boldsymbol{\rho}_{\beta}\right\rVert_{\mathcal{L}^{\infty}}$
$\displaystyle\leq\frac{2}{Z_{\mu}}\max\left(\sqrt{\beta},\beta\sqrt{\hbar}\right).$

Prior to giving the proof of the above proposition, let us make the following
observation. Since $\boldsymbol{\nabla}_{\\!x}H=x$ and
$\boldsymbol{\nabla}_{\\!\xi}H=\boldsymbol{p}$, we deduce that
$\boldsymbol{\nabla}_{z}H=\boldsymbol{z}:=(x,\boldsymbol{p})$ where
$z=(x,\xi)$. By Identity (16a), we have that

$\boldsymbol{\nabla}\boldsymbol{\rho}_{\beta}=-\beta\int_{0}^{1}(\mathbf{1}+\lambda
G_{\mu})^{-1}G_{\mu}^{1-s}\,\boldsymbol{z}\,G_{\mu}^{s}(\mathbf{1}+\lambda
G_{\mu})^{-1}\,\mathrm{d}s$

which implies the estimate

$\left\lVert\boldsymbol{\nabla}\boldsymbol{\rho}_{\beta}\right\rVert_{\mathcal{L}^{p}}\leq\beta\,\int_{0}^{1}\left\lVert
G_{\mu}^{1-s}\boldsymbol{z}\,G_{\mu}^{s}\right\rVert_{\mathcal{L}^{p}}\mathrm{d}s=\beta\,Z_{\mu}^{-1}\,\int_{0}^{1}\left\lVert
G^{1-s}\boldsymbol{z}\,G^{s}\right\rVert_{\mathcal{L}^{p}}\mathrm{d}s.$

Since $\left\lVert G\right\rVert_{\mathcal{L}^{\infty}}\leq 1$, then we have
the estimate

(26)
$\displaystyle\left\lVert\boldsymbol{\nabla}\boldsymbol{\rho}_{\beta}\right\rVert_{\mathcal{L}^{p}}$
$\displaystyle\leq\beta\,Z_{\mu}^{-1}\int_{0}^{1/2}\left\lVert
G^{1-s}\boldsymbol{z}\right\rVert_{\mathcal{L}^{p}}\mathrm{d}s+\beta\,Z_{\mu}^{-1}\int_{1/2}^{1}\left\lVert\boldsymbol{z}\,G^{s}\right\rVert_{\mathcal{L}^{p}}\mathrm{d}s$
$\displaystyle\leq
2\,\beta\,Z_{\mu}^{-1}\int_{1/2}^{1}\left\lVert\boldsymbol{z}\,G^{s}\right\rVert_{\mathcal{L}^{p}}\mathrm{d}s.$

Hence to estimate $\boldsymbol{\nabla}\boldsymbol{\rho}_{\beta}$ and
$\boldsymbol{\nabla}\boldsymbol{g}_{\beta}$, it remains to estimate the value
of $\left\lVert x\,e^{-\beta sH}\right\rVert_{\mathcal{L}^{p}}$ and
$\left\lVert\boldsymbol{p}\,e^{-\beta sH}\right\rVert_{\mathcal{L}^{p}}$ for
$s\in[\frac{1}{2},1]$. Let start with the case $p=\infty$.

######  Lemma 2.3.

Let $\beta>0$, then we have the estimate

$\displaystyle\left\lVert\left\lvert x\right\rvert^{n}e^{-\beta
H}\right\rVert_{\mathcal{L}^{\infty}}^{2/n}=\left\lVert\left\lvert\boldsymbol{p}\right\rvert^{n}e^{-\beta
H}\right\rVert_{\mathcal{L}^{\infty}}^{2/n}\leq
n\max\left(\tfrac{2}{\beta},\sqrt{2}\,\hbar\right).$

######  Remark 2.1.

In the classical case, it is not difficult to prove that the maximum of the
function $x\mapsto\left\lvert x\right\rvert^{n}e^{-\beta\left\lvert
x\right\rvert^{2}}$ is $\left(\frac{n}{2e\beta}\right)^{n/2}$, and more
generally,

$\left\lVert\left\lvert x\right\rvert^{n}e^{-\beta\left\lvert
x\right\rvert^{2}}\right\rVert_{L^{p}}^{p}=\omega_{d}\,\Gamma(\tfrac{d+np}{2})\left(\beta
p\right)^{-\frac{d+np}{2}}$

where $\omega_{d}$ is the volume of the $d$-dimensional unit ball and $\Gamma$
is the gamma function.

###### Proof of Lemma 2.3.

It is sufficient to prove a bound on the first quantity because $H$ is
symmetric in $x$ and $\boldsymbol{p}$. Let $\varphi\in L^{2}$ and
$\psi=e^{-tH}\varphi$. Notice that for any $t\geq 0$,
$\left\lVert\psi\right\rVert_{L^{2}}\leq\left\lVert\varphi\right\rVert_{L^{2}}$.
Let $y:=\left\lVert\left\lvert
x\right\rvert^{n}\psi\right\rVert_{L^{2}}^{2/n}$. Since
$2\,\partial_{t}\psi=-\big{(}\left\lvert
x\right\rvert^{2}-\hbar^{2}\Delta\big{)}\psi$, then integrating by parts
yields

$\frac{\mathrm{d}}{\mathrm{d}t}y^{n}=-\int_{\mathbb{R}^{d}}\left\lvert\psi\right\rvert^{2}\left\lvert
x\right\rvert^{2\left(n+1\right)}\mathrm{d}x-\hbar^{2}\operatorname{Re}\\!\left(\int_{\mathbb{R}^{d}}\nabla\left(\overline{\psi}\left\lvert
x\right\rvert^{2n}\right)\cdot\nabla\psi\,\mathrm{d}x\right)\\\
\leq-\int_{\mathbb{R}^{d}}\left\lvert\psi\right\rvert^{2}\left\lvert
x\right\rvert^{2\left(n+1\right)}+\hbar^{2}\left\lvert\nabla\psi\right\rvert^{2}\left\lvert
x\right\rvert^{2n}\mathrm{d}x+2\,n\,\hbar^{2}\int_{\mathbb{R}^{d}}\left\lvert\psi\right\rvert\left\lvert
x\right\rvert^{2n-1}\left\lvert\nabla\psi\right\rvert\,\mathrm{d}x.$

Applying Young’s inequality for the product, we get

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}y^{n}$
$\displaystyle\leq-\int_{\mathbb{R}^{d}}\left\lvert\psi\right\rvert^{2}\left\lvert
x\right\rvert^{2(n+1)}\mathrm{d}x+\left(n\hbar\right)^{2}\int_{\mathbb{R}^{d}}\left\lvert\psi\right\rvert^{2}\left\lvert
x\right\rvert^{2\left(n-1\right)}\mathrm{d}x$
$\displaystyle\leq-c^{-1}\,y^{n+1}+\left(n\hbar\right)^{2}c\,y^{n-1}$

where $c=\left\lVert\varphi\right\rVert_{L^{2}}^{2/n}$. This yields the
differential inequality

(27)
$y^{\prime}\leq-\frac{1}{n\,c}\left(y^{2}-\left(n\hbar\,c\right)^{2}\right).$

This ordinary differential equation has a fixed point at $y=n\hbar\,c$. If
initially, $y\leq\hbar\,c$, then $y^{\prime}\geq 0$ but $y$ remains smaller
than $n\hbar\,c$. If not, then at any time $y^{\prime}<0$ and
$\hbar\,c<y(t)<y(0)$. If initially $y>\sqrt{2}\,n\,\hbar\,c$, then as long as
it remains true, it holds

$y^{\prime}\leq-\frac{1}{2\,n\,c}\,y^{2}$

which implies

$y(t)\leq\max\left(\frac{2}{t},\sqrt{2}\,\hbar\right)n\,c$

and proves the result by taking $t=\beta$. ∎

To complete the proof of Proposition 2.2, we use the above lemma to get that
for any $s\in[\tfrac{1}{2},1]$,

$\left\lVert x\,e^{-\beta sH}\right\rVert_{\mathcal{L}^{\infty}}\leq
2\max\left(\tfrac{1}{\sqrt{\beta}},\sqrt{\hbar}\right).$

Finally, we conclude using Inequality (26).

###  2.4. $\mathcal{L}^{p}$ bounds for $2\leq p<\infty$

######  Proposition 2.3.

Let $\beta>0$ and $\hbar\in(0,1)$. Suppose $p\in[2,\infty]$ then there exists
$C_{d,p}>0$ such that

$\displaystyle\left\lVert\boldsymbol{\nabla}\boldsymbol{g}_{\beta}\right\rVert_{\mathcal{L}^{p}}$
$\displaystyle
C_{d,p}\,\frac{\beta^{\frac{1}{2}-\frac{d}{p}}}{Z_{\beta}}\frac{\max\left(2\sqrt{2},\beta\hbar\right)^{\frac{1}{2}-\frac{1}{p}}}{\theta(\beta\hbar)^{\frac{1}{p}}},$
$\displaystyle\left\lVert\boldsymbol{\nabla}\boldsymbol{\rho}_{\beta}\right\rVert_{\mathcal{L}^{p}}$
$\displaystyle\leq
C_{d,p}\,\frac{\beta^{\frac{1}{2}-\frac{d}{p}}}{Z_{\mu}}\frac{\max\left(2\sqrt{2},\beta\hbar\right)^{\frac{1}{2}-\frac{1}{p}}}{\theta(\beta\hbar)^{\frac{1}{p}}}.$

where where $\theta(x)=\operatorname{th}(x)/x$ with
$\operatorname{th}(x)=\tfrac{e^{x}-e^{-x}}{e^{x}+e^{-x}}$.

######  Lemma 2.4 (Moments of the thermal state).

Let $\beta>0$ and $n>-d$. Then

$\displaystyle h^{d}\operatorname{Tr}\\!\left(\left\lvert
x\right\rvert^{n}\boldsymbol{g}_{\beta}\right)=h^{d}\operatorname{Tr}\\!\left(\left\lvert\boldsymbol{p}\right\rvert^{n}\boldsymbol{g}_{\beta}\right)=\frac{C_{d,n}}{\left(\beta\,\theta(\beta\hbar/2)/2\right)^{n/2}}$

where $C_{d,n}=\Gamma(\tfrac{d+n}{2})/\Gamma(\tfrac{d}{2})$.

###### Proof.

Let $H_{\circ}=\frac{{x}^{2}+\left\lvert\boldsymbol{p}\right\rvert^{2}}{2}$ be
the one-dimensional Hamiltonian of the harmonic oscillator, and $\psi_{n}$ be
its eigenvalues verifying

$H_{\circ}\psi_{n}=\left(n+\tfrac{1}{2}\right)\hbar\,\psi_{n}.$

The Wigner transform of the corresponding density operator
$\ket{\psi_{n}}\bra{\psi_{n}}$, is classically given by (see e.g. [12, Section
5.04] or [7, Theorem 1.105])

$f_{n}(z)=2\left(-1\right)^{n}e^{-\left\lvert
z\right\rvert^{2}/\hbar}\,L_{n}\\!\left(\frac{2\left\lvert
z\right\rvert^{2}}{\hbar}\right)$

where $z=(x,\xi)$ and $L_{n}$ is the Laguerre polynomial of order $n$ defined
by

$L_{n}(z)=\frac{e^{x}}{n!}\,\partial_{x}^{n}\\!\left(x^{n}e^{-x}\right)=\sum_{k=0}^{n}\binom{n}{k}\frac{(-x)^{k}}{k!}$

By the formula of the generating function of the Laguerre polynomials, we
deduce that

$\sum_{n=0}^{\infty}t^{n}f_{n}(z)=\frac{2}{1+t}\,e^{-\frac{\left\lvert
z\right\rvert^{2}}{\hbar}\frac{1-t}{1+t}}$

Taking $t=e^{-\beta\hbar}$, we obtain the Wigner transform of the thermal
state. Since in dimension $d$ it is factorized, it yields to the following
formula for the Wigner transform of $Z_{\beta}^{-1}\,e^{-\beta H}$

$f_{\beta}(z)=\left(\frac{\beta\,\theta(\beta\hbar/2)}{2\pi}\right)^{d}e^{-\beta\left\lvert
z\right\rvert^{2}\theta(\beta\hbar/2)}$

Its spatial moments are given by

$\displaystyle h^{d}\operatorname{Tr}\\!\left(\left\lvert
x\right\rvert^{n}\boldsymbol{g}_{\beta}\right)$
$\displaystyle=\iint_{\mathbb{R}^{2d}}f_{\beta}(z)\left\lvert
x\right\rvert^{n}\mathrm{d}z=\left(\frac{\beta\,\theta(\beta\hbar/2)}{2\pi}\right)^{d/2}\int_{\mathbb{R}^{d}}\left\lvert
x\right\rvert^{n}e^{-\beta\left\lvert
x\right\rvert^{2}\theta(\beta\hbar/2)}\,\mathrm{d}x$

which yields the result. ∎

###### Proof of Proposition 2.3.

Since $h^{d}\operatorname{Tr}\\!\left(\left\lvert\left\lvert
x\right\rvert^{n}G\right\rvert^{2}\right)=Z_{2\beta}\,h^{d}\operatorname{Tr}\\!\left(\left\lvert
x\right\rvert^{2n}\boldsymbol{g}_{2\beta}\right)$, then, by Lemma 2.4, we have
the identity

$\left\lVert\left\lvert
x\right\rvert^{n}G\right\rVert_{\mathcal{L}^{2}}^{2}=\left\lVert\left\lvert\boldsymbol{p}\right\rvert^{n}G\right\rVert_{\mathcal{L}^{2}}^{2}=\frac{C_{d,2n}\,Z_{2\beta}}{\left(\beta\,\theta(\beta\hbar)\right)^{n}}$

whenever $n>-d/2$. Now applying linear interpolation of Schatten norms and
Lemma 2.3, we obtain the intermediate Schatten norm bounds

$\left\lVert\left\lvert
x\right\rvert^{n}G\right\rVert_{\mathcal{L}^{p}}=\left\lVert\left\lvert\boldsymbol{p}\right\rvert^{n}G\right\rVert_{\mathcal{L}^{p}}\leq
C_{d,n,p}\,Z_{2\beta}^{\frac{1}{p}}\,\frac{\max\left(\tfrac{2}{\beta},\sqrt{2}\,\hbar\right)^{n\left(\frac{1}{2}-\frac{1}{p}\right)}}{\left(\beta\,\theta(\beta\hbar)\right)^{\frac{n}{p}}}$

where
$C_{d,n,p}=\,C_{d,2n}^{\frac{1}{p}}\,n^{n\left(\frac{1}{2}-\frac{1}{p}\right)}$.
Finally, by Identity (16a), we arrive at

$\left\lVert\boldsymbol{\nabla}\boldsymbol{g}_{\beta}\right\rVert_{\mathcal{L}^{p}}\leq\frac{2\beta}{Z_{\beta}}\int_{1/2}^{1}\left\lVert\boldsymbol{z}\,G^{s}\right\rVert_{\mathcal{L}^{p}}\mathrm{d}s\\\
\leq
C_{d,1,p}\frac{2\beta}{Z_{\beta}}\sup_{s\in[\frac{1}{2},1]}\,Z_{2s\beta}^{\frac{1}{p}}\,\frac{\max\left(\tfrac{2}{s\beta},\sqrt{2}\,\hbar\right)^{\frac{1}{2}-\frac{1}{p}}}{\left(s\beta\,\theta(s\beta\hbar)\right)^{\frac{1}{p}}}\\\
\leq
C_{d,p}\,\frac{\beta^{\frac{1}{2}-\frac{d}{p}}}{Z_{\beta}}\frac{\max\left(2\sqrt{2},\beta\hbar\right)^{\frac{1}{2}-\frac{1}{p}}}{\left(\theta(\beta\hbar)\right)^{\frac{1}{p}}}$

where we used the fact that $Z_{\beta}\leq(2\pi/\beta)^{d}$ and
$C_{d,p}=2^{\frac{5}{4}+\frac{2d+1}{p}}\,C_{d,1,p}\,\pi^{\frac{d}{p}}$.
Similarly, Inequality (26) implies the same bound for
$\boldsymbol{\rho}_{\beta}$ with $Z_{\beta}$ is replaced by $Z_{\mu}$. ∎

Acknowledgements. J.C. was supported by the NSF through the RTG grant DMS- RTG
184031. C.S. acknowledges the support of the Swiss National Science Foundation
through the Eccellenza project PCEFP2_181153 and of the NCCR SwissMAP.

## References

  * [1] N. Benedikter.  Effective Dynamics of Interacting Fermions from Semiclassical Theory to the Random Phase Approximation.  Journal of Mathematical Physics, 63(8):081101, Aug. 2022.
  * [2] N. Benedikter, V. Jakšić, M. Porta, C. Saffirio, and B. Schlein.  Mean-Field Evolution of Fermionic Mixed States.  Communications on Pure and Applied Mathematics, 69(12):2250–2303, Dec. 2016.
  * [3] N. Benedikter, M. Porta, and B. Schlein.  Mean-field Evolution of Fermionic Systems.  Communications in Mathematical Physics, 331(3):1087–1131, Nov. 2014\.
  * [4] O. Bratteli and D. W. Robinson.  Operator Algebras and Quantum Statistical Mechanics II.  Springer-Verlag, Berlin, Heidelberg, 1981.
  * [5] S. Chapman and T. G. Cowling.  The Mathematical Theory of Non-uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases.  Cambridge Mathematical Library Series. Cambridge University Press, Cambridge ; New York, 3rd edition edition, 1990.
  * [6] J. J. W. Chong, L. Lafleche, and C. Saffirio.  From Many-Body Quantum Dynamics to the Hartree–Fock and Vlasov Equations with Singular Potentials.  arXiv:2103.10946, pages 1–74, Mar. 2021.
  * [7] G. B. Folland.  Harmonic Analysis in Phase Space, volume 122 of The Annals of Mathematics Studies.  Princeton University Press, Princeton, N.J., Mar. 1989.
  * [8] S. Fournais and S. Mikkelsen.  An optimal semiclassical bound on certain commutators.  arXiv:1912.08467, pages 1–24, Dec. 2019.
  * [9] F. Golse, S. Jin, and T. Paul.  On the Convergence of Time Splitting Methods for Quantum Dynamics in the Semiclassical Regime.  Foundations of Computational Mathematics, 21(3):613–647, June 2021\.
  * [10] F. Golse and T. Paul.  Empirical Measures and Quantum Mechanics: Application to the Mean-Field Limit.  Communications in Mathematical Physics, 369(3):1021–1053, Aug. 2019\.
  * [11] F. Golse and T. Paul.  Semiclassical Evolution With Low Regularity.  arXiv:2011.14884 [math], pages 1–55, Nov. 2020.
  * [12] H. J. Groenewold.  On the principles of elementary quantum mechanics.  Physica, 12(7):405–460, Oct. 1946.
  * [13] L. Lafleche and C. Saffirio.  Strong Semiclassical Limit from Hartree and Hartree–Fock to Vlasov–Poisson Equation.  Analysis & PDE, to appear:1–35, Oct. 2021.
  * [14] M. Porta, S. Rademacher, C. Saffirio, and B. Schlein.  Mean Field Evolution of Fermions with Coulomb Interaction.  Journal of Statistical Physics, 166(6):1345–1364, Mar. 2017.
  * [15] J. P. Solovej.  Many Body Quantum Mechanics.  Course notes, University of Vienna, Wien, 2014.


# Identification of time scales of the violation of the Stokes-Einstein
relation in Yukawa liquids

Zahra Ghannad<EMAIL_ADDRESS>Department of Physics, Alzahra
University, P. O. Box 19938-93973, Tehran, Iran

###### Abstract

We investigate the origin of the violation of the Stokes-Einstein (SE)
relation in two-dimensional Yukawa liquids. Using comprehensive molecular
dynamics simulations, we identify the time scales supporting the violation of
the SE relation $D\propto(\eta/T)^{-1}$, where $D$ is the self-diffusion
coefficient and $\eta$ is the shear viscosity. We first compute the self-
intermediate scattering function $F_{s}(k,t)$, the non-Gaussian parameter
$\alpha_{2}$, and the autocorrelation function of the shear stress
$C_{\eta}(t)$. The timescales obtained from these functions are included the
structural relaxation time $\tau_{\alpha}$, the peak time of the non-Gaussian
parameter $\tau_{\alpha_{2}}$, and the shear stress relaxation time
$\tau_{\eta}$. We find that $\tau_{\eta}$ is coupled with $D$ for all
temperatures indicating the SE preservation, however, $\tau_{\alpha}$ and
$\tau_{\alpha_{2}}$ are decoupled with $D$ at low temperatures indicating the
SE violation. Surprisingly, we find that the origins of this violation are
related to the non-exponential behavior of the autocorrelation function of the
shear stress and non-Gaussian behavior of the distribution function of
particle displacements. These results confirm dynamic heterogeneity that
occurs in two-dimensional Yukawa liquids that reflects the presence of regions
in which dust particles move faster than the rest when the liquid cools to
below the phase transition temperature.

††preprint: AIP/123-QED

##  I Introduction

The well-known Stokes-Einstein (SE) relation relates the diffusion coefficient
$D$ of a Brownian particle immersed in a liquid to the shear viscosity $\eta$
of a liquid at temperature $T$. For three-dimension liquids, it is given by
Einstein (1956); Hansen and McDonald (2006)

$D=\frac{k_{B}T}{c\pi\eta R}.$ (1)

Here, $R$ is the effective radius of the Brownian particle, $k_{B}$ is the
Boltzman constant, and $c$ is a constant that depends on the boundary
condition at the particle surface. For two-dimensional (2D) liquids, due to
the different dimensionality of $\eta$ in 2D, the SE relation is given by Liu,
Goree, and Vaulina (2006); Pan, Garrahan, and Chandler (2005)

$D=\frac{k_{B}T}{c^{\prime}\pi\eta}.$ (2)

Generally, the SE relation can be written as $D\eta/T=constant$, and any
deviations from it reveal a SE violation.

The violation of the Stokes-Einstein relation is a significant anomaly that
occurs in liquids Dubey  _et al._ (2019); Kawasaki and Kim (2017, 2019a,
2019b); Becker, Poole, and Starr (2006); Jeong  _et al._ (2010); Köddermann,
Ludwig, and Paschek (2008); Tsimpanogiannis  _et al._ (2019); Fernandez-Alonso
_et al._ (2007); Mallamace  _et al._ (2010); Kob  _et al._ (1997). A two-
dimensional Yukawa liquid (2DYL) is one of these liquids for which the SE
violation has been reported when the liquid is cooled to near its phase
transition temperature $T_{\rm{m}}$ Liu, Goree, and Vaulina (2006). 2DYL is a
liquid dusty plasma composed of charged dust particles that are strongly
coupled together and interact through the Yukawa potential, i.e.,
$\phi(r)=Q^{2}\rm{exp}(-\it{r}/\lambda_{D})/\rm{4}\pi\epsilon_{0}\it{r}$,
where $\lambda_{D}$ is the Debye shielding length and $Q$ is the dust charge
Yukawa (1935); Konopka, Morfill, and Ratke (2000). Due to the electric field
in the plasma sheath, dust particles can be confined and floated in a
monolayer, with an ignorable out-of-plane motion to form a 2D dusty plasma
liquid Ghannad (2019); Feng  _et al._ (2016); Wang, Huang, and Feng (2018);
Feng, lin, and Murillo (2017).

The violation of the SE relation in liquids is attributed to the dynamic
heterogeneity, that is, the presence of regions in which particles move faster
than the rest, i.e., they are more mobile than other particles. Although the
dynamic heterogeneity is cited as a commonly proposed reason for the SE
violation in literature Tarjus and Kivelson (1995); Pan  _et al._ (2017);
Sengupta  _et al._ (2013), the true origin of the violation remains unclear.
In this study, we quantify the dynamic heterogeneity in 2DYL to gain more
insight into the origin of the SE violation. We provide time scales that are
the signatures of the dynamic heterogeneity. We investigate the time of
structural relaxation, the non-Gaussian parameter, and the non-exponentiality
of the autocorrelation function of the shear stress and we show how each leads
to violation of the SE relation. This paper with identifying time scales
enhances our understanding of the SE violation mechanism and reports a
significant advance in understanding the SE violation.

The paper is organized as follows. In Section II, we describe the fundamental
features of the simulation technique to mimic 2DYL. In Section III, we compute
transport coefficients including the self-diffusion coefficient and the shear
viscosity. In Section IV, we compute correlation functions including the self-
part of the intermediate scattering function and the autocorrelation function
of the shear stress to identify time scales and interpret the physical origins
of the SE violation. In Section V, we present the conclusions.

##  II MODEL AND SIMULATION TECHNIQUE

We performed extensive molecular dynamics (MD) simulations to mimic 2DYL
Frenkel and Smit (2002). We integrated the equation of the motion
$m\ddot{\bm{\mathrm{r}}}_{i}=-\nabla\Sigma_{j}\phi_{ij}$ for $N$= 1024 dust
particles, where $m$ is the mass of the dust particle, and $\phi_{ij}$ is the
Yukawa pair interaction potential Yukawa (1935); Konopka, Morfill, and Ratke
(2000). An equilibrium Yukawa system can be characterized by two dimensionless
parameters Ohta and Hamaguchi (2000); Hartmann  _et al._ (2019). The first is
the screening parameter $\kappa=a/\lambda_{D}$. The second is the Coulomb
coupling parameter $\Gamma=Q^{2}/4\pi\epsilon_{0}ak_{B}T$, where
$\epsilon_{0}$ is the dielectric constant, $Q$ is the charge of a dust
particle, $T$ is the kinetic temperature of dust particles, $a=(n\pi)^{-1/2}$
is the Wigner-Seitz radius for 2D systems, and $n$ is the surface number
density of dust particles. We chose $\kappa=0.56$, which is a common value in
dusty plasma experiments Donkó  _et al._ (2006); Liu and Goree (2005); Liu,
Goree, and Vaulina (2006). For this value of $\kappa$, the phase transition is
at $\Gamma_{\rm{m}}\approx$ 142, so that for $\Gamma$ $<$
$\Gamma_{\rm{m}}$($T>T_{\rm{m}}$), the 2D Yukawa system is in the liquid phase
Hartmann  _et al._ (2005).

We applied normalized units in this work. The normalized temperature is
$\Gamma^{-1}$. The time is normalized by $\omega_{pd}^{-1}$, where
$\omega_{pd}=(Q^{2}/2\pi\epsilon_{0}ma^{3})^{1/2}$ is the nominal 2D dusty
plasma frequency Kalman  _et al._ (2004). The length is normalized by $a$, the
wave number by $a^{-1}$, the shear viscosity by $nm\omega_{pd}a^{2}$, and the
self-diffusion coefficient by $\omega_{pd}a^{2}$. Initially, the particles
were placed randomly into a rectangular box, and the periodic boundary
condition was applied to eliminate boundary effects due to the finite size of
the box. The size of the box is chosen so that the surface number density is
consistent with the definition of the Wigner-Seitz radius. Thus, we chose the
box with the sizes of $56.99a\times 49.08a$ so that $n\approx 1/(\pi a^{2})$.
The algorithm for integrating the equation motion is Verlet algorithm Swope
_et al._ (1982) with the integration time step of 0.037 $\omega^{-1}_{pd}$
Feng  _et al._ (2016), which is adequately small to conserve energy. A
Nos$\acute{e}$-Hoover thermostat Nosé (1984); Hoover (1985) was used to
maintain constant temperature $T$. The equilibration period was $10^{5}$ time
steps and the length of the production runs was $10^{6}$ time steps.

According to a standard test, if the simulations model a canonical ensemble in
thermal equilibrium, the ratio of the variance of the temperature obtained
from them to the variance of the temperature in the canonical ensemble,
$2\langle T\rangle^{2}/(dN)$($d$ is the dimension), must be equal to unity
Holian, Voter, and Ravelo (1995); Ghannad (2019). We reached a satisfactory
result for this ratio, i.e., 0.997.

##  III Transport properties

To investigate the Stokes-Einstein relationship, we first need to calculate
transport properties including the self-diffusion coefficient $D$ and the
shear viscosity $\eta$ in a 2DYL.

###  III.1 Self-diffusion coefficient

The self-diffusion coefficient can be computed from the linear fit of the
mean-squared displacement (MSD) versus time as Maginn  _et al._ (2019)

$D=\mathrm{\lim}_{t\to\infty}\frac{\langle|\bm{\mathrm{r}}_{i}(t)-\bm{\mathrm{r}}_{i}(0)|^{2}\rangle}{4t},$
(3)

where $\bm{\mathrm{r}}_{i}(t)$ is the position of particle $i$ at time $t$ and
$\langle...\rangle$ denotes an ensemble average. It is important to choose the
time interval to fit in a way that the motion of dust particles is diffusive
and $D$ is meaningful. Generally, the MSD obeys a power-law, i.e., MSD($t$)
$\propto t^{\mu}$, where for $\mu$ = 1, the motion is diffusive and for
$\mu\neq 1$, it is anomalous. A very short time, $\omega_{pd}t$ $\lesssim$ 5,
must be excluded from the time interval of the fit because it shows a
ballistic motion due to the trapping of dust particles in a cage created by
neighboring dust particles (Fig. 1(a)). Also, very long times must be excluded
due to the statistical error. We have selected 100 $<$ $\omega_{pd}t$ $<$ 1000
to sample $D$ and we can be sure that the motion of dust particles is
diffusive in this regime (Fig. 1(a)).

The results of the computed values of $D$ for $\Gamma>88$ are shown in Fig.
1(b). For $\Gamma\lesssim 88$, the exponent of MSD($t$) is $\mu\neq$ 1 at all
times (Fig. 1(a)), which indicates anomalous diffusion, as a result, $D$ is
meaningless and the SE relation cannot be examined. Therefore, in the
following, we examine the SE relation for the temperature range
$88<\Gamma<\Gamma_{\rm{m}}.$

Figure 1:  (a) Mean-squared displacements
$\langle(\Delta\mathit{r}(t))^{2}\rangle$ as a function of time for different
temperatures. For short times, the motion is ballistic, MSD($t$)$\propto
t^{2}$. At later time, For $\Gamma\lesssim 88$, $\mu$ is not 1.00 at all times
indicating anomalous diffusion, but for $\Gamma>88$, the motion is diffusive,
MSD($t$)$\propto t$. (b) Temperature dependence of self-diffusion coefficient
$D$. The error bars for $D$ are much smaller than the symbol size.
$T_{\rm{onset}}$ specifies the temperature threshold in which the SE relation
is violated.

Figure 2: Autocorrelation function of the shear stress $C_{\eta}(t)$ and its
time integration for $\Gamma=128.5$. Results are shown three independet
simulation runs (a) Run 1, (b) Run2, (c) Run3 with different random number
seeds for the initial velocities. The upper limit of the Green-Kubo integral
in Eq. (4), $\omega_{pd}t_{I}$, is chosen as the time when this noisy
$C_{\eta}(t)$ first crosses zero, indicated by the circle and
arrow.$\eta/nm\omega_{pd}a^{2}$ is the normalized viscosity obtaine from the
integral. Note that for better clarity, the autocorrelation function of the
shear stress $C_{\eta}(t)$ is multiplied by a factor of 10.

###  III.2 Shear viscosity

In dusty plasma studies, the viscosity has been quantified using two methods,
with and without macroscopic velocity gradient. In dusty plasma experiments,
macroscopic shear stress is applied to generate macroscopic velocity gradient,
and viscosity is calculated using the hydrodynamic approach Gavrikov  _et al._
(2005); Vorona  _et al._ (2007); Nosenko and Goree (2004). In equilibrium
simulations, viscosity is calculated with no macroscopic velocity gradient.
Here, we calculate the viscosity using the equilibrium MD simulations and will
show that our results for viscosity values are in good agreement with those
obtained by non-equilibrium simulations. For equilibrium systems, without
gradients, the shear viscosity $\eta$ can be calculated by the Green-Kubo
relation Feng  _et al._ (2011)

$\eta=\int_{0}^{\infty}C_{\eta}(t)dt=\frac{1}{Ak_{B}T}\int_{0}^{\infty}\langle
P_{xy}(t)P_{xy}(0)\rangle dt,$ (4)

where $C_{\eta}(t)$ is the autocorrelation function of the shear stress, and
$A$ is the area of the 2D system. The off-diagonal element of the stress
tensor $P_{xy}(t)$ is given by

$P_{xy}(t)=\sum_{i=1}^{N}\left(m_{i}v_{xi}v_{yi}-\sum_{i}\sum_{j>i}\frac{x_{ij}y_{ij}}{r_{ij}}\frac{\partial\phi_{ij}}{\partial
r_{ij}}\right),$ (5)

where $r_{ij}=|\bm{\mathrm{r}}_{i}-\bm{\mathrm{r}}_{j}|$, $m_{i}$ is the mass
of the particle $i$, and $v_{xi}$ and $v_{yi}$ are the components $x$ and $y$
of the velocity of the particle $i$, respectively.

The upper limit of the time integral in Eq. (4) must be chosen cautiously. A
very long integration time due to the noise in the correlation function
$C_{\eta}(t)$ leads to unreliable values for viscosity. Therefore, in
practice, the infinite upper limit is replaced by a finite time $t_{I}$ to
avoid the statistical errors in the viscosity estimation. In agreement with
Refs. Feng  _et al._ (2011); Donkó, Goree, and Hartmann (2010); Feng, Goree,
and Liu (2012, 2013), we choose the $t_{I}$ as the time when this noisy
$C_{\eta}(t)$ first crosses zero, as shown in Figs. 2(a), (b), (c). We
calculate numerically the integral of Eq. (4) by the well-known trapezoidal
rule Press  _et al._ (1992), and the results for viscosity are shown in Fig.
3. To improve statistical accuracy, we carry out several independent
simulation runs with different random number seeds for the initial velocities,
and calculate the standard deviation from the mean (i.e., error bar) of these
simulations (Table 1).

Table 1: Results for the mean shear viscosity $\bar{\eta}$ from three independent simulation runs normalized by $nm\omega_{pd}a^{2}$ along with standard deviation from the mean $\sigma_{M}$. $\Gamma$ | $\bar{\eta}$ | $\sigma_{M}$
---|---|---
3.62 | 0.250 | 0.009
8.77 | 0.156 | 0.013
18.12 | 0.152 | 0.007
54.46 | 0.271 | 0.010
89.39 | 0.406 | 0.017
106.50 | 0.52 | 0.03
123.42 | 0.65 | 0.04
128.50 | 0.76 | 0.04
135.50 | 1.18 | 0.09

In Figure. 3, we compare the values of the viscosity coefficient obtained from
our simulation with the results obtained from the equilibrium molecular
dynamics simulation (EMD) under the same conditions in Ref. [29] and the
results obtained from the non-equilibrium molecular dynamics simulation (NEMD)
in Ref. [28]. In NEMD, an appropriate perturbation $E$ is applied, then, the
ensemble average of the resulting flux $\langle J\rangle$ is measured and the
ratio of flux and field $E$ gives the viscosity coefficient Müller-Plathe
(1999). In the method used in Ref. [28], cause and effect are reversed in a
NEMD simulation, i.e., the flux is imposed and the corresponding field is
measured. Our computed values of $\eta$ are in good agreement with the
previously simulated shear viscosity values obtained from EMD and NEMD for
2DYL Liu and Goree (2005); Donkó  _et al._ (2006).

Figure 3:  Temperature dependence of the shear viscosity $\eta$. The error
bars are calculated from three independent simulation runs. For comparison,
data are shown from the equilibrium molecular dynamics simulation (EMD) at
$\kappa=0.56$ and $N=1024$ Liu and Goree (2005), and from the non-equilibrium
molecular dynamics simulations (NEMD) at $\kappa=0.5$ and $N=1600$ Donkó  _et
al._ (2006).

Note that in addition to the equilibrium molecular dynamics method, there are
other methods for determining the viscosity coefficient, all of which are
acceptable within the framework of their assumptions. For example, Haralson
and Goree obtained the viscosity by the hydrodynamic method for low
temperatures, but as discussed in these papers Haralson and Goree (2016,
2017), they compared their results with Refs. [28,29] and found in addition to
the similar trend, their viscosity data also show similar values to
simulations. In particular, their experimental viscosity values for
$\kappa\approx 0.75$ lie between the simulation values for $\kappa=0.5$ and
$\kappa=1$ in Ref. [28].

Figure 4: Self-diffusion coefficeint $D$ as the function of shear viscosity
scaled by the temperature
$\eta\Gamma/nm\omega_{pd}a^{2}(=\eta/nm\omega_{pd}a^{2}T)$. The arrow
indicates the position of $T_{\rm{onset}}$. The SE relation is obtained as a
power-law $D\propto(\eta/T)^{-\xi}$ by the fitting of the data, and $\xi$ is
unity for $T>T_{\rm{onset}}$ corresponding to $\Gamma<\Gamma_{\rm{onset}}$
(because $\Gamma=1/T$), satisfying the SE relation. With decreasing temperatre
below $T_{\rm{onset}}$ (corresponding to increasing $\Gamma$ above
$\Gamma_{\rm{onset}}$ ), the exponent reduces to $\xi=0.7$ indicating the SE
violation.

The relationship between $D$ and $\eta/T$ is shown in Fig. 4. Data are fitted
to a power-law of the form $D\propto(\eta/T)^{-\xi}$. For the temperatures
above the onset temperature $T_{\rm{onset}}\approx 0.0087$ (corresponding to
$\Gamma_{\rm{onset}}\approx 115$), the exponent $\xi$ is unity indicating the
preservation of the SE relation. For temperatures below $T_{\rm{onset}}$ , the
exponent is $\xi=0.7$, which is different from unity indicating the violation
of the SE relation. The reason for this violation is explained as follows:
With decreasing temperature, the shear viscosity increases but the self-
diffusion coefficient decreases (Fig. 1(b)). At temperatures above the
$T_{\rm{onset}}$, $\eta$ and $D$ are coupled together, which means that rates
for increasing viscosity and decreasing self-diffusion coefficient compensate
each other, and the SE relation is preserved. However, with decreasing
temperature below the $T_{\rm{onset}}$, the rate of increase in the shear
viscosity is greater than the rate of decrease in the self-diffusion
coefficient, i.e., the increase in the shear viscosity is faster than the
decrease in the self-diffusion coefficient, consequently, this lack of
coupling between $\eta$ and $D$ leads to the violation of the SE relation in
the 2DYL near the melting point.

##  IV Time scales

To answer the question as to what is the origin of the SE violation in 2DYL,
we study the time scales that support the violation. We first focus on the
structural relaxation time of the incoherent density-density correlation
function or the self-intermediate scattering function
$F_{s}(\bm{\mathrm{k}},t)$. This function is the Fourier transform of the
distribution of the particle displacement, defined as Hansen and McDonald
(2006)

$F_{s}(\bm{\mathrm{k}},t)=\frac{1}{N}\Bigg{\langle}\sum_{\mathit{i}=1}^{N}\rm{exp}\left(-i\bm{\mathrm{k}}\cdot[\bm{\mathrm{r}}_{\mathit{i}}(t)-\bm{\mathrm{r}}_{\mathit{i}}(0)]\right)\Bigg{\rangle},$
(6)

where $\bm{\mathrm{k}}$ is the wave vector. For an isotropic system,
$F_{s}(\bm{\mathrm{k}},t)$ depends only on the magnitude
$k=|\bm{\mathrm{k}}|$, therefore, averaging over all directions yields Ghannad
(2019)

$\big{\langle}\rm{exp}\left(-i\bm{\mathrm{k}}.\bm{\mathrm{r}}\right)\big{\rangle}_{\phi}=\frac{1}{2\pi}\int_{0}^{2\pi}\rm{exp}\left(-i\mathit{kr}cos\phi\right)\mathit{d}\phi=\mathit{J}_{0}(\mathit{kr})$
(7)

where $\phi$ is the angle between the vectors $\bm{\mathrm{k}}$ and
$\bm{\mathrm{r}}$, and $\mathit{J}_{0}(\mathit{kr})$ = sin($kr$)/($kr$) is the
ordinary Bessel function of order zero. Therefore, for an isotropic system,
$F_{s}(\bm{\mathrm{k}},t)$ reduces to

$F_{s}(k,t)=\frac{1}{N}\Bigg{\langle}\sum_{\mathit{i}=1}^{N}\frac{\mathrm{sin}(k|\bm{\mathrm{r}}_{\mathit{i}}(t)-\bm{\mathrm{r}}_{\mathit{i}}(0)|)}{k|\bm{\mathrm{r}}_{\mathit{i}}(t)-\bm{\mathrm{r}}_{\mathit{i}}(0)|}\Bigg{\rangle}.$
(8)

We choose the wave number $k=|\bm{\mathrm{k}}|=4.23a^{-1}$, which is the
position of the first peak in the static structure factor $S(k)$ of the 2DYL.
For an isotropic system in two dimensions, the static structure factor is
calculated by Pathria and Beale (2011)

$S(k)=1+2\pi n\int_{0}^{\infty}r(g(r)-1)J_{0}(kr)dr$ (9)

where $g(r)$ is the radial distribution function and for a homogeneous uniform
system is defined by Haile (1992)

$g(r)=\frac{1}{n}\Big{\langle}\frac{1}{N}\sum_{i=1}^{N}\sum_{j\neq
i}^{N}\delta(r-|r_{i}-r_{j}|)\Big{\rangle},$ (10)

where $\delta$ is the Dirac delta function. The functions $g(r)$ and $S(k)$
for the 2DYL are shown in Figs. 5(a), (b).

Figure 5:  (a) Radial distribution functions $g(r)$ vs distance in units of
the Wigner-Seitz radius $a$ for the 2D Yukawa liquid in various temperatures
($\Gamma$ = 106.5, 118, 128.5, 136). (b) Static structure factor $S(k)$ vs
reduced wave number $ka$ in various temperatures. The first peak occurs at
$ka$=4.23.

In a diffusive regime with linear MSD, $F_{s}(k,t)$ is obtained theoretically
by

$F^{\mathrm{Gauss}}_{s}(k,t)=\mathrm{exp}\left(-\frac{k^{2}\langle(\Delta\mathrm{r}(t))^{2}\rangle}{4}\right)=\mathrm{exp}\left(-k^{2}Dt\right),$
(11)

which is known as the Gaussian approximation Hansen and McDonald (2006). Here,
MSD = $\langle(\Delta\it{r}(t\rm))^{2}\rangle$. Deviations of the particle
displacements from a Gaussian distribution are quantified by kurtosis. We used
the excess kurtosis or non-Gaussian parameter $\alpha_{2}(t)$, as Rahman
(1964); Charbonneau  _et al._ (2012)

$\alpha_{2}(t)=\frac{d}{d+2}\frac{\langle(\Delta\mathit{r}(t))^{4}\rangle}{\langle(\Delta\mathit{r}(t))^{2}\rangle^{2}}-1,$
(12)

where $d$ denotes spatial dimension, and For two dimensions, it reduces to

$\alpha_{2}(t)=\frac{1}{2}\frac{\langle(\Delta\mathit{r}(t))^{4}\rangle}{\langle(\Delta\mathit{r}(t))^{2}\rangle^{2}}-1.$
(13)

For a Gaussian distribution, $\alpha_{2}(t)$ approaches zero.

Figure 6:  (a) The self-intermediate scattering function $F_{s}(k,t)$ (solid
lines) as well as the non-Gaussian parameter $\alpha_{2}(t)$ (dash-dotted
lines). (b) Relationship between $\eta/T(=\eta\Gamma)$ and $\tau_{\alpha}$.
Solid line is a linear fit to the simulation results. (c) Temperature
dependence of the peak height $\alpha^{\rm{peak}}_{2}$ of the non-Gaussian
parameter $\alpha_{2}(t)$. The peak hights of $\alpha_{2}(t)$ increase with
decreasing temperature.The solid line a fit to guide the eye.

As shown in Fig. 6(a), for $\Gamma=92.5$, $\alpha_{2}(t)$ is nearly zero and
it increases with decreasing temperature ($=$ increasing $\Gamma$) indicating
the non-Gaussian behavior of the function $F_{s}(k,t)$ in the 2DYL when cools
to near the melting point.

The structure relaxation time $\tau_{\alpha}$ is defined as Becker, Poole, and
Starr (2006); Jeong  _et al._ (2010) the time where the function $F_{s}(k,t)$
decays to a value of 1/e. If $\tau_{\alpha}$ is proportional to $\eta/T$, we
have another form for SE relation as $D\propto\tau_{\alpha}^{-1}$. The
justification for replacing $\tau_{\alpha}$ with $\eta/T$ results from Eq.
(11), i.e,
$F^{\rm{Gauss}}_{s}(k,t)=\mathrm{exp}(-k^{2}Dt)\equiv\mathrm{exp}(-t/\tau)$,
consequently,
$\tau^{-1}$=$k^{2}D$=$\frac{k_{B}k^{2}}{c^{\prime}\pi}\frac{T}{\eta}$. As
shown in Fig. 6(b) , the proportional relationship
$\eta/T\propto\tau_{\alpha}$ holds in 2DYL. Therefore, the relation
$D\tau_{\alpha}=constant$, can be applied as an alternative for the SE
relation. With decreasing temperature, the SE relation is violated and the
peak hights of $\alpha_{2}(t)$ increase (Fig. 6(c)). We denote the peak time
of $\alpha_{2}(t)$ with $\tau_{\alpha_{2}}$, i.e., when the non-Gaussian
parameter reaches its maximum value. The results obtained from our simulation
data show a linear relationship between the $\tau_{\alpha_{2}}$ and
$\tau_{\alpha}$ (Fig. 7) indicating that the relation
$D\tau_{\alpha_{2}}=costant$ is another acceptable form of the SE relation.

The temperature dependence of $D\tau_{\alpha}$ and $D\tau_{\alpha_{2}}$ along
with that of $D\eta/T$ is shown in Fig. 8. At temperatures much higher than
the phase transition temperature, the SE relations, $D\eta/T,D\tau_{\alpha}$,
and $D\tau_{\alpha_{2}}$, are preserved. However, with decreasing temperature
to near the phase transition temperature, the SE relations are violated. As a
result, the SE violation is linked with deviations of the particle
displacements from Gaussian.

Figure 7: Relationship between $\tau_{\alpha}$ and $\tau_{\alpha_{2}}$. Solid
line is a linear fit to the simulation results. Figure 8: Temperature
dependence of the SE relations $D\eta/T,D\tau_{\alpha},D\tau_{\alpha_{2}},$
and $C_{0}\Gamma(1+1/\beta)/T$ normalized to their values at $T=0.0106$
$(\Gamma=94.3)$. All the SE relations indicate the SE violation with
desreasing temperature to near the phase transition temperature $T_{m}$
corresponding to increasing $\Gamma$ to near $\Gamma_{m}$.

Figure 9:  (a) Time dependence of the autocorrelation function of the shear
stress $C_{\eta}(t)$ fitted by a stretched exponential function (solid line).
(b) Relationship between $D$ and $\tau_{\eta}^{-1}$. Solid line is a linear
fit to the simulation results. (c) Temperature dependence of the stretched
exponent $\beta$. Decreasing $\beta$ with decreasing temperature implies
increasing the non-exponentiality of the shear stress autocorrelation
function. The dashed line is a guide to the eye.

Next, we look at the time dependence of the autocorrelation function of the
shear stress. It is well fitted by a stretched exponential function
$C_{\eta}(t)\simeq C_{0}\mathrm{exp}(-(t/\tau_{\eta})^{\beta})$ where $C_{0}$
is the amplitude of the stretched exponential, $\tau_{\eta}$ is the stress
decay time, and the exponent $\beta$, ranging from 0 to 1, is the degree of
non-exponentiality, which determines the degree of deviations from an
exponential (Fig. 9(a)). All three parameters are given by fitting. Then, the
shear viscosity can be approximated by substituting the stretched exponential
into Eq. (4) as

$\eta\approx\int_{0}^{\infty}C_{0}e^{-(t/\tau_{\eta})^{\beta}}dt=\frac{C_{0}\tau_{\eta}}{\beta}\Gamma(1/\beta)=C_{0}\tau_{\eta}\Gamma(1+1/\beta),$
(14)

where $\Gamma$ is the Gamma function.

Figure 10:  Relative differences between $P(r,t)$ and Gaussian distribution
$P^{\rm{Gauss}}(r,t)$ vs normalized displacement $r/a$ at time
$\omega_{pd}\tau_{\alpha_{2}}$. For $r<r_{m}$, the relative differences are
nearly zero.

The relationship between $D$ and inverse $\tau_{\eta}$ is shown in Fig. 9(b).
Our simulation data show a proportional relationship between $D$ and
$\tau_{\eta}^{-1}$ for 2DYL. As a result, by using $D\propto\tau_{\eta}^{-1}$
and Eq. (4) for $\eta$, we obtain another acceptable form for the SE relation
as $D\eta/T\propto C_{0}\Gamma(1+1/\beta)/T$ shown in Fig. 8, that confirms
the SE violation with decreasing temperature. The temperature dependence of
$\beta$ is shown in Fig. 9(c). With decreasing temperature, $\beta$ decreases,
i.e., deviations from exponentiality of the autocorrelation function of the
shear stress increases. Therefore, the SE violation is linked with deviations
of the shear stress autocorrelation function $C_{\eta}(t)$ from exponential
shape.

Overall, we find that the physical origins of the SE violation in 2DYL are
attributed to the non-exponentiality of the stress correlation function and
the non-Gaussianity of the distribution of particle displacements. The smaller
exponent $\beta$, the greater the deviation from the exponential, which
results in stronger heterogeneous dynamics. The relationship between the non-
Gaussian behavior of the distribution of particle displacements and dynamic
heterogeneity can also be explained by calculating the relative difference
between the distribution of particle displacements $P(r,t)$ and that obtained
from the Gaussian approximation $P^{\rm{Gauss}}(r,t)$ given by
$P^{\rm{Gauss}}\rm{(}\mathit{r},\it{t}\rm{)}=(1/\pi\langle(\Delta\mathit{r}(\mathit{t}))^{2}\rangle)\rm{exp}(-\it{r}\rm^{2}/\langle(\Delta\mathit{r}(\it{t}\rm))^{2}\rangle)$,
where $\langle(\Delta\mathit{r}(t))^{2}\rangle$ is obtained from the
simulation. The distribution of particle displacements $P(\bm{\mathrm{r}},t)$
is given by Haile (1992)

$P(\bm{\mathrm{r}},t)=\frac{1}{N}\Bigg{\langle}\sum_{\mathit{i}=1}^{N}\delta(\bm{\mathrm{r}}-\bm{\mathrm{r}}_{i}(t)+\bm{\mathrm{r}}_{i}(0))\Bigg{\rangle},$
(15)

For isotropic liquids, $P(\bm{\mathrm{r}},t)$ depends only on the scalar
distance, $r=|\bm{\mathrm{r}}|$. Thus, Eq. (15) reduces to

$P(r,t)=\frac{1}{N}\Bigg{\langle}\sum_{\mathit{i}=1}^{N}\delta(r-|r_{i}(t)-r_{i}(0)|)\Bigg{\rangle}.$
(16)

Physically, 2$\pi rP(r,t)dr$ measures the probability of finding a dust
particle at distance $r$ from an origin at time $t$ given that the same dust
particle was at the origin at the initial time $t=0$ Haile (1992). Thus,
$P(r,t)$ is normalized by

$\int P(r,t)d\bm{\mathrm{r}}=1.$ (17)

The relative differencs between $P\rm{(}\it{r},t\rm{)}$ and
$P^{\rm{Gauss}}\rm{(}\it{r},t\rm{)}$ are shown in Fig. 10, for
$t=\tau_{\alpha_{2}}$. We first note that with decreasing temperature, the
relative difference increases, which supports the violation of the SE with
decreasing temperature corresponding to increasing $\Gamma$. For displacements
smaller than $r_{m}\simeq 1.7a$, the relative difference between $P(r,t)$ and
$P^{\rm{Gauss}}(r,t)$ is nearly zero, i.e.,
$P(r_{m},\tau_{\alpha_{2}})=P^{\rm{Gauss}}(r_{m},\tau_{\alpha_{2}})$. However,
for $r>r_{m}$, the relative difference increases and becomes as large as
$10^{3}$, which belongs to the lowest temperature, i.e., $T=0.00735$ or
$\Gamma=136$. It means that a significant number of particles travel farther
than $r_{m}$ at time $\omega_{pd}\tau_{\alpha_{2}}$. It reflects the existence
of dynamic heterogeneity with decreasing temperature, i.e., the regions are
formed in 2DYL that dust particles are more mobile than expected from a
Gaussian approximation.

##  V CONCLUSIONS

We have reported extensive molecular dynamics results to recognize the origin
of the Stokes-Einstein violation in 2D Yukawa liquids. First, we have computed
the mean-squared displacement MSD($t$) at different temperature ranges in the
Yukawa liquid. Then, for the temperature ranges that MSD($t$) is linear and
therefore the self-diffusion coefficient $D$ is meaningful, we have computed
$D$ from the linear fit of the MSD($t$) versus time. Then, we have calculated
the shear viscosity using the Green-Kubo integral with the assumption that
there is no macroscopic velocity gradient in the equilibrium. We have found
that there is a temperature threshold $T_{\rm{onset}}$ in which the Stokes-
Einstein relation, $D\propto(\eta/T)^{-1}$, is preserved for
$T>T_{\rm{onset}}$ but is violated for $T<T_{\rm{onset}}$.

To identify the time scales supporting this violation, we have calculated the
structural relaxation time $\tau_{\alpha}$ from the incoherent density-density
correlation function or the self-intermediate scattering function $F_{s}(k,t)$
when decays to a value of $\rm{e^{-1}}$. We have shown that $\tau_{\alpha}$ is
proportional to $\eta/T$. As a result, We have defined that an alternative
relationship to the Stokes-Einstein relation as $D\propto\tau_{\alpha}^{-1}$.
Then, by calculating the peak time of deviations of the particle displacements
from a Gaussian distribution, $\tau_{\alpha_{2}}$, we have shown that there is
a linear relation between $\tau_{\alpha_{2}}$ and $\tau_{\alpha}$ indicating
that the relation $D\propto\tau_{\alpha_{2}}^{-1}$ is another acceptable form
of the SE relation. Then, we have shown that with decreasing temperature to
below the phase transition temperature $\Gamma_{m}$, $D\tau_{\alpha}$ and
$D\tau_{\alpha_{2}}$ are violated which implies that the SE violation is
linked with the non-Gaussian behavior of the particle displacements
distribution. With calculation $\eta$ from the autocorrelation function of the
shear stress approximated by the stretched exponential,
$C_{\eta}\propto\mathrm{exp}(-(t/\tau_{\eta})^{\beta})$, we have found another
SE relation, and have shown that with decreasing temperature to below
$\Gamma_{m}$, $\beta$ decreases, i.e., deviation from exponentiality of
$C_{\eta}(t)$ increases, indicating a link between deviations of the
autocorrelation function of the shear stress from the exponential shape and
the SE relation.

Generally, our results provide a deep insight into understanding the breakdown
of the SE relation, not only in 2D Yukawa liquids but also in other strongly
coupled systems such as one-component Coulomb liquids and 3D Yukawa liquids.
For these 3D systems, failure of the SE relation has been reported at high
temperatures, in contrast to 2D Yukawa liquids Daligault (2006); Donkó and
Hartmann (2008). Therefore, it will be interesting to examine the effect of
the spatial dimension on the origin of the SE violation.

## DATA AVAILABILITY

The data that support the findings of this study are available from the author
upon reasonable request.

## REFERENCES

## References

  * Einstein (1956) A. Einstein,  _Investigations on the Theory of the Brownian Movement_ (Dover, New York, 1956).
  * Hansen and McDonald (2006) J. Hansen and I. R. McDonald,  _Theory of Simple Liquids_ (Academic, London, 2006).
  * Liu, Goree, and Vaulina (2006) B. Liu, J. Goree, and O. S. Vaulina, Phys. Rev. Lett. 96, 015005 (2006).
  * Pan, Garrahan, and Chandler (2005) A. C. Pan, J. P. Garrahan, and D. Chandler, Phys. Rev. E. 72, 041106 (2005).
  * Dubey  _et al._ (2019) V. Dubey, S. Erimban, S. Indra, and S. Daschakraborty, J. Phys. Chem. B 123, 10089 (2019).
  * Kawasaki and Kim (2017) T. Kawasaki and K. Kim, Sci. Adv. 3, e1700399 (2017).
  * Kawasaki and Kim (2019a) T. Kawasaki and K. Kim, Sci. Rep. 9, 8118 (2019a).
  * Kawasaki and Kim (2019b) T. Kawasaki and K. Kim, J. Stat. Mech. 2019, 084004 (2019b).
  * Becker, Poole, and Starr (2006) S. R. Becker, P. H. Poole, and F. W. Starr, Phys. Rev. Lett. 97, 055901 (2006).
  * Jeong  _et al._ (2010) D. Jeong, M. Y. Choi, H. J. Kimwcd, and Y. J. Jung, Phys. Chem. Chem. Phys. 12, 2001 (2010).
  * Köddermann, Ludwig, and Paschek (2008) T. Köddermann, R. Ludwig, and D. Paschek, ChemPhysChem 9, 1851 (2008).
  * Tsimpanogiannis  _et al._ (2019) I. N. Tsimpanogiannis, S. H. Jamali, I. G. Economou, T. J. H. Vlugt, and O. A. Moultos, Mol. Phys. 118, 1 (2019).
  * Fernandez-Alonso  _et al._ (2007) F. Fernandez-Alonso, F. J. Bermejo, S. E. McLain, J. F. C. Turner, J. J. Molaison, and K. W. Herwig, Phys. Rev. Lett. 98, 077801 (2007).
  * Mallamace  _et al._ (2010) F. Mallamace, C. Branca, C. Corsaro, N. Leone, J. Spooren, H. E. Stanley, and S.-H. Chen, ChemPhysChem 114, 1870 (2010).
  * Kob  _et al._ (1997) W. Kob, C. Donati, S. J. Plimpton, P. H. Poole, and S. C. Glotzer, Phys. Rev. Lett. 79, 2827 (1997).
  * Yukawa (1935) H. Yukawa, Proc. Phys. Math. Soc. Jpn. 17, 48 (1935).
  * Konopka, Morfill, and Ratke (2000) U. Konopka, G. E. Morfill, and L. Ratke, Phys. Rev. Lett. 84, 891 (2000).
  * Ghannad (2019) Z. Ghannad, Phys. Rev. E 100, 033211 (2019).
  * Feng  _et al._ (2016) Y. Feng, J. Goree, B. Liu, L. Wang, and W. Tian, J. Phys. D: Appl. Phys. 49, 235203 (2016).
  * Wang, Huang, and Feng (2018) K. Wang, D. Huang, and Y. Feng, J. Phys. D: Appl. Phys. 51, 245201 (2018).
  * Feng, lin, and Murillo (2017) Y. Feng, W. lin, and M. S. Murillo, Phys. Rev. E 96, 053208 (2017).
  * Tarjus and Kivelson (1995) G. Tarjus and D. Kivelson, J. chem. Phys. 103, 3071 (1995).
  * Pan  _et al._ (2017) S. Pan, Z. W. Wu1, W. H. Wang, M. Z. Li, and L. Xu, Sci. Rep. 7, 39938 (2017).
  * Sengupta  _et al._ (2013) S. Sengupta, S. Karmakar, C. Dasgupta, and S. Sastry, J. chem. Phys. 138, 12A548 (2013).
  * Frenkel and Smit (2002) D. Frenkel and B. Smit,  _Understanding Molecular Dynamics Simulation_ (Academic, San Diego, 2002).
  * Ohta and Hamaguchi (2000) H. Ohta and S. Hamaguchi, Phys. Plasmas 7, 4506 (2000).
  * Hartmann  _et al._ (2019) P. Hartmann, J. C. Reyes, E. G. Kostadinova, L. S. Matthews, T. W. Hyde, R. U. Masheyeva, K. N. Dzhumagulova, T. S. Ramazanov, T. Ott, H. Kählert, M. Bonitz, I. Korolov, and Z. Donkó, Phys. Rev. E. 99, 013203 (2019).
  * Donkó  _et al._ (2006) Z. Donkó, J. Goree, P. Hartmann, and K. Kutasi, Phys. Rev. Lett. 96, 145003 (2006).
  * Liu and Goree (2005) B. Liu and J. Goree, Phys. Rev. Lett. 94, 185002 (2005).
  * Hartmann  _et al._ (2005) P. Hartmann, G. J. Kalman, Z. Donkó, and K. Kutasi, Phys. Rev. E 72, 026409 (2005).
  * Kalman  _et al._ (2004) G. J. Kalman, P. Hartmann, Z. Donkó, and M. Rosenberg, Phys. Rev. Lett 92, 065001 (2004).
  * Swope  _et al._ (1982) W. C. Swope, H. C. Andersen, P. H. Berens, and K. R. Wilson, J. Chem. Phys. 76, 637 (1982).
  * Nosé (1984) S. Nosé, J. Chem. Phys.. 81, 511 (1984).
  * Hoover (1985) W. G. Hoover, Phys. Rev. A 31, 1695 (1985).
  * Holian, Voter, and Ravelo (1995) B. L. Holian, A. F. Voter, and R. Ravelo, Phys. Rev. E 52, 2338 (1995).
  * Maginn  _et al._ (2019) E. J. Maginn, R. A. Messerly, D. J. Carlson, D. R. Roe, and J. R. Elliott, Living J. Comput. Mol. Sci. 1, 6324 (2019).
  * Gavrikov  _et al._ (2005) A. Gavrikov, I. Shakhova, A. Ivanov, O. Petrov, N. Vorona, and V. Fortov, Phys. Lett. A 336, 378 (2005).
  * Vorona  _et al._ (2007) N. A. Vorona, A. V. Gavrikov, A. S. Ivanov, O. F. Petrov, V. E. Fortov, and I. A. Shakhova, J. Exp. Theor. Phys. 105, 824 (2007).
  * Nosenko and Goree (2004) V. Nosenko and J. Goree, Phys. Rev. Lett 93, 155004 (2004).
  * Feng  _et al._ (2011) Y. Feng, J. Goree, B. Liu, and E. G. D. Cohen, Phys. Rev. E 84, 046412 (2011).
  * Donkó, Goree, and Hartmann (2010) Z. Donkó, J. Goree, and P. Hartmann, Phys. Rev. E. 81, 056404 (2010).
  * Feng, Goree, and Liu (2012) Y. Feng, J. Goree, and B. Liu, Phys. Rev. E 85, 066402 (2012).
  * Feng, Goree, and Liu (2013) Y. Feng, J. Goree, and B. Liu, Phys. Rev. E. 87, 013106 (2013).
  * Press  _et al._ (1992) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,  _Numerical Recipes in C: The Art of Scientific Computing_ (Cambridge University Press, Cambridge, 1992).
  * Müller-Plathe (1999) Müller-Plathe, Phys. Rev. E. 59, 4894 (1999).
  * Haralson and Goree (2016) Z. Haralson and J. Goree, Phys. Plasmas 23, 093703 (2016).
  * Haralson and Goree (2017) Z. Haralson and J. Goree, Phys. Rev. Lett 118, 195001 (2017).
  * Pathria and Beale (2011) R. K. Pathria and P. D. Beale,  _Statistical Mechanics_ (Elsevier, Amesterdam, 2011).
  * Haile (1992) J. M. Haile,  _Molecular Dynamics Simulation: Elementary Methods_ (John Wiley and Sons, New York, 1992).
  * Rahman (1964) A. Rahman, Phys. Rev. 136, A405 (1964).
  * Charbonneau  _et al._ (2012) P. Charbonneau, A. Ikedac, G. Parisid, and F. Zamponi, Proc. Natl.Acad. Sci. USA 109, 13939 (2012).
  * Daligault (2006) J. Daligault, Phys. Rev. Lett. 96, 065003 (2006).
  * Donkó and Hartmann (2008) Z. Donkó and P. Hartmann, Phys. Rev. E. 78, 026408 (2008).