# Spatio-temporal modeling of saltatory conduction in neurons using Poisson-
Nernst-Planck treatment and estimation of conduction velocity

Rahul Gulati  Shiva Rudraraju<EMAIL_ADDRESS>Department of
Mechanical Engineering, University of Wisconsin-Madison, WI, USA

###### Abstract

Action potential propagation along the axons and across the dendrites is the
foundation of the electrical activity observed in the brain and the rest of
the central nervous system. Theoretical and numerical modeling of this action
potential activity has long been a key focus area of electro-chemical neuronal
modeling, and over the years, electrical network models of varying complexity
have been proposed. Specifically, considering the presence of nodes of Ranvier
along the myelinated axon, single-cable models of the propagation of action
potential have been popular. Building on these models, and considering a
secondary electrical conduction pathway below the myelin sheath, the double-
cable model has been proposed. Such cable theory based treatments, including
the classical Hodgkin-Huxley model, single-cable model, and double-cable model
have been extensively studied in the literature. But these have inherent
limitations in their lack of a representation of the spatio-temporal evolution
of the neuronal electro-chemistry. In contrast, a Poisson-Nernst-Planck (PNP)
based electro-diffusive framework accounts for the underlying spatio-temporal
ionic concentration dynamics and is a more general and comprehensive
treatment. In this work, a high-fidelity implementation of the PNP model is
demonstrated. This electro-diffusive model is shown to produce results similar
to the cable theory based electrical network models, and in addition, the rich
spatio-temporal evolution of the underlying ionic transport is captured. Novel
to this work is the extension of PNP model to axonal geometries with multiple
nodes of Ranvier, its correlation with cable theory based models, and multiple
variants of the electro-diffusive model - PNP without myelin, PNP with myelin,
and PNP with the myelin sheath and peri-axonal space. Further, we apply this
spatio-temporal model to numerically estimate conduction velocity in a rat
axon using the three model variants. Specifically, spatial saltatory
conduction due to the presence of myelin sheath and the peri-axonal space is
investigated.

###### keywords:

Action potential , saltatory conduction , signal velocity , Poisson-Nernst-
Planck , electro-diffusion , neuronal electrophysiology , myelin , peri-axonal
space , cable theory

††journal: Brain Multiphysics

##  1 Introduction

Electrical activity in nerve cells, enabled through the propagation of action
potentials, is critical to the entire signaling and communication cascade of
the nervous system. Disruption of this requisite signaling can lead to a
number of neurological disorders such as the motor neuron diseases and is
often linked with traumatic brain injury (TBI), Alzheimers, cognitive
impairment, depression etc [1, 2, 3, 4, 5]. To gain insight into the neuronal
electrophysiology, varied experimental investigations using the patch-clamp
technique, electroencephalograms (EEG), electrocardiogram (ECG), MRI, calcium
imaging, voltage imaging etc have been reported in literature. Despite the
wealth of information achieved by these investigations, there is a need to
supplement these studies with a robust numerical implementation which has the
potential to represent the electrophysiology to a far greater resolution as
recorded by the experiments.

Based on numerous voltage clamp experiments on the giant squid, Hodgkin-Huxley
came up with a first mathematical model in the form of an electrical circuit
to describe the current through the neuronal membrane [6]. They quantitatively
detailed the respective ionic conductance in respect to the membrane voltage.
Huxley delineated that the action potential propagation along the axon closely
follows Ohm’s law [7]. Using cable theory, they arrived at a one dimensional
model of the action potential propagation. Based on the physiology of the
neuron, the importance of myelin sheath that surrounds the axon, has been
emphasized. Degradation of this protective covering for example with age can
lead to slowdown of the signal or even signal disruption resulting in various
diseases [8, 9]. To be able to study the effect of myelin on action potential
propagation, the Hodgkin-Huxley model was modified to incorporate the myelin
sheath and to obtain the single-cable model [10, 11]. Experimental
investigation of the action potential by Barrett and Blight in the 1980’s lead
to the discovery of an after-potential [12, 13, 14]. With the advent of
advanced microscopy techniques, the existence of a secondary electrical
conduction pathway in the peri-axonal space under the myelin sheath has been
uncovered [10].

The cable theory based models namely Hodgkin-Huxley, single-cable, double-
cable, etc., have greatly contributed to our understanding of the neuronal
electrophysiology. They provide us an excellent insight into the membrane
potential, action potential propagation along the neuron, the electric current
propagating along the axon, the effect of saltatory conduction due to the
presence of myelin sheath, the effect of the submyelin peri-axonal space
dictated by the double-cable model, the conduction velocity implied by each of
these models etc. The cable theory based models, however, have some
limitations. First, these are a one dimensional reduction of the complex
heterogeneous spatial propagation of the action potential. Secondly, the cable
theory based models fail to describe the underlying spatial ionic diffusion
and the generated electric field during the electrical conduction propagation.
Therefore, these models are not able to accurately describe the dynamics after
a prolonged electrical activity or when the diameter is relatively smaller, as
in the case of dendrites. Further, the cable theory based models cannot be
easily extended to account for the membrane microenvironment, such as
incorporating the membrane-glia interaction.

The Poisson-Nernst-Planck based electro-diffusive model has the capability to
overcome the limitations of the cable theory based models and provide a
spatio-temporal representation of the electrical potential along with the
ionic distributions [15]. The PNP model is a more generalized model which can
be reduced to the electroneutral model and can subsequently be reduced to the
one dimensional cable theory based model [16]. Qian and Sejnowski modelled one
of the first intracellular dynamics incorporating one dimensional PNP theory
[17]. PNP model has also been applied to neuron-ECM-astrocyte interations
[18]. Assuming electroneutrality, ionic dynamics have also been represented
using Kirchoff-Nernst-Planck [19, 20]. However, the assumption of
electroneutrality for nonuniform geometries is invalid [21]. Using PNP model,
the neuronal intracellular-extracellular dynamics have been represented for a
single node of Ranvier [21, 22]. It has been demonstrated that the dynamics of
the PNP model resemble to that of the cable theory at higher ion channel
density [21]. The PNP model can also be coupled with mechanics to represent
the complex neuronal mechano-electrophysiology interactions [23, 24, 25, 26].

In this work, we extend the electro-diffusive PNP model to multiple nodes of
Ranvier, enhancing our ability to study the electrophysiology in a full length
neuronal axon. We present novel variants of the PNP model based on the
discrete cable theory based models, i.e. PNP model, PNP model with
myelination, and PNP model with myelin and peri-axonal space. As an example,
we demonstrate these models by simulating action potential conduction in a rat
neuron. Spatial saltatory conduction due to the presence of myelin sheath and
the peri-axonal space is demonstrated. Finally, we provide a detailed insight
into the numerically estimated conduction velocity for a rat and squid neuron
using various representative PNP electro-diffusive models. The Finite element
(FE) method is used to discretize the set of PDEs underlying the PNP model. As
suggested by an earlier work, non-homogeneous adaptive mesh is employed [22].
Results indicate that the conduction velocity (CV) of the action potential
increases with the presence of myelin sheath and the peri-axonal space. The CV
for the PNP with myelin is comparable to the single cable network but the CV
for the PNP with myelin and peri-axonal space model does not increase
drastically as compared to the double cable model. We observe that the action
potential amplitude is lower for the PNP model when the myelin sheath is
present.

In section 2, we briefly review and illustrate the well known models based on
the one dimensional cable theory. The electro-diffusive PNP model is presented
in section 3. The mathematical formulation of the numerical framework is
elaborated in section 4. The simulation results of the various models of the
PNP are detailed in section 5. Finally, a discussion of the conduction
velocity for a rat and a squid neuron is in section 6, followed by conclusion
in section 7.

##  2 Review of electrical network models of action potential propagation

Figure 1: Electrical circuit of cable theory based models. (A) Schematics of
cable theory based Hodgkin-Huxley circuit consists of membrane capacitance,
resistance offered by the ion channels. Action potential propagates like a
soliton through the axon and is depicted on the right. Extracellular voltage
is taken as 0 mV. (B) Electrical circuit of the single-cable model considers
the presence of myelin sheath. The capacitance and resistance offered by the
myelin are explicitly modeled. (C) Electrical network of the double-cable
representation incorporates an additional cable pathway for the submyelin
peri-axonal space. Action potential jumps from one node of Ranvier to the
next, resulting in a higher conduction velocity.

###  2.1 Hodgkin-Huxley model

The classical work of Hodgkin and Huxley was a landmark model in terms of
providing deep insights into the ionic basis of action potential propagation
in nerve cells. Based on the voltage clamp experiments on the Giant Squid, the
physiology of the initiation and propagation of action potential in a neuron
was posed as a coupled set of ordinary differential equations. This electrical
network model takes into account the membrane capacitance and the ionic
currents due to the sodium ions, potassium ions and some leak current through
the respective ion channels in the neuronal membrane. The influence of the
conductance of the respective ion channels, or conversely the resistance, on
the action potential was quantitatively estimated using experimental data. The
resulting electrical circuit is depicted in Figure 1(A), and the corresponding
governing equation linking the evolution of the membrane potential with the
underlying ionic transport is the following:

$C_{m}\frac{\partial V_{m}}{\partial
t}+\bar{G}_{Na}m^{3}h(V_{m}-E_{Na})+\bar{G}_{K}n^{4}(V_{m}-E_{K})+G_{L}(V_{m}-V_{rest})=I_{inj}$
(1)

where $C_{m}$ is the membrane capacitance of the lipid bilayer, $V_{m}$ is the
membrane potential, $V_{rest}$ is the resting potential of the nerve cells and
$I_{inj}$ is the injected current through the voltage clamp experiments to
initiate the action potential in the cell. $\bar{G_{i}}$ and $E_{i}$ are the
peak conductance and the Nernst potential of the sodium or potassium ions, and
${G_{L}}$ is the conductance of the membrane leak channel.

Numerous experiments have pointed towards the effective behaviour of the
intra-cellular region to that of a resistor along the axon [27]. Utilizing
cable theory, one can then arrive at:

$\frac{1}{R_{i}}\frac{\partial^{2}V_{m}}{\partial x^{2}}=I_{inj}$ (2)

where $R_{i}$ is the resistance per unit length of the intra-cellular region
along the axon. The partial differential equation (PDE) form of the Hodgkin-
Huxley equation can then be written as:

