# Self-Consistent Determination of Long-Range Electrostatics in Neural Network
Potentials

Ang Gao<EMAIL_ADDRESS>Department of Physics, Beijing University of Posts
and Telecommunications, Beijing, China 100876  Richard C. Remsing
<EMAIL_ADDRESS>Department of Chemistry and Chemical Biology, Rutgers
University, Piscataway, NJ, USA 08854

###### Abstract

Machine learning has the potential to revolutionize the field of molecular
simulation through the development of efficient and accurate models of
interatomic interactions. In particular, neural network models can describe
interactions at the level of accuracy of quantum mechanics-based calculations,
but with a fraction of the cost, enabling the simulation of large systems over
long timescales with ab initio accuracy. However, implicit in the construction
of neural network potentials is an assumption of locality, wherein atomic
arrangements on the scale of about a nanometer are used to learn interatomic
interactions. Because of this assumption, the resulting neural network models
cannot describe long-range interactions that play critical roles in dielectric
screening and chemical reactivity. To address this issue, we introduce the
self-consistent field neural network (SCFNN) model — a general approach for
learning the long-range response of molecular systems in neural network
potentials. The SCFNN model relies on a physically meaningful separation of
the interatomic interactions into short- and long-range components, with a
separate network to handle each component. We demonstrate the success of the
SCFNN approach in modeling the dielectric properties of bulk liquid water, and
show that the SCFNN model accurately predicts long-range polarization
correlations and the response of water to applied electrostatic fields.
Importantly, because of the separation of interactions inherent in our
approach, the SCFNN model can be combined with many existing approaches for
building neural network potentials. Therefore, we expect the SCFNN model to
facilitate the proper description of long-range interactions in a wide-variety
of machine learning-based force fields.

Computer simulations have transformed our understanding of molecular systems
by providing atomic-level insights phenomena of wide importance. The earliest
models used efficient empirical descriptions of interatomic interactions, and
similar force field-based simulations form the foundation of molecular
simulations today Allen and Tildesley (2017). However, it is difficult to
describe processes like chemical reactions that involve bond breakage and
formation, as well as electronic polarization effects within empirical force
fields. The development of quantum mechanics-based ab initio simulations
enabled the description of these complex processes, leading to profound
insights across scientific disciplines Tuckerman  _et al._ (1996); Car and
Parrinello (1985); Chen  _et al._ (2018); Geissler  _et al._ (2001); Lee  _et
al._ (2008); Walker, Crowley, and Case (2008); Senn and Thiel (2006); Dal
Peraro  _et al._ (2007). The vast majority of these first principles
approaches rely on density functional theory (DFT), and the development of
increasingly accurate density functionals has greatly improved the reliability
of ab initio predictions Sun, Ruzsinszky, and Perdew (2015); Sun  _et al._
(2016); Chen  _et al._ (2017); Furness  _et al._ (2020); Zhang  _et al._
(2021a); Adamo and Barone (1999). But, performing electronic structure
calculations are expensive, and first principles simulations are limited to
small system sizes and short time scales.

The prohibitive expense of ab initio simulations can be overcome through
machine learning. Armed with a set of ab initio data, machine learning can be
used to train neural network (NN) potentials that describe interatomic
interactions at the same level of accuracy as the ab initio methods, but with
a fraction of the cost. Consequently, NN potentials enable ab initio quality
simulations to reach the large system sizes and long time scales needed to
model complex phenomena, such as phase diagrams Gartner  _et al._ (2020);
Zhang  _et al._ (2021b); Deringer  _et al._ (2018); Niu  _et al._ (2020);
Deringer  _et al._ (2021) and nucleation Khaliullin  _et al._ (2011); Bonati
and Parrinello (2018).

Despite the significant advances made in this area, there are still practical
and conceptual difficulties with NN potential development, especially with
regard to long-range electrostatics. To make NN potential construction
computationally feasible, most approaches learn only local arrangements of
atoms around a central particle, where the meaning of “local” is defined by a
distance cutoff usually less than 1 nm. Because of this locality, the
resulting NN potentials are inherently short-ranged. The lack of long-range
interactions in NN potentials can lead to both quantitative and qualitative
errors, especially when describing polar and charged species Yue  _et al._
(2021); Grisafi and Ceriotti (2019); Niblett, Galib, and Limmer (2021).

The need for incorporating long-range electrostatics into NN potentials has
led to the development of several new approaches Xie, Persson, and Small
(2020); Yao  _et al._ (2018); Grisafi and Ceriotti (2019); Yue  _et al._
(2021); Ko  _et al._ (2021). Many of these approaches exclude all or some of
the electrostatic interactions from training and then assign effective partial
charges to each atomic nucleus that are used to calculate long-range
electrostatic interactions using traditional methods Xie, Persson, and Small
(2020); Yao  _et al._ (2018); Yue  _et al._ (2021); Ko  _et al._ (2021);
Niblett, Galib, and Limmer (2021). The values of these effective charges can
be determined using machine learning methods. For example, 4G-HDNNP Ko  _et
al._ (2021) employs deep neural networks to predict the electronegativities of
each nucleus, which are subsequently used within a charge equilibration
process to determine the effective charges. These approaches can predict
binding energies and charge transfer between molecules, but they also
introduce quantities that are not direct physical observables, such as the
effective charges and electronegativities. Another approach introduced feature
functions to explicitly incorporate nonlocal geometric information into the
construction of NN potentials Grisafi and Ceriotti (2019). However, these
feature functions depend on the system size. The resulting NN models cannot be
used to accurately model systems that are larger than the original training
set. This size restriction severely limits applicability by making this
approach unable to model extended system sizes.

The difficulties that current approaches to NN potentials have when treating
long-range interactions can be resolved by a purely ab initio strategy that
uses no effective quantities. Such a strategy can be informed by our
understanding of the roles of short- and long-range interactions in condensed
phases Widom (1967); Weeks, Chandler, and Andersen (1971); Chandler, Weeks,
and Andersen (1983); Rodgers and Weeks (2008a). In uniform liquids,
appropriately chosen uniformly slowly-varying components of the long-range
forces — van der Waals attractions and long-range Coulomb interactions —
cancel to a good approximation in every relevant configuration. As a result,
the local structure is determined almost entirely by short-range interactions.
In water, these short-range interactions correspond to hydrogen-bonding and
packing Rodgers and Weeks (2009); Remsing, Rodgers, and Weeks (2011); Rodgers
and Weeks (2008b); Rodgers, Hu, and Weeks (2011). Therefore, short-range
models, including current NN potentials, can describe the structure of uniform
systems. This idea, that short-range forces determine the structure of uniform
systems, forms the foundation for the modern theory of bulk liquids Widom
(1967); Weeks, Chandler, and Andersen (1971); Chandler, Weeks, and Andersen
(1983), in which the averaged effects of long-range interactions can be
treated as a small correction to the purely short-range system.

In contrast, the effects of long-range interactions are more subtle and play a
role in collective effects that are important for dielectric screening.
Moreover, long-range forces do not cancel at extended interfaces and instead
play a key role in interfacial physics. As a result, short-range systems
cannot describe interfacial structure and thermodynamics, as they do in the
bulk, and standard NN models fail to describe even the simplest liquid-vapor
interfaces Niblett, Galib, and Limmer (2021). The local molecular field (LMF)
theory of Weeks and coworkers provides a framework for capturing the average
effects of long-range interactions at interfaces through an effective external
field Rodgers and Weeks (2008a); Remsing, Liu, and Weeks (2016); Gao  _et al._
(2018); Gao, Remsing, and Weeks (2020); Remsing (2019). LMF theory also
provides physically intuitive insights into the roles of short- and long-range
forces at interfaces that can be leveraged to model nonuniform systems.

Here, we exploit the physical picture provided by liquid-state theory to
develop a general approach for learning long-range interactions in NN
potentials from ab initio calculations. We separate the atomic interactions
into appropriate short-range and long-range components and construct a
separate network to handle each part. Importantly, the short-range model is
isolated from the long-range interactions, such that each component is treated
independently. This separation also isolates the long-range response of the
system, enabling it to be learned. Short-range interactions can be learned
using established approaches. The short- and long-range components of the
potential are then connected through a rapidly-converging self-consistent
loop. The resulting self-consistent field neural network (SCFNN) model is able
to describe the effects of long-range interactions without the use of
effective charges or similar artificial quantities. We illustrate this point
through the development of a SCFNN model of liquid water. In addition to
capturing the local structure of liquid water, as evidenced by the radial
distribution function, the SCFNN model accurately describes long-range
structural correlations connected to dielectric screening, as well as the
response of liquid water to electrostatic fields.

##  I Results

###  I.1 Workflow of the Self-Consistent Field Neural Network Model

The SCFNN model consists of two modules that each target a specific response
of the system (Fig. 1). Module 1 predicts the electronic response via the
position of the maximally localized Wannier function centers (MLWFCs). Module
2 predicts the forces on the nuclear sites. In turn, each module consists of
two networks: one to describe the short-range interactions and one to describe
perturbations to the short-range system from long-range electric fields.
Together, these two modules (four networks) enable the model to predict the
total electrostatic properties of the system.

In the short-range system, the $v(r)=1/r$ portion of the Coulomb potential is
replaced by the short-range potential $v_{0}(r)=\mathop{\rm
erfc}(r/\sigma)/r$. Physically, $v_{0}(r)$ corresponds to screening the charge
distributions in the system through the addition of neutralizing Gaussian
charge distributions of opposite sign - the interactions are truncated by
Gaussians. Therefore, we refer to this system as the Gaussian-truncated (GT)
system Rodgers and Weeks (2008b); Remsing, Rodgers, and Weeks (2011); Rodgers
and Weeks (2009); Rodgers, Hu, and Weeks (2011); Rodgers and Weeks (2008a). By
making a physically meaningful choice for $\sigma$, the GT system can describe
the structure of bulk liquids with high accuracy but with a fraction of the
computational cost. Moreover, the GT system has served as a useful short-range
component system when modeling the effects of long-range fields Rodgers and
Weeks (2008b); Remsing, Liu, and Weeks (2016); Gao, Remsing, and Weeks (2020);
Baker III, Rodgers, and Weeks (2020); Cox (2020). Here, we choose $\sigma$ to
be 4.2 Å (8 Bohr), which is large enough for the GT system to accurately
describe hydrogen bonding and the local structure of liquid water Rodgers and
Weeks (2008b); Remsing, Rodgers, and Weeks (2011); Rodgers and Weeks (2009);
Rodgers, Hu, and Weeks (2011); Rodgers and Weeks (2008a).

Figure 1: Schematic of the self-consistent field neural network (SCFNN). The
SCFNN consists of two modules, each with two networks. One networks learns the
short-range interactions (S) and the other learns the effects of long-range
interactions (L). Module 1 learns the positions of maximally localized Wannier
function centers, $\mathbf{r}_{w}$, and Module 2 learns the forces,
$\mathbf{F}$, on the atomic nuclei, the positions of which are indicated by
$\mathbf{R}$.

The remaining part of the Coulomb interaction,
$v_{1}(r)=v(r)-v_{0}(r)=\mathop{\rm erf}(r/\sigma)/r$, is long ranged, but
varies slowly over the scale of $\sigma$. Because $v_{1}(r)$ is uniformly
slowly-varying, the effective field produced by $v_{1}(r)$ usually induces a
linear response in the GT system. The linear nature of the response makes the
effects of $v_{1}(r)$ able to captured by linear models. In the context of
neural networks, we demonstrate below that a linear network is sufficient to
learn the linear response induced by long-range interactions.

####  I.1.1 Module 1

