The following article has been submitted to The Journal of Chemical Physics.

©2021 Hiroki Kusudo. This article is distributed under a Creative Commons
Attribution (CC BY) License.

#  Local stress tensor calculation by the Method-of-Plane in microscopic
systems with macroscopic flow: a formulation based on the velocity
distribution function

Hiroki Kusudo<EMAIL_ADDRESS>Department of Mechanical
Engineering, Osaka University, 2-1 Yamadaoka, Suita 565-0871, Japan  Takeshi
Omori<EMAIL_ADDRESS>Department of Mechanical and Physical Engineering,
Osaka City University, 3-3-138 Sugimoto, Sumiyoshi-ku, Osaka 558-8585, Japan
Yasutaka Yamaguchi<EMAIL_ADDRESS>Department of Mechanical
Engineering, Osaka University, 2-1 Yamadaoka, Suita 565-0871, Japan  Water
Frontier Research Center (WaTUS), Research Institute for Science & Technology,
Tokyo University of Science, 1-3 Kagurazaka, Shinjuku-ku, Tokyo, 162-8601,
Japan

###### Abstract

In this work, we showed a calculation method of local stress tensor applicable
to non-equilibrium MD systems based on the Method of Plane (MoP). From the
relation between the macroscopic velocity distribution function and the
microscopic molecular passage across a fixed control plane, we derived a
method to calculate the basic properties of the macroscopic momentum
conservation law including the density, the velocity, the momentum flux, the
interaction and kinetic terms of the stress tensor defined on a surface with a
finite area. Any component of the streaming velocity can be obtained on a
control surface, which enables the separation of the kinetic momentum flux
into the advection and stress terms in the framework of MoP. We verified the
present method through the extraction of the density, velocity and stress
distributions in a quasi-1D steady-state Couette flow system and in a quasi-2D
steady-state system with a moving contact line. In our method, as opposed to
volume average method, the density, mass and momentum fluxes are defined on a
surface, which is essential to be consistent with the mass and momentum
conservation laws in dynamic systems.

##  I Introduction

With the increasing interest in microfluidic devices and nanotechnologies,
molecular dynamics (MD) simulations have become a powerful computational tool
to examine the fluid behavior for small scale systems. For the understanding
in terms of flow fields, the microscopic motion of individual molecules must
be averaged, and the stress tensor plays a key role in such macroscopic flow.
Within the framework of fluid mechanics, the stress tensor is determined from
the velocity fields through the constitutive equation typically including the
viscosity, and the local acceleration of the fluid is given by the gradient of
the stress tensor as well as the external field to satisfy the momentum
conservation. On the other hand, the molecular motion is governed by the
intermolecular interaction, and the stress tensor should be defined through
the average of the molecular motion and interaction.

More concretely, the macroscopic equation of continuity is given with the
density $\rho(\mbox{\boldmath$x$},t)$ and velocity vector
$\mbox{\boldmath$u$}(\mbox{\boldmath$x$},t)$ both as functions of position $x$
and time $t$ by

$\displaystyle\frac{\partial\rho}{\partial t}+\frac{\partial\rho
u_{k}}{\partial x_{k}}=0,$ (1)

where $x_{k}$ and $u_{k}$ are the $k$-direction components of the position and
the velocity vector, respectively. The Einstein notation is used with a dummy
index $k$ for the vectors. The second term of the LHS is called the advection
or streaming term. Equation (1) comes from the mass conservation

$\displaystyle\int\\!\\!\\!\int\\!\\!\\!\int_{V}\mathrm{d}V\frac{\partial\rho}{\partial
t}=-\int\\!\\!\\!\int_{S}\mathrm{d}S\rho u_{k}n_{k}$ (2)

satisfied for an arbitrary volume $V$ in an enclosing surface $S$, where
$n_{k}$ is the $k$-direction component of the outward unit normal vector $n$
with respect to the infinitesimal surface element $\mathrm{d}S$. This integral
form in Eq. (2) means that the mass in $V$ can be changed only by the mass
flux passing the surface $S$. By applying Gauss’ divergence theorem for the
RHS as

$\displaystyle\int\\!\\!\\!\int_{S}\mathrm{d}S\rho
u_{k}n_{k}=\int\\!\\!\\!\int\\!\\!\\!\int_{V}\mathrm{d}V\frac{\partial\rho
u_{k}}{\partial x_{k}},$ (3)

it follows for Eq. (2) that

$\displaystyle\int\\!\\!\\!\int\\!\\!\\!\int_{V}\mathrm{d}V\left(\frac{\partial\rho}{\partial
t}+\frac{\partial\rho u_{k}}{\partial x_{k}}\right)=0,$ (4)

for an arbitrary $V$, which corresponds to Eq. (1).

Similarly, the macroscopic momentum equation, the Navier-Stokes equation, is
given by

$\displaystyle\frac{\partial\rho u_{l}}{\partial t}+\frac{\partial\rho
u_{l}u_{k}}{\partial x_{k}}=\frac{\partial\tau_{kl}}{\partial x_{k}}+\rho
F_{l},$ (5)

where the fluid stress tensor component $\tau_{kl}$ expresses the stress in
the $l$-direction exerted on a surface element with an outward normal in the
$k$-direction, and $F_{l}$ denotes the external force per mass. This is also
derived from the momentum conservation for an arbitrary volume $V$ enclosed by
$S$:

$\displaystyle\int\\!\\!\\!\int\\!\\!\\!\int_{V}\mathrm{d}V\frac{\partial\rho
u_{l}}{\partial t}=-\int\\!\\!\\!\int_{S}\mathrm{d}S\rho
u_{l}u_{k}n_{k}+\int\\!\\!\\!\int_{S}\mathrm{d}S\tau_{kl}n_{k}+\int\\!\\!\\!\int\\!\\!\\!\int_{V}\mathrm{d}V\rho
F_{l},$ (6)

meaning that the total momentum in $V$ can be changed by the momentum flux
passing the surface $S$ as well as the impulse due to the stress exerted on
the surface $S$ and external force exerted on the volume $V$. Specifically
note that the advection $\rho u_{l}u_{k}$ and stress $\tau_{kl}$ in the 1st
and 2nd terms of the RHS, respectively, are separated.

In contrast to the above-mentioned macroscopic feature, microscopic molecules
have their own velocities, and their local average and variance correspond to
the macroscopic velocity $u$ and temperature under the assumption of local
equilibrium in space and time, respectively. Hence, an instantaneous
microscopic mass or momentum transfer across an arbitrary surface $S$ exists
due to the passage of the constituent molecules even in macroscopically static
system without mean flow with $\mbox{\boldmath$u$}=\mbox{\boldmath$0$}$. In
the kinetic theory of gases based on the Boltzmann equation with respect to
the velocity distribution function (VDF) of the constituent molecules, the
macroscopic velocity is subtracted from the microscopic molecular velocity
upon the definition of the stress. Before the establishment of the MD method,
Irving and Kirkwood [1] (IK) put forward the connection between the
macroscopic conservation laws and microscopic molecular motion governed by the
inter-molecular interaction through the statistical mechanical theory based on
the distribution function in the phase space, and derived an expression of the
local pointwise stress comprised of kinetic and interaction parts including
Taylor series expansion of differences in delta functions to express the
microscopic particle feature.

After the introduction of numerical MD simulations, [2, 3] the calculation of
average stress in homogeneous bulk systems, equivalent to the average bulk
pressure with its sign inverted, was enabled based on the virial theorem,
which indeed corresponds to the average of the IK form integrated in space and
time. From this bulk stress, the viscosity as a transport coefficient can also
be obtained based on the Green-Kubo relation from the time-fluctuation of the
off-diagonal stress component averaged in space. [4, 5] Regarding the local
stress implemented for MD simulations with a discrete time-step, Tsai [6]
proposed a pragmatic scheme to calculate the averaged stress defined on a flat
plane in quasi-1D planer systems, and Thompson  _et al._ [7] extended the
approach toward a spherical curved surface for the analysis of the surface
tension of a spherical droplet. In both cases, all the momentum flux and
intermolecular force across the plane or the sphere, which divide the
computational domain, were summed up during the time integration in
macroscopically static systems. This type of stress definition is usually
called the Method of Plane (MoP), [8] or Hardy [9] stress, where the momentum
conservation law with the above-mentioned MoP for quasi-1D systems was proved
to be an exact consequences of Newton’s laws, [9] i.e., the MoP meets Gauss’
divergence theorem for any control volume (CV) surrounded by enclosing
surface(s) irrespective of whether the local system in the control volume is
homogeneous or not. In addition, this momentum conservation is not restricted
to quasi-1D systems but also applicable to quasi-2D systems, e.g., the present
authors adopted the MoP for a CV with an rectangular enclosing surface set
around the contact line of an equilibrium droplet to examine the nanoscale
wetting behavior through the mechanical balance exerted on the fluid in the
CV. [10] On the other hand, the volume average (VA) is proposed as another
approach that gives the local mean stress in space, where weighted average of
the pair interaction in local CVs is included in the formulation. [11, 12, 13,
14, 15, 16, 17, 18, 19] This VA should in principle be applied only for
homogeneous CVs, and is advantageous especially to explore the link with the
local fluctuation and thermodynamic limit, [16, 17] while stress integral can
also be obtained by the VA, and this enables the calculation of the surface
tension based on Bakker’s equation using this stress integral, [20] which is
known as a mechanical route to evaluate the surface tension. [21, 12, 22, 23,
5, 19] The momentum conservation is satisfied for the whole system if the VA
is properly summed up; however, special care is needed to consider the
momentum conservation for local CVs because the VA originally was not designed
to examine local momentum conservation to be satisfied through the link to
Gauss’ divergence theorem. [9, 14, 18] This feature is similar to the atomic
stress, [24] for instance provided as stress/atom command in LAMMPS package,
[25] often used to simply visualize the stress field.

Going back to the momentum conservation in Eq. (6), for systems with a non-
zero local flow, i.e., macroscopically dynamic systems with
$\mbox{\boldmath$u$}\neq\mbox{\boldmath$0$}$, the local macroscopic velocity
$u$ must be defined on a surface $S$ enclosing a control volume $V$ in MD
systems so that the stress calculated in MD simulations may be consistent with
the macroscopic momentum Eq. (5). In other words, the momentum transfer due to
the microscopic molecular passage across a control surface should be separated
into advection and stress contributions to examine the local flow from a
macroscopic point of view.

In this paper, we show a calculation method of the MoP-based local stress
tensor applicable to non-equilibrium molecular dynamics (NEMD) systems. We
provide the formulation for systems consisting of single-component mono-atomic
fluid molecules for simplicity while the present framework is also applicable
to systems of multi-component or poly-atomic fluid molecules. For the
derivation, we introduced the VDF to clearly give the average of physical
properties defined on a fixed control plane as we shall see later. With this
procedure, we provide not only the expression of the kinetic term of stress
tensor but also those of the components in the advection term, i.e., the
density and macroscopic velocity as an extension of the MoP. To check its
validity, we performed test calculation in two systems: 1) a quasi-1D Couette
flow system and 2) a quasi-2D system with liquid-solid-vapor contact lines,
both consisted of a Lennard-Jones fluid between parallel solid walls moving in
the opposite directions tangential to the walls. In the first system, we
compared the density and velocity distributions obtained by the present method
and the VA, and we calculated the distributions of the stress components and
advection term. Furthermore, we showed that the same velocity distribution was
obtained on bin faces with different normal directions, which is essential to
determine the advection term. In the second system, the density, velocity and
stress distributions are calculated in the complex flow with liquid-vapor
interfaces and contact lines.

