††thanks: Current address: IBM Research - Haifa, Haifa University Campus,
Mount Carmel, Haifa 31905, Israel

# Static kinks in chains of interacting atoms

Haggai Landa<EMAIL_ADDRESS>Institut de Physique Théorique, Université
Paris-Saclay, CEA, CNRS, 91191 Gif-sur-Yvette, France  Cecilia Cormick  IFEG,
CONICET and Universidad Nacional de Córdoba, X5016LAE Córdoba, Argentina
Giovanna Morigi  Department of Physics, Universität des Saarlandes, D-66123
Saarbrücken, Germany

###### Abstract

We theoretically analyse the equation of topological solitons in a chain of
particles interacting via a repulsive power-law potential and confined by a
periodic lattice. Starting from the discrete model, we perform a gradient
expansion and obtain the kink equation in the continuum limit for a power-law
exponent $n\geq 1$. The power-law interaction modifies the sine-Gordon
equation, giving rise to a rescaling of the coefficient multiplying the second
derivative (the kink width) and to an additional integral term. We argue that
the integral term does not affect the local properties of the kink, but it
governs the behaviour at the asymptotics. The kink behaviour at the center is
dominated by a sine-Gordon equation and its width tends to increase with the
power law exponent. When the interaction is the Coulomb repulsion, in
particular, the kink width depends logarithmically on the chain size. We
define an appropriate thermodynamic limit and compare our results with
existing studies performed for infinite chains. Our formalism allows one to
systematically take into account the finite-size effects and also slowly
varying external potentials, such as for instance the curvature in an ion
trap.

††preprint: APS/123-QED

##  I Introduction

The Frenkel-Kontorova model reproduces in one dimension the essential features
of stick-slip motion between two surfaces Frenkel-Kontorova ; Bak:RPP ;
Tosatti:RMP . The ground state is expected to describe the structure of a one-
dimensional crystal monolayer growing on top of a substrate crystal. In one
dimension the elastic crystal is modelled by a periodic chain of classical
particles with uniform equilibrium distance $a$, which interact with a
sinusoidal potential with periodicity $b$ Pokrovski:Book ; Kivshar:Book .
Frustration emerges from the competition between the two characteristic
lengths: Depending on the mismatch between $a$ and $b$ and on the strength of
the substrate potential, a continuous transition occurs between a structure
with the substrate’s period (commensurate) and an incommensurate structure
Talapov . The transition is characterized by proliferation of kinks, namely,
of local distributions of excess particles (or holes) in the substrate
potential. When the interactions of the elastic crystal are nearest-neighbour,
in the long-wavelength limit the dynamics of a single kink is governed by the
integrable sine-Gordon equation Merwe ; Talapov ; Tabor .

The experimental realizations of crystals of interacting atoms, such as ions
Dubin:RMP , dipolar gases Ferlaino , and Rydberg excitons Rydberg:Review ,
offer unique platforms for analysing the Frenkel-Kontorova model dynamics
Garcia-Mata ; Pruttivarasin ; Vanossi . The substrate potential can be
realised by means of optical lattices Linnet:2012 ; Enderlein:2013 ; Ferlaino
; Cetina:2013  or of a second atomic crystal Drewsen ; Mehlstaeubler .
Periodic boundary conditions can be implemented in ring traps Haeffner . Kinks
and dislocations can be imaged Pyka ; Ulm ; Mielenz ; Monroe ; Atom:microscope
; Rydberg:Arimondo  and spectroscopically analysed Brox . Differing from
textbook models, the particles’ interaction is a power-law potential, whose
exponent could be engineered by means of lasers Porras:RPP . Nano-friction
have been experimentally investigated in small ion chains in periodic
potentials Vuletic:1 ; Vuletic:2 ; Vuletic:3 .

The study of kinks and of nanofriction in these systems requires one to
analyse the effects of the tails of the interactions on the sine-Gordon
equation. Specifically, in one dimension the energy is non additive in Coulomb
systems Kac . Yet, long-range effects are marginal: the exponent of the
Coulomb interaction formally separates two regimes, such that for slower
power-law decays the dynamical equations are characterized by fractional
spatial derivatives, while for faster decays the spatial derivatives are of
integer order Zavlasky ; Staffilani . The effect of long-range interactions on
the commensurate-incommensurate transition have been discussed Bak:PRL ;
Pokrovsky:LR , the kink solutions in a periodic potentials have been analysed
numerically for the long-range Kac-Baker interactions Kac-Baker  and for
dipolar and Coulomb potentials Braun:1990 . Analytic studies of the kink
solutions have been performed in the thermodynamic limit Kac-Baker ;
Pokrovsky:LR ; Braun:1990 .

The aim of the present work is to review the analytical derivation of the kink
equation for power-law interacting potentials by means of a gradient
expansion, which is implemented following the lines of the study of Ref.
Morigi:2004 ; Morigi:2004:PRL . This derivation allows one to determine the
local properties of the kink as a function of the interaction range for
integer exponents when this decays with the distance as the Coulomb repulsion
or faster. Moreover, it allows one to determine its asymptotic behaviour, as
well as to systematically take into account the finite-size effects, thus
setting the basis of a study where these effects can be included in a
perturbative fashion.

This manuscript is organised as follows. In Sec. II we introduce the
Lagrangian of an atomic chain in a periodic potential, where the atoms
interact via a power-law potential. In Sec. III we consider the long-
wavelength limit and derive the equation for the static kink. We discuss
separately the case of Coulomb interactions. We analyse then the thermodynamic
limit and compare our results with the ones of Ref. Braun:1990 . The
conclusions are drawn in Sec. IV.

##  II An atomic chain in a periodic substrate potential

A chain of $N$ interacting atoms with mass $m$ is confined in a finite volume
and is parallel to the $x$ axis. Their atomic positions and canonically-
conjugated momenta are $x_{j}$ and $p_{j}$, with $j=-N/2,N/2+1\ldots,N/2-1$
and $x_{j}<x_{j+1}$. Their Lagrangian $\mathcal{L}$ reads

$\mathcal{L}=\sum_{i=1}^{N}\frac{m\dot{x}_{j}^{2}}{2}-V_{\rm pot}\,,$ (1)

where $V_{\rm pot}$ is the potential energy and is thus the sum of a periodic
substrate potential and of the harmonic interaction between pairs of
particles, $V_{\rm pot}=V_{\rm opt}+V_{n}$. The periodic substrate potential
$V_{\rm opt}$ is a sinusoidal lattice with periodicity $b$ and amplitude
$V_{0}$:

$V_{\rm opt}=\sum_{i=1}^{N}V_{0}\left[1-\cos\left(\frac{2\pi
x_{i}}{b}\right)\right]\,,$ (2)

The atomic interaction $V_{n}$ couples atoms at distance $ra$ with strength
scaling as $1/r^{n+2}$:

$V_{n}=\frac{1}{2}\sum_{i}\sum_{r>0}\frac{K_{n}}{r^{n+2}}(x_{i+r}-x_{i}-ra)^{2}\,,$
(3)

Here, $n$ an integer number, $n\geq 1$, and $K_{n}$ is the spring constant
between nearest neighbour. The interaction term vanishes when the atoms are at
the equilibrium positions $x_{j}^{(0)}=ja$.

We note that interactions of this form are obtained by expanding the
interaction potential till second order in the displacements about the
equilibrium positions of the interaction forces, and discarding
anharmonicities. Let the interaction between two particles at distance $x$ be
given by

$W_{\rm int}(x)=W_{n}/x^{n}\,,$ (4)

where $W_{n}$ is a constant which depends on $n$. Then, the spring constant
takes the form

$K_{n}=\partial_{x}^{2}W_{\rm int}(x)|_{x=a}=\frac{n(n+1)W_{n}}{a^{n+2}}\,.$
(5)

The textbook limit, where the particles interact via nearest-neighbour
interactions, is recovered by letting $n\to\infty$. In this case $W_{n}$ shall
be appropriately rescaled in order to warrant that the spring constant
$K_{\infty}$ remains finite, with $K_{\infty}=\lim_{n\to\infty}K_{n}$. In the
following we discuss the cases of $n>1$ (with dipolar and van der Waals
interactions corresponding to the particular values $n=3$ and $n=6$
respectively), as well as the Coulomb interaction $n=1$. For Coulomb
interaction between particles of charge $q$, the spring constant is
$K=2q^{2}/(4\pi\varepsilon_{0}a^{3})$, with $q$ the charge of the particles
and $\varepsilon_{0}$ the vacuum’s permittivity Morigi:2004 .