The separation of interactions into short- and long-range components is
crucial to the SCFNN model. In particular, the two networks of each module are
used to handle this separation. Network 1S of Module 1 predicts the positions
of the MLWFCs in the short-range GT system, while Network 1L predicts the
perturbations to the MLWFC positions induced by the effective long-range
field. Networks 1S and 1L leverage Kohn’s theory on the nearsightedness of
electronic matter (NEM) Kohn (1996); Prodans and Kohn (2005). The NEM states
that Prodans and Kohn (2005) “local electronic properties, such as the density
$n(r)$, depend significantly on the effective external potential only at
nearby points.” Here, the effective external potential refers to the Kohn-Sham
effective potential, which includes the external potential and the self-
consistently determined long-range electric fields. Therefore, the NEM
suggests that the electronic density, and consequently the positions of the
MLWFCs, are ‘nearsighted’ with respect to the effective potential, but not to
the atomic coordinates, contrary to what has been assumed in previous work
that also uses local geometric information of atoms as input to neural
networks Krishnamoorthy  _et al._ (2021); Zhang  _et al._ (2020). An atom
located at $\mathbf{r}^{\prime}$ will affect the effective potential at
$\mathbf{r}$, even if $\mathbf{r}^{\prime}$ is far from $\mathbf{r}$, through
long-range electrostatic interactions. Consequently, current approaches to
generating NN models can only predict the position of MLWFCs for a purely
short-range system without long-range electrostatics, such as the GT system
Krishnamoorthy  _et al._ (2021); Zhang  _et al._ (2020). We exploit this fact
and use established NNs to predict the locations of the MLWFCs in the GT
system Krishnamoorthy  _et al._ (2021). To do so, we create a local reference
frame around each water molecule (Fig. 2) and use the coordinates of the
surrounding atoms as inputs to the neural network. The local reference system
preserves the rotational and translational symmetry of the system. The network
outputs the positions of the four MLWFCs around the central water, which are
then transformed to the laboratory frame of reference.

Network 1L predicts the response of the MLWFC positions to the effective field
$\mathbf{E}(\mathbf{r})$, defined as the sum of the external field,
$\mathbf{E}_{\rm ext}(\mathbf{r})$, and the long-range field from $v_{1}(r)$:

$\mathbf{E}(\mathbf{r})=\mathbf{E}_{\rm ext}(\mathbf{r})+\int
d\mathbf{r}\mkern 2.0mu\vphantom{r}^{\prime}\rho(\mathbf{r}\mkern
2.0mu\vphantom{r}^{\prime})\nabla v_{1}(|\mathbf{r}-\mathbf{r}\mkern
2.0mu\vphantom{r}^{\prime}|)\,,$ (1)

where $\rho(\mathbf{r}\mkern 2.0mu\vphantom{r}^{\prime})$ is the
_instantaneous_ charge density of the system, including nuclear and electronic
charges. Network 1L also introduces a local reference frame for each water
molecule. However, Network 1L takes as input both the local coordinates _and_
local effective electric fields. The NEM suggests that this local information
is sufficient to determine the perturbation in the MLWFC positions. Network 1L
outputs this change in the positions of the water molecule’s four MLWFCs, and
this perturbation is added to the MLWFC position determined in the GT system
to obtain the MLWFCs in the full system. We note that $\mathbf{E}(\mathbf{r})$
is a slowly-varying long-range field, such that the MLWFCs respond linearly to
this field. Therefore, Network 1L is constructed to be linear in
$\mathbf{E}(\mathbf{r})$. Table 1 demonstrates that the linear response
embodied by Network 1L predicts the perturbation of the MLWFCs with reasonable
accuracy.

Figure 2: Local frame around the central water. The $y$-axis is along the OH
bond. The $z$-axis is perpendicular to the plane of the molecule. The $x$-axis
is perpendicular to the $x$ and $z$ axis.

We now need to determine the effective field $\mathbf{E}(\mathbf{r})$. This
effective field depends on the electron density distribution, but evaluating
and including the full three-dimensional electron density for every
configuration in a training set requires a prohibitively large amount of
storage space. Instead, we approximate the electron density by the charge
density of the MLWFCs, assuming each MLWFC is a point charge of magnitude
$-2e_{0}$. This approximation is often used when computing molecular
multipoles, as needed to predict vibrational spectra, for example Zhang  _et
al._ (2020, 2021a). Here, it is important to note that the MLWFs of water are
highly localized, so that the center gives a reasonable representation of the
location of the MLWF. Moreover, the electron density is essentially smeared
over the scale of $\sigma$ through a convolution with $v_{1}(r)$, which makes
the resulting fields relatively insensitive to small-wavelength variations in
the charge density. As a result, the electron density can be accurately
approximated by the MLWFC charge density within our approach.

The effective field is a functional of the set of MLWFC positions,
$\mathbf{E}[\left\\{\mathbf{r}_{w}\right\\}]$, and the positions of the MLWFCs
themselves depend on the field, $\mathbf{r}_{w}[\mathbf{E}]$. Therefore, we
determine $\mathbf{E}$ and $\left\\{\mathbf{r}_{w}\right\\}$ through self-
consistent iteration. Our initial guess for $\mathbf{E}$ is obtained from the
positions of the MLWFCs in the GT system. We then iterate this self-consistent
loop until the MLWFC positions no longer change, within a tolerance of
$2.6\times 10^{-4}$ Å. In practice, we find that self-consistency is achieved
quickly.

####  I.1.2 Module 2

After Module 1 predicts the positions of the MLWFCs, Module 2 predicts the
forces on the atomic sites. As with the first module, Module 2 consists of two
networks: one that predicts the forces of the GT system and another that
predicts the forces produced by $\mathbf{E}(\mathbf{r})$. To predict the
forces in the GT system, we adopt the network used by Behler and coworkers
Morawietz  _et al._ (2016). This network, Network 2S, takes local geometric
information of the atoms as inputs and, consequently, cannot capture long-
range interactions. To describe long-range interactions, we introduce a second
network (Network 2L in Fig. 1). This additional network predicts the forces on
atomic sites due to the effective field $\mathbf{E}(\mathbf{r})$, which
properly accounts for long-range interactions in the system. In practice, we
again introduce a local reference frame for each water molecule and use local
atomic coordinates and local electric fields as inputs. In this case, we also
find that a network that is linear in $\mathbf{E}(\mathbf{r})$ accurately
predicts the resulting long-range forces, consistent with the linear response
of the system to a slowly-varying field.

In practice, separating the data obtained from standard DFT calculations into
the GT system and the long-range effective field is not straightforward. To
solve this problem, we apply homogeneous electric fields of varying strength
while keeping the atomic coordinates fixed. The fields only perturb the
positions of the MLWFCs and the forces on the atoms — these perturbations are
not related to the GT system. The changes induced by these electric fields are
directly obtained from DFT calculations and are used to train Networks 1L and
2L, which learn the response to long-range effective fields. The remaining
part of the DFT data is used to train Networks 1S and 2S, which learn the
response of the short-ranged GT system. See the Methods section for a more
detailed discussion of the networks and the training procedure.

We emphasize that our approach to partitioning the system into a short-range
GT piece and a long-range perturbation piece is different from other machine
learning approaches for handling long-range electrostatics. The standard
approaches usually partition the total energy into two parts, a short-ranged
energy and an Ewald energy that is used to evaluate the long-range
interactions. However, this partitioning results in a coupling between the
short- and long-range interactions. For example, the short-range part of the
energy in the 4G-HDNNP model depends on the effective charges that are
assigned to the atoms, but these effective charges depend on long-range
electrostatic interactions through the global charge equilibration process
used to determine their values Ko  _et al._ (2021). In contrast, the approach
we propose here isolates the short-range and long-range physics. The GT system
does not depend on long-range electrostatics even implicitly; it is completely
uncoupled from the long-range interactions. The effects of long-range
electrostatic interactions are completely isolated within the second network
of each module, Network 1L and Network 2L in Fig. 1. This separation of short-
and long-ranged effects is similar in spirit to the principles underlying LMF
theory Rodgers and Weeks (2008a); Remsing, Liu, and Weeks (2016); Gao,
Remsing, and Weeks (2020) and related theories of uniform liquids Weeks,
Chandler, and Andersen (1971); Chandler, Weeks, and Andersen (1983); Rodgers,
Hu, and Weeks (2011); Rodgers and Weeks (2009); Remsing, Rodgers, and Weeks
(2011).

|  $0.1$ V/Å |  $0.2$ V/Å
---|---|---
MLWFC | $F_{\rm O}$ | $F_{\rm H}$ | MLWFC | $F_{\rm O}$ | $F_{\rm H}$
MAE ($\times 100$) | $0.028$ | $1.4$ | $0.98$ | $0.056$ | $2.8$ | $2.0$
Table 1: When homogeneous external fields are applied, the location of the
maximally localized Wanner function center (MLWFC) and the forces on the
oxygen and hydrogen nuclei, $F_{\rm O}$ and $F_{\rm H}$, respectively, are
changed. This table shows the Mean Absolute Error (MAE) of Network 1L and 2L
in predicting the changes in the MLWFC positions (Å) and the forces (eV/Å)
along the $z$-direction when fields of strength 0.1 V/Å and 0.2 V/Å are
applied along $z$-direction. The predictions are made for the test sets and
the error is computed with respect to the DFT results. Figure 3: Comparison of
the radial distribution functions for (a) O-O, (b) O-H and (c) H-H
correlations in liquid water, as predicted by molecular dynamics simulations
of the self-consistent field neural network (SCFNN) and Behler-Parrinello (BP)
models.

###  I.2 Water’s Local Structure is Insensitive to Long-Range Interactions

We demonstrate the success of the SCFNN approach by modeling liquid water.
Water is the most important liquid on Earth. Yet, the importance of both
short- and long-range interactions makes it difficult to model. Short-range
interactions are responsible for water’s hydrogen bond network that is
essential to its structure and unusual but important thermodynamic properties
Ball (2008); Remsing, Rodgers, and Weeks (2011). Long-range interactions play
key roles in water’s dielectric response, interfacial structure, and can even
influence water-mediated interactions Prelesnik  _et al._ (2021). Because of
this broad importance, liquid water has served as a prototypical test system
for many machine learning-based models Morawietz  _et al._ (2016); Cheng  _et
al._ (2019); Zhang  _et al._ (2020, 2021b); Grisafi and Ceriotti (2019) Here,
we test our SCFNN model on a system of bulk liquid water by performing
molecular dynamics simulations of 1000 molecules in the canonical ensemble
under periodic boundary conditions.

One conventional test on the validity of a NN potential is to compare the
radial distribution function, $g(r)$, between atomic sites for the different
models. The $g(r)$ predicted by the SCFNN model is the same as that predicted
by the Behler-Parrinello model Morawietz  _et al._ (2016) for all three site-
site correlations in water (Fig. 3). This level of agreement may be expected,
based on previous work examining the structure of bulk water Rodgers and Weeks
(2008b); Rodgers, Hu, and Weeks (2011); Remsing, Rodgers, and Weeks (2011);
Gao  _et al._ (2018); Gao, Remsing, and Weeks (2020). The radial distribution
functions of water are determined mainly by short-range, nearest-neighbor
interactions, which arise from packing and hydrogen bonding; long-range
interactions have little effect on the main features of $g(r)$. Consequently,
purely short-range models, like the GT system, can quantitatively reproduce
the $g(r)$ of water Rodgers and Weeks (2008b); Rodgers, Hu, and Weeks (2011);
Remsing, Rodgers, and Weeks (2011); Gao  _et al._ (2018); Gao, Remsing, and
Weeks (2020). Similarly, the short-range Behler-Parrinello model accurately
describes the radial distribution functions, as does the SCFNN model, which
includes long-range interactions.

###  I.3 Long-Range Electrostatics and Dielectric Response

Though the short-range structure exemplified by the radial distribution
function is insensitive to long-range interactions, long-range correlations
are not. For example, the longitudinal component of the dipole density or
polarization correlation function evaluated in reciprocal space,
$\chi^{0}_{zz}(\mathbf{k})$, was recently shown to be sensitive to long-range
interactions Cox (2020). This correlation function is defined according to

$\chi^{0}_{zz}(\mathbf{k})=\frac{1}{V}\sum_{l,j}\frac{(\mathbf{k}\cdot\mathbf{p}_{l})\,(\mathbf{k}\cdot\mathbf{p}_{j})}{k^{2}}\,e^{-i\mathbf{k}\cdot\left(\mathbf{r}_{l}-\mathbf{r}_{j}\right)}\,,\,\,\mathrm{with}\,\,\mathbf{k}=k\hat{\mathbf{z}}\,.$
(2)

Here $\mathbf{p}_{j}$ is the dipole moment of water molecule $j$ and
$\mathbf{r}_{j}$ is the position of the oxygen atom of water molecule $j$.