##  II Theory

We show the derivation of the stress averaged on a finite bin face in a
Cartesian coordinate system for single-component mono-atomic fluid in the
following for simplicity. Note that the Einstein notation with dummy indices
used in Sec. I is not applied hereafter. The fluid stress tensor component
$\tau_{kl}$ which expresses the stress in the $l$-direction exerted on a
surface element with an outward normal in the $k$-direction, is given by the
kinetic term $\tau_{kl}^{\mathrm{kin}}$ and the inter-molecular interaction
term $\tau_{kl}^{\mathrm{int}}$ as

$\displaystyle\tau_{kl}=\tau_{kl}^{\mathrm{kin}}+\tau_{kl}^{\mathrm{int}}.$
(7)

In the standard MoP for equilibrium MD systems without mean flow consisting of
single-component mono-atomic fluid molecules, [8, 10] the kinetic term
$\tau_{kl}^{\mathrm{kin}}$ in Eq. (7) on a bin face of area $S_{k}$ with its
normal vector pointing to the $k$-th Cartesian direction is calculated by

$\tau_{kl}^{\mathrm{kin}}\equiv-\frac{1}{S_{k}\delta
t}\left<\sum_{i\in\mathrm{fluid},\delta t}^{\text{across\
}S_{k}}m^{i}v_{l}^{i}\frac{v_{k}^{i}}{\left|v_{k}^{i}\right|}\right>,$ (8)

where $m^{i}$ and $v_{l}^{i}$ denote the mass and $l$-component of the
velocity vector $\mbox{\boldmath$v$}^{i}$ of fluid particle $i$, respectively.
We also denote the bin face by $S_{k}$ hereafter. The angular brackets denote
the ensemble average, and the summation $\sum_{i\in\mathrm{fluid},\delta
t}^{\text{across\ }S_{k}}$ is taken for every fluid particle $i$ passing
through $S_{k}$ within a time interval of $\delta t$, which is equal to the
time increment for the numerical integration. Considering that we deal with a
single-component fluid molecules of an identical mass $m$, we substitute
$m^{i}$ with $m$ hereafter. A sign function $\frac{v_{k}^{i}}{|v_{k}^{i}|}$
equal to $\pm 1$ is multiplied to the momentum transfer $mv_{l}^{i}$ across
$S_{k}$ to evaluate the kinetic effect on the stress depending on the passing
direction. Note that in static equilibrium systems, i.e., systems without
macroscopic local mean flow, the advection term is zero in the whole system.

On the other hand, the intermolecular interaction term
$\tau_{kl}^{\mathrm{int}}$ in Eq. (7) in the case of simple two body potential
is calculated by

$\displaystyle\tau_{kl}^{\mathrm{int}}=-\frac{1}{S_{k}}\left<\sum_{(i,j)\in\mathrm{fluid}}^{\text{across\
}S_{k}}F^{ij}_{l}\frac{r^{ij}_{k}}{|r^{ij}_{k}|}\right>,$ (9)

where $r^{ij}_{k}$ and $F^{ij}_{l}$ denote the $k$-component of the relative
position vector
$\mbox{\boldmath$r$}^{ij}\equiv\mbox{\boldmath$x$}^{j}-\mbox{\boldmath$x$}^{i}$
and the $l$-component of the force vector $\mbox{\boldmath$F$}^{ij}$ on
particle $j$ at position $\mbox{\boldmath$x$}^{j}$ from particle $i$ at
position $\mbox{\boldmath$x$}^{i}$, respectively. The summation
$\sum_{(i,j)\in\mathrm{fluid}}^{\text{across\ }S_{k}}$ is taken for all line
segments of the inter-particle interaction between $\mbox{\boldmath$x$}^{i}$
and $\mbox{\boldmath$x$}^{j}$ which cross $S_{k}$. A sign function
$\frac{r^{ij}_{k}}{|r^{ij}_{k}|}$ is multiplied for this interaction term to
evaluate the force effect depending on the force direction. Note that
technically the fluid-solid interaction can also be included as $i$-$j$ pair
in the summation $\sum_{(i,j)}^{\mathrm{across\ }S_{k}}$ in Eq. (9), but only
the fluid-fluid interaction was taken into account as the fluid stress, and
fluid-solid contribution was considered as an external force field. [10, 26,
27] Also note that for multi-component systems or systems with poly-atomic
molecules, difficulties mainly arise to treat the interaction force between
different kind of molecules or the constraint force [28] of the polyatomic
molecules, where the interaction forces should be properly implemented into
the stress calculation to satisfy the conservation laws. [29]

To extend the standard MoP to steady-state NEMD systems with a non-zero
macroscopic mean local flow, the mean velocity should be properly subtracted
from the kinetic term $\tau_{kl}^{\mathrm{kin}}$ in Eq. (8) so that the
macroscopic momentum flux as the advection term due to the mean velocity $u$
may be included not in the stress term but in the advection term within the
macroscopic description of the momentum conservation, i.e., in the Navier-
Stokes equation (5). In the following, we provide a general framework to
connect a microscopic variable $\xi^{i}$ of particles and a macroscopic field
value $\xi(\mbox{\boldmath$x$},t)$ averaged on $S_{k}$ under non-zero mean
velocity based on the local VDF in the Cartesian $xyz$-coordinate system.

At first, we define the VDF $f(\mbox{\boldmath$x$},\mbox{\boldmath$v$},t)$ for
the mass with a velocity $\mbox{\boldmath$v$}=(v_{x},v_{y},v_{z})$ at position
$\mbox{\boldmath$x$}=(x,y,z)$ at time $t$, which gives the local density
$\rho(\mbox{\boldmath$x$},t)$ by

$\displaystyle\rho(\mbox{\boldmath$x$},t)=\int_{-\infty}^{\infty}\mathrm{d}v_{x}\int_{-\infty}^{\infty}\mathrm{d}v_{y}\int_{-\infty}^{\infty}\mathrm{d}v_{z}f(\mbox{\boldmath$x$},\mbox{\boldmath$v$},t)\equiv\int\\!\\!\\!\int\\!\\!\\!\int_{-\infty}^{\infty}\mathrm{d}\mbox{\boldmath$v$}f(\mbox{\boldmath$x$},\mbox{\boldmath$v$},t),$
(10)

where we rewrite
$\int_{-\infty}^{\infty}\mathrm{d}v_{x}\int_{-\infty}^{\infty}\mathrm{d}v_{y}\int_{-\infty}^{\infty}\mathrm{d}v_{z}$
by
$\int\\!\\!\\!\int\\!\\!\\!\int_{-\infty}^{\infty}\mathrm{d}\mbox{\boldmath$v$}$.
Then a microscopic variable $\xi^{i}$ per mass of particle $i$ can be related
to a corresponding macroscopic field variable $\xi(\mbox{\boldmath$x$},t)$ as

$\lim_{\delta t\to 0}\left<\sum_{i\in\mathrm{fluid},\delta t}^{\text{crossing
}S_{k}}m\xi^{i}\right>\equiv\lim_{\delta t\to
0}\int\\!\\!\\!\int\\!\\!\\!\int_{-\infty}^{\infty}\mathrm{d}\mbox{\boldmath$v$}\int_{0}^{\left|v_{k}\right|\delta
t}\mathrm{d}x_{k}S_{k}f(\mbox{\boldmath$x$},\mbox{\boldmath$v$},t)\xi(\mbox{\boldmath$x$},t).$
(11)

The RHS denotes the integral weighted with VDF considering an oblique pillar
of a base area $S_{k}$ and a height $|v_{k}|\delta t$ with its central axis
parallel to $v$ , which is typically assumed upon the derivation of the
equilibrium pressure in the kinetic theory of gases.

With the limit $\delta t\to 0$, and by rewriting the average of
$f(\mbox{\boldmath$x$},\mbox{\boldmath$v$},t)$ and
$\xi(\mbox{\boldmath$x$},t)$ in the oblique pillar by
$f(S_{k},\mbox{\boldmath$v$},t)$ and $\xi(S_{k},t)$, respectively, the
integral with respect to $x_{k}$ in the RHS of Eq. (11) writes

$\lim_{\delta t\to 0}\int_{0}^{\left|v_{k}\right|\delta
t}\mathrm{d}x_{k}S_{k}f(\mbox{\boldmath$x$},\mbox{\boldmath$v$},t)\xi(\mbox{\boldmath$x$},t)=\lim_{\delta
t\to 0}S_{k}f(S_{k},\mbox{\boldmath$v$})\xi(S_{k})\left|v_{k}\right|\delta t,$
(12)

and it follows for Eq. (11) that

$\lim_{\delta t\to 0}\left<\sum_{i\in\mathrm{fluid},\delta t}^{\text{crossing
}S_{k}}m\xi^{i}\right>=\lim_{\delta t\to
0}S_{k}\int\\!\\!\\!\int\\!\\!\\!\int_{-\infty}^{\infty}\mathrm{d}\mbox{\boldmath$v$}f(S_{k},\mbox{\boldmath$v$},t)\xi(S_{k},t)\left|v_{k}\right|\delta
t.$ (13)

Hence, by dividing both sides by $S_{k}\delta t$,

$\int\\!\\!\\!\int\\!\\!\\!\int_{-\infty}^{\infty}\mathrm{d}\mbox{\boldmath$v$}f(S_{k},\mbox{\boldmath$v$},t)\xi(S_{k},t)\left|v_{k}\right|=\lim_{\delta
t\to 0}\frac{1}{S_{k}\delta t}\left<\sum_{i\in\mathrm{fluid},\delta
t}^{\text{crossing }S_{k}}m\xi^{i}\right>$ (14)

is derived as a basic equation for the connection between the macroscopic
field variable $\xi(S_{k},t)$ and microscopic variable $\xi^{i}$ which belongs
to the constituent particle $i$ upon crossing $S_{k}$.

Now, we proceed to the expressions of the macroscopic field variables averaged
on $S_{k}$. By substituting $\xi(S_{k},t)$ and $\xi^{i}$ in Eq. (14) with
$\frac{1}{|v_{k}|}$ and $\frac{1}{|v^{i}_{k}|}$, respectively, and using Eq.
(10), it follows

$\rho(S_{k},t)=\lim_{\delta t\to 0}\frac{1}{S_{k}\delta
t}\left<\sum_{i\in\mathrm{fluid},\delta t}^{\text{crossing
}S_{k}}\frac{m}{\left|v_{k}^{i}\right|}\right>,$ (15)

where $v_{k}^{i}$ denotes the velocity component in the $k$-direction of
particle $i$. Similarly, regarding the macroscopic mass flux $\rho u_{l}$
given by

$\displaystyle\rho u_{l}(\mbox{\boldmath$x$},t)=$
$\displaystyle\int\\!\\!\\!\int\\!\\!\\!\int_{-\infty}^{\infty}\mathrm{d}\mbox{\boldmath$v$}f(\mbox{\boldmath$x$},\mbox{\boldmath$v$},t)v_{l},$
(16)

substituting $\xi(S_{k},t)$ and $\xi^{i}$ in Eq. (14) with
$\frac{v_{l}}{|v_{k}|}$ and $\frac{v_{l}^{i}}{|v_{k}^{i}|}$, respectively,
leads to

