Duality and hidden equilibrium in transport models

Rouven Frassek1$\star$, Cristian Giardinà2$\dagger$ and Jorge
Kurchan1$\ddagger$

1 Laboratoire de Physique de l’École Normale Supérieure, ENS, Université PSL,
CNRS, Sorbonne Université, Université de Paris, 75005 Paris, France

2 University of Modena and Reggio Emilia, FIM, Via G. Campi 213/b, 41125
Modena, Italy

$\star$<EMAIL_ADDRESS>$\dagger$<EMAIL_ADDRESS>$\ddagger$<EMAIL_ADDRESS>
## Abstract

A large family of diffusive models of transport that have been considered in
the past years admit a transformation into the same model in contact with an
equilibrium bath. This mapping holds at the full dynamical level, and is
independent of dimension or topology. It provides a good opportunity to
discuss questions of time reversal in out of equilibrium contexts. In
particular, thanks to the mapping one may define the free energy in the non-
equilibrium states very naturally as the (usual) free energy of the mapped
system.

Copyright R. Frassek et al.
This work is licensed under the Creative Commons
Attribution 4.0 International License.
Published by the SciPost Foundation.  |  Received 04-06-2020
Accepted 14-10-2020
Published 20-10-2020  Check for
updates

10.21468/SciPostPhys.9.4.054
---|---

###### Contents

  1. 1 Introduction
  2. 2 Group-based diffusive models
    1. 2.1 Particle transport models
    2. 2.2 Energy transport models
    3. 2.3 Other models
  3. 3 Probability evolution generators
    1. 3.1 Symmetric partially excluded processes
    2. 3.2 KMP processes
    3. 3.3 Integrable version of the energy transport model
  4. 4 Duality transformations
    1. 4.1 Microscopic mapping
    2. 4.2 Perturbative approach
    3. 4.3 Integrable models
  5. 5 Fluctuating hydrodynamics
  6. 6 Discussion

##  1 Introduction

The purpose of this paper is to show that a large class of models of transport
that have been studied by the mathematical physics community in the past years
are, in fact, hidden equilibrium models. Stochastic energy transport models
were already introduced by Kac [2] for a fully connected graph. On a one
dimensional lattice with leads at the ends connected to two reservoirs, two
extensively studied models are: the Kipnis-Marchioro-Presutti (KMP) model [3]
with thermal baths, where energies are randomly redistributed among nearest-
neighbor sites and the symmetric exclusion process (SEP) [4, 5] with particle
baths, where particles jump stochastically between neighboring sites.

In a remarkable series of papers (see [6] for a review), Bertini et al.
constructed the coarse-grained ‘fluctuating hydrodynamic’ limit of such
diffusive systems. They showed how to construct a (WKB / Friedlin-Wentzel)
theory where the amplitude of the effective noise is the small parameter. The
problem maps then onto a Hamiltonian (or Hamilton-Jacobi) field theory in one
space dimension, expressed in terms of a density field $\rho(x)$ and its
conjugate $\hat{\rho}(x)$. The papers contain a second development that comes
as a surprise: there is an explicit transformation that maps the evolution
between a configuration $\rho_{1}(x,t^{\prime})$ and $\rho_{2}(x,t)$ into one
between $\rho_{2}(x,t^{\prime})$ and $\rho_{1}(x,t)$ even when the system is
driven out of equilibrium by the baths. Such transformations are available
when there is a time-reversal symmetry - typically detailed balance - in
systems with equilibrium baths. In order to explain this miracle, Tailleur et
al. [7] showed that in fact, the large deviations around the hydrodynamic
limit of SEP and of KMP could be transformed, via non-local canonical
transformations, into equilibrium. The mapping from ‘downhill’ to ‘uphill’ is
then nothing but the usual one, once the system is transformed into an
undriven system with detailed balance. It was however never clear which
systems enjoy such exceptional feature.

We shall show here that, for a large class of systems, there exists a mapping
from the non-equilibrium process to the equilibrium process already at
microscopic level, and in any dimensionality. The main result of this paper is
then that diffusive systems driven in a non-equilibrium state by (multiple)
external reservoirs are in direct relation with the corresponding equilibrium
systems. See Figure 1 for a pictorial representation of the mapping between a
system with three reservoirs at temperatures $T_{1},T_{2},T_{3}$ and a system
with reservoirs all having the same temperature $T$.

Before going on, it is useful to mention that the particle transport models
like SEP, and the energy transport models like KMP, may be treated in a
unified group-theoretical way that uncovers their similarities – and indeed
opens a whole set of connections with the work on quantum chains, especially
those related to AdS/CFT and high-energy QCD [8, 9, 10, 11, 12].

It turns out that the mapping to equilibrium is not dependent on
integrability, but it derives from a property extensively studied by
probabilists: duality [13]. This property in turn is a consequence of a group
of invariance of the generator of the evolution, that is by no means obvious
when looking at the transition rules.

Although integrability is inessential, when the systems are integrable the
mapping may be constructed explicitly using the conserved charges. As an
example, we will work out in detail (cf. Section 4.3) the explicit formulas
for the mapping for the SEP model following the recent results of [14].

The results of this paper may be extended to quantum stochastic systems driven
in a non-equilibrium steady state by Lindblad reservoirs [15]. For the sake of
clarity we restrict here to classical interacting particle systems and we
shall discuss the quantum counterpart [16] in a forthcoming paper.

$\Gamma_{3},T_{3}$$\Gamma_{2},T_{2}$$\Gamma_{1},T_{1}$

$\Gamma_{3},T$$\Gamma_{2},T$$\Gamma_{1},T$

Figure 1: A probability distribution in the system of the left is mapped onto
one in the system of the right. The mapping is preserved by subsequent
evolution and may be inverted. The stationary state attained in the system of
the left is mapped onto the equilibrium state of the system to the right. The
couplings to the reservoirs $\Gamma_{i}$ remain the same.

##  2 Group-based diffusive models

Our results apply to a whole class of transport models, in particular the
spatial structure does not matter. The system is thus defined on a general
graph $G=(V,E)$, not necessarily planar. The graph has vertices
$V=\\{1,2,\ldots,N\\}$ and edges $E$. On the vertices (labelled by $i\in V$)
we have, when the models are discrete, a given number of particles which may
or may not be limited; otherwise, when the models are continuous, we have a
continuous quantity (‘energy’). In the edge set $E$, we distinguish:

  * •

Internal edges, along which these quantities are transported, according to a
certain probabilistic rule which may depend on the edge. We assume this rule
is symmetric, i.e. the edges are unoriented.

  * •

External edges (‘leads’), connecting vertices to different reservoirs, which
may be particle or thermal baths. The baths are defined by their thermodynamic
properties (chemical potential $\mu_{i}$ or temperature $T_{i}$), and by their
coupling strength $\Gamma_{i}\geq 0$. If the $\mu_{i}$ (or the $T_{i}$) are
all the same, the system eventually reaches thermodynamic equilibrium.
Otherwise, a stationary state with transport is eventually reached.

The most studied graphs of this kind are one-dimensional chains, but we shall
not restrict to these here. The class of models that we can treat is
identified by the following property: they are diffusive models and they
posses a symmetry group, yielding an algebraic description of the process
generator. On the surface, this looks quite restrictive and abstract, but the
class actually contains several natural models considered so far in the
mathematical study of transport. We defer the algebraic description, together
with the identification of the symmetry group, to the next section. Here we
introduce a few models representative of this class; we also recall their
duals.

###  2.1 Particle transport models

#### Symmetric partially excluded processes: SEP$(n)$.

These are systems of particles with repulsive interactions. In the most
general version [17, 18] one assumes maximal occupancy $n\in\mathbb{N}$, which
justifies the name SEP$(n)$. Particles evolve stochastically, following a
sequence of jumps at (exponentially distributed) random times, as follows:

  * •

Rule for transport from a vertex $k$ to a vertex $\ell$: the internal edge
$(k,\ell)$ is chosen with rate $A_{k\ell}>0$. The jumping rate to move one
particle from $k$ to $\ell$ is proportional to the number of particles on the
departure site $n_{k}$, and to the number of ‘holes’ ($n-n_{\ell}$) in the
arrival site $\ell$. Similarly for a jump from $\ell$ to $k$.

  * •

Rule for particle reservoirs: a lead connected to site $i$ is chosen with rate
$\Gamma_{i}>0$. A particle is injected in $i$ with rate $\rho_{i}$ times the
number of holes in $i$, and extracted from $i$ with rate $(1-\rho_{i})$ times
the number of particles in $i$.

#### Dual process.

For any maximal occupation number $n$ the model has a dual process [19, 17,
20] which is obtained by adding extra sites (one for each lead connected to
reservoirs). The dual particles move on the internal edges like the particles
in the original process did (self-duality). Furthermore, a dual particle
sitting on a site that in the original process was connected to a lead, is
absorbed at rate 1 in the corresponding extra site. As a consequence, the dual
dynamics voids the sites of the graphs and all particles are eventually
absorbed at the extra sites.

In this paper, to keep working always on the same Hilbert space, we will avoid
the introduction of the extra sites. As we shall see, the price we will need
to pay is a transformed Hamiltonian which is not stochastic, but only as an
intermediate step.

###  2.2 Energy transport models

#### KMP processes.

The first microscopic stochastic model of energy transport, in the context of
Fourier’s law, was introduced in 1982 by Kipnis, Marchioro and Presutti (KMP)
[3], who studied a chain of oscillators randomly exchanging their energies
between themselves and with two energy reservoirs at different temperatures.
In the same spirit, we consider here a family of systems called (generalized)
KMP processes [21]. The family is labelled by a positive real $s>0$. In the
course of time, the random jumps occur as follows:

  * •

Rule for transport from a vertex $k$ to a vertex $\ell$: the internal edge
$(k,\ell)$ is chosen with rate proportional to $A_{k\ell}$. Energies $z_{k}$
and $z_{\ell}$ of vertices $k$ and $\ell$ are randomly redistributed by
allotting a fraction $B(z_{k}+z_{\ell})$ on site $k$, and the remaining
fraction $(1-B)(z_{k}+z_{\ell})$ on site $l$. Here $B$ is a beta-distributed
random variable with parameters $(2s,2s)$, i.e. it is a law on in the interval
$[0,1]$ with probability density $f(b)\propto b^{2s}(1-b)^{2s}$. Thus the case
$s=1/2$ gives the uniform redistribution rule of the original KMP model.

  * •

A lead $i$ is chosen with rate $\Gamma_{i}$. The energy of the site $i$ is
refreshed with an equilibrium rule of temperature $T_{i}$ (heat bath or Monte
Carlo).

#### Dual process.

Similarly to the previous example, KMP processes have dual processes made of
particles and absorbing extra sites [3, 22, 21, 20]. In the dual process, an
internal edge $(k,\ell)$ is chosen with rate proportional to $A_{k\ell}$ and
then, out of the available $n_{k}+n_{\ell}$ particles, a random number $M$ of
them is put on site $k$ and the remaining number on site $\ell$. The random
number $M$ has a beta-binomial distribution 111The beta-binomial distribution
with parameters $(n,a,b)$ gives the number of heads in $n$ Bernoulli
experiments with a coin whose head probability is fixed but randomly drawn
from a beta distribution with parameters $(a,b)$. with parameters
$(n_{k}+n_{\ell},2s,2s)$. The case $s=1/2$ gives a uniform discrete
redistribution of dual particles. Furthermore, a dual particle sitting on a
site that in the original process was connected to a lead, is absorbed in the
corresponding extra site at rate 1.

###  2.3 Other models

Other group-based diffusive models, that have physical properties similar to
those discussed above, are introduced.

#### Integrable models.

Models that guarantee integrability play a special role as they offer the
possibility of getting explicit formulas. In the framework of partially
excluded processes of Section 2.1, the integrable model is obtained for $n=1$,
which gives the Symmetric Exclusion Process made of particles with hard core
exclusion [13].  The non-equilibrium stationary state is constructed by the
Matrix Product Ansatz [4], from which large deviations of the density and of
the current may be computed, see [5] for a review of the exact solution. An
alternative construction of the stationary measure was recently obtained in
[14] using duality within the algebraic framework of the quantum inverse
scattering method [23]. For the SEP model, the mapping from non-equilibrium to
equilibrium has been described in [7] at the hydrodynamic level, by exploiting
results from the macroscopic fluctuation theory [6]. As we will discuss in
Section 5 the results of [14] allow for such mapping on the microscopic level.

As for the energy transport model of Section 2.2, none of the KMP processes is
integrable. A family of boundary-driven integrable energy transport model has
been recently introduced in [24] (see the discussion on Section 3.3 for the
relation of these models to other models in the probability theory/field
theory literature). For spin $s=1/2$ the model is a Lévy process described as
follows:

  * •

