headers [0][Authors][]Computational inverse scattering with internal sources:
a RKHS approach0  propositionlemma proposition corollarylemma corollary
theoremlemma theorem definitionlemma definition assumptionlemma assumption

# Computational inverse scattering with internal sources: a reproducing kernel
Hilbert space approach

Yakun Dong1
<EMAIL_ADDRESS>Kamran Sadiq2
<EMAIL_ADDRESS>Otmar Scherzer1,2,3
<EMAIL_ADDRESS>John C. Schotland4
<EMAIL_ADDRESS>
1Faculty of Mathematics
University of Vienna
Oskar-Morgenstern-Platz 1
A-1090 Vienna, Austria

2Johann Radon Institute for Computational
and Applied Mathematics (RICAM)
Altenbergerstraße 69
A-4040 Linz, Austria

3Christian Doppler Laboratory for Mathematical
Modeling and Simulation of Next Generations
of Ultrasound Devices (MaMSi)
Oskar-Morgenstern-Platz 1
A-1090 Vienna, Austria

4Department of Mathematics and
Department of Physics
Yale University
New Haven, CT 06520-8283, USA

## Abstract

We present a method to reconstruct the dielectric susceptibility (scattering
potential) of an inhomogeneous scattering medium, based on the solution to the
inverse scattering problem with internal sources. We employ the theory of
reproducing kernel Hilbert spaces, together with regularization to recover the
susceptibility of two- and three-dimensional scattering media. Numerical
examples illustrate the effectiveness of the proposed reconstruction method.

## 1\. Introduction

The inverse scattering problem is concerned with determining the structure of
a scattering medium from measurements of the scattered field. Such problems
have numerous applications including medical imaging, nondestructive testing,
remote and radar sensing, ocean acoustics and geophysical exploration; see
e.g. [13, 7, 22, 4], and the references therein.

In this paper we investigate an inverse scattering problem with internal
sources, and show that it is possible to recover the dielectric susceptibility
(scattering potenial) of an inhomogeneous scattering medium. Suitable internal
sources take the form of photoactivated fluorophores, which are used in
superresolution microscopy [2, 10], and subwavelength bubbles, which are used
in ultrasound localization microscopy [8]. The considered inverse problem is
ill-posed. We apply the theory of reproducing kernel Hilbert spaces (RKHS) [1,
12] accompanied by tools from regularization theory (see e.g. [21, 6]) to
reconstruct the susceptibility.

In [9], a similar inverse problem is considered, also using an integral
equation approach. The integral equation is converted into a linear system of
algebraic equations, which is solved by means of a regularized pseudoinverse
to recover the susceptibility. Instead, here we use the RKHS approach, which
yields reconstructed images of higher quality compared to [9].

This paper is organized as follows. In Section 2, we formulate the forward
problem and the associated inverse scattering problem. We also briefly recall
some tools from the theory of regularization and RKHS, which we need for our
reconstruction method. In addition, we present the method for both two- and
three-dimensional inverse problems. In Section 3, we present the results
obtained by applying the RKHS reconstruction method to three numerical
experiments. We also present some comparison results to [9], as shown in Table
LABEL:tb:comp-john. Finally, in Section 4, we present our concluding remarks.

## 2\. Inverse scattering problem with internal sources

We consider an experiment in which light from a monochromatic source
propagates in an inhomogeneous medium in two dimensions. The three-dimensional
case is taken up in Section 2.3. The source is taken to be located in the
interior of the medium and to consist of a single photoactivated fluorescent
molecule. The scattered light is then registered by a detector in the far
field for each source. We ignore the effects of polarization and consider a
scalar model of the optical field, which satisfies the time-harmonic wave
equation,

$\displaystyle\mathop{}\\!\mathbin{\bigtriangleup}U(\operatorname{\mathbf{r}})+k^{2}(1+\eta(\operatorname{\mathbf{r}}))U(\operatorname{\mathbf{r}})=-a(\operatorname{\mathbf{r}}_{1})\delta(\operatorname{\mathbf{r}}-\operatorname{\mathbf{r}}_{1}),$
(2.1)

where $U$ is the total field, $\eta(\operatorname{\mathbf{r}})$ is the
susceptibility of the medium at the position $\operatorname{\mathbf{r}}$, and
$k$ is the wave number. Since the source consists of a single molecule, it is
taken to be a point source with position $\operatorname{\mathbf{r}}_{1}$ and
amplitude $a(\operatorname{\mathbf{r}}_{1})$. The total field may be
decomposed into the incident field $U_{i}$ and scattered field $U_{s}$.
Following standard procedure [4, 3], the total field satisfies the integral
equation

$U(\operatorname{\mathbf{r}})=U_{i}(\operatorname{\mathbf{r}})+k^{2}\int
G(\operatorname{\mathbf{r}},\operatorname{\mathbf{r}}^{\prime})\eta(\operatorname{\mathbf{r}}^{\prime})U(\operatorname{\mathbf{r}}^{\prime})d\operatorname{\mathbf{r}}^{\prime},$
(2.2)

where

${G}(\operatorname{\mathbf{r}},\operatorname{\mathbf{r}}^{\prime}):=\frac{\mathrm{i}}{4}H_{0}^{(1)}(k|\operatorname{\mathbf{r}}-\operatorname{\mathbf{r}}^{\prime}|)$
(2.3)

is the two-dimensional Green’s function, with $H_{0}^{(1)}$ being the Hankel
function of first kind [5]. The Green’s function satisfies the equation

$\displaystyle\mathop{}\\!\mathbin{\bigtriangleup}_{\mathbf{r}}G(\mathbf{r},\mathbf{r^{\prime}})+k^{2}G(\mathbf{r},\mathbf{r^{\prime}})=-\delta(\mathbf{r}-\mathbf{r^{\prime}}),$
(2.4)

together with the Sommerfeld radiation condition [5].

In the far-zone, the Hankel function has the asymptotic form, see e.g. [5,
Section 3.4 (3.59)]:

$\displaystyle\frac{\mathrm{i}}{4}H_{0}^{(1)}(kR)$
$\displaystyle\sim\frac{e^{\mathrm{i}\frac{\pi}{2}}}{4}\sqrt{\frac{2}{\pi
kR}}e^{-\frac{\pi}{4}\mathrm{i}}e^{\mathrm{i}kR}=\frac{e^{\mathrm{i}\frac{\pi}{4}}}{\sqrt{8\pi
k}}\frac{e^{\mathrm{i}kR}}{\sqrt{R}},$ (2.5)

which propagates as a cylindrical wave with amplitude $R^{-1/2}$. Thus, in the
far field, we obtain the asymptotic formula for the two-dimensional Green’s
function,

$\displaystyle G\left(\mathbf{r},\mathbf{r}^{\prime}\right)$
$\displaystyle=\frac{\mathrm{i}}{4}H_{0}^{(1)}\left(k\left|\mathbf{r}-\mathbf{r}^{\prime}\right|\right)\sim\frac{e^{\mathrm{i}\frac{\pi}{4}}}{\sqrt{8\pi
k}}\frac{e^{\mathrm{i}k\left|\mathbf{r}-\mathbf{r}^{\prime}\right|}}{\sqrt{\left|\mathbf{r}-\mathbf{r}^{\prime}\right|}}$
$\displaystyle\simeq\frac{e^{\mathrm{i}\frac{\pi}{4}}}{\sqrt{8\pi
k}}\frac{e^{\mathrm{i}k|\mathbf{r}|}}{\sqrt{|\mathbf{r}|}}e^{-\mathrm{i}k\hat{\mathbf{r}}\cdot\mathbf{r}^{\prime}}\left(1+\frac{1}{2}\frac{\widehat{\mathbf{r}}\cdot\mathbf{r}^{\prime}}{|\mathbf{r}|}\right).$
(2.6)

Since the incident field is given by
$U_{i}(\mathbf{r})=a\left(\mathbf{r}_{1}\right)G\left(\mathbf{r},\mathbf{r}_{1}\right)$,
we derive by using (2.6) the asymptotic behavior of the total field,

$\displaystyle U(\mathbf{r})$
$\displaystyle=a\left(\mathbf{r}_{1}\right)G\left(\mathbf{r},\mathbf{r}_{1}\right)+k^{2}\int
G\left(\mathbf{r},\mathbf{r}^{\prime}\right)\eta\left(\mathbf{r}^{\prime}\right)U\left(\mathbf{r}^{\prime}\right)d\mathbf{r}^{\prime}$
(2.7)
$\displaystyle\sim\frac{e^{\mathrm{i}\frac{\pi}{4}}}{\sqrt{8\pi}}\frac{e^{\mathrm{i}k|\mathbf{r}|}}{\sqrt{|\mathbf{r}|}}\left(e^{-\mathrm{i}k\hat{\mathbf{r}}\cdot\mathbf{r}_{1}}\frac{a\left(\mathbf{r}_{1}\right)}{\sqrt{k}}+k^{\frac{3}{2}}\int
e^{-\mathrm{i}k\hat{\mathbf{r}}\cdot\mathbf{r}^{\prime}}\eta\left(\mathbf{r}^{\prime}\right)U\left(\mathbf{r}^{\prime}\right)d\mathbf{r}^{\prime}\right)$

Thus, in the far zone of the scatterer, the field behaves as an outgoing
cylindrical wave:

$U(\mathbf{r})\sim\frac{e^{\mathrm{i}\frac{\pi}{4}}}{\sqrt{8\pi}}\frac{e^{\mathrm{i}k|\mathbf{r}|}}{\sqrt{|\mathbf{r}|}}\left(\frac{a(\operatorname{\mathbf{r}}_{1})}{\sqrt{k}}e^{-\mathrm{i}k\hat{\mathbf{r}\cdot\mathbf{r}_{1}}}+A(\mathbf{r})\right),$
(2.8)

where the scattering amplitude $A$ is defined by

$A(\mathbf{r}):=k^{\frac{3}{2}}\int_{\Omega}e^{-\mathrm{i}k\hat{\mathbf{r}}\cdot\mathbf{r}^{\prime}}\eta\left(\mathbf{r}^{\prime}\right)U\left(\mathbf{r}^{\prime}\right)d\mathbf{r}^{\prime},$
(2.9)

and $\Omega$ is the volume of the scatterer.

In the weak-scattering approximation, the scattering amplitude can be
calculated by replacing $U$ with $U_{i}$. Since the scattering amplitude $A$
is recorded in the far zone at at the point $\operatorname{\mathbf{r}}_{2}$,
the scattering amplitude can be rewritten as:

$\displaystyle A(\operatorname{\mathbf{r}}_{1},\operatorname{\mathbf{r}}_{2})$
$\displaystyle=a(\operatorname{\mathbf{r}}_{1})k^{\frac{3}{2}}\int_{\Omega}e^{-\mathrm{i}k\widehat{\operatorname{\mathbf{r}}}_{2}\cdot\operatorname{\mathbf{r}}}\eta(\operatorname{\mathbf{r}})G(\operatorname{\mathbf{r}},\operatorname{\mathbf{r}}_{1})d\operatorname{\mathbf{r}},$
(2.10)

where the dependence of the scattering amplitude
$A(\operatorname{\mathbf{r}}_{1},\operatorname{\mathbf{r}}_{2})$ on the source
position $\operatorname{\mathbf{r}}_{1}$ and the detector position
$\operatorname{\mathbf{r}}_{2}$ has been made explicit. The inverse problem
thus consists of solving the integral equation (2.10) to recover the
susceptibility $\eta$.

If we differentiate (2.10) with respect to $\operatorname{\mathbf{r}}_{1}$,
use the fact that
$\displaystyle\mathop{}\\!\mathbin{\bigtriangleup}_{\mathbf{r}}G(\mathbf{r},\mathbf{r^{\prime}})+k^{2}G(\mathbf{r},\mathbf{r^{\prime}})=-\delta(\mathbf{r}-\mathbf{r^{\prime}})$
for the Green’s function, and assume unit amplitude for all sources in (2.10),
we obtain

$\displaystyle\mathop{}\\!\mathbin{\bigtriangleup}_{\operatorname{\mathbf{r}}_{1}}A(\operatorname{\mathbf{r}}_{1},\operatorname{\mathbf{r}}_{2})+k^{2}A(\operatorname{\mathbf{r}}_{1},\operatorname{\mathbf{r}}_{2})$
$\displaystyle=k^{\frac{3}{2}}\int_{\Omega}e^{-\mathrm{i}k\widehat{\operatorname{\mathbf{r}}}_{2}\cdot\operatorname{\mathbf{r}}}\eta(\operatorname{\mathbf{r}})\left(\mathop{}\\!\mathbin{\bigtriangleup}_{\operatorname{\mathbf{r}}_{1}}G(\operatorname{\mathbf{r}},\operatorname{\mathbf{r}}_{1})+k^{2}G(\operatorname{\mathbf{r}},\operatorname{\mathbf{r}}_{1})\right)d\operatorname{\mathbf{r}}$
(2.11)
$\displaystyle=k^{\frac{3}{2}}\int_{\Omega}e^{-\mathrm{i}k\widehat{\operatorname{\mathbf{r}}}_{2}\cdot\operatorname{\mathbf{r}}}\eta(\operatorname{\mathbf{r}})\left(-\delta(\operatorname{\mathbf{r}}-\operatorname{\mathbf{r}}_{1})\right)d\operatorname{\mathbf{r}}=-k^{\frac{3}{2}}e^{-\mathrm{i}k\widehat{\operatorname{\mathbf{r}}}_{2}\cdot\operatorname{\mathbf{r}}_{1}}\eta(\operatorname{\mathbf{r}}_{1}).$