Here and in what follows we assume motion along one axis and open boundary
conditions. This can be realised by means of anisotropic traps and
sufficiently cold atoms, where one can assume that the motion in the
transverse direction is frozen out. Typically the trap curvature gives rise to
inhomogeneity in the equilibrium particle distribution, which leads to a
position-dependent spring constant $K_{n}$. Below we assume that the atoms are
uniformly distributed. We note, however, that this formalism can
systematically include the trap inhomogeneity, as shown in Ref. Morigi:2004 ,
as long as the trap curvature can be treated in the continuum limit Dubin:1997
.

###  II.1 Equilibrium configuration in the discrete chain

The equilibrium configurations of the discrete chain are solutions of the
equation of motion

$\displaystyle m\ddot{x}_{j}$ $\displaystyle=$
$\displaystyle\sum_{r>0}\frac{K_{n}}{r^{n+2}}\left[x_{j+r}-x_{j}-(x_{j}-x_{j-r})\right]$
$\displaystyle-\frac{2\pi}{b}V_{0}\sin\left(\frac{2\pi x_{i}}{b}\right)\,,$

which satisfy $m\ddot{x}_{j}=0$ for all $j$. For $V_{0}=0$, namely, in the
absence of the substrate potential, the ground state is the uniform chain with
equilibrium positions

$x_{j}^{(0)}=ja\,.$

For $V_{0}\neq 0$ the equilibrium configuration $\\{\bar{x}_{j}^{(0)}\\}$
reads

$\bar{x}_{j}^{(0)}=x_{j}^{(0)}+\bar{u}_{j}=ja+\bar{u}_{j}\,,$

where $\bar{u}_{j}$ are static displacements that solve the set of equations:

$\displaystyle\frac{2\pi
V_{0}}{b}\sin\left[2\pi(x_{i}^{(0)}+\bar{u}_{i})/b)\right]=\sum_{r\neq
0}\frac{K_{n}}{|r|^{n+2}}(\bar{u}_{i+r}-\bar{u}_{i})\,.$ (7)

A static kink describes particles displacements which are localized in a
region of the chain, such that the chain is uniform at the edges. The solution
interpolates between the two boundary values $u_{j}\to 0$ for $j\to-N/2$ and
$u_{j}\to b$ for $j\to+N/2$ and describes a topological soliton. The antikink
is the topological soliton of opposite charge and interpolates between
$u_{j}\to 0$ for $j\to-N/2$ and $u_{j}\to-b$ for $j\to+N/2$.

It is convenient to introduce the phase variable $\theta_{j}$ of particle $j$,
which is a dimensionless scalar and is defined as

$\theta_{j}=2\pi(\bar{x}_{j}^{(0)}-\ell bj)/b\,,$ (8)

or, alternatively:

$\bar{x}_{j}^{(0)}=j\ell b+\frac{b}{2\pi}\theta_{j}\,.$ (9)

The phase function $\theta_{j}$ gives the shift of the ion $j$ from the
commensurate configuration, such that for every $2\pi$ change in the phase the
respective ion position is shifted by one period $b$ of the periodic
potential. The equation of the phase function is then given by

$\displaystyle-\sum_{r\neq
0}\frac{(\theta_{i+r}-\theta_{i}-(\delta_{i+r}-\delta_{i}))}{|r|^{n+2}}+m_{K}^{2}\sin(\theta_{i})=0\,,$
(10)

where

$\delta_{j}=2\pi(x_{j}^{(0)}-\ell bj)/b\,,$ (11)

and represents the mismatch between the former equilibrium positions of the
crystal and the nodes of the periodic potential. The equation depends now on
the power-law exponent, on the mismatches $\delta_{j}$, and more specifically,
on the ratio $a/b$, and on the dimensionless ratio

$m_{K}=2\pi\sqrt{\frac{V_{0}}{Kb^{2}}}\,,$ (12)

which scales the weight of the interactions with the respect to the localizing
potential. We now conveniently rewrite the mismatch by introducing the ratio
$\ell=\lfloor a/b\rceil$, with $\lfloor x\rceil$ the nearest integer to $x$.
The configuration is commensurate when the ratio $a/b=\ell$, then the ground-
state equilibrium positions are $x_{j}^{(0)}=ja=\ell jb$. When instead
$\ell\neq a/b$ the displacement is generally $\bar{u}_{j}\neq 0$. For
$x_{j}^{(0)}=ja$ in Eq. (11), one obtains $\delta_{j}=j\delta$ with

$\delta=2\pi(a-\ell b)/b\,,$ (13)

which vanishes when both length scales are commensurate with each other. In
general, $\delta$ is a periodic function of the ratio $a/b$, in fact
$\delta(a/b+p)=\delta(a/b)$, $p\in\mathbb{Z}$, and its value ranges in the
interval $-\pi<\delta<\pi$. We use this expression and rewrite Eq. (10) in
terms of the phase function:

$\displaystyle-\sum_{r\neq
0}\frac{(\theta_{i+r}-\theta_{i}-r\delta)}{|r|^{n+2}}+m_{K}^{2}\sin(\theta_{i})=0\,.$
(14)

In the rest of this manuscript we will consider structures whose ground state
is commensurate, and thus take $\delta=0$ in Eq. (14). We analyse the equation
of motion of single topological solitons in the long-wavelength limit when the
periodic substrate potential is a small perturbation to the elastic crystal,
$m_{K}\ll 1$.

##  III Equation of the static kink in the long-wavelength limit

In this section we derive the continuum limit of Eq. (10) extending the
procedure of Ref. Morigi:2004  to kink’s equation. In our treatment we keep
$n$ finite and consider the commensurate case, $a=\ell b$, when the mismatch
$\delta$ vanishes. For completeness, we first shortly review the derivation of
the sine-Gordon equation in the case of nearest-neighbour interactions and
briefly discuss the properties of its exact solution.

###  III.1 Nearest-neighbour interactions

The sine-Gordon equation of the Frenkel-Kontorova model is found in our model
by considering the limit $n\to\infty$ in Eq. (14), thus keeping only the
nearest neighbour terms ($r=1$) and taking $K_{\infty}=\lim_{n\to\infty}K_{n}$
to be constant. The long-wavelength limit is determined by first taking the
Fourier transform of the phase function,
$\theta_{j}\to\tilde{\theta}(k)=\sum_{j}{\rm e}^{{\rm
i}kx_{j}}\theta_{j}/\sqrt{N}$ and expanding the equation of motion for small
$k$ Zavlasky . In $k$ space the interaction term takes the form
$-2K_{\infty}(1-\cos(ka))\tilde{\theta}(k)\sim-k^{2}a^{2}K_{\infty}\tilde{\theta}(k)$.
Going back to position space, this term is cast in terms of a second-order
derivative of the phase function. Then, Eq. (14) takes the form of a sine-
Gordon equation Tabor :

$\displaystyle-a^{2}\frac{\partial^{2}}{\partial
x^{2}}\theta(x)+m_{K}^{2}\sin(\theta(x))=0\,.$ (15)

One can rewrite the equation as

$\displaystyle-d^{2}\frac{\partial^{2}}{\partial
x^{2}}\theta(x)+\sin(\theta(x))=0\,,$ (16)

where $d$ is the width of the kink Braun:1990 ,

$d=a\,\sqrt{\frac{K_{\infty}}{V_{0}}\left(\frac{b^{2}}{4\pi^{2}}\right)}=\frac{a}{m_{K}}\,.$
(17)

Equation (16) admits the solution Tabor

$\theta(x)=4\tan^{-1}[e^{\sigma(x-x_{0})/d}]\,,$ (18)

with $\sigma=\pm 1$ the topological charge and $x_{0}$ the position of the
kink. Specifically, for $\sigma=1$, the solution is a single kink with
$\theta(x)\to 2\pi$ for $x\to\infty$, and $\theta(x)\to 0$ for $x\to-\infty$.
The antikink corresponds to $\sigma=-1$ and is thus the mirror reflection of
the kink at $x_{0}$. We note that the continuum limit is consistent when the
width of the kink is much larger than the interparticle distance, $d\gg a$,
and the potential is a small perturbation to the elastic crystal.