On the internal edge $(k,\ell)$ with energies $z_{k}$ and $z_{\ell}$, an
amount $\alpha$ of energy (with $0<\alpha<z_{k}$) is moved from site $k$ to
site $\ell$ as a Poisson process with intensity $\frac{d\alpha}{\alpha}$.
Similarly for a jump from $\ell$ to $k$.

  * •

At the sites $i$ connected to leads, jumps decreasing the energy by an amount
$\alpha$ (with $0<\alpha<z_{i}$) occur with intensity $\frac{d\alpha}{\alpha}$
and jumps increasing the energy by any amount $\alpha>0$ occur with intensity
$\frac{d\alpha}{\alpha}e^{-\lambda_{i}\alpha}$ with $\lambda_{i}>0$.

The absorbing dual process of this model is made of particles with clustered
jumps. Namely, if $n_{i}$ particles are sitting on site $i$ then a jump of $r$
particles (with $1\leq r\leq n_{i}$) occurs at rate $1/r$. From a site $i$,
the particles are moved, with the same probability, either to a neighbouring
site (connected by an internal edge) or to the extra site (connected by lead).
The general integrable model, labelled by spin value $s>0$, is analyzed in
[24].

#### Diffusions.

For energy transport models one can replace the evolution through Markovian
jumps of the generalized KMP processes with a diffusions having continuous
paths. This gives the so-called Brownian energy processes introduced in [21]
or the Brownian momentum processes [25, 22]. The dual processes of those are
called inclusion processes [26], they are the bosonic counterpart of partially
excluded processes. More precisely, in the model with spin $s>0$, the
particles jump from vertex $k$ to a vertex $\ell$ with rate
$n_{k}(2s+n_{\ell})$.

#### Inhomogeneiteis and/or multispecies.

One can add inhomogeneities to allow different maximal occupancy at each site
keeping the key duality property [27]. Several kinds of particles or
continuous quantities, with mutual exclusion properties, are also possible,
leading to higher rank algebras [28, 29].

##  3 Probability evolution generators

In this section we develop the algebraic description of the Markov processes
discussed above. Such algebraic approach relies on applying to the equation
describing the evolution of the probability for the system of interest – the
so-called forward Kolmogorov equation – tools that are standard for the
Schrödinger equation. It is well-known [17] that exclusion processes can be
expressed in terms of a Hamiltonian (in fact, the generator of the probability
evolution) with spin operators of the $\mathfrak{su}(2)$ Lie algebra. For the
energy transport models it was realized in [22, 21] that also KMP processes
can be described in this way, just replacing the compact spin algebra by the
non-compact $\mathfrak{su}(1,1)$ algebra. See [30] for an extended
presentation of the algebraic approach, and its use in finding dual processes.

###  3.1 Symmetric partially excluded processes

The evolution operator of these processes can be written as the Hamiltonian
$H$ of the half-integer spin $j=n/2$ ferromagnet [17]

$\displaystyle H$ $\displaystyle=$ $\displaystyle\sum_{k,l\in
V^{2}}A_{kl}\left(J^{+}_{k}J^{-}_{l}+J^{-}_{k}J^{+}_{l}+2J^{0}_{k}J^{0}_{l}-2j^{2}\right)$
$\displaystyle+$ $\displaystyle\sum_{i\in
V}\Gamma_{i}\left[\rho_{i}(J^{+}_{i}+J^{0}_{i}-j)+(1-\rho_{i})(J^{-}_{i}-J^{0}_{i}-j)\right]\,.$

The operators $J^{+}_{i},J^{-}_{i},J^{0}_{i}$ act on the Hilbert space
corresponding to $0\leq r\leq n$ particles per site $i$ as follows:

$\displaystyle J^{+}_{i}|r\rangle_{i}$ $\displaystyle=$
$\displaystyle(2j-r)|r+1\rangle_{i}$ $\displaystyle J^{-}_{i}|r\rangle_{i}$
$\displaystyle=$ $\displaystyle r|r-1\rangle_{i}$ $\displaystyle
J^{0}_{i}|r\rangle_{i}$ $\displaystyle=$ $\displaystyle(r-j)|r\rangle_{i}\,.$
(2)

We order the states like $|2j\rangle,|2j-1\rangle,\ldots,|0\rangle$ so that
the $J^{-}_{i}$ operators are lower-triangular. The Hamiltonian (3.1) acts on
the tensor product space with states $\otimes_{i\in V}|r\rangle_{i}$ and it is
easy to read off the (negative) transition rates. For instance the first term
$-J^{+}_{k}J^{-}_{l}$ gives the rate $n_{\ell}(2j-n_{k})$ for the jump of a
particle from site $\ell$ with $n_{\ell}$ particles to site $k$ with $n_{k}$
particles. The operators $J^{+}_{i},J^{-}_{i},J^{0}_{i}$ satisfy the
commutation relations of the $\mathfrak{su}(2)$ Lie algebra:

$[J_{i}^{0},J_{i}^{\pm}]=\pm
J_{i}^{\pm}\;\;\;\;\;\;\;\;\;\;[J_{i}^{-},J_{i}^{+}]=-2J_{i}^{0}\,.$ (3)

Representations are labeled by the spin value $j$ related to the squared
angular momentum operator via:

$\displaystyle{(\vec{J}_{i})^{2}}|j,m\rangle_{i}$ $\displaystyle=$
$\displaystyle\left([J_{i}^{0}]^{2}+\frac{1}{2}[J_{i}^{+}J_{i}^{-}+J_{i}^{-}J_{i}^{+}]\right)|j,m\rangle_{i}$
$\displaystyle=$ $\displaystyle j(j+1)|j,m\rangle_{i}\,.$

Here $j=n/2$, so that the ordinary SEP with $(0,1)$ occupation corresponds to
a representation of spin $j=1/2$. For a given $j$ half-integer, there are
$2j+1$ eigenstates of the $J^{0}_{i}$ operator

$J^{0}_{i}|j,m\rangle_{i}=m|j,m\rangle_{i}\,,$ (5)

with $m=r-j$ and $m\in\\{-j,-(j-1),\ldots,j\\}$. In other words, for a given
half-integer $j$ we identify $|r\rangle_{i}=|j,r-j\rangle_{i}$, cf. (3.1).

###  3.2 KMP processes

We discuss here the diffusion processes associated to KMP, i.e. having the
same algebraic structure. Moreover, following Kac’s idea of energy transfer
via random collisions [2], we think of kinetic energies, and thus discuss a
velocity diffusion process where energy transfer is obtained by ‘random
rotations’ of velocity vectors [31].

Thus we consider, on each site $i\in V$, a velocity vector with
$M\in\mathbb{N}$ components, i.e. $v_{i,\alpha}$ with $\alpha=1,\ldots,M$.
Suppose that they evolve as a diffusion process with Markov generator (acting
on the core of $\mathscr{C}^{\infty}$ functions with compact support)

${L}=\sum_{k,\ell\in V^{2}}A_{k\ell}{L}_{k\ell}+\sum_{i\in
V}\Gamma_{i}{L}_{i}\,,$ (6)

where

${L}_{k\ell}=\sum_{\alpha,\beta=1}^{M}\left(v_{k,\alpha}\frac{\partial}{\partial
v_{\ell,\beta}}-v_{\ell,\beta}\frac{\partial}{\partial
v_{k,\alpha}}\right)^{2}\,,\qquad\qquad{L}_{i}=\sum_{\alpha=1}^{M}\left(T_{i}\frac{\partial^{2}}{\partial
v_{i,\alpha}^{2}}-v_{i,\alpha}\frac{\partial}{\partial
v_{i,\alpha}}\right)\,.$

Each term in ${L}_{k\ell}$ represents a rotation in the plane
$(v_{k,\alpha},v_{\ell,\beta})$, therefore it conserves the total kinetic
energy of the edge $v^{2}_{k,\alpha}+v^{2}_{\ell,\beta}$; on top of this,
${L}_{i}$ gives a Langevin bath with temperature $T_{i}$ and strength of
coupling $\Gamma_{i}$.

The Fokker-Planck equation, yielding the evolution of the time-dependent
probability density $p(v,t)$, reads

$\frac{\partial p(v,t)}{\partial t}=L^{*}p(v,t)\,,$ (7)

where $L^{*}$ is the adjoint (in $\mathscr{L}^{2}(dv)$) of $L$. Expanding the
brackets in (6) and defining

${S}_{i}^{+}=\sum_{\alpha=1}^{M}\frac{1}{2}v_{i,\alpha}^{2}\qquad\qquad{S}_{i}^{-}=\sum_{\alpha=1}^{M}\frac{1}{2}\frac{\partial^{2}}{\partial
v_{i,\alpha}^{2}}\qquad\qquad{S}_{i}^{0}=\sum_{\alpha=1}^{M}\frac{1}{4}\left(v_{i,\alpha}\frac{\partial}{\partial
v_{i,\alpha}}+\frac{\partial}{\partial v_{i,\alpha}}v_{i,\alpha}\right)$

we can rewrite this as an imaginary time Schrödinger equation with Hamiltonian
$H=L^{*}$ and

$\begin{split}L&=\sum_{k,\ell\in
V^{2}}A_{k\ell}\left(S^{+}_{k}S^{-}_{l}+S^{-}_{k}S^{+}_{l}-2S^{0}_{k}S^{0}_{l}+\frac{M^{2}}{8}\right)+\sum_{i\in
V}\Gamma_{i}\left(T_{i}S^{-}_{i}-S^{0}_{i}+\frac{M}{4}\right)\,.\end{split}$
(8)

The generators $S^{a}$ satisfy the $\mathfrak{su}(1,1)$ Lie algebra relations:

$[S_{i}^{0},S_{i}^{\pm}]=\pm
S_{i}^{\pm}\;\;\;\;\;\;\;\;\;\;[S_{i}^{-},S_{i}^{+}]=2S_{i}^{0}\,.$ (9)

Representations are labeled in an analogous manner by the spin value $s>0$
with the Casimir operator:

$\displaystyle(\vec{S}_{i})^{2}|s,m\rangle_{i}$ $\displaystyle=$
$\displaystyle\left([S_{i}^{0}]^{2}-\frac{1}{2}[S_{i}^{+}S_{i}^{-}+S_{i}^{-}S_{i}^{+}]\right)|s,m\rangle_{i}$
$\displaystyle=$ $\displaystyle s(s-1)|s,m\rangle_{i}\,.$

Note the sign difference with respect to $\mathfrak{su}(2)$. To identify the
representation (i.e. the value of $s$), we compute the eigenvalue of $S^{0}$
and $(\vec{S})^{2}_{i}$ when acting on the (constant) lowest weight state
$|-\rangle_{i}$, which is the eigenvector with zero eigenvalue of $L^{*}$:

$S^{0}_{i}|-\rangle_{i}=\frac{M}{4}|-\rangle_{i}\,,$ (11)
$(\vec{S}_{i})^{2}|-\rangle_{i}=\frac{M}{4}\left(\frac{M}{4}-1\right)|-\rangle_{i}\,.$
(12)

Hence $s=M/4$, and in particular $s=1/4$ for the process with one velocity per
site.

###  3.3 Integrable version of the energy transport model

For integrable transport models one may perform explicitly the mapping from
non-equilibrium to equilibrium. However, in principle integrability is not
necessary for the mapping to exist.

Once one identifies the KMP energy diffusion model with an
$\mathfrak{su}(1,1)$ chain, the question immediately arises as to its
integrability. While the chain with Hamiltonian (8) is not integrable, its
integrable cousin has been known for a long time [32]. It remains to check
that such a model may indeed be interpreted as a stochastic system, which it
is indeed [24]. It turns out that the particle version of such a system has
also been introduced in the probability literature, in the asymmetric form, by
Sasamoto-Wadati [33] for spin $s=1/2$ and by Povolotsky [34] and Barraquand-
Corwin [35] for higher spin, just constructing it on the basis of the Bethe
ansatz properties – without making a connection with integrable spin chain
explicit. For the latter we refer the reader to [24] and [36].

The boundary driven integrable version of the energy transport model recently
introduced in [24] is the Lévy process described in Section 2.3. In the
corresponding particle version, the probability evolves with a generator of
the form (6) where now

$\displaystyle L_{k,\ell}$
$\displaystyle=2(\psi(\mathbb{S}_{k,\ell})-\psi(2s))\,,$ (13) $\displaystyle
L_{i}$
$\displaystyle=e^{\beta_{i}S^{+}_{i}}e^{\frac{1}{\beta_{i}-1}S^{-}_{i}}\psi(S^{0}_{i}+s)e^{-\frac{1}{\beta_{i}-1}S^{-}_{i}}e^{-\beta_{i}S^{+}_{i}}-\psi(2s)\,.$