The longitudinal polarization correlation function predicted by our SCFNN
model and the Behler-Parrinello agree everywhere except at small $k$,
indicating that long-range correlations are different in the two models (Fig.
4a). The long-wavelength behavior of the polarization correlation function is
related to the dielectric constant via

$\lim_{k\rightarrow 0}\chi^{0}_{zz}(\mathbf{k})=\epsilon_{0}k_{\rm
B}T\frac{\epsilon-1}{\epsilon}\,,$ (3)

where $\epsilon=78.4$ is the dielectric constant of water. The
$\chi^{0}_{zz}(\mathbf{k})$ predicted by our SCFNN model is consistent with
the expected behavior at small $k$. In contrast, short-range models, like the
GT system Cox (2020) and the Behler-Parrinello model, significantly deviate
from the expected asymptotic value. Consequently, these short-range models are
expected to have difficulties describing the dielectric screening that is
important in nonuniform systems Rodgers and Weeks (2008b); Remsing, Liu, and
Weeks (2016); Gao, Remsing, and Weeks (2020); Cox (2020); Niblett, Galib, and
Limmer (2021), for example.

To further examine the dielectric properties of the NN models, we can apply
homogeneous fields of varying strength to the system and examine its response.
To do so, we performed finite-field simulations at constant displacement
field, $\mathbf{D}$. These finite-$\mathbf{D}$ simulations Zhang and Sprik
(2016) can be naturally combined with our SCFNN model, unlike many other
neural network models. Following previous work Cox (2020), we use
$\mathbf{D}=D\hat{\mathbf{z}}$, vary the magnitude of the displacement field
from $D=0$ V/Å to $D=0.4$ V/Å, and examine the polarization, $P$, induced in
water. As shown in Fig. 4b, the polarization response of water to the external
field is accurately predicted by dielectric continuum theory, as expected,
further suggesting that the SCFNN model properly describes the dielectric
response of water. To the best of our knowledge, this is the first NN model
that can accurately describe the response of a system to external fields. We
emphasize that this response is achieved by learning the long-range response
via Networks 1L and 2L.

Figure 4: (a) The longitudinal polarization correlation function in reciprocal
space, $\chi^{0}_{zz}(\mathbf{k})$, shows differences between the self-
consistent field neural network (SCFNN) and Behler-Parrinello (BP) models at
low $k$. In particular, the SCFNN model plateaus as $k\rightarrow 0$ in a
manner consistent with the theoretical prediction (green line), while the BP
(short-range) model does not. (b) The polarization, $P$, induced by a
homogeneous displacement field of magnitude $D$ along the $z$-axis is
accurately predicted by the SCFNN model, evidenced by the agreement with
dielectric continuum theory (DCT) predictions.

##  II Discussion

In this work, we have presented a general strategy to construct NN potentials
that can properly account for the long-range response of molecular systems
that is responsible for dielectric screening and related phenomena. We
demonstrated that this model produces the correct long-range polarization
correlations in liquid water, as well as the correct response of liquid water
to external electrostatic fields. Both of these quantities are related to the
dielectric constant and require a proper description of long-range
interactions. In contrast, current derivations of NN potentials result in
short-range models that cannot capture these effects.

We anticipate that this approach will be of broad use to the molecular machine
learning and simulation community for modeling electrostatic and dielectric
properties of molecular systems. In contrast to short-range interactions that
must be properly learned to describe the different local environments
encountered at extended interfaces and at solute surfaces, the response of the
system to long-range, slowly-varying fields is quite general. Learning the
long-range response (through Networks 1L and 2L) is analogous to learning a
linear response, and we expect the resulting model to be relatively
transferable. As such, our resulting SCFNN model can make _predictions_ about
conditions on which it was not trained. For example, we trained the model for
_electric fields_ of magnitude 0, 0.1, and 0.2 V/Å, and then used this model
to successfully predict the response of the system to _displacement fields_
with magnitudes between 0 and 0.4 V/Å. This suggests that our approach can be
used to train NN models in more complex environments, like water at extended
interfaces, and then accurately predict the response of water to long-range
fields in those environments.

Finally, we note that our SCFNN approach is complementary to many established
methods for creating NN potentials. Learning the short-range, GT system
interactions can be accomplished with any method that uses local geometric
information. In this case, the precise form of Networks 1S and 2S can be
replaced with an alternative NN. Then, Networks 1L and 2L can be used as
defined here, within the general SCFNN workflow, resulting in a variant of the
desired NN potential that can describe the effects of long-range interactions.
Because of this, we expect our SCFNN approach to be transferable and readily
interfaced with current and future machine learning methods for modeling
short-range molecular interactions.

##  III Methods

###  III.1 Neural Network Potentials.

Our training and test set consists of 1571 configurations of 64 water
molecules. Homogeneous electric fields were applied to the system, as
described further in the next section. We used two thirds of the
configurations for training and one third to test the training of the network.

To train the networks we need to separate the DFT data into the GT system and
the long-range effective field. However, that separation is not
straightforward in practice. To achieve this, we use the differences in the
MLWFC locations and forces induced by different fields to fit Networks 1L and
2L. We now describe this procedure in detail for fitting Network 1L, and
Network 2L was fit following a similar approach.

To learn the effects of long-range interactions, we consider perturbations to
the positions of the MLWFCs induced by external electric fields of different
magnitudes. Consider applying two fields of strength $\left|\mathbf{E}\right|$
and $\left|\mathbf{E^{\prime}}\right|$. These fields will alter the MLWFC
positions by $\Delta\mathbf{r}_{w}[\mathbf{R},\mathbf{E}]$ and
$\Delta\mathbf{r}\mkern
2.0mu\vphantom{r}^{\prime}_{w}[\mathbf{R},\mathbf{E^{\prime}}]$, respectively.
However, both $\Delta\mathbf{r}_{w}$ and $\Delta\mathbf{r}\mkern
2.0mu\vphantom{r}^{\prime}_{w}$ are not directly obtainable from a single DFT
calculation. Instead, we can readily compute the difference in perturbations,
$\Delta\mathbf{r}_{w}-\Delta\mathbf{r}\mkern 2.0mu\vphantom{r}^{\prime}_{w}$,
directly from the DFT data, because

$\Delta\mathbf{r}_{w}-\Delta\mathbf{r}\mkern
2.0mu\vphantom{r}^{\prime}_{w}=\mathbf{r}_{w}-\mathbf{r}\mkern
2.0mu\vphantom{r}^{\prime}_{w}\,.$ (4)

Here, $\mathbf{r}_{w}$ and $\mathbf{r}\mkern 2.0mu\vphantom{r}^{\prime}_{w}$
are the locations of the MLWFCs in the full system in the presence of the
field $\mathbf{E}$ and $\mathbf{E^{\prime}}$, respectively, and these
positions can be readily computed in the simulations. These differences in the
MLWFC positions are used to fit Network 1L. In addition, we also exploit the
fact that $\Delta\mathbf{r}_{w}=0$ when $\mathbf{E}=0$. This allows us to fix
the zero point of Network 1L.

After fitting Networks 1L and 2L, we use them to predict the contribution of
the effective field to the MLWFC locations and forces. We then subtract that
part from the DFT data. What remains corresponds to the short-range GT system,
and this is used to train Networks 1S and 2S.

We now describe the detailed structure of the four networks used here.

Network 1S. In the local frame of water molecule $i$, we construct two types
of symmetry functions as inputs to Network 1S. The first type is the type 2
Behler-Parrinello symmetry function Behler (2011),

$G_{i}^{2}=\sum_{j\neq
i}\exp\left(-\eta(r_{ij}-r_{s})^{2}\right)f_{c}(r_{ij}).$ (5)

Here, $\eta$ and $r_{s}$ are parameters that adjust the width and center of
the Gaussian, and $f_{c}$ is a cutoff function whose value and slope goes to
zero at the radial cutoff $r_{c}$. We adopted the same cutoff function as
previous work Morawietz  _et al._ (2016), and the cutoff $r_{c}$ is set equal
to 12 Bohr.

The second type of symmetry function is similar to the type 4 Behler-
Parrinello symmetry functionBehler (2011). This symmetry function depends on
the angle between $\mathbf{r}_{ij}$ and the axis of the local frame,

$\mathbf{G}_{i}^{4}=\sum_{j\neq
i}2^{1-\zeta}(1+\lambda\frac{\mathbf{r}_{ij}}{r_{ij}})^{\zeta}\exp(-\eta
r_{ij}^{2})f_{c}(r_{ij}).$ (6)

Here, $\zeta$ and $\lambda$ are parameters that adjust the dependence of the
angular term.

We use 36 symmetry functions as input to Network 1S. Network 1S itself
consists of two hidden layers that contain 24 and 16 nodes. The output layer
consists of 12 nodes, corresponding to the three-dimensional coordinates of
the four MLWFCs of a central water molecule. Network 1S is a fully connected
feed-forward network, and we use $\tanh(x)$ as its activation function.

Network 1L. In the local frame of water molecule $i$, we construct one type of
symmetry function as input to Network 1L,

$\mathbf{EG}_{i}^{2}=\sum_{j}\mathbf{E}_{j}\exp(-\eta(r_{ij}-r_{s})^{2})f_{c}(r_{ij}).$
(7)

Here, $\mathbf{E}_{j}$ is the effective field exerted on atom $j$. We use 36
symmetry functions as inputs to Network 1L. Network 1L has no hidden layers.
The output layer consists of 12 nodes, corresponding to the three-dimensional
coordinates of the _perturbations_ of a water molecule’s four MLWFCs induced
by the external field.

Network 2S. Network 2S is exactly the same as the Behler-Parrinello Network
employed in previous work Morawietz  _et al._ (2016). In brief, the network
contains 2 hidden layers, each containing 25 nodes. Type 2 and 4 Behler-
Parrinello symmetry functions are used as inputs to the network. The network
for oxygen takes 30 symmetry functions as inputs, while the network for
hydrogen takes 27 symmetry functions as inputs. A hyperbolic tangent is used
as the activation function.

Network 2L. Network 2L uses the same type of symmetry function as Network 1L.
The network for the force on the oxygen and for the force on hydrogen are
trained independently. To predict the force on the oxygen, we center the local
frame on the oxygen atom. When the force on a hydrogen atom is the target, we
center the local frame on a hydrogen atom. We use 36 symmetry functions as
inputs to Network 2L. Network 2L has no hidden layers. The inputs map linearly
onto the forces on the atoms.

###  III.2 DFT Calculations.