$C_{m}\frac{\partial V_{m}}{\partial
t}+\bar{G}_{Na}m^{3}h(V_{m}-E_{Na})+\bar{G}_{K}n^{4}(V_{m}-E_{K})+G_{m}(V_{m}-V_{rest})=\frac{1}{R_{i}}\frac{\partial^{2}V_{m}}{\partial
x^{2}}$ (3)

Using the above one-dimensional PDE, one can model the propagation of action
potential along the length of the axon.

###  2.2 Single-cable model

From the physiology of nerve cells, it is well known that the glial cells
provide a protective covering to the axonal membrane. This protective covering
consists of multiple layers of myelin sheath along the axon, with gaps in
between. These axonal gaps, void of any myelin covering, are identified as the
nodes of Ranvier. These myelin lamellae play an indispensable role in the
rapid movement of the action potential, since the direct ionic exchange with
the extra-cellular medium occurs only at the nodes of Ranvier. This leads to a
local current and the action potential jumping from one node of Ranvier to the
next, more commonly referred to as saltatory conduction. The degradation of
the myelin layers is known to lead to a decline in the conduction velocity and
is also linked with various neuronal disease conditions [8, 9].

To accommodate this spatial heterogeneity of the myelin sheath, the standard
Hodgkin-Huxley treatment needs to be altered to model myelinated-axons. A
simple electrical circuit to realize this, known as the single-cable model,
has been proposed [11, 28]. The ionic exchange between the intra-cellular
regions and the extra-cellular regions takes place at the membrane, but only
at the location of the nodes of Ranvier. As an example, for a rat axon, the
span of the nodes of Ranvier is around $2.3\mu m$ and the nodes are separated
by distance of $70\mu m-100\mu m$ [10]. The single-cable model with the myelin
sheaths having their respective capacitance and resistance is shown in Figure
1(B) [10].

###  2.3 Double-cable model

The presence of the submyelin peri-axonal region has been proposed as a
potential pathway for rapid electrical conduction along the axon [10, 29]. Due
to the presence of two conduction pathways, each modeled using cable theory,
this treatment is referred to as the double-cable model. An electrical circuit
representing such a double-cable model can be seen in Figure 1(C) . The basic
governing equations for the electrical circuit at the node of Ranvier are
similar to the set of equations from the single cable model/Hodgkin-Huxley
model. In addition, using cable theory, the partial differential equation
modeling the peri-axonal space takes the form:

$\frac{1}{R_{i}}\frac{\partial^{2}V_{m}}{\partial
x^{2}}+\frac{1}{R_{pa}}\frac{\partial^{2}V_{my}}{\partial
x^{2}}=C_{my}\frac{\partial V_{my}}{\partial t}+\frac{V_{my}}{R_{my}}$ (4)

where $R_{pa}$ is the resistance per unit length in the peri-axonal space,
$C_{my}$ is the cumulative capacitance of the myelin sheath, $R_{my}$ is the
myelin resistance and $V_{my}$ is the potential in the peri-axonal region.

###  2.4 Comparison of action potential profiles for electrical network models

The main results of this manuscript will be discussed in Section 5, but it is
worthwhile to present here a brief comparison of the relative differences
between the action potential (voltage) profiles predicted by modeling the
electrical network models described above. Figure 1 presents such a comparison
of the voltage profiles. All these plots were generated by solving the
governing equations listed earlier in this section using an in-house 1D Finite
Element Method implementation, and the model input parameters available in the
literature for a giant squid axon [6] and rat axon [10] were used.

Shown as part of Figure 1(A) is a typical voltage profile predicted by the
Hodgkin-Huxley model. The voltage profile is a single spike that propagates
along the length of the axon, in this case from left to right. After the
refractory period of the propagating spike, if a potential/current
perturbation exceeding the threshold value is injected at the left boundary,
then it would result in a voltage wave traveling with the same amplitude.

The voltage profile produced by a single-cable model is shown in Figure 1(B).
As can be seen, the potential difference between the membrane potential and
the myelin potential, ($\text{V}_{m}-\text{V}_{my}$) jumps across the nodes of
Ranvier. This saltatory conduction leads to fast propagation of the voltage
envelope, i.e. the membrane potential, which travels like a soliton-like wave.
The decreased capacitance and increased resistance due to the myelin sheath
enables this mode of propagation. In this simulation, the nodes of Ranvier are
assumed to be equidistant. However, if the inter-nodal distance is varying
along the length of an axon, it will lead to an axially varying conduction
speed.

The voltage profile produced by a double-cable model, that incorporates a
secondary conduction pathway, is shown in Figure 1(C). It can be observed that
the profile of the voltage envelope is flatter than the voltage profiles
predicted by the previous two models. Like in the case of the single-cable
model, the action potential jumps from one node of Ranvier to the next, but
through the peri-axonal space, and the propagation is relatively much faster.
The peak amplitude is comparable to the peak amplitude for the Hodgkin-Huxley
model and single-cable model.

##  3 A field theoretic model of action potential propagation

Cable theory based electrical network models have proven to be fundamental to
our current understanding of the dynamics involved in action potential
propagation. However, at the core, the network models are reduced order
representations that try to capture the complex spatio-temporal variations of
the ionic transport, voltage distribution, and most importantly the membrane
structural heterogeneity into effective electrical properties like capacitance
and resistance of the membrane and the channels. As mentioned in the
introduction, a deeper investigation into spatial and temporal interactions of
the ionic transport with the membrane microstructure (ion channel and pump
distributions, myelin, glial environment, etc.) is desired, but the resulting
coupled evolution of the voltage distributions becomes necessary when more
complex phenomena like neuronal injury effects on action potential propagation
[30, 31, 32], and conditions like neuronal hyperexcitability observed with
Alzheimer’s disease [33, 34, 35] are to be modeled.
An electro-diffusive framework for modeling spatio-temporal ionic charge
distribution and the resulting voltage evolution using high-fidelity partial
differential equations (PDE) modeling coupled electrostatics and
electrochemistry can be capable of faithfully representing spatio-temporally
heterogenous evolution of the ionic and voltage distributions leading to
generation, propagation and potentially disruption of the neuronal action
potential. While not addressed in this current manuscript, such a capability
to model spatial heterogeneity and membrane geometry-action potential
interactions is evidently more important in the dendrites than the axons, due
to the complex morphology of neuronal dendritic structures and synapses. We
now propose a PDE based field theoretic implementation of the Poisson-Nerst-
Planck (PNP) framework that can model 2D/3D ionic and voltage field
distributions and their interactions with the membrane microstructure.
In this model, the transport of the respective ionic species, $c_{i}$, due to
the corresponding diffusive and electromigration flux, $\bm{F}_{i}$, is
modeled using the classical Nernst-Planck equation:

$\frac{\partial c_{i}}{\partial t}=-\nabla\cdot\bm{F}_{i}$ (5)

The ionic flux comprises of the diffusion term and the electro-migration term
as presented in Eq. 6. Here $D_{i}$ is the diffusion coefficient of the
$i^{th}$ ion, $R$ is the gas constant, $T$ is the temperature, $F$ is the
faraday constant and $z_{i}$ is the valency of the $i^{th}$ ionic species. It
is to be noted that the cable theory based models only account for the
electro-migration term i.e. the voltage gradient term. While the diffusion
term may be neglected for axons with larger diameter, for finer geometries
like dendrites, the diffusion term takes over the electro-migration term [17].
Hence, the PNP treatment has a higher-fidelity in representing the underlying
action potential propagation. The total flux is given by:

$\bm{F}_{i}=-[D_{i}(\nabla c_{i}+\frac{c_{i}Fz_{i}}{RT}\nabla V)]$ (6)

The coupling between the voltage and local ionic distribution is modeled with
Poisson’s equation [36], shown below.

$-\nabla\cdot(\epsilon\nabla V)=F\sum_{i=1}z_{i}c_{i}$ (7)

where $\epsilon$ is the permittivity of the medium. In this work, we consider
three ionic species: $Na^{+}$, $K^{+}$ and $Cl^{-}$. The presence of an anion
species ensures regulation of the required potential - for instance at the
resting state. The sodium and potassium ions, as is well understood,
depolarise and re-polarise the neuron. While the above PDE formulation by
itself is well known, the novelty of its implementation in this work is from
its application to model field variations and interactions in the presence of
geometric heterogeneity of the nodes of Ranvier, myelin distribution, and the
peri-axonal space.

### Cable theory models as special cases of the general PNP framework

As advancements in high resolution imaging techniques lead to an improved
understanding of the neuronal membrane microenvironment, various electrical
network models incorporating the membrane spatial heterogeneity have been
proposed over the years. Of these models, a primary classification in terms of
increasing complexity of the membrane heterogeneity leads to the three
treatments described above: Hodgkin-Huxley model, single-cable model and
double-cable model. The PNP model is capable of representing the spatial
heterogeneity of each of these models, and thus reproduce the action potential
conduction profiles observed with these electrical network models. To the best
of our knowledge, numerical demonstration of this equivalence for each of
these three models, and the extraction of equivalent 1D electrical network
results (effective capacitance, effective resistance, action potential
profile, etc.) from a more general electro-diffusion framework, has not been
shown earlier in the literature. Figure 2 depicts this equivalence between:

  * 1.

PNP model and the classical Hodgkin-Huxley model,

  * 2.

PNP model with myelination and the single-cable model,

  * 3.

PNP model with myelin plus peri-axonal space and the double-cable model.

The axonal membrane has ion channels and the embedded membrane capacitance in
the PNP model similar to the Hodgkin-Huxley model. The existence of the myelin
sheath leads to the presence of nodes of Ranvier and decreased internodal
capacitance in the PNP model with myelin. Note that there is no direct
correlation of $V_{my}$ as in the single cable model. Finally, the PNP model
with myelin and peri-axonal space incorporates the sub-myelin peri-axonal
space, leading to lower net resistance and faster propagation speed of the
action potential. Here, the axial resistance and the peri-axonal resistance
act in parallel, and hence the net axial resistance is lower. While the peri-
axonal resistance does not correlate to the axial resistance in the double
cable model, this is not the case for the equivalent PNP model as the
resistance is computed using the underlying ionic constants. Here, the
diffusion constants of the underlying ions in the peri-axonal space are
assumed to be the same as in the cytoskeleton region.

