# Violation of the magnonic Wiedemann-Franz law in the strong nonlinear regime

Kouki Nakata  Advanced Science Research Center, Japan Atomic Energy Agency,
Tokai, Ibaraki 319-1195, Japan  Yuichi Ohnuma  Research Center for Advanced
Science and Technology (RCAST), The University of Tokyo, Meguro, Tokyo
153-8904, Japan  Se Kwon Kim  Department of Physics, Korea Advanced Institute
of Science and Technology, Daejeon 34141, Republic of Korea

###### Abstract

The celebrated Wiedemann-Franz (WF) law which governs the relation between
charge and heat transport traces back to the experimental discovery in 1853 by
Wiedemann and Franz. Despite the fundamental difference of the quantum-
statistical properties between fermions and bosons, the linear-in-$T$ behavior
of the WF law at low temperatures has recently been found to be the universal
property by the discovery of the WF law for magnon transport. However, the WF
law is for the linear response, and whether or not the universal law is valid
even in the nonlinear regime of Bose systems remains an open issue. Here we
provide a solution to this fundamental challenge. We show that the ratio of
the thermal to spin transport coefficient of magnons in topologically trivial
insulating magnets exhibits a different behavior from the linear response and
the universal law breaks down in the strong nonlinear regime. This finding is
within experimental reach with current device and measurement technologies.
Our discovery is the key ingredient in magnon-based spintronics, in the
evaluation of the figure of merit for thermomagnetic conversion elements of
spintronics devices.

##  I Introduction

The research on thermoelectric properties of materials started more than two
centuries ago and it has a long history. The celebrated Wiedemann-Franz (WF)
law which dictates the linear relation between charge and heat transport
traces back to the experimental discovery in 1853 by Wiedemann and Franz Franz
and Wiedemann (1853) that the ratio of the thermal $\kappa$ to electrical
conductivity $\sigma$ of several metals reduces to approximately the same
value for a fixed temperature. Lorenz established that this ratio is linear in
the absolute temperature $T$ and the proportionality constant takes a
material-independent value Lorenz (1881). Using quantum theory on solids,
Sommerfeld appropriately derived the universal constant, dubbed the Lorenz
number ${\cal{L}}_{\text{e}}$, which is independent of any material parameters
Sommerfeld (1928). Thus the WF law has been formulated Lifshitz and Pitaevskii
(1981); Ashcroft and Mermin (1976) at low temperature in the form of

$\displaystyle\frac{\kappa}{\sigma}$
$\displaystyle\stackrel{{\scriptstyle\rightarrow}}{{=}}{\cal{L}}_{\text{e}}T,$
$\displaystyle{\cal{L}}_{\text{e}}$
$\displaystyle:=\frac{{\pi}^{2}}{3}\Big{(}\frac{{k_{\text{B}}}}{e}\Big{)}^{2},$

where $e$ is the elementary charge and $k_{\text{B}}$ the Boltzmann constant.
This universal law characterizes the figure of merit for thermoelectric
conversion elements and has been playing a central role in electronics.

Toward efficient transmission of information that goes beyond what is offered
by conventional electronics, the last two decades have seen a rapid
development of spintronics, aiming at utilizing another degrees of freedom of
electrons, spins, by means of spin transport Chumak  _et al._ (2015); Nakata
_et al._ (2017a). For this holy grail, it is desirable to formulate the spin
analog of the WF law because the law is expected to be a promising building
block in spintronics, in the evaluation of the figure of merit for
thermomagnetic conversion elements of spintronics devices 111 We refer to the
spin analog of the figure of merit for thermoelectric conversion elements as
that for thermomagnetic ones $Z_{\text{s}}$, i.e.,
$Z_{\text{s}}:={\mathcal{S}}^{2}G/K$, where ${\mathcal{S}}$, $G$, and $K$ are
the spin Seebeck coefficient, the spin conductivity, and the thermal
conductivity of magnons, respectively.. Then we have established the magnon
analog of the WF law, namely, the WF law for magnon transport, in ferromagnets
and antiferromagnets Nakata  _et al._ (2015, 2017b, 2017c). Magnons are
bosonic magnetic excitations, i.e., the quantized spin-waves. Since magnons
carry the spin angular momentum, spin currents are generated in the absence of
charge currents when magnons propagate in insulating magnets. Thus the WF law
for magnon transport, dubbed the magnonic WF law, is defined as the ratio of
the thermal conductivity $K_{1}$ within the linear response regime to spin
conductivity $G$ of magnons. In the bulk of topologically trivial insulating
magnets at sufficiently low temperatures compared to the magnon energy gap,
the ratio reduces to Nakata  _et al._ (2017b)

$\displaystyle\frac{K_{1}}{G}$
$\displaystyle\stackrel{{\scriptstyle\rightarrow}}{{=}}{\cal{L}}_{1}T,$ (1a)
$\displaystyle{\cal{L}}_{1}$
$\displaystyle:=\frac{5}{2}\Big{(}\frac{{k_{\text{B}}}}{g\mu_{\text{B}}}\Big{)}^{2},$
(1b)

where $g$ is the $g$-factor of the constituent spins and $\mu_{\text{B}}$ the
Bohr magneton. The thermomagnetic ratio is linear in temperature. In analogy
to charge transport in metals, we refer to this behavior as the magnonic WF
law. The constant ${\cal{L}}_{1}$ analogous to the Lorenz number, i.e., the
magnetic Lorenz number of magnons, is independent of any material parameters
except the $g$-factor, which is material specific. The role of the charge $e$
is played by $g\mu_{\text{B}}$.

Magnons are bosonic excitations, while electrons are fermions. Still,
remarkably, the magnonic WF law exhibits the same linear-in-$T$ behavior at
low temperatures as the one for electronic transport despite the fact that the
quantum-statistical properties of bosons and fermions are fundamentally
different, particularly in the low temperature regime where quantum effects
dominate. In that sense, the linear-in-$T$ behavior of the WF law is found to
be the universal property.

As spintronics technologies develop, there has recently been a growing
interest in the properties of the nonlinear response Vera-Marun  _et al._
(2012); Sakimura  _et al._ (2014); Omar  _et al._ (2020); Zhang  _et al._
(2021); Tanikawa  _et al._ (2021); Nakai and Nagaosa (2019); Zeng  _et al._
(2019, 2020); Watanabe and Oshikawa (2020). However, the existing theory for
the magnonic WF law is intended to be applied only for the linear response
regime. Whether or not the universal law is valid even in the nonlinear regime
remains an open issue. In this paper, we provide a solution to this
fundamental challenge by using the Boltzmann equation. This is the main aim of
this paper.

This paper is organized as follows. In Sec. II we investigate longitudinal
thermal transport of magnons in the bulk of insulating magnets, and study the
validity and violation of the magnonic WF law in the nonlinear response
regime. Then, we remark on several issues in Sec. III. Finally, we give an
estimate for the experimental feasibility in Sec. IV and summarize in Sec. V.
Technical details are described in the Appendix.

##  II Nonlinear thermal transport

We consider longitudinal transport of magnons in the bulk of a topologically
trivial three-dimensional insulating magnet 222 We refer to the magnet where
Berry curvatures are zero as the topologically trivial magnet., subjected to a
temperature gradient, where the magnon of the energy dispersion relation
$\epsilon_{{{\mathbf{k}}}}=Dk^{2}+\Delta$ with the group velocity
${\mathbf{v}}_{{\mathbf{k}}}=\partial\epsilon_{{\mathbf{k}}}/(\partial\hbar{\mathbf{k}})$
carries the spin angular momentum $-1$ in units of the reduced Planck constant
$\hbar$: In which $k:=|{\mathbf{k}}|$ denotes the magnitude of the wavenumber
${\mathbf{k}}=(k_{x},k_{y},k_{z})$, $D$ represents the spin stiffness
constant, and $\Delta$ is the magnon energy gap, e.g., due to an external
magnetic field, a spin anisotropy, etc. Throughout this paper assuming that
the magnon energy gap takes a nonzero value $\Delta\neq 0$ and that the
nonequilibrium Bose distribution function of magnons $f_{{\mathbf{k}}}$ is
described by the Boltzmann equation within the quasiparticle approximation, we
study magnon transport at low temperatures using a relaxation time
approximation see . Note that if one assumed a magnon energy dispersion
including the $k$-linear term, nonreciprocal responses Tokura and Nagaosa
(2018), e.g., the Doppler shift of spin-waves Vlaminck and Bailleul (2008),
could arise in certain materials with broken inversion symmetry, which is
outside the scope of this paper: See Ref. Go  _et al._ (2022) for the
generation of the magnon nonreciprocity in the absence of a finite energy gap,
i.e., nonreciprocal transport of gapless spin-waves.

The applied temperature gradient, ${\mathbf{\nabla}}T=({\text{const.}})$,
induces a magnonic spin current along the longitudinal direction, which leads
to an accumulation of magnons at the boundaries and builds up a nonuniform
magnetization in the magnet. This magnetization gradient plays a role of an
effective magnetic field gradient ${\mathbf{\nabla}}B$, which works as the
gradient of a nonequilibrium spin chemical potential Basso  _et al._ (2016);
Cornelissen  _et al._ (2016a); Du  _et al._ (2017); Demokritov  _et al._
(2006), and drives magnon currents Haldane and Arovas (1995); Meier and Loss
(2003); Fujimoto (2009); Nakata  _et al._ (2017b). Therefore, magnon transport
subjected to the temperature gradient along the $x$ axis,
$\partial_{x}T=({\text{const.}})$, is characterized including the nonlinear
response by the longitudinal transport coefficient $L_{ij}$ for
$i\in{\mathbb{N}}$ and $j\in{\mathbb{N}}$ as

$\displaystyle\begin{pmatrix}j_{x}^{\text{s}}\\\
j_{x}^{\text{E}}\end{pmatrix}=\begin{pmatrix}L_{11}&L_{12}&L_{13}&L_{14}&L_{15}&L_{16}\\\
L_{21}&L_{22}&L_{23}&L_{24}&L_{25}&L_{26}\end{pmatrix}\begin{pmatrix}\partial_{x}B\\\
-\frac{\partial_{x}T}{T}\\\ (\partial_{x}B)^{2}\\\
-(\frac{\partial_{x}T}{T})(\partial_{x}B)\\\ (\frac{\partial_{x}T}{T})^{2}\\\
-(\frac{\partial_{x}T}{T})^{3}\end{pmatrix},$

where the spin current density
${\mathbf{j}}^{\text{s}}=(j_{x}^{\text{s}},j_{y}^{\text{s}},j_{z}^{\text{s}})$
and the energy current density
${\mathbf{j}}^{\text{E}}=(j_{x}^{\text{E}},j_{y}^{\text{E}},j_{z}^{\text{E}})$
are defined as Nakata  _et al._ (2017b)

$\displaystyle{\mathbf{j}}^{\text{s}}:=$
$\displaystyle-\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}g\mu_{\text{B}}{\mathbf{v}}_{{\mathbf{k}}}f_{{\mathbf{k}}},$
(3a) $\displaystyle{\mathbf{j}}^{\text{E}}:=$
$\displaystyle\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\epsilon_{{\mathbf{k}}}{\mathbf{v}}_{{\mathbf{k}}}f_{{\mathbf{k}}},$
(3b)