The DFT calculations followed previous work Cheng  _et al._ (2019); Marsalek
and Markland (2017) and used published configurations of water as the training
set Cheng  _et al._ (2019). In short, all calculations were performed with
CP2K (version 7) Kühne  _et al._ (2020); VandeVondele  _et al._ (2005), using
the revPBE0 hybrid functional with 25% exact exchange Zhang and Yang (1998);
Adamo and Barone (1999); Goerigk and Grimme (2011), the D3 dispersion
correction of Grimme Grimme  _et al._ (2010), Goedecker-Tetter-Hutter
pseudopotentials Goedecker, Teter, and Hutter (1996), and TZV2P basis sets
VandeVondele and Hutter (2007), with a plane wave cutoff of 400 Ry. Maximally
localized Wannier function centers Marzari  _et al._ (2012) were evaluated
with CP2K, using the LOCALIZE option. The maximally localized Wannier function
spreads were minimized according to previous work Berghold  _et al._ (2000). A
homogeneous, external electric field was applied to the system using the Berry
phase approach, with the PERIODIC_EFIELD option in CP2K Souza, Íniguez, and
Vanderbilt (2002); Umari and Pasquarello (2002); Zhang, Hutter, and Sprik
(2016). Electric fields of magnitude 0, 0.1, and 0.2 V/Å were applied along
the $z$-direction of the simulation cell. Sample input files are given at
Zenodo (https://doi.org/10.5281/zenodo.5521328).

###  III.3 Molecular Dynamics Simulations

MD simulations are performed in the canonical (NVT) ensemble, with a constant
temperature of 300 K maintained using a Berendsen thermostat. The system
consisted of 1000 water molecules in a cubic box 31.2 Å in length. The
equations of motion were integrated with a timestep of 0.5 fs. Radial
distribution functions and longitudinal polarization correlations functions
were computed from 100 independent trajectories that were each 50 ps in
length. Finite-$\mathbf{D}$ simulations were performed under the same
simulation conditions, and each trajectory was 50 ps long at each magnitude of
$D$.

###  III.4 Data Availability.

Datasets used to train and test the NNP can be found at Zenodo
(https://doi.org/10.5281/zenodo.5521328).

###  III.5 Code Availability.

All DFT calculations were performed with CP2K version 7. In-house code was
used to construct the NN potentials and perform the MD simulations. These
codes are available at Github (https://github.com/andy90/SCFNN).

###### Acknowledgements.

We acknowledge the Office of Advanced Research Computing (OARC) at Rutgers,
The State University of New Jersey for providing access to the Amarel cluster
and associated research computing resources that have contributed to the
results reported here. We thank John Weeks for helpful comments on the
manuscript.

## References

  * Allen and Tildesley (2017) M. P. Allen and D. J. Tildesley,  _Computer simulation of liquids_ (Oxford university press, 2017).
  * Tuckerman  _et al._ (1996) M. E. Tuckerman, P. J. Ungar, T. Von Rosenvinge, and M. L. Klein, The Journal of Physical Chemistry 100, 12878 (1996).
  * Car and Parrinello (1985) R. Car and M. Parrinello, Physical review letters 55, 2471 (1985).
  * Chen  _et al._ (2018) M. Chen, L. Zheng, B. Santra, H.-Y. Ko, R. A. DiStasio Jr, M. L. Klein, R. Car, and X. Wu, Nature chemistry 10, 413 (2018).
  * Geissler  _et al._ (2001) P. L. Geissler, C. Dellago, D. Chandler, J. Hutter, and M. Parrinello, Science 291, 2121 (2001).
  * Lee  _et al._ (2008) T.-S. Lee, C. S. López, G. M. Giambaşu, M. Martick, W. G. Scott, and D. M. York, Journal of the American Chemical Society 130, 3053 (2008).
  * Walker, Crowley, and Case (2008) R. C. Walker, M. F. Crowley, and D. A. Case, Journal of computational chemistry 29, 1019 (2008).
  * Senn and Thiel (2006) H. M. Senn and W. Thiel, Atomistic approaches in modern biology , 173 (2006).
  * Dal Peraro  _et al._ (2007) M. Dal Peraro, P. Ruggerone, S. Raugei, F. L. Gervasio, and P. Carloni, Current opinion in structural biology 17, 149 (2007).
  * Sun, Ruzsinszky, and Perdew (2015) J. Sun, A. Ruzsinszky, and J. P. Perdew, Physical review letters 115, 036402 (2015).
  * Sun  _et al._ (2016) J. Sun, R. C. Remsing, Y. Zhang, Z. Sun, A. Ruzsinszky, H. Peng, Z. Yang, A. Paul, U. Waghmare, X. Wu, M. L. Klein, and J. P. Perdew, Nat. Chem. 8, 831 (2016).
  * Chen  _et al._ (2017) M. Chen, H.-H. Ko, R. C. Remsing, M. F. C. Andrade, B. Santra, Z. Sun, A. Selloni, R. Car, M. L. Klein, J. P. Perdew, and X. Wu, Proc. Natl. Acad. Sci. U.S.A. 114, 10846 (2017).
  * Furness  _et al._ (2020) J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew, and J. Sun, The journal of physical chemistry letters 11, 8208 (2020).
  * Zhang  _et al._ (2021a) C. Zhang, F. Tang, M. Chen, J. Xu, L. Zhang, D. Y. Qiu, J. P. Perdew, M. L. Klein, and X. Wu, The Journal of Physical Chemistry B (2021a).
  * Adamo and Barone (1999) C. Adamo and V. Barone, The Journal of chemical physics 110, 6158 (1999).
  * Gartner  _et al._ (2020) T. E. Gartner, L. Zhang, P. M. Piaggi, R. Car, A. Z. Panagiotopoulos, and P. G. Debenedetti, Proceedings of the National Academy of Sciences 117, 26040 (2020).
  * Zhang  _et al._ (2021b) L. Zhang, H. Wang, R. Car, and E. Weinan, Physical Review Letters 126, 236001 (2021b).
  * Deringer  _et al._ (2018) V. L. Deringer, N. Bernstein, A. P. Bartók, M. J. Cliffe, R. N. Kerber, L. E. Marbella, C. P. Grey, S. R. Elliott, and G. Csányi, The journal of physical chemistry letters 9, 2879 (2018).
  * Niu  _et al._ (2020) H. Niu, L. Bonati, P. M. Piaggi, and M. Parrinello, Nature communications 11, 1 (2020).
  * Deringer  _et al._ (2021) V. L. Deringer, N. Bernstein, G. Csányi, C. Ben Mahmoud, M. Ceriotti, M. Wilson, D. A. Drabold, and S. R. Elliott, Nature 589, 59 (2021).
  * Khaliullin  _et al._ (2011) R. Z. Khaliullin, H. Eshet, T. D. Kühne, J. Behler, and M. Parrinello, Nature materials 10, 693 (2011).
  * Bonati and Parrinello (2018) L. Bonati and M. Parrinello, Physical review letters 121, 265701 (2018).
  * Yue  _et al._ (2021) S. Yue, M. C. Muniz, M. F. Andrade, L. Zhang, R. Car, and A. Z. Panagiotopoulos, Journal of Chemical Physics 154 (2021), 10.1063/5.0031215.
  * Grisafi and Ceriotti (2019) A. Grisafi and M. Ceriotti, Journal of Chemical Physics 151 (2019), 10.1063/1.5128375, 1909.04512 .
  * Niblett, Galib, and Limmer (2021) S. P. Niblett, M. Galib, and D. T. Limmer, arXiv preprint arXiv:2107.06208 (2021).
  * Xie, Persson, and Small (2020) X. Xie, K. A. Persson, and D. W. Small, Journal of Chemical Theory and Computation 16, 4256 (2020).
  * Yao  _et al._ (2018) K. Yao, J. E. Herr, D. W. Toth, R. McKintyre, and J. Parkhill, Chemical Science 9, 2261 (2018).
  * Ko  _et al._ (2021) T. W. Ko, J. A. Finkler, S. Goedecker, and J. Behler, Nature Communications 12, 1 (2021).
  * Widom (1967) B. Widom, Science 157, 375 (1967).
  * Weeks, Chandler, and Andersen (1971) J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chem. Phys. 54, 5237 (1971).
  * Chandler, Weeks, and Andersen (1983) D. Chandler, J. D. Weeks, and H. C. Andersen, Science 220, 787 (1983).
  * Rodgers and Weeks (2008a) J. M. Rodgers and J. D. Weeks, Journal of Physics: Condensed Matter 20, 494206 (2008a).
  * Rodgers and Weeks (2009) J. M. Rodgers and J. D. Weeks, The Journal of chemical physics 131, 244108 (2009).
  * Remsing, Rodgers, and Weeks (2011) R. C. Remsing, J. M. Rodgers, and J. D. Weeks, J. Stat. Phys. 145, 313 (2011).
  * Rodgers and Weeks (2008b) J. M. Rodgers and J. D. Weeks, Proceedings of the National Academy of Sciences 105, 19136 (2008b), https://www.pnas.org/content/105/49/19136.full.pdf .
  * Rodgers, Hu, and Weeks (2011) J. M. Rodgers, Z. Hu, and J. D. Weeks, Molecular Physics 109, 1195 (2011).
  * Remsing, Liu, and Weeks (2016) R. C. Remsing, S. Liu, and J. D. Weeks, Proc. Natl. Acad. Sci. USA 113, 2819 (2016).
  * Gao  _et al._ (2018) A. Gao, L. Tan, M. I. Chaudhari, D. Asthagiri, L. R. Pratt, S. B. Rempe, and J. D. Weeks, The Journal of Physical Chemistry B 122, 6272 (2018).
  * Gao, Remsing, and Weeks (2020) A. Gao, R. C. Remsing, and J. D. Weeks, Proc. Natl. Acad. Sci. U.S.A. 117, 1293 (2020).
  * Remsing (2019) R. C. Remsing, Proceedings of the National Academy of Sciences 116, 23874 (2019), https://www.pnas.org/content/116/48/23874.full.pdf .
  * Baker III, Rodgers, and Weeks (2020) E. B. Baker III, J. M. Rodgers, and J. D. Weeks, The Journal of Physical Chemistry B 124, 5676 (2020).
  * Cox (2020) S. J. Cox, Proceedings of the National Academy of Sciences of the United States of America 117, 19746 (2020), 2008.01682 .
  * Kohn (1996) W. Kohn, Physical Review Letters 76, 3168 (1996).
  * Prodans and Kohn (2005) E. Prodans and W. Kohn, Proceedings of the National Academy of Sciences of the United States of America 102, 11635 (2005).
  * Krishnamoorthy  _et al._ (2021) A. Krishnamoorthy, K. I. Nomura, N. Baradwaj, K. Shimamura, P. Rajak, A. Mishra, S. Fukushima, F. Shimojo, R. Kalia, A. Nakano, and P. Vashishta, Physical Review Letters 126, 216403 (2021).
  * Zhang  _et al._ (2020) L. Zhang, M. Chen, X. Wu, H. Wang, E. Weinan, and R. Car, Physical Review B 102, 1 (2020), 1906.11434 .
  * Morawietz  _et al._ (2016) T. Morawietz, A. Singraber, C. Dellago, and J. Behler, Proceedings of the National Academy of Sciences of the United States of America 113, 8368 (2016).
  * Ball (2008) P. Ball,  _Chemical Reviews_ , Chemical Reviews 108, 74 (2008).
  * Prelesnik  _et al._ (2021) J. L. Prelesnik, R. G. Alberstein, S. Zhang, H. Pyles, D. Baker, J. Pfaendtner, J. J. De Yoreo, F. A. Tezcan, R. C. Remsing, and C. J. Mundy, Proceedings of the National Academy of Sciences 118 (2021).
  * Cheng  _et al._ (2019) B. Cheng, E. A. Engel, J. Behler, C. Dellago, and M. Ceriotti, Proceedings of the National Academy of Sciences 116, 1110 (2019).
  * Zhang and Sprik (2016) C. Zhang and M. Sprik, Physical Review B 93, 1 (2016), 1510.00643 .
  * Behler (2011) J. Behler, Journal of Chemical Physics 134 (2011), 10.1063/1.3553717.
  * Marsalek and Markland (2017) O. Marsalek and T. E. Markland, The journal of physical chemistry letters 8, 1545 (2017).
  * Kühne  _et al._ (2020) T. D. Kühne, M. Iannuzzi, M. Del Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, _et al._ , The Journal of Chemical Physics 152, 194103 (2020).
  * VandeVondele  _et al._ (2005) J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing, and J. Hutter, Comput. Phys. Commun. 167, 103 (2005).
  * Zhang and Yang (1998) Y. Zhang and W. Yang, Physical Review Letters 80, 890 (1998).
  * Goerigk and Grimme (2011) L. Goerigk and S. Grimme, Physical Chemistry Chemical Physics 13, 6670 (2011).
  * Grimme  _et al._ (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010).
  * Goedecker, Teter, and Hutter (1996) S. Goedecker, M. Teter, and J. Hutter, Phys. Rev. B 54, 1703 (1996).
  * VandeVondele and Hutter (2007) J. VandeVondele and J. Hutter, J. Chem. Phys. 127, 114105 (2007).
  * Marzari  _et al._ (2012) N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Rev. Mod. Phys. 84, 1419 (2012).
  * Berghold  _et al._ (2000) G. Berghold, C. J. Mundy, A. H. Romero, J. Hutter, and M. Parrinello, Phys. Rev. B 61, 10040 (2000).
  * Souza, Íniguez, and Vanderbilt (2002) I. Souza, J. Íniguez, and D. Vanderbilt, Physical review letters 89, 117602 (2002).
  * Umari and Pasquarello (2002) P. Umari and A. Pasquarello, Physical Review Letters 89, 157602 (2002).
  * Zhang, Hutter, and Sprik (2016) C. Zhang, J. Hutter, and M. Sprik, The journal of physical chemistry letters 7, 2696 (2016).


# Absence of anomalous dissipation for weak solutions of the Maxwell–Stefan
system

Luigi C. Berselli, Stefanos Georgiadis and Athanasios E. Tzavaras  Luigi C.
Berselli, Dipartimento di Matematica, Università di Pisa, Via F. Buonarroti
1/c, I56127, Pisa, Italy<EMAIL_ADDRESS>Stefanos Georgiadis,
Computer, Electrical and Mathematical Science and Engineering Division, King
Abdullah University of Science and Technology (KAUST), Thuwal 23955-6900,
Saudi Arabia and Institute for Analysis and Scientific Computing, Vienna
University of Technology, Wiedner Hauptstraße 8-10, 1040 Wien, Austria
<EMAIL_ADDRESS>Athanasios E. Tzavaras, Computer, Electrical
and Mathematical Science and Engineering Division, King Abdullah University of
Science and Technology (KAUST), Thuwal 23955-6900, Saudi Arabia
<EMAIL_ADDRESS>
###### Abstract.

In this paper we give a short and self-contained proof of the fact that weak
solutions to the Maxwell-Stefan system automatically satisfy an entropy
equality, establishing the absence of anomalous dissipation.

###### Key words and phrases:

Gas mixture, Maxwell–Stefan equations, isothermal model, nonequilibrium
thermodynamics, anomalous dissipation.

###### 2000 Mathematics Subject Classification:

35D30, 35Q35, 76N15, 76R50, 76T30.

##  1\. Introduction

Energy balance is an important tool in the understanding of fluid systems.
Already from the work of Kolmogorov in the 1940s, it became evident that the
evolution of energy in fluids plays a crucial role in the mechanism of
creation of turbulence. To this extent, establishing an energy equality
implies that the flow is not chaotic, as there is no turbulence (with some
caveat due to the variety of phenomena related with fluid flows). It was
formalized in the famous Onsager conjecture, stating that Hölder continuous
solutions of the incompessible Euler equations flows with exponent
$\alpha>1/3$ conserve the kinetic energy. The positive part of the conjecture
was proved in Constantin, E, and Titi [5], according to which weak solutions
to the Euler equations with some additional regularity conserve their energy.
The negative part of the conjecture was settled in an important work by Isett
[10] and Buckmaster et al. [3] building on several previous works (see
references therein). Following the work [5], several studies investigated
regularity criteria under which an energy conservation (or dissipation)
identity can be attained, including the incompressible [7] or compressible
[13] Navier-Stokes equations, the system of magnetohydrodynamics [4], the
Euler-Korteweg system [6] and general conservation laws [9]. See [1] for a
review of results on incompressible flows under various regularity assumptions
that induce conservation of energy.

The study of multicomponent systems is an active research field of
mathematical fluid mechanics for compressible flows. The Maxwell-Stefan
system, see (1.1) below, is one of the simplest models describing the
diffusive transport of the components of a gaseous mixture with zero mean
flow. It consists of a system of continuity equations, one for the mass
fraction of each component, together with an algebraic system describing
frictional balance among the components, and coupling the molar fluxes with
the mass fractions. The system is equipped with an entropy dissipation
structure obtained from thermodynamic considerations which plays an important
role in its analysis. It induces a gradient flow structure [11] which renders
it formally a parabolic system. Local existence of smooth solutions is shown
in [2] while global existence of weak solutions is shown in [12]. At present
it is not known if smooth solutions break down and in what norm.

The solutions constructed in [12] inherit from construction an energy
dissipation inequality. The objective of the present work is to show that weak
solutions to the Maxwell-Stefan system satisfy an energy dissipation identity,
making a first qualitative connection with similar results about standard
turbulence theory. The main result we prove is the absence of anomalous
dissipation for weak solutions to the Maxwell-Stefan system in the same class
where we have global existence, see Theorem 2. This property is obtained as a
consequence of the regularity of the weak solutions, and of the structure of
the equations, without additional assumptions.

###  1.1. Setting

We consider the system of Maxwell-Stefan equations which describes
multicomponent diffusive phenomena in a gas mixture. The system consists of
the partial differential equations

(1.1) $\partial_{t}c_{i}+\operatorname{div}J_{i}=0,\quad\textnormal{in
}\Omega\times(0,T),$

describing the evolution of the mass fractions $c_{i}$ for each component
$i=1,\dots,n$, where $J_{i}$ is the molar flux of the $i$-th component and the
evolution takes place in a bounded domain $\Omega\subset{\mathbb{R}}^{3}$ with
smooth boundary $\partial\Omega$.

Initially, the mass fractions are given by

(1.2) $c_{i}(x,0)=c_{i}^{0}(x)\quad\textnormal{ in }\Omega,$

where $c_{i}^{0}(x)\geq 0$ and $\sum_{i=1}^{n}c_{i}^{0}(x)=1$ and we assume
that no mass can enter or leave the domain, i.e., that

(1.3) $J_{i}\cdot\nu=0,\quad\textnormal{ on }\partial\Omega\times(0,T).$

where $\nu$ is the exterior normal unit vector to $\partial\Omega$.

The molar fluxes $J_{i}$ are solutions of the linear system

(1.4) $-\sum_{j\not=i,j=1}^{n}\frac{c_{j}J_{i}-c_{i}J_{j}}{D_{ij}}=\nabla
c_{i}\quad\textnormal{in }\Omega\times(0,T),$

where $D_{ij}$, for $i\not=j$, are positive and symmetric coefficients
describing binary interactions between the $i$-th and $j$-th components, and
the molar fluxes are subject to the constraint

(1.5) $\sum_{i=1}^{n}J_{i}=0.$

In the case of nonvanishing concentrations, the molar fluxes $J_{i}$ are
defined as the product of the mass fractions $c_{i}$ and the diffusional
velocities $u_{i}$, i.e. $J_{i}:=c_{i}u_{i}$. In [11] it was shown that system
(1.4)-(1.5) can be rewritten in terms of the quantities

$m_{i}:=\sqrt{c_{i}}u_{i},$

which make sense even for $c_{i}=0$, so that $J_{i}$ should be understood as
$J_{i}:=\sqrt{c_{i}}m_{i}$ and $m_{i}$ can be found by inverting the algebraic
linear system

(1.6) $-\sum_{j=1}^{n}A_{ij}m_{j}=2\nabla\sqrt{c_{i}},$

with

$A_{ij}=\begin{cases}-\frac{\sqrt{c_{i}c_{j}}}{D_{ij}},&i\not=j\\\
\sum_{k\not=i}\frac{c_{k}}{D_{ik}},&i=j\end{cases},$

subject to the constraint

(1.7) $\sum_{i=1}^{n}\sqrt{c}_{i}m_{i}=0.$

The main difficulty of the Maxwell-Stefan system is that system (1.4) (and
hence also when reformulated as (1.6)) is singular and thus yields no unique
solution. Nevertheless, under the additional constraint (1.5) (equivalently
(1.7)), which ensures conservation of total mass, the linear system can be
inverted on its image by means of the Bott-Duffin inverse (for more details we
refer to [11]).

To properly define the class of solution which we will handle, we define a
weak solution as follows:

######  Definition 1.

We call the vector field $\bm{c}=(c_{1},\dots,c_{n})$, such that $c_{i}\in
L^{2}(0,T;H^{1}(\Omega))$ and $\partial_{t}c_{i}\in
L^{2}(0,T;(H^{1}(\Omega))^{*})$, a weak solution of (1.1)-(1.4) if:

  * (i)

it satisfies (1.1) in the weak sense, i.e.

(1.8)
$\int_{0}^{T}\int_{\Omega}c_{i}\partial_{t}\varphi_{i}\,\textnormal{d}x\textnormal{d}t+\int_{0}^{T}\int_{\Omega}\sqrt{c_{i}}m_{i}\cdot\nabla\varphi_{i}\,\textnormal{d}x\textnormal{d}t=0,$

for all $i=1,\dots,n$ and for all $\varphi_{i}\in
C^{1}_{c}(\overline{\Omega}\times(0,T))$, where $m_{i}=\sqrt{c_{i}}u_{i}$ is
the solution of (1.6)-(1.7);

  * (ii)

$c_{i}\in L^{\infty}(0,T;L^{\infty}(\Omega))$, for all $i=1,\dots,n$,;

  * (iii)

$c_{i}\geq 0$ with $\sum_{i=1}^{n}c_{i}=1$, a.e. in $\Omega\times[0,T]$ for
all $i=1,\dots,n$,;

  * (iv)

$\sqrt{c_{i}}\in L^{2}(0,T;H^{1}(\Omega))$ and $c_{i}\in
C([0,T];L^{2}(\Omega))$ for all $i=1,\dots,n$, so that the initial conditions
are satisfied in the strong $L^{2}(\Omega)$ sense.

######  Remark 1.

Note that from the weak formulation (1.8) one can obtain the one considered in
[11], namely

$\int_{\Omega}c_{i}(T)\psi_{i}(T)\,\textnormal{d}x-\int_{\Omega}c_{i}^{0}\psi_{i}(0)\,\textnormal{d}x-\int_{0}^{T}\int_{\Omega}c_{i}\partial_{t}\psi_{i}\,\textnormal{d}x\textnormal{d}t-\int_{0}^{T}\int_{\Omega}\sqrt{c_{i}}m_{i}\cdot\nabla\psi_{i}\,\textnormal{d}x\textnormal{d}t=0,$

for all $\psi_{i}\in C^{1}([0,T];C^{1}(\overline{\Omega}))$, by choosing in
(1.8) $\varphi_{i}(x,t)=\eta_{\sigma}(t)\psi_{i}(x,t)$, where $\eta_{\sigma}$
is the cut-off in time defined by

(1.9) $\eta_{\sigma}(t)=\begin{cases}\frac{t-\sigma}{\sigma},&\sigma\leq t\leq
2\sigma\\\ 1,&2\sigma\leq t\leq T-2\sigma\\\
\frac{T-t-\sigma}{\sigma},&T-2\sigma\leq t\leq T-\sigma\\\
0,&t\leq\sigma,~{}t\geq T-\sigma,\end{cases}$

for $0<\sigma<\frac{T}{4}$ and $\psi_{i}\in
C^{1}([0,T];C^{1}(\overline{\Omega}))$ and then pass to the limit $\sigma\to
0$.

We note, here, that since weak solutions satisfy $c_{i}\in
L^{\infty}(0,T;L^{\infty}(\Omega))$ and $\nabla\sqrt{c_{i}}\in
L^{2}(0,T;L^{2}(\Omega))$, the inversion of system (1.6) yields $m_{i}\in
L^{2}(0,T;L^{2}(\Omega))$.

The existence of such a weak solution, satisfying all the above conditions was
proven in [12]. In addition to the properties of Definition 1, weak solutions
satisfy an entropy inequality. Let us define the quantity

(1.10)
$H(\bm{c}):=\sum_{i=1}^{n}\int_{\Omega}c_{i}(\ln{c_{i}}-1)\textnormal{d}x,$

which is the entropy functional associated to the system. Notice that
$H(\bm{c})\leq 0$, since $0\leq c_{i}\leq 1$ almost everywhere. Moreover, the
entropy is well-defined even for zero concentrations, since the limit of the
function $s\ln s$, as $s\to 0^{+}$, is equal to zero. Previous works,
especially [12, 11], have shown that weak solutions, as in Definition 1,
satisfy the following entropy inequality:

(1.11)
$H(\bm{c}(T))+\frac{1}{2}\sum_{i,j=1}^{n}\int_{0}^{T}\int_{\Omega}\frac{c_{i}c_{j}}{D_{ij}}|u_{i}-u_{j}|^{2}\,\textnormal{d}x\textnormal{d}t\leq
H(\bm{c}^{0}),$

where the dissipation term can be rewritten in terms of $\sqrt{c_{i}}$ and
$m_{i}$ as follows:

$\frac{1}{2}\sum_{i,j=1}^{n}\int_{0}^{T}\int_{\Omega}\frac{1}{D_{ij}}|\sqrt{c_{j}}m_{i}-\sqrt{c_{i}}m_{j}|^{2}\,\textnormal{d}x\textnormal{d}t,$

and thus is well-defined, since $c_{i}\in L^{\infty}(0,T;L^{\infty}(\Omega))$
and $m_{i}\in L^{2}(0,T;L^{2}(\Omega))$. In this work, we will show that, in
fact, weak solutions to the Maxwell-Stefan system as in Definition 1 satisfy
(1.11) as an identity. It should be noted that, according to the results from
[8], weak solutions are in fact renormalized solutions, meaning that they
further satisfy the identity

(1.12)
$\partial_{t}\beta(c_{i})+\operatorname{div}(\beta^{\prime}(c_{i})c_{i}u_{i})=c_{i}u_{i}\cdot\nabla\beta^{\prime}(c_{i}),$

in the weak sense, for all functions $\beta\in C^{2}([0,+\infty))$. Even
though a formal computation suggests that the entropy identity could be
obtained by choosing in (1.12) $\beta(s)=s\ln{s}-s$, such a $\beta$ is not
admissible, due to the fact that $s\ln{s}-s$ is not a $C^{1}$ function up to
$s=0$. To this degree, the aim of this paper is to make this formal
computation rigorous.

###  1.2. Main result

The fact that weak solutions satisfy instead of (1.11) the entropy equality
(1.13), is known in the literature as “absence of anomalous dissipation” and
it is generally a basic, but additional, assumptions used in the proofs of
weak-strong uniqueness results (see Theorem 1 in [11]). In particular, this
assumption is at the basis of the relative entropy method, which assumes that
the strong solution satisfies an entropy identity, allowing for the derivation
of the relative entropy inequality.

In this work we prove that all weak solutions satisfy the entropy equality,
implying that there is no anomalous dissipation. This can be interpreted in
analogy with incompressible fluids as the fact that the flow is not turbulent.
Since the entropy equality is “formally” obtained by testing the equations by
$\ln{c_{i}}$ and integrating by parts over the domain, the essential steps
will be the “rigorous” justification of these computations. As we will see,
this is not just a matter of regularity of the solution: the fact that the
solution may not be bounded away from zero (it is not enough simply to assume
that $c_{i}>0$ since estimates involving the logarithm are not uniform near
zero and one should ask the unrealistic $c_{i}\geq\overline{c}>0$) is a
peculiar issue and makes the justification non trivial.

The main result of this paper is, then, the following:

######  Theorem 2.

Let $\bm{c}$ be a weak solution to the Maxwell-Stefan system (1.1)-(1.4) and
let the initial data be such that $c_{i}^{0}(x)\geq 0$ and
$\sum_{i=1}^{n}c_{i}^{0}(x)=1$ a.e. in $\Omega$. Then $\bm{c}$ satisfies the
entropy identity

(1.13)
$H(\bm{c}(t))+\frac{1}{2}\sum_{i,j=1}^{n}\int_{0}^{t}\int_{\Omega}\frac{c_{i}c_{j}}{D_{ij}}|u_{i}-u_{j}|^{2}\,\textnormal{d}x\textnormal{d}t=H(\bm{c}^{0})\qquad\forall\,t\in(0,T].$

##  2\. Proof of Theorem 2

We first notice that even if we use the weak formulation in terms of
$\sqrt{c_{i}}$ and $m_{i}$ as in [11], the mass fractions $c_{i}$ of a weak
solution satisfy the regularity stated in Definition 1 and especially the
strong $L^{2}$-continuity.

######  Lemma 3.

Let $\bm{c}$ be a weak solution such that $c_{i}\in
L^{\infty}(0,T;L^{\infty}(\Omega))$ and $m_{i}\in L^{2}(0,T;L^{2}(\Omega))$.
Then, it follows that

(2.1) $\partial_{t}c_{i}\in L^{2}(0,T;H^{1}(\Omega)^{*}),$

and thus

(2.2) $c_{i}\in C([0,T];L^{2}(\Omega)).$

###### Proof.

Let $\langle\cdot,\cdot\rangle$ denote the duality pairing between the spaces
$L^{2}(0,T;H^{1}(\Omega)^{*})$ and $L^{2}(0,T;H^{1}(\Omega))$. Then, the
computation

$\begin{split}\int_{0}^{T}|\langle\partial_{t}c_{i},\varphi\rangle|\,\textnormal{d}t&=\int_{0}^{T}|\langle-\operatorname{div}(c_{i}u_{i}),\varphi\rangle|\,\textnormal{d}t\\\
&=\int_{0}^{T}|\langle c_{i}u_{i},\nabla\varphi\rangle|\,\textnormal{d}t\\\
&\leq\|\sqrt{c_{i}}\|_{L^{\infty}(0,T;L^{\infty}(\Omega))}\|m_{i}\|_{L^{2}(0,T;L^{2}(\Omega))}\|\nabla\varphi\|_{L^{2}(0,T;L^{2}(\Omega))},\end{split}$

shows that due to $m_{i}$ being in $L^{2}(0,T;L^{2}(\Omega))$,

$\partial_{t}c_{i}\in L^{2}(0,T;(H^{1}(\Omega)^{*}).$

Moreover, since

$c_{i}\in L^{2}(0,T;H^{1}(\Omega))\quad\textnormal{ and
}\quad\partial_{t}c_{i}\in L^{2}(0,T;H^{1}(\Omega)^{*}),$

and since

$H^{1}(\Omega)\hookrightarrow L^{2}(\Omega)\hookrightarrow H^{1}(\Omega)^{*}$

is a Gelfand evolution triple (by continuous and dense embedding), (2.2)
follows by Lions-Magenes interpolation theory. ∎

We can now give the proof of the main result of this paper.

###### Proof of Theorem 2.

In the weak formulation (1.8) choose
$\varphi_{i}(x,t)=(\eta_{\sigma}(t)\ln(c_{i}^{\epsilon}(x,t)+\delta))^{\epsilon}$,
where $\eta_{\sigma}$ is the cut-off in time defined in (1.9) for
$0<\sigma<\frac{T}{4}$ and $0<\delta<1$, a constant. Moreover, for
$0<\epsilon<\sigma$ we define $v^{\epsilon}$ to be the mollification of $v$ in
time, i.e.

$v^{\epsilon}(x,t):=\int_{0}^{T}\rho_{\epsilon}(t-\tau)v(x,\tau)\,\textnormal{d}\tau,$

where $\rho$ is a standard symmetric mollifier. Defining mollification this
way, we assume that the function $v$ is extended by zero outside $[0,T]$. This
gives the usual convergence $v^{\epsilon}\to v$ in the interior of $[0,T]$ and
eventually problems on the boundary; this is why we work with the cut-off in
time. We note that admissible test functions belong to the space
$C^{1}_{c}(\overline{\Omega}\times(0,T))$, but by a density argument it is
enough to consider them in $H^{1}(0,T;H^{1}(\Omega))$.

We first notice that

(2.3)
$\int_{0}^{T}\int_{\Omega}v\,\partial_{t}\phi^{\epsilon}\,\textnormal{d}x\textnormal{d}t=\int_{0}^{T}\int_{\Omega}v^{\epsilon}\,\partial_{t}\phi\,\textnormal{d}x\textnormal{d}t,$

for $\phi$ such that $\phi(x,0)=\phi(x,T)=0$.

Indeed, we have:

$\displaystyle\int_{0}^{T}\int_{\Omega}v(x,t)\partial_{t}\phi^{\epsilon}(x,t)\,\textnormal{d}x\textnormal{d}t$
$\displaystyle=\int_{0}^{T}\int_{\Omega}v(x,t)\left(\int_{0}^{T}\partial_{t}\rho_{\epsilon}(t-\tau)\phi(x,\tau)\,\textnormal{d}\tau\right)\,\textnormal{d}x\textnormal{d}t$
$\displaystyle=-\int_{0}^{T}\int_{\Omega}v(x,t)\left(\int_{0}^{T}\partial_{\tau}\rho_{\epsilon}(t-\tau)\phi(x,\tau)\,\textnormal{d}\tau\right)\,\textnormal{d}x\textnormal{d}t$
$\displaystyle=\int_{0}^{T}\int_{\Omega}v(x,t)\left(\int_{0}^{T}\rho_{\epsilon}(t-\tau)\partial_{\tau}\phi(x,\tau)\,\textnormal{d}\tau\right)\,\textnormal{d}x\textnormal{d}t$
$\displaystyle=\int_{0}^{T}\int_{\Omega}\partial_{\tau}\phi(x,\tau)\left(\int_{0}^{T}\rho_{\epsilon}(\tau-t)v(x,t)\,\textnormal{d}t\right)\,\textnormal{d}x\textnormal{d}\tau$
$\displaystyle=\int_{0}^{T}\int_{\Omega}v^{\epsilon}(x,t)\partial_{t}\phi(x,t)\,\textnormal{d}x\textnormal{d}t,$

where in the second equality we used
$\partial_{t}\rho_{\epsilon}(t-\tau)=-\partial_{\tau}\rho_{\epsilon}(t-\tau)$,
in the third one an integration by parts combined with $\phi(x,0)=\phi(x,T)=0$
and in the fourth equality we used Fubini-Tonelli and the fact that $\rho$ is
a symmetric mollifier.

Then, because of (2.3) and since
$\eta_{\sigma}(T)\ln(c_{i}^{\epsilon}(T)+\delta)=\eta_{\sigma}(0)\ln(c_{i}^{\epsilon}(0)+\delta)=0$,
the first term of (1.8) yields:

$\begin{split}\int_{0}^{T}\int_{\Omega}c_{i}\partial_{t}\varphi_{i}\,\textnormal{d}x\textnormal{d}t&=\int_{0}^{T}\int_{\Omega}c_{i}^{\epsilon}\partial_{t}\left[\eta_{\sigma}(t)\ln(c_{i}^{\epsilon}+\delta)\right]\,\textnormal{d}x\textnormal{d}t\\\
&=\int_{0}^{T}\int_{\Omega}c_{i}^{\epsilon}\eta_{\sigma}^{\prime}(t)\ln(c_{i}^{\epsilon}+\delta)\,\textnormal{d}x\textnormal{d}t+\int_{0}^{T}\int_{\Omega}c_{i}^{\epsilon}\eta_{\sigma}(t)\partial_{t}\ln(c_{i}^{\epsilon}+\delta)\,\textnormal{d}x\textnormal{d}t,\end{split}$

where

$\begin{split}\int_{0}^{T}\int_{\Omega}c_{i}^{\epsilon}\eta_{\sigma}^{\prime}(t)\ln(c_{i}^{\epsilon}+\delta)\,\textnormal{d}x\textnormal{d}t&=\frac{1}{\sigma}\int_{\sigma}^{2\sigma}\int_{\Omega}c_{i}^{\epsilon}\ln{(c_{i}^{\epsilon}+\delta)}\,\textnormal{d}x\textnormal{d}t-\frac{1}{\sigma}\int_{T-2\sigma}^{T-\sigma}\int_{\Omega}c_{i}^{\epsilon}\ln{(c_{i}^{\epsilon}+\delta)}\,\textnormal{d}x\textnormal{d}t,\end{split}$

and

$\begin{split}\int_{0}^{T}\int_{\Omega}c_{i}^{\epsilon}\eta_{\sigma}(t)\partial_{t}\ln(c_{i}^{\epsilon}+\delta)\,\textnormal{d}x\textnormal{d}t&=\int_{0}^{T}\int_{\Omega}c_{i}^{\epsilon}\eta_{\sigma}(t)\frac{\partial_{t}c_{i}^{\epsilon}}{c_{i}^{\epsilon}+\delta}\,\textnormal{d}x\textnormal{d}t\\\
&=\int_{0}^{T}\int_{\Omega}\eta_{\sigma}(t)\partial_{t}c_{i}^{\epsilon}\,\textnormal{d}x\textnormal{d}t-\delta\int_{0}^{T}\int_{\Omega}\eta_{\sigma}(t)\partial_{t}\ln(c_{i}^{\epsilon}+\delta)\,\textnormal{d}x\textnormal{d}t\\\
&=-\int_{0}^{T}\int_{\Omega}\eta_{\sigma}^{\prime}(t)c_{i}^{\epsilon}\,\textnormal{d}x\textnormal{d}t+\delta\int_{0}^{T}\int_{\Omega}\eta_{\sigma}^{\prime}(t)\ln(c_{i}^{\epsilon}+\delta)\,\textnormal{d}x\textnormal{d}t\\\
&=-\frac{1}{\sigma}\int_{\sigma}^{2\sigma}\int_{\Omega}(c_{i}^{\epsilon}-\delta\ln{(c_{i}^{\epsilon}+\delta)})\,\textnormal{d}x\textnormal{d}t\\\
&\phantom{xx}+\frac{1}{\sigma}\int_{T-2\sigma}^{T-\sigma}\int_{\Omega}(c_{i}^{\epsilon}-\delta\ln{(c_{i}^{\epsilon}+\delta)})\,\textnormal{d}x\textnormal{d}t.\end{split}$

The second term of (1.8) reads

$\begin{split}\int_{0}^{T}\int_{\Omega}c_{i}u_{i}\cdot\nabla\varphi_{i}\,\textnormal{d}x\textnormal{d}t&=\int_{0}^{T}\int_{\Omega}(c_{i}u_{i})^{\epsilon}\eta_{\sigma}(t)\cdot\nabla\ln{(c_{i}^{\epsilon}+\delta)}\,\textnormal{d}x\textnormal{d}t\\\
&=\int_{0}^{T}\int_{\Omega}(c_{i}u_{i})^{\epsilon}\eta_{\sigma}(t)\cdot\frac{\nabla
c_{i}^{\epsilon}}{c_{i}^{\epsilon}+\delta}\,\textnormal{d}x\textnormal{d}t,\end{split}$

and putting everything together we arrive at

(2.4)
$\begin{split}&\frac{1}{\sigma}\int_{\sigma}^{2\sigma}\int_{\Omega}c_{i}^{\epsilon}\ln{(c_{i}^{\epsilon}+\delta)}\,\textnormal{d}x\textnormal{d}t-\frac{1}{\sigma}\int_{T-2\sigma}^{T-\sigma}\int_{\Omega}c_{i}^{\epsilon}\ln{(c_{i}^{\epsilon}+\delta)}\,\textnormal{d}x\textnormal{d}t\\\
&\phantom{x}-\frac{1}{\sigma}\int_{\sigma}^{2\sigma}\int_{\Omega}(c_{i}^{\epsilon}-\delta\ln{(c_{i}^{\epsilon}+\delta)})\,\textnormal{d}x\textnormal{d}t+\frac{1}{\sigma}\int_{T-2\sigma}^{T-\sigma}\int_{\Omega}(c_{i}^{\epsilon}-\delta\ln{(c_{i}^{\epsilon}+\delta)})\,\textnormal{d}x\textnormal{d}t\\\
&\phantom{xx}+\int_{0}^{T}\int_{\Omega}(c_{i}u_{i})^{\epsilon}\eta_{\sigma}(t)\cdot\frac{\nabla
c_{i}^{\epsilon}}{c_{i}^{\epsilon}+\delta}\,\textnormal{d}x\textnormal{d}t=0.\end{split}$

The plan is to first let $\epsilon\to 0^{+}$ keeping $\sigma>0$ and $\delta>0$
fixed, then keeping $\delta>0$ fixed take the limit $\sigma\to 0^{+}$ and
finally let $\delta\to 0^{+}$.

Step 1: Keeping $\sigma>0$ and $\delta>0$ fixed, we pass to the limit
$\epsilon\to 0^{+}$. For the first two terms, due to mollification properties,
$0\leq c_{i}\leq 1$ and $\Omega$ being bounded, we have

$\begin{split}&c_{i}^{\epsilon}\ln{(c_{i}^{\epsilon}+\delta)}\to
c_{i}\ln{(c_{i}+\delta)}~{}\textnormal{almost everywhere in
}\Omega\times[0,T]\\\
&|c_{i}^{\epsilon}\ln{(c_{i}^{\epsilon}+\delta)}|\leq\max\\{\ln{(1+\delta)},|\ln{\delta}|\\}\in
L^{1}(0,T;L^{1}(\Omega)),\end{split}$

so that by the Dominated Convergence Theorem we can pass to the limit. The
next two terms are treated in a similar way.

For the last term, we have

$\begin{split}&\int_{0}^{T}\int_{\Omega}(c_{i}u_{i})^{\epsilon}\eta_{\sigma}\cdot\frac{\nabla
c_{i}^{\epsilon}}{c_{i}^{\epsilon}+\delta}\,\textnormal{d}x\textnormal{d}t-\int_{0}^{T}\int_{\Omega}c_{i}u_{i}\eta_{\sigma}\cdot\frac{\nabla
c_{i}}{c_{i}+\delta}\,\textnormal{d}x\textnormal{d}t\\\
&\phantom{xx}=\int_{0}^{T}\int_{\Omega}[(c_{i}u_{i})^{\epsilon}-c_{i}u_{i}]\eta_{\sigma}\cdot\frac{\nabla
c_{i}^{\epsilon}}{c_{i}^{\epsilon}+\delta}\,\textnormal{d}x\textnormal{d}t+\int_{0}^{T}\int_{\Omega}c_{i}u_{i}\eta_{\sigma}\cdot\frac{\nabla
c_{i}^{\epsilon}-\nabla
c_{i}}{c_{i}^{\epsilon}+\delta}\,\textnormal{d}x\textnormal{d}t\\\
&\phantom{xxxx}+\int_{0}^{T}\int_{\Omega}c_{i}u_{i}\eta_{\sigma}\cdot\nabla
c_{i}\left(\frac{1}{c_{i}^{\epsilon}+\delta}-\frac{1}{c_{i}+\delta}\right)\,\textnormal{d}x\textnormal{d}t,\end{split}$

and notice that $c_{i}u_{i}$ and $\nabla c_{i}$ are both in
$L^{2}(0,T;L^{2}(\Omega))$, due to $\sqrt{c_{i}}u_{i}$ and
$\nabla\sqrt{c_{i}}$ being in $L^{2}(0,T;L^{2}(\Omega))$ and $c_{i}\in
L^{\infty}(0,T;L^{\infty}(\Omega))$. This implies that
$(c_{i}u_{i})^{\epsilon}\to c_{i}u_{i}$ and $\nabla c_{i}^{\epsilon}=(\nabla
c_{i})^{\epsilon}\to\nabla c_{i}$ as $\epsilon\to 0$, both in
$L^{2}(0,T;L^{2}(\Omega))$. Then, due to the bounds $|\eta_{\sigma}(t)|\leq 1$
and $\frac{1}{c_{i}^{\epsilon}+\delta}\leq\frac{1}{\delta}$, we have

$\begin{split}\left|\int_{0}^{T}\int_{\Omega}[(c_{i}u_{i})^{\epsilon}-c_{i}u_{i}]\eta_{\sigma}\cdot\frac{\nabla
c_{i}^{\epsilon}}{c_{i}^{\epsilon}+\delta}\,\textnormal{d}x\textnormal{d}t\right|&\leq\frac{1}{\delta}\int_{0}^{T}\int_{\Omega}|(c_{i}u_{i})^{\epsilon}-c_{i}u_{i}||\nabla
c_{i}^{\epsilon}|\,\textnormal{d}x\textnormal{d}t\\\
&\leq\frac{1}{\delta}\|(c_{i}u_{i})^{\epsilon}-c_{i}u_{i}\|_{L^{2}(0,T;L^{2}(\Omega))}\|\nabla
c_{i}\|_{L^{2}(0,T;L^{2}(\Omega))}\\\ &\to 0,\textnormal{ as $\epsilon\to
0$,}\end{split}$

and similarly

$\begin{split}\left|\int_{0}^{T}\int_{\Omega}c_{i}u_{i}\eta_{\sigma}\cdot\frac{\nabla
c_{i}^{\epsilon}-\nabla
c_{i}}{c_{i}^{\epsilon}+\delta}\,\textnormal{d}x\textnormal{d}t\right|&\leq\frac{1}{\delta}\int_{0}^{T}\int_{\Omega}|c_{i}u_{i}||\nabla
c_{i}^{\epsilon}-\nabla c_{i}|\,\textnormal{d}x\textnormal{d}t\\\
&\leq\frac{1}{\delta}\|c_{i}u_{i}\|_{L^{2}(0,T;L^{2}(\Omega))}\|(\nabla
c_{i})^{\epsilon}-\nabla c_{i}\|_{L^{2}(0,T;L^{2}(\Omega))}\\\ &\to
0,\textnormal{ as $\epsilon\to 0$}.\end{split}$

For the third term, we apply the Dominated Convergence Theorem, since

$c_{i}u_{i}\eta_{\sigma}\cdot\nabla
c_{i}\left(\frac{1}{c_{i}^{\epsilon}+\delta}-\frac{1}{c_{i}+\delta}\right)\to
0\textnormal{ a.e. in $\Omega\times(0,T)$, as $\epsilon\to 0$,}$

and

$\begin{split}\left|c_{i}u_{i}\eta_{\sigma}\cdot\nabla
c_{i}\left(\frac{1}{c_{i}^{\epsilon}+\delta}-\frac{1}{c_{i}+\delta}\right)\right|&=|c_{i}u_{i}||\eta_{\sigma}||\nabla
c_{i}|\frac{|c_{i}-c_{i}^{\epsilon}|}{|(c_{i}+\delta)(c_{i}^{\epsilon}+\delta)|}\\\
&\leq\frac{2}{\delta^{2}}|c_{i}u_{i}||\nabla c_{i}|\in
L^{1}(0,T;L^{1}(\Omega)),\end{split}$

thus establishing

$\int_{0}^{T}\int_{\Omega}(c_{i}u_{i})^{\epsilon}\eta_{\sigma}\cdot\frac{\nabla
c_{i}^{\epsilon}}{c_{i}^{\epsilon}+\delta}\,\textnormal{d}x\textnormal{d}t\to\int_{0}^{T}\int_{\Omega}c_{i}u_{i}\eta_{\sigma}\cdot\frac{\nabla
c_{i}}{c_{i}+\delta}\,\textnormal{d}x\textnormal{d}t.$

Therefore, letting $\epsilon\to 0^{+}$ in (2.4), we obtain

(2.5)
$\begin{split}&\frac{1}{\sigma}\int_{\sigma}^{2\sigma}\int_{\Omega}(c_{i}+\delta)\ln{(c_{i}+\delta)}\,\textnormal{d}x\textnormal{d}t-\frac{1}{\sigma}\int_{T-2\sigma}^{T-\sigma}\int_{\Omega}(c_{i}+\delta)\ln{(c_{i}+\delta)}\,\textnormal{d}x\textnormal{d}t\\\
&\phantom{x}-\frac{1}{\sigma}\int_{\sigma}^{2\sigma}\int_{\Omega}c_{i}\,\textnormal{d}x\textnormal{d}t+\frac{1}{\sigma}\int_{T-2\sigma}^{T-\sigma}\int_{\Omega}c_{i}\,\textnormal{d}x\textnormal{d}t+\int_{0}^{T}\int_{\Omega}c_{i}u_{i}\eta_{\sigma}(t)\cdot\frac{\nabla
c_{i}}{c_{i}+\delta}\,\textnormal{d}x\textnormal{d}t=0.\end{split}$

Step 2: Now keeping $\delta>0$ fixed, let $\sigma\to 0^{+}$. We start with the
third term and see that

$\frac{1}{\sigma}\int_{\sigma}^{2\sigma}\int_{\Omega}c_{i}\,\textnormal{d}x\textnormal{d}t\to\int_{\Omega}c_{i}^{0}\,\textnormal{d}x.$

The reason is that

$\displaystyle\left|\frac{1}{\sigma}\int_{\sigma}^{2\sigma}\int_{\Omega}c_{i}\,\textnormal{d}x\textnormal{d}t-\int_{\Omega}c_{i}^{0}\,\textnormal{d}x\right|$
$\displaystyle=\left|\frac{1}{\sigma}\int_{\sigma}^{2\sigma}\int_{\Omega}(c_{i}-c_{i}^{0})\,\textnormal{d}x\textnormal{d}t\right|$
$\displaystyle\leq\frac{1}{\sigma}\int_{\sigma}^{2\sigma}\int_{\Omega}|c_{i}-c_{i}^{0}|\,\textnormal{d}x\textnormal{d}t$
$\displaystyle=\frac{1}{\sigma}\int_{\sigma}^{2\sigma}\|c_{i}-c_{i}^{0}\|_{L^{1}(\Omega)}\,\textnormal{d}t,$

which converges to zero as $\sigma\to 0$ since the integrand is continuous due
to $c_{i}\in C([0,T];L^{2}(\Omega))$. The fourth term is treated the exact
same way.

For the first term of (2.5), we use the Mean Value Theorem to write:

$\displaystyle\left|\frac{1}{\sigma}\int_{\sigma}^{2\sigma}\int_{\Omega}(c_{i}+\delta)\ln{(c_{i}+\delta)}\,\textnormal{d}x\textnormal{d}t-\int_{\Omega}(c_{i}^{0}+\delta)\ln{(c_{i}^{0}+\delta)}\,\textnormal{d}x\right|$
$\displaystyle\phantom{xxxxxx}\leq\frac{1}{\sigma}\int_{\sigma}^{2\sigma}\int_{\Omega}\left|(c_{i}+\delta)\ln{(c_{i}+\delta)}-(c_{i}^{0}+\delta)\ln{(c_{i}^{0}+\delta)}\right|\,\textnormal{d}x\textnormal{d}t$
$\displaystyle\phantom{xxxxxx}=\frac{1}{\sigma}\int_{\sigma}^{2\sigma}\int_{\Omega}|\ln{(\xi+\delta)}+1||c_{i}-c_{i}^{0}|\,\textnormal{d}x\textnormal{d}t$
$\displaystyle\phantom{xxxxxx}\leq
M(\delta)\frac{1}{\sigma}\int_{\sigma}^{2\sigma}\|c_{i}-c_{i}^{0}\|_{L^{1}(\Omega)}\,\textnormal{d}x\textnormal{d}t,$

where $\xi$ lies between $c_{i}$ and $c_{i}^{0}$ and thus $0<\xi<1$ and
$M(\delta)$ is a positive constant which only depends on $\delta$. Then,
taking the limit $\sigma\to 0^{+}$ we conclude that

$\frac{1}{\sigma}\int_{\sigma}^{2\sigma}\int_{\Omega}(c_{i}+\delta)\ln{(c_{i}+\delta)}\,\textnormal{d}x\textnormal{d}t\to\int_{\Omega}(c_{i}^{0}+\delta)\ln{(c_{i}^{0}+\delta)}\,\textnormal{d}x,$

and similarly for the second term of (2.5).

For the last term of (2.5), since $\eta_{\sigma}(t)\to 1$ almost everywhere in
$(0,T)$ we see that

$c_{i}u_{i}\eta_{\sigma}(t)\cdot\frac{\nabla c_{i}}{c_{i}+\delta}\to
c_{i}u_{i}\cdot\frac{\nabla c_{i}}{c_{i}+\delta}~{}\textnormal{almost
everywhere in }\Omega\times(0,T),$

and again due to the bounds $|\eta_{\sigma}(t)|\leq 1$ and
$\frac{1}{c_{i}^{\epsilon}+\delta}\leq\frac{1}{\delta}$ and the fact that
$c_{i}u_{i}$ and $\nabla c_{i}$ are in $L^{2}(0,T;L^{2}(\Omega))$ we obtain

$\left|c_{i}u_{i}\eta_{\sigma}(t)\cdot\frac{\nabla
c_{i}}{c_{i}+\delta}\right|\leq\frac{1}{\delta}|c_{i}u_{i}||\nabla c_{i}|\in
L^{1}(0,T;L^{1}(\Omega)).$

and Lebesgue’s Dominated Convergence Theorem can be applied.

Then, passing to the limit $\sigma\to 0^{+}$ in (2.5) gives

(2.6)
$\begin{split}&\int_{\Omega}(c_{i}^{0}+\delta)\ln{(c_{i}^{0}+\delta)}\,\textnormal{d}x-\int_{\Omega}(c_{i}(T)+\delta)\ln{(c_{i}(T)+\delta)}\,\textnormal{d}x\\\
&\phantom{xx}-\int_{\Omega}c_{i}^{0}\,\textnormal{d}x+\int_{\Omega}c_{i}(T)\,\textnormal{d}x+\int_{0}^{T}\int_{\Omega}c_{i}u_{i}\cdot\frac{\nabla
c_{i}}{c_{i}+\delta}\,\textnormal{d}x\textnormal{d}t=0.\end{split}$

Step 3: We remove $\delta$, needed to regularize the logarithm, by letting
$\delta\to 0^{+}$. The first term gives

$\int_{\Omega}(c_{i}^{0}+\delta)\ln(c_{i}^{0}+\delta)\,\textnormal{d}x\to\int_{\Omega}c_{i}^{0}\ln
c_{i}^{0}\,\textnormal{d}x,$

by the Dominated Convergence Theorem, since
$(c_{i}^{0}+\delta)\ln(c_{i}^{0}+\delta)\to c_{i}^{0}\ln c_{i}^{0}$ almost
everywhere as $\delta\to 0^{+}$ and
$|(c_{i}^{0}+\delta)\ln(c_{i}^{0}+\delta)|\leq C$ uniformly in $\delta$ for
some constant $C>0$ independent of $\delta$, since $0\leq c_{i}^{0}\leq 1$ and
$0<\delta<1$ and constants are integrable in $\Omega$. The second term is
treated similarly.

Finally, for the last term we obtain

$\int_{0}^{T}\int_{\Omega}c_{i}u_{i}\cdot\frac{\nabla
c_{i}}{c_{i}+\delta}\,\textnormal{d}x\textnormal{d}t\to
2\int_{0}^{T}\int_{\Omega}m_{i}\cdot\nabla\sqrt{c_{i}}\,\textnormal{d}x\textnormal{d}t,$

by the Dominated Convergence Theorem, since for $m_{i}=\sqrt{c_{i}}u_{i}$

$c_{i}u_{i}\cdot\frac{\nabla c_{i}}{c_{i}+\delta}\to
2m_{i}\cdot\nabla\sqrt{c_{i}}\quad\textnormal{ almost everywhere as $\delta\to
0^{+}$},$

and

$\begin{split}\left|c_{i}u_{i}\cdot\frac{\nabla
c_{i}}{c_{i}+\delta}\right|&=|m_{i}||2\nabla\sqrt{c_{i}}|\frac{c_{i}}{c_{i}+\delta}\\\
&\leq|m_{i}||2\nabla\sqrt{c_{i}}|\in L^{1}(0,T;L^{1}(\Omega)),\end{split}$

since by the regularity of weak solutions $m_{i}$ and $\nabla\sqrt{c_{i}}$ are
in $L^{2}(0,T;L^{2}(\Omega))$ and the bound is uniform in $\delta$.

Therefore, letting $\delta\to 0^{+}$ we obtain the entropy equality

(2.7) $\int_{\Omega}c_{i}(T)(\ln
c_{i}(T)-1)\,\textnormal{d}x-2\int_{0}^{T}\int_{\Omega}m_{i}\cdot\nabla\sqrt{c_{i}}\,\textnormal{d}x\textnormal{d}t=\int_{\Omega}c_{i}^{0}(\ln
c_{i}^{0}-1)\,\textnormal{d}x.$

We sum (2.7) over all $i\in\\{1,\dots,n\\}$ and use the definition of the
entropy (1.10), to arrive at

$H(\bm{c}(T))-2\sum_{i=1}^{n}\int_{0}^{T}\int_{\Omega}m_{i}\cdot\nabla\sqrt{c_{i}}\,\textnormal{d}x\textnormal{d}t=H(\bm{c}^{0}).$

Performing the following computation, which derives from (1.6) and the
symmetry of $D_{ij}$

$\begin{split}-\sum_{i=1}^{n}2m_{i}\cdot\nabla\sqrt{c_{i}}&=\sum_{i=1}^{n}\sum_{j=1}^{n}A_{ij}m_{j}m_{i}\\\
&=-\sum_{i=1}^{n}\sum_{j\neq i}\frac{1}{D_{ij}}c_{i}u_{i}\cdot
c_{j}u_{j}+\sum_{i=1}^{n}\sum_{j\neq i}\frac{1}{D_{ij}}c_{j}m_{i}^{2}\\\
&=\frac{1}{2}\sum_{i=1}^{n}\sum_{j\neq
i}\frac{c_{i}c_{j}}{D_{ij}}|u_{i}-u_{j}|^{2},\end{split}$

we obtain

(2.8)
$H(\bm{c}(T))+\frac{1}{2}\int_{0}^{T}\int_{\Omega}\sum_{i=1}^{n}\sum_{j\neq
i}\frac{c_{i}c_{j}}{D_{ij}}|u_{i}-u_{j}|^{2}\,\textnormal{d}x\textnormal{d}t=H(\bm{c}^{0}).$

The fact that $T>0$ is arbitrary, implies that (2.8) holds for all times
$t\in(0,T]$ and the proof is complete. ∎

## Acknowledgments

LCB acknowledges support by INdAM GNAMPA and also by MIUR, within project
PRIN20204NT8W4-004: Nonlinear evolution PDEs, fluid dynamics and transport
equations: theoretical foundations and applications. LCB also thanks King
Abdullah University of Science and Technology (KAUST) for the support and
hospitality during the preparation of the paper. The research of SG and AET
was partially supported by King Abdullah University of Science and Technology
(KAUST) baseline funds. SG also acknowledges partial support from the Austrian
Science Fund (FWF), grants P33010 and F65. This work has received funding from
the European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme, ERC Advanced Grant no. 101018153.

## References

  * [1] L.C. Berselli, _Three-dimensional Navier–Stokes equations for turbulence_ , Mathematics in Science and Engineering, Academic Press, London, [2021].
  * [2] D. Bothe. On the Maxwell–Stefan equations to multicomponent diffusion. In: J. Escher et al. (eds). Parabolic Problems. Progress in Nonlinear Differential Equations and their Applications, pp. 81–93. Springer, Basel, 2011.
  * [3] T. Buckmaster, C. De Lellis, L. Székelyhidi Jr., and V. Vicol, Onsager’s conjecture for admissible weak solutions. Comm. Pure Appl. Math. 72 (2019), 229–274.
  * [4] R. Caflisch, I. Klapper, G. Steele. Remarks on Singularities, Dimension and Energy Dissipation for Ideal Hydrodynamics and MHD. Commun. Math. Phys. 184 (1997), 443–455.
  * [5] P. Constantin, W. E, E. Titi. Onsager’s conjecture on the energy conservation for solutions of Euler’s equation. Commun. Math. Phys. 165 (1994), 207–209.
  * [6] T. Debiec, P. Gwiazda, A. Swierczewska-Gwiazda, A.E. Tzavaras. Conservation of energy for the Euler–Korteweg equations. Calc. Var. 57 (2018), 160.
  * [7] J. Duchon, R. Robert. Inertial energy dissipation for weak solutions of incompressible Euler and Navier–Stokes equations. Nonlinearity 13 (2000), 249-255.
  * [8] S. Georgiadis, H. Kim, A.E. Tzavaras. Uniqueness of renormalized solutions for the Maxwell–Stefan system. Under review, arXiv: 2311.10465 (2023).
  * [9] P. Gwiazda, M. Michalek, A. Swierczewska-Gwiazda. A Note on Weak Solutions of Conservation Laws and Energy/Entropy Conservation. Arch. Rational Mech. Anal. 229 (2018), 1223–1238.
  * [10] P. Isett. A proof of Onsager’s conjecture. Ann. of Math. (2) 188 (2018), 871–963.
  * [11] X. Huo, A. Jüngel, A.E. Tzavaras, Weak-Strong Uniqueness for Maxwell–Stefan Systems, SIAM J. Math. Anal. 54(3) (2022), 3215-–3252.
  * [12] A. Jüngel, I.V. Stelzer, Existence analysis of Maxwell-Stefan systems for multicomponent mixtures, SIAM J. Math. Anal. 45 (2013), 2421–2440.
  * [13] C. Yu. Energy Conservation for the Weak Solutions of the Compressible Navier–Stokes Equations. Arch. Rational Mech. Anal. 225 (2017), 1073–1087.