Here $\psi$ is the digamma function (the logarithmic derivative of the gamma
function), $\beta_{i}$ is a parameter tuning the density of the reservoirs and
the operator $\mathbb{S}_{k,\ell}$ is related to the two-site Casimir acting
on two sites $k$ and $\ell$ via
$(\vec{S}_{k}+\vec{S}_{\ell})^{2}=\mathbb{S}_{k,\ell}(\mathbb{S}_{k,\ell}-1)$.

##  4 Duality transformations

As discussed in previous sections, dual Markov processes of transport models
are usually obtained by introducing extra sites, one for each lead with
$\Gamma_{i}>0$. The key property of the dual process is the absorbing property
of the extra sites, which is reflected in the fact that the dual generator has
a triangular structure. Here we fully exploit the consequences of this crucial
property. To keep working always on the same Hilbert space we avoid the extra
sites, yielding a non-stochastic transformed Hamiltonian.

###  4.1 Microscopic mapping

Let us for definiteness concentrate on the SEP$(n)$ case (the same reasoning
can be repeated in the other cases). As pointed out by Schütz and Sandow [17],
one can make a “duality” transformation of the Hamiltonian $H$ in (3.1) as
follows. From the Hadamard formula for the $\mathfrak{su}(2)$ algebra one
obtains

$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{e}^{\mu
J^{\pm}}J^{0}e^{-\mu J^{\pm}}}$ $\displaystyle=$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{J}^{0}\mp\mu
J^{\pm}}$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{e}^{\mu
J^{\pm}}J^{\mp}e^{-\mu J^{\pm}}}$ $\displaystyle=$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{J}^{\mp}\pm
2\mu J^{0}-\mu^{2}J^{\pm}}\,.$ (14)

Setting $\mu=1$, we get that the boundary terms of (3.1) transform into

$e^{J^{+}_{i}}\Gamma_{i}\left[\rho_{i}(J^{+}_{i}+J^{0}_{i}-j)+(1-\rho_{i})(J^{-}_{i}-J^{0}_{i}-j)\right]e^{-J^{+}_{i}}=\Gamma_{i}\left[(J^{0}_{i}-j)+(1-\rho_{i})J^{-}_{i}\right]\,.$
(15)

This is essentially equivalent to the general notion of duality of Markov
processes [21, 30]. Indeed (15) can be rewritten as

$\left[\rho_{i}(J^{+}_{i}+J^{0}_{i}-j)+(1-\rho_{i})(J^{-}_{i}-J^{0}_{i}-j)\right]^{tr}R_{i}e^{-J^{+}_{i}}=R_{i}e^{-J^{+}_{i}}\left[\rho_{i}J^{-}_{i}-J^{0}_{i}-j\right]\,,$
(16)

where tr denotes transposition and $R_{i}$ is the diagonal matrix (with
entries the inverse of the Binomial distribution) such that
$R_{i}J_{i}^{\pm}R_{i}^{-1}=J_{i}^{\mp}.$ If we introduce an extra site,
called ${c(i)}$ and associated to site $i$, this can be further rewritten as

$H_{i}^{tr}D_{i}=D_{i}H_{i}^{dual}\,,$ (17)

where the duality function $D_{i}=\sum_{m}{\rho_{i}^{m}\langle
m|R_{i}e^{J_{i}^{+}}}$, and $H_{i}^{dual}=a^{+}_{c(i)}J^{-}_{i}-J^{0}_{i}-j$
is the stochastic Hamiltonian that describes the dual absorbing process at
site $i$ (here $a^{+}_{c(i)}$ is a bosonic creation operator in the extra
site). Equation (17) is the standard form of a duality relation between two
Markov processes [13, 30], that allows to connect expectations of the original
process to expectation of the dual process. For a general discussion on the
relation between duality relation and symmetries we refer the reader to [21]
and reference therein; for a constructive approach aiming to introducing
group-based models with duality we refer to [37].

Transforming further,  with an arbitrary $0<\bar{\rho}<1$, we have

$\displaystyle
e^{-(1-\bar{\rho})J^{-}_{i}}e^{J^{+}_{i}}\Gamma_{i}\left[\rho_{i}(J^{+}_{i}+J^{0}_{i}-j)+(1-\rho_{i})(J^{-}_{i}-J^{0}_{i}-j)\right]e^{-J^{+}_{i}}e^{(1-\bar{\rho})J^{-}_{i}}=$
$\displaystyle=\Gamma_{i}\left[(J^{0}_{i}-j)+(\bar{\rho}-\rho_{i})J^{-}_{i}\right]\,.$
(18)

All in all, defining
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{J}^{\pm}_{tot}=\sum_{i\in
V}J^{\pm}_{i}}$ we get that the original Hamiltonian (3.1) is mapped into:

$\displaystyle
H_{B}=e^{-(1-\bar{\rho})J^{-}_{tot}}e^{J^{+}_{tot}}He^{-J^{+}_{tot}}e^{(1-\bar{\rho})J^{-}_{tot}}$
$\displaystyle=$ $\displaystyle H_{0}+\underbrace{\sum_{i\in
V}\Gamma_{i}(\bar{\rho}-\rho_{i})J^{-}_{i}}_{B}\,.$ (19)

This defines $B$ as the underbraced term and $H_{0}$ as

$H_{0}=\sum_{k,l\in
V^{2}}A_{kl}\left(J^{+}_{k}J^{-}_{l}+J^{-}_{k}J^{+}_{l}+2J^{0}_{k}J^{0}_{l}-2j^{2}\right)+\sum_{i\in
V}\Gamma_{i}(J^{0}_{i}-j)\,.$ (20)

In $H_{0}$, the bulk term of (3.1) is left invariant, precisely because of the
global $\mathfrak{su}(2)$ symmetry. If all the $\rho_{i}=\bar{\rho}$ are the
same, we may eliminate the $B$ term. Note that $H_{B}$ does not correspond to
a probability conserving process.

Since the Hamiltonian $H_{0}$ commutes with $J^{0}_{tot}$, we may
simultaneously diagonalize the two operators so that
$J^{0}_{tot}|\Lambda_{k},m_{tot}\rangle=m_{tot}|\Lambda_{k},m_{tot}\rangle$
and $H_{0}|\Lambda_{k},m_{tot}\rangle=\Lambda_{k}|\Lambda_{k},m_{tot}\rangle$
where $k=1,\ldots,(n+1)^{|V|}$ and $m_{tot}$ labels the corresponding number
of particles. If we order this set of vectors by decreasing values of
$m_{tot}$, and since each $J^{-}_{i}$ lowers the value of $m_{tot}$ by one,
the matrix elements of
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{B}}$
are lower diagonal in this base. Hence it becomes apparent that the
$J^{-}_{i}$ in $B$ do not modify the spectrum. We conclude that the spectra of
$H_{B}$ and $H$ coincide with that of $H_{0}$. We remark that such
isospectrality property for spin chain Hamiltonians differing by operators
that are lower triangular in the right basis was also observed in [38].  By
inspection of (20) we see that the spectrum depends only on the
$\Gamma_{i}$’s, but not on the $\rho_{i}$’s. We conclude that there exists an
operator $W$ such that

$H_{B}=WH_{0}W^{-1}.$ (21)

In particular, denoting $|R_{0}^{k}\rangle\equiv|\Lambda_{k},m_{tot}\rangle$,
$\langle L_{0}^{k}|$ and $|R_{B}^{k}\rangle$, $\langle L_{B}^{k}|$ the right
and left eigenvectors of $H_{0}$ and $H_{B}$, and $\Lambda_{k}$ the
corresponding eigenvalues labelled by $k$, we have the relation

$W|R_{0}^{k}\rangle=|R_{B}^{k}\rangle\,.$ (22)

The transformation $W$ can be written as

$W=\lim_{\epsilon\rightarrow
0}\epsilon\sum_{k}\left[\Lambda_{k}-H_{B}+\epsilon\right]^{-1}|R_{0}^{k}\rangle\langle
L_{0}^{k}|\,.$ (23)

Alternatively, one may also write an equivalent expression that, introducing
an additional integral, does not require the knowledge of the spectrum of
$H_{B}$:

$W=\lim_{\epsilon\rightarrow 0}\epsilon\lim_{\bar{\epsilon}\rightarrow
0}\Im\;\oint d\Lambda\left[\Lambda-H_{B}+\epsilon\right]^{-1}\left[\Lambda-
H_{0}+i\bar{\epsilon}\right]^{-1}\,.$ (24)

Here $\Im$ denotes the imaginary part. This expression in fact is essentially
formal and it does not seem to be useful at this stage.

To recap, we have up to now a similarity transformation relating $H$ to
$H_{0}$:

$W^{-1}e^{-(1-\bar{\rho})J^{-}_{tot}}e^{J^{+}_{tot}}He^{-J^{+}_{tot}}e^{(1-\bar{\rho})J^{-}_{tot}}W=H_{0}.$
(25)

One further transformation gives the final result

$P\;H\;P^{-1}=H_{eq}\,,$ (26)

with

$P=e^{-J^{+}_{tot}}e^{(1-\bar{\rho})J^{-}_{tot}}W^{-1}e^{-(1-\bar{\rho})J^{-}_{tot}}e^{J^{+}_{tot}}$
(27)

and $H_{eq}$ denoting the Hamiltonian (3.1) in equilibrium with
$\rho_{i}=\bar{\rho}$ for all $i\in V$. Thus we found a mapping of the non-
equilibrium process $H$ (with chemical potentials $\rho_{i}$) to the
equilibrium process with Hamiltonian $H_{eq}$ with reservoirs having the same
chemical potential $\bar{\rho}$ at all leads! The coupling intensities
$\Gamma_{i}$ of $H$ and $H_{eq}$ remain the same. If $|\psi\rangle$ is the ket
vector which encodes the probability distribution of the original process as
$\frac{d}{dt}|\psi\rangle=H|\psi\rangle$, and $|\psi\rangle_{eq}$ is the ket
vector which encodes the probability distribution of the transformed process
as $\frac{d}{dt}|\psi\rangle_{eq}=H_{eq}|\psi\rangle_{eq}$ then, as
consequence of (26), one has $|\psi\rangle=P^{-1}|\psi\rangle_{eq}$. It would
be interesting to study the relation of our approach to the matrix product
ansatz [4] and [45, 46, 47].

For later use, it is convenient to define $A$ by

$W=e^{A}\,.$ (28)

The transformation (26) may be written in terms of $A$ as

$e^{-[J^{+}_{tot},\;\;]}e^{(1-\bar{\rho})[J^{-}_{tot},\;\;]}e^{-[{A},\;\;]}e^{-(1-\bar{\rho})[J^{-}_{tot},\;\;]}e^{[J^{+}_{tot},\;\;]}{H}={H}_{eq}\,,$
(29)

where we used the notation $e^{[X,\;\;]}Y=e^{X}Ye^{-X}$.

###  4.2 Perturbative approach

For all systems, the expression for $W$ we introduced in the previous
paragraph may be evaluated perturbatively.  To this aim we insert $\Delta$, a
bookkeeping parameter we shall set to one at the end, to order the
perturbation series. We start by writing

$\displaystyle H_{\Delta}=H_{0}+\Delta B\,,$ (30)

and define $A_{\Delta}$ by

$W_{\Delta}=e^{A_{\Delta}}\,,$ (31)

which reduce to $H_{B}$ and $W$ for $\Delta=1$ respectively. The relation
between $H_{B}$ and $H_{0}$ in (21) becomes

$e^{-A_{\Delta}}\left[H_{0}+\Delta
B\right]e^{A_{\Delta}}=e^{-[A_{\Delta},\;\;]}H_{\Delta}=H_{0}\,.$ (32)

Then we may write a power series for $A_{\Delta}$ by putting
$A_{\Delta}=\Delta A^{(1)}+\Delta^{2}A^{(2)}+...$, so that

$\displaystyle B$ $\displaystyle=$ $\displaystyle[A^{(1)},H_{0}]$
$\displaystyle\frac{1}{2}[B,A^{(1)}]$ $\displaystyle=$
$\displaystyle[A^{(2)},H_{0}]$ $\displaystyle\vdots$

or, in the basis where $H_{0}$ is diagonal