$\rho u_{l}(S_{k},t)=\lim_{\delta t\to 0}\frac{1}{S_{k}\delta
t}\left<\sum_{i\in\mathrm{fluid},\delta t}^{\text{crossing
}S_{k}}\frac{mv_{l}^{i}}{\left|v_{k}^{i}\right|}\right>.$ (17)

From Eqs. (15) and (17), the macroscopic velocity $u_{l}$ results in

$u_{l}(S_{k},t)=\frac{\rho u_{l}(S_{k},t)}{\rho(S_{k},t)}=\lim_{\delta t\to
0}\frac{\displaystyle\left<\sum_{i\in\mathrm{fluid},\delta t}^{\text{crossing
}S_{k}}\frac{mv_{l}^{i}}{\left|v_{k}^{i}\right|}\right>}{\displaystyle\left<\sum_{i\in\mathrm{fluid},\delta
t}^{\text{crossing }S_{k}}\frac{m}{\left|v_{k}^{i}\right|}\right>}.$ (18)

Finally, to write the kinetic contribution of the stress
$\tau_{kl}^{\text{kin}}$, we use the expression in the kinetic theory of gases
given by

$\displaystyle\tau_{kl}^{\text{kin}}(\mbox{\boldmath$x$},t)=-\int\\!\\!\\!\int\\!\\!\\!\int_{-\infty}^{\infty}\mathrm{d}\mbox{\boldmath$v$}f(\mbox{\boldmath$x$},\mbox{\boldmath$v$},t)\left(v_{k}-u_{k}(\mbox{\boldmath$x$},t)\right)\left(v_{l}-u_{l}(\mbox{\boldmath$x$},t)\right).$
(19)

By expanding Eq. (19), it follows

$\displaystyle\tau_{kl}^{\text{kin}}(\mbox{\boldmath$x$},t)=$
$\displaystyle-\int\\!\\!\\!\int\\!\\!\\!\int_{-\infty}^{\infty}\mathrm{d}\mbox{\boldmath$v$}f(\mbox{\boldmath$x$},\mbox{\boldmath$v$},t)v_{k}(v_{l}-u_{l}(\mbox{\boldmath$x$},t))+u_{k}(\mbox{\boldmath$x$},t)\int\\!\\!\\!\int\\!\\!\\!\int_{-\infty}^{\infty}\mathrm{d}\mbox{\boldmath$v$}f(\mbox{\boldmath$x$},\mbox{\boldmath$v$},t)(v_{l}-u_{l}(\mbox{\boldmath$x$},t))$
$\displaystyle=$
$\displaystyle-\int\\!\\!\\!\int\\!\\!\\!\int_{-\infty}^{\infty}\mathrm{d}\mbox{\boldmath$v$}f(\mbox{\boldmath$x$},\mbox{\boldmath$v$},t)v_{k}(v_{l}-u_{l}(\mbox{\boldmath$x$},t))+u_{k}\rho
u_{l}(\mbox{\boldmath$x$},t)-u_{k}\rho u_{l}(\mbox{\boldmath$x$},t)$
$\displaystyle=$
$\displaystyle-\int\\!\\!\\!\int\\!\\!\\!\int_{-\infty}^{\infty}\mathrm{d}\mbox{\boldmath$v$}f(\mbox{\boldmath$x$},\mbox{\boldmath$v$},t)v_{k}(v_{l}-u_{l}(\mbox{\boldmath$x$},t)).$
(20)

Hence, by substituting $\xi(S_{k},t)$ and $\xi^{i}$ in Eq. (14) with
$-\frac{v_{k}(v_{l}-u_{l})}{|v_{k}|}$ and
$-\frac{v_{k}^{i}(v_{l}^{i}-u_{l})}{|v_{k}^{i}|}$, respectively, it follows

$\displaystyle\tau^{\text{kin}}_{kl}(S_{k},t)$ $\displaystyle=\lim_{\delta
t\to 0}\left[-\frac{1}{S_{k}\delta t}\left<\sum_{i\in\mathrm{fluid},\delta
t}^{\text{crossing
}S_{k}}\frac{mv_{k}^{i}\left(v_{l}^{i}-u_{l}(S_{k},t)\right)}{\left|v_{k}^{i}\right|}\right>\right]$
$\displaystyle=\lim_{\delta t\to 0}\left(-\frac{1}{S_{k}\delta
t}\left<\sum_{i\in\mathrm{fluid},\delta t}^{\text{crossing
}S_{k}}\frac{mv_{k}^{i}v_{l}^{i}}{\left|v_{k}^{i}\right|}\right>+u_{l}(S_{k},t)\frac{1}{S_{k}\delta
t}\left<\sum_{i\in\mathrm{fluid},\delta t}^{\text{crossing
}S_{k}}\frac{mv_{k}^{i}}{\left|v_{k}^{i}\right|}\right>\right)$
$\displaystyle=-\lim_{\delta t\to 0}\frac{1}{S_{k}\delta
t}\left<\sum_{i\in\mathrm{fluid},\delta t}^{\text{crossing
}S_{k}}\frac{mv_{k}^{i}v_{l}^{i}}{\left|v_{k}^{i}\right|}\right>+\rho
u_{l}u_{k}(S_{k},t),$ (21)

where Eq. (17) is used in the final equality. Note that the second term in the
rightmost-HS can be obtained by

$\rho u_{l}u_{k}(S_{k},t)=\frac{\rho u_{l}(S_{k},t)\cdot\rho
u_{k}(S_{k},t)}{\rho(S_{k},t)}$ (22)

using Eqs. (17) and (18), which correspond to the advection term in the
macroscopic momentum conservation in the Navier-Stokes equation (5). By
subtracting $\rho u_{l}u_{k}(S_{k},t)$ from the rightmost-HS and leftmost-HS
of Eq. (21), it follows

$\tau^{\text{kin}}_{kl}(S_{k},t)-\rho u_{l}u_{k}(S_{k},t)=-\lim_{\delta t\to
0}\frac{1}{S_{k}\delta t}\left<\sum_{i\in\mathrm{fluid},\delta
t}^{\text{crossing
}S_{k}}\frac{mv_{k}^{i}v_{l}^{i}}{\left|v_{k}^{i}\right|}\right>,$ (23)

meaning that the microscopic total momentum transfer in the RHS corresponds to
the stress minus the advection term in the LHS. Technically, the summation in
the RHS of Eq. (23) is calculated during the MD simulation, and as the post
process, the stress $\tau^{\text{kin}}_{kl}(S_{k},t)$ is obtained by adding
the advection term $\rho u_{l}u_{k}$ to the total microscopic momentum
transfer as

$\displaystyle\tau^{\text{kin}}_{kl}(S_{k},t)$
$\displaystyle=\left[\tau^{\text{kin}}_{kl}(S_{k},t)-\rho
u_{l}u_{k}(S_{k},t)\right]+\rho u_{l}u_{k}(S_{k},t),$ (24)

where the advection term is calculated by the dividing $\rho
u_{l}(S_{k},t)\cdot\rho u_{k}(S_{k},t)$ by the density $\rho(S_{k},t)$ as in
Eq. (22): all obtained also as the post process.

The relation between the macroscopic variables in Eqs. (1) and (5)
corresponding microscopic expressions are summarized in TABLE 1.

Table 1:  Microscopic expressions for the calculation of the corresponding macroscopic properties defined as the average on bin face $S_{k}$ in steady-state systems. The top four properties can be directly calculated from steady-state systems through the MoP procedure, whereas the others below are derived from the four.  macroscopic property | microscopic expression | corresponding equation(s)
---|---|---
$\rho(S_{k},t)$ | $\displaystyle\lim_{\delta t\to 0}\frac{1}{S_{k}\delta t}\left<\sum_{i\in\mathrm{fluid},\delta t}^{\text{crossing }S_{k}}\frac{m}{\left|v_{k}^{i}\right|}\right>$ | Eq. (15)
$\rho u_{l}(S_{k},t)$ | $\displaystyle\lim_{\delta t\to 0}\frac{1}{S_{k}\delta t}\left<\sum_{i\in\mathrm{fluid},\delta t}^{\text{crossing }S_{k}}\frac{mv_{l}^{i}}{\left|v_{k}^{i}\right|}\right>$ | Eq. (17)
$\tau_{kl}^{\mathrm{int}}(S_{k},t)$ | $\displaystyle-\frac{1}{S_{k}}\left<\sum_{(i,j)\in\mathrm{fluid}}^{\text{across\ }S_{k}}F^{ij}_{l}\frac{r^{ij}_{k}}{|r^{ij}_{k}|}\right>$ | Eq. (9)
$\tau_{kl}^{\mathrm{kin}}(S_{k},t)-\rho u_{l}u_{k}(S_{k},t)$ | $\displaystyle-\lim_{\delta t\to 0}\frac{1}{S_{k}\delta t}\left<\sum_{i\in\mathrm{fluid},\delta t}^{\text{crossing }S_{k}}\frac{mv_{k}^{i}v_{l}^{i}}{\left|v_{k}^{i}\right|}\right>$ | Eq. (23)
$\displaystyle u_{l}=\frac{\rho u_{l}}{\rho}$ | - | Eq. (18)
$\rho u_{l}u_{k}$ | - | Eq. (22)
$\tau_{kl}^{\mathrm{kin}}=\left(\tau_{kl}^{\mathrm{kin}}-\rho u_{l}u_{k}\right)+\rho u_{l}u_{k}$ | - | Eq. (24)
$\tau_{kl}=\tau_{kl}^{\mathrm{kin}}+\tau_{kl}^{\mathrm{int}}$ | - | Eq. (7)

In practice, within the framework of MD, $\delta t\ (\to 0)$ must be replaced
by a small non-zero time step of $\Delta t$ for the numerical integration.
Upon this procedure without this limit, we have to assume the following: 1)
the change of the distribution function
$f(\mbox{\boldmath$x$},\mbox{\boldmath$v$},t)$ within the distance range of
$|\mbox{\boldmath$v$}|\delta t$ is negligibly small, and 2) the values of
$v^{i}_{k}$ and $v^{i}_{l}$ upon ‘crossing’ should be properly evaluated based
on the position update procedure of particles depending on the time
integration scheme. For the velocity Verlet method, which is applied in the
numerical test in Sec. III, we adopted
$\mbox{\boldmath$v$}^{i}\equiv\frac{\mbox{\boldmath$x$}^{i}(t+\Delta
t)-\mbox{\boldmath$x$}^{i}(t)}{\Delta t}$ using the positions
$\mbox{\boldmath$x$}^{i}(t)$ and $\mbox{\boldmath$x$}^{i}(t+\Delta t)$ of
fluid particle $i$ at time $t$ and $t+\Delta t$ before and after crossing the
bin face to avoid the discrepancy of the mass flux by the MoP calculation and
by the position update.