Figure 2 also shows, as an overlaid line plot in each case, the 1D action
potential profile (obtained as a line-out along the axonal axis from the 2D
PNP simulations) for the three PNP model variants listed above. It can be
noticed that the curvature of the voltage potential is lower when myelin is
present as it leads to faster action potential propagation. Further, the
presence of myelin leads to a lower action potential amplitude. This maybe due
to the fact that the action potential jumps node to node as soon the threshold
potential is attained.

Figure 2: Schematic of adaptations of the PNP model based on various cable
theory based models. Normalised profile of the action potential propagation
computed in this work using the corresponding PNP model is depicted for the
three variants: A) The PNP model comprises of ionic exchange between the
intra-cellular and the extra-cellular region through the ionic channels
present uniformly on the neuronal membrane. B) PNP with myelin model unifies
the PNP model and the presence of myelin sheath. The ion channels are only
present at the nodes of Ranvier. The low capacitance of the myelin sheath
steers the action potential with higher velocity. C) The presence of the
submyelin peri-axonal space is accounted for in the PNP model with myelin and
submyelin peri-axonal space.

##  4 Numerical implementation of the PNP model

The coupled nonlinear system of PDE’s for ionic concentration and electrical
potential in the electro-diffusive PNP model are solved using the standard
Finite Element Method (FEM). The primal fields that are solved for are the
voltage and the concentration of $Na^{+}$/ $K^{+}$/ $Cl^{-}$ ions. The
electric field and the net charge are derived fields. The salient features of
the computational implementation are: adaptive mesh refinement near the nodes
of Ranvier, adaptive time-stepping schemes, support for parallel direct and
iterative (Krylov-subspace) solvers with Jacobi/SOR preconditioning. The weak
formulation of the governing equations solved with FEM are given in Section
4.1. The computational framework is made available to the wider research
community as an open source library [37], and we hope it serves as a platform
for wider adoption of the high-fidelity PNP framework by neuronal modeling
researchers.

###  4.1 Weak formulation of the PNP model

The Nernst-Planck and Poisson equations, expressed in their weak (integral)
formulation that is suitable for the FEM implementation, following standard
notation, are as follows:

Find the primal fields $\\{V,c_{Na},c_{K},c_{Cl}\\}$, where,

$\displaystyle V$
$\displaystyle\in\mathscr{S}_{V},\quad\mathscr{S}_{V}=\\{V~{}|V=~{}\bar{V}~{}\forall~{}\textbf{X}\in\Gamma^{V}_{g}\\},$
$\displaystyle c_{i}$
$\displaystyle\in\mathscr{S}_{c_{i}},\quad\mathscr{S}_{c_{i}}=\\{c_{i}~{}|c_{i}=~{}\bar{c}_{i}~{}\forall~{}\textbf{X}\in\Gamma^{c_{i}}_{g}\\}$

and $i\in\\{Na^{+},K^{+},Cl^{-}\\}$, such that,

$\displaystyle\forall~{}w^{V}$
$\displaystyle\in\mathscr{V}_{V},\quad\mathscr{V}_{V}=\\{V~{}|V=~{}0~{}\forall~{}\textbf{X}\in\Gamma^{V}_{g}\\},$
$\displaystyle\forall~{}w^{c_{i}}$
$\displaystyle\in\mathscr{V}_{c_{i}},\quad\mathscr{V}_{c_{i}}=\\{c_{i}~{}|c_{i}=~{}0~{}\forall~{}\textbf{X}\in\Gamma^{c_{i}}_{g}\\}$

we have,

$-\frac{F}{\epsilon}\int_{\Omega}w^{V}~{}(c_{Na}+c_{K}-c_{Cl})~{}dV+\int_{\Omega}\nabla
w^{V}\cdot\nabla V~{}dV-\int_{\Gamma^{V}_{h}}w^{V}~{}(\nabla
V\cdot\bm{n})~{}dS=0$ (8)

and,

$\int_{\Omega}w^{c_{i}}~{}\frac{\partial c_{i}}{\partial
t}~{}dV+\int_{\Omega}\nabla w^{c_{i}}\cdot~{}D_{i}\nabla
c_{i}~{}dV+\int_{\Omega}\nabla
w^{c_{i}}\cdot~{}D_{i}\frac{c_{i}z_{i}F}{RT}\nabla
V~{}dV+\int_{\Gamma^{V}_{h}}w^{c_{i}}~{}(\bm{F}_{i}\cdot\bm{n})~{}dS=0$ (9)

where $w^{V}$ is the variation for the voltage field, $w^{c_{i}}$ are the
variations for the ionic fields. $\Omega$ is the problem geometry,
$\\{\Gamma^{V}_{g},\Gamma^{c_{i}}_{g}\\}$ are the Dirichlet boundaries and
$\\{\Gamma^{V}_{h},\Gamma^{c_{i}}_{h}\\}$ are the Neumann (Flux) boundaries of
the voltage and ionic fields, respectively. $\bm{n}$ is the unit normal vector
. In this work, there is no voltage flux ($\nabla V\cdot\bm{n}=0$) at all the
boundaries. Eq. 9 and Eq. 8 are the governing equations solved using FEM.

The PNP framework models the evolution of the voltage field and the ionic
concentrations of $Na^{+}$, $K^{+}$ and $Cl^{-}$. The initial ionic
concentrations in the various regions such as the extra-cellular region,
membrane, myelin, cytoskeleton, etc., are mentioned in the Supplemental
Information. The initial voltage in the ECM is taken to be 0 mV and in the
cytoskeleton region to be the resting value of -70 mV. The boundary conditions
on the top and bottom surface of the domain are applied so that the fields,
i.e. voltage and the ionic concentrations, represent their bulk value in the
extra-cellular region as in Eq. 10. The ionic exchange at the nodes of Ranvier
is incorporated as an ionic flux as given by Eq. 11. Here $I_{i}$ is the
current of each ionic concentration computed using their respective Hodgkin-
Huxley ionic conductance.

At top and bottom boundaries;

$c_{i}=c_{i}^{e},\qquad V=0$ (10)

At the Node of Ranvier:

$\bm{F}_{i}\cdot\bm{n}=f_{i},\qquad f_{i}=\frac{I_{i}}{z_{i}F}$ (11)

where,

$I_{Na}=\bar{G}_{Na}m^{3}h(V_{m}-E_{Na})$
$I_{K}=\bar{G}_{K}n^{4}(V_{m}-E_{K})$ $I_{Cl}=0$

###  4.2 Estimating effective electrical properties from ionic distributions

In order to compare the action potential propagation modelled using the cable
theory and the PNP model, the embedded electrical properties of the PNP model
should be the same as the input values of the cable theory based models. The
capacitance of the PNP model can be computed using the membrane thickness and
myelin thickness to obtain $C_{m}$ and $C_{my}$, respectively. The axial
resistance offered to the conduction can be computed using the ionic constants
as presented in Figure 3. Initially, the membrane potential is at its resting
value of $-70$ mV with higher concentration of sodium ions in the
extracellular region and a higher concentration of potassium ions in the
intra-cellular region. The model is first equilibrated for a few timesteps
until there are no fluctuations in the field variables. Next, to initiate the
action potential, we assume a sodium ion influx in a region of length $5\mu m$
near the left end of the neuronal membrane until the local potential reaches
the threshold potential. Thereafter, the activation parameters regulate the
depolarization/repolarization.

Figure 3: Electrical connection between the cable theory based model and the
Poisson-Nernst-Planck model. The spatial electro-diffusive PNP model comprises
of multiple fields including the voltage and the concentration of the
respective ions. The capacitance of the membrane and the myelin sheath are
embedded in the PNP model in the form of membrane/myelin thickness. The ion
channel currents are treated as a flux (Neumann boundary condition) in the
spatial model. The resistance offered to the action potential along the axon
can be computed using the ionic constants and the expression given in the
figure. The derived fields namely, the electric field and net charge are also
plotted. Note that the plot of electric field is zoomed in around the membrane
as there is no appreciable electric field in the ECM.

##  5 Results

This section seeks to establish the significance of the PNP numerical
framework by analyzing the results of the three variants of the electro-
diffusive model discussed in the preceding sections. First and foremost, an
unmyelinated PNP model for a squid axon is investigated and distinctions with
the cable theory based Hodgkin-Huxley model are probed. Subsequently, the
simulation results for the PNP with myelin and PNP with myelin plus peri-
axonal space model are carefully studied and a detailed comparison with the
corresponding cable theory based model is presented. Neuronal axons of a rat
and squid are considered for the numerical simulations. The parameters of the
respective axons are collected from [6, 10, 22, 21] and are listed in the
Supplemental Information.

###  5.1 PNP model

Figure 4 compares the action potential resulting from the cable theory based
Hodgkin-Huxley model to the PNP model. For this simulation, the neuronal axon
of a squid having a length of $20,000$ $\mu$m, a radius of $238$ $\mu$m, and a
rat neuron having a length of $1500$ $\mu$m, a diameter of $1.1$ $\mu$m are
considered. Note that a larger length for the squid axon is selected so that
we can clearly observe the propagation profile of the action potential in the
simulations. From the plots, we can clearly visualize that the one-dimensional
lineout of the spatial voltage resulting from the PNP model is similar to the
voltage from the cable theory based Hodgkin-Huxley model. The charge is
accumulated near the membrane due to the property of the membrane to act as a
capacitor. This leads to the formation of the Debye layer near to the
membrane. As reflected from the plot of net charge, the variation of ionic
concentration is appreciable near the membrane. The electric field is
prominent throughout the membrane due to the potential difference between the
intracellular and the extracellular region, however, the magnitude of this
electric field varies with the propagation of the action potential along the
axon. The initial equivalent resistance offered to the action potential
propagation along the squid axon is computed to be $38.18$ $\Omega$ cm using
the underlying ionic constants. The conduction velocity for the squid and the
rat axon is computed to be 15.38 m/s and 0.36 m/s respectively.

###  5.2 PNP model with myelin