$\displaystyle B_{ij}(\Lambda_{i}-\Lambda_{j})^{-1}$ $\displaystyle=$
$\displaystyle A^{(1)}_{ij}$
$\displaystyle\frac{1}{2}[B,A^{(1)}]_{ij}(\Lambda_{i}-\Lambda_{j})^{-1}$
$\displaystyle=$ $\displaystyle A^{(2)}_{ij}$ $\displaystyle\vdots$

This allows to compute the coefficients $A^{(i)}$ of $W_{\Delta}$ recursively.

###  4.3 Integrable models

As we have seen in the previous sections, of all the transformations we need
to perform, the only one that is non-local in the sites, and does not have an
explicit expression, is $W$. Consider now the expressions for the matrix $W$
(23): it depends on $H_{0}$ and $H_{B}$. It is easy to see that if we have an
operator $Q_{B}$ (conserved charge) that commutes with $H_{B}$, and an
operator $Q_{0}$ that commutes with $H_{0}$, with $Q_{0}$ and $Q_{B}$ iso-
spectral, we may just as well substitute the $H$’s by the $Q$’s in (23). The
change of basis induced by $W$ will be the same. If, furthermore, the spectrum
of $Q_{0}$ is known explicitly, we may obtain an important simplification and
sum up the perturbative series. This is what integrability gives us.

We illustrate this procedure in the prototype case of boundary-driven SEP on a
1D chain of $N$ sites. In the bulk, particles jumps at rate $1$ to their
nearest neighbors, provided there is space in the arrival site. The reservoir
connected to the first site (resp. last site) inject particles at rate
$\alpha$ (resp. $\delta$) and remove particles at rate $\gamma$ (resp.
$\beta$). We consider the Hamiltonian (3.1) with $j=1/2$ and define
$J_{i}^{\pm}=\sigma_{i}^{\pm}\;,J_{i}^{0}=\sigma^{0}/2$ with
$\sigma^{\pm}_{i}=(\sigma^{x}_{i}\pm\sigma^{y}_{i})/2$ and
$\sigma^{0}_{i}=\sigma^{z}_{i}$ where $\sigma_{i}$ are the Pauli matrices at
site $i$, so that

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N-1}\left(\sigma^{+}_{i}\sigma^{-}_{i+1}+\sigma^{-}_{i}\sigma^{+}_{i+1}+\frac{1}{2}\sigma^{0}_{i}\sigma^{0}_{i+1}-\frac{1}{2}\right)$
$\displaystyle+$
$\displaystyle\sum_{i\in\\{1,N\\}}\Gamma_{i}\left[\rho_{i}(\sigma^{+}_{i}+\sigma^{0}_{i}/2-1/2)+(1-\rho_{i})(\sigma^{-}_{i}-\sigma^{0}_{i}/2-1/2)\right]\,,$

with
$\Gamma_{1}=\alpha+\gamma\,,\Gamma_{N}=\delta+\beta\,,\rho_{1}=\alpha/(\alpha+\gamma)\,,\rho_{N}=\delta/(\delta+\beta)$.
It is convenient to set $\bar{\rho}=\rho_{N}$ so that, after performing the
transformation in (19), one gets $H_{B}=H_{0}+B$ with

$H_{0}=\sum_{i=1}^{N-1}\left(\sigma^{+}_{i}\sigma^{-}_{i+1}+\sigma^{-}_{i}\sigma^{+}_{i+1}+\frac{1}{2}\sigma^{0}_{i}\sigma^{0}_{i+1}-\frac{1}{2}\right)+\sum_{i\in\\{1,N\\}}\frac{\Gamma_{i}}{2}\left(\sigma^{0}_{i}-1\right)$
(36)

and

$B=\Gamma_{1}(\rho_{N}-\rho_{1})\sigma_{1}^{-}\,.$ (37)

The study of the Hamiltonian $H_{0}$ using the coordinate Bethe ansatz goes
back to [39, 40] and to Sklyanin [23] using the quantum inverse scattering
method. For its relation to the Hamiltonian of the SEP within the latter
framework we refer the reader to [41, 42, 43]. The isospectrality of the
Hamiltonians (4.3) and (36) is straight-forwardly obtained. The transformation
$W$, cf. Section 4.1, however was only obtained recently [14]. Without going
into the details of the quantum inverse scattering method [32], which allows
to obtain the conserved charge, we take the operator $Q_{B}$ as given and show
how the results of [14] fit into our perturbative framework. We further
present some alternative ways of writing the resulting similarity
transformation, see in particular (47) and (51).

As mentioned before, the explicit computation of the similarity transformation
$W$ relies on the existence of a “charge” operator

$Q_{\Delta}=Q_{0}+\Delta Q_{-}\,,$ (38)

that commutes with the Hamiltonian $H_{\Delta}$ for any $\Delta$, i.e.

$[H_{\Delta},Q_{\Delta}]=0\,.$ (39)

The operator $Q_{\Delta}$ has the same eigenvectors as the Hamiltonian
$H_{\Delta}$ which reduce to the eigenvectors of $H_{0}$ (or $Q_{0}$) for
$\Delta=0$ and to the ones of $H_{B}$ for $\Delta=1$. Thus the transformation
$W_{\Delta}$ maps $Q_{0}$ to $Q_{\Delta}$:

$Q_{\Delta}=W_{\Delta}Q_{0}W_{\Delta}^{-1}\,,$ (40)

cf. (32). The operator $Q_{\Delta}$ follows from the transfer matrix of the
integrable spin chain, see [14] for further details. The unperturbed part only
depends on the total spin

$Q_{0}=\frac{\sigma^{0}_{tot}}{2}\Gamma_{1}\Gamma_{N}\left(\frac{\sigma^{0}_{tot}}{2}+\Gamma_{1}^{-1}+\Gamma_{N}^{-1}\right)\,,$
(41)

while the lower triangular part is given by the non-local expression

$Q_{-}=\Gamma_{1}(\rho_{N}-\rho_{1})\left(\sigma_{tot}^{-}+\Gamma_{N}\sum_{i=1}^{N}\sigma^{-}_{i}\left(\frac{\sigma^{0}_{i}-1}{2}+\sum_{k=i+1}^{N}\sigma^{0}_{k}\right)\right)\,.$
(42)

The similarity transformation can then be computed using the algebraically
more simple operator $Q_{\Delta}$:

$W_{\Delta}=\lim_{\epsilon\rightarrow
0}\epsilon\sum_{m}\left[\Lambda^{Q}_{m}-Q_{\Delta}+\epsilon\right]^{-1}|R_{0}^{m}\rangle\langle
L_{0}^{m}|\,.$ (43)

The eigenvalues $\Lambda^{Q}_{m}$ can be read off immediately from the
explicit form of $Q_{0}$ in (41). One way to determine the $W_{\Delta}$ is by
recursively solving (40) in powers of $\Delta$ as discussed in Section 4.2.
Writing

$W_{\Delta}=1+\sum_{k=1}^{N}\Delta^{k}g_{k}$ (44)

we get

$g_{k}Q_{0}=Q_{0}g_{k}+Q_{-}\,g_{k-1}\qquad\qquad k=1,\ldots N\,,$ (45)

with $g_{0}=1$. From this it follows that $g_{k}$ must annihilate $k$
particles which implies its exchange relation with $Q_{0}$:

$g_{k}Q_{0}(\sigma_{tot}^{0})=Q_{0}(\sigma_{tot}^{0}+2k)g_{k}\,.$

This allows us to solve the recursion and explicitly write $g_{k}$ in terms of
$Q_{0}$ and $Q_{-}$. By using then the definition of $Q_{0}$ (41) we find

$W_{\Delta}=\sum_{k=0}^{N}\frac{\Delta^{k}}{k!}\left(\frac{Q_{-}}{\Gamma_{1}\Gamma_{N}}\right)^{k}\frac{\Gamma\left(\sigma^{0}_{tot}+\Gamma_{1}^{-1}+\Gamma_{N}^{-1}-k\right)}{\Gamma\left(\sigma^{0}_{tot}+\Gamma_{1}^{-1}+\Gamma_{N}^{-1}\right)}\,.$
(46)

The final expression for $W$ can be obtained by taking $\Delta=1$. From this
expression, the non-local character of the $W$-transform is clearly seen, cf.
(42).

Curiously, $W$ can be written as an exponential function

$W=\exp\left(\sum_{k=1}^{N}A^{(k)}\right)\,,$ (47)

with

$A^{(k)}=\gamma_{k}\frac{(-1)^{k+1}}{k!\,\Gamma_{1}^{k}\,\Gamma_{N}^{k}}\frac{\Gamma\left(1+\sigma_{tot}^{0}+\Gamma_{1}^{-1}+\Gamma_{N}^{-1}\right)}{\Gamma\left(2k+\sigma_{tot}^{0}+\Gamma_{1}^{-1}+\Gamma_{N}^{-1}\right)}Q_{-}^{k}\,,$
(48)

where the function $\gamma_{k}$ are related to the Clebsch-Gordan coefficients
(compare [44]),
$\gamma_{k}=(1,1,3,14,80,468,2268,10224,313632,9849600,\ldots)$. This
exponential form of $W$ allows to read off its inverse.

The expression for $W$ may be expressed in terms of Bessel function which may
be more suitable to take the hydrodynamic limit. To see this, we define
${\mathrm{P}}_{j}$ the projector on the subspace of total spin $j$. Then we
may rewrite (46) for $\Delta=1$ as

$W=\sum_{j}\sum_{k=0}^{N}\frac{1}{k!}\left(\frac{Q_{-}}{\Gamma_{1}\Gamma_{N}}\right)^{k}\frac{\Gamma\left(2j+\Gamma_{1}^{-1}+\Gamma_{N}^{-1}-k\right)}{\Gamma\left(2j+\Gamma_{1}^{-1}+\Gamma_{N}^{-1}\right)}\;{\mathrm{P}}_{j}\,.$
(49)

Now, writing the gamma function in the numerator as $\Gamma(z)=\int
dxx^{z-1}e^{-x}$ we may perform the sum over $k$, to obtain:

$W=\sum_{j}\int
dx\;e^{-x+\frac{Q_{-}}{\Gamma_{1}\Gamma_{N}}x^{-1}}\;x^{\left(2j+\Gamma_{1}^{-1}+\Gamma_{N}^{-1}-1\right)}\frac{1}{\Gamma\left(2j+\Gamma_{1}^{-1}+\Gamma_{N}^{-1}\right)}\;{\mathrm{P}}_{j}\,.$
(50)

We recognize the generator of the Bessel function $K_{n}(z)$, so that we may
write:

$W=2\sum_{j}\;\frac{K_{-2j-\Gamma_{1}^{-1}-\Gamma_{N}^{-1}}\left[2\sqrt{\frac{-Q_{-}}{\Gamma_{1}\Gamma_{N}}}\,\right]}{\Gamma\left(2j+\Gamma_{1}^{-1}+\Gamma_{N}^{-1}\right)}\;\left(\sqrt{\frac{-Q_{-}}{\Gamma_{1}\Gamma_{N}}}\right)^{2j+\Gamma_{1}^{-1}+\Gamma_{N}^{-1}}\;{\mathrm{P}}_{j}\,.$
(51)

This is an explicit expression for $W$ acting on a function of given
$\sigma^{0}_{tot}$.

One final remark: when we map a probability distribution via this
transformation, we will not necessarily obtain a function that is positive
definite, and may be interpreted as a probability. This is no real problem if
the system is ergodic: one may add to the function an equilibrium distribution
multiplied by a sufficiently large factor, and normalize the result to one.
The evolution of this linear combination, which may be interpreted as a bona-
fide probability, is just a linear combination of the one of the original
vector and the equilibrium one, and the latter does not evolve.

##  5 Fluctuating hydrodynamics

Fluctuating hydrodynamics is obtained by coarse graining a system: the
quantities obtained depend on space and have fluctuations that are the
smaller, the larger the coarse graining. The hydrodynamic limit consists
formally of two parts: one first goes to a continuous chain, and then performs
the long-wavelength (infrared) limit. These cases amount to a ‘semiclassical’
approximation in the following sense: just as the probability evolution can be
seen as a Schrödinger equation with imaginary time, its coarse-grained limit
correspond to the semiclassical limit of it – the role of $\hbar$ being played
by the intensity of fluctuations. The (WKB) treatment of this with the usual
classical tools is variously known as ‘Freidlin-Wentzel’ or ‘Hamilton-Jacobi’
approach in the mathematical physics literature. Exactly the same problem
arises in the AdS/CFT literature, where the long-wavelength limit becomes the
high angular momentum regime of the string [48]. Let us stress again the role
of group theoretical nature of these problems: it is because of this structure
that ‘coarse graining’ amounts to ‘large spin representation’, and through it
the semiclassical limit.