respectively. The transport coefficient $L_{11}$ is identified with the spin
conductivity of magnons Meier and Loss (2003) $G:=L_{11}$. In contrast to the
junction system Nakata  _et al._ (2015, 2018), the second-order response
vanishes 333 This result changes in general if one assumes a magnon energy
dispersion with the $k$-linear term. in the bulk of topologically trivial
insulating magnets due to the property of the odd function in $k_{x}$ see :

$\displaystyle L_{13}=L_{14}=L_{15}=L_{23}=L_{24}=L_{25}=0.$ (4)

Note that since the third-order response to the temperature gradient
$O\big{(}(\partial_{x}T)^{3}\big{)}$ can become dominant at sufficiently low
temperatures compared to the magnon gap see , we neglect the other third-order
terms such as $O\big{(}(\partial_{x}T)^{2}(\partial_{x}B)^{1}\big{)}$,
$O\big{(}(\partial_{x}T)^{1}(\partial_{x}B)^{2}\big{)}$, and
$O\big{(}(\partial_{x}B)^{3}\big{)}$.

In analogy to charge transport Ashcroft and Mermin (1976) and using Eq. (II),
we formulate thermal transport of magnons in the nonlinear response regime.
Under the applied temperature gradient, the magnonic spin current is generated
and this leads to an accumulation of magnons at the boundaries. Consequently,
the nonuniform magnetization is developed and this effective magnetic field
gradient $\partial_{x}B$ generates a counter-current of magnons. Then, the
system reaches a stationary state such that in- and out-flowing magnonic spin
currents balance each other $j_{x}^{\text{s}}=0$, which results in
$\partial_{x}B=\partial_{x}B^{\ast}$ with

$\displaystyle\partial_{x}B^{\ast}:=\frac{L_{12}}{L_{11}}\frac{\partial_{x}T}{T}+\frac{L_{16}}{L_{11}}\Big{(}\frac{\partial_{x}T}{T}\Big{)}^{3}.$
(5)

Thus the thermal conductivity is measured. This effective magnetic field
gradient in the new quasiequilibrium state $\partial_{x}B^{\ast}$ brings the
nonequilibrium spin chemical potential Basso  _et al._ (2016), being peculiar
to the system out of equilibrium, and contributes to the thermal conductivity
associated with the heat current density Nakata  _et al._ (2017b); Basso  _et
al._ (2016),
$j_{x}^{Q}=L_{21}\partial_{x}B^{\ast}-L_{22}{\partial_{x}T}/{T}-L_{26}({\partial_{x}T}/{T})^{3}$,
as

$\displaystyle
j_{x}^{Q}=-K_{1}\partial_{x}T-K_{2}(\partial_{x}T)^{2}-K_{3}(\partial_{x}T)^{3},$
(6)

where

$\displaystyle K_{1}$
$\displaystyle:=\frac{1}{T}\Big{(}L_{22}-\frac{L_{12}L_{21}}{L_{11}}\Big{)},$
(7a) $\displaystyle K_{2}$ $\displaystyle=0,$ (7b) $\displaystyle K_{3}$
$\displaystyle:=\frac{1}{T^{3}}\Big{(}L_{26}-\frac{L_{16}L_{21}}{L_{11}}\Big{)},$
(7c)

and $K_{1}$ represents the thermal conductivity in the linear response regime
and $K_{2(3)}$ is the thermal transport coefficient of the second-order
(third-order) nonlinear response. Note that in the stationary state under the
applied temperature gradient, the heat current density $j_{x}^{Q}$ is
different from the energy current density $j_{x}^{\text{E}}$ Nakata  _et al._
(2017b); Basso  _et al._ (2016); Nakata and Ohnuma (2021),

$\displaystyle j_{x}^{Q}\neq j_{x}^{\text{E}},$ (8)

in that

$\displaystyle j_{x}^{Q}=$ $\displaystyle
L_{21}\partial_{x}B^{\ast}-L_{22}\frac{\partial_{x}T}{T}-L_{26}\Big{(}\frac{\partial_{x}T}{T}\Big{)}^{3},$
(9a) $\displaystyle j_{x}^{\text{E}}=$ $\displaystyle\ \ \ \ \ 0\ \ \ \ \
-L_{22}\frac{\partial_{x}T}{T}-L_{26}\Big{(}\frac{\partial_{x}T}{T}\Big{)}^{3}.$
(9b)

We remark that if one wrongly omits the contribution of $\partial_{x}B^{\ast}$
associated with the counter-current and identifies $j_{x}^{\text{E}}$ as the
heat current density in theoretical calculation, the ratio of the thermal to
spin conductivity would not obey the magnonic WF law, breaking the linear-
in-$T$ behavior, even in the linear response regime Nakata  _et al._ (2017b).
Note that for thermal transport of electrons in metals, the contribution of
the counter-current is strongly suppressed by the sharp Fermi surface of
fermions at temperatures $k_{\text{B}}T$, which is much smaller than the Fermi
energy even at room temperature. This is the crucial difference in the thermal
conductivity between magnons and electrons, i.e., bosons and fermions,
respectively.

We evaluate the thermal transport coefficient of the third-order nonlinear
response $K_{3}$. The Onsager relation holds $L_{12}=L_{21}$ and at low
temperatures $k_{\text{B}}T\ll\Delta$, it reduces to Nakata  _et al._ (2017b)
$L_{12}/L_{11}=L_{21}/L_{11}\stackrel{{\scriptstyle\rightarrow}}{{=}}-{\Delta}/({g\mu_{\text{B}}})$.
Thus, the thermal transport coefficient of the third-order nonlinear response
at low temperature is recast into

$\displaystyle
K_{3}\stackrel{{\scriptstyle\rightarrow}}{{=}}\frac{1}{T^{3}}\Big{(}L_{26}+\frac{{\Delta}}{{g\mu_{\text{B}}}}L_{16}\Big{)}.$
(10)

A straightforward calculation using the Boltzmann equation provides the
transport coefficients of the nonlinear response at low temperatures as see

$\displaystyle L_{16}$
$\displaystyle\stackrel{{\scriptstyle\rightarrow}}{{=}}-g\mu_{\text{B}}F{\text{e}}^{-b}\Big{[}\frac{{\Delta}^{3}}{(\beta
D)^{7/2}}A_{3}+\frac{3D{\Delta}^{2}}{(\beta
D)^{9/2}}A_{4}+O(T^{11/2})\Big{]},$ (11a) $\displaystyle L_{26}$
$\displaystyle\stackrel{{\scriptstyle\rightarrow}}{{=}}F{\text{e}}^{-b}\Big{[}\frac{{\Delta}^{4}}{(\beta
D)^{7/2}}A_{3}+\frac{4D{\Delta}^{3}}{(\beta
D)^{9/2}}A_{4}+O(T^{11/2})\Big{]},$ (11b)

where $A_{n}:={\sqrt{\pi}}(2n)!/[2^{2n+1}(n!)]$ is the Gaussian integral for
$n\in{\mathbb{N}}$,
$F:=({2D}/{\hbar})^{4}[{{\tau}^{3}{\beta}^{3}}/({{10}{\pi}^{2}})]$, the
inverse temperature $\beta:=1/(k_{\text{B}}T)$, $b:=\beta\Delta$, and the
relaxation time $\tau$. Note that at low temperatures the relaxation time
takes a constant value of being temperature-independent: At sufficiently low
temperatures, the effect of magnon-magnon interactions and that of phonons are
negligibly small, and impurity scattering makes a major contribution to the
relaxation. Under the assumption that impurities are dilute and scattering is
elastic and spatially isotopic with the impurity potential localized in space,
the relaxation time at low temperatures reduces to Nakata and Ohnuma (2021)
$\tau\stackrel{{\scriptstyle\rightarrow}}{{=}}{\hbar}/({2\alpha\Delta})$,
where $\alpha$ is the Gilbert damping constant. Since the Gilbert damping
constant is little influenced by temperature Tserkovnyak  _et al._ (2005)
(i.e., the temperature dependence is negligibly small), the relaxation time at
low temperatures takes a constant value of being temperature-independent.

At low temperatures the thermal transport coefficient of the third-order
nonlinear response reduces to

$\displaystyle
K_{3}\stackrel{{\scriptstyle\rightarrow}}{{=}}\frac{1}{T^{3}}F{\text{e}}^{-b}\frac{D{\Delta}^{3}}{(\beta
D)^{9/2}}A_{4},$ (12)

and the spin conductivity to
$G\stackrel{{\scriptstyle\rightarrow}}{{=}}(g\mu_{\text{B}})^{2}{\text{e}}^{-b}\tau{(k_{\text{B}}T)^{3/2}}/({4{\pi}^{3/2}{\hbar}^{2}\sqrt{D}})$.
Thus, we find at low temperatures $k_{\text{B}}T\ll\Delta$ that the ratio of
the thermal transport coefficient of the third-order nonlinear response to the
spin conductivity is given as

$\displaystyle\frac{K_{3}}{G}$
$\displaystyle\stackrel{{\scriptstyle\rightarrow}}{{=}}{\cal{L}}_{3}\frac{1}{T^{3}},$
(13a) $\displaystyle{\cal{L}}_{3}$
$\displaystyle:=\frac{32{\tau}^{2}D{\Delta}^{3}A_{4}}{{5}{\sqrt{\pi}}{\hbar}^{2}(g\mu_{\text{B}})^{2}}.$
(13b)

The thermomagnetic ratio in the nonlinear regime is proportional to $1/T^{3}$
[Eq. (13a)]. This is in contrast to the one in the linear response regime
$K_{1}/G$, which exhibits the linear-in-$T$ behavior [Eq. (1a)]. The
proportionality constant ${\cal{L}}_{3}$, which is independent of temperature,
is less universal than the magnetic Lorenz number ${\cal{L}}_{1}$ [Eq. (1b)]
in that the constant depends on other material parameters as well as the
$g$-factor, such as the spin exchange interaction, the spin anisotropy, the
spin length, the lattice constant, etc. Note that each transport coefficient
does not diverge even at low temperature
$T\stackrel{{\scriptstyle\rightarrow}}{{=}}0$ as
$K_{3}\propto{\text{e}}^{-b}/T^{3/2}\stackrel{{\scriptstyle\rightarrow}}{{=}}0$
and $G\propto
T^{3/2}{\text{e}}^{-b}\stackrel{{\scriptstyle\rightarrow}}{{=}}0$.

Finally, we discuss the validity and violation of the magnonic WF law in the
nonlinear response regime. The magnonic WF law is originally for the linear
response [Eq. (1a)]. Instead of the thermal conductivity $K_{1}$, it is recast
in terms of the heat current density $j_{x}^{Q}=-K_{1}\partial_{x}T$ as
${j_{x}^{Q}}/{G}=-[{\cal{L}}_{1}(\partial_{x}T)]T$, which states that the
ratio of the heat current density to the magnonic spin conductivity $G$ is
linear in temperature for a fixed temperature gradient
$(\partial_{x}T)=(\text{const.})$. This is the magnonic WF law in terms of the
heat current density. Since the heat current density includes the nonlinear
response [Eq. (6)], it can be concluded that the magnonic WF law does hold
even in the nonlinear response regime if the ratio of the heat current density
to the magnonic spin conductivity $G$ exhibits the linear-in-$T$ behavior for
the fixed temperature gradient: This is the criterion for the magnonic WF law
in the nonlinear response regime.

From Eq. (6) the ratio of the heat current density including the nonlinear
response to the magnonic spin conductivity $G$ becomes