The constant $m_{K}$ scales the mass $m_{\rm SG}$ of the sine-Gordon kink. We
recall that the kink’s mass $m_{\rm kink}$ is defined as $m_{\rm
kink}=m\sum_{j}(\partial\bar{u}_{j}/\partial x_{0})^{2}$. In the continuum
limit and using Eq. (18), one obtains Braun:1990

$m_{\rm SG}=\frac{8m_{K}}{a}\,m\left(\frac{b}{2\pi}\right)^{2}\,.$ (19)

Anharmonicities of the potential give rise to higher-order derivatives and
thus to an asymmetry between kink and antikink Willis ; Braun:1990 . In this
treatment we discard these terms, restricting the expansion of the interaction
potential to the harmonic terms. In this limit the kink and antikink in the
continuous limit just differ because of the topological charge $\sigma$. We
further note that, by performing the continuum limit, we discarded higher
order terms in the gradient expansion. These terms account for discreteness
effects and give rise to an effective kink narrowing Braun:1990 .

###  III.2 Power-law interactions

We now perform a gradient expansion for power-law interactions. This is done
by means of a manipulation of Eq. (14) which is equivalent to a Taylor
expansion for low momenta in Fourier space, and consists in singling out the
terms contributing to the second derivative. The equation of motion for the
static kink then becomes an integro-differential equation. In the following we
will assume $m_{K}\ll 1$, unless otherwise specified.

In order to study the long-wavelength limit we use the prescription
$2i/N\to\xi$ and $\theta_{j}\to\theta(\xi)$. For $N\gg 1$, then $\xi$ can be
treated as a continuous variable defined in the interval $[-1/2,+1/2]$. With
this prescription Eq. (14) takes the form:

$-\frac{(d_{n}/a)^{2}}{N^{n+1}}\mathcal{I}_{n}+\sin[\theta(\xi)]=0\,,$ (20)

where $\mathcal{I}_{n}$ is dimensionless:

$\displaystyle\mathcal{I}_{n}=$
$\displaystyle\int_{\bar{a}}^{1/2}d\xi^{\prime}\frac{1}{\xi^{\prime
n+2}}[\theta(\xi+\xi^{\prime})-2\theta(\xi)+\theta(\xi-\xi^{\prime})]\,,$ (21)

and the discrete nature of the chain at atomic distances enters through the
(high-frequency) cutoff

$\bar{a}=1/N\,.$

In Eq. (20) we have also introduced the characteristic length

$d_{n}=a\,\sqrt{\frac{K_{n}}{V_{0}}\left(\frac{b}{2\pi}\right)^{2}}\,,$ (22)

which would correspond to the kink’s width when simply truncating the sum in
Eq. (14) till the nearest-neighbours. We note that in this formalism the chain
length is a low-frequency cutoff, and in particular the limit $N\gg 1$ is
equivalent to small $k$ in Fourier space. We will first keep $N$ constant,
thus consider a finite chain, and take the thermodynamic limit only after
performing the gradient expansion.

We now integrate Eq. (21) by parts applying the procedure as in Ref.
Morigi:2004  and rewrite it as the sum of three terms Morigi:2004 :

$\mathcal{I}_{n}=I_{\rm edge}+I_{\bar{a}}+I_{n}\,,$ (23)

where the first term on the RHS contains the contributions of the chain’s
edges:

$\displaystyle I_{\rm edge}=-$
$\displaystyle\frac{1}{(1+n)}\frac{[\theta(\xi+1/2)+\theta(\xi-1/2)-2\theta(\xi)]}{(1/2)^{n+1}}$
$\displaystyle-\frac{1}{2n(1+n)}\frac{[\theta^{\prime}(\xi+1/2)-\theta^{\prime}(\xi-1/2)]}{(1/2)^{n+1}}\,,$
(24)

and $\theta^{\prime}(\xi)\equiv\partial\theta(\xi)/\partial\xi$. The second
integral describes the contribution of the nearest-neighbour terms:

$\displaystyle I_{\bar{a}}=$
$\displaystyle\frac{1}{(1+n)}\frac{[\theta(\xi+\bar{a})-2\theta(\xi)+\theta(\xi-\bar{a})]}{\bar{a}^{n+1}}$
$\displaystyle+\frac{\bar{a}}{n(1+n)}\frac{[\theta^{\prime}(\xi+\bar{a})-\theta^{\prime}(\xi-\bar{a})]}{\bar{a}^{n+1}}\,.$
(25)

By making a Taylor expansion about $\xi$, $I_{\bar{a}}$ can be cast into the
form :

$\displaystyle I_{\bar{a}}=$
$\displaystyle\frac{\bar{a}^{2}}{\bar{a}^{n+1}}\left(\frac{1}{1+n}+\frac{2}{n(1+n)}\right)\left(\theta^{\prime\prime}(\xi)+{\rm
O}(\bar{a}^{2}\theta^{(4)}(\xi))\right)\,,$ (26)

with $\theta^{\prime\prime}(\xi)\equiv\partial^{2}\theta(\xi)/\partial\xi^{2}$
and $\theta^{j}(\xi)\equiv\partial^{j}\theta(\xi)/\partial\xi^{j}$. Finally,
the term $I_{n}$ reads:

$\displaystyle I_{n}$ $\displaystyle=$
$\displaystyle\frac{1}{n(n+1)}\int_{\bar{a}}^{1/2}d\xi^{\prime}\frac{\theta^{\prime\prime}(\xi+\xi^{\prime})+\theta^{\prime\prime}(\xi-\xi^{\prime})}{\xi^{\prime
n}}\,,$ (27)

and still contains second derivatives of the phase function. We note that it
can be also rewritten in the form:

$\displaystyle I_{n}$ $\displaystyle=$
$\displaystyle\frac{1}{n(n+1)}\frac{\partial}{\partial\xi}\int_{\bar{a}}^{1/2}d\xi^{\prime}\frac{\theta^{\prime}(\xi+\xi^{\prime})+\theta^{\prime}(\xi-\xi^{\prime})}{\xi^{\prime
n}}\,,$ (28)

so that the equation for the static kink is the integro-differential equation:

$\displaystyle-$ $\displaystyle
d_{n}^{2}\frac{n+2}{n(n+1)}\frac{\partial^{2}\theta}{\partial
x^{2}}+\sin\theta$
$\displaystyle=\frac{1}{n(n+1)}\frac{d_{n}^{2}}{a^{2}}\frac{1}{N^{n+1}}\frac{\partial}{\partial\xi}\int_{\bar{a}}^{1/2}d\xi^{\prime}\frac{\theta^{\prime}(\xi+\xi^{\prime})+\theta^{\prime}(\xi-\xi^{\prime})}{\xi^{\prime
n}}\,.$

Here, the contribution due to $I_{\rm edge}$ has been discarded, and we have
used $N\bar{a}=1$ and $x=Na\xi$. Moreover, we have discarded higher order
derivatives $\theta^{(2j)}(a)$, which account for discreteness effects. If we
would take now the thermodynamic limit, and thus let $N\to\infty$, then this
expression would coincide with the one reported in Ref. Braun:1990 , apart for
the definition of the scaling coefficients and for higher order local
derivatives. We note, however, that the integral term in Eq. (III.2) still
contains terms which can be of the same order as $k^{2}a^{2}$. In order to
single them out, we perform a further step of the partial integration. We
identify two qualitatively different cases, the case $n>1$ and the Coulomb
case $n=1$, which we discuss individually below.

####  III.2.1 Power-law interactions with $n>1$

Performing partial integration of $I_{n}$ for $n>1$ we obtain

$\displaystyle I_{n}$ $\displaystyle=$
$\displaystyle-\frac{1}{n(n^{2}-1)}\frac{\theta^{\prime\prime}(\xi+\xi^{\prime})+\theta^{\prime\prime}(\xi-\xi^{\prime})}{\xi^{\prime
n-1}}\Bigl{|}_{\bar{a}}^{1/2}+I_{n}^{\prime}\,,$ (30)

and

$\displaystyle I_{n}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{1}{n(n^{2}-1)}\int_{\bar{a}}^{1/2}d\xi^{\prime}\frac{\theta^{(3)}(\xi+\xi^{\prime})-\theta^{(3)}(\xi-\xi^{\prime})}{\xi^{\prime
n-1}}\,,$ (31)