Note that Eq. (23) without the limit $\delta t\to 0$ is the same as the RHS of
Eq. (8), which simply sums up the momentum transfer across the bin face
$S_{k}$ with a sign function $\frac{v_{k}^{i}}{|v_{k}^{i}|}$. Hence, if one
locates a control volume with a closed surface consisting of the MoP bin
faces, then the momentum conservation is strictly satisfied with Eq. (23).
Different choices are indeed possible to determine the advection term $\rho
u_{l}u_{k}(S_{k},t)$ in Eq. (22) to separate the stress
$\tau_{kl}^{\mathrm{kin}}(S_{k},t)$ from
$\tau_{kl}^{\mathrm{kin}}(S_{k},t)-\rho u_{l}u_{k}(S_{k},t)$ by Eq. (24), and
this may sound that the definition of $\tau_{kl}^{\mathrm{kin}}(S_{k},t)$ is
not unique. However; by setting $l=k$ in Eq. (17), the surface normal mass
flux is evaluated as the simple sum of the mass passage with a sign function
$\frac{v^{i}_{k}}{|v^{i}_{k}|}$, and this strictly satisfies the mass
conservation in Eq. (2), meaning that one can choose a unique definition of
$\rho u_{l}u_{k}(S_{k},t)$ that simultaneously satisfies the macroscopic mass
and momentum conservation.

Another point to be noted is that the final forms in Eqs. (15), (17) and (23)
are formally equivalent to the MoP expressions by Daivis, Travis, and Todd
[30], which were derived for a quasi-1-dimensional flow through the
expressions of the time derivative of the fluxes in a control volume with the
Fourier transform, and were in principle applicable for the average on an
infinite plane under a periodic boundary condition. On the other hand, our
non-flux-based derivation with a definition of physical properties averaged on
a face through the VDF enables the calculation of physical properties on a
finite area. In addition, taking advantage of this non-flux-based feature, one
can calculate, for instance, the velocity component $u_{l}$ on a bin face
$S_{k}\ (l\neq k)$ tangential to the velocity component by Eq. (18). This
point will be discussed more in detail with a quasi-1D Couette-type flow in
Sec. III.1 as an example.

##  III Numerical Test

The extended MoP was tested through the calculation of the density,
macroscopic mean velocity and stress distributions in two systems with a
Lennard-Jones (LJ) fluid: a quasi-1D Couette-type flow and a quasi-2D shear
flow with solid-liquid-vapor contact lines. Note that both systems are in
steady state and we applied time average instead of ensemble average.

###  III.1 Quasi-1D Couette-type flow

Figure 1 (a) shows the MD simulation system of a quasi-1D Couette-type flow,
where the basic setup is a standard one similar to our previous study.[31, 32]
The two parallel solid walls were fcc crystals and every pair of the nearest
neighbors in the walls was bound through a harmonic potential
$\Phi_{\text{h}}(r)=\frac{k}{2}(r-r_{\text{eq}})^{2}$, with $r$ being the
interparticle distance, $r_{\text{eq}}=0.277$ nm, and $k=46.8$ N/m.
Interactions between fluid particles and between fluid and solid particles
were modeled by a 12-6 LJ potential
$\Phi^{\mathrm{LJ}}(r_{ij})=4\epsilon_{ij}\left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12}-\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right]$,
where $r_{ij}$ was the distance between the particles $i$ and $j$, while
$\epsilon$ and $\sigma$ denoted the LJ energy and length parameters,
respectively. This LJ interaction was truncated at a cut-off distance of
$r_{\mathrm{c}}=3.5\sigma$ and quadratic functions were added so that the
potential and interaction force smoothly vanished at $r_{\mathrm{c}}$. [22] We
used the following parameters for fluid-fluid (ff) and solid-fluid (sf)
interactions: $\sigma_{\text{ff}}=0.340$ nm, $\epsilon_{\text{ff}}=1.67\times
10^{-21}$ J, $\sigma_{\text{sf}}=0.345$ nm, $\epsilon_{\text{sf}}=0.646\times
10^{-21}$ J. The atomic masses of fluid and solid particles were
$m_{\text{f}}=39.95$ u and $m_{\text{s}}=195.1$ u, respectively. Finally, the
equations of motion were integrated using the velocity-Verlet algorithm, with
a time step $\Delta t$ of 5 fs.

The periodic boundary condition was set in the $x$ and $y$-directions, and
4000 LJ particles were confined between two parallel solid walls consisting of
the fcc crystal located on the bottom and top sides of the calculation cell,
which directed (001) and (00$-1$) planes normal to the $z$-direction. Both had
eight layers so that the possible minimum distance between the fluid particle
and the solid particle in the outmost layer was longer than the cutoff
distance. The relative positions of the solid particles in the outmost layers
of each base crystal were fixed and the temperature of those in the second
outermost layers was controlled at a control temperature of 100 K by using the
standard Langevin thermostat. [33] The system was first equilibrated for 10 ns
using the top wall as a piston with a control pressure of 4 MPa without shear
so that a quasi-1D system with a LJ liquid confined between fcc solid walls
was achieved. After the equilibration, further relaxation run to achieve a
steady shear flow was carried out for 10 ns by moving the particles in the
outmost layers of both walls with opposite velocities of $\pm$100 m/s in the
$x$-direction, using the top wall as a piston with a control pressure of 4
MPa. Finally, steady shear flow simulation was carried out, keeping their
$z$-position constant at the average position during the 2nd relaxation run,
where the system pressure resulted in 3.61 MPa.

We tested the MoP expression (TABLE 1) in the steady state, where the local
density, velocity, advection term, and stress were obtained as the time-
average of 200 ns on a grid with $x$-normal bin faces with a height $\Delta
z=0.150$ nm and $z$-normal ones with a width $\Delta x=0.145$ nm. Assuming
that the system is quasi-1D, the distribution in the $x$-direction was
averaged for bins with identical $z$-positions.

Figure 1:  (a) Quasi-1D Couette-type flow system of a Lennard-Jones liquid
confined between two solid walls. (b) Distributions of density $\rho$ and
velocity $u_{x}$ calculated by the proposed Method of Plane (MoP) and the
volume average (VA). Solid and dashed lines denote the results of MoP and VA,
respectively while the two lines almost overlap in this scale. (c) Difference
between MoP and VA regarding density $\rho^{\text{MoP}}-\rho^{\text{VA}}$ and
velocity $u_{x}^{\text{MoP}}-u_{x}^{\text{VA}}$ with their error bars depicted
with semi-transparent areas around the average.

For a comparison, we also obtained the density and velocity distributions
based on a standard volume average (VA), where the time-average in equally
divided bin volumes parallel to the solid wall with a height of 0.15 nm were
calculated.

Figure 1 (b) shows the distributions of density $\rho$ and macroscopic
velocity in the $x$-direction $u_{x}$ calculated by the proposed MoP and
standard VA as a reference. Note that these distributions by the MoP can be
calculated both on $x$-normal bin faces and on $z$-normal ones as shown later,
while only the distributions obtained on $x$-normal bin faces are shown as the
MoP results here. Overall, the MoP well reproduced the results by the VA, and
the two lines almost overlap in this scale. Regarding the density
distribution, except near the walls where layered structures are observed,
bulk liquid with almost constant density was formed. Note that the bulk
density was not completely constant because the temperature was not constant
due to the viscous heat dissipation induced by the extreme shear imposed on
this system. The shear velocity profiles are linear throughout almost all the
liquid part except in layered structures, which can be understood by the
change of local viscosity there. The density and velocity differences between
MoP and VA are shown in Fig. 1 (c). The density difference was within 10
$\text{kg}/\text{m}^{3}$, which is less than 1 % of the bulk density, and the
velocity difference was also within 0.5 m/s, showing that proposed MoP can
extract the density and velocity distribution consistent with VA.

Figure 2:  Distributions of (a) the diagonal stress component
$\tau_{zz}(S_{z})(\equiv\tau_{zz}^{\mathrm{int}}+\tau_{zz}^{\mathrm{kin}})$,
advection term $\rho u_{x}u_{x}(S_{x})$ and $\tau_{xx}(S_{x})-\rho
u_{x}u_{x}(S_{x})[\equiv\tau_{xx}^{\mathrm{int}}+(\tau_{xx}^{\mathrm{kin}}-\rho
u_{x}u_{x})]$, and (b) diagonal and off-diagonal stress components
$\tau_{xx}(S_{x})$, $\tau_{zz}(S_{z})$, $\tau_{zx}(S_{z})$ and
$\tau_{xz}(S_{x})$.

Figure 2 (a) shows the distributions of
$\tau_{zz}(\equiv\tau_{zz}^{\mathrm{int}}+\tau_{zz}^{\mathrm{kin}})$,
$\tau_{xx}-\rho
u_{x}u_{x}[\equiv\tau_{xx}^{\mathrm{int}}+(\tau_{xx}^{\mathrm{kin}}-\rho
u_{x}u_{x})]$ and $\rho u_{x}u_{x}$, where the first two were directly
obtained with simple addition based on Eqs. (9) and (21) as also listed in the
top part of TABLE 1, while $\rho u_{x}u_{x}$ was obtained from the density
$\rho$ and velocity $u_{x}$. Note that $\tau_{zz}-\rho u_{z}u_{z}$ is shown as
$\tau_{zz}$ because $u_{z}$ is equal to zero in the whole area of the present
system. Also note that the calculation of $\tau_{kl}^{\mathrm{int}}$ in Eq.
(9) was the same as in equilibrium systems without macroscopic flow. As
clearly observed, $\tau_{xx}-\rho u_{x}u_{x}$ including the advection and
$\tau_{zz}$ are different away from the solid walls, indicating that the flow
effect should be removed to properly evaluate the fluid stress. Figure 2 (b)
displays the distributions of the stress component $\tau_{xx}$, $\tau_{zz}$,
$\tau_{zx}$ and $\tau_{xz}$, where $\tau_{zz}$ and $\tau_{zx}$ were calculated
on $z$-normal bins whereas the others were obtained on $x$-normal bins. As
explained above, the stress value $\tau_{xx}$ was calculated by adding $\rho
u_{x}u_{x}$ to $\tau_{xx}-\rho u_{x}u_{x}$ whereas the advection terms for the
others can be neglected considering $u_{z}=0$. In the bulk region sufficiently
away from the walls, $\tau_{xx}=\tau_{zz}$ and $\tau_{zx}=\tau_{xz}$ are
satisfied as expected from the solution of a laminar Couette flow, and the
former indicates that the stress $\tau_{xx}$ is adequately calculated by the
proposed MoP with the resulting value $-\tau_{xx}(=-\tau_{zz})$ equal to the
external pressure value of 3.61 MPa. The wall-tangential diagonal stress
$\tau_{xx}$ fluctuates near the walls as typically observed also in
equilibrium systems, [22, 10] because of the layered structure of the liquid
as displayed in the density distribution in Fig. 1 (b). On the other hand,
$\tau_{zz}$ was constant except near the walls, where the solid-liquid (SL)
interaction acts as the external force on the liquid. Regarding the off-
diagonal components $\tau_{zx}(=\tau_{xz})$, they were constant except just
around the walls, where friction from the solid is included in the force
balance even in the laminar flow.