$\displaystyle\frac{j_{x}^{Q}}{G}=-\frac{K_{1}}{G}(\partial_{x}T)-\frac{K_{2}}{G}(\partial_{x}T)^{2}-\frac{K_{3}}{G}(\partial_{x}T)^{3}.$
(14)

Since the thermal transport coefficient of the second-order nonlinear response
vanishes $K_{2}=0$ [Eq. (7b)], the magnonic WF law holds even in the nonlinear
regime if the ratio of the thermal transport coefficient of the third-order
nonlinear response to the spin conductivity, $K_{3}/G$, exhibits the linear-
in-$T$ behavior. However, we find from Eq. (13a) that the ratio $K_{3}/G$ is
proportional to $1/T^{3}$ and does not exhibit the linear-in-$T$ behavior.
Thus, it is concluded that the magnonic WF law violates in the nonlinear
regime.

We remark that the magnonic WF law breaks down in the strong nonlinear regime
where the third-order nonlinear response to the temperature gradient
$O\big{(}(\partial_{x}T)^{3}\big{)}$ becomes relevant by the large temperature
gradient. It is not until the third-order response contributes that the law
violates. Since the thermal transport coefficient of the second-order
nonlinear response vanishes $K_{2}=0$, the universal law remains valid even in
the nonlinear regime where the temperature gradient is large but not enough
for the third-order response to become relevant: We refer to this region as
the weak nonlinear regime for convenience. In conclusion, in the bulk of
topologically trivial insulating magnets, the magnonic WF law remains valid
even in the weak nonlinear regime but breaks down in the strong nonlinear
regime.

##  III Discussion

In contrast to the bulk of topologically trivial materials studied in this
paper [Eq. (4)], the second-order nonlinear response does not vanish in
junction systems Nakata  _et al._ (2015, 2018), including quantum dot systems,
and it contributes to thermal transport. Therefore, in the quantum dot system
Sánchez and López (2013); López and Sánchez (2013) the WF law for electronic
transport violates in the weak nonlinear regime due to the second-order
response. Thus, we find that in the bulk of topologically trivial materials
the WF law is more robust against the nonlinear effect compared with in
junction systems, in that the law breaks down in the strong nonlinear regime
for the bulk of topologically trivial materials, while it violates in the weak
nonlinear regime for junction systems.

Note that throughout this paper, we focus on longitudinal thermal transport in
the bulk of topologically trivial insulating magnets where Berry curvatures
are zero, and find that the WF law breaks down in the strong nonlinear regime.
In the bulk of topological materials 444See Ref. Kondo and Akagi (2022) for
topological magnon systems., however, the WF law for electronic Hall transport
violates in the weak nonlinear regime due to the second-order response arising
from nonzero Berry curvatures Zeng  _et al._ (2020). Thus, it is concluded
that in the bulk of topologically trivial materials the WF law is more robust
against the nonlinear effect compared with in the bulk of topological
materials.

In this paper, we have studied magnon transport in the nonlinear response
regime under the assumption that the energy dispersion of magnons is gapped
and parabolic in terms of $k$. We remark that the second-order nonlinear
response [Eq. (4)] does not vanish in general if one assumes a magnon energy
dispersion including the $k$-linear term. In that case, the magnon
nonreciprocity Tokura and Nagaosa (2018), e.g., the spin-wave Doppler shift
Vlaminck and Bailleul (2008), could arise in certain materials with broken
inversion symmetry. See Ref. Go  _et al._ (2022) for nonreciprocal transport
in a gapless spin-wave system, i.e., the magnon nonreciprocity in the absence
of a finite energy gap.

In general, there exists a spin anisotropy in magnets, which causes the magnon
energy gap. Therefore, toward the development of various functions of
spintronics devices, it is of importance to establish the fundamental
principle of magnon transport even in the gapped systems. Hence, we studied
the gapped magnonic systems. In the gapped system, the magnonic WF law holds
only at low temperature $k_{\text{B}}T\ll\Delta$ and the linear-in-$T$
behavior violates at higher temperatures even within the linear response
Nakata  _et al._ (2015, 2017b). Therefore, in this paper focusing on the
gapped magnonic system at such low temperatures, we have studied the effect of
the nonlinear response on the magnonic WF law (i.e., the linear-in-$T$
behavior). We also remark that if the magnon gap is much smaller than the
thermal energy, then there would be a large number of low-energy magnons and
therefore frequent interactions between magnons. In this case, we would need
to consider the magnon-magnon interaction to capture the transport properties
well. In this work, we focused on the systems with the magnon gap larger than
the thermal energy in part so that the magnon-magnon interactions can be
neglected. Studying the effect of magnon-magnon interactions is beyond the
scope of our current work. For these reasons, in this paper, we have focused
on the gapped magnonic system at low temperatures. Still, it will be of
significance to develop our work into the gapless magnon mode that possesses
appropriate symmetry and include the magnon-magnon interaction. We leave the
advanced study for future work.

By manipulating the applied magnetic field, magnon and phonon thermal
conductivities can be distinguished experimentally Prasai  _et al._ (2017)
because the former does depend strongly on the field whereas the latter does
not. In fact, in Ref. Prasai  _et al._ (2017), it has been experimentally
shown that thermal conductivities of phonons and magnons in a magnetic
insulator can be separately characterized at low temperatures $T\leq O(1)$K.
For this reason, we focus only on magnon thermal conductivity in our work.

##  IV Estimates for experiments

In the bulk of topologically trivial insulating magnets, the heat current
density [Eq. (6)] consists of the linear response $K_{1}\partial_{x}T$ and the
third-order nonlinear response $K_{3}(\partial_{x}T)^{3}$. When the large
temperature gradient is applied enough that the third-order nonlinear response
begins to contribute in that the ratio of the nonlinear to linear response
amounts to $1/10$, $K_{3}(\partial_{x}T)^{3}/(K_{1}\partial_{x}T)=1/10$, we
identify it with the strong nonlinear regime. Thus, for observation of the
violation of the magnonic WF law in the bulk of insulating magnets, the
temperature gradient needs to reach
$|\partial_{x}T|=\sqrt{(1/10)(K_{1}/K_{3})}=:\mathcal{T}$, where
$K_{1}/K_{3}={25}{\sqrt{\pi}}{\hbar}^{2}k_{\text{B}}^{2}T^{4}/(64{\tau}^{2}D{\Delta}^{3}A_{4})$.
The criterion for observation of the violation is whether or not the
temperature gradient exceeds the value of $\mathcal{T}$.

For an estimate, we assume the following experimental parameter values Artman
_et al._ (1965); Kota and Imamura (2017) for ${\text{Cr}}_{2}{\text{O}}_{3}$:
$D=10$ meV(nm)2, ${\Delta}\sim 4$ meV, $\alpha=O(10^{-3})$, and $T=4$ K. This
results in ${\mathcal{T}}=O(10)$ K/mm. In addition, since the value of
${\mathcal{T}}$ is proportional to the Gilbert damping constant and it takes
Tserkovnyak  _et al._ (2005); Heinrich  _et al._ (2011) $\alpha=O(10^{-4})$ or
even $\alpha<O(10^{-4})$ for YIG depending on the shape, it is roughly
estimated as ${\mathcal{T}}=O(1)$ K/mm or even ${\mathcal{T}}<O(1)$ K/mm for
YIG, respectively. In both cases, those temperature gradients are
experimentally realizable Chudo .

Observation of long-distance transport of spin-wave spin currents Cornelissen
_et al._ (2015) and measurement of the nonequilibrium spin chemical potential
Du  _et al._ (2017), the magnonic spin conductivity Cornelissen  _et al._
(2016b, a), and the thermal conductivity Onose  _et al._ (2010) have been
reported. Refs. Tabuchi  _et al._ (2014, 2015); Kosen  _et al._ (2018) develop
magnonics technologies at low temperatures. Given these estimates, we expect
that observation of the magnonic WF law and the violation, while being
challenging, seem within experimental reach with current device and
measurement techniques. The key is to decrease the overall temperature while
the temperature gradient is maintained constant.

##  V Summary

Focusing on longitudinal transport of magnons at low temperatures in the bulk
of topologically trivial insulating magnets, we have studied the validity and
have found the violation of the magnonic Wiedemann-Franz law in the nonlinear
response regime. In terms of the heat current density, the magnonic Wiedemann-
Franz law is recast into the form that the ratio of the heat current density
to the magnonic spin conductivity is linear in temperature for a fixed
temperature gradient. Then we have shown that the universality of the
Wiedemann-Franz law, the linear-in-$T$ behavior, breaks down in the strong
nonlinear regime. In contrast to the linear response, the ratio of the thermal
transport coefficient of the third-order nonlinear response to the spin
conductivity is proportional to $1/T^{3}$, and the proportionality constant is
less universal in that it depends on other material parameters as well as the
$g$-factor. Thus, the universal law violates in the strong nonlinear regime
where the third-order nonlinear response becomes relevant by the large
temperature gradient. Since the second-order nonlinear response vanishes in
the bulk of topologically trivial insulating magnets, the magnonic Wiedemann-
Franz law remains valid even in the weak nonlinear regime but breaks down in
the strong nonlinear regime. Those findings are within experimental reach with
current device and measurement technologies. Toward efficient transmission of
information that goes beyond what is offered by conventional electronics, our
discovery is a promising building block in magnon-based spintronics, in the
evaluation of the figure of merit for thermomagnetic conversion elements of
spintronics devices.

###### Acknowledgements.

The author (K. N.) would like to thank D. Loss for the collaborative work on
the related study through fruitful discussions. We are grateful also to Y.
Araki for useful discussions and H. Chudo for helpful feedback on the
experimental feasibility. We acknowledge support by JSPS KAKENHI Grant Number
JP20K14420 (K. N.) and JP22K03519 (K. N.), by Leading Initiative for Excellent
Young Researchers, MEXT, Japan (K. N.), and by JST ERATO Grant No. JPMJER1601
(Y. O.). The author (S. K. K.) is supported by Brain Pool Plus Program through
the National Research Foundation of Korea funded by the Ministry of Science
and ICT (Grant No. NRF-2020H1D3A2A03099291) and by the National Research
Foundation of Korea funded by the Korea Government via the SRC Center for
Quantum Coherence in Condensed Matter (Grant No. NRF-2016R1A5A1008184). *

##  Appendix A Transport coefficients of magnons in the nonlinear response
regime

In this Appendix, we provide some details of the straightforward calculation
for the transport coefficient $L_{ij}$ of the nonlinear response in the bulk
of topologically trivial insulating magnets, and make a few remarks. Under the
relaxation time approximation, the Boltzmann equation of the quasiparticle
approximation describes the transport property of a steady state in terms of
time as Mahan (2000); Haug and Jauho (2007); Haldane and Arovas (1995); Meier
and Loss (2003); Fujimoto (2009); Basso  _et al._ (2016)

$\displaystyle\Big{(}{\mathbf{v}}_{{\mathbf{k}}}\cdot{\mathbf{\nabla}}T\frac{\partial}{\partial
T}-g\mu_{\text{B}}{\mathbf{v}}_{{\mathbf{k}}}\cdot{\mathbf{\nabla}}B\frac{\partial}{\partial\epsilon_{{{\mathbf{k}}}}}\Big{)}f_{{\mathbf{k}}}=-\frac{f_{{\mathbf{k}}}-f^{0}_{{\mathbf{k}}}}{{\tau}},$
(15)