We now collect the edge contributions and verify that they scale like
$1/N^{n-1}$, thus we neglect them under the reasonable assumption that the
kink derivatives vanish at the edges. Using that $N\bar{a}=1$ and going back
to dimensional coordinates ($x=Na\xi$), we obtain the integro-differential
equation for a static kink in a chain of atoms interacting via power-law
interactions:

$\displaystyle-$
$\displaystyle\frac{d_{n}^{2}}{n-1}\frac{\partial^{2}}{\partial
x^{2}}\theta(x)+\sin\theta(x)=$
$\displaystyle=\frac{d_{n}^{2}}{n-1}\frac{a^{n-1}}{n(n+1)}\int_{a}^{L/2}dx^{\prime}\frac{\theta^{(3)}(x+x^{\prime})-\theta^{(3)}(x-x^{\prime})}{x^{\prime
n-1}}\,,$

where $L=Na$ is the chain’s length, and now the third-order derivative is
taken with respect to the dimensional coordinates.

####  III.2.2 Coulomb interactions

Fore $n=1$ partial integration of Eq. (27) leads to the expression

$\displaystyle I_{n=1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\theta^{\prime\prime}(\xi+\xi^{\prime})+\theta^{\prime\prime}(\xi-\xi^{\prime}))\log\xi^{\prime}\Bigl{|}_{\bar{a}}^{1/2}+I_{1}^{\prime}\,,$
(33)

where

$\displaystyle I_{1}^{\prime}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\int_{\bar{a}}^{1/2}d\xi^{\prime}(\theta^{(3)}(\xi+\xi^{\prime})-\theta^{(3)}(\xi-\xi^{\prime}))\log\xi^{\prime}\,.$
(34)

Neglecting the contributions from the edges we obtain the equation for a
static kink in a sufficiently long chain of single-component charges:

$\displaystyle-$ $\displaystyle d_{1}^{2}\left(\frac{3}{2}+\log
N\right)\frac{\partial^{2}}{\partial x^{2}}\theta+\sin\theta(x)$
$\displaystyle=-\frac{d_{1}^{2}}{2}\int_{a}^{L/2}dx^{\prime}(\theta^{(3)}(x+x^{\prime})-\theta^{(3)}(x-x^{\prime}))\log(x^{\prime}/L)\,.$

###  III.3 Discussion

(a) $N=101$

(b) $N=301$

(c) $N=1001$

Figure 1: The normalized phase variable $\theta_{j}/2\pi$ of particle $j$ [Eq.
(8)], measuring the particle’s deviation in the kink solution from its ground
state position along the chain. The circles are the result of numerically
finding an equilibrium solution of Eq. (II.1), with
$K_{n}=12,n=1,m=1,V_{0}=1,a=2\pi,b=2\pi$. The dotted red line is obtained by
plotting the analytic sine-Gordon kink with nearest-neighbor (NN) coupling
(with the kink width $d$ given in Eq. (17) taking $K_{\infty}=K_{1}=12$). The
dashed black line is obtained using the same analytic expression with the kink
width $\bar{d}_{1}\propto({\frac{3}{2}+\log N})^{1/2}$ in Eq. (37), setting
$\alpha=0$. The number of simulated particles in the numerical solution is
increased for successive panels with (a) 101, (b) 301, and (c) 1001 particles,
with open boundary conditions. In (b) and (c), only the 100 particle at the
center are plotted. The convergence of the numerical solution toward the
analytic formula is visible.

The integro-differential equations (III.2.1) and (III.2.2) are characterised
by a left-hand side (LHS), which is a SG equation, and a the right-hand side
(RHS), which is an integral term depending on the kink. We first observe that
a priori the integral term on the RHS cannot be discarded. This term, in
particular, is responsible for the behaviour at the edges of the chain, far
away from the kink’s core: Simple considerations show that the tail of the
kink decays algebraically with $1/x^{n+1}$. This behaviour is in contrast to
the exponential decay of the sine-Gordon kink at the asymptotics in Eq. (18).
It is recovered by inspecting the behaviour of the integral term at distances
much larger than the kink core, where one can replace the kink’s first
derivative with a Dirac-delta function Braun:1990 , or can be derived from the
general properties of the interactions Pokrovsky:LR . In this limit, in
particular, the second derivative on the LHS can be neglected. Thus, the
integral term is majorly responsible for the non-local properties of the kink
and gives rise to a power-law interactions between distant kinks Pokrovsky:LR
; Braun:1990 .

At the kink’s core the integral term is negligible for $n>1$ as long as the
ratio $a/d$ is sufficiently small. Analytical estimates and numerical
calculations indicate that the RHS of Eq. (III.2.1) scales approximately like
$a/d$ for $n=2$ and has a sharper decay for larger $n$. Therefore, for $n>1$
the integral term scales like discreteness effects at the kink’s core, and in
the continuum limit the kink’s core is determined to good approximation by a
sine-Gordon equation with a rescaled kink’s width

$\bar{d}_{n>1}=\frac{d_{n}}{\sqrt{n-1}}=\frac{a}{\sqrt{n-1}}\sqrt{\frac{K_{n}}{V_{0}}\left(\frac{b}{2\pi}\right)^{2}}\,.$
(36)

It is also interesting to analyse the scaling of $d_{n}$ with $n$: For $n$
finite it increases as $n$ decreases: In fact, the repulsive interactions
become locally increasingly strong. The nearest-neighbour case can be
recovered by appropriately rescaling the elastic constant $K_{n}$, such as
$K_{n}\sim nK_{\infty}$. In this limit the term on the right-hand side of Eq.
(III.2.1) tends to zero and one recovers the sine-Gordon equation.

The behaviour at the kink’s core for $n=1$, corresponding to repulsive Coulomb
interactions, shall be discussed apart. For this case we have considered a SG
kink and numerically verified that this slowly converges to the solution of
Eq. (III.2.2) as $a/d$ becomes smaller. For finite but large chain the
integral term can be approximated by the function $-\alpha\sin\theta$, plus a
correction which is negligible at the kink’s core. The kink’s equation at the
core is then given by a SG equation, with the kink’s width

$\bar{d}_{n=1}=\frac{a}{\sqrt{1+\alpha}}\sqrt{\frac{3}{2}+\log
N}\sqrt{\frac{K_{1}}{V_{0}}\left(\frac{b}{2\pi}\right)^{2}}\,.$ (37)

where $\alpha=\alpha(a/d)$ and $1>\alpha>0$. This coefficient monotonously
decreases with $a/d$ for the values we checked. These predictions are in
excellent agreement with the numerical results with a discrete chain of ions.
Figure 1 displays the kink’s solution for $N=100,300,500$. The solid blue line
corresponds to a Sine-Gordon kink whose width is given by Eq. (37), the
convergence of the behaviour at the kink’s center with the SG kink is visible.
On the basis of these numerical analysis we conclude that the kink at the core
is described by a SG kink whose width is proportional to $\sqrt{\log N}$ and
thus weakly depends on the chain’s length.

This dependence of the kink’s width on the ions number is a consequence of the
weak non-additivity of Coulomb systems in one dimension. Extensivity can be
formally re-established for instance, by rescaling the spring constant as

$K_{1}=K/\log N\,,$ (38)

which is an implementation of Kac’s scaling for one-dimensional Coulomb
interactions Kac ; Morigi:2004:PRL . For a static kink, this is equivalent to
increasing the depth of the potential as $V_{0}\to V_{0}\log N$. This
rescaling leads to the definition:

$\bar{d}_{1}^{\rm
Kac}=a\sqrt{\frac{K}{V_{0}}\left(\frac{b}{2\pi}\right)^{2}}\,.$ (39)

Thus, in the thermodynamic limit the core of the static kink is described by a
SG kink with width (39), while at the tails the kink decays as $1/x^{2}$.

##  IV Conclusions

We have shown a procedure that permits to determine the equation for a static
kink in a chain of atoms interacting with a repulsive, power-law potential,
and to infer the properties of the solution. Starting from the discrete
equation we have taken a continuum limit and cast the equation into the sum of
a local term, which has the form of a SG equation, and a integral term. We
have argued that in the continuum limit the SG equation determines the
properties at the kink’s core, while the integral gives the behaviour at the
tails. The correction of the integral term to the behaviour at the kink’s
core, in particular, are of the same order of the discreteness effects. These
effects modify the kink’s width and form. Moreover, they significantly modify
the dynamical properties Willis ; Vladan:preprint .