In addition to the normal velocity component $u_{k}$ on the MoP plane $S_{k}$,
the calculation of $u_{l}\ (l\neq k)$ tangentially to $S_{k}$ is needed for
the separation of $\tau_{kl}^{\mathrm{kin}}(S_{k})-\rho u_{l}u_{k}(S_{k})$ in
Eq. (24) to properly define the stress in general flows with $u_{l}\neq 0$ and
$u_{k}\neq 0$. Including this tangential velocity, we compared the
distributions of the density $\rho$, the mass flux $\rho u_{x}$ and the
velocity $u_{x}$ averaged on $x$-normal and $z$-normal bin faces as another
numerical test in the present system in Fig. 1. More concretely, the density
$\rho(S_{x})$ and $\rho(S_{z})$ averaged on $x$-normal and $z$-normal bin
faces $S_{x}$ and $S_{z}$, respectively were calculated by Eq. (15) with
setting $k=x$ and $k=z$, whereas the macroscopic mass flux $\rho u_{x}(S_{x})$
and $\rho u_{x}(S_{z})$ were obtained by Eq. (17) with $l=x$ on $S_{x}$ and
$S_{z}$, respectively. With these definitions, $u_{x}(S_{x})\equiv\frac{\rho
u_{x}(S_{x})}{\rho(S_{x})}$ on $S_{x}$ as well as
$u_{x}(S_{z})\equiv\frac{\rho u_{x}(S_{z})}{\rho(S_{z})}$ on $S_{z}$ can be
obtained as in Eq. (18). Note that in the present laminar flow system with
$u_{z}=0$, the calculation of $u_{x}$ is practically not needed for the stress
separation for $\tau_{zz}^{\mathrm{kin}}$, $\tau_{zx}^{\mathrm{kin}}$ and
$\tau_{xz}^{\mathrm{kin}}$ in Eq. (24). Figure 3 shows the distributions of
the (a) density $\rho$, (b) mass flux $\rho u_{x}$, and (c) velocity $u_{x}$
defined on $x$-normal and $z$-normal bin faces, in the system in Fig. 1, where
the values averaged on each bin face of $x$-normal and $z$-normal are plotted
with setting the $z$-position at the center of each bin face, respectively,
i.e., they are staggered by $\Delta z/2$. In the bulk, $\rho$, $\rho u_{x}$
and resulting $u_{x}$ averaged on bin faces with different normal directions
agreed well, indicating that the separation of the stress and advection terms
in Eq. (24) is possible with the velocity values properly evaluated by the
proposed method. The difference seen around the top and bottom is due to the
layered structures around the two walls, i.e., the values on $S_{z}$ are the
average on a surface parallel to the layered structure whereas those on
$S_{x}$ are the average across the layers.

Figure 3:  Comparison of the time-averaged distributions of the (a) density
$\rho$ (b) mass flux $\rho u_{x}$ and (c) velocity $u_{x}$ averaged on
$x$-normal and $z$-normal bin faces $S_{x}$ and $S_{z}$, respectively.

###  III.2 Quasi-2D shear flow with solid-liquid-vapor contact lines

Figure 4:  Top: Quasi-2D Couette-type flow system of a Lennard-Jones liquid
confined between two solid walls. Middle: Distributions of density $\rho$,
velocity $u$, and off-diagonal stress component $\tau_{zx}$. Black arrow
denotes the macroscopic velocity calculated by the proposed Method of Plane.
Bottom: Distributions of diagonal stress components $\tau_{xx}$ and
$\tau_{zz}$.

The top panel of Fig. 4 shows the MD simulation system of a quasi-2D Couette-
type flow, where the basic setups are the same as in the quasi-1D system. The
periodic boundary condition was set in the $x$\- and $y$-directions, and
20,000 LJ particles were confined between two parallel solid walls with a
distance about 10.4 nm and a dimension of $x\times y=3.92\times 39.2$ nm2 so
that the LJ fluid may form two quasi-2D menisci with contact lines on the
walls upon the preliminary equilibration at a control temperature $T=85$ K
without shear. The static contact angle on both top and bottom walls were
about 57 degrees. After the equilibration, further relaxation run to achieve a
steady shear flow with asymmetric menisci were carried out for 10 ns by moving
the particles in the outmost layers of both walls with opposite velocities of
$\pm$10 m/s in the $x$-direction.

After the relaxation run, the density, velocity and stress distributions were
obtained by the present MoP expression in the steady state with the time-
average of 500 ns on $x$-normal bin faces with a length of $\Delta z=0.149$ nm
and $z$-normal ones with a length of $\Delta x=0.150$ nm.

The middle panel of Fig. 4 shows the distributions of density $\rho$
calculated on the $x$-normal bin faces, velocity vector with components
calculated on each bin face corresponding the component direction, and a
stress component $\tau_{zx}$, where those for $\rho$ and $\tau_{zx}$ are
displayed only for the half of the system with respect to the center of mass
of the fluid considering the symmetry. A clockwise caterpillar like flow is
clearly captured by the present method, where the shear stress $\tau_{zx}$
distribution in the liquid phase shows the non-uniformity of the viscous
stress. The strong tensile stress seen in the $\tau_{zx}$ distribution around
the LV interfaces is due to the LV interfacial tension. The bottom panel of
Fig. 4 shows the distributions of diagonal stress components $\tau_{xx}$ and
$\tau_{zz}$. Layered structures are observed for $\tau_{xx}$ near the SL
interfaces due to the adsorption layers in the density distribution. The
relation between the density layers and stress distribution near the solid
walls are qualitatively the same as in our previous study of static droplet:
[10] negative stress was seen in the adsorption layers, i.e., the adsorption
layers were compressed whereas tensile stress appeared between the layers.

Indeed, these apparent flow features can be qualitatively visualized by
another methods such as atomic stress [24], but the present method provides
the distributions of physical properties defined on a surface establishing a
direct link with the conservation laws for arbitrary local volume as described
in Introduction, and is generally applicable to a wide range of nanoscale
systems with liquid flow. One of our future research targets is dynamic
wetting [34, 35, 36], for which we plan to examine the mechanical balance
exerted on the fluid around a CV set around the moving contact line in Fig. 4.
Through the comparison with the static case, [10] this would enable the
analysis of advancing and receding contact angle from a mechanical point of
view.

##  IV Concluding remarks

In this work, we showed a calculation method of local stress tensor applicable
to non-equilibrium MD systems based on the Method of Plane (MoP). From the
relation between the macroscopic velocity distribution function and the
microscopic molecular passage across a fixed control plane, we derived a basic
equation to connect the macroscopic field variable and the microscopic
molecular variable. Based on the connection, we derived a method to calculate
the basic properties of the macroscopic momentum conservation law including
the density, velocity and momentum flux as well as the interaction and kinetic
terms of the stress tensor defined on a surface with a finite area. Any
component of the streaming velocity can be obtained on a control surface,
which enables the separation of the kinetic momentum flux into the advection
and stress terms in the framework of the MoP. We verified the present method
through the extraction of the density and velocity distributions by volume
average (VA) and the MoP in a quasi-1D steady-state Couette flow system,
seeing that the stress tensor distribution by the MoP satisfies the solution
of a laminar Couette flow in the bulk, indicating that the flow effect, i.e.,
the advection term, was removed to evaluate stress properly. Furthermore, we
showed the density, velocity, and stress tensor distributions by the MoP even
in a quasi-2D steady-state system with a moving contact line. In our method as
opposed to VA, the density, mass and momentum fluxes are defined on a surface,
which is essential to have consistency with the conservation laws in dynamic
systems.

###### Acknowledgements.

H.K., T.O. and Y.Y. are supported by JSPS KAKENHI Grant Nos. JP20J20251,
JP18K03929, and JP18K03978, Japan, respectively. Y.Y. is also supported by JST
CREST Grant No. JPMJCR18I1, Japan.

DATA AVAILABILITY

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