where $\tau$ is the relaxation time, the Bose distribution function of magnons
out of equilibrium $f_{{\mathbf{k}}}$, and the one in equilibrium
$f^{0}_{{\mathbf{k}}}:=({\text{e}}^{\beta\epsilon_{{{\mathbf{k}}}}}-1)^{-1}$.
Defining the deviation from equilibrium
$g_{{\mathbf{k}}}:=f_{{\mathbf{k}}}-f^{0}_{{\mathbf{k}}}$, it is described as

$\displaystyle
g_{{\mathbf{k}}}=-\tau\Big{(}{\mathbf{v}}_{{\mathbf{k}}}\cdot{\mathbf{\nabla}}T\frac{\partial}{\partial
T}-g\mu_{\text{B}}{\mathbf{v}}_{{\mathbf{k}}}\cdot{\mathbf{\nabla}}B\frac{\partial}{\partial\epsilon_{{{\mathbf{k}}}}}\Big{)}(f^{0}_{{\mathbf{k}}}+g_{{\mathbf{k}}}).$
(16)

Using the method of successive substitution for $g_{{\mathbf{k}}}$, the
deviation up to $O({\tau}^{3})$ is given as

$\displaystyle g_{{\mathbf{k}}}=$ $\displaystyle
g_{1}+g_{2}+g_{3}+O({\tau}^{4}),$ (17a) $\displaystyle g_{1}:=$
$\displaystyle-\tau\Big{(}{\mathbf{v}}_{{\mathbf{k}}}\cdot{\mathbf{\nabla}}T\frac{\partial}{\partial
T}-g\mu_{\text{B}}{\mathbf{v}}_{{\mathbf{k}}}\cdot{\mathbf{\nabla}}B\frac{\partial}{\partial\epsilon_{{{\mathbf{k}}}}}\Big{)}f_{{\mathbf{k}}}^{0},$
(17b) $\displaystyle g_{2}:=$
$\displaystyle{\tau}^{2}\Big{(}{\mathbf{v}}_{{\mathbf{k}}}\cdot{\mathbf{\nabla}}T\frac{\partial}{\partial
T}-g\mu_{\text{B}}{\mathbf{v}}_{{\mathbf{k}}}\cdot{\mathbf{\nabla}}B\frac{\partial}{\partial\epsilon_{{{\mathbf{k}}}}}\Big{)}^{2}f_{{\mathbf{k}}}^{0},$
(17c) $\displaystyle g_{3}:=$
$\displaystyle-{\tau}^{3}\Big{(}{\mathbf{v}}_{{\mathbf{k}}}\cdot{\mathbf{\nabla}}T\frac{\partial}{\partial
T}-g\mu_{\text{B}}{\mathbf{v}}_{{\mathbf{k}}}\cdot{\mathbf{\nabla}}B\frac{\partial}{\partial\epsilon_{{{\mathbf{k}}}}}\Big{)}^{3}f_{{\mathbf{k}}}^{0},$
(17d)

where $g_{1}=O({\tau})$, $g_{2}=O({\tau}^{2})$, and $g_{3}=O({\tau}^{3})$.
From $\partial/(\partial
T)=-(\epsilon_{{{\mathbf{k}}}}/T)[\partial/(\partial\epsilon_{{{\mathbf{k}}}})]$,
the component $g_{3}$ is recast into

$\displaystyle
g_{3}={\tau}^{3}\Big{(}{\mathbf{v}}_{{\mathbf{k}}}\cdot{\mathbf{\nabla}}T\frac{\epsilon_{{{\mathbf{k}}}}}{T}+g\mu_{\text{B}}{\mathbf{v}}_{{\mathbf{k}}}\cdot{\mathbf{\nabla}}B\Big{)}^{3}\frac{{\partial}^{3}f_{{\mathbf{k}}}^{0}}{\partial{\epsilon_{{{\mathbf{k}}}}}^{3}}.$
(18)

Focusing on sufficiently low temperatures compared to the magnon gap and thus
assuming $\Delta/(k_{\text{B}}T)\gg
g\mu_{\text{B}}\partial_{x}B/(k_{\text{B}}\partial_{x}T)$, the third-order
response to the temperature gradient $O\big{(}({\mathbf{\nabla}}T)^{3}\big{)}$
becomes dominant. This condition can be met when the system has a sufficiently
large magnetic anisotropy and the resulting magnon gap, such as $\Delta\sim 4$
meV $=O(10)$ T for Artman  _et al._ (1965); Kota and Imamura (2017); Nakata
_et al._ (2017b) ${\text{Cr}}_{2}{\text{O}}_{3}$, $\Delta\sim 2$ meV $=O(10)$
T for Itoh  _et al._ (2013) ${\text{SrRu}}{\text{O}}_{3}$, and $\Delta\sim
0.4$ meV $=O(1)$ T for Lado and Fernández-Rossier (2017)
${\text{Cr}}{\text{I}}_{3}$. Then, the component reduces to

$\displaystyle
g_{3}\stackrel{{\scriptstyle\rightarrow}}{{=}}-{\tau}^{3}({\mathbf{v}}_{{\mathbf{k}}}\cdot{\mathbf{\nabla}}T)^{3}\frac{{\partial}^{3}f_{{\mathbf{k}}}^{0}}{\partial
T^{3}}+O\big{(}({\mathbf{\nabla}}T)^{2}({\mathbf{\nabla}}B)^{1}\big{)}+O\big{(}({\mathbf{\nabla}}T)^{1}({\mathbf{\nabla}}B)^{2}\big{)}+O\big{(}({\mathbf{\nabla}}B)^{3}\big{)},$
(19)

and we neglect the other third-order terms such as
$O\big{(}({\mathbf{\nabla}}T)^{2}({\mathbf{\nabla}}B)^{1}\big{)}$,
$O\big{(}({\mathbf{\nabla}}T)^{1}({\mathbf{\nabla}}B)^{2}\big{)}$, and
$O\big{(}({\mathbf{\nabla}}B)^{3}\big{)}$. Using the relation at low
temperatures $k_{\text{B}}T\ll\Delta$,

$\displaystyle\frac{{\partial}^{3}f_{{\mathbf{k}}}^{0}}{\partial
T^{3}}\stackrel{{\scriptstyle\rightarrow}}{{=}}{\text{e}}^{-\beta\epsilon_{{{\mathbf{k}}}}}\frac{(\beta\epsilon_{{{\mathbf{k}}}})^{3}}{T^{3}},$
(20)

we obtain the component as

$\displaystyle
g_{3}\stackrel{{\scriptstyle\rightarrow}}{{=}}-{\tau}^{3}({\mathbf{v}}_{{\mathbf{k}}}\cdot{\mathbf{\nabla}}T)^{3}{\text{e}}^{-\beta\epsilon_{{{\mathbf{k}}}}}\frac{(\beta\epsilon_{{{\mathbf{k}}}})^{3}}{T^{3}}.$
(21)

Assuming magnons which carry the spin angular momentum $-1$ in units of
$\hbar$, we define the spin current density ${\mathbf{j}}^{\text{s}}$ and the
energy current density ${\mathbf{j}}^{\text{E}}$ as

$\displaystyle{\mathbf{j}}^{\text{s}}=$
$\displaystyle-\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}g\mu_{\text{B}}{\mathbf{v}}_{{\mathbf{k}}}g_{{\mathbf{k}}},$
(22a) $\displaystyle{\mathbf{j}}^{\text{E}}=$
$\displaystyle\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\epsilon_{{\mathbf{k}}}{\mathbf{v}}_{{\mathbf{k}}}g_{{\mathbf{k}}}.$
(22b)

Substituting the function $g_{{\mathbf{k}}}=g_{1}+g_{2}+g_{3}$ into each
current density and performing the Gaussian integrals, we obtain the transport
coefficient of the nonlinear response at low temperatures as

$\displaystyle L_{13}$ $\displaystyle=L_{14}=L_{15}=L_{23}=L_{24}=L_{25}=0,$
(23a) $\displaystyle L_{16}$
$\displaystyle\stackrel{{\scriptstyle\rightarrow}}{{=}}-g\mu_{\text{B}}F{\text{e}}^{-b}\Big{[}\frac{D^{3}}{(\beta
D)^{13/2}}A_{6}+\frac{3D^{2}\Delta}{(\beta
D)^{11/2}}A_{5}+\frac{3D{\Delta}^{2}}{(\beta
D)^{9/2}}A_{4}+\frac{{\Delta}^{3}}{(\beta D)^{7/2}}A_{3}\Big{]},$ (23b)
$\displaystyle L_{26}$
$\displaystyle\stackrel{{\scriptstyle\rightarrow}}{{=}}F{\text{e}}^{-b}\Big{[}\frac{D^{4}}{(\beta
D)^{15/2}}A_{7}+\frac{4D^{3}\Delta}{(\beta
D)^{13/2}}A_{6}+\frac{6D^{2}{\Delta}^{2}}{(\beta
D)^{11/2}}A_{5}+\frac{4D{\Delta}^{3}}{(\beta
D)^{9/2}}A_{4}+\frac{{\Delta}^{4}}{(\beta D)^{7/2}}A_{3}\Big{]},$ (23c)
$\displaystyle\frac{L_{26}}{L_{16}}$
$\displaystyle\stackrel{{\scriptstyle\rightarrow}}{{=}}-\frac{\Delta}{g\mu_{\text{B}}},$
(23d)

where $A_{n}:={\sqrt{\pi}}(2n)!/[2^{2n+1}(n!)]$ for $n\in{\mathbb{N}}$,
$F:=({2D}/{\hbar})^{4}[{{\tau}^{3}{\beta}^{3}}/({{10}{\pi}^{2}})]$, and
$b:=\beta\Delta$. The second-order response vanishes due to the property of
the odd function, and this results in $K_{2}=0$. See Ref. Nakata  _et al._
(2017b) for the transport coefficient of the linear response, where the
Onsager relation holds $L_{12}=L_{21}$. At low temperatures, it reduces to
Nakata  _et al._ (2017b)
$L_{12}/L_{11}=L_{21}/L_{11}\stackrel{{\scriptstyle\rightarrow}}{{=}}-{\Delta}/({g\mu_{\text{B}}})$.
Therefore, the thermal transport coefficient of the third-order nonlinear
response $K_{3}$ at low temperature is given as

$\displaystyle
K_{3}\stackrel{{\scriptstyle\rightarrow}}{{=}}\frac{1}{T^{3}}\Big{(}L_{26}+\frac{{\Delta}}{{g\mu_{\text{B}}}}L_{16}\Big{)}.$
(24)

Finally, substituting $L_{26}$ and $L_{16}$ into Eq. (24) we obtain the
thermal transport coefficient of the third-order nonlinear response $K_{3}$ at
low temperatures in the main text.

Next, we remark on the relaxation time. The relaxation time is different from
the lifetime of magnons in general. Those are distinct quantities. However,
under the assumption that impurities are dilute and scattering is elastic and
spatially isotopic with the impurity potential localized in space (i.e., the
relaxation time depends solely on the magnitude of the wavenumber), the
relaxation time coincides with the lifetime. At low temperatures
$k_{\text{B}}T\ll\Delta$, the relaxation time reduces to Nakata and Ohnuma
(2021)