We thus obtain an inversion formula for $\eta$ of the form

$\displaystyle\eta(\operatorname{\mathbf{r}}_{1})=-\frac{e^{\mathrm{i}k\widehat{\operatorname{\mathbf{r}}}_{2}\cdot\operatorname{\mathbf{r}}_{1}}}{k^{\frac{3}{2}}}\left(\mathop{}\\!\mathbin{\bigtriangleup}_{\operatorname{\mathbf{r}}_{1}}A(\operatorname{\mathbf{r}}_{1},\operatorname{\mathbf{r}}_{2})+k^{2}A(\operatorname{\mathbf{r}}_{1},\operatorname{\mathbf{r}}_{2})\right).$
(2.12)

Theoretically, reconstruction of the susceptibility $\eta$ using (2.12) only
requires the position of the sources $\operatorname{\mathbf{r}}_{1}$. However,
numerically, reconstructions of $\eta$ based on different observations of the
scattering amplitude from different detectors $\operatorname{\mathbf{r}}_{2}$
are obtained. Thus we can determine the susceptibility
$\eta(\operatorname{\mathbf{r}}_{1})$ by averaging the reconstructions:

$\eta(\operatorname{\mathbf{r}}_{1})=\frac{\sum_{\operatorname{\mathbf{r}}_{2}}\eta(\operatorname{\mathbf{r}}_{1};\operatorname{\mathbf{r}}_{2})}{n_{d}},$
(2.13)

where the dependence of $\eta$ on the detector positions is made explicit and
$n_{d}$ represents the number of detectors $\operatorname{\mathbf{r}}_{2}$.

We illustrate the proposed reconstruction method in both two-dimensional and
three dimensional scattering media (see subsection 2.3 below). To avoid any
ambiguity, we call the inversion of the susceptibility in a two-dimensional
medium as the two-dimensional reconstruction method, while the inversion in a
three-dimensional medium is called the three-dimensional reconstruction
method.

### 2.1. Numerical analysis of the forward problem

We begin with the numerical analysis of the forward problem. Since the
scattering amplitude is computed from (2.10) and there is a logarithmic
singularity in the two-dimensional Green’s function, a careful evaluation of
the integral is required. To proceed, we discretize $\Omega$ into pixels
$C_{k}\left(k=1,2,\cdots N_{\nu}\right)$ of area $h^{2}$, where $N_{\nu}$ is
the total number of pixels. Assuming unit amplitude $a=1$ for all sources, we
split the integral of (2.10) into singular- and non-singular parts:

$\displaystyle
k^{\frac{3}{2}}\int_{\Omega}e^{-\mathrm{i}k\widehat{\operatorname{\mathbf{r}}}_{2}\cdot\operatorname{\mathbf{r}}}\eta(\operatorname{\mathbf{r}})G(\operatorname{\mathbf{r}},\operatorname{\mathbf{r}}_{1})d\operatorname{\mathbf{r}}=k^{\frac{3}{2}}\sum_{k=1}^{N_{\nu}}\int_{\Omega_{k}\backslash\Omega_{s}}e^{-\mathrm{i}k\widehat{\operatorname{\mathbf{r}}}_{2}\cdot\operatorname{\mathbf{r}}}\eta(\operatorname{\mathbf{r}})G(\operatorname{\mathbf{r}},\operatorname{\mathbf{r}}_{1})d\operatorname{\mathbf{r}}+k^{\frac{3}{2}}\int_{\Omega_{s}}e^{-\mathrm{i}k\widehat{\operatorname{\mathbf{r}}}_{2}\cdot\operatorname{\mathbf{r}}}\eta(\operatorname{\mathbf{r}})G(\operatorname{\mathbf{r}},\operatorname{\mathbf{r}}_{1})d\operatorname{\mathbf{r}},$
(2.14)

where $\Omega_{s}$ is the pixel containing the singularity. We assume that
$kh\ll 1$. As the non-singular part can be computed directly, the singular
part remains to be estimated.

The Green’s function $G\left(\mathbf{r},\mathbf{r}^{\prime}\right)$ is weakly
singular at $\mathbf{r}=\mathbf{r}^{\prime}$. We will approximate the two-
dimensional Green’s function near the singularity using its asymptotic form:

$\displaystyle\begin{aligned}
G\left(\mathbf{r},\mathbf{r}^{\prime}\right)=&\frac{1}{2\pi}\ln\frac{1}{\left|\mathbf{r}-\mathbf{r}^{\prime}\right|}+\frac{\mathrm{i}}{4}-\frac{\gamma}{2\pi}-\frac{1}{2\pi}\ln\frac{k}{2}+O\left(\left|\mathbf{r}-\mathbf{r}^{\prime}\right|^{2}\ln\frac{1}{\left|\mathbf{r}-\mathbf{r}^{\prime}\right|}\right),\end{aligned}$
(2.15)

for $\left|\mathbf{r}-\mathbf{r}^{\prime}\right|\rightarrow 0$, where
$\gamma=0.5772156$ is the Euler-Mascheroni constant. Thus, the singular part
in (2.14) can be computed as

$\displaystyle\begin{aligned}
k^{\frac{3}{2}}\int_{\Omega_{s}}\eta(\mathbf{r})e^{-\mathrm{i}k\widehat{\mathbf{r}}_{2}\cdot\mathbf{r}}G\left(\mathbf{r},\mathbf{r}_{1}\right)d\mathbf{r}&\approx\frac{k^{\frac{3}{2}}}{2\pi}\eta(\widetilde{\mathbf{r}})e^{-\mathrm{i}k\hat{\mathbf{r}}_{2}\cdot\tilde{\mathbf{r}}}\left[\left(\frac{\mathrm{i}\pi}{2}-\gamma-\ln\frac{k}{2}\right)h^{2}+\int_{-\frac{h}{2}}^{\frac{h}{2}}\int_{-\frac{h}{2}}^{\frac{h}{2}}\left(\ln\frac{1}{\left|\mathbf{r}-\mathbf{r}_{1}\right|}\right)d\mathbf{r}\right]\\\
&=\eta(\widetilde{\mathbf{r}})e^{-\mathrm{i}k\hat{\mathbf{r}}_{2}\cdot\widetilde{\mathbf{r}}}\frac{k^{\frac{3}{2}}h^{2}}{2\pi}\left(\xi-\ln\frac{hk}{2}\right),\end{aligned}$

where $\displaystyle\xi=\frac{1}{2}(3+\ln
2)-\frac{\pi}{4}-\gamma+\frac{\mathrm{i}\pi}{2}$ and $\widetilde{\mathbf{r}}$
is the central point in $\Omega_{s}$.

### 2.2. Numerical analysis of the inverse problem

Next, we focus on estimating the susceptibility through differentiation of the
measured scattering amplitude $A(\mathbf{r}_{1},\mathbf{r}_{2})$ using (2.12).
We apply the theory of RKHS [1, 12], to carry out the indicated
differentiation, which allows for an explicit formulation of a functional by
means of a reproducing kernel [20, 18]. In the setting of our problem, we
identify a functional $f$ within a space $\mathcal{W}$ that approximates the
scattering amplitude $A(\mathbf{x})$. Here
$\mathbf{x}=(\mathbf{r}_{1},\mathbf{r}_{2})=\mathcal{D}\in\mathbb{R}^{4}$,
where $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$ are the position of the sources
and detectors, respectively. This leads to the following optimization problem:

$\operatorname*{argmin}_{f}\frac{1}{n}\sum_{i=1}^{n}|\bar{A}_{i}-f(\mathbf{x}_{i})|^{2}+\lambda\|f\|^{2}_{\mathcal{W}},\quad\bar{A}_{i}=A(\mathbf{x}_{i}),$
(2.16)

where $n=n_{s}\times n_{d}$, with $n_{s}$ representing the number of sources
and $n_{d}$ representing the number of detectors.

The following result (commonly referred as the representer theorem [12, 11])
demonstrates that the minimizer of problem (2.16) still lies in the space
$\mathcal{W}$ under certain conditions. This result effectively reduces the
original infinite-dimensional minimization problem (2.16) to a finite-
dimensional one for the scalar coefficients. More precisely, if the space
$\mathcal{W}$ has a reproducing kernel $\kappa:X\times
X\rightarrow\mathbb{R}$, then the minimizer $f^{*}$ of the problem (2.16) lies
in a finite-dimensional subspace of the infinite-dimensional space
$\mathcal{W}$, where

$f^{*}(\cdot)=\sum_{i=1}^{n}c^{*}_{i}\kappa(\mathbf{x}_{i},\cdot),\ \
c^{*}_{i}\in\mathbb{R},\ \ i=1,2,\dots,n,$ (2.17)

and $\mathbf{x}_{i}$ is the known input.

Let $\mathcal{W}$ be a reproducing kernel Hilbert space, $\kappa$ be its
reproducing kernel, and $K$ be the Gram matrix, with
${K}_{i,j}=\kappa(\mathbf{x}_{i},\mathbf{x}_{j})$, where
$\mathbf{x}_{i},\mathbf{x}_{j}\in\mathcal{D}$. Then any functional
$f\in\text{span}\\{\kappa(\mathbf{x}_{i},\cdot),i=1,2,\dots,n\\}$, can be
written as $f=\sum^{n}_{i=1}c_{i}\kappa(\mathbf{x}_{i},\cdot)$. Thus the
vector $(f(\mathbf{x}_{1}),\dots,f(\mathbf{x}_{n}))=\bar{c}{K}$, where
$\bar{c}=(c_{1},\cdots,c_{n})$. The components of (2.16) can be rewritten in
terms of the Gram matrix ${K}$ and the coefficient $\bar{c}$ as

$n^{-1}\sum_{i=1}^{n}|\bar{A}_{i}-f(\mathbf{x}_{i})|^{2}=n^{-1}\|\bar{A}-\bar{c}{K}\|^{2}_{\mathbb{R}^{n}},$

where $\bar{A}=(\bar{A}_{1},\dots,\bar{A}_{n})$ and

$\displaystyle\lambda\|f\|^{2}_{\mathcal{W}}$
$\displaystyle=\lambda\langle\sum^{n}_{i=1}c_{i}\kappa(\mathbf{x}_{i},\cdot),\sum^{n}_{j=1}c_{j}\kappa(\mathbf{x}_{j},\cdot)\rangle_{\mathcal{W}}=\lambda\sum^{n}_{i,j=1}c_{i}c_{j}\kappa(\mathbf{x}_{i},\mathbf{x}_{j})=\lambda\langle\bar{c}K,\bar{c}\rangle_{\mathbb{R}^{n}},$

where we have used the symmetry property of the kernel function. This allows
us to consider the optimization problem (2.16) as a function of $\bar{c}$
admitting a representation,

$\displaystyle\mathcal{M}(\bar{c})$
$\displaystyle:=n^{-1}\sum_{i=1}^{n}|\bar{A}_{i}-f(\mathbf{x}_{i})|^{2}+\lambda\|f\|^{2}_{\mathcal{W}}=n^{-1}\|\bar{A}-\bar{c}{K}\|^{2}_{\mathbb{R}^{n}}+\lambda\langle\bar{c}K,\bar{c}\rangle_{\mathbb{R}^{n}}$
(2.18)
$\displaystyle=n^{-1}\|\bar{A}\|^{2}_{\mathbb{R}^{n}}-2n^{-1}\langle\bar{c},\bar{A}K\rangle_{\mathbb{R}^{n}}+\langle\bar{c},\bar{c}(\lambda
K+n^{-1}K^{2})\rangle_{\mathbb{R}^{n}},$

where in the last equality we use the fact that $K$ is a symmetric matrix. By
the representer theorem, if the minimizer of the problem (2.16) has the form
of the equation (2.17), then the vector
$\bar{c}^{*}=(c_{1}^{*},\dots,c_{n}^{*})$ satisfies

$\nabla_{\bar{c}}\mathcal{M}(\bar{c})=-2n^{-1}\bar{A}K+2\bar{c}K(\lambda
I+n^{-1}K)=0,$ (2.19)