## References

  * Irving and Kirkwood [1950] J. H. Irving and J. G. Kirkwood, “The statistical mechanical theory of transport processes. IV. The equations of hydrodynamics,” J. Chem. Phys. 18, 817–829 (1950).
  * Alder and Wainwright [1957] B. J. Alder and T. E. Wainwright, “Phase Transition for a Hard Sphere System,” J. Chem. Phys. 27, 1208–1209 (1957).
  * Alder and Wainwright [1959] B. J. Alder and T. E. a. Wainwright, “Studies in Molecular Dynamics. I. General Method,” J. Chem. Phys. 31, 459–466 (1959).
  * Evans and Morriss [2008] D. Evans and G. Morriss,  _Statistical Mechanics of Nonequilibrium Liquids_, 2nd ed. (Cambridge University Press, 2008) pp. 71–72.
  * P. and J. [2017] A. M. P. and T. D. J.,  _Computer Simulation of Liquids_ , 2nd ed. (Oxford, 2017).
  * Tsai [1979] D. H. Tsai, “The virial theorem and stress calculation in molecular dynamics,” J. Chem. Phys. 70, 1375–1382 (1979).
  * Thompson  _et al._ [1984] S. M. Thompson, K. E. Gubbins, J. P. R. B. Walton, R. A. R. Chantry, and J. S. Rowlinson, “A molecular dynamics study of liquid drops,” J. Chem. Phys. 81, 530–542 (1984).
  * Todd, Evans, and Daivis [1995] B. D. Todd, D. J. Evans, and P. J. Daivis, “Pressure tensor for inhomogeneous fluids,” Phys. Rev. E 52, 1627–1638 (1995).
  * Hardy [1982] R. J. Hardy, “Formulas for determining local properties in molecular‐dynamics simulations: Shock waves,” J. Chem. Phys. 76, 622–628 (1982).
  * Yamaguchi  _et al._ [2019] Y. Yamaguchi, H. Kusudo, D. Surblys, T. Omori, and G. Kikugawa, “Interpretation of Young’s equation for a liquid droplet on a flat and smooth solid surface: Mechanical and thermodynamic routes with a simple Lennard-Jones liquid,” J. Chem. Phys. 150, 044701 (2019).
  * Lutsko [1988] J. F. a. Lutsko, “Stress and elastic constants in anisotropic solids: Molecular dynamics techniques,” Journal of Applied Physics 64, 1152–1154 (1988).
  * Weng  _et al._ [2000] J.-G. Weng, S. Park, J. R. Lukes, and C.-L. Tien, “Molecular dynamics investigation of thickness effect on liquid films,” J. Chem. Phys. 113, 5917–5923 (2000).
  * Cormier, Rickman, and Delph [2001] J. Cormier, J. M. Rickman, and T. J. Delph, “Stress calculation in atomistic simulations of perfect and imperfect solids,” Journal of Applied Physics 89, 99–104 (2001).
  * Heyes  _et al._ [2011] D. M. Heyes, E. R. Smith, D. Dini, and T. A. Zaki, “The equivalence between volume averaging and method of planes definitions of the pressure tensor at a plane,” J. Chem. Phys. 135 (2011), 10.1063/1.3605692.
  * Yang, Wu, and Li [2012] J. Z. Yang, X. Wu, and X. Li, “A generalized Irving-Kirkwood formula for the calculation of stress in molecular dynamics models,” J. Chem. Phys. 137 (2012), 10.1063/1.4755946, 84870207982 .
  * Heyes, Dini, and Smith [2016] D. M. Heyes, D. Dini, and E. R. Smith, “Equilibrium fluctuations of liquid state static properties in a subvolume by molecular dynamics,” J. Chem. Phys. 145, 104504 (2016).
  * Smith, Heyes, and Dini [2017] E. R. Smith, D. M. Heyes, and D. Dini, “Towards the Irving-Kirkwood limit of the mechanical stress tensor,” J. Chem. Phys. 146, 224109 (2017).
  * Shi, Santiso, and Gubbins [2021] K. Shi, E. E. Santiso, and K. E. Gubbins, “Can we define a unique microscopic pressure in inhomogeneous fluids?” J. Chem. Phys. 154, 84502 (2021).
  * Smith and Braga [2020] E. R. Smith and C. Braga, “Hydrodynamics across a fluctuating interface,” J. Chem. Phys. 153, 134705 (2020).
  * Bakker [1928] G. Bakker,  _Kapillarität und Oberflächenspannung_ , Vol. 6 (Wien-Harms, 1928).
  * Rowlinson and Widom [1982] J. S. Rowlinson and B. Widom,  _Molecular Theory of Capillarity_ (Dover, 1982).
  * Nishida  _et al._ [2014] S. Nishida, D. Surblys, Y. Yamaguchi, K. Kuroda, M. Kagawa, T. Nakajima, and H. Fujimura, “Molecular dynamics analysis of multiphase interfaces based on in situ extraction of the pressure distribution of a liquid droplet on a solid surface,” J. Chem. Phys. 140, 074707 (2014).
  * Surblys  _et al._ [2014] D. Surblys, Y. Yamaguchi, K. Kuroda, M. Kagawa, T. Nakajima, and H. Fujimura, “Molecular dynamics analysis on wetting and interfacial properties of water-alcohol mixture droplets on a solid surface,” J. Chem. Phys. 140, 034505 (2014).
  * Thompson, Plimpton, and Mattson [2009] A. P. Thompson, S. J. Plimpton, and W. Mattson, “General formulation of pressure and stress tensor for arbitrary many-body interaction potentials under periodic boundary conditions,” J. Chem. Phys. 131, 1–6 (2009).
  * Plimpton [1995] S. Plimpton, “Fast Parallel Algorithms for Short-Range Molecular Dynamics,” J. Comp. Phys. 117, 1–19 (1995).
  * Kusudo, Omori, and Yamaguchi [2019] H. Kusudo, T. Omori, and Y. Yamaguchi, “Extraction of the equilibrium pinning force on a contact line exerted from a wettability boundary of a solid surface through the connection between mechanical and thermodynamic routes,” J. Chem. Phys. 151, 154501 (2019).
  * Imaizumi  _et al._ [2020] Y. Imaizumi, T. Omori, H. Kusudo, C. Bistafa, and Y. Yamaguchi, “Wilhelmy equation revisited: A lightweight method to measure liquid–vapor, solid–liquid, and solid–vapor interfacial tensions from a single molecular dynamics simulation,” J. Chem. Phys. 153, 034701 (2020).
  * Andersen [1983] H. C. Andersen, “Rattle: A ”velocity” version of the shake algorithm for molecular dynamics calculations,” J. Comput. Phys. 52, 24–34 (1983).
  * Surblys  _et al._ [2019] D. Surblys, H. Matsubara, G. Kikugawa, and T. Ohara, “Application of atomic stress to compute heat flux via molecular dynamics for systems with many-body interactions,” Phys. Rev. E 99, 1–6 (2019).
  * Daivis, Travis, and Todd [1996] P. J. Daivis, K. P. Travis, and B. D. Todd, “A technique for the calculation of mass, energy, and momentum densities at planes in molecular dynamics simulations,” J. Chem. Phys. 104, 9651–9653 (1996).
  * Ogawa  _et al._ [2019] K. Ogawa, H. Oga, H. Kusudo, Y. Yamaguchi, T. Omori, S. Merabia, and L. Joly, “Large effect of lateral box size in molecular dynamics simulations of liquid-solid friction,” Phys. Rev. E 100, 023101 (2019).
  * Oga  _et al._ [2019] H. Oga, Y. Yamaguchi, T. Omori, S. Merabia, and L. Joly, “Green-Kubo measurement of liquid-solid friction in finite-size systems,” J. Chem. Phys. 151, 054502 (2019).
  * Blömer and Beylich [1999] J. Blömer and A. E. Beylich, “Molecular dynamics simulation of energy accommodation of internal and translational degrees of freedom at gas–surface interfaces,” Surf. Sci. 423, 127–133 (1999).
  * Hizumi  _et al._ [2015] Y. Hizumi, T. Omori, Y. Yamaguchi, and T. Kajishima, “Study on the Navier boundary condition for flows with a moving contact line by means of molecular dynamics simulation,” Trans. JSME (in Japanese) 81, 15–00409 (2015).
  * Omori  _et al._ [2019] T. Omori, Y. Kobayashi, Y. Yamaguchi, and T. Kajishima, “Understanding the asymmetry between advancing and receding microscopic contact angles,” Soft Matter 15, 3923–3928 (2019).
  * Thalakkttor and Mohseni [2020] J. J. Thalakkttor and K. Mohseni, “Role of the rate of surface dilatation in determining microscopic dynamic contact angle,” Phys. Fluids 32, 012111 (2020).


# Topological waves guided by a glide-reflection symmetric crystal interface

Julio Andrés Iglesias Martínez, Nicolas Laforge, Muamer Kadic, and Vincent
Laude  Institut FEMTO-ST, CNRS UMR 6174, Université Bourgogne Franche-Comté,
Besançon, France

###### Abstract

A domain wall separating two different topological phases of the same crystal
can support the propagation of backscattering-immune guided waves. In valley-
Hall and quantum-Hall crystal waveguides, this property stems from symmetry
protection and results from a topological transition at a Dirac point. Since
an initially closed band gap has to open, the guidance bandwidth remains
limited compared to that of wide band gap crystals. When a glide-symmetric
dislocation is introduced in a 2D crystal, we show that a pair of wide-
bandwidth, single-mode, and symmetry-protected guided waves appear in the bulk
band gap. The 2D Zak phase changes by $\pi$ on either side of the interface,
providing a topological invariant protected by glide-reflection symmetry at
the X point of the Brillouin zone. A demonstration experiment is performed
with acoustic waves in water, at ultrasonic frequencies, and shows the
continuous tuning of transmission as a function of the glide parameter. The
concept further extends to other types of waves, including the case of elastic
waves in solids, but also of optical and electromagnetic waves.

##  I Introduction

Topological phononics promises unprecedented wave properties inspired by the
concepts of topological insulators Hasan and Kane (2010); Fleury  _et al._
(2014); Yang  _et al._ (2015); Zhang  _et al._ (2018); Ma  _et al._ (2019);
Gao and Christensen (2021). One promising direction is the achievement of uni-
directional and backscattering-free guided wave propagation along a boundary
of a crystal or a domain wall between two crystal phases. Passive topological
waveguides, for instance of the valley-Hall Lu  _et al._ (2017); Yan  _et al._
(2018) and the quantum-Hall type Wu and Hu (2015); He  _et al._ (2016);
Miniaci  _et al._ (2018), lead to symmetry-protected, single-mode guided waves
along a domain wall separating two phases of the same crystal with different
topological invariants. The topological properties of the waveguide are
inherited from those of the two-dimensional (2D) bulk crystal according to the
bulk-boundary correspondence principle Mong and Shivamoggi (2011). The one-
dimensional (1D) domain wall hetero-structure is formed without breaking the
periodicity of the 2D lattice, by tuning a continuous geometrical parameter
that controls a topological transition Lu  _et al._ (2014); Chen  _et al._
(2019). For instance, in valley-Hall crystals, a triangular inclusion is
rotated continuously to reduce the symmetry of a 2D crystal possessing a band
structure with a Dirac point at the K point of the first Brillouin zone,
causing a gap to open there Lu  _et al._ (2016, 2017); Zhu  _et al._ (2018).
In quantum-Hall crystals with C6v symmetry, a double Dirac at the $\Gamma$
point undergoes a topological transition under the continuous tuning of the
internal structure of the unit cell of the crystal Wu and Hu (2015); He  _et
al._ (2016); Pal and Ruzzene (2017). In both cases, however, the available
bandwidth for the dispersion of the guided wave is limited by the effective
opening of the band gap that the control parameter allows Laude (2021). In
contrast, artificial crystals have long been designed to present very wide
complete band gaps Vasseur  _et al._ (2001); Laude (2020); Iglesias Martínez
_et al._ (2021) that the guided bands could in principle cover. Phononic
crystal waveguides formed by coupling a sequence of crystal defects, however,
lack topological protection and are generally multimodal, leading to a
competition of the guided bands inside the complete band gap that can severely
flatten the guided bands Laude (2021).

Can we obtain topological crystal waveguides that make full use of a wide
complete band gap crystal? We propose in this Letter to start from a
2-periodic crystal and to introduce a glide-reflection (GR) symmetric
dislocation running all along one of the periodicity axes. The resulting
structure loses one periodicity, along the directional orthogonal to the glide
operation, but gains a glide-reflection symmetry that the initial 2-periodic
crystal did not possess. The 2D Zak phase of bulk bands, measured along the
interface direction, changes by $\pi$ on either side, providing a topological
invariant protected by GR symmetry. All pairs of bands at the boundary of the
first Brillouin zone (X point) are degenerate, leading to the appearance of
pairs of left- and right-propagating guided waves in all Bragg band gaps of
the crystal. The pairs of guided Bloch waves are protected by GR symmetry at
the X point and their smooth dispersion covers most of the band gap of the
bulk crystal.

The Letter is organized as follows. We first discuss the topology of the band
structure of a square-lattice crystal and its transformation under a glide
dislocation. We show how the dispersion of waves guided along the glide
dislocation closes the complete band gap exactly for a half-lattice glide
dislocation. Tuning the glide parameter, the spectral transmission can be
changed continuously from no transmission at all to full transmission through
the phononic band gap. The glide-reflection symmetric crystal waveguide offers
wide bandwidth, single mode operation, and symmetry-protected backscattering
immunity. An experiment performed with ultrasonic acoustic waves around 0.5
MHz and a crystal of steel rods in water demonstrates the operation of the
glide-reflection symmetric phononic crystal waveguide.

Figure 1: Principle of the glide-reflection symmetric topological phononic
crystal waveguide. (a) A 2-periodic square-lattice phononic crystal is
composed of steel rods in water (lattice constant $a$, diameter $d=0.9a$). For
every frequency within the complete band gap, transmission of incident
acoustic waves is forbidden. (b) Two pieces of the same square-lattice crystal
are now glided along the $x$-axis. The glide parameter $g$ is periodic with
period $a$. For a half-lattice glide parameter ($g=a/2$), waves are guided
along the glide dislocation, for frequencies within the complete band gap. In
numerical simulations, $P$ is the normalized pressure field, frequency is
taken at the center of the band gap, and waves are incident from the left.

For demonstration purposes, we consider in the following a 2-periodic square-
lattice phononic crystal of circular inclusions, as shown in Figure 1.
Numerical simulations are performed in this Letter considering a phononic
crystal of steel rods in water, but the results extend naturally to other
material systems [SM]See Supplemental Material at link_to_be_defined for more
information on other material systems, for a discussion of modal symmetry, for
experimental details, and for exploration of resilience to disorder.. In the
glide dislocation, one half of the crystal is spatially shifted by an amount
$g$, along direction $x$. The glide operation creates an interface between two
crystal phases that remain identical except for the spatial shift. The initial
crystal ($g=0$) possesses a complete band gap within which transmission
decreases exponentially with crystal thickness. Fig. 1(a) illustrates
numerically, for a frequency at the center of the complete band gap, that
total reflection of incident waves results. When $g=a/2$, guided waves appear
along the dislocation and transmission is obtained, as Fig. 1(b) shows.