At the fluctuating hydrodynamics level, the transformations connecting non-
equilibrium to equilibrium can be read from the previous Section 4.1; just
changing commutators into Poisson brackets:

$\\{{\cal J}_{k}^{0},{\cal J}_{k}^{\pm}\\}=\pm{\cal
J}_{k}^{\pm}\;\;\;\;\;\;\;\;\;\;\\{{\cal J}_{k}^{-},{\cal
J}_{k}^{+}\\}=-2{\cal J}_{k}^{0}\,.$ (52)

For the one-dimensional SEP(2j), a realisation of (52) can be obtained
starting from a representation of the $\mathfrak{su}(2)$ algebra on the
lattice, cf. [7]:

$J_{k}^{+}=(2j-\rho_{k})e^{\hat{\rho}_{k}}\>,\qquad
J_{k}^{-}=\rho_{k}\,e^{-\hat{\rho}_{k}}\>,\qquad J_{k}^{0}=-(\rho_{k}+j)\,,$
(53)

where $\hat{\rho}_{k}=\frac{\partial}{\partial\rho_{k}}$. Redefining
$\rho_{k}\rightarrow\frac{\rho_{k}}{2j}$, so that the commutator is
$\frac{1}{2j}$, and going to the continuous limit, in terms of a density
variable $\rho(x)$ and its conjugate field $\hat{\rho}(x)$, then (53) becomes:

${\cal J}^{+}=2j(1-\rho(x))e^{\hat{\rho}(x)}\>,\qquad{\cal
J}^{-}=2j\rho(x)\,e^{-\hat{\rho}(x)}\>,\qquad{\cal J}^{0}=-2j(\rho(x)+1/2)\,.$
(54)

We obtain the small-noise, hydrodynamic limit as a ‘semiclassical’ one in the
usual way. The transition probability between $\rho(x,t_{1})$ and
$\rho(x,t_{2})$ reads (see [19, 6, 49])

$\displaystyle P(\rho(x,t),\rho(x,t^{\prime}))$ $\displaystyle=$
$\displaystyle\int{\mathcal{D}}{\hat{\rho}}{\mathcal{D}}\rho\,e^{-2j\,N\,S[\rho,\hat{\rho}]}\,,$
(55) $\displaystyle S[\rho,\hat{\rho}]$ $\displaystyle=$
$\displaystyle\int_{t_{1}}^{t_{2}}\int_{0}^{1}dxdt\big{\\{}\hat{\rho}\partial_{t}\rho-{\mathcal{H}}[\rho,\hat{\rho}]\big{\\}}\,.$
(56)

This is a path integral with boundaries $\rho(x,t_{1})$ and $\rho(x,t_{2})$
with ‘classical’ Hamiltonian

${\mathcal{H}}[\rho,\hat{\rho}]=\frac{1}{2}\sigma(\rho)(\nabla\hat{\rho})^{2}-\frac{1}{2}\nabla\rho\nabla\hat{\rho}\,,$
(57)

where $\sigma(\rho)=\rho(1-\rho)$. The fields $\rho$ are constrained, in the
hydrodynamic limit, to satisfy the spatio-temporal boundary conditions

$\displaystyle\forall t,$
$\displaystyle\quad\rho(0,t)=\rho_{0},\quad\rho(1,t)=\rho_{1}\,.$ (58)

The quantum model suggests that, at the fluctuating hydrodynamics level, the
mapping (29) between non-equilibrium and equilibrium is replaced by

$e^{-\\{{\cal J}^{+}_{tot},\;\\}}e^{(1-\bar{\rho})\\{{\cal
J}^{-}_{tot},\;\\}}e^{-\\{{\cal A},\;\\}}e^{-(1-\bar{\rho})\\{{\cal
J}^{-}_{tot},\;\\}}e^{\\{{\cal J}^{+}_{tot},\;\\}}{\cal H}={\cal H}_{eq}\,,$
(59)

where ${\cal H}_{eq}$ coincides in form with ${\cal H}$ but the field
$\rho^{\prime}$ satisfies the boundary conditions

$\displaystyle\forall t,$
$\displaystyle\quad\rho^{\prime}(0,t)=\rho^{\prime}(1,t)=\bar{\rho}.$ (60)

This is a canonical transformation that may be seen as mapping
$(\rho(x),\hat{\rho}(x))$ into another pair
$(\rho^{\prime}(x),\hat{\rho}^{\prime}(x))$, but not a simple contact
transformation mapping $\rho(x)$ into another pair $\rho^{\prime}(x)$.

Now, in principle ${\cal A}$ can be obtained with the classical perturbation
equations:

$\displaystyle\mathcal{B}+\\{{\cal H}_{0},{\cal A}^{(1)}\\}$ $\displaystyle=$
$\displaystyle 0$ $\displaystyle\frac{1}{2}\\{\mathcal{B},{\cal
A}^{(1)}\\}+\\{{\cal H}_{0},{\cal A}^{(2)}\\}$ $\displaystyle=$ $\displaystyle
0$ $\displaystyle\vdots$

Here $\mathcal{A},\mathcal{H}_{0}$ and $\mathcal{B}$ denote the
’semiclassical’ counterparts of $A,H_{0}$ and $B$ introduced in the
microscopic model.

One has to solve these equations under the assumption that ${\cal A}^{(r)}$
are functions of the spins, or of $\rho(x),\hat{\rho}(x)$. This will have a
solution to the extent that ${\cal H}_{0}$ and ${\cal H}_{\cal B}$ are
classically integrable. In that case, one may take advantage of the fact that
${\cal H}_{0}$ may be written in action-angle variables. It turns out that the
fluctuating hydrodynamic limit equations for all the models considered here
are classical integrable models. One may understand this by arguing that
several models (e.g. all the partial exclusion models) share a common
hydrodynamic limit with an integrable one (in this case the SEP). The same is
true, for example, with the KMP model, which is not integrable, but shares the
same hydrodynamic limit as the model of (13). In models that are not
classically integrable, the classical perturbation would be applied on a
trajectory – the instanton in this case.

##  6 Discussion

The results in this paper have to be seen in the light of the work of Graham
[50], who long time ago observed, in the case of Langevin dynamics, that the
distinction between processes that do satisfy reversibility and processes that
do not is somewhat artificial. Our results are in this vein, although here the
mapping is of the system driven out of equilibrium to the same system with
equilibrium baths.

As mentioned above, the transformation mapping the driven problem into the
undriven one is not a function between particle occupations of the two
respective systems, but it involves the whole set including the auxiliary
‘spin’ variables. For example, in a problem where occupation numbers are
associated with the spins $s^{0}_{i}$, the transformation is not only non-
local in the site $i$ of the form
$\\{s^{0}_{i}\\}\rightarrow\\{{s^{\prime}}^{0}_{j}\\}$, but it involves all
the spin operators on a site: in other words, the mapped occupation numbers
are a function, not only of the original occupation numbers, but also of the
current operators on each link. This explains an apparent paradox about
Onsager’s principle: in the driven problem closed trajectories form a loop,
while in the undriven one they do not. A pure ‘contact’ transformation can
only transform loops into loops.

The mapping $P$ in (26) allows one to obtain the non-equilibrium steady state
by transforming the equilibrium measure. At the hydrodynamic level, the
mapping reproduces the non-equilibrium free energy found as the density large
deviation function in MFT [6, 7, 51].

A recurrent question in driven out of equilibrium problems concerns the
existence of an ‘out of equilibrium free energy’. Here we have, for the driven
system, a distinguished candidate: it is the (usual) free energy of the mapped
system. This may in principle be written as a function of the old density
fields and, from the discussion above, their time derivatives. All in all, we
believe that this angle of study of these problems provides an excellent
framework in which one may discuss what properties one can expect, and which
ones one cannot expect, from problems driven out of equilibrium.

## Acknowledgments

We thank the anonymous referees for their comments.

#### Funding information

J.K. is supported by the Simons Foundation Grant No 454943. R.F. is supported
by the German research foundation (Deutsche Forschungsgemeinschaft – DFG)
Research Fellowships Programme 41652715.