The action potential and the respective ionic concentrations modeled using the
PNP with myelination model are presented in Figure 5. The one-dimensional
lineout of the primary fields is plotted on the right. We observe that the
action potential propagates as a soliton-like wave as in the single-cable
model. While the saltatory conduction observed in the single-cable model is
due to the potential jump from one node of Ranvier to the next, the saltatory
conduction in the electro-diffusive model is observed in the respective ionic
concentrations. As in the single-cable model, the saltatory conduction is due
to the presence of nodes of Ranvier, and therefore the myelin sheath. Low
capacitance of the myelin sheath which quickly charges and discharges, leads
to fast propagation of the charge. The envelope of the one-dimensional lineout
of the respective ionic concentration propagates with the action potential. We
also note that there is no correlation of the myelin voltage in the
electrodiffusive model as in the single cable model. The electric field,
prominent at the nodes of Ranvier and the membrane is also depicted. As in the
previous model, the accumulation of the net charge closer to the membrane can
be observed.

Here, we discuss the mechanism of action potential propagation using the PNP
electrodiffusive model. As compared to the bulk, the net charge is
considerable near the membrane due to the membrane capacitance. This leads to
lower net resistance in the intra-cellular region near the membrane. To
initiate the action potential, there is sodium ion influx from the first node
of Ranvier. The high electric field/ high voltage gradient at this node of
Ranvier steers the charge to propagate at a faster velocity along the membrane
as this offers a low resistance conduction path. These ions reach the adjacent
node of Ranvier, resulting in the membrane potential attaining the threshold
potential. At this point, the Hodgkin-Huxley activation/inactivation
parameters govern the ionic flux and the propagation of action potential along
the entire length of the axon. One must realize that the ionic diffusion is
immensely slow but the combination of high voltage gradient/electric field at
the respective nodes of Ranvier, the ionic flux through the membrane leads to
the local current and thus faster propagation of the action potential along
the axon. This ability of the PNP framework to model the electric field and
net charge, in addition to the action potential and the ionic distributions,
make it a high-fidelity and more comprehensive model compared to the cable
theory models.

###  5.3 PNP model with myelin and peri-axonal space

The propagation of the action potential using the PNP model with myelin and
peri-axonal space is depicted in Figure 6. The profile of the action potential
propagation is similar to the earlier presented PNP model. However, the one-
dimensional lineout plot of the respective ionic concentration depicts a lot
of surges due to the presence of the peri-axonal space in addition to the
nodes of Ranvier. The envelope of the sodium ion concentration propagates with
the action potential. However, the envelope of the potassium and the chloride
ion travel relatively fast with respect to the action potential. Again, we can
visualize the spatial saltatory conduction from the respective ionic
concentration which jumps from one node of Ranvier to the next. As in the
previous PNP model, we note the lower amplitude of action potential
propagation, the dominant electric field at the nodes of Ranvier, and the
accumulation of net charge near the membrane. Note that the magnitude of the
electric field varies with the neuronal conduction along the axon. Simulation
video of the myelinated PNP models showing the spatio-temporal evolution of
all the fields is provided as part of the Supplemental Information.

###  5.4 Spatial resolution of the primal fields

Based on the simulation results presented in the previous sections, one can
appreciate the potential of the PNP framework to augment the voltage imaging
experiments where the neuron is stained with a voltage-sensitive or ion-
sensitive dye to carry out an analysis of the action potential propagation and
its properties such as conduction velocity. Figure 7 demonstrates the
capability of the PNP framework to determine the voltage fields and the
respective ionic concentration to a fine resolution. This fine resolution is
important, as it’s the small variations in the respective ionic concentrations
that eventually result in the quintessential propagation of the action
potential. At any point P in the domain, one can visualize the temporal
variation of the respective primary or derived fields. Thus, in addition to
the voltage imaging experiments, this model can potentially act as a digital
twin for the patch-clamp/voltage-clamp experiments where the profile of the
voltage or the current is plotted at a specific point in the neuron where the
electrode is inserted.

Figure 4: Comparison of the cable theory based Hodgkin-Huxley to the PNP model
for a squid axon. The models have identical inputs of the membrane capacitance
and resistance. The schematics depict the capability of the PNP model to
visualize the spatial distribution of the voltage, ionic concentration, net
charge and the electric field. Net charge and electric field are dominant at
the membrane. The profile of the one dimensional lineout of the action
potential propagation from the PNP model resembles to the action potential
propagation modelled using the cable theory based classical Hodgkin-Huxley
model. Figure 5: Schematics of the various fields modelled using the PNP model
with myelin for a rat axon. A lineout of the one dimensional profile of the
primary fields extracted from the spatial PNP model is depicted in the right.
The action potential propagates as a soliton-like wave as in the HH model. The
ionic concentrations depict a saltatory conduction where it appears that the
concentration jumps from one node of Ranvier to another. This is unlike the
saltatory conduction observed using the one dimensional single-cable model.
This is observed for all the ionic species present. Net charge is accumulated
near the membrane. Electric field is dominant at the nodes of Ranvier. Figure
6: Schematics of the fields modelled using the PNP model with myelin and
periaxonal space for a rat axon. A lineout of the one dimensional profile of
the primary fields extracted from the spatial PNP model is depicted in the
right. The profile of the action potential is similar to earlier models. As in
the PNP with myelin model, the ionic concentrations depict a saltatory
conduction jumping from one node from Ranvier to the next. However, this
conduction is faster here as compared to the PNP with myelin model. Figure 7:
Spatial resolution of the primary and derived fields in the PNP with myelin
model of a rat axon. The framework of the PNP based models can provide insight
into the various spatial fields present at a material point on the neuronal
membrane or its vicinity. The time evolution of various field values at a
point ‘P’ in the intra-cellular region is depicted in the plots.

##  6 Discussion on estimated velocity of action potential propagation

The conduction velocity for a rat axon computed using the various cable theory
based models and the PNP models are compared in Figure 8a. For each of the
models based on cable theory or PNP, the conduction velocity increases due to
myelination or the presence of submyelin peri-axonal space. While the
conduction velocity computed using the PNP and the PNP with myelin are in
close proximity to that computed from the classical Hodgkin-Huxley and the
single-cable model respectively, the double cable model has nearly a nine fold
increase in the speed compared to the PNP model with myelin and periaxonal
space. This may be due to the fact that, while the double-cable model employs
the peri-axonal resistance and the para-nodal resistance, these values for the
PNP model depend on the underlying ionic constants and are therefore of the
same order of magnitude as in the axonal cytoplasm. It has been assumed that
the ionic diffusion coefficient is the same in the cytoskeleton of the
neuronal axon and the sub-myelin peri-axonal space which may not be the case.
However, once these values are obtained experimentally, we can easily
incorporate this in the PNP framework. The slight difference in the conduction
velocity computed using the PNP model, PNP with myelin model as compared to
the cable theory based Hodgkin-Huxley, single cable model is discussed next.
Initially the equivalent axial resistance of the PNP model is comparable to
the cable theory based model. With the propagation of action potential and
because of the nature of the membrane to act as a capacitor, the net charge is
appreciable near the membrane in the intracellular region. This leads to lower
net resistance near the membrane in the PNP models and the differences in the
conduction velocity when compared to their corresponding cable theory based
model.

For a squid, there is a two fold increase in the computed conduction velocity
when the electro-diffusive model includes myelination. The conduction speed of
the Hodgkin-Huxley model is close to the value computed using the PNP model.
The primary difference in neuronal conduction velocity of a squid and a rat is
due to the larger diameter of the squid axon. The net axial resistance in the
rat axon is therefore relatively higher.

(a)  (b)

Figure 8: Estimated conduction velocity (CV) of action potential propagation.
(a) Numerical estimates of the CV of the action potential propagating through
the rat neuron modeled using variants of cable theory and PNP based models.
Using the cable theory based models, the CV increases due to the presence of
the myelin sheath and the peri-axonal space. The same trend is observed for
the PNP models but the increase in the conduction speed when the peri-axonal
space is present in the PNP model, is comparatively lower than the double-
cable model. (b) CV for the action potential propagating through the squid
neuron modeled using cable theory and PNP models. The increase in the CV is
significant when myelin is present. Here, HH, SC, and DC abbreviations are
used for the Hodgkin-Huxley, Single-Cable, and Double-Cable models,
respectively, to distinguish between the three 1D electrical network models.

##  7 Conclusion

Cable theory based Hodgkin-Huxley, single-cable and double cable models have
been extensively used to model the action potential propagation along the
neuron. These models are essentially a one dimensional reduced order
representation of the complex action potential dynamics, and ignore the
spatial ionic diffusion which can be crucial to study propagation in non-
trivial geometries and in the heterogeneous neuronal microenvironment, such as
in action potential propagating through thin dimensions in dendrites, and of
multiple action potentials. The Poisson-Nernst-Planck model provides a robust
framework to represent the action potential propagation due to the spatio-
temporal variation of the ionic concentrations at a much finer resolution.
Also, the derived fields, such as the electric field and the net charge, can
be easily obtained and analyzed. Moreover, as shown in this work, this model
is able to emulate the various variants of the cable theory based models. This
framework can also be extended to model the heterogeneous microenvironment
around the neuronal membrane, such as the neuron-glia interactions.

To the best of our knowledge, the framework presented here for numerically
estimating the conduction velocity of the spatial action potential propagating
using the PNP model, is a first-of-its-kind. Further, the spatio-temporal
manifestation of the well-known saltatory conduction mechanism has been
demonstrated. This electro-diffusive PNP framework can augment experiments to
yield an elevated understanding of the action potential propagation, and thus
assisting us to get a better insight into the functioning of our
neurophysiology, and into the action potential related underpinnings of
various neurological diseases.

## Acknowledgement

The authors thank the U.S. Office of Naval Research (ONR) for supporting this
work under the “PANTHER” program (award number: N00014-21-1-2918) through Dr.
Timothy Bentley. The authors would also like to thank Prof. Christian Franck
at the University of Wisconsin-Madison and Prof. Ashfaq Adnan at the
University of Texas at Arlington for insightful discussions on the subject.

## References