The formalism discussed here can be extended by including a non-homogeneous
density distribution, as it is the case in the presence of a harmonic trap. In
this case the gradient expansion shall be performed including the density in
the continuum limit and will generally give rise to first-order derivatives of
the phase function Morigi:2004 .

The experimental study of our findings could be pursued with several physical
systems. With trapped ions, for instance, recent years saw significant
advances in the trapping of long chains PhysRevA.95.013413 , and subjecting
trapped ions to optical lattices (though still with smaller numbers of ions
PhysRevA.99.031401 ). These findings could also be observed in dipolar gases
Ferlaino  and Rydberg excitons Rydberg:Review  and in the future perhaps also
in magnetically repelling colloids straube2011pattern .

###### Acknowledgements.

The authors are deeply grateful for the privilege they had to know and to work
with Shmuel Fishman. Shmuel Fishman was a great scientist, a person of high
moral standards, and a generous friend. This work is dedicated to his memory.
We thank Eugene Demler, Thomas Fogarty, David Mukamel, and Vladimir
Stojanovic, for stimulating discussions and helpful comments. We acknowledge
the contribution of Andreas A. Buchheit in the initial stage of this project.
This project was supported by the German Research Foundation (the priority
program No. 1929 GiRyd), the European Commission (ITN “ColOpt”), and the
German Ministry of Education and Research (BMBF, QuantERA project “NAQUAS”).
Project NAQUAS has received funding from the QuantERA ERA-NET Cofund in
Quantum Technologies implemented within the European Union’s Horizon 2020
program. G.M. thanks P. Eschner-Morigi for insightful comments. C.C.
acknowledges funding from grant BID-PICT 2015-2236. H.L. thanks Roni Geffen
for fruitful discussions, and acknowledges support by IRS-IQUPS of Université
Paris-Saclay and by LabEx PALM under grant number ANR-10-LABX-0039-PALM.

## References

  * (1) Ya. I. Frenkel and T. A. Kontorova, The model of dislocation in solid body, Zh. Eksp. Teor. Fiz. 8, 1340 (1938) [Phys. Z. Sowjetunion 13, 1 (1938)].
  * (2) P. Bak, Rep. Prog. Phys. 45, 587 (1982).
  * (3) A. Vanossi, N. Manini, M. Urbakh, S. Zapperi, and E. Tosatti, Rev. Mod. Phys. 85, 529 (2013).
  * (4) O. M. Braun and Yu. S. Kivshar, The Frenkel-Kontorova Model: Concepts, Methods, and Applications (Springer, New York, 2004).
  * (5) V. L. Pokrovskij and A. L. Talapov, Theory of incommensurate crystals (Harwood Academic Publishers, New York, 1984).
  * (6) V. L. Pokrovskij and A. L. Talapov, Sov. Phys. JETP 48, 579 (1978).
  * (7) F. C. Frank and J. H. van der Merwe, Proc. R. Soc. Lond. A 198, 205 (1949).
  * (8) C. Rebbi and G. Soliani, Solitons and particles (World Scientific, Singapore, 1984).
  * (9) D. H. E. Dubin and T. M. O’Neil, Rev. Mod. Phys. 71, 87 (1999).
  * (10) S. Baier, M. J. Mark, D. Petter, K. Aikawa, L. Chomaz, Zi Cai, M. Baranov, P. Zoller, and F. Ferlaino, Science 352, 201 (2016).
  * (11) P Schauss, Quantum Sci. Technol 3, 023001 (2018).
  * (12) I. Garcia-Mata, O.V. Zhirov, and D. L. Shepelyansky, Eur. Phys. J. D 41, 325 (2007).
  * (13) T. Pruttivarasin, M. Ramm, I. Talukdar, A. Kreuter, and H. Häffner, New J. Phys. 13, 075012 (2011).
  * (14) A. Benassi, A. Vanossi, and E. Tosatti, Nature Commun. 2, 236 (2011).
  * (15) R. B. Linnet, I. D. Leroux, M. Marciante, A. Dantan, and M. Drewsen, Phys. Rev. Lett. 109, 233005 (2012).
  * (16) M. Enderlein, T. Huber, C. Schneider, and T. Schaetz. Phys. Rev. Lett. 109, 233004 (2012).
  * (17) M. Cetina, A. Bylinskii, L. Karpa, D. Gangloff, K. M. Beck, Y. Ge, M. Scholz, A. T. Grier, I. Chuang, and V. Vuletić, New J. Phys. 15, 053001 (2013).
  * (18) L. Hornekær, N. Kjærgaard, A. M. Thommesen, and M. Drewsen, Phys. Rev. Lett. 86, 1994 (2001).
  * (19) J. Kiethe, R. Nigmatullin, D. Kalincev, T. Schmirander, T. E. Mehlstäubler, Nature Commun. 8, 15364 (2017).
  * (20) H.-K. Li, E.Urban, C. Noel, A. Chuang, Y. Xia, A. Ransford, B. Hemmerling, Y. Wang, T. Li, H. Häffner, X. Zhang, Phys. Rev. Lett. 118, 053001 (2017).
  * (21) K. Pyka, J. Keller, H. L. Partner, R. Nigmatullin, T. Burgermeister, D. M. Meier, K. Kuhlmann, A. Retzker, M. B Plenio, W. H. Zurek, A. Del Campo, and T. E. Mehlstäubler, Nature Commun. 4, 2291 (2013).
  * (22) S. Ulm, J. Rossnagel, G. Jacob, C. Degunther, S. T. Dawkins, U. G. Poschinger, R. Nigmatullin, A. Retzker, M. B. Plenio, F. Schmidt-Kaler, and K. Singer, Nature Commun. 4, 2290 (2013).
  * (23) M. Mielenz, J. Brox, S. Kahra, G. Leschhorn, M. Albert, T. Schaetz, H. Landa, and B. Reznik Phys. Rev. Lett. 110, 133004 (2013).
  * (24) C. Senko, J. Smith, P. Richerme, A. Lee, W.C. Campbell, and C. Monroe, Science 345, 430 (2014).
  * (25) W. S. Bakr, J. Gillen, A. Peng, S. Fölling, and M. Greiner, Nature 462, 74 (2009).
  * (26) M. Viteau, M. G. Bason, J. Radogostowicz, N. Malossi, D. Ciampini, O. Morsch, and E. Arimondo, Phys. Rev. Lett. 107, 060402 (2011).
  * (27) J. Brox, P. Kiefer, M. Bujak, H. Landa, T. Schaetz, Phys. Rev. Lett. 119, 153602 (2017)
  * (28) Ch. Schneider, D. Porras, and T. Schaetz, Rep. Prog. Phys. 75, 024401 (2012).
  * (29) A. Bylinskii, D. Gangloff, and V. Vuletić. Science 348 6239, 1115-1118 (2015).
  * (30) D. Gangloff, A. Bylinskii, I. Counts, W. Jhe, and V. Vuletić. Nature Phys. 11 11, 915-919 (2015).
  * (31) I. Counts, D.Gangloff, A. Bylinskii, J. Hur, R. Islam, and V. Vuletić. Phys. Rev. Lett., 119 043601 (2017).
  * (32) A. Campa, T. Dauxois, and S. Ruffo, Phys. Rep. 480, 57 (2009).
  * (33) N. Laskin and G. Zaslavsky, Physica A: Stat. Mech. and App. 368, 38 (2006).
  * (34) K. Kirkpatrick, E. Lenzmann, and G. Staffilani. Communications in Mathematical Physics 317, 563 (2013).
  * (35) P. Bak and R. Bruinsma, Phys. Rev. Lett. 49, 249 (1982).
  * (36) S. F. Mingaleev, Y. B. Gaididei, E. Majernikova, S. Shpyrko, Phys. Rev. E 61, 4454 (2000).
  * (37) V L Pokrovsky and A Virosztek, Journal of Physics C: Solid State Physics 16, 4513 (1983).
  * (38) O. M. Braun, Yu. S. Kivshar, and I. I. Zelenskaya, Phys. Rev. B 41, 7118 (1990).
  * (39) G. Morigi and S. Fishman, Phys. Rev. E 70, 066141 (2004).
  * (40) G. Morigi and S. Fishman, Phys. Rev. Lett. 93, 170602 (2004).
  * (41) D. H. E. Dubin, Phys. Rev. E 55, 4017 (1997).
  * (42) C. Willis, M. El-Batanouny, and P. Stancioff, Phys. Rev. B 33, 1904 (1986).
  * (43) D. A. Gangloff, A. Bylinskii, and V. Vuletić, Physical Review Research 2, 013380 (2020).
  * (44) M. R. Kamsap, C. Champenois, J. Pedregosa-Gutierrez, S. Mahler, M. Houssin, and M. Knoop, Phys. Rev. A 95, 013413 (2017).
  * (45) T. Lauprêtre, R. B. Linnet, I. D. Leroux, H. Landa, A. Dantan, and M. Drewsen, Phys. Rev. A 99, 031401 (2019).
  * (46) A. V. Straube, A. A. Louis, J. Baumgartl, C. Bechinger, and R. P. A. Dullens, EPL (Europhysics Letters) 94, 48008 (2011).