for $\bar{c}=\bar{c}^{*}$. Since the kernel function $\kappa$ is a positive
semidefinite function, then the Gram matrix $K$ is a positive semidefinite
matrix singular matrix. Thus, we obtain the formula for $\bar{c}^{*}$:

$\bar{c}^{*}=\bar{A}K[K(n\lambda I+K)]^{\dagger},$ (2.20)

where, $[K(n\lambda I+K)]^{\dagger}$ is the generalized inverse of the matrix
$[K(n\lambda I+K)]$.

Next, we specify the space $\mathcal{W}$. To ensure the necessary smoothness,
we choose the standard Sobolev space $W_{2}^{s}(\mathbb{R}^{d})$, where $d=4$
and $s=3$. In this case, the Sobolev embedding theorem is satsified for
$s>d/2$. This condition implies that $\mathcal{W}$ is a space of continuous
functions, and the function values are continuous linear functionals,
indicating that it is a reproducing kernel Hilbert space. The reproducing
kernel of the Sobolev space $W_{2}^{3}(\mathbb{R}^{4})$ is given by [17]:

$\kappa(\mathbf{x},\mathbf{t})=\int_{\mathbb{R}^{4}}\frac{\prod_{j=1}^{4}\cos\left(2\pi\left(x_{j}-t_{j}\right)u_{j}\right)}{(1+\sum_{0<|\alpha|_{1}\leq
3}\prod_{j=1}^{4}(2\pi u_{j})^{2\mathbf{\alpha}_{j}})}\mathrm{d}\bf u,$ (2.21)

for all
$\mathbf{x}=(x_{1},\dots,x_{4}),\mathbf{t}=(t_{1},\dots,t_{4})\in\mathbb{R}^{4}$,
where $|\mathbf{\alpha}|_{1}=\sum_{j=1}^{4}\alpha_{j}$ with non-negative
integers $\mathbf{\alpha}_{j}$.

We have now gathered all the ingredients needed to approximate the scattering
amplitude $A$. We approximate $A$ by

$f(\mathbf{t})=\sum^{n}_{i=1}\bar{c}^{*}_{i}\kappa(\mathbf{x}_{i},\mathbf{t}),$

where $\mathbf{t}=(t_{1},t_{2},t_{3},t_{4})\in\mathbb{R}^{4}$,
$\mathbf{x}_{i}\in\mathcal{D}$ and $\bar{c}^{*}$ is given by (2.20). The
reconstruction of the susceptibility $\eta$ via (2.12) requires
differentiation of
$A(\operatorname{\mathbf{r}}_{1},\operatorname{\mathbf{r}}_{2})$ with respect
to $\operatorname{\mathbf{r}}_{1}$, which in turn requires differentiation of
$f(t_{1},t_{2},t_{3},t_{4})$ with respect to the first two components $t_{1}$
and $t_{2}$, yielding

$\mathop{}\\!\mathbin{\bigtriangleup}_{t_{1},t_{2}}f(\mathbf{t})=\sum^{n}_{i=1}\bar{c}^{*}_{i}\kappa^{\prime\prime}(\mathbf{x}_{i},\mathbf{t}),$
(2.22)

where $\kappa^{\prime\prime}(\mathbf{x},\mathbf{t})$ is given by

$\kappa^{\prime\prime}(\mathbf{x},\mathbf{t})=-4\pi^{2}\int_{\mathbb{R}^{4}}\frac{(u_{1}^{2}+u_{2}^{2})\prod_{j=1}^{4}\cos\left(2\pi\left(x_{j}-t_{j}\right)u_{j}\right)}{(1+\sum_{0<|\alpha|_{1}\leq
3}\prod_{j=1}^{4}(2\pi u_{j})^{2\mathbf{\alpha}_{j}})}\mathrm{d}\bf u.$

Thus, we have identified a functional $f(\mathbf{t})$ within a RKHS
$W_{2}^{3}(\mathbb{R}^{4})$ that approximates the scattering amplitude
$A(\operatorname{\mathbf{r}}_{1},\operatorname{\mathbf{r}}_{2})$. Moreover,
the differentiation of $f$ (in turn $A$) is explicitly obtained in terms of
(2.22), which is needed for the reconstruction of $\eta$ in (2.12).

### 2.3. Three dimensional problem

We consider next the case of a three-dimensional medium. The total field $U$
still satisfies the time-harmonic wave equation (2.1) and the integral
equation (2.2), but the three-dimensional Green’s function is given by

$\displaystyle
G(\operatorname{\mathbf{r}},\operatorname{\mathbf{r}}^{\prime})=\frac{e^{\mathrm{i}k|\operatorname{\mathbf{r}}-\operatorname{\mathbf{r}}^{\prime}|}}{4\pi|\operatorname{\mathbf{r}}-\operatorname{\mathbf{r}}^{\prime}|}.$

In the far field, the Green’s function takes the asymptotic form for
$\operatorname{\mathbf{r}}\gg\operatorname{\mathbf{r}}^{\prime}$ given by,

$G\left(\mathbf{r},\mathbf{r}^{\prime}\right)=\frac{e^{\mathrm{i}k\left|\operatorname{\mathbf{r}}\right|}}{4\pi\left|\operatorname{\mathbf{r}}\right|}\left(e^{-\mathrm{i}k\hat{\operatorname{\mathbf{r}}}\cdot\operatorname{\mathbf{r}}^{\prime}}+\mathcal{O}\left(\frac{1}{\left|\operatorname{\mathbf{r}}\right|}\right)\right).$
(2.23)

It follow that the asymptotic behavior of the total field is given by,

$\displaystyle U(\mathbf{r})$
$\displaystyle=a\left(\mathbf{r}_{1}\right)G\left(\mathbf{r},\mathbf{r}_{1}\right)+k^{2}\int
G\left(\mathbf{r},\mathbf{r}^{\prime}\right)\eta\left(\mathbf{r}^{\prime}\right)U\left(\mathbf{r}^{\prime}\right)d\mathbf{r}^{\prime}$
(2.24)
$\displaystyle\sim\frac{e^{\mathrm{i}k\left|\operatorname{\mathbf{r}}\right|}}{4\pi\left|\operatorname{\mathbf{r}}\right|}e^{-\mathrm{i}k\hat{\mathbf{r}}\cdot\mathbf{r}_{1}}a\left(\mathbf{r}_{1}\right)+\frac{e^{\mathrm{i}k\left|\operatorname{\mathbf{r}}\right|}}{4\pi\left|\operatorname{\mathbf{r}}\right|}k^{2}\int
e^{-\mathrm{i}k\hat{\mathbf{r}}\cdot\mathbf{r}^{\prime}}\eta\left(\mathbf{r}^{\prime}\right)U\left(\mathbf{r}^{\prime}\right)d\mathbf{r}^{\prime}$

In the far-zone of the scatterer, the field behaves as an outgoing spherical
wave:

$U(\mathbf{r})\sim\frac{e^{\mathrm{i}k\left|\operatorname{\mathbf{r}}\right|}}{4\pi\left|\operatorname{\mathbf{r}}\right|}\left(a(\operatorname{\mathbf{r}}_{1})e^{-\mathrm{i}k\hat{\mathbf{r}\cdot\mathbf{r}_{1}}}+A(\mathbf{r})\right),$
(2.25)

where the scattering amplitude $A$ is defined by

$A(\mathbf{r})=k^{2}\int_{\Omega}e^{-\mathrm{i}k\hat{\mathbf{r}}\cdot\mathbf{r}^{\prime}}\eta\left(\mathbf{r}^{\prime}\right)U\left(\mathbf{r}^{\prime}\right)d\mathbf{r}^{\prime}.$
(2.26)

Following the same procedure as before, the scattering amplitude in three
dimensions is rewritten as,

$\displaystyle A(\operatorname{\mathbf{r}}_{1},\operatorname{\mathbf{r}}_{2})$
$\displaystyle=a(\operatorname{\mathbf{r}}_{1})k^{2}\int_{\Omega}e^{-\mathrm{i}k\widehat{\operatorname{\mathbf{r}}}_{2}\cdot\operatorname{\mathbf{r}}}\eta(\operatorname{\mathbf{r}})G(\operatorname{\mathbf{r}},\operatorname{\mathbf{r}}_{1})d\operatorname{\mathbf{r}},$
(2.27)

where we assume the sources have unit amplitude. The inversion formula for
$\eta$ is thus given by,

$\displaystyle\eta(\operatorname{\mathbf{r}}_{1})=-\frac{e^{\mathrm{i}k\widehat{\operatorname{\mathbf{r}}}_{2}\cdot\operatorname{\mathbf{r}}_{1}}}{k^{2}}\left(\mathop{}\\!\mathbin{\bigtriangleup}_{\operatorname{\mathbf{r}}_{1}}A(\operatorname{\mathbf{r}}_{1},\operatorname{\mathbf{r}}_{2})+k^{2}A(\operatorname{\mathbf{r}}_{1},\operatorname{\mathbf{r}}_{2})\right).$
(2.28)

The three-dimensional inverse problems is also considered in [9], where (2.27)
is converted into a linear system and then a regularized pseudoinverse
solution is used to recover the susceptibility. In addition, the number of
sources, detectors, and volume elements in $\Omega$ was chosen so that the
resulting linear system is overdetermined. Moreover, the two-dimensional case
was not considered in [9].

The forward problem of solving (2.26) was carried out by using the coupled-
dipole method [14, 19]. We discretized the three-dimensional domain $\Omega$
into voxels $C_{k}\left(k=1,2,\cdots N_{\nu}\right)$ of volume $h^{3}$, where
$N_{\nu}$ is the total number of voxels. We split the integral of (2.27) into
singular- and non-singular parts according to

$\displaystyle
k^{2}\int_{\Omega}e^{-\mathrm{i}k\widehat{\operatorname{\mathbf{r}}}_{2}\cdot\operatorname{\mathbf{r}}}\eta(\operatorname{\mathbf{r}})G(\operatorname{\mathbf{r}},\operatorname{\mathbf{r}}_{1})d\operatorname{\mathbf{r}}=k^{2}\sum_{k=1}^{N_{\nu}}\int_{\Omega_{k}\backslash\Omega_{s}}e^{-\mathrm{i}k\widehat{\operatorname{\mathbf{r}}}_{2}\cdot\operatorname{\mathbf{r}}}\eta(\operatorname{\mathbf{r}})G(\operatorname{\mathbf{r}},\operatorname{\mathbf{r}}_{1})d\operatorname{\mathbf{r}}+k^{2}\int_{\Omega_{s}}e^{-\mathrm{i}k\widehat{\operatorname{\mathbf{r}}}_{2}\cdot\operatorname{\mathbf{r}}}\eta(\operatorname{\mathbf{r}})G(\operatorname{\mathbf{r}},\operatorname{\mathbf{r}}_{1})d\operatorname{\mathbf{r}},$
(2.29)

where we assume that $kh\ll 1$. Combining the fact that
$e^{\mathrm{i}k\left|\operatorname{\mathbf{r}}-\operatorname{\mathbf{r}}_{1}\right|}\approx
1+\mathrm{i}k\left|\operatorname{\mathbf{r}}-\operatorname{\mathbf{r}}_{1}\right|$,
the singular part of (2.29) is computed as

$\displaystyle\begin{aligned}
\int_{\Omega_{s}}e^{-\mathrm{i}k\widehat{\mathbf{r}}_{2}\cdot\mathbf{r}}\eta(\mathbf{r})G\left(\mathbf{r},\mathbf{r}_{1}\right)d\mathbf{r}&\approx\frac{1}{4\pi}\eta(\widetilde{\mathbf{r}})e^{-\mathrm{i}k\widehat{\mathbf{r}}_{2}\cdot\widetilde{\mathbf{r}}}\int_{-\frac{h}{2}}^{\frac{h}{2}}\int_{-\frac{h}{2}}^{\frac{h}{2}}\int_{-\frac{h}{2}}^{\frac{h}{2}}\left(\frac{1}{\left|\operatorname{\mathbf{r}}-\operatorname{\mathbf{r}}_{1}\right|}+\mathrm{i}k\right)d\mathbf{r}=\frac{1}{4\pi}\eta(\widetilde{\mathbf{r}})e^{-\mathrm{i}k\widehat{\mathbf{r}}_{2}\cdot\widetilde{\mathbf{r}}}\zeta,\end{aligned}$

where $\zeta=h^{2}(\xi+\mathrm{i}kh)$ with
$\xi=\text{ln}(26+15\sqrt{3})-\frac{\pi}{2}\approx 2.38,$ and
$\widetilde{\mathbf{r}}$ is the central point in $\Omega_{s}$.

Following along similar lines as in the previous subsection, we reconstruct
the three-dimensional susceptibility $\eta$ based on (2.28). We consider the
optimization problem (2.16) in the Sobolev space $W_{2}^{s}(\mathbb{R}^{d})$,
where $d=6$ and $s=4$, which is a RKHS. The space $W_{2}^{4}(\mathbb{R}^{6})$
has the reproducing kernel given by [17]:

$\kappa(\mathbf{x},\mathbf{t})=\int_{\mathbb{R}^{6}}\frac{\prod_{j=1}^{6}\cos\left(2\pi\left(x_{j}-t_{j}\right)u_{j}\right)}{(1+\sum_{0<|\alpha|_{1}\leq
4}\prod_{j=1}^{6}(2\pi u_{j})^{2\mathbf{\alpha}_{j}})}\mathrm{d}\bf u,$ (2.30)

where
$\mathbf{x}=(\mathbf{r}_{1},\mathbf{r}_{2})=\mathcal{D}\in\mathbb{R}^{6}$ with
$\mathbf{r}_{1}$ and $\mathbf{r}_{2}$ are the source and detector positions,
respectively. Similarly, we identify a functional $f$ within the RKHS
$W_{2}^{4}(\mathbb{R}^{6})$ that approximates the scattering amplitude $A$.
The functional

$f(\mathbf{t})=\sum^{n}_{i=1}\bar{c}^{*}_{i}\kappa(\mathbf{x}_{i},\mathbf{t}),$

where $\mathbf{t}=(t_{1},t_{2},t_{3},t_{4},t_{5},t_{6})\in\mathbb{R}^{6}$,
$\mathbf{x}_{i}\in\mathcal{D}$ and $\bar{c}^{*}$ is as in (2.20). The
reconstruction of $\eta$ from (2.28) requires differentiation of
$A(\operatorname{\mathbf{r}}_{1},\operatorname{\mathbf{r}}_{2})$ with respect
to $\operatorname{\mathbf{r}}_{1}$, which in turn requires differentiation of
$f(t_{1},t_{2},t_{3},t_{4},t_{5},t_{6})$ with respect to $t_{1}$, $t_{2}$ and
$t_{3}$, yielding

$\mathop{}\\!\mathbin{\bigtriangleup}_{t_{1},t_{2},t_{3}}f(\mathbf{t})=\sum^{n}_{i=1}\bar{c}^{*}_{i}\kappa^{\prime\prime}(\mathbf{x}_{i},\mathbf{t}),$
(2.31)

where $\kappa^{\prime\prime}(\mathbf{x},\mathbf{t})$ is given by

$\displaystyle\displaystyle\kappa^{\prime\prime}(\mathbf{x},\mathbf{t})=-4\pi^{2}\int_{\mathbb{R}^{6}}\frac{(u_{1}^{2}+u_{2}^{2}+u_{3}^{2})\prod_{j=1}^{6}\cos\left(2\pi\left(x_{j}-t_{j}\right)u_{j}\right)}{(1+\sum_{0<|\alpha|_{1}\leq
4}\prod_{j=1}^{6}(2\pi u_{j})^{2\mathbf{\alpha}_{j}})}\mathrm{d}\bf u.$

Finally, following the same procedure as in the previous subsection, we
reconstruct $\eta$ using (2.28).

For three-dimensional reconstructions, we present two methods. First we
perform two-dimensional reconstructions using (2.12) for each slice, and then
make a three-dimensional reconstruction by assembling these cross-sectional
reconstructions together. The second, is to directly use (2.28) for three-
dimensional reconstruction.

## 3\. numerical experiments

To illustrate the reconstruction method, we present three numerical
experiments.

### 3.1. Reconstructions of a three-ball model

We consider reconstructing a three-ball model [9]. The model system consists
of a volume of dimensions 70 nm $\times$ 70 nm $\times$ 40 nm. Three spherical
scatterers of radius 12 nm are placed so their lowest point is 3 nm from the
bottom of the sample. The minimum distance between any two spheres is 5 nm.
The susceptibility of two of the spheres is set to 1.275, while the remaining
sphere has a susceptibility of 1.885. The wave number $k=2\pi/500=0.0126$.
Gaussian noise at the 1$\%$ level is added to the scattering amplitude and the
positions of the sources. We divide the system into 17 layers, thus each layer
is a two-dimensional medium, we solve the forward problem on each layer with a
discretization of 35$\times$35, while the susceptibility for each layer is
reconstructed on a grid of size 31$\times$31 by (2.12). This model system is
similar to what is considered in [9]. We use 150 randomly placed sources and 7
detectors of $L$ shape in this experiment, for each layer, to show the
feasibility of the two-dimensional reconstruction method.

(a) Positions of 150 sources (b) layer 1-5 and 13-17  (c) layer 6 and 12

(d) layer 7 and 11  (e) layer 8 and 10  (f) layer 9

Figure 1: Susceptibility reconstruction of three spherical scatterers on each
layer by the two-dimensional method.

Figure 1 shows the position of the sources and the reconstructions for each
layer. Due to the symmetry of the medium and the source positions for each
layer being the same, some layer reconstructions are also identical. In our
three-ball model, layer 9 is the central layer, so that layer 8 $\&$ 10, and
layer 7 $\&$ 11 (see Figure 1 (d) and (e) ) yield identical layer
reconstructions. Thus, one can reconstruct on either of the two layers, hence
reducing the reconstruction time. It is also worth noting that the Gram matrix
$K$ for calculating the vector $\bar{c}^{*}$ in (2.20) for each layer is the
same, since the position of the sources and the reconstruction grid are the
same for each layer. This also applies to subsection 3.2.1.

### 3.2. Reconstructions of a neuron

We consider a neuron from a rat brain, of size 3.928$\mu m$$\times$5.44$\mu m$
$\times$2.248$\mu m$. The image was obtained from experiment and manually
segmented [15, 16]. The main trunk of the neuron has a susceptibility of
$\eta=0.25$ and the branches spines have a susceptibility of $\eta=0.38$,
while the background has susceptibility $\eta=0$. Once again, 1$\%$ Gaussian
noise is added to the scattering amplitude and the positions of the sources.
Both two-dimensional method and three-dimensional methods were tested.

####  3.2.1 Two-dimensional reconstruction method

We manually segmented the neuron into a grid of size 33$\times$46$\times$18,
then we calculated the forward problem on 18 layers with each layer of size
33$\times$46\. Subsequently, we reconstructed each of these 18 layers onto a
31$\times$44 grid, denoted as $\mathcal{G}:=\\{(g_{i},g_{j})\\}$, with
$i=1,\dots,31$, and $j=1,\dots,44$. For each layer, 500 randomly spaced
sources and 7 detectors were employed.

Since we approximate the scattering amplitude by
$f(\mathbf{t})=\sum^{n}_{i=1}\bar{c}^{*}_{i}\kappa(\mathbf{x}_{i},\mathbf{t})$,
$\mathbf{t}=(t_{1},t_{2},t_{3},t_{4})\in\mathbb{R}^{4}$,
$\mathbf{x}_{i}\in\mathcal{D}$, we can approximate the value at grid points by
letting $(t_{1},t_{2})=(g_{i},g_{j})$ and randomly set $(t_{3},t_{4})$. In
Figure 2, we only set one $(t_{3},t_{4})$, which is the coordinate of one
detector, yielding one susceptibility reconstruction for each grid point.

(a) Positions of 500 sources (b) layer 1 and 18 (c) layer 2  (d) layer 3  (e)
layer 4  (f) layer 5

(g) layer 6 (h) layer 7 (i) layer 8  (j) layer 9  (k) layer 10  (l) layer 11

(m) layer 12 (n) layer 13 (o) layer 14  (p) layer 15  (q) layer 16  (r) layer
17

Figure 2: Neuron model and susceptibility reconstructions on each layer by the
two-dimensional method.

Then putting all of the layers together, we obtain the reconstruction shown in
Figure 3. To better display the three dimensional colour map, we add
transparency to the points of value less than 0.085.

Figure 3: Neuron reconstructions by two-dimensional method compared to the
groundtruth.

Let us define the relative error for the scattering amplitude:

$\chi^{2}=\frac{\sum_{{\operatorname{\mathbf{r}}}_{1},{\operatorname{\mathbf{r}}}_{2}}\left|A\left({\operatorname{\mathbf{r}}}_{1},{\operatorname{\mathbf{r}}}_{2}\right)-A_{\mathrm{rec}}\left({\operatorname{\mathbf{r}}}_{1},{\operatorname{\mathbf{r}}}_{2}\right)\right|^{2}}{\sum_{{\operatorname{\mathbf{r}}}_{1},{\operatorname{\mathbf{r}}}_{2}}\left|A\left({\operatorname{\mathbf{r}}}_{1},{\operatorname{\mathbf{r}}}_{2}\right)\right|^{2}},$

where $A_{\mathrm{rec}}$ is the reconstructed scattering amplitude. We define
the mean relative error of the susceptibility at grid $\mathcal{G}$ as

$\delta=\frac{\sum^{31}_{i=1}\sum^{44}_{j=1}|\eta\\{g_{i},g_{j}\\}-\eta_{rec}\\{g_{i},g_{j}\\}|}{31\times
44},$

where, $\eta\\{g_{i},g_{j}\\}$ is the groundtruth susceptibility at grid point
$(g_{i},g_{j})$, and $\eta_{rec}$ is the reconstructed susceptibility.

It is important to determine the optimal number of detectors, sources and
$(t_{3},t_{4})$ for a given experiment. In Figure 4 we employ 7 detectors and
show the behavior of the error with different numbers of sources on the
reconstructions of layer 9. We set $(t_{3},t_{4})$ equal to the coordinates of
7 detectors, which leads to multiple reconstructions of susceptibility. Then
the mean of the multiple reconstructions compared to the single
$(t_{3},t_{4})$ reconstruction is plotted.

Figure 4: Error behavior with different numbers of sources

We see that 500 sources almost reaches the minimum of the susceptibility
error. Moreover, the reconstructed susceptibility error by a single
$(t_{3},t_{4})$ is comparable to multiple detectors. Thus, Figure 4 shows that
one can choose 500 sources and single $(t_{3},t_{4})$.

Next, in Figure 5, we fix 500 sources and choose a single $(t_{3},t_{4})$
reconstruction for each grid point, we then test the impact of different
numbers of detectors on the reconstructions. We find that the number of
detectors does not significantly affect the susceptibility reconstructions. In
all settings (Figure 4 and Figure 5), the scattering amplitude error is on the
order of $10^{-6}$. This shows good performance of using the RKHS method. Note
that the reconstructions depicted in Figure 2 have been obtained using the
same configuration of 500 sources, 7 detectors, and single $(t_{3},t_{4})$ as
those also utilized in Figure 4 and Figure 5.

Figure 5: Error behavior using different numbers of detectors

####  3.2.2 Three-dimensional reconstruction method

Here we show the feasibility of our three-dimensional reconstruction method
and compare the reconstructions with known results. We segment the neuron into
a grid of size 33$\times$46$\times$19, and then we solve the forward problem
on this grid with 500 randomly spaced sources and 7 detectors. Subsequently,
we reconstruct the model on the grid 33$\times$44$\times$18, thus, avoiding an
inverse crime. The result is shown in Figure 6, where the purple dots indicate
the positions of the sources. Once again, we add transparency to the points of
value less than 0.08 to better display the colour map.

Figure 6: Neuron reconstructions by the three-dimensional method compared to
groundtruth.

We next compared our RKHS method with the method in [9]. We consider the same
problem setting as in [9] and compute the corresponding scattering amplitude
error results in Table LABEL:tb:comp-john. We see that for setup in [9], the
scattering amplitude error by the proposed method is less than reported in
[9]. Moreover, one can reduce the number of detectors from 225 to 7 and still
obtain the same scattering amplitude error. Thus the proposed RKHS method is
superior to the method presented in [9].

Table 1: Comparison between our RKHS method and the one reported in [9].

| Linear system [9] | RKHS method | RKHS method
---|---|---|---
Number of sources | 500 | 500 | 500
Number of detectors | 225 | 225 | 7
Discretization of the forward problem | 33$\times$46$\times$19 | 33$\times$46$\times$19 | 33$\times$46$\times$19
Scattering amplitude error | $1*10^{-2}$ | $1.9*10^{-4}$ | $1*10^{-4}$

## 4\. Conclusion

We have developed a method to reconstruct the dielectric susceptibility of an
inhomogeneous scattering medium from internal sources using a RKHS approach.
Our numerical results in two- and three dimensions problems exhibit better
performance compared to previous work in [9]. We note that the approach we
have presented is quite general and applicable to imaging with any scalar wave
field with internal sources.

### Acknowledgements