## References

  * Marion et al. [2018] C. M. Marion, K. L. Radomski, N. P. Cramer, Z. Galdzicki, R. C. Armstrong, Experimental Traumatic Brain Injury Identifies Distinct Early and Late Phase Axonal Conduction Deficits of White Matter Pathophysiology, and Reveals Intervening Recovery, The Journal of Neuroscience 38(41) (2018) 8723–8736.
  * Palop et al. [2007] J. J. Palop, J. Chin, E. D. Roberson, J. Wang, M. T. Thwin, N. Bien-Ly, J. Yoo, K. O. Ho, G.-Q. Yu, A. Kreitzer, S. Finkbeiner, J. L. Noebels, L. Mucke, Aberrant Excitatory Neuronal Activity and Compensatory Remodeling of Inhibitory Hippocampal Circuits in Mouse Models of Alzheimer’s Disease, Neuron 55 (2007) 697–711.
  * Palop and Mucke [2010] J. J. Palop, L. Mucke, Amyloid-$\beta$ Induced Neuronal Dysfunction in Alzheimer’s Disease: From Synapses toward Neural Networks, Nature Neuroscience 13 (2010) 812–818.
  * Ghatak et al. [2019a] S. Ghatak, N. Dolatabadi, D. Trudler, X. Zhang, Y. Wu, M. Mohata, R. Ambasudhan, M. Talantova, S. A. Lipton, Mechanisms of hyperexcitability in Alzheimer’s disease hiPSC-derived neurons and cerebral organoids vs isogenic controls, eLife 8 (2019a) e50333.
  * Chaudhury et al. [2015] D. Chaudhury, H. Liu, M.-H. Han, Neuronal Correlates of Depression, Cellular and molecular life sciences 72(24) (2015) 4825–4848.
  * Hodgkin and Huxley [1952] A. Hodgkin, A. Huxley, A quantitative description of membrane current and its application to conduction and excitation in nerve, The Journal of Physiology 117(4) (1952) 500–544.
  * Huxley [1959a] A. F. Huxley, Ion movements during nerve activity, Annals New York Academy of Sciences 81 (1959a) 221–246.
  * Duncan and Radcliff [2016] I. D. Duncan, A. B. Radcliff, Inherited and acquired disorders of myelin: The underlying myelin pathology, Experimental Neurology 283 (2016) 452–475.
  * Zhan et al. [2014] X. Zhan, G. Jickling, B. Ander, D. Liu, B. Stamova, C. Cox, L. Jin, C. DeCarli, F. Sharp, Myelin injury and degraded myelin vesicles in Alzheimer’s disease, Current Alzheimer research 11(3) (2014) 232–238.
  * Cohen et al. [2020] C. C. Cohen, M. A. Popovic, J. Klooster, M.-T. Weil, W. Möbius, K.-A. Nave, M. H. Kole, Saltatory conduction along myelinated axons involves a periaxonal nanocircuit, Cell 180 (2020) 311–322.
  * Richardson et al. [2000] A. Richardson, C. McIntyre, W. Grill, Modelling the effects of electric fields on nerve fibres: influence of the myelin sheath, Medical and Biological Engineering and Computing 38 (2000) 438–446.
  * Barett and Barrett [1982] E. F. Barett, J. N. Barrett, Intracellular recording from vertebrate myelinated axons: mechanism of the depolarizing afterpotential, Journal of Physiology 323 (1982) 117–144.
  * Blight [1985] A. Blight, Computer simulation of action potentials and afterpotentials in mammalian myelinated axons: the case for a lower resistance myelin sheath, Neuroscience 15 (1985) 13–31.
  * Blight and Someya [1985] A. Blight, S. Someya, Depolarizing afterpotentials in myelinated axons of mammalian spinal cord, Neuroscience 15 (1985) 1–12.
  * Pods et al. [2013] J. Pods, J. Schonke, P. Bastian, Electrodiffusion Models of Neurons and Extracellular Space Using the Poisson-Nernst-Planck Equations—Numerical Simulation of the Intra- and Extracellular Potential for an Axon Model, Biophysical Journal 105 (2013) 242–254.
  * Mori [2009] Y. Mori, From Three-Dimensional Electrophysiology to the Cable Model: an Asymptotic Study, arXiv preprint 0901. (2009) 3914.
  * Qian and Sejnowski [1989] N. Qian, T. Sejnowski, An Electro-Diffusion Model for Computing Membrane Potentials and Ionic Concentrations in Branching Dendrites, Spines and Axons, Biological Cybernetics 62 (1989) 1–15.
  * Halnes et al. [2013] G. Halnes, I. Ostby, K. H. Pettersen, S. W. Omholt, G. T. Einevoll, Electrodiffusive Model for Astrocytic and Neuronal Ion Concentration Dynamics, PLoS Computational Biology 9(12) (2013) e1003386.
  * Mori and Peskin [2009] Y. Mori, C. S. Peskin, A Numerical method for cellular electrophysiology based on the electrodiffusion equations with internal boundary conditions at membranes, Communications in Applied Mathematics and Computational Science 4 (2009) 85–134.
  * Halnes et al. [2016] G. Halnes, T. Maki-Marttunen, D. Keller, K. H. Pettersen, O. A. Andreassen, G. T. Einevoll, Effect of Ionic Diffusion on Extracellular Potentials in Neural Tissue, PLOS Computational Biology 12(11) (2016) e1005193.
  * Lopreore et al. [2008] C. L. Lopreore, T. M. Bartol, J. S. Coggan, D. X. Keller, G. E. Sosinsky, M. H. Ellisman, T. J. Sejnowski, Computational Modeling of Three-Dimensional Electrodiffusion in Biological Systems: Application to the Node of Ranvier, Biophysical Journal 95 (2008) 2624–2635.
  * Dione et al. [2016] I. Dione, J. Deteix, T. Briffard, E. Chamberland, Improved Simulation of Electrodiffusion in the Node of Ranvier by Mesh Adaptation, PLOS ONE 11(8) (2016) e0161318.
  * Goriely et al. [2015] A. Goriely, M. G. Geers, G. A. Holzapfel, J. Jayamohan, A. J$\acute{e}$rusalem, S. Sivaloganathan, W. Squier, J. A. W. v. Dommelen, S. Waters, E. Kuhl, Mechanics of the brain: perspectives, challenges, and opportunities, Biomechanics and Modeling in Mechanobiology 14 (2015) 931–965.
  * Garcia-Gonzalez and Jerusalem [2019] D. Garcia-Gonzalez, A. Jerusalem, Energy based mechano-electrophysiological model of CNS damage at the tissue scale, Journal of the Mechanics and Physics of Solids 125 (2019) 22–37.
  * Kwong et al. [2019] M. T. Kwong, F. Bianchi, M. Malboubi, J. A. Garc$\acute{i}$a-Grajales, L. Homsi, M. Thompson, H. Ye, L. Noels, A. J$\acute{e}$rusalem, 3D finite element formulation for mechanical–electrophysiological coupling in axonopathy, Computer Methods in Applied Mechanics and Engineering 346 (2019) 1025–1050.
  * Auddya et al. [2021] D. Auddya, X. Zhang, R. Gulati, R. Vasan, K. Garikipati, P. Rangamani, S. Rudraraju, Biomembranes undergo complex, non-axisymmetric deformations governed by Kirchhoff-Love kinematics and revealed by a three dimensional computational framework, Proceedings of the Royal Society A 477(2255) (2021) 20210246.
  * Huxley [1959b] A. F. Huxley, Ion movements during nerve activity, Second conference on physicochemical mechanism of nerve activity and second conference on muscular contraction 81(2) (1959b) 221–246.
  * Koch [2004] C. Koch, Biophysics of Computation: Information Processing in Single Neurons, vol. 1, Oxford University Press, 2004.
  * Stephanova and Bostock [1995] D. Stephanova, H. Bostock, A distributed-parameter model of the myelinated human motor nerve fibre: temporal and spatial distributions of action potentials and ionic currents, Biological Cybernetics 73 (1995) 275–280.
  * Bar-Kochba et al. [2016] E. Bar-Kochba, M. T. Scimone, J. B. Estrada, C. Franck, Strain and rate-dependent neuronal injury in a 3D in vitro compression model of traumatic brain injury, Scientific Reports 6 (1) (2016) 1–11.
  * Li et al. [2011] W. Li, Q. Fan, Z. Ji, X. Qiu, Z. Li, The effects of irreversible electroporation (IRE) on nerves, PloS one 6 (4) (2011) e18831.
  * Estrada et al. [2021] J. B. Estrada, H. C. Cramer III, M. T. Scimone, S. Buyukozturk, C. Franck, Neural cell injury pathology due to high-rate mechanical loading, Brain Multiphysics 2 (2021) 100034.
  * Palop and Mucke [2016] J. J. Palop, L. Mucke, Network abnormalities and interneuron dysfunction in Alzheimer disease, Nature Reviews Neuroscience 17 (12) (2016) 777–792.
  * Ghatak et al. [2019b] S. Ghatak, N. Dolatabadi, D. Trudler, X. Zhang, Y. Wu, M. Mohata, R. Ambasudhan, M. Talantova, S. A. Lipton, Mechanisms of hyperexcitability in Alzheimer’s disease hiPSC-derived neurons and cerebral organoids vs isogenic controls, Elife 8 (2019b) e50333.
  * Kim et al. [2007] D. Y. Kim, B. W. Carey, H. Wang, L. A. Ingano, A. M. Binshtok, M. H. Wertz, W. H. Pettingell, P. He, V. M.-Y. Lee, C. J. Woolf, et al., BACE1 regulates voltage-gated sodium channels and neuronal activity, Nature cell biology 9 (7) (2007) 755–764.
  * Freidberg [2008] J. P. Freidberg, Plasma Physics and Fusion Energy, vol. 1, Cambridge University Press, 2008.
  * Git [2022] Code repository for Action potential propagation along the neuronal axon and spatio-temporal saltatory conduction, https://github.com/cmmg/neuronalActionPotential, 2022\.


# Nematic Spin Liquid State in the $S=\frac{5}{2}$ Heisenberg Kagomé
Antiferromagnet Li9Fe3(P2O7)3(PO4)2

E. Kermarrec  R. Kumar  R. Hénaff  Université Paris-Saclay, CNRS, Laboratoire
de Physique des Solides, 91405, Orsay, France  B. Koteswararao  Faculty in
Physics, Indian Institute of Technology Tirupati  P. L. Paulose  Department of
Condensed Matter Physics and Materials Science, TIFR, Mumbai  P. Mendels  F.
Bert  Université Paris-Saclay, CNRS, Laboratoire de Physique des Solides,
91405, Orsay, France

###### Abstract