$\displaystyle\tau\stackrel{{\scriptstyle\rightarrow}}{{=}}\frac{1}{2\alpha}\frac{\hbar}{\Delta},$
(25)

where $\alpha$ is the Gilbert damping constant. Since the Gilbert damping
constant is little influenced by temperature Tserkovnyak  _et al._ (2005)
(i.e., the temperature dependence is negligibly small), it is concluded that
at low temperatures the relaxation time takes a constant value of being
temperature-independent. Note that at sufficiently low temperatures, the
effect of magnon-magnon interactions and that of phonons are negligibly small,
and impurity scattering makes a major contribution to the relaxation. See Ref.
Nakata and Ohnuma (2021) for details.

Lastly, we comment on the Boltzmann equation. Throughout this paper, we study
magnon transport using the conventional Boltzmann equation of the
quasiparticle approximation. From the viewpoint of quantum field theory, under
the assumption that the variation of the center-of-mass coordinates is slow
compared with that of the relative coordinates, the quantum kinetic equation
of the lowest order gradient approximation becomes the quantum Boltzmann
equation, which reduces to the conventional Boltzmann equation in the limit of
the quasiparticle approximation Mahan (2000); Haug and Jauho (2007).
Therefore, the criterion for magnon transport to be described by the
conventional Boltzmann equation is whether or not the quasiparticle
approximation is applicable to the system. From Ref. Nakata and Ohnuma (2021),
this results in the condition $\Delta\gg\hbar/(2\tau)$. Since
$\hbar/(2\tau)\stackrel{{\scriptstyle\rightarrow}}{{=}}\alpha\Delta$ at low
temperatures and $\alpha\leq O(10^{-3})$ for insulating magnets Tserkovnyak
_et al._ (2005); Heinrich  _et al._ (2011), the condition is satisfied. Thus,
it is concluded that magnon transport we study in this paper is described by
the conventional Boltzmann equation.

## References

  * Franz and Wiedemann (1853) R. Franz and G. Wiedemann, Annalen der Physik 165, 497 (1853).
  * Lorenz (1881) L. Lorenz, Ann. Phys. 249, 422 (1881).
  * Sommerfeld (1928) A. Sommerfeld, Z. Phys. 47, 1 (1928).
  * Lifshitz and Pitaevskii (1981) E. M. Lifshitz and L. P. Pitaevskii,  _Physical Kinetics: Landau and Lifshitz Course of Theoretical Physics Volume 10_ (Pergamon Press, 1981).
  * Ashcroft and Mermin (1976) N. W. Ashcroft and N. D. Mermin,  _Solid State Physics_ (Brooks Cole, Belmont, CA, 1976).
  * Chumak  _et al._ (2015) A. V. Chumak, V. I. Vasyuchka, A. A. Serga, and B. Hillebrands, Nat. Phys. 11, 453 (2015).
  * Nakata  _et al._ (2017a) K. Nakata, P. Simon, and D. Loss, J. Phys. D: Appl. Phys. 50, 114004 (2017a).
  * Note (1) We refer to the spin analog of the figure of merit for thermoelectric conversion elements as that for thermomagnetic ones $Z_{\text{s}}$, i.e., $Z_{\text{s}}:={\mathcal{S}}^{2}G/K$, where ${\mathcal{S}}$, $G$, and $K$ are the spin Seebeck coefficient, the spin conductivity, and the thermal conductivity of magnons, respectively.
  * Nakata  _et al._ (2015) K. Nakata, P. Simon, and D. Loss, Phys. Rev. B 92, 134425 (2015).
  * Nakata  _et al._ (2017b) K. Nakata, S. K. Kim, J. Klinovaja, and D. Loss, Phys. Rev. B 96, 224414 (2017b).
  * Nakata  _et al._ (2017c) K. Nakata, J. Klinovaja, and D. Loss, Phys. Rev. B 95, 125429 (2017c).
  * Vera-Marun  _et al._ (2012) I. J. Vera-Marun, V. Ranjan, and B. J. van Wees, Nat. Phys. 8, 313 (2012).
  * Sakimura  _et al._ (2014) H. Sakimura, T. Tashiro, and K. Ando, Nat. Commun. 5, 5730 (2014).
  * Omar  _et al._ (2020) S. Omar, M. Gurram, K. Watanabe, T. Taniguchi, M. Guimarães, and B. van Wees, Phys. Rev. Applied 14, 064053 (2020).
  * Zhang  _et al._ (2021) Z.-F. Zhang, Z.-G. Zhu, and G. Su, Phys. Rev. B 104, 115140 (2021).
  * Tanikawa  _et al._ (2021) Y. Tanikawa, K. Takasan, and H. Katsura, Phys. Rev. B 103, L201120 (2021).
  * Nakai and Nagaosa (2019) R. Nakai and N. Nagaosa, Phys. Rev. B 99, 115201 (2019).
  * Zeng  _et al._ (2019) C. Zeng, S. Nandy, A. Taraphder, and S. Tewari, Phys. Rev. B 100, 245102 (2019).
  * Zeng  _et al._ (2020) C. Zeng, S. Nandy, and S. Tewari, Phys. Rev. Research 2, 032066 (2020).
  * Watanabe and Oshikawa (2020) H. Watanabe and M. Oshikawa, Phys. Rev. B 102, 165137 (2020).
  * Note (2) We refer to the magnet where Berry curvatures are zero as the topologically trivial magnet.
  * (22) See the Appendix for details.
  * Tokura and Nagaosa (2018) Y. Tokura and N. Nagaosa, Nat. Commun. 9, 3740 (2018).
  * Vlaminck and Bailleul (2008) V. Vlaminck and M. Bailleul, Science 322, 410 (2008).
  * Go  _et al._ (2022) G. Go, S. Lee, and S. K. Kim, Phys. Rev. B 105, 134401 (2022).
  * Basso  _et al._ (2016) V. Basso, E. Ferraro, and M. Piazzi, Phys. Rev. B 94, 144422 (2016).
  * Cornelissen  _et al._ (2016a) L. J. Cornelissen, K. J. H. Peters, G. E. W. Bauer, R. A. Duine, and B. J. van Wees, Phys. Rev. B 94, 014412 (2016a).
  * Du  _et al._ (2017) C. Du, T. V. der Sar, T. X. Zhou, P. Upadhyaya, F. Casola, H. Zhang, M. C. Onbasli, C. A. Ross, R. L. Walsworth, Y. Tserkovnyak, and A. Yacoby, Science 357, 195 (2017).
  * Demokritov  _et al._ (2006) S. O. Demokritov, V. E. Demidov, O. Dzyapko, G. A. Melkov, A. A. Serga, B. Hillebrands, and A. N. Slavin, Nature (London) 443, 430 (2006).
  * Haldane and Arovas (1995) F. D. M. Haldane and D. P. Arovas, Phys. Rev. B 52, 4223 (1995).
  * Meier and Loss (2003) F. Meier and D. Loss, Phys. Rev. Lett. 90, 167204 (2003).
  * Fujimoto (2009) S. Fujimoto, Phys. Rev. Lett. 103, 047203 (2009).
  * Nakata  _et al._ (2018) K. Nakata, Y. Ohnuma, and M. Matsuo, Phys. Rev. B 98, 094430 (2018).
  * Note (3) This result changes in general if one assumes a magnon energy dispersion with the $k$-linear term.
  * Nakata and Ohnuma (2021) K. Nakata and Y. Ohnuma, Phys. Rev. B 104, 064408 (2021).
  * Tserkovnyak  _et al._ (2005) Y. Tserkovnyak, A. Brataas, G. E. W. Bauer, and B. I. Halperin, Rev. Mod. Phys. 77, 1375 (2005).
  * Sánchez and López (2013) D. Sánchez and R. López, Phys. Rev. Lett. 110, 026804 (2013).
  * López and Sánchez (2013) R. López and D. Sánchez, Phys. Rev. B 88, 045129 (2013).
  * Note (4) See Ref. Kondo and Akagi (2022) for topological magnon systems.
  * Prasai  _et al._ (2017) N. Prasai, B. A. Trump, G. G. Marcus, A. Akopyan, S. X. Huang, T. M. McQueen, and J. L. Cohn, Phys. Rev. B 95, 224407 (2017).
  * Artman  _et al._ (1965) J. O. Artman, J. C. Murphy, and S. Foner, Phys. Rev. 138, A912 (1965).
  * Kota and Imamura (2017) Y. Kota and H. Imamura, Appl. Phys. Express 10, 013002 (2017).
  * Heinrich  _et al._ (2011) B. Heinrich, C. Burrowes, E. Montoya, B. Kardasz, E. Girt, Y.-Y. Song, Y. Sun, and M. Wu, Phys. Rev. Lett. 107, 066604 (2011).
  * (44) H. Chudo, Private communication.
  * Cornelissen  _et al._ (2015) L. J. Cornelissen, J. Liu, R. A. Duine, J. B. Youssef, and B. J. van Wees, Nat. Phys. 11, 1022 (2015).
  * Cornelissen  _et al._ (2016b) L. J. Cornelissen, J. Shan, and B. J. van Wees, Phys. Rev. B 94, 180402 (2016b).
  * Onose  _et al._ (2010) Y. Onose, T. Ideue, H. Katsura, Y. Shiomi, N. Nagaosa, and Y. Tokura, Science 329, 297 (2010).
  * Tabuchi  _et al._ (2014) Y. Tabuchi, S. Ichino, T. Ishikawa, R. Yamazaki, K. Usami, and Y. Nakamura, Phys. Rev. Lett. 113, 083603 (2014).
  * Tabuchi  _et al._ (2015) Y. Tabuchi, S. Ichino, A. Noguchi, T. Ishikawa, R. Yamazaki, K. Usami, and Y. Nakamura, Science 349, 405 (2015).
  * Kosen  _et al._ (2018) S. Kosen, R. G. E. Morris, A. F. van Loo, and A. D. Karenowska, Appl. Phys. Lett. 112, 012402 (2018).
  * Mahan (2000) G. D. Mahan,  _Many-Particle Physics_ (Kluwer Academic, Plenum Publishers, New York, 2000).
  * Haug and Jauho (2007) H. Haug and A. Jauho,  _Quantum Kinetics in Transport and Optics of Semiconductors_ (Springer New York, 2007).
  * Itoh  _et al._ (2013) S. Itoh, Y. Endoh, T. Yokoo, D. Kawana, Y. Kaneko, Y. Tokura, and M. Fujita, Journal of the Physical Society of Japan 82, 043001 (2013).
  * Lado and Fernández-Rossier (2017) J. L. Lado and J. Fernández-Rossier, 2D Materials 4, 035002 (2017).
  * Kondo and Akagi (2022) H. Kondo and Y. Akagi, Phys. Rev. Research 4, 013186 (2022).


# Why does dissolving salt in water decrease its dielectric permittivity

Chunyi Zhang  Department of Physics, Temple University, Philadelphia,
Pennsylvania 19122, USA  Shuwen Yue  Department of Chemical and Biological
Engineering, Princeton University, Princeton, New Jersey 08544, USA
Athanassios Z. Panagiotopoulos  Department of Chemical and Biological
Engineering, Princeton University, Princeton, New Jersey 08544, USA  Michael
L. Klein<EMAIL_ADDRESS>Department of Physics, Temple University,
Philadelphia, Pennsylvania 19122, USA  Institute for Computational Molecular
Science, Temple University, Philadelphia, Pennsylvania 19122, USA  Department
of Chemistry, Temple University, Philadelphia, Pennsylvania 19122, USA  Xifan
Wu<EMAIL_ADDRESS>Department of Physics, Temple University, Philadelphia,
Pennsylvania 19122, USA  Institute for Computational Molecular Science, Temple
University, Philadelphia, Pennsylvania 19122, USA