This research was funded in whole, or in part, by the Austrian Science Fund
(FWF) 10.55776/P34981 (OS $\&$ YD) – New Inverse Problems of Super-Resolved
Microscopy (NIPSUM), SFB 10.55776/F68 (OS) “Tomography Across the Scales,”
project F6807-N36 (Tomography with Uncertainties). The financial support by
the Austrian Federal Ministry for Digital and Economic Affairs, the National
Foundation for Research, Technology and Development and the Christian Doppler
Research Association is gratefully acknowledged. The work of JCS was supported
by NSF grant DMS-1912821 and AFOSR grant FA9550-19-1-0320. For the purpose of
open access, the author has applied a CC BY public copyright license to any
Author Accepted Manuscript version arising from this submission. The
computational results presented have been achieved using the Vienna Scientific
Cluster (VSC).

## References

## References

  * [1] N. Aronszajn  “Theory of reproducing kernels”  In _Transactions of the American Mathematical Society_ 68.3, 1950, pp. 337–404  DOI: 10.1090/s0002-9947-1950-0051437-7
  * [2] E. Betzig et al.  “Imaging intracellular fluorescent proteins at nanometer resolution”  In _Science_ 313.5793, 2006, pp. 1642–1645  DOI: 10.1126/science.1127344
  * [3] M. Born and E. Wolf  “Principles of Optics”  Cambridge: Cambridge University Press, 1999
  * [4] R. Carminati and J.. Schotland  “Principles of Scattering and Transport of Light”  Cambridge University Press, 2021  DOI: 10.1017/9781316544693
  * [5] D. Colton and R. Kress  “Inverse acoustic and electromagnetic scattering theory” 93, Applied Mathematical Sciences  Berlin: Springer-Verlag, 1998  URL: https://www.springer.com/gp/book/9783662035375
  * [6] H.. Engl, M. Hanke and A. Neubauer  “Regularization of inverse problems”, Mathematics and its Applications 375  Dordrecht: Kluwer Academic Publishers Group, 1996  URL: https://link.springer.com/book/9780792341574
  * [7] C.. Epstein  “Introduction to the Mathematics of Medical Imaging”  New York: Prentice Hall, 2003
  * [8] C. Errico et al.  “Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging”  In _Nature_ 527.7579, 2015, pp. 499–502  DOI: 10.1038/nature16066
  * [9] A.. Gilbert, H.. Levinson and J.. Schotland  “Imaging from the inside out: inverse scattering with photoactivated internal sources”  In _Optics Letters_ 43.12, 2018, pp. 3005  DOI: 10.1364/ol.43.003005
  * [10] S.. Hell et al.  “The 2015 super-resolution microscopy roadmap”  In _Journal of Physics D: Applied Physics_ 48.44, 2015, pp. 443001  DOI: 10.1088/0022-3727/48/44/443001
  * [11] G. Kimeldorf and G. Wahba  “Some results on Tchebycheffian spline functions”  In _Journal of Mathematical Analysis and Applications_ 33.1, 1971, pp. 82–95  DOI: 10.1016/0022-247x(71)90184-3
  * [12] G.. Kimeldorf and G. Wahba  “A Correspondence Between Bayesian Estimation on Stochastic Processes and Smoothing by Splines”  In _The Annals of Mathematical Statistics_ 41.2, 1970, pp. 495–502  DOI: 10.1214/aoms/1177697089
  * [13] P. Kuchment  “The Radon Transform and Medical Imaging”, CBMS-NSF Regional Conference Series in Applied Mathematics  Philadelphia: SIAM, 2013  DOI: 10.1137/1.9781611973297
  * [14] H.. Levinson and V.. Markel  “Solution of the nonlinear inverse scattering problem by T-matrix completion. I. Theory”  In _Physical Review E_ 94.4, 2016  DOI: 10.1103/physreve.94.043317
  * [15] M.. Martone et al.  “Correlated 3D Light and Electron Microscopy: Use of High Voltage Electron Microscopy and Electron Tomography for Imaging Large Biological Structures”  In _Journal of Histotechnology_ 23.3, 2000, pp. 261–270  DOI: 10.1179/his.2000.23.3.261
  * [16] M.. Martone et al.  “A cell-centered database for electron tomographic data”  In _Journal of Structural Biology_ 138.1-2, 2002, pp. 145–155  DOI: 10.1016/s1047-8477(02)00006-0
  * [17] E. Novak, M. Ullrich, H. Woźniakowski and S. Zhang  “Reproducing kernels of Sobolev spaces on $\mathbb{R}^{d}$ and applications to embedding constants and tractability”  In _Analysis and Applications_ 16.5, 2018, pp. 693–715  DOI: 10.1142/s0219530518500094
  * [18] S. Pereverzyev  “An Introduction to Artificial Intelligence Based on Reproducing Kernel Hilbert Spaces”, 2022  DOI: 10.1007/978-3-030-98316-1
  * [19] E.. Purcell and C.. Pennypacker  “Scattering and Absorption of Light by Nonspherical Dielectric Grains”  In _The Astrophysical Journal_ 186, 1973, pp. 705–714  DOI: 10.1086/152538
  * [20] S. Saitoh and Y. Sawano  “Theory of Reproducing Kernels and Applications”, 2016  DOI: 10.1007/978-981-10-0530-5
  * [21] A.. Tikhonov and V.. Arsenin  “Solutions of Ill-Posed Problems”  Washington, D.C.: John Wiley & Sons, 1977
  * [22] A.. Weglein et al.  “Inverse scattering series and seismic exploration”  In _Inverse Problems_ 19.6, 2003, pp. R27–R83  DOI: 10.1088/0266-5611/19/6/r01


# Temperature-induced reversal effects of kink dynamics in carbon nanotube on
flat substrate

Alexander V. Savin  N.N. Semenov Federal Research Center for Chemical Physics,
Russian Academy of Science (FRCCP RAS), Moscow, 119991, Russia  Plekhanov
Russian University of Economics, Moscow, 117997 Russia  Margarita Kovaleva
N.N. Semenov Federal Research Center for Chemical Physics, Russian Academy of
Science (FRCCP RAS), Moscow, 119991, Russia  HSE University, Moscow, 109028
Russia

###### Abstract

Carbon nanotubes are nano-objects with quite anisotropic properties, for
example the mechanical properties in longitudinal and radial directions differ
significantly. This feature of the carbon nanotubes yields many interesting
phenomena investigated in last decades. One of them is the ability to form
both hollow and collapsed states if the radius of the nanotube is large
enough. The transitions between the two states have been also reported. In our
study we present single-walled carbon nanotube interacting with a plane
substrate and characterize the energy of interaction with the substrate using
effective Lennard-Jones-type potential. We show energy of the homogeneous open
and collapsed states depending on the radius of the carbon nanotube and report
on the bi-stability in some range of the nanotube diameters. Using the
molecular-dynamical simulations we look at the evolution of the initial half-
opened, half-collapsed state and demonstrate that the transition area from one
state to another is spatially localized having features of topological soliton
(kink or anti-kink). We show that the value and the direction of the kink
propagation speed depend significantly on the nanotube diameter as well as on
the temperature of the system. We also discuss the mechanism of the process
using a simplified model with asymmetric double-well potential and show the
entropic nature of the transition.

###### pacs:

44.10.+i, 05.45.-a, 05.60.-k, 05.70.Ln

##  I Introduction

Carbon nanotubes (CNTs) are cylindrical macromolecules with a diameter varying
from half a nanometer up to 20 nanometers. They are long, hollow tubule
structures made of graphene sheets. Similar structures were obtained firstly
more than 70 years ago during the thermal decomposition of carbon monoxide on
an iron contact Radushkevich52 . However, CNTs themselves were synthesized
only 30 years ago as by-products of fullerene C60 synthesis Iijima91 . CNTs
are promising engineering nano-materials with increasing usage and
significance in nanotechnology Eleckii02 . CNTs attract interest of
researchers in physics, material science, electronics and biotechnology and
nanotechnology due to their unique thermal, mechanical, optical and biological
properties Ferreira16 ; Menezes19 ; Venkataraman19 . From the mechanical point
of view nanotube is a quasi-one-dimensional molecular structure with
pronounced nonlinear properties Dresselhaus01 ; Astakhova01 ; Savin04 ;
Savin08 .

Nanotubes have high longitudinal (axial) and relatively weak transverse
(radial) stiffness. Because of this, at sufficiently large diameters,
nanotubes due to the weak non-covalent interaction of atoms can transform from
a hollow cylindrical shape to a collapsed state Chopra95 ; Gao98 ; Xiao07 ;
Baimova15 ; Xu2019 ; Grande2020 ; Magnin2021 ; Umeno . Non-valent interaction
with the substrate can also change the cylindrical shape of the nanotube
Hertel98 ; Xie10 ; Yuan18 . It has been shown that long multi-walled narrow
graphene nanoribbons can be created by squashing carbon nanotubes using a
thermally assisted high-pressure process Chen21 ; Toh21 . Such collapsed
nanotubes can be used as semiconducting graphene nanoribbons.

We report on the collapse of the carbon nanotube on the substrate, and
underline that this process is controlled by the temperature of the system.
Moreover, we show that the effect is reversal, which means that the nanotube
that collapsed due to temperature decrease can open back if the temperature is
increased again. The front of the opening or collapsing in the semi-collapsed
nanotube has the soliton-like profile, which can diffusively move in the
longitudinal direction. Using a simplified (biparabolic) effective model for
the evolution of the front, we show that Brownian motion of a kink in
asymmetric double-well potential can describe the process. Using this approach
we compare the theoretically obtained speed of the solitary wave with the
speed of the front from the molecular-dynamical investigation and obtain a
fairly good agreement in the low-temperature range.

##  II Model

We consider a CNT with chirality indices $(m,m)$. The cylindrical structure of
such a nanotube is formed by periodic repetition along the $x$ axis of
transverse cyclic zigzag chains of $K=4m$ carbon atoms:

$\displaystyle x_{n,(j-1)m+i}$ $\displaystyle=$ $\displaystyle h_{i}+a(n-1),$
$\displaystyle y_{n,(j-1)m+i}$ $\displaystyle=$ $\displaystyle
R\cos(\phi_{i}+(j-1)\Delta\phi),$ $\displaystyle z_{n,(j-1)m+i}$
$\displaystyle=$ $\displaystyle R\sin(\phi_{i}+(j-1)\Delta\phi),$

where the first index $n=0,\pm 1,\pm 2,...$ numbers the transverse rings of
atoms (unit cells), the second index $k=(j-1)m+i$, $j=1,...,m$, $i=1,2,3,4$, –
atoms in these rings. Here the angular pitch $\Delta\phi=2\pi/m$, the radius
of the nanotube $R=r_{0}/2\sin(\Delta\phi/6)$ ($r_{0}=1.418$Å is the
equilibrium length of the C–C valence bond in a graphene sheet), the
longitudinal pitch $a=r_{0}\sqrt{3}$, the longitudinal displacements
$h_{1}=h_{4}=0$, $h_{2}=h_{3}=a/2$, the angular displacements $\phi_{1}=0$,
$\phi_{2}=\Delta\phi/6$, $\phi_{3}=3\Delta\phi/6$, $\phi_{4}=4\Delta\phi/6$.

Let the index $n$ number the transverse ring elements of the CNT with the
armchair structure, then the positions of the atoms of each ring can be
described by a vector of $3K$ dimension ${\bf u}_{n}=\\{{\bf
v}_{n,k}=(x_{n,k},y_{n,k},z_{n,k})\\}_{k=1}^{K}$. The Hamiltonian of a
nanotube has the following form

$H=\sum_{n=-\infty}^{+\infty}\sum_{k=1}^{K}\left\\{\frac{1}{2}M(\dot{\bf
v}_{n,k},\dot{\bf v}_{n,k})+P_{n,k}+V(z_{n,k})\right\\},$ (1)

where $M$ is the mass of a carbon atom, ${\bf
v}_{\alpha}=\left(x_{\alpha}(t),y_{\alpha}(t),z_{\alpha}(t)\right)$ is a
vector defining the position of the atom with two-dimensional index
$\alpha=(n,k)$ at time $t$. Term $P_{\alpha}$ describes the interaction of the
atom with the index $\alpha$ with the other atoms of the nanotube, the last
term of $V(z_{\alpha})$ sets the energy of interaction of the atom with flat
substrate.

Figure 1: View of a carbon nanotube with a chirality index (31,31) located on
a flat substrate formed by the surface of a silicon carbide crystal 6HSiC
(0001). The left end of the nanotube is in an open stationary state (the
cross-section of the nanotube has the shape of a convex drop), the right end
is in a collapsed stationary state (the cross-section has the shape of an
asymmetric dumbbell with a flat two-layer central part). In the middle part of
the nanotube, a localized region of its smooth transition from one stable
stationary state to another is formed.

To describe the carbon-carbon valence interactions we use a standard set of
molecular dynamics potentials Savin10 . The valence bond between two
neighboring carbon atoms $\alpha$ and $\beta$ can be described by the Morse
potential