We investigate the low temperature magnetic properties of a $S=\frac{5}{2}$
Heisenberg kagomé antiferromagnet, the layered monodiphosphate
Li9Fe3(P2O7)3(PO4)2, using magnetization measurements and 31P nuclear magnetic
resonance. An antiferromagnetic-type order sets in at $T_{\rm N}=1.3$ K and a
characteristic magnetization plateau is observed at 1/3 of the saturation
magnetization below $T^{*}\sim 5$ K. A moderate 31P NMR line broadening
reveals the development of anisotropic short-range correlations within the
plateau phase concomitantly with a gapless spin-lattice relaxation time
$T_{1}\sim k_{B}T/\hbar S$, which both point to the presence of a
semiclassical nematic spin liquid state predicted for the Heisenberg kagomé
antiferromagnetic model.

###### pacs:

75.10.Jm, 76.60.-k, 75.30.Cr, 75.10.Kt, 75.40.Gb,

Magnetic frustation sets as a playground for the discovery of new exotic
phases of quantum matter. The Heisenberg antiferromagnetic model on the kagomé
lattice (KHAF) constitutes the archetype of the frustration in two dimensions,
whose ground-state and exotic spin dynamics still remain topical questions in
quantum magnetism even after 30 years of research [1, 2, 3, 4].

Early on, Chalker et al.[5] highlighted the distinctive magnetic properties of
the classical ($S=\infty$) KHAF, for which thermal fluctuations select a spin
nematic, or coplanar, ground state via the order-by-disorder mechanism in zero
field [6, 7], with spins on a single triangle oriented at 120∘ from each other
(Fig.1e). Its spin dynamics is governed by zero-energy soft modes, i.e.
special lines of weathervane defects [8]. The order-by-disorder mechanism is
also responsible for the remarkable existence of semiclassical fractional
plateaus of magnetization in frustrated magnets [1, 9, 10, 11, 12, 13]. In
particular, the 1/3 plateau –a hallmark of frustration[1, 14]–, arising from
collinear uud spin configurations (Fig.1d), has been observed in the
triangular compounds Cs2CuBr4[15] and RbFe(MoO4)2[16, 17, 18] but remains
debated and/or difficult to observe for quantum kagomé antiferromagnets with
large interaction [21] or with deviations from the pure KHAF model [19, 20].

More recently, the phase diagram of the KHAF in an external field was
revisited with state-of-the art numerical methods, notably tensor network
formalism[22, 23], DMRG[13] or exact diagonalization[12]. It encompasses fully
“quantum” phases and semiclassical plateaus, including the 1/3 plateau that
appears for all values of spin $S$. The specific aspect of the 1/3 plateau of
the classical kagomé lattice under an applied field was pinpointed by finite-
temperature Monte-Carlo studies which showed that thermal fluctuations favor a
a classical collinear spin liquid phase [24, 25] at $M=\frac{1}{3}M_{\rm
sat}$, made of a macroscopic number of degenerate uud states.

The comparison with classical experimental candidates has remained challenging
for long as naturally occuring kagomé lattice are rare. The jarosites family
AFe3(OH)6(SO4)2 (with A+ typically an univalent cation) is held up as a
paragon of the classical KHAF[1, 26, 27]. They all (but A+=D3O+[28]) display
long range ordering into a $q=0$ structure, as a result of a deficient kagomé
network and/or of perturbative anisotropic terms. Furthermore, their large
antiferromagnetic interaction $J\sim 600$ K requires unattainable applied
magnetic fields in order to reach the 1/3 plateau.

Figure 1:  (Color online) (a) Schematic structure of LiFePO illustrating the
corner-sharing triangular lattice made of Fe3+, $S=\frac{5}{2}$ (brown). Li
and O atoms are omitted for clarity. (b), (c) Fe atoms are magnetically
coupled through PO4 chemical groups with two different phosphorus sites: P1
(yellow) and P2 (blue). Illustration of the up-up-down uud (d) and of the 120∘
ordered structure (e).

In this Letter, we investigate the low temperature magnetic properties of the
layered iron monodiphosphate Li9Fe3(P2O7)3(PO4)2 (LiFePO), a $S=\frac{5}{2}$
Heisenberg kagomé antiferromagnet material with a moderate exchange $J\sim 1$
K. LiFePO thus offers the opportunity to probe and finely explore
experimentally the phase diagram of the classical KHAF under applied field.
Our measurements reveal a clear 1/3 magnetization plateau with an extended
spin-liquid regime characterized by exotic spin dynamics.

LiFePO belongs to the family of monodiphosphates, which were primarily studied
for their electronic transport properties in the context of battery materials
[29]. LiFePO crystallizes in the centrosymmetric trigonal space group
P$\overline{3}c1$ where Fe3+ ($S=\frac{5}{2}$) atoms form corner sharing
triangles, topologically equivalent to a regular kagomé lattice (Fig. 1), with
a strong Heisenberg character ($L\simeq 0$).

The macroscopic magnetic susceptibility, measured by a DC SQUID magnetometer,
is shown as a function of temperature in Fig. 2a. Between 100 and 300K the
susceptibility obeys a Curie-Weiss law, $C/(T-\theta_{\rm CW})$, with a Curie-
Weiss temperature $\theta_{\rm CW}=-11(2)$ K, indicative of antiferromagnetic
interactions. For a simple nearest-neighbour Heisenberg model on a kagomé
lattice, this leads to an exchange interaction $J\sim 1$ K. The extracted
value of the Curie constant $C=N_{A}\mu_{0}\mu_{\rm eff}^{2}/3k_{B}$ yields an
effective moment $\mu_{\rm eff}=5.85(2)$ $\mu_{B}$, in very good agreement
with the $5.92$ $\mu_{B}$ value expected for a fully occupied kagomé lattice
made of Fe3+ ions. At $T=T_{\rm N}=1.3$ K a sharp peak in the magnetic
susceptibility signals the onset of an antiferromagnetic order (Fig.2a). The
transition temperature $T_{\rm N}$ is reduced as compared to $|\theta_{\rm
CW}|=11$ K because of the frustration and the bidimensionality inherent to the
kagomé geometry. Li9Fe3(P2O7)3(PO4)2 thus shows a correlated paramagnetic
regime between 1.3 and $\sim 10$ K.

Magnetization versus field M($B$) curve is shown in Fig. 2b. Above 10 K,
$S=5/2$ spins are in a paramagnetic regime and M($B$) evolves almost linearly
from 0 to 16 T. As $T$ is decreased, a pronounced curvature is observed
towards the saturated value $M_{\rm sat}=5$ $\mu_{B}$. For $T<T^{*}\simeq 5$ K
the slope of the magnetization clearly diminishes before increasing again, in
the vicinity of $\sim M_{\rm sat}/3$. Fig. 2c shows the derivative $dM/dH$
which firmly confirms the existence of a narrow 1/3 magnetization plateau
developing below $T^{*}$. The evolution of the 1/3 plateau is quantified by
the size of the dip, $\Delta(dM/dH)$, and is shown as a function of $T$ in
Fig. 2.d. It displays a smooth crossover from a paramagnet for $T>T^{*}$ to
the plateau phase, i.e. a classical collinear spin liquid phase formed by the
set of the $uud$ states that remain degenerate at $M=M_{\rm sat}/3$, with no
clear sign of symmetry-breaking like for a gas-liquid type transition. In the
$T\rightarrow 0$ limit, its width can be estimated theoretically from Ref.[24]
to $\sim 1.5$ T in agreement with our experimental value of $\sim 1.6$ T
(FWHM). The small $dM/dH$ singularity observed around 0.7 T at 0.5 K is close
to a 1/9 fraction but is perhaps linked to a field-induced magnetic transition
as found in RbFe(MoO4)2 with an incommensurate to commensurate transition
along the $c$ axis[18].

The most recent theoretical phase diagram under applied fields has been
computed using numerical tensor network methods. It shows three different
incompressible phases or magnetization plateaus besides the 1/3 plateau[23].
The $1/9$-plateau is only expected for the most quantum case
($S=\frac{1}{2}$). The others (one-magnon and two-magnon plateaus [23]) should
appear for all values of $S$ and lie very close to the saturation,
respectively at 43/45 and 41/45 for $S=\frac{5}{2}$. These plateaus remain
absent in our 1.8 K data but a clear conclusion is difficult and would require
applying magnetic fields higher than 16T. The absence of the $C_{6}$
rotational symmetry of our $S=\frac{5}{2}$ Fe3+ lattice may also be
detrimental to these semiclassical phases.

Figure 2:  (Color online) Magnetization and susceptibility measurements of
LiFePO. (a) Magnetic susceptibility $\chi$ measured under 100 Oe showing an
antiferromagnetic transition at $T_{\rm N}=1.3$ K and Curie-Weiss fit (red
line) of $1/\chi$ (inset). (b) M vs $B$ measured at $T=1.8$ K shows a
saturated value $M_{\rm sat}\sim 5$ $\mu_{B}$/Fe and a 1/3-plateau for $B\sim
5$ T. (c) Derivative $dM/dH$ as a function of temperature (with a constant
offset for clarity). The dip at $\sim 5$ T signals the onset of the 1/3
plateau for $T=T^{*}$. (d) Evolution of the dip of the derivative of the
magnetization, $\Delta(dM/dH)$, with temperature (orange squares). Line is
guide to the eye. The classical spin liquid phase emerges below $T^{*}\simeq
5$ K.

In order to gain more insights into the magnetic properties of LiFePO, we
performed 31P NMR measurements. Li9Fe3(P2O7)3(PO4)2 naturally contains two
very sensitive NMR nuclei: 31P ($I=1/2$, $\gamma/2\pi=17.23467$ MHz/T) and 7Li
($I=3/2$, $\gamma/2\pi=16.54607$ MHz/T). 7Li was found to be mainly coupled to
Fe3+ through dipolar magnetic field resulting in broad lines which yield
relatively poor accuracy on the determination of the NMR shift. The NMR
spectra were thus obtained with $\nu_{0}=\gamma B_{0}/2\pi=109.732$ MHz for
the 31P local probe with $I=1/2$, avoiding any quadrupolar effects on the NMR
resonance line, using the standard $\pi/2-\tau-\pi$ pulse sequence, with
$\tau=20$ $\mu$s.