###### Abstract

The dielectric permittivity of salt water decreases on dissolving more salt.
For nearly a century, this phenomenon has been explained by invoking
saturation in the dielectric response of the solvent water molecules. Herein,
we employ an advanced deep neural network (DNN), built using data from density
functional theory, to study the dielectric permittivity of sodium chloride
solutions. Notably, the decrease in the dielectric permittivity as a function
of concentration, computed using the DNN approach, agrees well with
experiments. Detailed analysis of the computations reveals that the dominant
effect, caused by the intrusion of ionic hydration shells into the solvent
hydrogen-bond network, is the disruption of dipolar correlations among water
molecules. Accordingly, the observed decrease in the dielectric permittivity
is mostly due to increasing suppression of the collective response of solvent
waters.

In chemistry and biology, water is widely referred to as the universal solvent
[1, 2]. As salts dissolve in water, the anomalously large dielectric
permittivity of water promotes the solubilization of salt by screening
interionic Coulomb interactions. At the same time, the dielectric response of
water is influenced by the presence of dissolved salts [3, 4, 5, 6, 7, 8, 9].
Almost 100 years ago, it was found that the static dielectric permittivity of
sodium chloride (NaCl) solution decreases as more salt is dissolved [3].
Later, more sophisticated experiments revealed a nonlinear behavior in which
dielectric decrement slows down at high solute concentrations [4, 8, 9]. A
theoretical explanation of this phenomenon was conceived soon after the first
experiment. As stated in their dielectric saturation theory, Debye [10] and
Sack [11] envisioned the formation of hydration shells due to the tendency of
water dipoles to be aligned along electric fields of dissociated ions. Debye
further estimated that ionic electric fields are strong enough to saturate the
polarizability of water molecules near the ions and therefore lower the
dielectric response [10]. Because of its built-in physical intuition,
dielectric saturation has been, to date, the most adopted theory to explain
dielectric decrement in salt water [11, 12, 13, 14, 15].

The past half-century has witnessed significant progress in understanding
water through principles of quantum mechanics and statistical physics [16, 17,
18, 19, 20]. This progress calls into question the dielectric saturation
explanation. Indeed, consensus has been reached that the high dielectric
permittivity of water is closely associated with correlated dipole
fluctuations of water molecules on the underlying hydrogen(H)-bond network
[21, 22, 23]. However, this collective dipolar response is missing in the
picture of dielectric saturation which mainly concerns the suppressed
dielectric response of individual water molecules [11, 10, 13]. More
disturbingly, based on classical electrodynamics, dielectric saturation is
estimated to occur on water molecules that are a few angstroms away from ions
[10]. The above length scale is comparable to the estimated de Broglie
wavelength of electrons at room temperature [24]. Physical interactions at
such length scales are governed by quantum mechanics rather than a classical
description. In this regard, density functional theory (DFT)-based [25, 26] ab
initio molecular dynamics (AIMD) [27] provides an ideal framework to predict
properties of liquids from quantum mechanical principles. Indeed, recent AIMD
simulations found that polarizabilities of water molecules in ionic first
hydration shells are only slightly different from that in neat water [28, 29],
which contradicts the dielectric saturation hypothesis.

Due to the long-range nature of the dipole-dipole interaction and the
disordered liquid structure, the prediction of dielectric response in water
demands both a spatially extensive model containing many hundreds of water
molecules and a simulation time beyond nanoseconds [30, 31]. However, AIMD
simulations of such large timescale and system size are simply not feasible
using current computer architectures. Thus, to date, dielectric decrement has
been mostly studied using molecular dynamics with classical force fields, and
the effect of electronic polarizability has been neglected [13, 14, 32, 33].

Herein, we overcome the challenge by studying dielectric decrement by
combining AIMD and deep neural networks (DNNs) [34, 35, 36]. The liquid
structures of NaCl solutions are simulated by a DNN that explicitly
incorporates long-range electrostatic interactions [34] with periodic
simulation cells containing about 4000 water molecules. Importantly, the
potential is trained on DFT calculations based on the strongly constrained
appropriately normed (SCAN) functional [37, 38]. In addition, a second DNN
[39] is trained separately for centers of electronic orbitals, in terms of
maximally localized Wannier functions [40]. Notably, this second DNN allows us
to rigorously partition the electronic charge density into contributions from
dipole moments of individual water molecules. The dual DNNs enable efficient
computations of dielectric permittivity at the DFT accuracy. (See Supplemental
Material [41] for more details on this methodology.)

Figure 1: Static dielectric permittivity of NaCl solutions,
$\varepsilon_{\rm{NaCl(aq)}}$, from this work and experiments [8, 9] [42]. All
results are normalized by the dielectric permittivity of neat water
$\varepsilon_{\rm{water}}$. The inset shows the decomposition of the computed
dielectric permittivity.

Based on linear response theory, the static dielectric permittivity of NaCl
solutions, $\text{\large{$\varepsilon$}}_{\mathrm{NaCl(aq)}}$, is related to
the fluctuation of the total dipole moment, $\bm{\mathrm{M}}$, by [43, 13]

$\displaystyle\text{\large{$\varepsilon$}}_{\mathrm{NaCl(aq)}}$
$\displaystyle=\frac{{\langle}\bm{\mathrm{M}}^{2}{\rangle}}{3Vk_{B}T\text{\large{$\varepsilon$}}_{0}}+\text{\large{$\varepsilon$}}_{\infty}$
$\displaystyle=\frac{{\langle}(\bm{\mathrm{M}}_{\rm{W(aq)}}+\bm{\mathrm{M}}_{\rm{I(aq)}})^{2}{\rangle}}{3Vk_{B}T\text{\large{$\varepsilon$}}_{0}}+\text{\large{$\varepsilon$}}_{\infty}$
$\displaystyle=\text{\large{$\varepsilon$}}_{\mathrm{W(aq)}}+\text{\large{$\varepsilon$}}_{\mathrm{W(aq)-I(aq)}}+\text{\large{$\varepsilon$}}_{\mathrm{I(aq)}}+\text{\large{$\varepsilon$}}_{\infty}$
(1)

where $V$, $k_{B}$, $T$, and $\text{\large{$\varepsilon$}}_{0}$ are the system
volume, Boltzmann constant, temperature, and vacuum permittivity,
respectively. $\text{\large{$\varepsilon$}}_{\infty}$ is the electronic
contribution in the high-frequency limit. As expected, the theoretical
$\text{\large{$\varepsilon$}}_{\infty}$ are small values around 1.88-1.99 at
concentrations under consideration. We report the computed dielectric
permittivity of NaCl solutions in Fig. 1 together with experimental data. Note
that both results have been normalized to enable a better comparison of
dielectric decrement behavior. There is good agreement between experiments and
present calculations. In particular, the nonlinear behavior in dielectric
decrement observed in experiment is well reproduced. The dielectric
permittivity drops steeply at low concentrations, but its slope becomes
gradually flattened as solute concentration increases. Notably, the
nonlinearity generates a bowing feature in dielectric decrement. Absolute
values of the computed dielectric permittivity are reported in Supplemental
Material Table 1 [41]. It should be noted that the predicted dielectric
permittivity of neat water by SCAN functional is 102.5, which is larger than
the experimental value of 78. The overestimation of the dielectric
permittivity is consistent with a previous study employing the SCAN functional
[30], and this overestimation is particularly attributed to the self-
interaction error in the SCAN functional that over-strengthens H-bonds. The
slightly overstructured liquid water has been widely reported in literature
[17, 44, 45] and its effects on observables can be approximated by the effects
of decreasing the temperature, which does not affect our conclusions.

In NaCl solutions, the fluctuation of the overall dipole moment,
$\bm{\mathrm{M}}$, involves contributions from both water molecules,
$\bm{\mathrm{M}}_{\rm{W(aq)}}$, and ions, $\bm{\mathrm{M}}_{\rm{I(aq)}}$.
Therefore, the dielectric permittivity,
$\text{\large{$\varepsilon$}}_{\mathrm{NaCl(aq)}}$ in Eq. Why does dissolving
salt in water decrease its dielectric permittivity is composed of the self-
terms, $\text{\large{$\varepsilon$}}_{\mathrm{W(aq)}}$ and
$\text{\large{$\varepsilon$}}_{\mathrm{I(aq)}}$ whose dipole fluctuations are
restricted to water molecules and solvated ions only, and the cross-coupling
term $\text{\large{$\varepsilon$}}_{\mathrm{W(aq)-I(aq)}}$ reflecting dipole
fluctuations in water induced by the movements of ions or vice versa. The
computed values of above terms are presented in the inset of Fig. 1. Notably,
$\text{\large{$\varepsilon$}}_{\mathrm{NaCl(aq)}}$ is dominated by
$\text{\large{$\varepsilon$}}_{\mathrm{W(aq)}}$ at all concentrations, which
agrees with previous findings [13, 46]. Thus, dielectric decrement observed in
NaCl solutions is due to the weakened dielectric response of solvent water
molecules.

The dielectric component $\text{\large{$\varepsilon$}}_{\mathrm{W(aq)}}$ due
to solvent water can be further evaluated via the dipolar correlation
formalism proposed by Kirkwood [21] as

$\text{\large{{\hbox{\varepsilon}}}}_{\mathrm{W(aq)}}=\frac{\rho\mu^{2}G_{K}}{3k_{B}T\text{\large{{\hbox{\varepsilon}}}}_{0}},$
(2)

where $\rho$ and $\mu$ denote water number density and average dipole moment
per water molecule respectively, and $G_{K}$ is the so-called correlation
factor that measures the total angular correlations among water dipoles. In
polar liquids, $G_{K}$ is obtained by the integration of the dipolar
correlation function as
$G_{K}=\int\mathcal{C}(\bm{r})d\bm{r}=\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{N}\bm{\hat{\mu}}_{i}\cdot\bm{\hat{\mu}}_{j}$,
where $\bm{\hat{\mu}}_{i}$ is the unit vector of the $i$th molecular dipole
and $N$ is the number of water molecules. The dipolar correlation is defined
as $\mathcal{C}(\bm{r})={\langle}\bm{d}(\bm{0})\cdot\bm{d}(\bm{r}){\rangle}$,
accounting for the spatial correlation between the dipolar density as a
function of distance, $\bm{r}$. Because of the discretized nature of water
molecules, the dipolar density is defined as
$\bm{d}(\bm{r})=\sum_{i=1}^{N}\bm{\hat{\mu}}_{i}\delta(\bm{r}-\bm{r}_{i})$
with $\bm{r}_{i}$ denoting the position vector of the $i$th water molecule.

In neat water, both the dipole moment, $\mu$, and the correlation factor,
$G_{K}$, are largely enhanced by the underlying H-bond network, leading to the
anomalously large dielectric permittivity [23]. In NaCl solutions, as shown in
Fig. 2 (a), the correlation factor, $G_{K}$, the water number density, $\rho$,
and the water dipole moment, $\mu$, all decrease as increasing amounts of salt
dissolved, which according to Eq. 2 leads to dielectric decrement.