#  Bayesian graph neural networks for strain–based crack localization

C. Mylonas1 G. Tsialiamanis2 K. Worden2 E. N. Chatzi1
1Institute of Structural Engineering, ETH Zürich,
Stefano-Franscini-Platz 5, CH-8093 Zürich, Switzerland
2 Dynamics Research Group, Department of Mechanical Engineering,
University of Sheffield Mappin Street, Sheffield S1 3JD, UK

## Abstract

A common shortcoming of vibration-based damage localization techniques is that
localized damages, i.e. small cracks, have a limited influence on the spectral
characteristics of a structure. In contrast, even the smallest of defects,
under particular loading conditions, cause localized strain concentrations
with predictable spatial configuration. However, the effect of a small defect
on strain decays quickly with distance from the defect, making strain-based
localization rather challenging. In this work, an attempt is made to
approximate, in a fully data-driven manner, the posterior distribution of a
crack location, given arbitrary dynamic strain measurements at arbitrary
discrete locations on a structure. The proposed technique leverages Graph
Neural Networks (GNNs) and recent developments in scalable learning for
Bayesian neural networks. The technique is demonstrated on the problem of
inferring the position of an unknown crack via patterns of dynamic strain
field measurements at discrete locations. The dataset consists of simulations
of a hollow tube under random time-dependent excitations with randomly sampled
crack geometry and orientation.

Keywords: GraphNets, Strain-based crack localization, Bayesian Graph Neural
Network

This is a preprint of a work submitted for possible publication on conference
proceedings. Copyright may be transferred without notice, after which this
version may no longer be accessible.

##  1 Introduction

When the region around a defect is loaded, strain concentration occurs. For
different local loading conditions, depending on the geometric configuration
of the defect, different strain concentration patterns arise. The technique
proposed in this work, allows for strain-based localization of cracks with
learned, rather than hand-engineered, features. In contrast to the approaches
based on 2D Convolutional Neural Networks (CNNs) already in the literature
[1], the approach presented herein can effectively deal with arbitrary
(potentially sparse) strain sensor positions, and not only sensing positions
arranged in a grid as CNN-based computation dictates. This feature is achieved
by employing a message-passing Graph Neural Network (GNN) [2, 3]. The GNN is
used for recognising spatial strain patterns, and learns these patterns
directly from data. These patterns depend non-linearly on the sensor position
relative to the crack and the sensor position in the structure.

In essence, as CNNs perform local aggregation defined by a spatial 2D grid,
message-passing graph neural networks allow for local aggregation on the
neighborhood of the nodes of an arbitrary graph.

Figure 1: Schematic of the construction of a strain sensing graph. Each node
corresponds to a full strain-sensor. The neighborhood of each sensor is
defined by the sensors closer to it. An ellipsoidal crack is depicted in red.

In contrast to other works in the literature [4, 5], a model that predicts
directly the crack position from data is learned and no optimization is
performed during prediction time.

#### Training Data And Strain-PCA contrasting:

In this study, the dynamics of 450 hollow tubes with randomly-positioned and
rotated ellipsoidal defects were simulated with the Finite Element Method
under random white-noise excitations in three dimensions, different for every
run. The full strain tensors on a $150\times 150$ grid of values on the
surface of the tube were extracted. All simulated tubes are 10m in length and
1m in diameter.

A basis is first obtained from a Principal Components Analysis (PCA) of strain
snapshots of a crack-free structure. In the presented results, 50 components
of the PCA basis are retained (variance explained $88.3\%$). A separate PCA
for each component of the strain tensor and the strain-invariants is learned.
The strain fields corresponding to structures that contain a defect are pre-
processed by projecting their strain fields to the corresponding PCA-derived
bases. The projected strains are removed from the observed strains. Assuming a
linear superposition of the strain related to the defect and the strain of the
defect-free structure, the strain associated with the defect should appear in
the residuals of the projection. Indeed, this effect is demonstrated in the
right-most plot of Figure 2. During prediction, only a sparse set of sensors
is expected to be available. Therefore, the estimates of the strain field are
going to be less accurate compared to the case where the full strain field is
available. The projection using a sparse mask of strain sensors, is equivalent
to performing a simple linear regression with the PCA bases. The PCA
approximation of $\hat{\varepsilon}$, for a strain component $\varepsilon$
based on a sparse set of sensors at $i$ positions $S:\\{s_{1},\dots s_{i}\\}$,
is computed from

$\displaystyle\bar{\Psi}_{\varepsilon}$
$\displaystyle\stackrel{{\scriptstyle\mathclap{\mbox{S}}}}{{\leftarrow}}\Psi_{\varepsilon}$
$\displaystyle\hat{c}$
$\displaystyle=(\bar{\Psi}_{\varepsilon}^{T}\bar{\Psi}_{\varepsilon})^{-1}\bar{\Psi}_{\varepsilon}^{T}\bar{\varepsilon}$
$\displaystyle\hat{\varepsilon}$
$\displaystyle=\bar{\Psi}_{\varepsilon}\hat{c}$ $\displaystyle
r^{PCA}_{\varepsilon}$ $\displaystyle=\hat{\varepsilon}-\varepsilon$

where $\Psi_{\varepsilon}$ is the matrix of PCA components,
$\bar{\Psi}_{\varepsilon}$ is the matrix of PCA components, sparsified at the
sensor positions, $\hat{c}$ is a vector of coefficients and
$r^{PCA}_{\varepsilon}$ is the PCA-computed residual for that strain
component. In what follows, for notational convenience the PCA-derived strain
residuals $r^{PCA}_{\varepsilon_{ij}}$ are referred to simply as strains
$\varepsilon_{ij}$.

A representative example of the quality of the approximation for the
investigated structure, and for the crack geometries used is depicted in
Figure 2.

Figure 2: Contrasting sparse measurements using strain field PCA. The
positions on the plots correspond to length-wise and angular cylinder
coordinates. The left-most plot is an approximation of the defect-dependent
residual based on the projection on a limited number of sensors with
$\bar{\Psi}_{\varepsilon}$. The middle plot is the actual residual based on a
projection using the dense measurements and full PCA components
$\Psi_{\varepsilon}$. The rightmost plot is the full residual strain field for
$\varepsilon_{22}$ for a time-instant where $\varepsilon_{22}$ reveals the
position of the strain.

##  2 Message-Passing Graph Neural Networks