Fig. 3 displays the evolution of the 31P NMR spectrum with temperature. Two
NMR lines, accounting for the two 31P sites (labelled as 1 and 2) expected
from the structure [29] (Fig.1a,b), are clearly visible below 10 K. The NMR
intensity relative to site (1) gradually decreases as the temperature is
increased, and finally escapes our NMR time window above 50 K, due to a fast
relaxation, i.e. a shortening of the $T_{2}$ transverse relaxation time
($T_{2}<5-10$ $\mu$s). The 31P NMR Knight shift $\textbf{K}=\boldsymbol{A_{\rm
hf}\chi}/\left(g\mu_{B}\hbar\gamma\right)$ provides a direct measurement of
the local magnetic susceptibility tensor $\boldsymbol{\chi}$. Following the
common convention[30] in spherical coordinates, the eigenvalues of K relate to
the shift of the NMR line $\Delta B=B_{0}-B$ in the applied field $B$ through
[31]:

$\frac{\Delta B}{B}=K_{\rm iso}+K_{\rm
ax}\left(3\cos^{2}\theta-1\right)+K_{\rm ani}\sin^{2}\theta\cos 2\phi\\\ $ (1)

The spectra for sites (1) and (2) could be well fitted using Eq.1 for an
isotropic powder distribution (Fig. 3). The mismatch between observed
(circles) and simulated (red line) intensities is due to an anisotropic
$T_{2}$ spin-spin relaxation time, as confirmed by measurements using
different duration $\tau$ between the pulses (not shown). While site (1) and
site (2) probe the same local susceptibility, the anisotropic components of
the Knight shift tensor are best revealed from site (2) and we thus focus on
that NMR line in the following. We then extract the temperature dependence of
the components of the Knight shift tensor for site (2), as shown in Fig. 4.

In the paramagnetic regime ($300-10$ K), the NMR lineshape of site (2) is
anisotropic and shows a pronounced shoulder on its right-hand side. However,
above 10 K, $K_{\rm iso}$ and $K_{\rm ani}$ are found to scale with each other
(Fig. 4) and thus the spectrum could be well fitted using a unique
$T$-dependent parameter –an isotropic magnetic susceptibility $\chi(T)$–, with
$K_{\rm iso}=A_{\rm hf}^{\rm iso}\chi(T)$, $K_{\rm ani}=A_{\rm hf}^{\rm
ani}\chi(T)$ and $K_{\rm ax}\simeq 0$. The $T$-independent components of the
hyperfine tensor per Fe atom are $A_{\rm hf}^{\rm iso}=2.21(3)$ kOe$/\mu_{B}$
and $A_{\rm hf}^{\rm ani}=0.591(5)$ kOe$/\mu_{B}$ for site (2). An anisotropic
hyperfine coupling tensor is a usual feature of 31P NMR.

Upon cooling below $T^{*}=5$ K the NMR lineshape evolves and reveals the
development of a strongly temperature-dependent anisotropy, which can thus
only be ascribed to the magnetic susceptibility. The decrease of $K_{\rm
iso,\rm 2}$ signals the development of antiferromagnetic correlations when
$T<\theta_{\rm CW}=-11$ K, and clearly indicates the onset of a correlated
magnetic phase. Both site (1) and site (2) reveal an anisotropic
susceptibility as $K_{\rm ani}$ and $K_{\rm ax}$ undergo an abrupt change, a
remarkable feature for Heisenberg spins. The correlated magnetic phase
observed for $B\sim 6$ T and $1.3<T<10$ K, for which $M\sim\frac{1}{3}M_{\rm
sat}$, seems thus consistent with the theoretically expected, strongly
anisotropic, collinear phase.

Figure 3:  (Color online) Temperature dependence of the 31P NMR spectra
(circles) measured as a function of applied magnetic field $B$ for LiFePO,
from 80 K to 1.3 K. The phosphorus located on site (1) is only observed at low
$T$ (magnified at 40 K ($\times 3$)), due to strong $T_{2}$ relaxation
effects. Red line is a fit to Eq.1 (see text) from which the tensor components
$K$ are extracted. Vertical dashed line indicates the reference field
$B_{0}=6.3669$ T.

We now use 31P NMR to investigate the magnetic fluctuations in this
1/3-plateau phase. The spin-lattice relaxation rate $1/T_{1}$ probes the
imaginary part of the dynamic susceptibility
$\chi^{\prime\prime}_{\perp}(\nu,q)$ at low energy ($\nu\sim 100$ MHz)
through:

$\frac{1}{T_{1}}=\frac{\gamma^{2}}{g^{2}\mu_{B}^{2}}k_{B}T\sum_{q}\left|A_{\rm
hf}(q)\right|^{2}\frac{\chi^{\prime\prime}_{\perp}(\nu,q)}{2\pi\nu}\\\ $ (2)

and was determined through the saturation recovery method, with a
$\pi/2-\Delta t-\pi/2-\tau-\pi$ pulse sequence, from 1.3 to 300 K, at the
maximal intensity of site (2) (Fig. 5). Four different regimes in temperature
can be identified. A peak is clearly observed for $T\sim 1.5$ K, close to the
transition temperature $T_{\rm N}$, which separates the quasi-static magnetic
order from a classical spin-liquid state for $T_{\rm N}<T<T^{*}$ for which
$T_{1}\sim T$ (red line in Fig. 5). A crossover is observed at $T=T^{*}\sim 5$
K, and the dynamics is now dominated by a higher spin-lattice relaxation rate
for $T^{*}<T<80$ K, with an atypical power law behavior $\alpha=0.35(3)$ (blue
line) which persists well above $|\theta_{\rm CW}|$, up to $\sim 80$ K. We
speculate that lattice vibrational modes[44] or dipole-dipole interactions are
responsible for the unusual power law behavior of the spin dynamics observed
here. Finally, $T_{1}$ levels off in the high-$T$ paramagnetic regime. The
Moriya paramagnetic limit gives the $T$-independent rate
$1/T_{1}^{\infty}=\gamma^{2}g^{2}A_{\rm
hf}^{2}S(S+1)/3z_{1}\sqrt{2\pi}\nu_{e}$ [42, 43] in magnetic insulators, where
$\nu_{e}=J\sqrt{2zS(S+1)/3}/h$ the exchange frequency, with $z_{1}=2$ the
number of probed Fe and $z=4$ the number of magnetic nearest-neighbours of a
Fe atom. Using $A_{\rm hf}=A_{\rm hf}^{\rm iso}$ and the experimental value
$1/T_{1}^{\infty}=6$ ms-1, the formula leads to $J\simeq 2.2$ K, in reasonable
agreement with $J\sim 1$ K. We note that lithium diffusion likely starts to
contribute to the relaxation above 250 K.

Figure 4:  (Color online) (a) Components of the NMR shift K as a function of
temperature for site 2. (b) In the paramagnetic regime, $K_{\rm iso}$ and
$K_{\rm ani}$ scale with each other showing that the magnetic susceptibility
is isotropic above $T^{*}$. The red vertical dashed line indicates $T^{*}$.

We now discuss more specifically the classical spin liquid phase observed in
the $T$-range $1.5-5$ K in the context of theoretical predictions. The low$-T$
dynamics of the classical Heisenberg kagomé antiferromagnet has been
investigated early on by Monte-Carlo simulations [5, 33, 32]. They gave
evidence for the role of thermal fluctuations in selecting a coplanar short-
range order via the order-by-disorder mechanism for $T\rightarrow 0$. Short-
range nematic spin correlations are predicted to survive at low $T$ with non-
dispersive soft magnon modes[34], a characteristic feature associated with the
large degeneracy inherent to the kagomé lattice. Numerical simulations of the
spin autocorrelation function yield
$\left\langle\bm{S}(0)\bm{S}(t)\right\rangle\sim e^{-\nu t}$ for $H=0$ [1, 33,
35, 36], with $\nu=ck_{B}T/\hbar S$ ($c\sim 1$)[1, 37]. The kagomé bilayer
compound SrCr8Ga4O19 was the first experimental candidate to hint at such spin
dynamics, yet with an important contribution from quantum fluctuations [38,
33, 39, 40] and spin vacancies [41]. Later on, inelastic neutron scattering
exposed a $T$-linear inverse relaxation rate in deuterium jarosite over the
extended temperature range $0.06<T/J_{\rm cl}<1$ [28], with $J_{\rm
cl}=JS(S+1)$. Here, the $T$-dependence of $1/T_{1}$ can be inferred from $\nu$
through the Redfield formula:

$\frac{1}{T_{1}}=\frac{2\gamma^{2}\Delta^{2}\nu}{\nu^{2}+\gamma^{2}B^{2}}\simeq\frac{2\gamma^{2}\Delta^{2}}{\nu}\\\
$ (3)

where $\Delta$ is the amplitude of the fluctuating field at the phosphorus
position. Using the constant high-$T$ value $T_{1}^{\infty}$ set by the
Moriya’s exchange frequency $\nu_{e}$[42], this finally leads to $T_{1}=\beta
T$ with $\beta=ck_{B}T_{1}^{\infty}/\nu_{e}\hbar S$. In LiFePO, the dynamical
behavior of $T_{1}$ could be well fitted to such a law, with no adjustable
parameter but $c=1.06(3)$ (red line in Fig. 5), giving firm evidence of the
spin liquid phase within the $T$-range $0.2<T/J_{\rm cl}<0.6$ with $J_{\rm
cl}=8.25$ K.

In a more recent work[36], Taillefumier et al. characterized the dynamics for
the long time scales appropriate for NMR using semiclassical numerical
techniques, and predicted two different low temperature regimes in the
$B\rightarrow 0$ limit: (i) a cooperative, classical, spin-liquid regime with
$T_{1}\sim T$ in line with what is discussed above and (ii) a coplanar
dynamical regime for $T/J_{\rm cl}<0.005$ with a strong anisotropic relaxation
dominated by weathervane defects modes. Since $J_{\rm cl}=8.25$ K, our
measurements down to 1.3 K in LiFePO give only access to the region $T/J_{\rm
cl}>0.2$, i.e. the high-$T$ paramagnetic and cooperative spin liquid phases,
before residual interactions apparently lift the degeneracy and produce the
antiferromagnetic transition. The observed thermal behavior of our spin-
lattice relaxation rate, $T_{1}\sim T$, seems to show that this dynamics
unexpectedly persists under moderate applied fields ($B\sim 6$ T). This can
perhaps be understood because similar zero-energy excitations exist within the
$uud$ manifold of the 1/3 plateau. Furthermore, the short-range anisotropic
correlations that develop above $T_{N}$ according to our NMR local
susceptibility measurements may indicate either that the nematic correlations
are still present for $T_{N}<T<T^{*}$, or that we lie close to a quasi-static
magnetic order. In the light of our results, a theoretical study of the spin
dynamics under applied fields, in the 1/3-phase, would be highly relevant to
confirm the origin of the dynamics observed in NMR.