Figure 2: (a) Decrease of the correlation factor, $G_{K}$, water number
density, $\rho$, and the square of water molecular dipole, $\mu{{}^{2}}$, as a
function of concentration. All three quantities are normalized by their
corresponding values at zero concentration. (b) Schematic diagram of the
intrusion of ionic hydration shells into the tetrahedral H-bond network.

The effect from the disrupted H-bond network

As seen in Fig. 2 (a), dielectric decrement of NaCl solutions is mostly
attributed to the decreased correlation factor, $G_{K}$, relative to that of
neat water. Thus, the strong correlation among dipole moments in neat water is
significantly suppressed in salt solutions. In neat water, the large $G_{K}$
is closely associated with the tetrahedral H-bond structure, in which a water
molecule at the center of a tetrahedron is H-bonded with four neighboring
water molecules. The directions of dipole moments of any two H-bonded water
molecules, therefore, point in a similar direction, resulting in a positive
$\bm{\hat{\mu}}_{i}\cdot\bm{\hat{\mu}}_{j}$, which gives rise to the first
positive sharp peak at 2.7 Å in the dipolar correlation function in Fig. 3
(b). Under the influence of the directional H-bonding, dipole moments on
vertices of a tetrahedron also prefer to be aligned in a similar direction to
some extent, which yields a second positive peak around 5.1 Å in Fig. 3 (b).
In the same fashion, the dipolar correlation propagates to the third
coordination shell and beyond.

Figure 3: (a) Decomposition of the loss of correlation factor, $G_{K}$, of
NaCl solutions relative to neat water into contributions from correlations
between “bulk water”, $G_{K}^{\rm{B}}$, “hydration water”, $G_{K}^{\rm{H}}$,
and cross-correlations between “bulk water” and “hydration water”,
$G_{K}^{\rm{BH}}$. (b) The dipolar correlation function $\mathcal{C}(r)$ of
neat water and 4.9 M NaCl solution. Decomposition of the difference of
correlation factor, $\Delta\mathcal{C}(r)$, between 4.9 M NaCl solution and
neat water into the contribution from (c) the correlation between “hydration
water” $\Delta\mathcal{C}^{\rm{H}}(r)$ and (d) the cross-correlation between
“bulk water” and “hydration water” $\Delta\mathcal{C}^{\rm{BH}}(r)$. Insets in
(b-d) schematically show molecular configurations with arrows representing
some representative molecular dipoles: H-bonded water dipoles in (b) point in
similar directions contributing to a positive dipolar correlation; the
presence of ions in (c) and (d) introduces negative dipolar correlations, as
shown by the opposite directions of water molecular dipoles.

The H-bond network is disrupted increasingly as more salt is dissolved. Salt
ions exert electrostatic fields that can attract water molecules by competing
with the H-bonding. In the close vicinity of ions, water molecules hydrate the
ions by orienting their electric dipole moments towards the ions, thereby
lowering the electrostatic energy of the system, as schematically shown in
Fig. 2 (b). For a sodium cation, the first hydration shell can be described as
a relatively tight sphere comprised of about 5 or 6 water molecules, whose
oxygen is attractive to the cation at the center [47]. On the other hand, the
first hydration shell of a chloride ion is a relatively large sphere composed
of as many as 6-8 water molecules whose protons are attracted to the chloride
lone pair electrons [47, 48].

Because of the intrusion of the hydration shells, water molecules in the
solvent are now divided into two distinct categories: the “hydration (H)
water” inside the ionic hydration shells and the “bulk (B) water” outside it.
As such, the pattern of dipolar correlation is fundamentally revised. As shown
in Eq. Why does dissolving salt in water decrease its dielectric permittivity,

$\displaystyle G_{K}$
$\displaystyle=\int[\mathcal{C}^{\rm{B}}(\bm{r})+\mathcal{C}^{\rm{H}}(\bm{r})+\mathcal{C}^{\rm{BH}}(\bm{r})]d\bm{r}$
$\displaystyle=G_{K}^{\rm{B}}+G_{K}^{\rm{H}}+G_{K}^{\rm{BH}},$ (3)

the total correlation factor $G_{K}$ involves the self-terms of
$G_{K}^{\rm{B}}$ ($G_{K}^{\rm{H}}$) by dipolar correlation restricted to “bulk
water” (“hydration water”) only, and the coupling term $G_{K}^{\rm{BH}}$ due
to the dipolar correlation between “bulk water” and “hydration water”. The
above components in correlation factors, relative to neat water, are presented
in Fig. 3 (a). (See: Supplemental Material [41] for more details.)

As seen in Fig. 3 (a), the reduction in the overall correlation factor,
$G_{K}$, is mostly from $G_{K}^{\mathrm{H}}$, which describes the correlation
among “hydration water”. This is because water molecules in hydration shells
are constricted by the ion-water attraction instead of H-bonding. Within a
hydration shell, the cation (anion)-water attraction reorientates the dipole
moments from an H-bonding direction to a central-force direction pointing
outwards (towards) ions. As such, the dipolar correlation between two
neighboring “hydration water” molecules is thereby significantly suppressed.
This is evidenced by the sharp negative peak at $\sim$ 2.7 Å in the dipolar
correlation function $\Delta\mathcal{C}^{\rm{H}}(r)$ as plotted relative to
neat water in Fig. 3 (c). Moreover, the absence of H-bonding even causes anti-
correlations between two “hydration water” molecules located on the opposite
sides of a single ion as schematically shown by opposite directions of water
molecular dipoles in the inset of Fig. 3 (c). Therefore, the aforementioned
positive peak of neat water in Fig. 3 (b) due to correlated dipole moments on
vertices of a tetrahedron at 5.1 Å disappears. Instead, it is replaced by two
negative peaks at 4.8 and 6.1 Å, which are caused by the anti-correlated water
dipoles in hydration shells of Na+ and Cl- ions, respectively. At long range,
water molecules in a hydration shell, in principle, should be correlated to
those in another hydration shell. However, such correlations are also weaker
than those in neat water as expected in Fig. 3 (c). As concentration
increases, the loss of $G_{K}^{\mathrm{H}}$ should accumulate linearly, which
is responsible for most of the linear dielectric decrement in salt water.

Of course, “hydration water” is H-bonded to “bulk water”, and in this way, the
H-bond network is partially restored. Nevertheless, the reconstructed H-bond
structure deviates from that found in neat water. Within a hydration shell,
two water molecules located on opposite sides of a single ion are anti-
correlated, as mentioned above. Because of the highly directional nature of
H-bonding, the anti-correlation extends to the correlation between one
“hydration water” molecule and one “bulk water” molecule that is H-bonded to
another “hydration water” molecule at the other side of the ion, as
schematically shown by the opposite direction of green arrows in the inset of
Fig. 3 (d). Again, these anticorrelations can be identified as a broad
negative peak centered at 8 Å, which weakens the dipolar correlation. As a
result, $G_{K}^{\mathrm{BH}}$ also contributes to the decreased overall
correlation factor of $G_{K}$ relative to neat water, as shown in Fig. 3 (a).
Moreover, $G_{K}^{\mathrm{BH}}$ plays a surprisingly key role in the nonlinear
dielectric decrement as evidenced by its arc shape in Fig. 3 (a). This
nonlinearity is an intrinsic property because $G_{K}^{\mathrm{BH}}$ describes
the correlation between the dipolar density of “bulk water”
$\bm{d}^{\mathrm{B}}(\bm{r})$ and the dipolar density of “hydration water”
$\bm{d}^{\mathrm{H}}(\bm{r})$, and its value depends on the existence of both
types of water, i.e.,
${\langle}\bm{d}^{\mathrm{B}}(\bm{0})\cdot\bm{d}^{\mathrm{H}}(\bm{r}){\rangle}$.
In neat water, $G_{K}^{\mathrm{BH}}=0$ since the dipolar density of “hydration
water” $\bm{d}^{\mathrm{H}}(\bm{r})$ is zero. As salt dissolve in water,
hydration shells appear in the solution, and the absolute value of
$G_{K}^{\mathrm{BH}}$ starts to increase, reaching its maximum at about 2.3 M,
in which the NaCl solution is roughly equally occupied by “bulk water” and
“hydration water”. After the maximum, $G_{K}^{\mathrm{BH}}$ decreases with
further elevated concentrations. In principle, it will vanish again at
$\bm{d}^{\mathrm{B}}(\bm{r})=0$, when the entire solution is completely
occupied by hydration shells.

The tetrahedral H-bond network is expected to recover in the “bulk water”
outside the hydration shell. The dipolar correlation among “bulk water”
molecules is captured by the $G_{K}^{\mathrm{B}}$ component of the correlation
factor. Indeed, the analysis in Fig. 3 (a) shows that $G_{K}^{\mathrm{B}}$ of
NaCl solutions at all concentrations is little different from neat water.
Thus, the large decrease in the correlation factor, $G_{K}$, in salt water is
mostly due to the disrupted H-bond network in the “hydration water”.

Excluded volume effect Due to short-range repulsion, ions and water molecules
are separated by 2-4 Å. This extra volume demanded by ions is no longer
accessible to water molecules, and the water number density is therefore
decreased. In the literature, this is referred to as the excluded volume
effect [49]. According to Eq. 2, this effect should lead to the decreased
dielectric permittivity. Indeed, the present computations show that the
excluded volume effect makes a small contribution to dielectric decrement, in
which the water number density decreases slightly with increasing solute
concentration as shown in Fig. 2 (a). Since the repelled volume by ions is
proportional to the salt concentration, dielectric decrement due to the
excluded volume effect is indeed linear as expected.

Local field effect Hydrated ions, like all charged defects, change the
electrostatic potential profile throughout the solution. As expected, water
molecules nearby an ion are polarized in a different manner from neat water.
In condensed matter physics, related phenomena have been already identified,
for example around defects in semiconductors or at interfaces in solid
materials, and they have long been recognized as the local field effect [24].
There is consensus that a proper description of local field effects,
particularly for regions close to charged defects, demands electronic
structure details computed from quantum mechanics. Based on DFT, the present
DNN simulations yield a dipole moment, $\mu$=2.85 (2.91) Debye for the
“hydration water” of the cation (anion), which is only slightly smaller than
the value of 2.99 Debye in neat water. This suggests that the capability of
ions to polarize the water dipole is comparable to that of H-bonding. Indeed,
it is also consistent with the recent theoretical discovery that molecular
polarizabilities of the “hydration water” are only marginally different from
that in neat water [28, 29]. Since H-bonding is mostly electrostatic in
nature, it strongly indicates that water molecules nearby ions are far from
being saturated by ions’ local fields. Nevertheless, the local field effect
also contributes slightly to dielectric decrement as indicated by Eq. 2.
Because the $\mu$ of the “hydration water” is only a little smaller than in
neat water, $\mu^{2}$ of NaCl solutions drops slowly as a function of
concentration, as shown in Fig. 2 (a).

In addition to the SCAN ab initio simulations, we also simulated the
dielectric permittivity using the classical OPC water model [50]. As shown in
Supplemental Material [41], the results obtained using the OPC model agree
well with those from the SCAN-DFT approach. A notable distinction between the
OPC model and the SCAN-DFT model is that the OPC model is a rigid model with a
fixed dipole moment of 2.48 D, indicating that the DFT approach is necessary
for accurately capturing the local field effect.