## References

  * [1]
  * [2] M. Kac, _Foundations of kinetic theory_ , In _Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability, Volume 3: Contributions to Astronomy and Physics_ , pp. 171–197. University of California Press, Berkeley, Calif. (1956).
  * [3] C. Kipnis, C. Marchioro and E. Presutti, Heat flow in an exactly solvable model, J. Stat. Phys. 27, 65 (1982), 10.1007/BF01011740.
  * [4] B. Derrida, M. R. Evans, V. Hakim and V. Pasquier, Exact solution of a 1D asymmetric exclusion model using a matrix formulation, J. Phys. A: Math. Gen. 26, 1493 (1993), 10.1088/0305-4470/26/7/011.
  * [5] B. Derrida, Non-equilibrium steady states: fluctuations and large deviations of the density and of the current, J. Stat. Mech. P07023 (2007), 10.1088/1742-5468/2007/07/P07023.
  * [6] L. Bertini, A. De Sole, D. Gabrielli, G. Jona-Lasinio and C. Landim, Macroscopic fluctuation theory, Rev. Mod. Phys. 87, 593 (2015), 10.1103/RevModPhys.87.593.
  * [7] J. Tailleur, J. Kurchan and V. Lecomte, Mapping out-of-equilibrium into equilibrium in one-dimensional transport models, J. Phys. A: Math. Theor. 41, 505001 (2008), 10.1088/1751-8113/41/50/505001.
  * [8] L. N. Lipatov, High energy asymptotics of multi–colour QCD and exactly solvable lattice models (1993), arXiv:hep-th/9311037.
  * [9] L. D. Faddeev and G. P. Korchemsky, High energy QCD as a completely integrable model, Phys. Lett. B 342, 311 (1995), 10.1016/0370-2693(94)01363-H.
  * [10] G. P. Korchemsky, Bethe ansatz for QCD pomeron, Nucl. Phys. B 443, 255 (1995), 10.1016/0550-3213(95)00099-E.
  * [11] V. Braun, S. Derkachov and A. Manashov, Integrability of three-particle evolution equations in QCD, Phys. Rev. Lett. 81, 2020 (1998), 10.1103/PhysRevLett.81.2020.
  * [12] N. Beisert et al., Review of AdS/CFT integrability: An overview, Lett. Math. Phys. 99, 3 (2012), 10.1007/s11005-011-0529-2.
  * [13] T. M. Liggett, Interacting particle systems, Springer Berlin Heidelberg, ISBN 9783540226178 (2005), 10.1007/b138374.
  * [14] R. Frassek, Eigenstates of triangularisable open XXX spin chains and closed-form solutions for the steady state of the open SSEP, J. Stat. Mech. 053104 (2020), 10.1088/1742-5468/ab7af3.
  * [15] R. Frassek, C. Giardina’ and J. Kurchan, Duality in quantum transport models (2020), arXiv:2008.03476.
  * [16] D. Bernard and T. Jin, Open quantum symmetric simple exclusion process, Phys. Rev. Lett. 123, 080601 (2019), 10.1103/PhysRevLett.123.080601.
  * [17] G. Schütz and S. Sandow, Non-Abelian symmetries of stochastic processes: Derivation of correlation functions for random-vertex models and disordered-interacting-particle systems, Phys. Rev. E 49, 2726 (1994), 10.1103/PhysRevE.49.2726.
  * [18] J. Keisling, _An ergodic theorem for the symmetric generalized exclusion process_ , Markov Process. Related Fields 4, 351 (1998).
  * [19] H. Spohn, Long range correlations for stochastic lattice gases in a non-equilibrium steady state, J. Phys. A: Math. Gen. 16, 4275 (1983), 10.1088/0305-4470/16/18/029.
  * [20] G. Carinci, C. Giardinà, C. Giberti and F. Redig, Duality for stochastic models of transport, J. Stat. Phys. 152, 657 (2013), 10.1007/s10955-013-0786-9.
  * [21] C. Giardinà, J. Kurchan, F. Redig and K. Vafayi, Duality and hidden symmetries in interacting particle systems, J. Stat. Phys. 135, 25 (2009), 10.1007/s10955-009-9716-2.
  * [22] C. Giardinà, J. Kurchan and F. Redig, Duality and exact correlations for a model of heat conduction, J. Math. Phys. 48, 033301 (2007), 10.1063/1.2711373.
  * [23] E. K. Sklyanin, Boundary conditions for integrable quantum systems, J. Phys. A: Math. Gen. 21, 2375 (1988), 10.1088/0305-4470/21/10/015.
  * [24] R. Frassek, C. Giardinà and J. Kurchan, Non-compact quantum spin chains as integrable stochastic particle processes, J. Stat. Phys. 180, 135 (2020), 10.1007/s10955-019-02375-4.
  * [25] C. Giardinà and J. Kurchan, The Fourier law in a momentum-conserving chain, J. Stat. Mech. P05009 (2005), 10.1088/1742-5468/2005/05/P05009.
  * [26] C. Giardinà, F. Redig and K. Vafayi, Correlation inequalities for interacting particle systems with duality, J. Stat. Phys. 141, 242 (2010), 10.1007/s10955-010-0055-0.
  * [27] S. Floreani, F. Redig and F. Sau, Hydrodynamics for the partial exclusion process in random environment (2019), arXiv:1911.12564.
  * [28] J. Kuan, An algebraic construction of duality functions for the stochastic $U_{q}(A_{n}^{(1)})$ vertex model and its degenerations, Commun. Math. Phys. 359, 121 (2018), 10.1007/s00220-018-3108-x.
  * [29] M. Vanicat, Exact solution to integrable open multi-species SSEP and macroscopic fluctuation theory, J. Stat. Phys. 166, 1129 (2017), 10.1007/s10955-016-1705-7.
  * [30] G. Carinci, C. Giardinà and F. Redig, _The algebraic approach to duality of Markov processes_ , Monograph in preparation (2020) .
  * [31] C. Bernardin and S. Olla, Fourier’s law for a microscopic model of heat conduction, J. Stat. Phys. 121, 271 (2005), 10.1007/s10955-005-7578-9.
  * [32] L. D. Faddeev, How algebraic Bethe ansatz works for integrable model (1996), arXiv:hep-th/9605187.
  * [33] T. Sasamoto and M. Wadati, One-dimensional asymmetric diffusion model without exclusion, Phys. Rev. E 58, 4181 (1998), 10.1103/PhysRevE.58.4181.
  * [34] A. M. Povolotsky, On the integrability of zero-range chipping models with factorized steady states, J. Phys. A: Math. Theor. 46, 465205 (2013), 10.1088/1751-8113/46/46/465205.
  * [35] G. Barraquand and I. Corwin, The $q$-Hahn asymmetric exclusion process, Ann. Appl. Probab. 26, 2304 (2016), 10.1214/15-AAP1148.
  * [36] A. Kuniba, V. V. Mangazeev, S. Maruyama and M. Okado, Stochastic R matrix for $U_{q}(A_{n}^{(1)})$, Nucl. Phys. B 913, 248 (2016), 10.1016/j.nuclphysb.2016.09.016.
  * [37] G. Carinci, C. Giardinà, F. Redig and T. Sasamoto, A generalized asymmetric exclusion process with $U_{q}(\mathfrak{sl}_{2})$ stochastic duality, Probab. Theory Relat. Fields 166, 887 (2016), 10.1007/s00440-015-0674-0.
  * [38] F. C. Alcaraz, M. Droz, M. Henkel and V. Rittenberg, Reaction-diffusion processes, critical dynamics, and quantum chains, Ann. Phys. 230, 250 (1994), 10.1006/aphy.1994.1026.
  * [39] M. Gaudin, _La fonction d’onde de Bethe_ , Commissariat à l’Energie Atomique: série scientifique. Masson, Paris (1983).
  * [40] F. C. Alcaraz, M. N. Barber, M. T. Batchelor, R. J. Baxter and G. R. W. Quispel, Surface exponents of the quantum XXZ, Ashkin-Teller and Potts models, J. Phys. A: Math. Gen. 20, 6397 (1987), 10.1088/0305-4470/20/18/038.
  * [41] C. S. Melo, G. A. P. Ribeiro and M. J. Martins, Bethe ansatz for the XXX-S chain with non-diagonal open boundaries, Nucl. Phys. B 711, 565 (2005), 10.1016/j.nuclphysb.2004.12.008.
  * [42] J. de Gier and F. H. L. Essler, Exact spectral gaps of the asymmetric exclusion process with open boundaries, J. Stat. Mech. P12011 (2006), 10.1088/1742-5468/2006/12/P12011.
  * [43] N. Crampe, E. Ragoucy and M. Vanicat, Integrable approach to simple exclusion processes with boundaries. Review and progress, J. Stat. Mech. P11032 (2014), 10.1088/1742-5468/2014/11/P11032.
  * [44] _The Online-Encyclopedia of Integer Sequences_ , A027614, https://oeis.org.
  * [45] R. Chetrite and S. Gupta, Two refreshing views of fluctuation theorems through kinematics elements and exponential martingale, J. Stat. Phys. 143, 543 (2011), 10.1007/s10955-011-0184-0.
  * [46] J. Haegeman, J. I. Cirac, T. J. Osborne and F. Verstraete, Calculus of continuous matrix product states, Phys. Rev. B 88, 085118 (2013), 10.1103/PhysRevB.88.085118.
  * [47] R. Chetrite and H. Touchette, Nonequilibrium Markov processes conditioned on large deviations, Ann. Henri Poincaré 16, 2005 (2015), 10.1007/s00023-014-0375-8.
  * [48] M. Kruczenski, Spin chains and string theory, Phys. Rev. Lett. 93, 161602 (2004), 10.1103/PhysRevLett.93.161602.
  * [49] V. Lecomte, C. Appert-Rolland and F. van Wijland, Thermodynamic formalism for systems with Markov dynamics, J. Stat. Phys. 127, 51 (2007), 10.1007/s10955-006-9254-0.
  * [50] R. Graham, Solution of Fokker Planck equations with and without manifest detailed balance, Z. Physik B - Condensed Matter 40, 149 (1980), 10.1007/BF01295084.
  * [51] B. Derrida, J. L. Lebowitz and E. R. Speer, Free energy functional for nonequilibrium systems: An exactly solvable case, Phys. Rev. Lett. 87, 150601 (2001), 10.1103/PhysRevLett.87.150601.


11institutetext: University of Leeds, Leeds LS2 9JT, UK 11email:
<EMAIL_ADDRESS>22institutetext: University College
London, London WC1E 6BT, UK
22email<EMAIL_ADDRESS>
# Investigating Guiding Information for Adaptive Collocation Point Sampling in
PINNs

Jose Florido(✉) 11 0009-0001-7291-5029 He Wang 22 0000-0002-2281-5679 Amirul
Khan 11 0000-0002-7521-5458 Peter K. Jimack 11 0000-0001-9463-7595

###### Abstract

Physics-informed neural networks (PINNs) provide a means of obtaining
approximate solutions of partial differential equations and systems through
the minimisation of an objective function which includes the evaluation of a
residual function at a set of collocation points within the domain. The
quality of a PINNs solution depends upon numerous parameters, including the
number and distribution of these collocation points. In this paper we consider
a number of strategies for selecting these points and investigate their impact
on the overall accuracy of the method. In particular, we suggest that no
single approach is likely to be “optimal” but we show how a number of
important metrics can have an impact in improving the quality of the results
obtained when using a fixed number of residual evaluations. We illustrate
these approaches through the use of two benchmark test problems: Burgers’
equation and the Allen-Cahn equation.

###### Keywords:

Partial differential equations Deep learning Physics-informed neural networks
Adaptivity.

##  1 Introduction

###  1.1 Context

This paper is concerned with the mechanism by which Physics-informed neural
networks (PINNs) [1, 2], allow us to introduce our prior knowledge of physics
into a machine learning (ML) algorithm. PINNs are an ML approach to the
solution of systems of partial differential equations (PDEs), where there is
limited or noisy ground truth data. In this context, the PDEs provide a
mathematical model of the physics, which is captured in the method through the
inclusion of loss terms that are evaluated by calculating PDE residuals at
collocation points throughout the domain. This point-based, meshless approach
offers greater flexibility compared to meshed methods in traditional numerical
models - but also poses an interesting problem regarding the number and
distribution of said collocation points.

The original implementation of PINNs uniformly distributes a fixed number of
collocation points randomly throughout the domain. However, in certain
problems some areas are intrinsically harder to learn than others - and this
purely random approach can lead to slow or inefficient training. Different
collocation points contain different amounts of information for learning,
which indicates there is an underlying distribution of locations for
collocation points that can maximise the learning - which is mostly likely
non-uniform and problem dependent. Such a proposition suggests that adjusting
the distribution of points, by biasing towards features of interest, may be
beneficial in terms of efficiency. It may also be necessary for practical
reasons to minimise the risk of the network getting trapped in local minima
during training [3]. However, manually choosing this distribution of
collocation points can be an arduous task and requires a priori knowledge of
the solution, or at least of which areas of the domain will be “most
important”. The alternative is to increase the number of points used globally
which, whilst potentially effective, adds significantly to the computational
cost.

To find the optimal or a good distribution of the collocation points, there
are two common approaches to explicitly or implicitly approximate this
distribution. The first approach is to have fixed collocation points but
weight them differently, i.e. adaptive weighting. Examples of this are [4],
[5]’s Self-Adaptive PINNs, a concept extended by [6]’s DASA point weighing,
and more recently [7]’s Loss-Attentional PINNs. The residual-based attention
scheme in [8] works similarly to the above, weighting the influence of
specific collocation points in the domain to ensure key collocation points
aren’t overlooked.

The other approach is to refine the locations of collocation points, i.e.
adaptive resampling. [9] first presents this type of adaptive refinement based
on residual information. Other sampling approaches include [10, 11, 12].
[13]’s formulation for resampling is a more general version of some of the
previous approaches, and systematically compares adaptive collocation
resampling methods to fixed approaches. Also noteworthy is [14]’s
implementation of a cosine-annealing strategy for restarting training from a
uniform distribution when optimisation stalls. Whilst the above focus on the
collocation points where the PDEs are evaluated, [15] also considers
optimising the selection of experimental points (where the data is available
for supervised learning problems) at the same time.

From the sampling point of view, existing methods essentially seek the ideal
distribution of collocation points, via either explicit parameterisation [9],
or implicit approximation [10]. The information exploited up to date is mostly
focused on the loss function [10, 12, 13], e.g. information such as local PDE
residual at collocation points, or gradients of the loss term. Deviating from
existing literature, we investigate into a different category of information
that proves to be useful in parametrising the ideal collocation points. This
information is the geometric information of the estimated solution, e.g. their
spatial and temporal derivatives. The intuition behind is the solution
geometry reveals intrinsic information about the PDE. For instance, stiffness
regions are harder to learn [7]. By introducing the solution geometry, our aim
is to shed a deeper light into the complexity behind seeking an “ideal”
distribution, and into the interplay between the effectiveness of different
point sampling methods, the complexity of the problem and the computational
cost.

###  1.2 PINNs

The process of training a PINN is similar to that for a regular NN: once a
topology and a set of weights is defined we seek to find values of these
weights that minimise a prescribed loss function. As discussed below, we will
consider a fully unsupervised learning task for which the loss function is
obtained from the sum of all of the residuals of the PDEs at all of the
selected collocation points. The significance of the number and location of
the collocation points is substantial therefore, since they directly impact on
the loss function which is to be minimised. The crux of adaptive resampling is
that a fixed number of points are moved based upon some criteria (that we
investigate), with a view to this modified loss function being a better
objective for the overall minimisation . Previous investigations of this
approach have redistributed the collocation points based upon a probability
density function (PDF) that was formed using the values of the residuals at
the current points or the derivatives of the loss function with respect to the
location of the collocation points. In this work, we consider alternative
criteria based upon mixed derivatives of the estimated solution and of the
residual.

The loss function $\mathcal{L}(\theta)$ of a PINN can be typically
characterised as the sum of the individual losses due to data fitting
($\mathcal{L}_{data}$), enforcing the governing PDE ($\mathcal{L}_{PDE}$), and
weakly enforcing boundary conditions ($\mathcal{L}_{BCs}$) as follows:

$\mathcal{L}(\theta)=w_{data}\mathcal{L}(\theta)_{data}+w_{PDE}\mathcal{L}(\theta)_{PDE}+w_{BC}\mathcal{L}(\theta)_{BCs}$

Here, the relative impact of each loss term can be controlled through the
weight, $w$, whilst $\theta$ represents the trainable network parameters. For
the one-dimensionsal time-dependent problems that we look at, we will consider
a spatial dimension $x$, a temporal dimension $t$ and solve for the dependent
variable $u(x,t)$. The PDE loss term is obtained by evaluating the PDE
residual ($f$, say) for every collocation point in the set of selected
collocation points $\mathrm{x}\in\mathcal{T}$:

$\mathcal{L}_{PDE}=\frac{1}{|\mathcal{T}|}\sum_{\mathrm{x}\in\mathcal{T}}|f(\mathrm{x};u;u_{x},u_{t},u_{xx};\theta)|^{2}$