$U_{1}({\bf v}_{\alpha},{\bf
v}_{\beta})=\epsilon_{1}\\{\exp[-\alpha_{0}(r_{\alpha\beta}-r_{0})]-1\\}^{2},$
(2)

where $r_{\alpha\beta}=|{\bf v}_{\alpha}-{\bf v}_{\beta}|$,
$\epsilon_{1}=4.9632$ eV is the valence bond energy and $r_{0}=1.418$ Å is the
equilibrium valence bond length.

Valence angle deformation energy between three adjacent carbon atoms $\alpha$,
$\beta$, and $\gamma$ can be described by the potential

$U_{2}({\bf v}_{\alpha},{\bf v}_{\beta},{\bf
v}_{\gamma})=\epsilon_{2}(\cos\varphi-\cos\varphi_{0})^{2},$ (3)

where $\cos\varphi=({\bf v}_{\alpha}-{\bf v}_{\beta},{\bf v}_{\gamma}-{\bf
v}_{\beta})/r_{\alpha\beta}r_{\gamma\beta}$, and $\varphi_{0}=2\pi/3$ is the
equilibrium valent angle. Parameters $\alpha_{0}=1.7889$ Å-1 and
$\epsilon_{2}=1.3143$ eV can be found from the small amplitude oscillations
spectrum of the graphene sheet Savin08 .

Valence bonds between four adjacent carbon atoms $\alpha$, $\beta$, $\gamma$,
and $\delta$ constitute the torsion angles, the potential energy of which can
be defined as

$U_{3}(\phi)=\epsilon_{3}(1-\cos\phi),$ (4)

where $\phi$ is the corresponding torsion angle ($\phi_{0}=0$ is the
equilibrium value of the angle) and $\epsilon_{3}=0.499$ eV Gunlycke08 .

More detailed discussion and motivation of our choice of the interaction
potentials (2), (3), (4) can be found in our earlier publication Savin10 .

Non-valent interactions of carbons atoms are described Setton96  by the
Lennard-Jones potential

$V_{cc}(r)=4\epsilon_{cc}[(\sigma_{cc}/r)^{12}-(\sigma_{cc}/r)^{6}],$ (5)

where $\epsilon_{cc}=0.002757$ eV is the binding energy and
$\sigma_{cc}=3.393$ Å. The potential attains a minimum value of
$-\epsilon_{cc}$ at $r_{cc}=2^{1/6}\sigma_{cc}=3.807$ Å (equilibrium
interatomic distance).

To simulate the dynamics of a nanotube located on a flat substrate formed by
the surface of a molecular crystal, it is necessary to find the interaction
potential of the carbon atom with the substrate. For this purpose, the
interaction energy of a finite flat sheet of graphene with a flat surface of a
6H-SiC(0001) silicon carbide crystal was found Sforzini15 . The calculations
used a graphene sheet size $1.985\times 1.643$ nm2, consisting of 160 carbon
atoms, located parallel to the surface of the crystal $z=0$ at a distance of
$h$. The interaction energy of each carbon atom with a silicon carbide crystal
was calculated as the sum of the Lennard-Jones potentials (5) with the values
of the parameters from Rappe92 . At each value of $h$, the energy of the
interaction of the sheet with the crystal was averaged along the shift along
the $x$ and $y$ axes, and then normalized by the number of atoms in the sheet.
As a result, the dependence of the interaction energy of one atom of the sheet
on its distance to the substrate plane $V(h)$ was obtained.

The calculations showed that the interaction energy with the substrate $V(h)$
can be described with high accuracy by the $(k,l)$ Lennard-Jones potential

$V(h)=\epsilon_{0}[k(h_{0}/h)^{l}-l(h_{0}/h)^{k}]/(l-k),$ (6)

where degree $k=3.75$, $l=17$, binding energy $\epsilon_{0}=0.073$ eV,
equilibrium distance to plane $h_{0}=4.19$ Å.

##  III Homogeneous stationary states of the nanotube

To find a uniform stationary state of the nanotube it is necessary to solve
the problem of the minimum potential energy

$E=\sum_{k=1}^{K}\\{P_{(0,k)}+V(z_{(0,k)})\\}\rightarrow\min:\\{{\bf
v}_{0,k}\\}_{k=1}^{K},\leavevmode\nobreak\ a$ (7)

by the coordinates of atoms in one unit cell $(n=0)$ and the value of the
longitudinal period. Problem (7) was solved numerically by the conjugate
gradient method. By choosing the initial configuration of atoms in the cell,
all stable stationary states of the nanotube can be obtained.

The solution of the problem (7) showed that at the index $m<17$ (at the radius
of the nanotube $R<11.517$ Å), the nanotube $(m,m)$ on a flat substrate has
only one stable stationary state (open ground state), in which the cross-
section of the nanotube has the form of a convex drop – see Fig. 1. At $m\geq
17$, in addition to the open state, there is also a second stable stationary
state (collapsed ground state), in which the cross-section has the form of an
asymmetric dumbbell with a flat two-layer central part.

Figure 2: Dependence of (a) normalized energy $E/K$ and (b) diameters $D_{y}$,
$D_{z}$ nanotubes with chirality index $(m,m)$ from the value of the index
$m$. Curves 1, 2, 3 give the dependence for an open primary state, and curves
4, 5, 6 – the represent the ground state of the collapsed nanotube located on
a flat substrate of silicon carbide 6H-SiC(0001). Line 7 gives the diameter
dependence for an isolated cylindrical nanotube.

The stationary state of the nanotube $(\\{{\bf v}_{0,k}^{0}\\}_{k=1}^{K},a)$
is characterized by its normalized energy $E/K$, transversal and vertical
diameters

$D_{y}=\max_{k_{1},k_{2}}|y_{0,k_{1}}^{0}-y_{0,k_{2}}^{0}|,\leavevmode\nobreak\
\leavevmode\nobreak\
D_{z}=\max_{k_{1},k_{2}}|z_{0,k_{1}}^{0}-z_{0,k_{2}}^{0}|.$

The dependencies of the normalized energy $E/K$ and diameters $D_{y}$, $D_{z}$
on the value of the index $m$ of the nanotube located on a flat substrate of
crystalline silicon carbide are presented in Fig. 2. As can be seen from the
figure, on the substrate the vertical diameter of the nanotube is always
smaller than the transverse diameter ($D_{z}<D_{y}$) – interaction with the
substrate leads to flattening of the nanotube. The transverse diameter of the
stationary state always grows in proportion to the index value: $D_{y}\sim m$
for $m\rightarrow\infty$. The vertical diameter of the open state of the
nanotube grows as the logarithm of the index: $D_{z}\sim\ln m$, and the
vertical diameter of the collapsed state remains constant $D_{z}\approx 8.07$
Å (the diameter value is determined by the height of the end curves of the
cross-section of the nanotube).

At low values of the index $17\leq m\leq 30$, the open state is more
advantageous in energy than the collapsed state of the nanotube. If $m\geq 31$
the ground state becomes the collapsed one. With an increase in the index $m$,
the difference in the energy of the states monotonically decreases at $m<31$,
reaches a minimum value $\Delta E=|E_{o}-E_{c}|=0.014$ eV at $m=31$ (open
state energy value is $E_{o}=-3.909$ eV and collapsed state energy value is
$E_{c}=-3.923$ eV) and monotonically increases at $m>31$ [see Fig. 2 (a)].

Figure 3: Dependence of the energy $E$ of the nanotube $(m,m)$ located on the
substrate from the flat surface of the silicon carbide crystal 6H-SiC (0001)
on the distance $h$ from the substrate to the center of the opposite side of
the nanotube at different values of the index $m$. The energy is counted from
the minimum value.

##  IV Energy profile of the transition between the stationary states

The results obtained allow us to conclude that a single-walled nanotube with
an index $m\geq 16$ is a bistable system with two stable (open and collapsed)
states. The energy profile of the transition between these two states can be
found numerically; for that we solve the problem on the minimum of (7) for
each fixed value of the distance $h$ between the substrate plane and the
center of the upper side of the nanotube. The numerical solution of this
problem allows us to obtain the dependence $E(h)$, describing the change in
the energy of the nanotube during its homogeneous collapse.

The view of potential $E(h)$ at different values of the index of the nanotube
$m$ is shown in Fig. 3. As can be seen from the figure at $m<17$, the function
$E(h)$ is a one-well potential with a minimum corresponding to the open state.
At $m=17$, the potential obtains a two-well form, it has a new narrow minimum
corresponding to the collapsed stationary state. With the increase of the
index $m$ (with the increase of the radius) of the nanotube, the depth of the
new narrow minimum monotonically increases. At $m\geq 31$ this minimum becomes
energetically preferred, i.e. the collapsed state becomes the main state of
the CNT. In our case the energy profile $E(h)$ has the form of a strongly
asymmetric double-well potential with a narrow first valley and a wide second
valley.

The transition of a nanotube from one stationary state to another can be
described qualitatively as the motion of a kink (topological soliton) in the
$\phi$–4 model with an asymmetric double-well potential $E(h)$ having one
narrow deep well corresponding to the collapsed state and a second wide well
with higher energy corresponding to the open state of the nanotube. This case
is considered in detail in Savin97 ; Costantini01 , where it is shown that the
direction of motion of the kink in such a chain depends on the temperature
value. The kink motion describes the sequential transition of a chain from a
non-ground state to a ground state. At low temperatures, the main state will
always be in a deeper narrow valley, at high temperatures – a higher energy
minimum but in wider second valley. Switching the ground state with an
increase in temperature leads to a change in the direction of motion of the
topological soliton. Let’s model this effect for motion along the nanotube
topological soliton describing its sequential transition from one ground state
to another – see Fig. 1. For an isolated nanotube, the effect of changing the
direction of motion of this topological soliton with temperature changes was
detected in Chang10 , but did not receive sufficient explanation.

Figure 4: View of the upper part of the longitudinal section of the nanotube
$(m,m)$ with a topological soliton at index (a) $m=31$ and (b) $m=32$. The
view of the section is shown at thermostat temperatures $T=0,30,...,930$ K at
time $t=400$ and $t=500$ ps. the dotted line shows the position of the center
of the soliton at the initial time value.

##  V Temperature dependence of the direction of motion of a topological
soliton

Figure 5: View of the upper part of the longitudinal section of the nanotube
$(m,m)$ with a topological soliton at index (a) $m=33$ and (b) $m=34$. The
view of the section is shown at thermostat temperatures $T=0,30,...,930$ K at
time $t=300$ and $t=200$ ps. The dotted line shows the position of the center
of the soliton at the initial time value.

To model the motion of a topological soliton in a thermalized nanotube, we
consider a finite nanotube with a chirality index $(m,m)$ of $N=420$ cross
elements. At the initial moment of time, the left part of the nanotube $1\leq
n\leq 180$ will be transferred to the hollow stationary state, the right part
$N-180<n\leq N$ – to the collapsed state, and in the central part $180<n\leq
N-180$ we will define a linear continuous transition from one state to
another. We fix the position of atoms in the end cells of the nanotube $n=1,2$
and $n=N-1,N$ and consider the dynamics of the nanotube immersed in the
Langevin thermostat.

The dynamics of a thermalized nanotube is described by a system of Langevin
equations

$\displaystyle M\ddot{\bf v}_{n,k}=-\frac{\partial H}{\partial{\bf
v}_{n,k}}-\Gamma M\dot{\bf v}_{n,k}+\Xi_{n,k},$ (8) $\displaystyle
n=3,4,...,N-2,\leavevmode\nobreak\ \leavevmode\nobreak\ k=1,2,...,K,$

where the vector ${\bf v}_{n,k}=(x_{n,k},y_{n,k},z_{n,k})$ specifies the
coordinates of the $k$-th atom in the $n$-th transversal element, $H$ is the
Hamiltonian (1), $\Gamma=1/t_{r}$ is the friction coefficient (velocity
relaxation time is $t_{r}=1$ ps), and
$\Xi_{n,k}=(\xi_{n,k,1},\xi_{n,k,2},\xi_{n,k,3})$ is three-dimensional vector
of the normally distributed random forces, normalized by the conditions

$\langle\xi_{n_{1}k_{1}j_{1}}\\!(t_{1})\xi_{n_{2}k_{2}j_{2}}\\!(t_{2})\rangle\\!=\\!2M\Gamma
k_{B}T\delta_{n_{1}n_{2}}\delta_{k_{1}k_{2}}\delta_{j_{1}j_{2}}\delta(t_{1}-t_{2})$

($k_{B}$ is the Boltzmann constant, $T$ is thermostat temperature).

If thermostat temperature is $T=0$, than initially given topological defect in
nanotube, with the center at the node $N/2$ will always move to the right,
translating the main part of the nanotube in the hollow state, when the index
$m\leq 30$. Same way, the defect will move to the left, shifting the nanotube
to close state if the condition $m>30$ holds. In this case, a localized region
is formed in the nanotube in which there is a smooth transition from one
stationary state to another – see Fig. 1. This region moves at a constant
speed while maintaining its shape, i.e. behaves like a topological soliton.
The upper part of the central longitudinal section of the nanotube has the
shape of an anti-kink $\\{z_{n}\\}_{n=1}^{N}$. The motion of a topological
soliton along a nanotube is conveniently described as the motion of this one-
dimensional anti-kink.