Let us analyze the topology of the band structure of the phononic crystal
structure and its change with the glide parameter. Figure 2(a) shows a super-
cell and the phononic band structure for the 2-periodic crystal for $g=0$;
Fig. 2(b) shows similar information for the glide-reflection symmetric
waveguide for $g=a/2$. The super-cell is the numerical device used to obtain
the dispersion relation of crystal waveguides. Periodic boundary conditions
are applied on the vertical boundaries of the supercell while the horizontal
boundaries are left free (Neumann boundary condition). In the figure, we
consider $N=10$ unit cells in the vertical direction and $1$ unit cell in the
horizontal direction. The band structures account for acousto-elastic coupling
between pressure waves in water and elastic waves in steel Laforge  _et al._
(2021). The band structure for $g=0$ in Fig. 2(a) shows the complete band gap
separated by groups of bands. Counting the bands, there are exactly $N+1$
bands below the band gap. Those bands are actually sampled from the original
Brillouin zone as $(k_{y}a/\pi=n/N,k_{x})$ with $n=0,\cdots,N$ Laude (2020).
Hence, when $n=N$, the Bloch wavevector varies along the YM edge of the first
Brillouin zone ($k_{y}a/\pi=1$).

Figure 2: Band structure topology of the glide-reflection waveguide, computed
for a supercell made of $N=10$ unit cells of the crystal. (a) For $g=0$, the
supercell simply repeats $N$ times vertically the primitive cell of the 2D
square-lattice crystal ($d=2$ mm, $a=2.22$ mm). The band structure of the
waveguide is obtained from the projected band structure of the 2D crystal. (b)
For $g=a/2$, the bands group by pairs of symmetric (S) / anti-symmetric (AS)
Bloch waves, with respect to the glide reflection symmetry. They are
degenerate at the Brillouin zone edge (the X point), causing a pair of guided
waves to appear inside the complete band gap. The modal distributions of the S
(red color band) and AS (yellow color band) guided waves are shown on the left
for $k_{x}a/\pi=0.8$.

The guided waves for $g\neq 0$ appearing inside the band gap originate from
the $N$-th and $(N+1)$-th bands. Actually, as Fig. 2(b) shows for glide
parameter $g=a/2$, all bands are degenerate by pairs at the X point of the
first Brillouin zone. This essential property is obtained only for a half-
lattice glide; it is shown later to signal a topological transition of the
band structure occurring at $g=a/2$. Since the $N$-th and the $(N+1)$-th bands
were repelling and thus sitting on opposite sides on the band gap for $g=0$
and they are degenerate at the X point for $g=a/2$, they have to move inside
the band gap as $g$ is tuned continuously between those two values. When
$g>a/2$ and is tuned toward $g=a$, the gap closes continuously and
symmetrically from the case $g<a/2$.

Why the $N$-th and $(N+1)$-th bands hold a pair a guided Bloch waves can be
understood based upon the transformation of the band structure under the
continuous change of glide parameter $g$ from $0$ to $a/2$. The dispersion of
the guided wave extends inside the complete band gap, with a real wavevector
$k_{x}$; along the $y$ direction the guided wave is evanescent, i.e. its
amplitude is decreasing exponentially. Fig. 2(b) illustrates that property for
one particular value of $k_{x}$. In contrast to the case $g=0$, the Bloch
waves of the 2-periodic phononic crystal do not translate directly into Bloch
waves of the waveguide structure for $g\neq 0$. However, they can still be
used as a functional basis to express the $1$-periodic guided mode. Bloch
waves of the 2-periodic crystal are all evanescent for frequencies inside the
complete band gap. Hence their wavevector satisfies $\Re(k_{y})a/\pi=1$: the
real part of the Bloch wavevector is restricted to the top edge of the 2D
first Brillouin zone. Furthermore, the imaginary part of the Bloch wavevector
can only be directed along the $k_{y}$ direction in reciprocal space, since
propagation is lossless along the $x$-axis. $\Im(k_{y})\neq 0$ hence provides
the necessary exponential decrease away from the glide interface such that the
guided wave is confined. It can be check visually in the modal shapes of Fig.
2(b) that the exponential decrease in amplitude along the $y$-axis is
accompanied by an alternation of the sign from one unit cell to the next, in
accordance with the condition $\Re(k_{y})a/\pi=1$.

Zak phase was originally introduced Zak (1989) for 1D crystals as the integral
of the Berry connection along the 1D Brillouin zone. The 2D Zak phase for 2D
crystals Liu and Wakabayashi (2017); Liu  _et al._ (2018) is a natural
generalization where the integral of the Berry connection is taken along a 1D
contour, chosen as the interface direction in reciprocal space Ma  _et al._
(2019). Namely, Zak phase for band $n$ is

$\gamma_{n}=\int_{C}\rm{d}{\bf R}\cdot{\cal A}_{n}({\bf R}),$ (1)

with ${\bf R}=k_{x}$ at fixed $k_{y}$ and $C=[-\pi/a;\pi/a]$. The Berry
connection is

${\cal A}_{n}({\bf R})=i\langle u_{n}({\bf R})|\nabla_{\bf{R}}|u_{n}({\bf
R})\rangle,$ (2)

with $u_{n}({\bf R})$ a Bloch wave defined over the 2D unit cell and
$\langle.\rangle$ denoting the scalar product in real space defined on this
unit cell. Note that the integration contour chosen, $C$, is different from
the one used to define Chern numbers, that is the boundary enclosing the first
2D Brillouin zone. By construction, the bottom crystal B is the glide-
reflection (GR) image of crystal A. The glide operation implies a phase change
for every Bloch wave $\phi(k_{x})=-gk_{x}$ (a translation of the origin by
$g$). Since the Berry connection changes as ${\cal A}_{n}({\bf
R})\rightarrow{\cal A}_{n}({\bf R})-\frac{\partial\phi}{\partial k_{x}}$ under
any phase change, we have

$\gamma_{n}(B)=\gamma_{n}(A)+2\pi g/a.$ (3)

Hence there is a $\pi$ change of the 2D Zak phase across the interface for
every band, for $g=a/2$ exactly. Since the Zak phase is $2\pi$-periodic, its
value alternates by $\pi$ between both crystal images.

Why degenerescence of Bloch waves occurs specifically for $g=a/2$ and at the X
point of the first Brillouin zone results from the combination of the space
group symmetry of the waveguide and of its periodicity along the glide
dislocation, as the Supplemental Material (SM) details [SM]. Let us consider
here a compact demonstration based on operators of the 1D crystal interface.
For any glide parameter $g$, the composition $G_{a-g}\circ G_{g}$ is the
translation by one lattice constant $T_{a}$ in direct space. In reciprocal
space, this implies $G_{a-g}(k)G_{g}(k)=\exp(ika)$. For $g=a/2$, we then have
$G^{2}_{g}(k)=-1$ at the X point of the 1D Brillouin zone ($ka=\pi$). Hence
the eigenvalues of $G_{g}(\pi/a)$ are $\pm i$. Its eigenvectors form complex
conjugate pairs, since $G_{g}(\pi/a)u=iu$ implies $G_{g}(\pi/a)u^{*}=-iu^{*}$.
Since the glide operator commutes with the dynamical operator for the wave
equation, they share common eigenvectors. Hence $Du=\omega^{2}u$ implies
$Du^{*}=\omega^{2}u^{*}$ since the wave equation has real coefficients because
of its time reversal invariance (TRI). Summarizing, each complex conjugate
eigenvector pair shares a degenerate eigenvalue at the X point.

Figure 3: Glide-reflection symmetric topological phononic crystal waveguide
experiment. (a) Two pieces of the same square-lattice crystal of steel rods in
water (diameter $d=2$ mm; lattice constant $a=2.22$ mm) are glided. The
transmission of acoustic waves guided along the dislocation is probed using an
ultrasonic emitter (E) of short pulses that are detected by an ultrasonic
receiver (R). For a half-lattice glide parameter ($g=a/2$), the (b) numerical
and (c) experimental acoustic wave transmission covers most of the complete
phononic band gap (highlighted with the grey color). The reference for
transmission is the measurement in the absence of the phononic crystal
waveguide. Figure 4: Gap opening as a function of the glide parameter. (a) For
$g\neq a/2$, a mini-gap for guided waves opens in the phononic band structure
at the X point of the Brillouin zone. (b) Experiment confirms the opening of
the mini-gap, for $g=a/4$ and $g=3a/8$. (c) The eigenfrequencies of the two
guided waves at the X point vary with the glide parameter (blue line: S waves;
red line: AS waves). For exactly $g=a/2$, the waveguide is glide-reflection
symmetric and the guided wave gap closes. This gap opens symmetrically on
either side of that value.

We now turn to the experimental demonstration of the glide-reflection
symmetric phononic crystal waveguide. The phononic crystal of steel rods in
water is depicted in Figure 3(a). A total of $24\times 16$ parallel rods are
aligned using perforated parallel plates [SM]. The rod diameter is $d=2$ mm
and the lattice constant is $a=2.22$ mm ($d/a=0.9$). Fig. 3(b) shows the
numerical transmission as a function of frequency. The experimental
transmission of Fig. 3(c) is obtained based on the ultrasonic pulse-echo
technique described for instance in Ref. Khelif  _et al._ (2003). The complete
band gap extends from 0.28 MHz to 0.46 MHz whereas the guided mode
transmission covers the range from 0.28 MHz to 0.43 MHz, in accordance with
theory.

For other glide parameter values, transmission is observed as well inside the
complete band gap but within a reduced frequency range. Actually, for $g\neq
a/2$ the left and right propagating guided Bloch waves interfere and form an
anti-crossing and thus a mini band gap at the X point, as Fig. 4(a) shows. The
experiments reported in Fig. 4(b) clearly show the mini band gap opening when
$g=a/4$ and $g=3a/8$ and the corresponding reduction of the transmission
range. Furthermore, Fig. 4(c) shows the variation with $g$ of the $N$-th and
the $(N+1)$-th band intersection with the X point. The phononic band gap is
completely opened for $g=0$ and $g=a$, closes midway for $g=a/2$, and varies
continuously and symmetrically between these points. Since the gap is fully
opened for $g=0$ and closed for $g=a/2$, and the glide parameter can be
continuously tuned with periodicity $a$, a continously-tunable transmission
filter is obtained. Note that symmetry protection against backscattering of
the guided waves is only achieved when the glide parameter $g=a/2$.

Summarizing the results above, the interface waves are protected by a class of
topology relying on spatial symmetries and thus belong to crystalline
topological insulators Prodan and Schulz-Baldes (2016). Whereas crystalline
topological phases generally induce interface waves that have a gapped
spectrum, because the interface breaks the corresponding spatial symmetries,
the glide symmetry of the interface ensures a gapless Dirac point at the X
point of the Brillouin zone.

The importance of glide-reflection symmetry is further verified in the SM [SM]
by considering the oblique lattice instead of the square lattice. It is
specifically found that inversion symmetry Li  _et al._ (2020) combined with
the glide operation leads to a gapped spectrum, unlike GRS. We note that
glide-reflection symmetric waveguides have been considered before, e.g. for
microwaves Quevedo-Teruel  _et al._ (2021) or acoustic waves Janković and Alù
(2021), but that waveguiding is in this case ensured by structural boundaries
rather than by a phononic band gap. The existence of a complete band gap
without glide ($g=0$) is indeed essential to our result. Defect photonic
crystal waveguides possessing glide-reflection symmetry have also been
discussed in terms of symmetry Mock  _et al._ (2010) or low group velocity
Patil  _et al._ (2022), but without consideration of their topological
properties.