It is important to distribute these points throughout the domain however
evidence suggests that it may also be beneficial to identify and place a
greater proportion of points in particular regions of importance, so that the
loss term especially includes these areas. This can ensure that optimising the
cost function improves how well the PDE is satisfied around these locations,
and therefore hopefully improves the solution in those areas.

For forward problems, PINNs can be trained without the need for data, and
boundary conditions can be enforced in a ‘hard’ manner. This is the approach
that we take here in order to focus exclusively on the PDE residual
contribution to the loss. We strongly impose the boundary conditions by
applying an output transformation to the results of the neural network and we
treat the PINN as an unsupervised problem with no input data. As a result,
$\mathcal{L}_{PDE}$ becomes the only term to optimise in our loss function.

Throughout this paper, for consistency with other studies in the literature
[13], a simple feedforward NN is used consisting of 3 intermediate layers of
64 nodes each, with a tanh activation function after each layer. The training
consists of 15,000 initial steps using the ADAM optimiser (with learning rate
of 0.001), followed by 1000 steps using L-BFGS [13] before beginning the
resampling process for adaptive methods. For these methods, the points are
then redistributed at this stage and then the training continues with 1000
steps of ADAM and 1000 steps of L-BFGS, repeating until the number of
resamples specified are completed.

##  2 Adaptive resampling

###  2.1 Probability density functions

As mentioned in the introduction, for many problems PINNs will give more
accurate results when the collocation points are distributed following an
appropriate, problem-dependent distribution. The process of rearranging the
collocation points can be automated by using information gathered during the
training process. This approach has been shown to work well in the literature
using the local residual as a guiding metric [9, 13]; but this paper also
proposes alternatives based on the spatial and temporal derivatives of both
the residuals and estimates of the solution. These proposed alternatives are
assessed against both uniform, fixed distributions and existing adaptive
redistribution methods, and shown to be beneficial in many cases.

The method that we use to resample the points is consistent with [13] for
comparison, which chooses the next set of collocation points $\mathcal{T}$
from a fine grid of random points $\mathcal{X}$. All points in $\mathcal{X}$
are assigned a probability ($P(\mathcal{X})$), and the prescribed number of
collocation points are chosen according to the normalised PDF
$\hat{P}(\mathcal{X})$. These are obtained from equations 1 and 2:

$P(\mathcal{X})=\frac{{Y(\mathcal{X})}^{k}}{\overline{Y(\mathcal{X})^{k}}}+c$
(1) $\hat{P}(\mathcal{X})=\frac{P(\mathcal{X})}{\|P(\mathcal{X})\|_{1}}$ (2)

Here $Y(\mathcal{X})$ is the vector containing the magnitude of the
information source chosen (e.g. the current residual or its derivatives, or
some derivatives of the current estimate of the solution),
$\overline{Y(\mathcal{X})^{k}}$ is the mean value across all points in
$\mathcal{X}$ and $k,c$ are constant hyper-parameters that affect the
resampling behaviour. A visualisation of the distribution of points obtained
using the PDE residual as $Y(\mathcal{X})$ for solving the Burgers’ Equation
is shown in figure 1

Figure 1: $Y(\mathcal{X})=$ PDE Residual shown at every point in $\mathcal{X}$
(left); and the corresponding points selected for $\mathcal{T}$ (right).

As a baseline, we also compare the performance of the adaptive methods that we
consider against two fixed distributions: a uniform random distribution and an
alternative uniform, pseudo-random, distribution made with a Hammersley
sampling algorithm [16]. For these fixed distributions, instead of resampling
after a given number of steps, the training is continued until an equivalent
amount of training steps are carried out.

We consider two cases for the adaptive method for which $Y$ is based upon the
current local residuals: one where the initial distribution of points before
any resamples is randomly generated ($PDE,R$), and the other where the initial
distribution follows Hammersley sampling ($PDE,H$). We also consider the use
of the current solution estimates as an information source, defining $Y$ based
upon the mixed (spatial and temporal) second derivative of the current
solution. We will refer to this as the local geometric curvature of $u$, and
denote this case as “$U_{xt}$” in the subsequent text. We also consider a
version of $Y$ that is based upon the mixed second derivative of the local
residual values, which will subsequently be denoted as “$PDE_{xt}$”.

###  2.2 Problem Definition

The main benchmark that we initially use in order to assess the different
sampling strategies that we consider is the 1D Burgers’ Equation, given by:

$uu_{x}+u_{t}=\nu\;u_{xx},\quad x\in[-1,1],\quad t\in[0,1],$ (3)

where the magnitude of $\nu$ determines the relative effect of diffusion.
Dirichlet boundary conditions are assigned:

$u(-1,t)=u(1,t)=0,$

as is the following initial condition:

$u(x,0)=-\sin(\pi x).$

Note that, since we solve this problem on a space-time domain (see figure 2)
the initial condition and the boundary conditions are treated identically by
the PINN. To illustrate the main features of the analytical solution to this
problem a colour map of $u$ is plotted in figure 2. Note that the smooth
initial solution (the left boundary of the domain shown) steepens to a very
sharp front by $t=0.25$ and this remains an important feature of the solution
for all subsequent times.

Figure 2: Burgers’ Equation contour plot of example solution with
$\nu=\frac{\pi}{100}$

As discussed in the previous section the boundary and initial conditions may
be imposed either weakly or strongly when using PINNs. Here we opt for the
latter by applying an output transformation to the $u$ that is output by the
network:

$output=-\sin(\pi*x)+(1-x^{2})ut$

For this case; the PDE residual when used as an information source is
therefore given by:

$PDE=uu_{x}+u_{t}-\nu\;u_{xx}$ (4)

Similarly,

$U_{xt}=\frac{\partial^{2}}{\partial x\partial t}u(x,t)$ (5)

and

$PDE_{xt}=\frac{\partial^{2}}{\partial x\partial t}(uu_{x}+u_{t}-\nu\;u_{xx})$
(6)

To measure and compare the accuracy of the different methods considered an
$L^{2}$ relative error metric is used at the end of training. This compares
the prediction $u$ to ground truth $u_{gt}$ as follows:

$L^{2}=\frac{\sqrt{\sum(u(i)-u_{gt}(i))^{2}}}{\sqrt{\sum u_{gt}(i)^{2}}}\;.$

Note that the sums in this last expression are over every point $i$ in a very
fine background grid which contains many more points than the number of
collocation points used to evaluate $u$. Furthermore, due to the stochastic
nature of ML training, a minimum of 20 repeats with different seeds were
carried out for each case investigated throughout this paper; with the quoted
$L^{2}$ errors always representing the average of these.

###  2.3 Results: default case