Figure 6: View of the bistable potential with index $m=31$ $E(h)$ obtained
from the energy form in the MD simulations (curve 1) and its approximation by
a bi-parabolic potential $W(h)$ (curve 2) with parameters
$h_{1},h_{0},h_{2}=0.3316$, 0.5058, 2.2473 nm, $K_{1},K_{2}=50$, 0.5 eV/nm2.

Simulation of nanotube dynamics shows that in the nanotube $(31,31)$ the
topological soliton moves to the left at $T<90$K and to the right at $T>90$K –
see Fig. 4 (a). Here, the temperature at which the direction of motion of the
topological soliton is switched $T_{0}=90$K (the collapsed state of the
nanotube is the main one at $T<T_{0}$, and the open state – at $T>T_{0}$). In
the nanotube (32,32), the switching of the direction of motion of the soliton
occurs already at a higher temperature $T_{0}=600$K – see Fig. 4 (b). For
wider nanotubes (33,33) and (34,34), switching must occur at very high
temperatures, so at $T\leq 930$K, the topological soliton always moves to the
left, bringing the nanotube into a collapsed state. The velocity of the
soliton thus decreases with increasing temperature – see Fig. 5.

##  VI Mechanism of the front speed temperature dependence

The amplitudes of the deformations of the nanotube during the collapse from
the hollow state are sufficiently high. Description of such evolution using,
for example, nonlinear Sanders-Koiter thin shell theory should give us
nonlinear equations which turn out to be too complicated for analysis, as too
many modes of deformations need to be taken into account Smirnov16 .
Therefore, we suggest an effective model representing the evolution of CNT to
describe the process of the collapse and a vice-versa process.

To describe the mechanism of the temperature dependence we consider the
nanotube as coupled unit-cells (rings of atoms) in a bistable potential
defined by the properties of the nanotube and by the interaction with the
substrate. Let us represent the motion of each element as an effective
particle in a non-degenerate double-well potential – see Fig. 3. Than we can
model the evolution of the CNT with a chain in a bistable on-site potential
$E(h)$. The energy difference between the two potential minima can be defined
as $\Delta E$. For nanotube with chirality index (31,31) the two-well
potential $E(h)$ can be approximated by a bi-parabolic potential

$W(h)=\begin{cases}\frac{1}{2}K_{1}(h-h_{1})^{2}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{if $h\leq h_{0}$,}\\\
\frac{1}{2}K_{2}(h-h_{2})^{2}\leavevmode\nobreak\ \leavevmode\nobreak\
\text{if $h>h_{0}$.}\end{cases}$ (9)

where potential minima are $h_{1}=3.316$ Å, $h_{2}=22.473$ Å, maximum point is
$h_{0}=5.058$ Å, stiffness values are $K_{1}=0.5$ eV/Å2, $K_{2}=0.005$ eV/Å2 –
see Fig. 6, $\Delta E=0$.

Half-openned (half-collapsed) state of the nanotube (31,31) may be described
as kink (antikink) in one-dimensional chain with double-well substrate
potential (9). Hamiltonian of this chain

$H=\sum_{n}\frac{1}{2}\mu\dot{h}_{n}^{2}+\frac{1}{2}\kappa(h_{n+1}-h_{n})^{2}+W(h_{n}),$
(10)

where effective mass of one chain cell $\mu\approx mM$, $m$ is chirality of
the nanotube, $M$ is mass of carbon atom, (only about 1/4 part of the carbon
atoms in each transversal cell participate in the collapse of the nanotube).

For further analysis it is convenient to introduce the dimensionless
displacements $u_{n}=-1+(h_{n}-h_{1})/h_{d}$, $h_{d}=(h_{2}-h_{1})/2$ and
dimensionless energy $E=H/E_{max}$, where $E_{max}=K_{1}(h_{0}-h_{1})^{2}/2$.
Then dimensionless double-well potential $V(u)=W(h_{1}+(u+1)h_{d})$ can be
represented in the form:

$V(u)=\begin{cases}\frac{1}{2}k_{1}(u+1)^{2}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{if $u\leq u_{0}$,}\\\
\frac{1}{2}k_{2}(u-1)^{2}\leavevmode\nobreak\ \leavevmode\nobreak\ \text{if
$u>u_{0}$,}\end{cases}$ (11)

where $u_{0}=-1+2(h_{0}-h_{1})/(h_{2}-h_{1})=-0.8181$,
$k_{i}=(h_{2}-h_{1})^{2}K_{i}/4E_{0}$, $i=1,2$, $k_{1}=60.468$,
$k_{2}=0.60468$, see Fig. 7.

We suppose that the evolution along the CNT is smooth enough to apply the
long-wave approximation:

$\displaystyle u_{n}=u,$ $\displaystyle
u_{n+1}=u+a_{0}u_{x}+\frac{a_{0}^{2}}{2}u_{xx},$ (12) $\displaystyle
u_{n-1}=u-a_{0}u_{x}+\frac{a_{0}^{2}}{2}u_{xx},$

where $a_{0}=1$ is dimensionless distance between the two unit-cells, $x$ is a
dimentionless longitudinal coordinate scaled by the longitudinal pitch
$x_{0}=a$. Under these assumptions we can model the dynamics of
circumferential motion of the CNT as an equation:

$u_{tt}-c_{0}^{2}u_{xx}+\omega_{0}^{2}\bar{V}^{\prime}(u)=-\gamma
u_{t}+\zeta(x,t),$ (13)

where the $u$ is an effective coordinate having sense of the state coordinate,
which can be defined by the coordinate of center of mass of the moving part of
nanotube, $c_{0}=\sqrt{\kappa a^{2}/\mu}$ is a longitudinal wave-speed. In
general case the potential is non-degenerate, $u_{1}$ and $u_{2}$ are the
local minima coordinates with the energy gap
$\Delta\epsilon=|\bar{V}(u_{1})-\bar{V}(u_{2})|$. For nanotube with chirality
(31,31) the potential in the dimensionless form $\bar{V}(u)$ can be taken from
(11) as $\bar{V}(u)=V(u)$, and $\omega_{0}=\sqrt{E_{max}/(\mu h_{d}^{2})}$ is
characteristic frequency of the motion, which yields:

$\bar{V}(u)=\begin{cases}\omega_{1}^{2}(u-u_{1})^{2}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{if $u\leq u_{0}$,}\\\
\omega_{2}^{2}(u-u_{2})^{2}\leavevmode\nobreak\ \leavevmode\nobreak\ \text{if
$u>u_{0}$,}\end{cases}$ (14)

$u_{1}=-1$, $u_{2}=1$, $\omega_{1}=\sqrt{k_{1}/2}$ and
$\omega_{2}=\sqrt{k_{2}/2}$ are local curvatures of the potential, which
define frequencies of the motion in the vicinity of each minimum. For nanotube
with chirality (31,31) the characteristic frequency of the motion
$\omega_{0}=1.148\cdot 10^{10}$ 1/s, $\omega_{1}=0.02273$,
$\omega_{2}=0.00642$, see Fig. 6. We neglect the energy gap between the two
minima, for this case we take $\Delta\epsilon=0$. Let us note, that the more
precise evaluation of the analytically obtained form of the potential $V(u)$
does not significantly change the results of the analysis. Therefore we keep
simple bi-parabolic approximation.

Figure 7: View of the bistable dimensionless potential V(u), $u_{0}=-0.8181$,
$k_{1}=60.468$, $k_{2}=0.60468$.

The r.h.s of (13) defines the coupling to the equilibrium heat bath with a
viscous damping constant $\gamma$ and the Gaussian noise $\zeta(x,t)$ with
zero mean and autocorrelation function:

$\langle\zeta(x,t)\zeta(x^{\prime},t^{\prime})\rangle=2\gamma
D\delta(x-x^{\prime})\delta(t-t^{\prime}),$

where $D=k_{B}T/A$, $A=h_{d}\mu$. In our calculations we used value
$\gamma=\Gamma$ from the previous section. The unperturbed kink for the
degenerate potential $V(u)$ when $\Delta\epsilon=0$ in its implicit form can
be found as:

$x-vt=\left(1-\frac{v^{2}}{c_{0}^{2}}\right)^{1/2}\frac{d}{\sqrt{2}}\int_{u_{0}}^{u_{1,2}(x-vt)}\frac{du}{\sqrt{\bar{V}(u)}},$
(15)

where $d=c_{0}/\omega_{0}$ is size of quasi-particle, constant speed is $v$,
$|v|<c_{0}$, $u_{0}$ – maximum point of the potential
[$\bar{V}(u_{0})=\max\limits_{u_{1}<u<u_{2}}V(u)$], the rest energy of the
soliton is:

$E_{0}=\omega_{0}c_{0}A\int_{u_{1}}^{u_{2}}\sqrt{2V(u)}du,$ (16)

rest mass is $M_{0}=E_{0}/c_{0}^{2}$. The anti-kink is defined from (15) with
reversed sign of the right hand size. The kink solution describes the
transition of the nanotube profile from the collapsed state to the open state,
i. e. $u_{1,2}(x\rightarrow-\infty)=u_{1}$,
$u_{1,2}(x\rightarrow+\infty)=u_{2}$, while the anti-kink solution corresponds
to transition from the open state to the collapsed one,
$u_{2,1}(x\rightarrow-\infty)=u_{2}$, $u_{2,1}(x\rightarrow+\infty)=u_{1}$.
The initial conditions of the numerical simulations correspond to the anti-
kink evolution.

The main factors defining the anti-kink velocity and its direction is the form
of the two wells of the potential and the energy gap $\Delta\epsilon$. The
phenomenon has an entropic nature and is defined by the difference of profile
of the both wells. To illustrate this phenomenon, let us suppose that each
element of our initial discrete chain is subject to the bistable potential.
The initial condition will correspond to the placement of the right side of
the chain to the left of the potential barrier, while the left side rests to
the right of the barrier. We can calculate the probability of the existence of
the element in each of potential wells. If the probability to find the element
in the r.h.s. well is higher, than wave-front of the switch between the two
states will move to the right, and vice versa. Following the approach
described in the Costantini01  we define the wave-speed for the almost
degenerate case, when there energy gap is negligibly small. let us denote
characteristic frequencies of the motion in the left and right valleys of the
potential $\Omega_{1}=\omega_{1}\omega_{0}$,
$\Omega_{2}=\omega_{2}\omega_{0}$. In the frequency range
$(\Omega_{2},\Omega_{1})$ phonons yield net pressure of the antikink to the
right, while the phonons from the left side of the antikink are reflected
back. The result can be characterized by effective thermal force, which plays
its role when the local curvatures of potential close to both minima are
different. Using the approach Currie80 , we can characterize the motion of the
kink or anti-kink as a particle which undergoes Brownian motion with
additional thermal force term in the Langevin equation:

$M_{0}\ddot{X}=-\gamma M_{0}\dot{X}\mp(\Delta\epsilon+A_{th})+\xi(t),$ (17)

where $X(t)$ is coordinate of center of mass of kink or antikink
quasiparticle, $M_{0}$ is its effective mass, $A_{th}$ represents effect of
the thermal forcing, and $\xi(t)$ is a Gaussian noise with zero mean and
autocorrelation function
$\langle\xi(t)\xi(0)\rangle=2M_{0}k_{B}T/x_{0}\delta(t)$.

To provide calculations we use approximate double-well potential $V(u)$ (see
Fig.( 6)) for the case of $m=31$. Supposing that the two minima of potential
correspond to the same energy value, we obtain the relation for the speed of
the anti-kink solution:

$v_{th}=c_{0}\frac{kT}{2E_{0}}\frac{(\omega_{1}-\omega_{2})\omega_{0}}{\gamma},$
(18)

where the energy of the quasi-particle is defined in (16). Let us note, that
the result is independent on $c_{0}$, as the value $E_{0}$ is proportional to
$c_{0}$. Using this fact let us simplify the relation as follows:

$v_{th}=\frac{kT}{2IA\omega_{0}}\frac{(\Omega_{1}-\Omega_{2})}{\gamma},$ (19)

where constant
$I=E_{0}/(A\omega_{0}c_{0})=\int_{u_{1}}^{u_{2}}\sqrt{2\bar{V}(u)}du$ is
defined by the potential profile, $\Omega_{1}$, $\Omega_{2}$ are two
characteristic frequencies of the motion in the vicinity of the stationary
states, $\gamma$ is the friction parameter, term $A=h_{d}\mu$ defines the
scale of the system.