Figure 5:  (Color online) $1/T_{1}$ (up) and $T_{1}$ (bottom) probed by 31P
NMR as a function of temperature in LiFePO. Four different regimes can be
identified: a quasi-static order (QSO, grey), a classical spin liquid phase
(orange), an atypical high $T$ regime (green) and high-$T$ paramagnetic phase.
The classical spin liquid phase displays spin dynamics with a linear law
$T_{1}=\beta T$ (red line, see text).

In conclusion we have identified a new classical $S=\frac{5}{2}$ Heisenberg
kagomé antiferromagnet that displays the 1/3 magnetization plateau predicted a
long time ago. Our 31P NMR measurements further identified a classical spin-
liquid regime for $T_{N}<T<T^{*}$, for which the time scale of the spin
dynamics is set by the temperature only, with $T_{1}\sim T$, in very good
agreement with numerical estimates. These fluctuations could naturally
originate from the $uud$ spin configurations manifold on each triangle that
define the collinear spin liquid of the 1/3 plateau. The correlated low-$T$
regime of LiFePO combines thermal fluctuations with pronounced anisotropic
spin correlations and hence can be coined as a semiclassical nematic spin
liquid. The discovery of archetypal magnetic frustration relevant to the
kagomé lattice in the layered iron monodisphophate LiFePO offers great
perspectives to improve our knowledge on classical and quantum fluctuations
effects in zero and applied fields in the KHAF. In LiFePO, further elastic and
inelastic neutron scattering studies could help revisiting the central
question of the selection of the $q=0$ or $q=\sqrt{3}\times\sqrt{3}$ ground-
state in such model and also confirming the relevance of a simple, nearest-
neighbour exchange Hamiltonian, a natural question for this low symmetry
system. LiFePO and other members of this family [Cr3+ ($S=\frac{3}{2}$), V3+
($S=1$)] will certainly open new avenues to search for quantum phases relevant
to the KHAF.

## Acknowledgments

It is a pleasure to acknowledge useful discussion with Prof. S. Greenbaum on
lithium diffusion, P. Berthet and C. Decorse for sharing their expertise on
samples synthesis of lithium oxides. EK and RS acknowledge financial support
from the labex PALM for the QuantumPyroMan project (ANR-10-LABX-0039-PALM).
This work was supported by the French Agence Nationale de la Recherche under
grant ANR-18-CE30-0022-04 ‘LINK’.

## References

  * [1] Introduction to Frustrated Magnetism, edited by C. Lacroix, P. Mendels, and F. Mila (Springer, Berlin, 2011).
  * [2] M. R. Norman, Rev. Mod. Phys. 88, 041002 (2016).
  * [3] C. Broholm,R. J. Cava, S. A. Kivelson, D. G. Nocera, M. R. Norman, T. Senthil, Science 367, 6475 (2020).
  * [4] P. Khuntia, M. Velazquez, Q. Barthélemy, F. Bert, E. Kermarrec, A. Legros, B. Bernu, L. Messio, A. Zorko and P. Mendels, Nat. Phys. 16, 469-474 (2020).
  * [5] J. T. Chalker, P. C. W. Holdsworth and E. F. Shender, Phys. Rev. Lett. 68, 855 (1992).
  * [6] J. Villain, R. Bidaux, J.-P. Carton and R. Conte, J. Phys. (Paris) 41, 1263 (1980)
  * [7] C. L. Henley, Phys. Rev. Lett. 62, 2056 (1989).
  * [8] I. Ritchey, P. Chandra and P. Coleman, Phys. Rev. B 47, 15 342 (1993).
  * [9] A. Honecker, J. Schulenburg and J. Richter, J. Phys.: Condens. Matter 16 S749 (2004) and references therein.
  * [10] H. Kageyama, K. Yoshimura, R. Stern, N. V. Mushnikov, K. Onizuka, M. Kato, K. Kosuge, C. P. Slichter, T. Goto, and Y. Ueda, Phys.Rev.Lett. 82, 3168 (1999).
  * [11] Y. H. Matsuda, N. Abe, S. Takeyama, H. Kageyama, P. Corboz, A. Honecker, S. R. Manmana, G. R. Foltin, K. P. Schmidt, and F. Mila, Phys.Rev.Lett. 111, 137204 (2013).
  * [12] S. Capponi, O. Derzhko, A. Honecker, A. M. Läuchli and J. Richter, Phys.Rev.B 88, 144416 (2013).
  * [13] S. Nishimoto, N. Shibata and C. Hotta, Nat. Commun. 4, 2287 (2013).
  * [14] D. C. Cabra, M. D. Grynberg, P. C. W. Holdsworth and P. Pujol, Phys. Rev. B 65, 094418 (2002).
  * [15] T. Ono, H. Tanaka, H. Aruga Katori, F. Ishikawa, H. Mitamura, T. Goto, Phys. Rev. B 67, 104431 (2003)
  * [16] T. Inami, Y. Ajiro, T. Goto, J. Phys. Soc. Jpn. 65, 2374 (1996)
  * [17] L. E. Svistov, A. I. Smirnov, and L. A. Prozorova, O. A. Petrenko, L. N. Demianets and A. Ya. Shapiro, Phys. Rev. B 67, 094434 (2003).
  * [18] L. E. Svistov, A. I. Smirnov, L. A. Prozorova, O. A. Petrenko, A. Micheler, N. Büttgen, A. Ya. Shapiro and L. N. Demianets, Phys. Rev. B 74, 024412 (2006).
  * [19] Y. Okamoto, M. Tokunaga, H. Yoshida, A. Matsuo, K. Kindo and Z. Hiroi, Phys. Rev. B 83, 180407(R) (2011).
  * [20] H. Ishikawa, M. Yoshida, K. Nawa, M. Jeong, S. Krämer, M. Horvatić, C. Berthier, M. Takigawa, M. Akaki, A. Miyake, M. Tokunaga, K. Kindo, J. Yamaura, Y. Okamoto, and Z. Hiroi, Phys. Rev. Lett 114, 227202 (2015).
  * [21] R. Okuma , D. Nakamura, and S. Takeyama, Phys. Rev. B 102, 104429 (2020).
  * [22] T. Picot and D. Poilblanc, Phys. Rev. B 91, 064415 (2015).
  * [23] T. Picot, M. Ziegler, R. Orús, and D. Poilblanc, Phys. Rev. B 93, 060407(R)(2016).
  * [24] M.E. Zhitomirsky, Phys. Rev. Lett. 88, 057204 (2002).
  * [25] M. V. Gvozdikova, P.-E. Melchy and M. E. Zhitomirsky, J. Phys.: Condens. Matter 23, 164209 (2011).
  * [26] J. E. Greedan, J. Mater. Chem. 11, 37-53 (2001)
  * [27] A. S. Wills, A. Harrison, C. Ritter and R. I. Smith , Phys. Rev. B 61, 6156 (2000).
  * [28] B. Fåk, F. C. Coomer, A. Harrison, D. Visser and M. E. Zhitomirsky, Eur. Phys. Lett. 81, 17006 (2008).
  * [29] S. Poisson, F. d’Yvoire, NGuyen-Huy-Dung, E. Bretey and P. Berthet, J. of Solid State Chem. 138, 32—40 (1998).
  * [30] If $(x,y,z)$ is the vector basis in which the Knight shift tensor K is diagonal, then the components used here are defined as $K_{\rm iso}=\frac{1}{3}\left(K_{xx}+K_{yy}+K_{zz}\right)$, $K_{\rm ax}=\frac{1}{6}\left(2K_{zz}-K_{xx}-K_{yy}\right)$ and $K_{\rm ani}=\frac{1}{2}\left(K_{yy}-K_{xx}\right)$
  * [31] C. P. Slichter, Principles of Magnetic Resonance, (Springer, 1996).
  * [32] J. N. Reimers and A. J. Berlinsky, Phys. Rev. B 48, 9539 (1993).
  * [33] A. Keren, Phys. Rev. Lett. 72, 3254 (1994).
  * [34] A. B. Harris, C. Kallin and A. J. Berlinsky, Phys. Rev. B 45, 2899 (1992).
  * [35] J. Robert, B. Canals, V. Simonet and R. Ballou, Phys. Rev. Lett. 101, 117207 (2008).
  * [36] M. Taillefumier, J. Robert, C. L. Henley, R. Moessner and B. Canals, Phys. Rev. B 90, 064419 (2014).
  * [37] R. Moessner, J.T. Chalker, Phys. Rev. B 58, 12049 (1998).
  * [38] Y.J. Uemura, A. Keren, K. Kojima, L.P. Le, G. M. Luke, W. D. wu, Y. Ajiro, T. Asano, Y. Kuriyama, M. Mekata, H. Kikuchi and K. Kakurai, Phys. Rev. Lett. 73, 3306 (1994).
  * [39] D. Bono, P. Mendels, G. Collin, N. Blanchard, F. Bert, A. Amato, C. Baines, and A. D. Hillier, Phys. Rev. Lett. 93, 187201 (2004).
  * [40] S.-H. Lee, C. Broholm, G. Aeppli, T. G. Perring, B. Hessen and A. Taylor, Phys. Rev. Lett. 76, 4424 (1996).
  * [41] L. Limot, P. Mendels, G. Collin, C. Mondelli, B. Ouladdiaf, H. Mutka, N. Blanchard and M. Mekata, Phys. Rev. B 65, 144447 (2002).
  * [42] T. Moriya, Prog. Theor. Phys. 16, 23 (1956).
  * [43] E. Kermarrec, A. Zorko, F. Bert, R. H. Colman et al., Phys. Rev. B 90, 205103 (2014).
  * [44] A. J. Vega, P. A. Beckmann, S. Bai and C. Dybowski, Phys. Rev. B 74, 214420 (2006).