The SM [SM] further explores the resilience of the interface waves to
disorder, a direct check of symmetry protection. It is observed numerically
that they survive a position disorder of at least 5% of the lattice constant
and an inclusion diameter disorder of 10%.

The discussion so far has been limited to the scalar case of acoustic waves
and to a square lattice crystal of steel rods in water. It is obvious,
however, that the symmetry principles involved extend the existence of glide-
reflection symmetric crystal waveguides to other material systems and
lattices. The SM specifically illustrates the cases of acoustic waves in a
fluid with rigid inclusions and of vector elastic waves in a solid perforated
with cylindrical holes or containing solid inclusions [SM]. By virtue of the
well-established analogies between acoustic/elastic waves and
optical/electromagnetic waves Kushwaha  _et al._ (1993), the transposition to
photonic crystals is straightforward. Other wave systems such as plasmonic
crystals, gravity-capillary waves at the surface of water Laforge  _et al._
(2019), or solutions of the Schrödinger equation are likely to present similar
properties too. Indeed, one can start from any wide band gap artificial
crystal and produce a glide-reflection symmetric interface within it. Then,
when the glide parameter is exactly one half of the lattice constant, a
gapless spectrum is obtained whereas the non-glided crystal has a completely
gapped spectrum. Hence, any efficient artificial crystal that has been
designed in the past, from seismic waves at the meter scale Brûlé  _et al._
(2014) to phononic crystals for thermal transport control at the nanometer
scale Yang  _et al._ (2014), can be used as a starting basis to design a
completely new topological glided structure supporting symmetry protected edge
states.

On the practical side, the glide operation offers the opportunity to design a
continuously varying transmission that can be changed from no transmission
(for $g=0$) to full transmission through the phononic band gap (for $g=a/2$).
As a waveguide for transmission of information, the glide-reflection symmetric
crystal waveguide offers wide bandwidth, single mode operation, and symmetry-
protected backscattering immunity.

## Acknowledgments

The authors are grateful for support by the EIPHI Graduate School (Contract
No. ANR-17-EURE-0002).

## References

  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, “Colloquium: Topological insulators,” Rev. Modern Phys. 82, 3045–3067 (2010).
  * Fleury  _et al._ (2014) Romain Fleury, Dimitrios L Sounas, Caleb F Sieck, Michael R Haberman, and Andrea Alù, “Sound isolation and giant linear nonreciprocity in a compact acoustic circulator,” Science 343, 516–519 (2014).
  * Yang  _et al._ (2015) Zhaoju Yang, Fei Gao, Xihang Shi, Xiao Lin, Zhen Gao, Yidong Chong, and Baile Zhang, “Topological acoustics,” Phys. Rev. Lett. 114, 114301 (2015).
  * Zhang  _et al._ (2018) Xiujuan Zhang, Meng Xiao, Ying Cheng, Ming-Hui Lu, and Johan Christensen, “Topological sound,” Commun. Phys. 1, 1–13 (2018).
  * Ma  _et al._ (2019) Guancong Ma, Meng Xiao, and Che Ting Chan, “Topological phases in acoustic and mechanical systems,” Nat. Rev. Phys. 1, 281–294 (2019).
  * Gao and Christensen (2021) Penglin Gao and Johan Christensen, “Topological vortices for sound and light,” Nat. Nanotechnol. 16, 487–489 (2021).
  * Lu  _et al._ (2017) Jiuyang Lu, Chunyin Qiu, Liping Ye, Xiying Fan, Manzhu Ke, Fan Zhang, and Zhengyou Liu, “Observation of topological valley transport of sound in sonic crystals,” Nat. Phys. 13, 369 (2017).
  * Yan  _et al._ (2018) Mou Yan, Jiuyang Lu, Feng Li, Weiyin Deng, Xueqin Huang, Jiahong Ma, and Zhengyou Liu, “On-chip valley topological materials for elastic wave manipulation,” Nat. Mater. 17, 993–998 (2018).
  * Wu and Hu (2015) L.-H. Wu and X. Hu, “Scheme for achieving a topological photonic crystal by using dielectric material,” Phys. Rev. Lett. 114, 223901 (2015).
  * He  _et al._ (2016) C. He, X. Ni, H. Ge, X.-C. Sun, Y.-B. Chen, M.-H. Lu, X.-P. Liu, and Y.-F. Chen, “Acoustic topological insulator and robust one-way transport,” Nat. Phys. 12, 3867 (2016).
  * Miniaci  _et al._ (2018) Marco Miniaci, RK Pal, B Morvan, and M Ruzzene, “Experimental observation of topologically protected helical edge modes in patterned elastic plates,” Phys. Rev. X 8, 031074 (2018).
  * Mong and Shivamoggi (2011) Roger S. K. Mong and Vasudha Shivamoggi, “Edge states and the bulk-boundary correspondence in Dirac Hamiltonians,” Phys. Rev. B 83, 125109 (2011).
  * Lu  _et al._ (2014) Jiuyang Lu, Chunyin Qiu, Shengjun Xu, Yangtao Ye, Manzhu Ke, and Zhengyou Liu, “Dirac cones in two-dimensional artificial crystals for classical waves,” Phys. Rev. B 89, 134302 (2014).
  * Chen  _et al._ (2019) Yi Chen, Xiaoning Liu, and Gengkai Hu, “Topological phase transition in mechanical honeycomb lattice,” J. Mech. Phys. Solids , 54 (2019).
  * Lu  _et al._ (2016) Jiuyang Lu, Chunyin Qiu, Manzhu Ke, and Zhengyou Liu, “Valley vortex states in sonic crystals,” Phys. Rev. Lett. 116, 093901 (2016).
  * Zhu  _et al._ (2018) Hongfei Zhu, Ting-Wei Liu, and Fabio Semperlotti, “Design and experimental observation of valley-Hall edge states in diatomic-graphene-like elastic waveguides,” Phys. Rev. B 97, 174301 (2018).
  * Pal and Ruzzene (2017) Raj Kumar Pal and Massimo Ruzzene, “Edge waves in plates with resonators: an elastic analogue of the quantum valley Hall effect,” New J. Phys. 19, 025001 (2017).
  * Laude (2021) Vincent Laude, “Principles and properties of phononic crystal waveguides,” APL Mater. 9, 080701 (2021).
  * Vasseur  _et al._ (2001) J. O. Vasseur, P. A. Deymier, B. Chenni, B. Djafari-Rouhani, L. Dobrzynski, and D. Prevost, “Experimental and theoretical evidence for the existence of absolute acoustic band gaps in two-dimensional solid phononic crystals,” Phys. Rev. Lett. 86, 3012–3015 (2001).
  * Laude (2020) Vincent Laude,  _Phononic Crystals: Artificial Crystals for Sonic, Acoustic, and Elastic Waves_, 2nd ed. (De Gruyter, Berlin, 2020).
  * Iglesias Martínez  _et al._ (2021) Julio Andrés Iglesias Martínez, Johnny Moughames, Gwenn Ulliac, Muamer Kadic, and Vincent Laude, “Three-dimensional phononic crystal with ultra-wide bandgap at megahertz frequencies,” Appl. Phys. Lett. 118, 063507 (2021).
  * (22) SM, See Supplemental Material at https://www.dropbox.com/s/q7887ep2lo4aepk/Glide_SI.pdf?dl=0 for more information on other material systems, for a discussion of modal symmetry, for experimental details, and for exploration of resilience to disorder.
  * Laforge  _et al._ (2021) Nicolas Laforge, Richard Wiltshaw, Richard V Craster, Vincent Laude, Julio Andrés Iglesias Martínez, Guillaume Dupont, Sébastien Guenneau, Muamer Kadic, and Mehul P Makwana, “Acoustic topological circuitry in square and rectangular phononic crystals,” Phys. Rev. Appl. 15, 054056 (2021).
  * Zak (1989) J Zak, “Berry’s phase for energy bands in solids,” Phys. Rev. Lett. 62, 2747 (1989).
  * Liu and Wakabayashi (2017) Feng Liu and Katsunori Wakabayashi, “Novel topological phase with a zero Berry curvature,” Phys. Rev. Lett. 118, 076803 (2017).
  * Liu  _et al._ (2018) Feng Liu, Hai-Yao Deng, and Katsunori Wakabayashi, “Topological photonic crystals with zero Berry curvature,” Phys. Rev. B 97, 035442 (2018).
  * Khelif  _et al._ (2003) A. Khelif, A. Choujaa, B. Djafari-Rouhani, M. Wilm, S. Ballandras, and V. Laude, “Trapping and guiding of acoustic waves by defect modes in a full-band-gap ultrasonic crystal,” Phys. Rev. B 68, 214301 (2003).
  * Prodan and Schulz-Baldes (2016) Emil Prodan and Hermann Schulz-Baldes,  _Bulk and Boundary Invariants for Complex Topological Insulators_ (Springer, 2016).
  * Li  _et al._ (2020) Zheng-wei Li, Xin-sheng Fang, Bin Liang, Yong Li, and Jian-chun Cheng, “Topological interface states in the low-frequency band gap of one-dimensional phononic crystals,” Phys. Rev. Appl. 14, 054028 (2020).
  * Quevedo-Teruel  _et al._ (2021) Oscar Quevedo-Teruel, Qiao Chen, Francisco Mesa, Nelson JG Fonseca, and Guido Valerio, “On the benefits of glide symmetries for microwave devices,” IEEE Journal of Microwaves 1, 457–469 (2021).
  * Janković and Alù (2021) Nikolina Janković and Andrea Alù, “Glide-symmetric acoustic waveguides for extreme sensing and isolation,” Phys. Rev. Appl. 15, 024004 (2021).
  * Mock  _et al._ (2010) Adam Mock, Ling Lu, and John O’Brien, “Space group theory and Fourier space analysis of two-dimensional photonic crystal waveguides,” Phys. Rev. B 81, 155115 (2010).
  * Patil  _et al._ (2022) Chirag Murendranath Patil, Guillermo Arregui, Morten Mechlenborg, Xiaoyan Zhou, Hadiseh Alaeian, Pedro David García, and Søren Stobbe, “Observation of slow light in glide-symmetric photonic-crystal waveguides,” Opt. Express 30, 12565–12575 (2022).
  * Kushwaha  _et al._ (1993) M. S. Kushwaha, P. Halevi, L. Dobrzynski, and B. Djafari-Rouhani, “Acoustic band structure of periodic elastic composites,” Phys. Rev. Lett. 71, 2022–2025 (1993).
  * Laforge  _et al._ (2019) Nicolas Laforge, Vincent Laude, Franck Chollet, Abdelkrim Khelif, Muamer Kadic, Yuning Guo, and Romain Fleury, “Observation of topological gravity-capillary waves in a water wave crystal,” New J. Phys. 21, 083031 (2019).
  * Brûlé  _et al._ (2014) Stéphane Brûlé, EH Javelaud, Stefan Enoch, and Sébastien Guenneau, “Experiments on seismic metamaterials: molding surface waves,” Phys. Rev. Lett. 112, 133901 (2014).
  * Yang  _et al._ (2014) Lina Yang, Nuo Yang, and Baowen Li, “Extreme low thermal conductivity in nanoscale 3D Si phononic crystal with spherical pores,” Nano Lett. 14, 1734–1738 (2014).