For this first test problem we have applied six different solution strategies
using $N=2000$ collocation points. In each case we follow the training regime
of [13], described in section 1.2 (15,000 initial steps of ADAM followed by
1000 steps of L-BFGS and then resampling every 2000 steps (1000 ADAM/1000
L-BFGS). Figure 3 plots the error against the number of resamples taken for
each of the approaches considered, including the two constant point sets based
upon a uniform random distribution and a distribution based upon the
Hammersley sampling algorithm.

Figure 3: Error against number of resamples solving Burgers’ Equation with
described default parameters $v=\frac{\pi}{100};N=2000$

The initial distribution is seen to have an effect in both the fixed and
residual-based methods, with the use of the Hammersley sampling algorithm
resulting in consistently lower error. The local residual method with random
initialisation ($PDE,R$) exhibited the least robustness and greatest
variability, which can be seen from the fact increasing the number of
resamples did not consistently lower the average error between runs.
Generally, more resamples does result in improving accuracy, with error
descending below 0.1% for above 30 resamples for the most of the adaptive
resampling methods.

For the same case we also look at the error achieved using different numbers
of collocation points. This was obtained for 100 rounds of resampling and is
plotted in figure 4. This figure suggests that the threshold for the number of
points required to obtain results of a given accuracy varies from method to
method. Even for the most accurate method ($PDE_{xt}$), the error seems to
plateau however between 0.01% and 0.1%, which the other adaptive methods
eventually reach at the default 2000 points. Nevertheless, this clearly
suggests that the $PDE_{xt}$ approach could allow similar accuracy to other
techniques but at a significantly lower cost.

Figure 4: Error against number of collocation points for different resampling
methods solving Burgers’ Equation with described default parameters
$v=\frac{\pi}{100};$ 100 resamples

In conclusion, based on this single test case, we observe that the achieved
error of a PINN is heavily influenced by the chosen collocation point
distribution, echoing what has been seen previously (see Introduction).
Furthermore, the adaptive resampling methods easily outperform the fixed
distributions, with both of the derivative-based redistribution methods
performing quite well (especially when considering the mixed second derivative
of the residual). In the next section we explore the extent to which this
behaviour generalises by changing multiple aspects of the problem under
consideration, including the initial conditions and the magnitude of the
diffusion term, and by considering a second benchmark PDE: the Allen-Cahn
Equation.

##  3 Results: other cases

###  3.1 Alternate Initial Conditions

In this example we contrast the same methods as in the previous section to
once again solve Burgers’ Equation, but this time we vary the initial
conditions, and therefore the entire evolution of the solution of the PDE. For
each of the four cases considered we compute a high resolution numerical
solution to represent the ground truth and use this to assess the error
against an increasing number of adaptive resamples of the collocation points.
The corresponding results are plotted in figure 5.

Figure 5: Error versus number of resamples. Initial conditions 2 and 3 are
combinations of randomly weighted sin curves of different frequency; with 2
being made up of three low frequency curves and 3 being a higher frequency
combination of four curves. Initial conditions 4 and 5 are $sin(2\pi)$ and
$1.5sin(\pi)$ respectively.

Inspection of figure 5 suggests that for three of the four initial conditions
(2, 4 and 5), $PDE_{xt}$ clearly outperforms the other methods (as it did for
initial condition 1 in the previous section). Not only does it deliver a
smaller error after the last resample (the maximum amount of training
investigated) but the error curve lies beneath the others at each stage in the
adaptive process. This gap is most noticeable in the 5th case, which
corresponds to an initial condition with a larger amplitude, thus leading to a
faster shock formation and a sharper shock. In this case it far outperforms
the next best method. Furthermore, in each case other than initial condition
3, the adaptive collocation point algorithms all performed significantly
better than using fixed collocation points. Unsurprisingly, of the two fixed
point approaches Hammersley always outperforms the uniform random
distribution.

The obvious outlier in figure 5 is initial condition 3. Here all of the
methods to quite poorly in terms of reducing the error and least worst result
is obtained using fixed collocation points based upon a Hammersley sampling.
Further investigation is required to better understand this result, which
shows that the PINN performance (or, at least, the ease with which the PINN
can be trained) is highly dependent upon the problem being solved.

###  3.2 Adjusting PDE parameters

For the following example, we investigate adjusting the value of the $\nu$
parameter in the PDE; which affects the magnitude of the diffusion term
$\nu\;u_{xx}$. Lower values increase the shock sharpness, increasing the
complexity of the problem without fundamentally changing the shape of the
solution. Increasing $\nu$ sees the opposite effect. We use values of $\nu$ of
$0.01$ and $0.001$, approximately multiplying and dividing by 3 from the
default $\nu=\frac{\pi}{100}$ used in the previous section. The ground truth
is again computed using high-resolution numerical simulations and the error
when using different information sources to control the adaptive resampling is
again compared against the fixed uniform methods: see figures 6 and 7.

Figure 6: Error against number of resamples, for Burgers’ Equation with
$\nu=0.01$ Figure 7: Error against number of resamples, for Burgers’ Equation
with $\nu=0.001$

For the simpler case, $\nu=0.01$ (figure 6), all adaptive methods behave
similarly, delivering a small error after relatively few resamples. Comparing
the fixed methods, the Hammerlsey distribution is again consistently better
than a random fixed distribution. Furthermore, in this specific example, it
again outperforms the adaptive sampling methods. Unlike case 3 of figure 5
however, the accuracy of all of these PINNs solutions is high. One possible
explanation for the Hammersley approach being the best in this case is that,
unlike the resampling methods, it uses the L-BFGS optimisation algorithm for
the majority of training steps (as opposed to an equal split between ADAM and
L-BFGS for the adaptive methods). Overall, It is clear that by making the
solution of the problem smoother (through increasing $\nu$) all of the methods
do well and therefore use of the simplest approach, based upon a fixed
distribution of collocation points, can be appropriate.

For the more challenging case of $\nu=0.001$ (figure 7), all methods struggle
to converge to an accurate solution using $N=2000$ collocation points. The
least worst approach, with an $L^{2}$ error of 14%, is again the adaptive
method based upon $PDE_{xt}$. Based upon this set of results we can make the
opposite proposition: that for complex problems, the choice of adaptive
sampling method may not be the limiting factor to accuracy, and an increase to
either the number of collocation points or the network complexity may be
needed to increase accuracy.

We test this proposition by undertaking a study to look at the impact that
changing the number of collocation points has on the solutions obtained in the
$\nu=0.01$ and $\nu=0.001$ cases. In table 1 we consider the simpler
$\nu=0.01$ case, observing the change in performance when the number of
collocation points is halved to $1000$. Conversely, in table 2, the more
complex case with the sharper shock ($\nu=0.001$), we double the number of
collocation points to $N=4000$ and observe whether there are any significant
increases in accuracy.

Table 1: Effect of changing $N$; $L^{2}$ Error at 100 resamples for $\nu=0.01$ $N$ | Fixed, R | Fixed, H | $PDE,R$ | $PDE,H$ | $U_{xt}$ | $PDE_{xt}$
---|---|---|---|---|---|---
2000 | 2.019E-04 | 8.95E-05 | 1.024E-04 | 1.152E-04 | 1.089E-04 | 1.055E-04
1000 | 3.089E-03 | 3.925E-04 | 1.503E-04 | 1.867E-04 | 1.458E-04 | 1.258E-04
Error Change | 1429.9% | 338.4% | 46.8% | 62.1% | 33.9% | 19.3%

In table 1, decreasing the number of collocation points from 2000 to 1000
increases the error as expected. However, there is a big discrepancy in how
much the different methods are affected. A positive % signifying increase in
error, the fixed methods are substantially less effective with fewer points.
The opposite is true for adaptive methods, especially for $PDE_{xt}$ which
only sees an increase in error of 19.3% compared to the 3-fold increase of the
fixed Hammersley distribution of points.

Table 2: Effect of changing $N$; $L^{2}$ Error at 100 resamples for $\nu=0.001$ $N$ | Fixed, R | Fixed, H | $PDE,R$ | $PDE,H$ | $U_{xt}$ | $PDE_{xt}$
---|---|---|---|---|---|---
2000 | 4.209E-01 | 5.093E-01 | 3.925E-01 | 4.780E-01 | 3.797E-01 | 1.425E-01
4000 | 4.049E-01 | 3.836E-01 | 1.858E-01 | 1.996E-01 | 2.591E-01 | 3.660E-02
Error Change | -3.9% | -32.8% | -111.2% | -139.5% | -46.6% | -289.3%

In table 2, we are this time looking at the decrease in error as a result of
doubling the number of collocation points in the case of $\nu=0.001$. As
anticipated, the error does decrease in each column (signified by a negative
increase) however this improvement is relatively small in most cases, meaning
that the overall errors are still not competitive, even with this increase in
resources. The biggest decreases are seen for the $PDE,H$ and especially
$PDE_{xt}$, which is the only method to achieve an error under 10%.

###  3.3 An alternative PDE: Allen-Cahn

Having presented several variations of the Burgers’ Equation problem, our next
step is to consider a different PDE entirely and observe the results. The
Allen-Cahn equation is chosen for this due to its complex transient behaviours
which are known to be challenging to capture computationally. The problem is
defined as follows:

$\frac{\partial u}{\partial t}=D\frac{\partial^{2}u}{\partial
x^{2}}+5(u-u^{3}),\quad x\in[-1,1],\quad t\in[0,1]$ (7)

We take D = 0.001 and use the initial condition

$u(x,0)=x^{2}\cos(\pi x)\;,$

and boundary conditions

$u(-1,t)=u(1,t)=-1\;.$

For this problem we again compute a high-resolution numerical solution for the
ground truth and consider the $L^{2}$ error as the number of resamples is
increased. Based on the observations of the previous subsections, we undertake
these tests for three distinct numbers of collocation points:
$N=500,1000,2000$. The same resampling methods as in section 2 are used, using
a minimum of 20 repeats per case. The results for $N=500$ are shown below in
figure 8. For this case, the adaptive methods are again consistently superior
to using fixed collocation points. Similar to the Burger’s equation solutions
with a large $\nu=0.01$, there are no clear advantages between $PDE_{xt}$, and
the approaches based upon the local residual (though using $U_{xt}$ is not so
good).

Figure 8: Error against number of resamples solving the Allen-Cahn Equation
with $N=500$

For $N\geq 1000$, a similar pattern is observed to that shown in figure 6,
where there is little difference between resampling and fixed methods. In this
case, the increased number of points again makes the problem simple enough to
solve with fixed collocation strategies, so that the benefits previously
observed from using adaptive resampling methods disappear. In fact a fixed
pseudo-random Hammersley distribution performs better than the other curves:
perhaps again benefiting from more L-BFGS optimiser steps.

For the Allen-Cahn Equation case, the behaviour and importance of sampling in
these methods varies significantly depending on the number of points used. The
less sharp shocks in this solution could also be influencing the effectiveness
of using derivatives of the solution as redistribution criteria. Nevertheless,
the adaptive methods clearly outperform the fixed distribution of collocation
points for the case $N=500$, again suggesting that results of a prescribed
accuracy can be obtained with fewer collocation points when adaptive sampling
is permitted.

##  4 Discussion and Conclusion

In this paper, to better explore the problem of collocation point sampling in
PINNs we have examined the two benchmark problems of Burgers’ Equation and the
Allen-Cahn Equation. We have used novel information sources based on
derivatives of the local residual and the solution estimate to guide the
resampling process, which has been kept consistent with other approaches in
literature, and compared the results to other resampling methods and fixed
random and pseudo-random distributions.

Overall, we have observed that using adaptive point distribution typically
allows more accurate PINNs solutions to be obtained than using fixed point
distributions for small numbers of collocation points. Using fewer collocation
points is potentially important since the training cost of a PINN grows in
proportion to this number (all other factors being equal). Furthermore, we
have observed that defining the probability distribution for the collocation
point locations based upon the mixed second derivative of the residual with
respect to the independent variables typically outperforms, often
significantly, using the residual alone, as in [13].

For Burgers’ equation with a smooth initial condition that leads to the
formation of a single “shock”, and a diffusion parameter such that this
features is sufficiently sharp, we obtained the smallest errors using this
approach (Fig. 3). Equally importantly, we observed that we could reach such
errors with significantly fewer collocation points than for any of the other
techniques that we compared against (Fig. 4). When considering other initial
conditions and other diffusion parameters the relative performance of the
adaptive approaches was seen to depend upon the complexity of the specific
solution in each case. Typically, there appears to be a value of $N$ above
which the use of a uniform (Hammersley) distribution of points is preferred -
however this value can be relatively large for the more complex problems (and
is relatively small when the solution is smoothly varying everywhere). When
considering numbers of points below this critical value of $N$ we invariably
find that using the mixed derivative of the residual is the best strategy for
defining the PDF for drawing the collocation point locations.

Consideration of a different PDE, the Allen-Cahn equation in this case,
demonstrates similar behaviour. Specifically, for a given problem there is a
value of $N$ below which the adaptive sampling approach is preferable (and
above which a uniform Hammersley distribution suffices). In the particular
examples considered here, defining the PDF based upon the mixed derivative of
the residual is still reliable and robust, however its improvement relative to
using the residual alone is less clear. Future work will extend the
exploration of these resampling methods to other PDEs beyond Allen-Cahn,
helping assess the generalizability of our findings. As our approach does not
explicitly constrain the type of domain we can study, we should also be able
to consider domains with higher spatial dimensions and more complex shapes.

In terms of computational cost, the PINNs network we have considered is not
competitive with standard discretization methods for solving individual
forward problems. This is primarily due to the high computational cost of
training with full gradient descent methods like L_BFGS. As for the cost of
the resampling strategies we have discussed, these only contribute to a small
overhead (estimated at 3% in [13]) to the overall computational cost of PINNs
training.

Considering other common problems in ML, the biggest challenge is determining
what constitutes an appropriate training regime. We did not observe
overfitting to occur, even with fixed sampling methods. However, the choice of
optimisation algorithm and training regime was still impactful. The use of
faster stochastic gradient descent methods like ADAM helped initialise
training but had a limited ceiling to accuracy, hence the need for the more
expensive L_BFGS. Whilst we were not able to fully explore the impact of
different training regimes within the scope of this paper, efficiency could
feasibly be improved by further optimising this via methods such as cosine
annealing [14]; or by optimising the architecture employed for solving
problems of varying complexity.

In summary, we have shown that a number of different strategies are viable for
selecting the location of PINNs collocation points with resampling methods. In
particular, we have shown that resampling can be guided not only by the local
residuals, but also by the spatial and temporal derivatives of the residuals
and even of the solution estimates themselves. Furthermore, we have observed
that a number of different factors impact the accuracy of the PINN method,
making it difficult to analyse any single aspect (such as collocation point
selection) in isolation. Whilst we believe that our contribution begins to
shed some light on the interplay between the location of the collocation
points and the ability of a PINN to learn the PDE solution it is clear that
further research is desirable in order to be able to propose robust
implementations that support reliable training or provide accuracy guarantees.

####  4.0.1 Acknowledgements

The first author gratefully acknowledges financial support via the EPSRC
Centre for Doctoral Training in Fluid Dynamics at Leeds (EP/S022732/1).

####  4.0.2

The authors have no competing interests to declare that are relevant to the
content of this article.

## References

  * [1] Lagaris, I., Likas. A., Fotiadis, D.: Artificial neural Networks for Solving Ordinary and Partial Differential Equations. IEEE Transactions on Neural Networks 9(5), 987-1000 (1998)
  * [2] Raissi, M., Perdikaris, P., Karniadakis, G.: Physics-Informed Neural Networks: A Deep Learning Framework for Solving Forward and Inverse Problems Involving Nonlinear Partial Differential Equations. Journal of Computational Physics 378, 686-707 (2019)
  * [3] Wang, Y., Lai, C., Gomez-Serrano, J., Buckmaster, T.:Asymptotic Self-Similar Blow-Up Profile for Three-Dimensional Axisymmetric Euler Equations Using Neural Networks. Physical Review Letters 130, 244002 (2023)
  * [4] Li, W., Zhang C., Wang, C., et al.: Revisiting PINNs: a Generative Adversarial Physics-Informed Neural Networks and Point-Weighting Method, arXiv Preprint, arXiv:2205.08754 (2022)
  * [5] McClenny, L. Braga-Neto, U.: Self-Adaptive Physics-Informed Neural Networks using a Soft Attention Mechanism. Journal of Computational Physics 474, 111722 (2022)
  * [6] Zhang, G., Yang, H., Zhu, F., et al.: Dasa-pinns: Differentiable adversarial self-adaptive pointwise weighting scheme for physicsinformed neural networks, SSRN (2023)
  * [7] Song, Y., Wang, H., Yang, H.: Loss-Attentional Physics-Informed Neural Networks. Journal of Computational Physics 501, 112781 (2024)
  * [8] Anagnostopoulos, S., Toscano, J., Karniadakis, G., et al.: Residual-based attention and connection to information bottleneck theory in PINNs, arXiv Preprint, arXiv:2307.00379 (2023)
  * [9] Lu, L., Meng, X., Karniadakis, G., et al.: DeepXDE: A Deep Learning Library for Solving Differential Equations. SIAM Review 63, 208-228 (2021)
  * [10] Nabian, M., Gladstone, R., Meidani, H.: Efficient training of physics-informed neural networks via importance sampling. Computer-Aided Civil and Infrastructure Engineering 36(8), 962-977 (2021)
  * [11] Hanna, J., Aguado, J., Borzacchiello, D., et al.: Residual-based adaptivity for two-phase flow simulation in porous media using Physics-informed Neural Networks. Computer Methods in Applied Mechanics and Engineering 396, 115100 (2022)
  * [12] Gao, Z., Yan, L., Zhou, T.: Failure-Informed Adaptive Sampling for PINNs. SIAM Journal on Scientific Computing 45, A1971-A1994 (2023)
  * [13] Wu, C., Zhu, M., Lu, L., et al.: A comprehensive study of non-adaptive and residual-based adaptive sampling for physics-informed neural networks. Computer Methods in Applied Mechanics and Engineering 403, 115671 (2023)
  * [14] Subramanian, S., Kirby, R., Gholami, A., et al/: Adaptive Self-Supervision Algorithms for Physics-Informed Neural Networks. ECAI 2023 372, (2023)
  * [15] Lau, G., Hemachandra, A., Low, B., et al.: PINNACLE: PINN Adaptive ColLocation and Experimental points selection. In: The Twelfth International Conference on Learning Representations (2024)
  * [16] Hammersley, J. M., and D. C. Handscomb. Monte Carlo methods, methuen & co. Ltd., London 40 (1964): 32.