In conclusion, dielectric decrement, as a century-old problem, has been
extensively studied over decades. However, a critical question remains
unresolved in the field regarding the main origin behind the dielectric
decrement—whether it is the dielectric saturation effect [10, 11] or the loss
of dipolar correlation on the H-bond network [13, 14]. To provide an
unambiguous answer, theoretical simulations must explicitly include both a
polarizable model of water molecules and an accurate model of H-bonding, which
can account for the dielectric saturation effect and correlation effect
simultaneously. Importantly, the polarizable models of water molecules should
be described from first principles at the quantum mechanics level, because the
length scale of dielectric saturation effect is about a few angstroms which is
comparable to the de Broglie wavelength of electrons at room temperature. In
this work, we achieve the above goal by reproducing dielectric decrement in
NaCl solutions on the DFT level using advanced DNNs. The results unambiguously
determine that the dielectric decrement in NaCl solutions is dominated by the
loss of correlations between water molecules due to the intrusion of ionic
hydration shells into the H-bond network, while the contribution from
dielectric saturation effect is small. Importantly, the present computations
provide a quantitative explanation of dielectric decrement in salt water; we
found that the linear dielectric decrement is due to the loss of correlation
within hydration shells, while nonlinear dielectric decrement is due to the
loss of correlation between water in hydration shells and bulk water.

###### Acknowledgements.

We thank Roberto Car, Linfeng Zhang, and Han Wang for fruitful discussions.
This work was supported by National Science Foundation through Awards No.
DMR-2053195. We also acknowledge support from the “Chemistry in Solution and
at Interfaces” (CSI) Center funded by the U.S. Department of Energy through
Award No. DE-SC0019394. This research used resources of the National Energy
Research Scientific Computing Center (NERSC), which is supported by the U.S.
Department of Energy (DOE), Office of Science under Contract No. DE-
AC02-05CH11231. This research includes calculations carried out on HPC
resources supported in part by the National Science Foundation through major
research instrumentation grant number 1625061 and by the U.S. Army Research
Laboratory under contract No. W911NF-16-2-0189. This research used resources
of the Oak Ridge Leadership Computing Facility at the Oak Ridge National
Laboratory, which is supported by the Office of Science of the U.S. Department
of Energy under Contract No. DE-AC05-00OR22725.

## References

  * Franks [2000] F. Franks,  _Water: a matrix of life_ , 2nd ed. (Royal Society of Chemistry, Cambridge, UK, 2000).
  * Eisenberg  _et al._ [2005] D. Eisenberg, W. Kauzmann, and W. Kauzmann,  _The Structure and Properties of Water_ (OUP Oxford, 2005).
  * Blüh [1924] O. Blüh, Z. Physik 25, 220 (1924).
  * Hasted  _et al._ [1948] J. B. Hasted, D. M. Ritson, and C. H. Collie, J. Chem. Phys. 16, 1 (1948).
  * Haggis  _et al._ [1952] G. Haggis, J. Hasted, and T. Buchanan, J. Chem. Phys. 20, 1452 (1952).
  * Harris and O’Konski [1957] F. E. Harris and C. T. O’Konski, J. Phys. Chem. 61, 310 (1957).
  * Hasted and Roderick [1958] J. B. Hasted and G. W. Roderick, J. Chem. Phys. 29, 17 (1958).
  * Christensen  _et al._ [1966] J. H. Christensen, A. J. Smith, R. B. Reed, and K. L. Elmore, J. Chem. Eng. Data 11, 60 (1966).
  * Buchner  _et al._ [1999] R. Buchner, G. T. Hefter, and P. M. May, J. Phys. Chem. A 103, 1 (1999).
  * Debye [1929] P. J. W. Debye,  _Polar molecules_ (Dover publications, 1929).
  * Sack [1926] H. Sack, Phys. Z. 27, 206 (1926).
  * Stiles  _et al._ [1982] P. Stiles, J. Hubbard, and R. Kayser, J. Chem. Phys. 77, 6189 (1982).
  * Rinne  _et al._ [2014] K. F. Rinne, S. Gekle, and R. R. Netz, J. Chem. Phys. 141, 214502 (2014).
  * Seal  _et al._ [2019] S. Seal, K. Doblhoff-Dier, and J. Meyer, J. Phys. Chem. B 123, 9912 (2019).
  * Gorji and Bowler [2020] A. Gorji and N. Bowler, J. Chem. Phys. 153, 014503 (2020).
  * Laasonen  _et al._ [1993] K. Laasonen, M. Sprik, M. Parrinello, and R. Car, J. Chem. Phys. 99, 9080 (1993).
  * Chen  _et al._ [2017] M. Chen, H.-Y. Ko, R. C. Remsing, M. F. Calegari 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).
  * Zhang  _et al._ [2021] C. Zhang, F. Tang, M. Chen, J. Xu, L. Zhang, D. Y. Qiu, J. P. Perdew, M. L. Klein, and X. Wu, J. Phys. Chem. B 125, 11444 (2021).
  * O’Neill  _et al._ [2022] N. O’Neill, C. Schran, S. J. Cox, and A. Michaelides, arXiv preprint arXiv:2211.04345 (2022).
  * Schran  _et al._ [2021] C. Schran, F. L. Thiemann, P. Rowe, E. A. Müller, O. Marsalek, and A. Michaelides, Proc. Natl. Acad. Sci. U.S.A. 118, e2110077118 (2021).
  * Kirkwood [1939] J. G. Kirkwood, J. Chem. Phys. 7, 911 (1939).
  * Pople [1951] J. A. Pople, Proc. R. Soc. Lond. A 205, 163 (1951).
  * Sharma  _et al._ [2007] M. Sharma, R. Resta, and R. Car, Phys. Rev. Lett. 98, 247401 (2007).
  * Milovzorov [2018] D. E. Milovzorov, in  _Electric Field_, edited by M. S. Kandelousi (InTech, 2018).
  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Car and Parrinello [1985] R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985).
  * Gaiduk and Galli [2017] A. P. Gaiduk and G. Galli, J. Phys. Chem. Lett. 8, 1496 (2017).
  * Rozsa  _et al._ [2020] V. Rozsa, T. A. Pham, and G. Galli, J. Chem. Phys. 152, 124501 (2020).
  * 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, Phys. Rev. Lett. 126, 216403 (2021).
  * Zhang and Sprik [2016] C. Zhang and M. Sprik, Phys. Rev. B 93, 144201 (2016).
  * Zasetsky and Svishchev [2001] A. Y. Zasetsky and I. M. Svishchev, J. Chem. Phys. 115, 1448 (2001).
  * Fuentes-Azcatl and Barbosa [2016] R. Fuentes-Azcatl and M. C. Barbosa, J. Phys. Chem. B 120, 2460 (2016).
  * Zhang  _et al._ [2022a] L. Zhang, H. Wang, M. C. Muniz, A. Z. Panagiotopoulos, R. Car, and W. E, J. Chem. Phys. 156, 124107 (2022a).
  * Zhang  _et al._ [2018a] L. Zhang, J. Han, H. Wang, R. Car, and W. E, Phys. Rev. Lett. 120, 143001 (2018a).
  * Zhang  _et al._ [2018b] L. Zhang, J. Han, H. Wang, W. Saidi, R. Car, and W. E, in  _Advances in Neural Information Processing Systems 31_, edited by S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (Curran Associates, Inc., 2018) pp. 4436–4446.
  * Sun  _et al._ [2015] J. Sun, A. Ruzsinszky, and J. P. Perdew, Phys. Rev. Lett. 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, _et al._ , Nat. Chem. 8, 831 (2016).
  * Zhang  _et al._ [2020a] L. Zhang, M. Chen, X. Wu, H. Wang, W. E, and R. Car, Phys. Rev. B 102, 041121 (2020a).
  * Marzari and Vanderbilt [1997] N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997).
  * [41] See Supplemental Material [URL] for methodology and numerical details and the dielectric permittivity simulated using the OPC water model, which includes Refs. [51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65].
  * Shcherbakov  _et al._ [2014] V. Shcherbakov, Y. M. Artemkina, and E. Korotkova, Russ. J. Inorg. Chem. 59, 922 (2014).
  * Sala  _et al._ [2010] J. Sala, E. Guàrdia, and J. Martí, J. Chem. Phys. 132, 214505 (2010).
  * Xu  _et al._ [2020] J. Xu, C. Zhang, L. Zhang, M. Chen, B. Santra, and X. Wu, Phys. Rev. B 102, 214113 (2020).
  * Piaggi  _et al._ [2021] P. M. Piaggi, A. Z. Panagiotopoulos, P. G. Debenedetti, and R. Car, J. Chem. Theory Comput. 17, 3065 (2021).
  * Chandra [2000] A. Chandra, J. Chem. Phys. 113, 903 (2000).
  * Zhang  _et al._ [2022b] C. Zhang, S. Yue, A. Z. Panagiotopoulos, M. L. Klein, and X. Wu, Nat. Commun. 13, 822 (2022b).
  * Xu  _et al._ [2021] J. Xu, Z. Sun, C. Zhang, M. DelloStritto, D. Lu, M. L. Klein, and X. Wu, Phys. Rev. Materials 5, L012801 (2021).
  * Adar  _et al._ [2018] R. M. Adar, T. Markovich, A. Levy, H. Orland, and D. Andelman, J. Chem. Phys. 149, 054504 (2018).
  * Izadi  _et al._ [2014] S. Izadi, R. Anandakrishnan, and A. V. Onufriev, J. Phys. Chem. Lett. 5, 3863 (2014).
  * Wang  _et al._ [2018] H. Wang, L. Zhang, J. Han, and W. E, Comput. Phys. Commun. 228, 178 (2018).
  * Zhang  _et al._ [2020b] Y. Zhang, H. Wang, W. Chen, J. Zeng, L. Zhang, H. Wang, and W. E, Comput. Phys. Commun. 253, 107206 (2020b).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, and et al., J. Phys.: Condens. Matter 21, 395502 (2009).
  * Hamann  _et al._ [1979] D. Hamann, M. Schlüter, and C. Chiang, Phys. Rev. Lett. 43, 1494 (1979).
  * Vanderbilt [1985] D. Vanderbilt, Phys. Rev. B 32, 8412 (1985).
  * Hamann [2013] D. R. Hamann, Phys. Rev. B 88, 085117 (2013).
  * Pizzi  _et al._ [2020] G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, and et al., J. Phys. Condens. Matter. 32, 165902 (2020).
  * Plimpton [1995] S. Plimpton, J. Chem. Phys. 117, 1 (1995).
  * Hockney and Eastwood [1989] R. W. Hockney and J. W. Eastwood,  _Computer Simulation Using Particles_ (Adam Hilger, NY, 1989).
  * Spaldin [2012] N. A. Spaldin, J. Solid State Chem. 195, 2 (2012).
  * Vanderbilt and King-Smith [1993] D. Vanderbilt and R. D. King-Smith, Phys. Rev. B 48, 4442 (1993).
  * Silvestrelli and Parrinello [1999] P. L. Silvestrelli and M. Parrinello, Phys. Rev. Lett. 82, 3308 (1999).
  * Baroni  _et al._ [2001] S. Baroni, S. De Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
  * Soper [1997] A. K. Soper, J. Phys.: Condens. Matter 9, 2399 (1997).
  * Sengupta  _et al._ [2021] A. Sengupta, Z. Li, L. F. Song, P. Li, and K. M. Merz, J. Chem. Inf. Model. 61, 869 (2021).