A directed attributed graph $G:(V,E)$, (referred to in the following simply as
“graph”) is defined as a set of vertices (nodes)
$V:\\{\mathbf{v}_{1},\dots\mathbf{v}_{N^{n}}\\}$ and a set of edges
$E:\\{(\mathbf{e_{1}},s_{1},r_{1}),\cdots(\mathbf{e}_{N^{e}},s_{N^{e}},r_{N^{e}}\\}$.
Nodes contain attributes, and edges are described by an attribute
$\mathbf{e}_{i}$, and sender $s_{i}$ and receiver $r_{i}$ indices. Indices
$s_{i},r_{i}$ are integers representing the sender and receiver nodes.
Attributes may be discrete or continuous vectors or tensors. In [3] the
GraphNet (GN) computational block is proposed, which is a generalisation of
several approaches to computation using attributed graphs, with a focus on
implementing inductive biases. GNs are algorithms that operate on graphs and
return graphs. The attributes of the vertices and edges, parametrize how the
entities in the graph interact. In addition, graph-level (global) attributes
may exist, which depend on the aggregation of all vertex and node information.
More concretely, for the problem considered, the strain sensors have the role
of vertices, their relative positions are edge features (which dictate which
nodes will interact to yield a prediction), and the position and orientation
of the crack is the global variable which needs to be inferred from the sensor
observations.

A GraphNet computational block consists of

  * •

an edge function $f^{e}:\mathbf{e}\rightarrow\mathbf{e}^{\prime}$ which
parametrizes what information is passed along an edge ;

  * •

an edge message aggregation function $\rho^{e\rightarrow v}$, which is a
permutation invariant function that summarizes the messages to be passed
according to the graph connectivity;

  * •

a vertex (node) function that transforms the vertex attributes
$f^{v}:\mathbf{v}\rightarrow\mathbf{v}^{\prime}$, optionally using the
aggregated incoming messages, the node attributes and the global attributes,
and finally;

  * •

a global function $f^{u}:\mathbf{u}\rightarrow\mathbf{u}^{\prime}$ ; which
uses global attributes. Node and edge aggregated information which are
aggregated by two additional aggregators:

  * •

the edge-to-global aggregator $\rho^{e\rightarrow u}$ and

  * •

the node-to-global aggregator $\rho^{v\rightarrow u}$.

The interaction of vertices in GraphNets occurs through a process that
resembles message passing [2], with a learned message passing function. In
addition to the message-passing function, transformations of the edge and node
attributes are performed, for instance by a Multi-Layer Perceptron (MLP).
Graph Neural Networks (GNNs) are GNs that use MLPs for one or several of the
aforementioned functions. The reader is referred to [3] for a more detailed
and complete exposition on GraphNets. A custom TensorFlow-based GraphNets
library was used 111 Code can be accessed at: https://github.com/mylonasc/tf-
gnns. The actual computation graph used in this work is given in Figure 4. The
construction of the graphs and the choices for the node, edge, global and
aggregation functions is, in general, problem specific. In the following, the
feature construction and some architectural choices are detailed.

#### Node feature construction:

Since the strain-input of the network is a time-series, it is natural to
process it with a 1D convolutional neural network. The raw strain time-series
is augmented by the two first rotational strain invariants, which read,

$\displaystyle I_{1}(t,\mathbf{p_{i}})$
$\displaystyle=\varepsilon_{11}(t,\mathbf{p_{i}})+\varepsilon_{22}(t,\mathbf{p_{i}})+\varepsilon_{33}(t,\mathbf{p_{i}})$
(1) $\displaystyle I_{2}(t,\mathbf{p_{i}})$
$\displaystyle=\varepsilon_{12}(t,\mathbf{p_{i}})^{2}+\varepsilon_{23}(t,\mathbf{p_{i}})^{2}+\varepsilon_{13}(t,\mathbf{p_{i}})^{2}-$
(2)
$\displaystyle\quad\varepsilon_{11}(t,\mathbf{p_{i}})\varepsilon_{22}(t,\mathbf{p_{i}})-\varepsilon_{22}(t,\mathbf{p_{i}})\varepsilon_{33}(t,\mathbf{p_{i}})-\varepsilon_{33}(t,\mathbf{p_{i}})\varepsilon_{11}(t,\mathbf{p_{i}})$
(3)

where $\mathbf{p_{i}}=[p_{L},p_{\phi}]$ is the position of sensor $i$. The
position of the sensors along the length of the tube and the angular
dimension, is appended as an additional filter dimension. The angular position
$p_{\phi}$ is encoded as $sin(p_{\phi}),cos(p_{\phi})$.

#### Edge features:

In order to localize the crack, information between several nearby sensors may
need to be taken into account. Groups of sensors should collectively detect
patterns, parametrized by their relative position. This is implemented by
endowing the edges of the created graph with features that encode this
information. Denoting $s_{i}$ and $r_{i}$ the sender and receiver node, the
features given in Figure 3 are used.

Figure 3: Edge features parametrize the manner in which the different nodes in
the GNN exchange information. To that end, non linear relative geometric
information is included in the edge function.

The 1D convolutional network captures temporal variations whereas the GNN
should capture parameterized spatial patterns that are associated with the
existence of a crack.

#### Construction of the strain-sensor graph:

Since local features of the strain field are more informative for the problem,
the spatial proximities of the randomly-placed strain sensors are used to
define the graphs employed in the computations, as depicted in Figure 1. For
each sensor graph, a set of 150 strain sensors are used ($\frac{150}{10\cdot
2\cdot\pi\cdot 0.5m^{2}}\approx 4.8[sensors]/[m^{2}]$). Sensors at a distance
less than 10cm from the crack are removed to avoid detecting the crack only
from the local strain concentration magnitude, which renders the problem much
less challenging. Different and unique strain sensor graphs and strain sensor
positions were used for each training and testing input. In essense, the crack
positions, the loading conditions and the strain sensor configurations
encountered at testing, are never encountered during training. The edges in
the network are bi-directional, i.e. for each directional edge from node $i$
to node $j$ there is a directional edge from node $j$ to node $i$.

#### Local reparametrization and uncertainty quantification

Since, in general, the defect localization problem may be ill-posed222Consider
the cases where the structure is unloaded, or when the structure is loaded in
a manner that does not deform the defect. In these two extreme cases, the
problem is not solvable. Moreover, there are multiple material configurations
that yield the same discrete strain observations for the same loading., the
uncertainty in the predictions and learned model needs to be represented.
Recent scalable approaches to representing model uncertainty in deep learning
include Monte-Carlo Dropout [6] and variational inference with stochastic
gradient descent [7, 8]. Variational inference relies on an approximation of
the maximum likelihood objective, namely the Evidence Lower Bound (ELBO). In
variational techniques, a distribution over the parameters of network layers
is inferred from data. The local reparametrization from [7] is used, since 1D
regression experiments showed faster and more consistent training. Recent
literature for Bayesian GNNs [9] supports these observations.

###  2.1 Model architecture

The node inputs are passed through a 1D-CNN that acts as the node function of
the input GNN $GN^{(input)}$. A simple 2-layer CNN with global max-pooling was
chosen, as it is expected that the dynamics play a minor role in identifying
the position of the crack from strain. Simple two-layer ReLU Multi-Layer
Perceptrons (MLPs) were used as an edge function for the input GNN. No message
passing is performed in the input GNN (Graph-independent GNN block). Both
$GN^{(input)}$ edge function outputs and $GN^{(input)}$ node-function outputs
have a size of 32 units. Networks with fewer units also demonstrated
relatively good performance. Following the input GNN, two full GraphNet (core)
layers are applied. Residual connections were used for all MLPs involved in
$GN^{(core)}$. In the core MLP, instead of the simple two-layer ReLU MLPs, for
both the edge and node MLPs an MLP that contains a Bayesian Neural Network
layer (BNN layer) is used. The tensorflow-distributions implementation [10] of
the local-reparametrization layer [7] was used. The inputs are first passed to
a fully-connected layer, a single Bayesian Neural Network (BNN) layer with
ReLU activation is then applied, and the MLPs conclude with dense layers with
no activations. A mean-aggregator was used for the edge-to-node aggregators
($\rho^{e\rightarrow v}$).

In the final computation step, a graph-to-global layer is employed
($GN^{(output)}$), which directly returns estimates on the position of the
crack. The nodes and edges output from the $GN^{(core)}$ are mean-aggregated
and concatenated. The same type of BNN-MLP used in the core network is used,
with different output of three units, corresponding to
$p_{\phi},p_{L},p_{\psi}$, which are the crack parameters.

Some preliminary experiments were performed on the performance of the proposed
architecture for different choices of hyper-parameters. Experiments were
performed with more $GN^{(core)}$ steps (message-passing through repeated
application of the $GN^{(core)}$ network) with and without BNN layers, and it
was found that two and three layers offered comparable results. More than
three $GN^{(core)}$ steps did not offer any improvements. Networks with a
single message-passing step, regardless of the size of the involved MLPs,
performed badly. This observation provides some preliminary evidence that
aggregating the information from several sensors is important for the crack
localization problem.

#### Training loss and train-test split

Denoting as $G^{in}$ the input graphs, as
$\mathbf{p}:[p_{\phi},p_{L},p_{\psi}]$ the corresponding crack positions,
$\mathcal{D}:\\{G^{in},\mathbf{p}\\}$ the dataset, as $p(w)$ a prior
distribution over the parameters of the network, the loss function (ELBO)
reads

$\displaystyle\mathcal{L(\phi,\mathbf{w},\mathcal{D})}=-D_{KL}(q_{\phi}(\mathbf{w})||p(\mathbf{w}))+L_{\mathcal{D}}(\phi)$

where $\mathbf{w}$ are all the parameters of the associated MLPs, $\phi$ are
the parameters of the variational approximation to the posterior
$q_{\phi}(\mathbf{w})$, $p(\mathbf{w})$ is a prior distribution (a factored
Gaussian). Moreover, $D_{KL}(q_{\phi}(\mathbf{w})||p(\mathbf{w}))$ is the
Kullback-Leibler divergence, and $L_{\mathcal{D}}$ is

$L_{\mathcal{D}}(\phi)=\sum_{(G^{in},\mathbf{p})\in\mathcal{D}}E_{q_{\phi}(\mathbf{w})}[\log
p(G^{in}|\mathbf{p},\mathbf{w})]$

which is the expected negative log-likelihood. In addition, we are assuming
$p(\mathbf{p})\sim\mathcal{N}(\mu_{\mathbf{p}},\sigma^{2}_{\mathbf{p}})$, and
$q_{\phi}(\mathbf{w})\sim\mathcal{N}(\mu_{\mathbf{w}},\sigma^{2}_{\mathbf{w}})$.
Low variance stochastic gradient estimators exist for the loss function [7,
11].

350 experiments were used as training set and the remaining 100 as testing
set. The loss on the testing set was monitored and early stopping was applied
after 50 epochs. Random segments of 150 timesteps from the 401 that are
available (0.5 seconds) were used from each experiment of the training set,
different at each epoch. Interestingly, the network did not suffer from over-
fitting, as evaluated in the test-set of unseen sensor configurations and
crack positions, even though the number of experiments in the raw dataset
number is relatively small (350). This effect is probably a result of the
combination of the small networks used and the fact that the learning signal
from the dataset is richer per datapoint. Specifically, each graph contains
around 150 nodes and several thousands of edges. Essentially, each graph
contributes to learning the edge and node functions in proportion to the
number of edges and nodes in the graph, not as a single datapoint333This
observation motivates augmenting the dataset with more random graphs even for
the same experiment, not a single one as it is the case in this work. This was
not attempted due to resource constraints but this investigation may form part
of future works..

Figure 4: The computational graph used. The $GN^{(core)}$ block is repeatedly
applied in order to allow for message passing between a larger neighborhood
for each node.

##  3 Results

The strain localization results for some random cases from the test-set are
given in Figure 5. The GNN predicts a set of points close to the actual node
position. The results on the angle of the crack were less accurate than the
position estimates. Because of the incorporation of Bayesian neural network
layers, every forward estimation of the GNN the outputs are different, and
each evaluation can be thought of as corresponding to a different model
sampled from the posterior distribution of GNN models that was learned from
the data. The Normalized Root Mean-Squared Error $L^{NRMSE}$, computed as,

$L^{NRMSE}(y_{pred},y_{actual})=\frac{\sqrt{\frac{\sum^{N}(y_{pred}-y_{actual})^{2}}{N}}}{\max(y_{act})-\min(y_{act})}$
(4)

is given together with all the training and test set examples in Figure 7. For
comparison, localization results on the test-set for a deterministic GNN
trained on least-squares loss, with the BNN layers replaced with simple FFNN
layers, is given in Figure 6. Although the results are relatively close to the
correct predictions, no estimates on the uncertainty of the predictions can be
obtained.

Figure 5: Strain-based localization results on test-set for some random cases
for the Bayesian GNN model. The localization is performed on time-windows of
200 observations (250ms). Predictions are denoted in blue. The scatter on the
predictions, is partially because of the Bayesian neural network layers and
partially because of the different inputs for each segment. It is expected
that for some time instants the localization is inaccurate since the crack
does not contribute to the strain field surrounding it and cannot be detected.
Figure 6: Strain-based localization results on test-set for some random cases
for the least-squares GNN model. The localization is performed on two non-
overlapping time-windows of 200 observations (250ms). Predictions are denoted
in blue. The network returns deterministic and incorrect predictions with
confidence. Figure 7: Strain-based localization results - all training and
validation sets.

##  4 Conclusions

In this work, graph neural networks with temporal convolution feature
extractors were employed for a defect localization problem, using arbitrary
configurations of sensors. Because of the end-to-end differentiability of
GNNs, future extensions of this work include optimal sensor placement
optimization problems [12]. The incorporation of data from multiple sensor
types, such as active piezo-electric sensors can also be naturally implemented
in the GNN framework. In addition, GNN-based computation may be a viable
alternative to Bayesian PDE inversion problems [13], where observations are
sparse, in place of the grid-based CNN computational substrate. Finally, GNNs
offer great potential for inference of data-driven Reduced Order Models
(ROMs), directly from physical observations. GNNs have already been
successfully applied to learning complex physics with particle methods [14]
and multi-body dynamics [15] on simulation data.

#### Ackgnowledgement

C. M. and E. C. would like to gratefully acknowledge the support of the
European Research Council via the ERC Starting Grant WINDMIL (ERC-2015- StG
#679843) on the topic of Smart Monitoring, Inspection and Life-Cycle
Assessment of Wind Turbines. G. T. and K. W. were supported by the EC Marie
Skłodowska-Curie ETN grant “DyVirt” (764547).

## References

  * [1] Nur Sila Gulgec, Martin Takáč, and Shamim N Pakzad.  Convolutional neural network approach for robust structural damage detection and localization.  Journal of Computing in Civil Engineering, 33(3):04019005, 2019\.
  * [2] Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl.  Neural message passing for quantum chemistry.  arXiv preprint arXiv:1704.01212, 2017.
  * [3] Peter W Battaglia, Jessica B Hamrick, Victor Bapst, Alvaro Sanchez-Gonzalez, Vinicius Zambaldi, Mateusz Malinowski, Andrea Tacchetti, David Raposo, Adam Santoro, Ryan Faulkner, et al.  Relational inductive biases, deep learning, and graph networks.  arXiv preprint arXiv:1806.01261, 2018.
  * [4] Konstantinos Agathos, Eleni Chatzi, and Stéphane PA Bordas.  Multiple crack detection in 3d using a stable xfem and global optimization.  Computational mechanics, 62(4):835–852, 2018.
  * [5] Konstantinos Agathos, Stéphane PA Bordas, and Eleni Chatzi.  Parametrized reduced order modeling for cracked solids.  International Journal for Numerical Methods in Engineering, 121(20):4537–4565, 2020.
  * [6] Yarin Gal and Zoubin Ghahramani.  Dropout as a bayesian approximation: Representing model uncertainty in deep learning.  In International Conference on Machine Learning, pages 1050–1059, 2016.
  * [7] Durk P Kingma, Tim Salimans, and Max Welling.  Variational dropout and the local reparameterization trick.  In Advances in neural information processing systems, pages 2575–2583, 2015.
  * [8] Charles Blundell, Julien Cornebise, Koray Kavukcuoglu, and Daan Wierstra.  Weight uncertainty in neural networks.  arXiv preprint arXiv:1505.05424, 2015.
  * [9] George Lamb and Brooks Paige.  Bayesian graph neural networks for molecular property prediction.  arXiv preprint arXiv:2012.02089, 2020.
  * [10] Joshua V Dillon, Ian Langmore, Dustin Tran, Eugene Brevdo, Srinivas Vasudevan, Dave Moore, Brian Patton, Alex Alemi, Matt Hoffman, and Rif A Saurous.  Tensorflow distributions.  arXiv preprint arXiv:1711.10604, 2017.
  * [11] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  arXiv preprint arXiv:1312.6114, 2013.
  * [12] Costas Argyris, Sharmistha Chowdhury, Volkmar Zabel, and Costas Papadimitriou.  Bayesian optimal sensor placement for crack identification in structures using strain measurements.  Structural Control and Health Monitoring, 25(5):e2137, 2018.
  * [13] Liu Yang, Xuhui Meng, and George Em Karniadakis.  B-pinns: Bayesian physics-informed neural networks for forward and inverse pde problems with noisy data.  arXiv preprint arXiv:2003.06097, 2020.
  * [14] Alvaro Sanchez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter W Battaglia.  Learning to simulate complex physics with graph networks.  arXiv preprint arXiv:2002.09405, 2020.
  * [15] Alvaro Sanchez-Gonzalez, Nicolas Heess, Jost Tobias Springenberg, Josh Merel, Martin Riedmiller, Raia Hadsell, and Peter Battaglia.  Graph networks as learnable physics engines for inference and control.  arXiv preprint arXiv:1806.01242, 2018.