Comparison of the analytical speed-temperature dependence with results
obtained numerically from the simulation of the nanotube with chirality
(31,31) is presented in Fig. 8. Details of the analysis are presented in
Appendix. We see a fairly good correspondence of the speed dependence,
however, the starting point is not correct. This is due to the difference
between the energy values at the two minima of the potential, which was
neglected in our calculations. However, such a good accordance of the obtained
analytically speed with the numerical results allows us to conclude that the
considered phenomenon has entropic nature.

Figure 8: Speed of the soliton for temperature for the CNT index $m=31$
obtained in MD simulation (markers 1) and comparison with the theoretical
value form the degenerated bi-parabolic approximation (line 2).

##  VII Conclusions

In our work we consider the carbon nanotube on a surface with interaction
realized via effective Lennard-Jones potential. We considered the evolution of
the nanotube initially in the state with one end collapsed and another end
opened. The profile of the nanotube demonstrates the soliton-like transition
area from the opened to the collapsed states. During the MD simulations we
show the evolution of the front as a soliton-like localized wave. We show that
the direction of the soliton (kink or anti-kink) depends on the radius of the
nanotube. We also demonstrate, that the temperature can affect significantly
the wave speed of the soliton. We find the energy profile of the nanotube
depending on the coordinate of the upper middle-point of the nanotube for
different index values of the nanotube.

Using the effective model of the nanotube on the plane substrate connected
with thermostat we explain the dependence of the nanotube evolution as
effective chain of elements in a bistable potential. We prove the entropic
nature of the speed dependence on the temperature of the system. The
analytical results obtained in our asymptotic analysis sufficiently well for
such a sketch model correspond to those obtained via MD simulation.

## Acknowledgements

This work was supported by Program of Fundamental Researches of the Russian
Academy of Sciences (project no. FFZE-2022-0009, state registration no.
122040500069-7). Computational facilities for the work were provided by the
Joint SuperComputer Center of the Russian Academy of Sciences.

##  Appendix A Thermal speed of the soliton

Here we present the detailed analysis of kink and anti-kink dynamics. The
small deformations of the soliton (15) can be studied in linear approximation
and expanded on the basis of eigenfunctions $\chi(x,t)=\psi(x)e^{-i\omega t}$:

$-c_{0}^{2}\psi_{xx}+\tilde{V}_{0}(x)\psi=\omega^{2}\psi.$ (20)

This equation can be represented as a one-dimensional Schrodinger equation
with asymmetric potential-well
$\tilde{V}_{0}(x)=\omega_{0}^{2}\bar{V}(u_{12}^{(0)}(x))$, where the
$u_{1,2}^{(0)}$ is a static kink solution. The phonon modes with
$\Omega_{1}^{2}>\omega^{2}>\Omega_{2}^{2}$ exert an effective pressure on the
kink to the left, while for the anti-kink the pressure acts to the right. To
estimate the effective thermal force $F_{th}$, acting to the kink or anti-kink
let us consider the diluted gas of kinks and anti-kinks with additional
stopping tilt $F$ Costantini01 . Using the transfer-matrix formalism Currie80
we obtain the free energy of the gas in thermal equilibrium as a corresponding
eigenvalue problem:

$\displaystyle\hat{H}(u)\eta_{n}(u)=\epsilon_{n}(u),$
$\displaystyle\hat{H}(u)=-\frac{\hbar^{2}}{2m_{1}}\frac{d^{2}}{du^{2}}+\bar{V}(u)+V_{1}-\frac{F}{\omega_{0}^{2}}u,$
(21)

where $m_{1}=(\hbar\omega_{0}c_{0}A/k_{B}T)^{2}$ is an effective mass and
$V_{1}$ is the energy offset depending on temperature. The difference between
the lowest eigenvalues, i.e. the frequencies defining the curvature in the two
wells of the potential $\omega_{1,2}$ can be compensated by the stopping tilt:

$\displaystyle\frac{d_{u}F}{\omega_{0}^{2}}=\frac{\hbar}{2\sqrt{m_{1}}}\left(\sqrt{\bar{V}^{\prime\prime}(u_{2})}-\sqrt{\bar{V}^{\prime\prime}(u_{1})}\right)$
$\displaystyle=-k_{B}T\omega_{0}(\omega_{1}-\omega_{2})/(2Ac_{0}\omega_{0}^{2}),$
(22)

where $d_{u}=u_{2}-u_{1}$ is distance between the two minima. Consequently,
the internal asymmetry of the double-well potential on the kink dynamics can
be estimated as the action of the external tilt $F_{th}$:

$A_{th}=d_{u}F_{th}=\frac{k_{B}T\Delta\omega}{2Ad},$ (23)

where $\Delta\omega=\omega_{2}-\omega_{1}$. Using this fact, we suppose that
the evolution of asymmetric kink or anti-kink can be described by the Langevin
equation as evolution of the Brownian particle.

$M_{0}\ddot{X}=-\gamma M_{0}\dot{X}\mp(\Delta\epsilon+A_{th})+\xi(t),$ (24)

where $X(t)$ is the coordinate of center of mass. In the degenerate case
energy difference between the two minima of potential $\Delta\epsilon=0$, and
only the thermal effective force drives the kink or antikink with the
stationary speed to the left or to the right, consequently, with the speed:

$\frac{v_{th}}{c_{0}}=\mp\frac{k_{B}T}{2E_{0}}\frac{\Delta\omega\omega_{0}}{\gamma}$
(25)

## References

  * (1) L. Radushkevich and V. Lukyanovich, “On the structure of carbon formed during the thermal decomposition of carbon monoxide on an iron contact(russian),” Russian Journal of Physical Chemistry, vol. 26, pp. 88–95, 1952.
  * (2) S. Iijima, “Helical microtubules of graphitic carbon,” Nature, vol. 354, pp. 56–58, 1991.
  * (3) A. Eletskii, “Carbon nanotubes and their emission properties,” Phys. Usp., vol. 45, pp. 369–402, 2002.
  * (4) F. Ferreira, L. Cividanes, F. Brito, B. Menezes, and W. Franceschi, ch. 1, pp. 1–29.  Springer, 2016.
  * (5) B. de Menezes, K. Rodrigues, B. da Silva Fonseca, R. Ribas, T. do Amaral Montanheiro, and G. Thim, “Recent advances in the use of carbon nanotubes as smart biomaterials,” J. Mater. Chem. B, vol. 7, pp. 1343–1360, 2019.
  * (6) A. Venkataraman, E. Amadi, and Y. e. a. Y. Chen and, “Carbon nanotube assembly and integration for applications,” Nanoscale Res Lett, vol. 14, pp. 220–241, 2019.
  * (7) M. Dresselhaus, G. Dresselhaus, P. Avouris, and R. Smalley, Carbon Nanotubes: Synthesis, Structure, Properties and Applications.  Springer, 2001.
  * (8) T. Astakhova, O. Gurin, M. Menon, and G. Vinogradov, “Longitudinal solitons in carbon nanotubes,” Phys. Rev. B, vol. 64, p. 035418, 2001.
  * (9) A. Savin and O. Savina, “Nonlinear dynamics of carbon molecular lattices: Soliton plane waves in graphite layers and supersonic acoustic solitons in nanotubes,” Phys. Solid State, vol. 46, pp. 383–391, 2004.
  * (10) A. Savin and Y. Kivshar, “Discrete breathers in carbon nanotubes,” Europhys. Lett., vol. 82, p. 66002, 2008.
  * (11) N. Chopra, L. Benedict, V. Crespi, M. Cohen, S. Louie, and A. Zettl, “Fully collapsed carbon nanotubes,” Nature, vol. 377, pp. 135–138, 1995.
  * (12) G. Gao, T. Çağin, and W. G. III, “Energetics, structure, mechanical and vibrational properties of single-walled carbon nanotubes,” Nanotechnology, vol. 9, pp. 184–191, 1998.
  * (13) J. Xiao, B. Liu, Y. Huang, J. Zuo, K.-C. Hwang, and M.-F. Yu., “Collapse and stability of single- and multi-wall carbon nanotube,” Nanotechnology, vol. 18, p. 395703, 2007.
  * (14) J. Baimova, Q. Fan, L. Zeng, Z. Wang, S. Dmitriev, X. Feng, and K. Zhou, “Atomic structure and energy distribution of collapsed carbon nanotubes of different chiralities,” Journal of Nanomaterials, vol. 2015, p. 186231, 2015.
  * (15) H. Xu, G. Drozdov, B. Hourahine, J.G. Park, R. Sweat, T. Frauenheim, and T. Dumitric, “Collapsed Carbon Nanotubes: From Nano to Mesoscale via Density Functional Theory-Based Tight-Binding Objective Molecular Modeling,” Carbon, vol. 143,pp. 786–792,2019.
  * (16) R.R. Del Grande, Al.F. Fonseca, and R.B. Capaz, “Energy barriers for collapsing large-diameter carbon nanotubes,” Carbon, vol. 159,pp. 161–165,2020.
  * (17) Y. Magnin, F. Rondepierre, W. Cui, D.J. Dunstan, A. San-Miguel “Collapse phase diagram of carbon nanotubes with arbitrary number of walls. Collapse modes and macroscopic analog,” Carbon, vol. 178, pp. 552–562,2021.
  * (18) Y. Umeno, Y. Yachi, M. Sato, H. Shima T. Hertel, R. Walkup, and P. Avouris, “On the atomistic energetics of carbon nanotube collapse from AIREBO potential,” Physica E , vol. 106, pp. 319–325, 2019.
  * (19) T. Hertel, R. Walkup, and P. Avouris, “Deformation of carbon nanotubes by surface van der waals forces,” Phys. Rev. B, vol. 58, pp. 13870–13873, 1998.
  * (20) J. Xie, Q. Xue, H. Chen, D. Xia, C. Lv, and M. Ma, “Influence of solid surface and functional group on the collapse of carbon nanotubes,” J. Phys. Chem. C, vol. 114, pp. 2100–2107, 2018.
  * (21) X. Yuan and Y. Wang, “Collapsed adhesion of carbon nanotubes on silicon substrates: continuum mechanics and atomistic simulations,” Nanotechnology, vol. 29, pp. 184–191, 2018.
  * (22) C. Chen, Y. Lin, W. Zhou, et al. ”Sub-10-nm graphene nanoribbons with atomically smooth edges from squashed carbon nanotubes.” Nat Electron 4, 653-663 (2021).
  * (23) CT. Toh, ”Squashing carbon nanotubes into nanoribbons.” Nat Electron 4, 633-634 (2021).
  * (24) D. Gunlycke, H. Lawler, and C. White, “Lattice vibrations in single-wall carbon nanotubes,” Phys. Rev. B, vol. 77, p. 195422, 2008.
  * (25) A. Savin, Y. Kivshar, and B. Hu, “Suppression of thermal conductivity in graphene nanoribbons with rough edges,” Phys. Rev. B, vol. 82, p. 195422, 2010.
  * (26) R. Setton, “Carbon nanotubes – ii. cohesion and formation energy of cylindrical nanotubes,” Carbon, vol. 34, p. 69, 1996.
  * (27) J. Sforzini, L. Nemec, T. Denig, B. Stadtmuller, T.-L. Lee, C. Kumpf, S. Soubatch, U. S. P. Rinke, V. Blum, F. Bocquet, and F. Tautz, “Approaching truly freestanding graphene: The structure of hydrogen-intercalated graphene on 6h-sic(0001),” Phys. Rev. Lett., vol. 114, p. 106804, 2015.
  * (28) A. Rappe, C. Casewit, K. Colwell, W. G. III, and W. Skiff, “Uff, a full periodic table force field for molecular mechanics and molecular dynamics simulations,” J. Am. Chem. Soc., vol. 114, pp. 10024–10035, 1992.
  * (29) A. Savin, G. Tsironis, and A. Zolotaryuk, “Reversal effects in stochastic kink dynamics,” Phys. Rev. E, vol. 56, pp. 2457–2466, 1997.
  * (30) G. Costantini and F. Marchesoni, “Asymmetric kinks: Stabilization by entropic forces,” Phys. Rev. Lett., vol. 87, p. 114102, 2001.
  * (31) T. Chang and Z. Guo, “Temperature-induced reversible dominoes in carbon nanotubes,” Nano Lett., vol. 10, pp. 3490–3493, 2010.
  * (32) L.I.Manevitch, V.V.Smirnov, M.Strozzi, and F.Pellicano, “Nonlinear optical vibrations of single-walled carbon nanotubes.,” Physica D, vol. 325, pp. 113–125, 2016.
  * (33) J. F. Currie, J. A. Krumhansl, A. R. Bishop, and S. E. Trullinger, “Statistical mechanics of one-dimensional solitary-wave-bearing scalar fields: Exact results and ideal-gas phenomenology,” Phys. Rev., vol. 2, pp. 477–496, 1980.