# Multi-soliton solutions of the sine-Gordon equation with elliptic-function
background

Daisuke A. Takahashi1,2777The current primary affiliation is 2. 1Research and
Education Center for Natural Sciences, Keio University, Hiyoshi 4-1-1,
Yokohama, Kanagawa 223-8521, Japan
2Department of Physics, Chuo University, 1-13-27 Kasuga, Bunkyo-ku, Tokyo
112-8551, Japan<EMAIL_ADDRESS>
###### Abstract

The multi-soliton solution of the sine-Gordon equation in the presence of
elliptic-function background is derived by the inverse scattering method. The
key tool in our formulation is the Lax pair written by $4\times 4$ matrix
differential operators given by Takhtadzhyan and Faddeev in 1974, which
enables us to use the conventional form of the integral representation of the
Jost solutions and Krichever’s theory of commuting differential operators. As
a by-product we also provide generalized orthogonality and completeness
relations for eigenfunctions associated with indefinite inner product. The
multi-soliton solution is expressed by a determinant of theta functions and
the shift of the background lattice due to solitons is also determined using
addition formula. One kink and one breather solutions are presented by
animated gifs.

##  1 Introduction

The sine-Gordon (SG) equation has a wide variety of applications in physics
and other natural sciences. For example, it appears as a continuum limit of
the Frenkel-Kontorowa model describing the commensurate-incommensurate
transition (Ref. [1], chapter 6 and Ref. [2], chapter 10), the effective model
for the chiral magnet and soliton propagation on it [3, 4], and kinks in the
long Josephson junction of superconductors [5, 6, 7]. The model where the
partial derivative is changed from hyperbolic to elliptic type has also been
solved and applied to vortices and dislocations in spatially two-dimensional
systems [8, 9, 10]. The rogue-wave solutions with unstable background has been
recently investigated [11]. From the view of the theory of classical
integrable systems, the SG equation is one of the earliest equations
formulated via the zero-curvature expression, known as the Ablowitz-Kaup-
Newell-Segur (AKNS) formalism [12, 13]. The analogical integrable models with
higher group symmetries have been also investigated (Ref. [14] and references
found in Ref. [15], part II, chapter I, §8).
In some of the above-mentioned works, the multi-soliton excited states not
only for uniform background but also for oscillating non-uniform background
attract considerable physical attention. While the multi-soliton solutions of
an integrable equation in the presence of the elliptic-function, or more
general quasiperiodic Riemann theta function background, can be in principle
obtained by taking a special limit of the finite-zone quasiperiodic solutions
[16], the procedure is often complicated. Furthermore, in physical
applications, the eigenfunctions of the Lax pair, which is necessary in
construction of soliton solutions, often play additional roles in, e.g.,
calculation of dispersion relations of linearized waves and linear stability
analysis. Therefore, constructing the soliton solutions with non-uniform
background not from a degenerate case of quasiperiodic solutions but by adding
solitons to the stationary background using eigenfunctions of the Lax pair via
various conventional methods provides helpful by-products in investigation of
physical phenomena.
Motivated by these circumstances, in this paper, we derive the multi-soliton
solution of the SG equation in the presence of elliptic-function background.
Though this problem itself has been solved in many references including the
above-mentioned ones, we believe that the following features (i)-(iii) will
bring some new light. (i) We use the Lax formalism based on $4\times 4$ matrix
differential operator by Takhtadzhyan and Faddeev [17], which is a variant of
the Lax pair written by an integral operator [18]. While the famous treatment
of the SG equation in soliton theory is the zero-curvature expression using
$2\times 2$ matrices which depend on the spectral parameter [12, 13], the use
of this method has the following advantages: (a) the theory of commuting
differential operators [19] can be applied and (b) the common integral
representation of Jost solutions can be used without modification so we do not
need trial and error to find suitable form. (ii) The final expression of the
multi-soliton solution, which is obtained by formulating the inverse
scattering method (ISM) and solving the Gelfand-Levitan-Marchenko (GLM)
equation, is compactly summarized as a determinant of theta functions, and its
asymptotic form is also derived using the addition formulas in Ref. [20].
(iii) The orthogonal and completeness relations of eigenfunctions arising from
the indefinite inner product is discussed in detail, which is necessary in
formulation of the ISM and determination of the possible emerging patterns of
discrete eigenvalues corresponding to solitons. This treatment is regarded as
a generalization of the corresponding finite-dimensional linear algebra [21]
to continuous space and differential operators.
The organization of this paper is as follows. In Sec. 2, we write down the Lax
pair of the SG equation. In Sec. 3, we point out that the Lax operators are
“$\sigma$-self-adjoint”, and introduce the associated indefinite inner product
which defines the generalized orthogonal relations to eigenfunctions. In Sec.
4, we derive simultaneous eigenfunctions of the Lax pair for the stationary
soliton lattice solution, and identify the corresponding algebraic curve,
based on the theory of commuting differential operators [19]. The eigenvalues
and eigenfunctions are parametrized by uniformization variable on torus, and
their symmetries are also presented. In Sec. 5, we define left and right Jost
solutions and the scattering matrix for the system where the background
potential asymptotically tends to the stationary soliton lattice at spatial
infinities. In Sec. 6, we introduce the integral representation of the Jost
solution, and present the GLM equation which determines the kernel function of
the integral representation from the scattering data. In Sec. 7, we solve the
GLM equation for reflectionless case and determine the multi-soliton solution,
which is expressed by determinant of theta functions. The phase shift of the
background lattice is also found. In Sec. 8, by solving the time-dependent Lax
equation, we determine the time-evolution of the scattering matrix, and using
this, we derive the time-dependent multi-soliton solution for the SG equation.
In Sec. 9, we provide the constraint that the discrete eigenvalues must
satisfy in order for the resultant solution to become real and bounded. In
Sec. 10, we address the gif animations of the soliton solutions generated by
Mathematica. The method of visualization and used numerical values are
presented. Sec. 11 is devoted to summary and discussion. The appendices
discuss several technical details. In A, the integration formula necessary to
obtain the eigenfunction is presented. In B, we show the detailed derivation
of the eigenfunctions given in Sec. 4. In C, we derive the completeness
relation for eigenfunctions of the soliton lattice potential. There, we also
discuss several technical important points on, e.g., the zero-norm
eigenfunction defined by indefinite inner product and the necessity of the
expression by meromorphic integrand. In D, we provide the detailed derivation
of the GLM equation. In E, we calculate the determinant of theta functions
appearing in the asymptotic form of the soliton solution using addition
formulas.

##  2 Lax pair

We write an $n\times n$ identity and zero matrix as $I_{n}$ and $O_{n}$. Let
$\sigma_{1},\sigma_{2},$ and $\sigma_{3}$ be the Pauli matrices. The Lax pair
and the Lax equation for the SG equation is given by [17]

$\displaystyle 4\mathrm{i}\frac{\partial\hat{L}}{\partial t}$
$\displaystyle=[\hat{L},\hat{B}],$ (2.1) $\displaystyle\hat{L}$
$\displaystyle=-4\mathrm{i}\partial_{x}\begin{pmatrix}\sigma_{3}&\\\
&O_{2}\end{pmatrix}+\begin{pmatrix}\mathrm{i}\sigma_{1}w&\mathrm{e}^{-\mathrm{i}\phi\sigma_{2}/2}\\\
\mathrm{e}^{-\mathrm{i}\phi\sigma_{2}/2}&O_{2}\end{pmatrix},\quad
w=\phi_{x}+\phi_{t},$ (2.2) $\displaystyle\hat{B}$
$\displaystyle=4\mathrm{i}\partial_{x}\begin{pmatrix}I_{2}&\\\
&-I_{2}\end{pmatrix}-2\begin{pmatrix}&\sigma_{3}\mathrm{e}^{-\mathrm{i}\phi\sigma_{2}/2}\\\
\mathrm{e}^{-\mathrm{i}\phi\sigma_{2}/2}\sigma_{3}&\end{pmatrix}.$ (2.3)

Equation (2.1) then reduces to the SG equation

$\displaystyle\phi_{tt}-\phi_{xx}+\sin\phi=0.$ (2.4)

Here we slightly changed prefactors of operators and choice of the SU(2) basis
in the way
$(\sigma_{1},\sigma_{2},\sigma_{3})\to(\sigma_{1},\sigma_{3},-\sigma_{2})$
from the original work Ref. [17]. This makes no essential difference in
formulation.

##  3 Symmetry of the Lax operator and orthogonality of eigenfunctions

Henceforth we simply write
$\sigma:=\sigma_{3}\oplus\sigma_{3}=I_{2}\otimes\sigma_{3}$. Then, $\hat{L}$
and $\hat{B}$ are both “$\sigma$-self-adjoint”[21], i.e.,

$\displaystyle\sigma\hat{L}^{\dagger}\sigma=\hat{L},\quad\sigma\hat{B}^{\dagger}\sigma=\hat{B}.$
(3.1)

Therefore the orthogonal and completeness relations of eigenfunctions of these
operators are defined through the indefinite inner product, called
“$\sigma$-inner product” in [21] : $(f_{1},f_{2})_{\sigma}:=\int
dxf_{1}^{\dagger}\sigma f_{2}$. The correct identification of the completeness
relation is important in derivation of the GLM equation (C and D).
Following the established procedure of the ISM, we first consider the
eigenvalue problem of $\hat{L}$. Since the highest-order coefficient matrix of
$\hat{L}$ is given by $\sigma_{3}\oplus O_{2}$, which obviously has rank 2,
there exist two linearly independent eigenfunctions for a given eigenvalue
$\lambda$. They can be written in the form $f=\left(\begin{smallmatrix}g\\\
\lambda^{-1}\mathrm{e}^{-\mathrm{i}\phi\sigma_{2}/2}g\end{smallmatrix}\right)$
with $g$ a two-component column vector. If we rewrite the equation with
respect to $g$ using the light-cone coordinate $x^{\prime}=\frac{x+t}{2},\
t^{\prime}=\frac{x-t}{2}$, the famous $2\times 2$ AKNS form [12] is
reproduced. Let $f_{1}=\left(\begin{smallmatrix}g_{1}\\\
\frac{1}{\lambda}\mathrm{e}^{-\mathrm{i}\phi\sigma_{2}/2}g_{1}\end{smallmatrix}\right)$
and $f_{2}=\left(\begin{smallmatrix}g_{2}\\\
\frac{1}{\lambda}\mathrm{e}^{-\mathrm{i}\phi\sigma_{2}/2}g_{2}\end{smallmatrix}\right)$
be eigenfunctions of $\hat{L}$ with eigenvalues $\lambda_{1}$ and
$\lambda_{2}$. Then, we can show

$\displaystyle 4\mathrm{i}\left(f_{1}^{\dagger}(I_{2}\oplus
O_{2})f_{2}\right)_{x}=(\lambda_{1}^{*}-\lambda_{2})f_{1}^{\dagger}\sigma
f_{2}.$ (3.2)

Therefore, if $\lambda_{1}^{*}\neq\lambda_{2}$, then we have $\int
dxf_{1}^{\dagger}\sigma f_{2}=0$, which is the above-mentioned orthogonality.
On the other hand, when $\lambda_{1}^{*}=\lambda_{2}$, we obtain
$f_{1}^{\dagger}(I_{2}\oplus
O_{2})f_{2}=g_{1}^{\dagger}g_{2}=(x\text{-independent})$. We can also show
that if $\lambda_{1}=\lambda_{2}$, $\det(g_{1},g_{2})$ does not depend on $x$,
which can be more easily shown using the AKNS form.
In addition to Eq. (3.1), $\hat{L}$ further has the following symmetry:

$\displaystyle(I_{2}\otimes\sigma_{2})\hat{L}(I_{2}\otimes\sigma_{2})=\hat{L}^{*},\quad(\sigma_{3}\otimes\sigma_{2})\hat{L}(\sigma_{3}\otimes\sigma_{2})=-\hat{L},$
(3.3)

which immediately means

$\displaystyle\hat{L}f$ $\displaystyle=\lambda f$
$\displaystyle\quad\leftrightarrow\quad\hat{L}[(I_{2}\otimes\sigma_{2})f^{*}]$
$\displaystyle=\lambda^{*}(I_{2}\otimes\sigma_{2})f^{*}$
$\displaystyle\leftrightarrow\quad$
$\displaystyle\hat{L}[(\sigma_{3}\otimes\sigma_{2})f]$
$\displaystyle=-\lambda(\sigma_{3}\otimes\sigma_{2})f$
$\displaystyle\quad\leftrightarrow\quad\hat{L}[(\sigma_{3}\otimes
I_{2})f^{*}]$ $\displaystyle=-\lambda^{*}(\sigma_{3}\otimes I_{2})f^{*}.$
(3.4)

Thus the eigenvalues $\lambda,\lambda^{*},-\lambda,-\lambda^{*}$ always appear
simultaneously.
$\hat{L}$ also has a little unfamiliar symmetry:

$\displaystyle\hat{L}f=\lambda
f\quad\leftrightarrow\quad\tilde{L}[(\sigma_{2}\otimes\sigma_{3})f]=-\lambda^{-1}(\sigma_{2}\otimes\sigma_{3})f,$
(3.5) $\displaystyle\tilde{L}:=\text{an operator such that $w$ in $\hat{L}$ is
replaced by $\tilde{w}=2\phi_{x}-w$.}$ (3.6)

We thus find the identity

$\displaystyle\hat{L}^{-1}=-(\sigma_{2}\otimes\sigma_{3})\tilde{L}(\sigma_{2}\otimes\sigma_{3}).$
(3.7)

It is unusual that the inverse of a matrix differential operator can be
written down explicitly, and furthermore, its expression does not include an
integral operator. Such situation seems to occur when the highest-order
coefficient matrix is not full-rank, which is $-4\mathrm{i}\sigma_{3}\oplus
O_{2}$ for the present $\hat{L}$.
In particular, for the stationary solution of the SG equation $\phi_{t}=0$, we
find $w=\tilde{w}=\phi_{x}$ and hence $\hat{L}=\tilde{L}$. In this case
$\hat{L}^{-1}=-(\sigma_{2}\otimes\sigma_{3})\hat{L}(\sigma_{2}\otimes\sigma_{3})$
holds and the eigenvalues $\lambda$ and $-\lambda^{-1}$ appear in pairs. We,
however, again emphasize that this relation does not hold for general
$w\neq\phi_{x}$, and therefore when we consider the scattering matrix $S$ in
the ISM in Sec. 5, we must not impose this symmetry to $\hat{L}$.

##  4 Soliton lattice solution, Riemann surface, and eigenfunctions

Let us consider the eigenfunction for the stationary solution, which we
henceforth write $\phi_{0}(x)$. Since the Lax pair commutes
$[\hat{L},\hat{B}]=0$ for the stationary solution, we consider the
simultaneous eigenfunction:

$\displaystyle\hat{L}f_{0}=\lambda f_{0},\quad\hat{B}f_{0}=\omega f_{0}.$
(4.1)

We first determine $\phi_{0}(x)$. The first integral for the stationary SG
equation is

$\displaystyle\phi_{0x}^{2}+2\cos\phi_{0}=\frac{4}{m}-2\quad\leftrightarrow\quad\frac{(\frac{1}{2}\phi_{0x})^{2}}{1-m\sin^{2}\frac{\phi_{0}-\pi}{2}}=\frac{1}{m}.$
(4.2)

Integrating this equation once again yields the soliton lattice solution as

$\displaystyle\phi_{0}(x)$
$\displaystyle=\pi+2\operatorname{am}\left(\tfrac{x}{\sqrt{m}}\big{|}m\right),$
(4.3) $\displaystyle\phi_{0}^{\prime}(x)$
$\displaystyle=\frac{2}{\sqrt{m}}\operatorname{dn}\left(\tfrac{x}{\sqrt{m}}\big{|}m\right),\quad\cos\tfrac{\phi_{0}(x)}{2}=-\operatorname{sn}\left(\tfrac{x}{\sqrt{m}}\big{|}m\right),\quad\sin\tfrac{\phi_{0}(x)}{2}=\operatorname{cn}\left(\tfrac{x}{\sqrt{m}}\big{|}m\right).$
(4.4)

Here and henceforth, we use the notations of elliptic functions in the
Abramowitz-Stegun book [22] unless otherwise noted. If $0<m<1$, it is the
rotating solution. If $m=1$, it represents the stationary one-kink solution,
and if $m>1$, it is an oscillating solution. In this paper we only consider
the rotating background $0<m<1$. Generally, the pair of commuting differential
operators satisfies an algebraic relation $P(\hat{L},\hat{B})=0$ [19].
Strictly speaking, Krichever’s original paper [19] only consider the case
where the highest-order coefficient matrices of the matrix differential
operators are invertible, but similar nature can be seen even if this
assumption is not satisfied and an algebraic curve can be defined in many
cases. In the present case, $P(\lambda,\omega)=0$ gives a genus-one curve,
i.e., an elliptic curve:

$\displaystyle\lambda^{2}\omega^{2}=\lambda^{4}+2\left(\frac{2}{m}-1\right)\lambda^{2}+1.$
(4.5)

It can be parametrized by elliptic functions:

$\displaystyle\lambda(z)$
$\displaystyle=\frac{-\mathrm{i}\sqrt{m}\operatorname{sn}(\mathrm{i}z|m)\operatorname{cn}(\mathrm{i}z|m)}{\operatorname{dn}(\mathrm{i}z|m)}=\mathrm{i}\sqrt{m}\operatorname{sn}(\mathrm{i}z|m)\operatorname{sn}(\mathrm{i}z-K|m),$
(4.6) $\displaystyle\omega(z)$
$\displaystyle=\frac{\operatorname{sn}^{2}(\mathrm{i}z|m)-\operatorname{sn}^{2}(\mathrm{i}z-K|m)}{\lambda(z)},$
(4.7)

where $K=K(m)$ and $K^{\prime}=K(1-m)$. The eigenfunction $f_{0}(x,z)$ for the
above parametrized eigenvalue $\lambda=\lambda(z)$ and the corresponding
crystal momentum $k(z)$ is

$\displaystyle k(z)$
$\displaystyle=\frac{Z(2\mathrm{i}z+\mathrm{i}K^{\prime}|m)+Z(2\mathrm{i}z-\mathrm{i}K^{\prime}|m)}{4\mathrm{i}\sqrt{m}},$
(4.8) $\displaystyle f_{0}(x,z)$
$\displaystyle=\frac{\mathrm{e}^{\mathrm{i}k(z)x}\Theta_{4}(0)}{2\Theta_{4}(\frac{x}{\sqrt{m}})}\begin{pmatrix}\Theta_{1}(\frac{x}{\sqrt{m}}-\mathrm{i}z)/\Theta_{4}(\mathrm{i}z)\\\
\Theta_{2}(\frac{x}{\sqrt{m}}-\mathrm{i}z)/\Theta_{3}(\mathrm{i}z)\\\
-\mathrm{i}\Theta_{4}(\frac{x}{\sqrt{m}}-\mathrm{i}z)/\Theta_{1}(\mathrm{i}z)\\\
-\mathrm{i}\Theta_{3}(\frac{x}{\sqrt{m}}-\mathrm{i}z)/\Theta_{2}(\mathrm{i}z)\end{pmatrix}.$
(4.9)

Here, we define the scaled theta functions by
$\Theta_{i}(u|m):=[\vartheta_{i}(\frac{\pi u}{2K},q)]_{\text{Abramowitz-
Stegun}}$ with the nome $q=\mathrm{e}^{-\pi K^{\prime}/K}$, and the second
argument $m$ is omitted. Using this scaled theta function, the Jacobi zeta
function is expressed as
$Z(u|m)=\frac{\mathrm{d}}{\mathrm{d}u}\ln\Theta_{4}(u|m)$. The derivation of
Eqs. (4.8) and (4.9) are given in B.
$\lambda(z),\omega(z),k(z)$ satisfy the following (quasi-)periodicity, parity,
and the complex conjugation relations:

$\displaystyle\lambda(z)$
$\displaystyle=(-1)^{l+n}\lambda(z+nK^{\prime}+\mathrm{i}lK)^{(-1)^{n}}=-\lambda(-z)=\lambda(z^{*})^{*},$
(4.10) $\displaystyle\omega(z)$
$\displaystyle=(-1)^{n}\omega(z+nK^{\prime}+\mathrm{i}lK)=-\omega(-z)=\omega(z^{*})^{*},$
(4.11) $\displaystyle k(z)$
$\displaystyle=k(z+nK^{\prime}+\mathrm{i}lK)+\frac{n\pi}{2K\sqrt{m}}=-k(-z)=k(z^{*})^{,}$
(4.12)

for $n,l\in\mathbb{Z}$. The eigenfunction $f_{0}(x,z)$ also has the symmetries

$\displaystyle f_{0}(x,z+nK^{\prime}+\mathrm{i}lK)$
$\displaystyle=\mathrm{i}^{l}\left[(\sigma_{3}^{l}\sigma_{2}^{n})\otimes(\sigma_{2}^{l}\sigma_{3}^{n})\right]f_{0}(x,z),$
(4.13) $\displaystyle f_{0}(x,-z^{*})$ $\displaystyle=(\sigma_{3}\otimes
I_{2})f_{0}(x,z)^{*}.$ (4.14)

The two linearly independent solution for a given eigenvalue
$\lambda=\lambda(z)$ are $f_{0}(x,z)$ and $f_{0}(x,-z-\mathrm{i}K)$, because
$\lambda(z)=\lambda(-z-\mathrm{i}K)$.
To cover all eigenfunctions, the region we must consider is the rectangle with
vertices $z=\pm K^{\prime}\pm\mathrm{i}K$, which is the fundamental period
parallelogram of $\lambda(z)$. The bounded eigenfunctions with
$k(z)\in\mathbb{R}$, i.e., the scattering eigenstates, lie on the lines
$z\in\mathbb{R}+\frac{\mathrm{i}K}{2}\mathbb{Z}$, and there are four
independent lines in one fundamental period parallelogram. Among them,
$\mathbb{R},\ \mathbb{R}-\mathrm{i}K$ gives monotonically increasing $k(z)$
and real $\lambda(z)$, while $\mathbb{R}\pm\frac{\mathrm{i}K}{2}$ gives
monotonically decreasing $k(z)$ and pure imaginary $\lambda(z)$. The region
satisfying $\operatorname{Im}k(z)>0$ is given by
$-K<\operatorname{Im}z<-\frac{K}{2},\ 0<\operatorname{Im}z<\frac{K}{2}$, and
the positions of zeros of $a(z)$ corresponding to bound eigenstates are chosen
from this region in Sec. 6.
The completeness relation satisfied for scattering eigenstates is

$\displaystyle\delta(x-y)I_{4}=\oint_{C_{1}+C_{2}}\frac{dz}{2\pi}f_{0}(x,z)f_{0}(y,z^{*})^{\dagger}\sigma.$
(4.15)

Here, $C_{1}$ is the rectangular contour passing through the vertices
$-K^{\prime}\to K^{\prime}\to
K^{\prime}+\frac{\mathrm{i}K}{2}\to-K^{\prime}+\frac{\mathrm{i}K}{2}\to-K^{\prime}$
and $C_{2}$ is the one $-K^{\prime}-\mathrm{i}K\to K^{\prime}-\mathrm{i}K\to
K^{\prime}-\frac{\mathrm{i}K}{2}\to-K^{\prime}-\frac{\mathrm{i}K}{2}\to-K^{\prime}-\mathrm{i}K$.
Note that the integrations along vertical lines cancel out due to the
periodicity. Therefore, the actual path with non-zero contribution is
$\oint_{C_{1}+C_{2}}=\int_{-K^{\prime}}^{K^{\prime}}-\int_{-K^{\prime}+\frac{\mathrm{i}K}{2}}^{K^{\prime}+\frac{\mathrm{i}K}{2}}+\int_{-K^{\prime}-\mathrm{i}K}^{K^{\prime}-\mathrm{i}K}-\int_{-K^{\prime}-\frac{\mathrm{i}K}{2}}^{K^{\prime}-\frac{\mathrm{i}K}{2}}$.
The derivation of Eq. (4.15) is given in C.

##  5 Scattering matrix

Henceforth, we consider the scattering problem under the condition that the
background potential tends to the stationary soliton lattice at spatial
infinities $x\to\pm\infty$. Let the potential $\phi(x)$ satisfy the boundary
condition

$\displaystyle\phi(x)\to\begin{cases}\phi_{0}(x)&(x\to-\infty),\\\
\phi_{0}(x-x_{0})&(x\to+\infty),\end{cases}$ (5.1)

where $x_{0}$ represents the phase shift of the background lattice caused by
solitons and ripple waves. We define the left and right Jost solutions by the
asymptotic form

$\displaystyle f_{-}(x,z)$ $\displaystyle\to f_{0}(x,z)\quad(x\to-\infty),$
(5.2) $\displaystyle f_{+}(x,z)$ $\displaystyle\to
f_{0}(x-x_{0},z)\quad(x\to+\infty).$ (5.3)

The scattering matrix $S(z)$ is then defined by the linear transformation
between these two:

$\displaystyle\begin{pmatrix}f_{+}(x,z)&f_{+}(x,-z-\mathrm{i}K)\end{pmatrix}=\begin{pmatrix}f_{-}(x,z)&f_{-}(x,-z-\mathrm{i}K)\end{pmatrix}S(z).$
(5.4)

The uniqueness of the definition of the Jost solution by its asymptotic form
and the properties of eigenfunctions (4.13) and (4.14) imply

$\displaystyle
f_{\pm}(x,z)=(-\mathrm{i}\sigma_{3}\otimes\sigma_{2})^{l}f_{\pm}(x,z+2nK^{\prime}+\mathrm{i}lK)=(\sigma_{3}\otimes
I_{2})f_{\pm}(x,-z^{*})^{*},\quad n,l\in\mathbb{Z}.$ (5.5)

Note that while $f_{0}(x,z)$ has an additional symmetry with respect to
translation of $z$ by $(2n+1)K^{\prime}$ (see Eq. (4.13)), the corresponding
symmetry does not exist for Jost solutions $f_{\pm}(x,z)$, because this
additional symmetry is coming from special nature of $\phi_{0}(x)$ satisfying
$\hat{L}=\tilde{L}$ in Eqs. (3.5)-(3.7). Combining Eqs. (5.4) and (5.5), we
have

$\displaystyle
S(z)=\sigma_{3}^{l}S(z+2nK^{\prime}+\mathrm{i}lK)\sigma_{3}^{l}=S(-z^{*})^{*}=\sigma_{1}S(-z-\mathrm{i}K)\sigma_{1},\quad
n,l\in\mathbb{Z}.$ (5.6)

Equating the integration constant obtained by integrating Eq. (3.2) at
$x\to\pm\infty$, we find

$\displaystyle S(z)^{-1}=S(z^{*})^{\dagger}.$ (5.7)

We can also show $\det S=1$ from the constancy of the Wronskian
$\det(g_{1},g_{2})$ (see the sentence after Eq. (3.2)). Therefore, $S$ and its
inverse matrix has the forms

$\displaystyle S(z)$ $\displaystyle=\begin{pmatrix}a(z)&-b(z^{*})^{*}\\\
b(z)&a(z^{*})^{*}\end{pmatrix},$ (5.8) $\displaystyle S(z)^{-1}$
$\displaystyle=\begin{pmatrix}a(z^{*})^{*}&b(z^{*})^{*}\\\
-b(z)&a(z)\end{pmatrix},$ (5.9)

with matrix elements satisfying

$\displaystyle a(z)a(z^{*})^{*}+b(z)b(z^{*})^{*}=1,$ (5.10) $\displaystyle
a(z)=a(-z^{*})^{*}=a(z-\mathrm{i}K),\quad
b(z)=b(-z^{*})^{*}=-b(z-\mathrm{i}K).$ (5.11)

From this, the zeros of $a(z)$ always appear simultaneously at four points
$z,-z^{*},z-\mathrm{i}K,-z^{*}-\mathrm{i}K$, and the orders of these zeros are
the same.

##  6 Integral representation of Jost solution and Gelfand-Levitan-Marchenko
equation

Let us introduce the integral representation for the left Jost solution

$\displaystyle
f_{-}(x,z)=f_{0}(x,z)+\int_{-\infty}^{x}dy\Gamma(x,y)f_{0}(y,z).$ (6.1)

The kernel function $\Gamma(x,y)$ is a $4\times 4$ matrix and assumed to
decrease exponentially in the limits $x,y\to-\infty$. In order for this
integral to converge, the integrand must decrease $y\to-\infty$, and this
condition is satisfied at least for eigenfunctions with
$\operatorname{Im}k\leq 0$, i.e., $-K/2\leq\operatorname{Im}z\leq 0,\
K/2\leq\operatorname{Im}z\leq K$. From Eq. (3.4), we immediately conclude that

$\displaystyle\Gamma(x,y)=(\sigma_{3}\otimes\sigma_{2})\Gamma(x,y)(\sigma_{3}\otimes\sigma_{2})=(\sigma_{3}\otimes
I_{2})\Gamma(x,y)^{*}(\sigma_{3}\otimes
I_{2})=(I_{2}\otimes\sigma_{2})\Gamma(x,y)^{*}(I_{2}\otimes\sigma_{2}).$ (6.2)

Let us now write

$\displaystyle\hat{L}_{0}(x)$ $\displaystyle=-4\mathrm{i}\sigma_{3}\oplus
O_{2}+U_{0}(x),\quad
U_{0}(x)=\begin{pmatrix}\mathrm{i}\sigma_{1}w_{0}&\mathrm{e}^{-\mathrm{i}\phi_{0}\sigma_{2}/2}\\\
\mathrm{e}^{-\mathrm{i}\phi_{0}\sigma_{2}/2}&O_{2}\end{pmatrix},$ (6.3)
$\displaystyle\hat{L}(x)$ $\displaystyle=-4\mathrm{i}\sigma_{3}\oplus
O_{2}+U(x),\quad
U_{0}(x)=\begin{pmatrix}\mathrm{i}\sigma_{1}w&\mathrm{e}^{-\mathrm{i}\phi\sigma_{2}/2}\\\
\mathrm{e}^{-\mathrm{i}\phi\sigma_{2}/2}&O_{2}\end{pmatrix}.$ (6.4)

Following the same derivation as Ref. [23], section 2.8, (In the present case,
the first-order coefficient matrix of $\hat{L}$ is not full-rank, but it makes
no difference to the proof.)

$\displaystyle U(x)-U_{0}(x)=4\mathrm{i}\big{[}\sigma_{3}\oplus
O_{2},\Gamma(x,x)\big{]},$ (6.5)
$\displaystyle\hat{L}(x)\Gamma(x,y)=\Gamma(x,y)\hat{L}_{0}(y).$ (6.6)

From Eq. (6.5),

$\displaystyle w$ $\displaystyle=w_{0}+2\Gamma_{12},$ (6.7)
$\displaystyle\mathrm{e}^{\mathrm{i}\phi/2}$
$\displaystyle=\mathrm{e}^{\mathrm{i}\phi_{0}/2}+4(\mathrm{i}\Gamma_{13}+\Gamma_{14}).\quad\left(\leftrightarrow\cos\tfrac{\phi}{2}=\cos\tfrac{\phi_{0}}{2}+4\mathrm{i}\Gamma_{13},\
\sin\tfrac{\phi}{2}=\sin\tfrac{\phi_{0}}{2}-4\mathrm{i}\Gamma_{14}.\right),$
(6.8)

where we briefly write the matrix element of $\Gamma(x,x)$ as
$\Gamma_{ij}=[\Gamma(x,x)]_{ij}$. Note that $\Gamma_{12},\
\mathrm{i}\Gamma_{13}$, and $\mathrm{i}\Gamma_{14}$ are real-valued functions
due to Eq. (6.2). Using these relations, we can re-construct the potential
from the kernel $\Gamma$, which is determined from scattering data by solving
the GLM equation shown below.
The integral kernel $\Gamma(x,y)$ is uniquely determined from the following
scattering data: (i) the values of the reflection coefficient for scattering
states $r(z)=b(z)/a(z)$, where scattering states mean the bounded
eigenfunction s.t. $k(z)\in\mathbb{R}$, appearing on lines
$z\in\mathbb{R}+\frac{\mathrm{i}K}{2}\mathbb{Z}$, and (ii) the list of the
zeros of $a(z)$, which we write $z_{1},\dots,z_{n}$, and the normalization
factor $c_{j}^{2}:=\frac{b(z_{j})}{\mathrm{i}\dot{a}(z_{j})}$, where $z_{j}$’s
are to be chosen from the region $\operatorname{Im}k(z_{j})>0$, and due to the
property (5.11), the zeros simultaneously appear at
$z,-z^{*},z-\mathrm{i}K,-\mathrm{z}^{*}-\mathrm{i}K$. The values of
$c_{j}^{2}$ also have some constraint, whose detail will be described in Sec.
9. Here, we only treat the case where all zeros of $a(z)$ are first order. For
uniform background, the solitons corresponding to higher-order zeros are
discussed in Ref. [24]. The case of unstable uniform background is considered
in Ref. [25].
The GLM equation that determines the kernel $\Gamma(x,y)$ from the above-
mentioned scattering data is given by

$\displaystyle\Gamma(x,y)+\Omega(x,y)+\int_{-\infty}^{x}\mathrm{d}w\Gamma(x,w)\Omega(w,y)=0,\quad
y\leq x.$ (6.9)

Here, $\Omega=\Omega_{\text{bd}}+\Omega_{\text{sc}}$ with

$\displaystyle\Omega_{\text{bd}}(x,y)$
$\displaystyle:=W(x)W(y)^{T}(I_{2}\otimes\sigma_{1}),$ (6.10) $\displaystyle
W(x)$
$\displaystyle:=-\mathrm{i}\sigma_{3}\otimes\sigma_{2}\left(f_{0}(x,-z_{1})c_{1},\dots,f_{0}(x,-z_{n})c_{n}\right),$
(6.11) $\displaystyle\Omega_{\text{sc}}(x,y)$
$\displaystyle:=\left[\int_{-K^{\prime}}^{K^{\prime}}-\int_{-K^{\prime}+\frac{\mathrm{i}K}{2}}^{K^{\prime}+\frac{\mathrm{i}K}{2}}+\int_{-K^{\prime}-\mathrm{i}K}^{K^{\prime}+\mathrm{i}K}-\int_{-K^{\prime}-\frac{\mathrm{i}K}{2}}^{K^{\prime}-\frac{\mathrm{i}K}{2}}\right]\frac{\mathrm{d}z}{2\pi}f_{0}(x,-z-\mathrm{i}K)r(z)f_{0}(y,z^{*})^{\dagger}\sigma.$
(6.12)

$\Omega_{\text{bd}}$ represents the contribution from bound states, and
$\Omega_{\text{sc}}$ is that from scattering states. The derivation of the GLM
equation (6.9) with (6.10)-(6.12) is given in D.

##  7 Reflectionless potentials

Let us solve the GLM equation (6.9) for the reflectionless case, i.e.,
$\Omega_{\text{sc}}=0$. The ansatz for the kernel $\Gamma(x,y)$ is as follows.
Let us introduce the notation $H(x)=(h_{1}(x),\dots,h_{n}(x))$, where each
$h_{i}(x)$ is four-component column vector. We then assume

$\displaystyle\Gamma(x,y)=H(x)W(y)^{T}(I_{2}\otimes\sigma_{1}).$ (7.1)

Substituting it to the GLM equation, we obtain the equation for $H(x)$:

$\displaystyle H(x)+W(x)+H(x)G(x)=0,$ (7.2)

where we define the $n\times n$ Gram matrix by

$\displaystyle G(x)=\int_{-\infty}^{x}dyW(y)^{T}(I_{2}\otimes\sigma_{1})W(y).$
(7.3)

From Eq. (7.2), $H=-W(I_{n}+G)^{-1}$. Let $G_{lj}(x)$ be the $(l,j)$-component
of $G(x)$. Using Eq. (3.2) and the symmetries of $\lambda(z)$ and
$f_{0}(x,z)$, Eqs. (4.10), (4.13), and (4.14), we find

$\displaystyle G_{lj}(x)=c_{l}c_{j}\frac{4f_{0}(x,-z_{l})^{T}(\sigma_{2}\oplus
O_{2})f_{0}(x,-z_{j})}{\lambda(z_{l})-\lambda(z_{j})}.$ (7.4)

It can be simplified using the Weierstrass three-term formula in Ref. [20],
Eq. (3.8). The result is

$\displaystyle
G_{lj}(x)=c_{l}c_{j}\mathrm{e}^{-\mathrm{i}[k(z_{l})+k(z_{j})]x}\frac{\Theta_{2}\Theta_{4}\Theta_{3}(\frac{x}{\sqrt{m}}+\mathrm{i}(z_{l}+z_{j}))}{\Theta_{3}\Theta_{2}(\mathrm{i}(z_{l}+z_{j}))\Theta_{4}(\frac{x}{\sqrt{m}})}.$
(7.5)

Thus,

$\displaystyle\Gamma(x,y)=-W(x)[I_{n}+G(x)]^{-1}W(y)^{T}(I_{2}\otimes\sigma_{1}).$
(7.6)

Let us write $W$ as an array of row vectors:
$W=\left(\begin{smallmatrix}W_{1}\\\ W_{2}\\\ W_{3}\\\
W_{4}\end{smallmatrix}\right)$, where each $W_{i}$ is a $1\times n$ row
vector. Then we get the expressions
$\Gamma_{13}=-W_{1}(I_{n}+G)^{-1}W_{4}^{T},\
\Gamma_{14}=-W_{1}(I_{n}+G)^{-1}W_{3}^{T}$. Using the relation (6.8) and the
Woodbury-type identity
$1+\vec{a}^{\dagger}A^{-1}\vec{b}=\frac{\det(A+\vec{b}\vec{a}^{\dagger})}{\det
A}$,

$\displaystyle\cos\tfrac{\phi}{2}$
$\displaystyle=\cos\tfrac{\phi_{0}}{2}+4\mathrm{i}\Gamma_{13}=\left(\cos\tfrac{\phi_{0}}{2}\right)\frac{\det(I+P)}{\det(I+G)},\quad
P_{ij}=G_{ij}-\frac{4\mathrm{i}}{\cos\frac{\phi_{0}}{2}}f_{03}(x,-z_{i})f_{02}(x,-z_{j}),$
(7.7) $\displaystyle\sin\tfrac{\phi}{2}$
$\displaystyle=\sin\tfrac{\phi_{0}}{2}-4\mathrm{i}\Gamma_{14}=\left(\sin\tfrac{\phi_{0}}{2}\right)\frac{\det(I+Q)}{\det(I+G)},\quad
Q_{ij}=G_{ij}-\frac{4\mathrm{i}}{\sin\frac{\phi_{0}}{2}}f_{04}(x,-z_{i})f_{02}(x,-z_{j}),$
(7.8)

where $f_{0j}(x,z)$ represents the $j$-th component of $f_{0}(x,z)$ defined by
(4.9). These can be simplified by using Ref. [20], Eq. (3.5b). Rewriting
$C_{j}=\frac{1}{c_{j}^{2}\sqrt{m}}$, and introducing new matrices
$\mathcal{G},\mathcal{P},\mathcal{Q}$ through the relations
$G_{ij}=m^{1/4}c_{i}c_{j}\mathrm{e}^{-\mathrm{i}[k(z_{i})+k(z_{j})]x}\mathcal{G}_{ij},\
P_{ij}=m^{1/4}c_{i}c_{j}\mathrm{e}^{-\mathrm{i}[k(z_{i})+k(z_{j})]x}\mathcal{P}_{ij},\
Q_{ij}=m^{1/4}c_{i}c_{j}\mathrm{e}^{-\mathrm{i}[k(z_{i})+k(z_{j})]x}\mathcal{Q}_{ij}$,
we obtain the final expressions

$\displaystyle\cos\tfrac{\phi}{2}$
$\displaystyle=\left(\cos\tfrac{\phi_{0}}{2}\right)\frac{\det(\mathcal{E}+\mathcal{P})}{\det(\mathcal{E}+\mathcal{G})},$
(7.9) $\displaystyle\sin\tfrac{\phi}{2}$
$\displaystyle=\left(\sin\tfrac{\phi_{0}}{2}\right)\frac{\det(\mathcal{E}+\mathcal{Q})}{\det(\mathcal{E}+\mathcal{G})},$
(7.10)

where the components of $n\times n$ matrices $\mathcal{E},\ \mathcal{G},\
\mathcal{P},\ $ and $\mathcal{Q}$ are given by

$\displaystyle\mathcal{E}_{ij}$
$\displaystyle=\delta_{ij}C_{j}\mathrm{e}^{2\mathrm{i}k(z_{j})x},$ (7.11)
$\displaystyle\mathcal{G}_{ij}$
$\displaystyle=\frac{\Theta_{4}\Theta_{3}(\frac{x}{\sqrt{m}}+\mathrm{i}(z_{i}+z_{j}))}{\Theta_{2}(\mathrm{i}(z_{i}+z_{j}))\Theta_{4}(\frac{x}{\sqrt{m}})},$
(7.12) $\displaystyle\mathcal{P}_{ij}$
$\displaystyle=\left(\frac{\Theta_{4}(\mathrm{i}z_{i})}{\Theta_{1}(\mathrm{i}z_{i})}\right)\left(\frac{-\Theta_{4}\Theta_{2}(\frac{x}{\sqrt{m}}+\mathrm{i}(z_{i}+z_{j}))}{\Theta_{2}(\mathrm{i}(z_{i}+z_{j}))\Theta_{1}(\frac{x}{\sqrt{m}})}\right)\left(\frac{\Theta_{2}(\mathrm{i}z_{j})}{\Theta_{3}(\mathrm{i}z_{j})}\right),$
(7.13) $\displaystyle\mathcal{Q}_{ij}$
$\displaystyle=\left(\frac{\Theta_{4}(\mathrm{i}z_{i})}{\Theta_{2}(\mathrm{i}z_{i})}\right)\left(\frac{\Theta_{4}\Theta_{1}(\frac{x}{\sqrt{m}}+\mathrm{i}(z_{i}+z_{j}))}{\Theta_{2}(\mathrm{i}(z_{i}+z_{j}))\Theta_{2}(\frac{x}{\sqrt{m}})}\right)\left(\frac{\Theta_{1}(\mathrm{i}z_{j})}{\Theta_{3}(\mathrm{i}z_{j})}\right).$
(7.14)

Since the wavenumbers used to make bound states all satisfy
$\operatorname{Im}k(z_{j})>0$,

$\displaystyle\mathcal{E}\to\begin{cases}\infty&(x\to-\infty),\\\
0&(x\to+\infty).\end{cases}$ (7.15)

Therefore, the asymptotic form of $\phi$ is given by

$\displaystyle\cos\tfrac{\phi}{2}$
$\displaystyle\to\begin{cases}\cos\tfrac{\phi_{0}}{2}&(x\to-\infty),\\\
\cos\tfrac{\phi_{0}}{2}\det\mathcal{P}\mathcal{G}^{-1}&(x\to+\infty),\end{cases}$
(7.16) $\displaystyle\sin\tfrac{\phi}{2}$
$\displaystyle\to\begin{cases}\sin\tfrac{\phi_{0}}{2}&(x\to-\infty),\\\
\sin\tfrac{\phi_{0}}{2}\det\mathcal{Q}\mathcal{G}^{-1}&(x\to+\infty).\end{cases}$
(7.17)

The determinant appearing in $x\to+\infty$ is calculated in E. The resultant
asymptotic form is

$\displaystyle\phi(x)\to\begin{cases}\phi(x)&(x\to-\infty),\\\
\phi(x-x_{0})&(x\to+\infty),\end{cases}$ (7.18) $\displaystyle
x_{0}:=-2\mathrm{i}\sqrt{m}\sum_{j=1}^{n}z_{n}.$ (7.19)

$x_{0}$ has the meaning of the phase shift of the background lattice caused by
solitons. The realness of $x_{0}$ follows from the fact that the zeros of
$a(z)$ always simultaneously emerge at four points $z,z-\mathrm{i}K,-z^{*},$
and $-z^{*}-\mathrm{i}K$.

##  8 Time evolution

Finally, let us determine the time dependence of scattering matrix if the
system obeys the Lax equation (2.1). Let $f$ be a time-dependent eigenfunction
of $\hat{L}$. Then, the Lax equation implies

$\displaystyle\hat{L}f$ $\displaystyle=\lambda f,$ (8.1) $\displaystyle
4\mathrm{i}f_{t}$ $\displaystyle=-\hat{B}f.$ (8.2)

We now define the time-dependent eigenfunction by the initial condition
$\tilde{f}_{+}(t=0,x,z)=f_{+}(x,z)$. Though $\hat{B}$ is now a time-dependent
operator for finite $x$, it is time-independent at $x\to\pm\infty$. Therefore,
the asymptotic forms of this $\tilde{f}_{+}(t,x,z)$ is given by

$\displaystyle(\tilde{f}_{+}(t,x,z),\tilde{f}_{+}(t,x,-z-\mathrm{i}K))\to\begin{cases}(f_{0}(x,z),f_{0}(x,-z-\mathrm{i}K))\mathrm{e}^{\mathrm{i}\omega(z)t\sigma_{3}/4}S(z)&(x\to-\infty),\\\
(f_{0}(x-x_{0},z),f_{0}(x-x_{0},-z-\mathrm{i}K))\mathrm{e}^{\mathrm{i}\omega(z)t\sigma_{3}/4}&(x\to+\infty).\end{cases}$
(8.3)

Here, we have used the relation $\omega(-z-\mathrm{i}K)=-\omega(z)$.
Therefore, if we define the time-dependent right Jost solution by the
asymptotic form $f_{+}(t,x,z)\to f_{0}(x-x_{0},z)\text{ for }x\to+\infty$,
which is different from $\tilde{f}_{+}(t,x,z)$, then the relation
$(\tilde{f}_{+}(t,x,z),\tilde{f}_{+}(t,x,-z-\mathrm{i}K))=(f_{+}(t,x,z),f_{+}(t,x,-z-\mathrm{i}K))\mathrm{e}^{\mathrm{i}\omega(z)t/4}$
holds. Then the time evolution of the scattering matrix is given by

$\displaystyle
S(t,z)=\mathrm{e}^{\mathrm{i}\omega(z)t\sigma_{3}/4}S(z)\mathrm{e}^{-\mathrm{i}\omega(z)t\sigma_{3}/4},$
(8.4)

or for each component,

$\displaystyle a(t,z)=a(z),\quad
b(t,z)=\mathrm{e}^{-\mathrm{i}\omega(z)t/2}b(z).$ (8.5)

Since $C_{j}$ is defined by
$C_{j}=\frac{1}{c_{j}^{2}\sqrt{m}}=\frac{\mathrm{i}\dot{a}(z_{j})}{\sqrt{m}b(z_{j})}$,
its time evolution is written as

$\displaystyle C_{j}(t)=C_{j}\mathrm{e}^{\mathrm{i}\omega(z_{j})t/2}.$ (8.6)

Therefore, we can define the time-dependent $\mathcal{E}_{ij}(t)$ as follows.
In Eq. (7.11), we replace $C_{j}$ by
$C_{j}(t)=C_{j}\mathrm{e}^{\mathrm{i}\omega(z_{j})t/2}$, i.e.,

$\displaystyle\mathcal{E}_{ij}(t)=\delta_{ij}C_{j}\mathrm{e}^{2\mathrm{i}k(z_{j})x+\mathrm{i}\omega(z_{j})t/2}.$
(8.7)

Then the time-dependent solution of the SG equation is given by replacing
$\mathcal{E}$ with $\mathcal{E}(t)$ in Eqs. (7.9) and (7.10):

$\displaystyle\cos\tfrac{\phi}{2}$
$\displaystyle=\left(\cos\tfrac{\phi_{0}}{2}\right)\frac{\det(\mathcal{E}(t)+\mathcal{P})}{\det(\mathcal{E}(t)+\mathcal{G})},$
(8.8) $\displaystyle\sin\tfrac{\phi}{2}$
$\displaystyle=\left(\sin\tfrac{\phi_{0}}{2}\right)\frac{\det(\mathcal{E}(t)+\mathcal{Q})}{\det(\mathcal{E}(t)+\mathcal{G})},$
(8.9)

where the definition of $\mathcal{G},\mathcal{P}$, and $\mathcal{Q}$ remains
unchanged (Eqs. (7.12)-(7.14)).

##  9 Limitation to eigenvalues $z_{j}$’s and normalization coefficients
$C_{j}$’s

In order for the function $\phi(x,t)$ to be real and bounded, $z_{j}$’s, which
are the zeros of $a(z)$ corresponding to the discrete eigenvalues for bound
states, and the normalization coefficients $C_{j}$’s must satisfy several
conditions.111If we consider application to scientific problems where the
complex-valued and/or divergent solutions have some appropriate physical
interpretations, the restriction stated here can be loosened. Here we describe
it.
First, depending on its value, $z_{j}$’s must satisfy the following (i) or
(ii):

  1. (i)

$z_{j},\ z_{j}-\mathrm{i}K,\ -z_{j}^{*},\ -z_{j}^{*}-\mathrm{i}K$ appear
simultaneously. Due to this symmetry, one of these four $z_{j}$ can be chosen
from $0\leq\operatorname{Re}z_{j}\leq K^{\prime},\
0<\operatorname{Im}z_{j}<\frac{K}{2}$ without loss of generality.

  2. (ii)

If, as a special case, $z_{j}=\mathrm{i}r_{j}$ or
$z_{j}=\mathrm{i}r_{j}+K^{\prime}$ with $r_{j}\in(0,\frac{K}{2}),$ then just
two zeros $z_{j}$ and $z_{j}-\mathrm{i}K$ appear simultaneously.

The case (i) corresponds to the breather solution and the case (ii) the
traveling one kink solution. The total number of $z_{j}$’s are always even. If
there are breathers and even number of kinks, the number is a multiple of
four, while if there exists odd number of kinks, the number is of the form
$4n+2$. The velocity of the soliton is given by
$v_{j}=-\frac{\operatorname{Im}\omega(z_{j})}{4\operatorname{Im}k(z_{j})}$. In
particular, for one kink solution, it reduces to
$v_{j}=-\frac{\omega(z_{j})}{4k(z_{j})}$, which is negative if
$z_{j}=\mathrm{i}r_{j}$ and positive if $z_{j}=\mathrm{i}r_{j}+K^{\prime}$.
Next let us consider the condition for $C_{j}$’s. By Eq. (5.11), the
derivative of $a(z)$ satisfies

$\displaystyle\dot{a}(z)=-\dot{a}(-z^{*})^{*}=\dot{a}(z-\mathrm{i}K),$ (9.1)

and thus if we define

$\displaystyle C(z):=\frac{\mathrm{i}\dot{a}(z)}{\sqrt{m}b(z)},$ (9.2)

then it has the symmetry

$\displaystyle C(z)=C(-z^{*})^{*}=-C(z-\mathrm{i}K).$ (9.3)

This symmetry is the same as that of $b(z)$ given in Eq. (5.11). Therefore,
for the breather solution where four zeros appear, if we temporarily write
these four as $z_{j},z_{j+1}=z_{j}-\mathrm{i}K,\ z_{j+2}=-z_{j}^{*},\
z_{j+3}=-z_{j}^{*}-\mathrm{i}K$, then we should choose the corresponding
coefficients as $C_{j+1}=-C_{j},\ C_{j+2}=C_{j}^{*},\ C_{j+3}=-C_{j}^{*}$. For
the kink solution, where two zeros $z_{j}$ and $z_{j+1}=z_{j}-\mathrm{i}K$
appear, $C_{j+1}=-C_{j}$ are both real and must have the opposite sign.
Whether the solution becomes kink or anti-kink, i.e., the phase rotation
becomes counterclockwise or clockwise, is determined by which one is chosen
positive.

##  10 Animation of soliton solutions

Here, we present a few examples of the soliton solutions, Eqs. (8.8) and
(8.9), by gif animations. We visualize the solution through the 3D plot
$(x,\,\cos\phi(x,t),\,\sin\phi(x,t))$. An example of snapshot is shown in Fig.
1. See attached files, testx.gif with x=1,2, and 3. Below, we provide the
parameters for each solution. We write $K=K(m)$ and $K^{\prime}=K(1-m)$.

  1. (i)

test1.gif: one kink solution.
The number of zeros is $n=2$. The parameters are $m=0.99,\
z_{1}=0.02\mathrm{i}K+K^{\prime},\ z_{2}=z_{1}-\mathrm{i}K,\ C_{1}=-1,\
C_{2}=-C_{1}$. The time range is $-10\leq t\leq 10$.

  2. (ii)

test2.gif: one anti-kink solution.
$C_{1}=1,\ C_{2}=-C_{1}$, and other parameters are the same as (i).

  3. (iii)

test3.gif: one breather solution.
The number of zeros is $n=4$. The parameters are $m=0.99,\
z_{1}=0.1\mathrm{i}K+0.5K^{\prime},\ z_{2}=z_{1}-\mathrm{i}K,\
z_{3}=-z_{1}^{*},\ z_{4}=-z_{2}^{*},\ C_{1}=-1,\ C_{2}=-C_{1},\
C_{3}=C_{1}^{*},\ C_{4}=C_{2}^{*}$. The time range is $-10\leq t\leq 10$.

Solutions with more number of solitons will be generated using the Mathematica
file in google drive.222In this file, since the default build-in Jacobi
elliptic functions are a little slow, the functions defined by the ratio of
theta functions are used alternatively.

Figure 1:  A snapshot of one kink solution. The number of zeros is $n=2$ and
the parameters are $m=0.95,\ z_{1}=0.02\mathrm{i}K+K^{\prime},\
z_{2}=z_{1}-\mathrm{i}K,\ C_{1}=-1,\ C_{2}=-C_{1}$. The snapshot at $t=1$.

##  11 Summary and discussion

In this paper, we have derived the multi-soliton solutions of the SG equation
with elliptic-function background by the ISM. The result is expressed by a
determinant of theta functions, i.e., Eqs. (8.8) and (8.9). The shift of the
background lattice caused by solitons (7.18) with (7.19) has also been found
using the addition formula of theta functions. One key tool in our work is the
Lax pair represented by $4\times 4$ matrix differential operators (2.2) and
(2.3), originally introduced in Ref. [17]. This most conventional but
seemingly outdated approach makes it possible to use the common form of the
integral representation of the Jost solution (6.1) without modification and
the formulation of the ISM is simplified. Also, the completeness relations in
an indefinite inner product space, which is necessary in the formulation of
the ISM, has been discussed in detail (C). Application to various physical
phenomena including the dynamics of defects with periodic background,
extension to unstable oscillating background ($m>1$) and studying multi rogue
waves, and the solitons associated with higher-order zeros of $a(z)$, are left
as possible future works.
Before concluding, we provide several technical remarks and perspectives. In
this work, we can introduce the concept of orthogonality between
eigenfunctions based on the indefinite inner product, since the Lax operators
have the symmetry (3.1). If the Lax pair has no such symmetry and any inner
product cannot be defined, the completeness relation will be constructed from
the set of left and right eigenfunctions, an analog of left and right
eigenvectors for finite dimensional matrix. Such basis is sometimes called a
bi-orthogonal basis.
Whether we can always reduce the classical integrable systems written by zero-
curvature condition (or a compatibility condition) and/or the Lax pair
including an integral operator to the Lax formalism with differential
operators seems to be unclear. (The inverse operation is easy; if one has an
integrable equation written by the Lax form using matrix differential
operators, one can immediately rewrite it in a zero-curvature expression.) If
we can do it, the integral representation of the Jost solution will be widely
applicable and multi-soliton solutions will be easily obtained by the dressing
method [26], where we can even omit the full formulation of the ISM if we do
not have an interest in general initial-value problem. We also note that, as
emphasized in Secs. 3 and 4, the highest-order coefficient matrix of the Lax
operator $\hat{L}$ is not full-rank and $\hat{L}^{-1}$ can be written down
explicitly without using an integral operator. Though we do not show detail
here, a similar property also emerges in the derivative nonlinear Schrödinger
equation, which is usually treated by the Kaup-Newell form. That is, if we
rewrite it to the Lax form, the highest-order coefficient matrix is not full-
rank. Even in these systems, the algebraic curve from the pair of commuting
operators can be appropriately defined (see Eq. (4.5)), though, strictly
speaking, these cases are not included in the seminal work [19]. Exploring
these examples more comprehensively and exhaustively might bring a slight
extension to the theory of commuting differential operators.

This work was supported by MEXT-Supported Program Grant No. S1511006 and JSPS
KAKENHI Grant No. JP19H05821.

##  Appendix A Integral often emerging in calculation of eigenfunctions for
elliptic potentials

Let us define
$\operatorname{scd}(z|m):=\operatorname{sn}(z|m)\operatorname{cn}(z|m)\operatorname{dn}(z|m)$.
When an elliptic-function potential is given in some hierarchy of integrable
systems, the calculation of its eigenfunctions often reduces to the following
integral:

$\displaystyle\int dx\frac{\alpha[\operatorname{scd}(\alpha
x)+\operatorname{scd}\beta]}{\operatorname{sn}^{2}(\alpha
x)-\operatorname{sn}^{2}\beta}=\alpha Z(\beta)x+\ln\frac{\Theta_{1}(\alpha
x-\beta)}{\Theta_{4}(\beta)\Theta_{4}(\alpha x)},$ (A.1)

even if the corresponding Riemann surface has higher genus $g>1$. This formula
can be derived by using [27], appendix B. If we interpret this formula using
Weierstrass’s functions, it reduces to the formula given in Ref. [28], II,
Kap. 6, §3, and is called the standard form of the elliptic integral of the
third kind. Using it, the solution to the first-order differential equation

$\displaystyle\frac{f_{x}}{f}=\gamma+\frac{\alpha[\operatorname{scd}\alpha
x+\operatorname{scd}\beta]}{\operatorname{sn}^{2}\alpha
x-\operatorname{sn}^{2}\beta}$ (A.2)

is given by

$\displaystyle f(x)=C\mathrm{e}^{\mathrm{i}kx}\frac{\Theta_{1}(\alpha
x-\beta)}{\Theta_{4}(\beta)\Theta_{4}(\alpha x)},\ k=-\mathrm{i}(\gamma+\alpha
Z(\beta)).$ (A.3)

If $m\in[0,1]$, $\alpha>0$, and $k$ is real, then $f$ becomes a twisted
periodic function

$\displaystyle
f\left(x+\tfrac{4Kl}{\alpha}\right)=\mathrm{e}^{\mathrm{i}k\frac{4Kl}{\alpha}}f(x),\quad
l\in\mathbb{Z}.$ (A.4)

Therefore, if $k$ can be physically interpreted to be crystal momentum of the
Bloch function, it is defined up to $\operatorname{mod}\frac{\alpha\pi}{2K}$
and the corresponding Brillouin zone is given by
$[-\frac{\alpha\pi}{4K},\frac{\alpha\pi}{4K}]$.

##  Appendix B Eigenfunctions for soliton lattice

In this appendix, we derive the simultaneous eigenfunction $f_{0}(x,z)$ with
its crystal momentum $k(z)$ in Eqs. (4.8) and (4.9) for the time-independent
Lax pair $\hat{L}$ and $\hat{B}$ with the soliton lattice potential (4.3),
when the eigenvalues are parametrized by Eqs. (4.6) and (4.7). The method is
in fact the special case of the one given by Krichever [19].
Let us write $f_{0}=\left(\begin{smallmatrix}g\\\
\frac{1}{\lambda}\mathrm{e}^{-\mathrm{i}\phi\sigma_{2}}g\end{smallmatrix}\right),$
where $g$ is a two-component vector $g=\left(\begin{smallmatrix}g_{1}\\\
g_{2}\end{smallmatrix}\right)$. Eliminating $g_{1}^{\prime}$ and
$g_{2}^{\prime}$ using the eigenequation of $\hat{L}$ in Eq. (4.1), and
substituting them to that of $\hat{B}$, one obtains two linear relations with
respect to $g_{1}$ and $g_{2}$. The determinant of this coefficient matrix
must vanish, which gives an elliptic curve (4.5). Using it, one can eliminate
either of $g_{1}$ and $g_{2}$ and obtain the first-order differential equation
only containing one function:

$\displaystyle\frac{g_{1}^{\prime}}{g_{1}}$
$\displaystyle=\frac{-\mathrm{i}(2\lambda-\omega)}{4}+\frac{\mathrm{i}(-1+\lambda^{4}-2\lambda^{3}\omega+\lambda^{2}\omega^{2})}{2\lambda(-1+\lambda^{2}-\lambda\omega+2\cos^{2}\frac{\phi_{0}}{2})}+\frac{-\phi_{0x}\cos\frac{\phi_{0}}{2}\sin\frac{\phi_{0}}{2}}{-1+\lambda^{2}-\lambda\omega+2\cos^{2}\frac{\phi_{0}}{2}},$
(B.1) $\displaystyle\frac{g_{2}^{\prime}}{g_{2}}$
$\displaystyle=\frac{\mathrm{i}(2\lambda+\omega)}{4}+\frac{-\mathrm{i}(-1+\lambda^{4}+2\lambda^{3}\omega+\lambda^{2}\omega^{2})}{2\lambda(-1+\lambda^{2}+\lambda\omega+2\cos^{2}\frac{\phi_{0}}{2})}+\frac{-\phi_{0x}\cos\frac{\phi_{0}}{2}\sin\frac{\phi_{0}}{2}}{-1+\lambda^{2}+\lambda\omega+2\cos^{2}\frac{\phi_{0}}{2}}.$
(B.2)

Using Eq. (4.4) and the parametrizations (4.6) and (4.7), the above equations
reduce to

$\displaystyle\frac{g_{1}^{\prime}}{g_{1}}$
$\displaystyle=\frac{-\mathrm{i}(2\lambda-\omega)}{4}+\frac{1}{\sqrt{m}}\frac{\operatorname{scd}(\frac{x}{\sqrt{m}})+\operatorname{scd}(\mathrm{i}z)}{\operatorname{sn}^{2}(\frac{x}{\sqrt{m}})-\operatorname{sn}^{2}(\mathrm{i}z)},$
(B.3) $\displaystyle\frac{g_{2}^{\prime}}{g_{2}}$
$\displaystyle=\frac{\mathrm{i}(2\lambda+\omega)}{4}+\frac{1}{\sqrt{m}}\frac{\operatorname{scd}(\frac{x}{\sqrt{m}})+\operatorname{scd}(\mathrm{i}z-K)}{\operatorname{sn}^{2}(\frac{x}{\sqrt{m}})-\operatorname{sn}^{2}(\mathrm{i}z-K)}.$
(B.4)

Therefore, using the formula in A, we find

$\displaystyle g_{1}$
$\displaystyle=C_{1}\mathrm{e}^{\mathrm{i}k_{1}x}\frac{\Theta_{1}(\frac{x}{\sqrt{m}}-\mathrm{i}z)}{\Theta_{4}(\mathrm{i}z)\Theta_{4}(\frac{x}{\sqrt{m}})},\
k_{1}=\frac{\omega-2\lambda}{4}-\frac{\mathrm{i}}{\sqrt{m}}Z(\mathrm{i}z),$
(B.5) $\displaystyle g_{2}$
$\displaystyle=C_{2}\mathrm{e}^{\mathrm{i}k_{2}x}\frac{\Theta_{2}(\frac{x}{\sqrt{m}}-\mathrm{i}z)}{\Theta_{3}(\mathrm{i}z)\Theta_{4}(\frac{x}{\sqrt{m}})},\
k_{2}=\frac{\omega+2\lambda}{4}-\frac{\mathrm{i}}{\sqrt{m}}Z(\mathrm{i}z-K).$
(B.6)

Due to the Bloch (Floquet) theorem, $g_{1}$ and $g_{2}$ must share the same
crystal momentum, i.e., $k_{1}$ and $k_{2}$ must be related by $k_{1}\equiv
k_{2}\mod\frac{\pi}{2\sqrt{m}K}$. In fact, using the periodicity of the Jacobi
elliptic and zeta functions, we can check $k_{1}=k_{2}$. Therefore, henceforth
we simply write $k_{1}=k_{2}=k(z)$. The ratio $C_{1}/C_{2}$ must also be
fixed. Rewriting $g_{1}/g_{2}$ using the Jacobi elliptic function, and
checking its consistency with the linear relation between $g_{1}$ and $g_{2}$
obtained from the eigenequation of $\hat{B}$, we find $C_{1}/C_{2}=1$.
The crystal momentum $k(z)$ can be further simplified as follows. The
derivative $k^{\prime}(z)$ can be simplified using the double-angle formula of
the $\operatorname{cs}$ function:

$\displaystyle
k^{\prime}(z)=\frac{1}{\sqrt{m}}\left(-\operatorname{cs}^{2}(2\mathrm{i}z)-\frac{E}{K}\right)=\frac{1}{2\sqrt{m}}\left(\operatorname{dn}^{2}(2\mathrm{i}z+\mathrm{i}K^{\prime})+\operatorname{dn}^{2}(2\mathrm{i}z-\mathrm{i}K^{\prime})-\frac{2E}{K}\right).$
(B.7)

Since $Z^{\prime}=\operatorname{dn}^{2}-\frac{E}{K}$, integrating this
expression soon gives the Jacobi zeta function. Furthermore, the constant of
integration can be fixed using some specific value at any point, e.g.,
$k_{1}(\frac{K^{\prime}}{2})=-\frac{\pi}{4K\sqrt{m}}$. Thus we obtain Eq.
(4.8).
The third and fourth components of the eigenfunction,
$g_{3}=\lambda^{-1}\left(-\operatorname{sn}(\frac{x}{\sqrt{m}})g_{1}-\operatorname{cn}(\frac{x}{\sqrt{m}})g_{2}\right)$
and
$g_{4}=\lambda^{-1}\left(\operatorname{cn}(\frac{x}{\sqrt{m}})g_{1}-\operatorname{sn}(\frac{x}{\sqrt{m}})g_{2}\right)$,
can be simplified using another expression of Eq. (4.6),
$\lambda(z)=-\mathrm{i}\frac{\Theta_{1}(\mathrm{i}z)\Theta_{2}(\mathrm{i}z)}{\Theta_{3}(\mathrm{i}z)\Theta_{4}(\mathrm{i}z)}$,
and the addition formula of theta functions, in particular, Ref. [20], Eqs.
(3.5b) and (3.8). The final expression is then given by Eq. (4.9).

##  Appendix C Completeness relation

In this appendix, we derive the completeness relation (4.15) for
eigenfunctions of the soliton lattice potential $\phi_{0}(x)$ of Eq. (4.3).333
In order to eliminate confusion, we should explain the difference of the usage
of the terminology “eigenvalues an eigenfunctions” between classical
integrable systems and other fields. In theoretical physics, particularly in
quantum mechanics, when one finds a solution of an eigenequation of a given
differential operator $\hat{L}f=\lambda f$ , the operand function $f$ is
called an eigenfunction only when it is a bounded function and only in this
case $\lambda$ is included to the set of eigenvalues. The eigenfunctions
possessing everywhere-bounded plane-wave type behavior are called a scattering
state and constitute the set of continuous eigenvalues, i.e., a band. The
eigenfunction with localized profile and finite norm is called a bound state
and make a discrete eigenvalue. The completeness relation, which is in bra-ket
notation of quantum mechanics often expressed as $1=\sum_{n}\ket{n}\bra{n}$,
is constructed by gathering all these scattering and bound eigenstates.
In classical integrable systems, the exponentially divergent solution of an
eigenequation, which is usually not counted as an eigenfunction, plays an
important role in generating multi-soliton solutions by various methods. An
elementary example is the divergent solution of the Schrödinger operator
$-f_{xx}=-\kappa^{2}f$ with eigenvalue $\lambda=-\kappa^{2}<0$ and
$f=e^{\pm\kappa x}$, which is indeed used to construct the multi-soliton
solution of the KdV equation. Thus, identification of all solutions of the
eigenequation for any complex $\lambda$ is important, and therefore, in the
context where no confusion occurs, these divergent solutions $f$ and
corresponding values $\lambda$ are also sometimes called eigenfunctions and
eigenvalues, without distinction. This manuscript also adopts this loose use
of terminology in several sections, e.g., in Sec. 4. To preserve the logical
unambiguity, the “genuine” eigenfunctions with no divergent behaviors are
explicitly referred to as scattering and bound eigenstates. The completeness
relation includes only these “genuine” eigenfunctions. The solution of
eigenequation for arbitrary complex eigenvalue is also often called the Baker-
Akhiezer function, whose original meaning is the single-valued function
defined on a Riemann surface and possessing finitely many essential singular
points but meromorphic excepting these points.
The first important point is that $\hat{L}$ is not self-adjoint, i.e.,
$\hat{L}^{\dagger}\neq\hat{L}$ and instead satisfies
$\hat{L}^{\dagger}=\sigma\hat{L}\sigma$ with $\sigma=I_{2}\otimes\sigma_{3}$
as stated in Sec. 3. For such operator, the eigenfunction satisfies the
orthogonal relation with respect to the indefinite inner product
$(f,g)_{\sigma}=\int dxf^{\dagger}\sigma g$. The normalization of
eigenfunction is also made based on $(f,f)_{\sigma}$ and all eigenfunctions
are classified into the ones possessing positive, negative, and zero norm. The
completeness relation for the case of finite-dimensional linear algebra is
shown in section 3 of Ref. [21], and we now need to consider its analog in
continuous space for a differential operator. (A specific example for a
differential operator in continuous space is found in Refs. [29, 30], though
not fully general.) In mathematical physics, a space equipped with indefinite
inner product is called the Krein space. Eigenfunctions with nonzero norm can
be almost analogously treated to the case of self-adjoint operators whose
eigenvalues are real and eigenfunctions are normalized by positive-definite
norm. On the other hand, we must carefully consider eigenfunctions with zero
norm. In particular, the eigenfunction with complex eigenvalue always has zero
norm due to Eq. (3.2). In this case, in order to prepare a $\sigma$-orthogonal
basis, we should construct positive- and negative-norm functions by linear
combination of a pair of zero-norm eigenfunctions with mutually complex
conjugate eigenvalues.
The second important point in writing down the completeness relation supposed
to be applied in the ISM is that the integrand must be expressed by
meromorphic function with respect to spectral variable $\lambda$ or its
parametrization variable $z$, because, in the derivation of the GLM equation
in D, we want to use the residue theorem. Therefore, we must eliminate $z^{*}$
from the integrand using the complex conjugation relations described in Eqs.
(4.10)-(4.14). The key relation is as follows. A little calculation using the
addition formula of theta functions shows

$\displaystyle f_{0}(x,z^{*})^{\dagger}\sigma
f_{0}(x,z)=\sum_{i=1}^{4}(-1)^{i-1}g_{i}(x,z^{*})^{*}g_{i}(x,z)=\frac{-1}{\sqrt{m}}\left(\operatorname{dn}^{2}(\tfrac{x}{\sqrt{m}})+\operatorname{cs}^{2}(2\mathrm{i}z)\right).$
(C.1)

This quantity is meromorphic, i.e., only including $z$ due to Eq. (4.14). The
spatial average of this quantity is, using
$\overline{\operatorname{dn}^{2}}=\frac{E}{K}$ and recalling Eq. (B.7),

$\displaystyle\overline{f_{0}(x,z^{*})^{\dagger}\sigma
f_{0}(x,z)}=\frac{1}{\sqrt{m}}\left(-\frac{E}{K}-\operatorname{cs}^{2}(2\mathrm{i}z)\right)=k^{\prime}(z),$
(C.2)

which is used to normalize scattering eigenstates and rewrite the $k$-integral
to $z$-integral. As shown below, the spatial average of the norm density
$\overline{f_{0}(x,z)^{\dagger}\sigma f_{0}(x,z)}$ for scattering states all
reduces to Eq. (C.2).
Thirdly, if we have a $\sigma$-orthogonal basis of the Bloch-type functions
parametrized by crystal momentum $k$ and band index $\alpha$, and all of them
have non-vanishing norm and are normalized so that the spatial average becomes
$\overline{\phi_{\pm}(k,\alpha)^{\dagger}\sigma\phi_{\pm}(k,\alpha)}=\pm 1$,
where $\pm$ represents the sign of norm, the completeness relation is given by
$\sum_{\alpha}\int\frac{dk}{2\pi}(\phi_{+}(k,\alpha)\phi_{+}(k,\alpha)^{\dagger}-\phi_{-}(k,\alpha)\phi_{-}(k,\alpha)^{\dagger})\sigma=\text{Id}$.
This fact can be proved by considering the finite-length system where the
scattering states can be explicitly normalizable and their spectrum is
discretized, and then taking the infinite-length limit. Here we omit the
detail of this argument, since it is tedious but straightforward. For readers’
reference, we note that a similar discussion in another physical problem for
discretization of scattering states of a self-adjoint operator and its
infinite limit can be found in Ref. [31], section 3.
Keeping in mind the above-mentioned things, let us now write down the
completeness relation. Since the potentials (4.4) are periodic, there are only
scattering states by the Bloch theorem and no bound state exist. Hence, the
completeness relation is written only by scattering states. As stated in Sec.
4, up to periodicity, there are four distinct lines in $z$-plane where crystal
momentum becomes real $k(z)\in\mathbb{R}$, representing the scattering states.
Here we discuss the scattering states on each line in detail.

  1. (i)

$z\in\mathbb{R}$.
We can check $\lambda(z)\in\mathbb{R}$ and $f_{0}(x,z)$ has positive norm, so
$f_{0}(x,z^{*})=f_{0}(x,z)$ and hence
$\overline{f_{0}(x,z^{*})^{\dagger}\sigma
f_{0}(x,z)}=\overline{f_{0}(x,z)^{\dagger}\sigma f_{0}(x,z)}=k^{\prime}(z)>0$.
Therefore, the contribution of these states to completeness relation is
$\int\frac{dk}{2\pi}\frac{f_{0}(x,z)f_{0}(y,z)^{\dagger}\sigma}{k^{\prime}(z)}=\int_{-K^{\prime}}^{K^{\prime}}\frac{dz}{2\pi}f_{0}(x,z)f_{0}(y,z^{*})^{\dagger}\sigma$.

  2. (ii)

$z\in\mathbb{R}-\mathrm{i}K$.
We can check $\lambda(z)\in\mathbb{R}$ and $f_{0}(x,z)$ has negative norm, and
$f_{0}(x,z^{*})=-f_{0}(x,z)$ using periodicity (4.13). Thus,
$\overline{f_{0}(x,z)^{\dagger}\sigma
f_{0}(x,z)}=-\overline{f_{0}(x,z^{*})^{\dagger}\sigma
f_{0}(x,z)}=-k^{\prime}(z)<0$. Therefore, the contribution to these states to
completeness relation is $\int
dk\frac{f_{0}(x,z)f_{0}(y,z)^{\dagger}\sigma}{-k^{\prime}(z)}=\int_{-K^{\prime}-\mathrm{i}K}^{K^{\prime}+\mathrm{i}K}\frac{dz}{2\pi}f_{0}(x,z)f_{0}(y,z^{*})^{\dagger}\sigma.$

  3. (iii)

$z\in\mathbb{R}\pm\frac{\mathrm{i}K}{2}$.
We can check $\lambda(z)\in\mathrm{i}\mathbb{R}$ and hence $f_{0}(x,z)$ has
zero norm. Therefore we must make nonzero-norm functions by linear combination
of a pair of eigenfunctions with mutually complex conjugate eigenvalues. Since
$\lambda(z^{*})=\lambda(z)^{*}$ and $k(z^{*})=k(z)^{*}$, this pair share the
same wavenumber. Now let us assume that $f_{1}$ and $f_{2}$ are zero-norm
eigenfunctions of the Bloch type with complex eigenvalues $\lambda$ and
$\lambda^{*}$ and share the same and real-valued wavenumber $k$. By definition
the spatial average of norm density is zero: $\overline{f_{i}^{\dagger}\sigma
f_{i}}=0,\ i=1,2$, but $\overline{f_{1}^{\dagger}\sigma f_{2}}\neq 0$ unless
they belong to nontrivial Jordan blocks. Then, if an overall prefactors of
$f_{1}$ and $f_{2}$ are adjusted to satisfy the relation
$\overline{f_{1}^{\dagger}\sigma f_{2}}=\overline{f_{2}^{\dagger}\sigma
f_{1}}>0$, then $f_{\pm}=\frac{f_{1}\pm f_{2}}{\sqrt{2}}$ are positive- and
negative-norm functions and $\sigma$-orthogonal to each other. Thus these two
can be used as a $\sigma$-orthogonal basis, and the contribution of these
states to completeness relation is written as
$\int\frac{dk}{2\pi}\frac{f_{+}f_{+}^{\dagger}\sigma-
f_{-}f_{-}^{\dagger}\sigma}{\overline{f_{2}^{\dagger}\sigma
f_{1}}}=\int\frac{dk}{2\pi}\frac{(f_{1}f_{2}^{\dagger}+f_{2}f_{1}^{\dagger})\sigma}{\overline{f_{2}^{\dagger}\sigma
f_{1}}}.$ Now let us take $f_{1}=f_{0}(x,z),\
z\in\mathbb{R}+\frac{\mathrm{i}K}{2}$ and $f_{2}=-f_{0}(x,z^{*}),\
z^{*}\in\mathbb{R}-\frac{\mathrm{i}K}{2}$. By Eq. (C.2),
$\overline{f_{2}^{\dagger}\sigma f_{1}}=-k^{\prime}(z)>0$ is real and
positive. (Note: $k(z)$ is a decreasing function on
$z\in\mathbb{R}\pm\frac{\mathrm{i}K}{2}$.) Therefore, the contribution of
these states to completeness relation is given by
$\int\frac{dk}{2\pi}\frac{(-f_{0}(x,z)f_{0}(y,z^{*})^{\dagger}-f_{0}(x,z^{*})f_{0}(y,z)^{\dagger})\sigma}{-k^{\prime}(z)}=-\int_{-K^{\prime}+\frac{\mathrm{i}K}{2}}^{K^{\prime}+\frac{\mathrm{i}K}{2}}\frac{dz}{2\pi}(f_{0}(x,z)f_{0}(y,z^{*})^{\dagger}+f_{0}(x,z^{*})f_{0}(y,z)^{\dagger})\sigma$.
Changing the dummy variable $z\to z^{*}$, the latter term can be rewritten as
$-\int_{-K^{\prime}-\frac{\mathrm{i}K}{2}}^{K^{\prime}-\frac{\mathrm{i}K}{2}}\frac{dz}{2\pi}f_{0}(x,z)f_{0}(y,z^{*})^{\dagger}\sigma$.

Summarizing all contributions from (i)-(iii), the completeness relation is
given by

$\displaystyle\delta(x-y)I_{4}$
$\displaystyle=\int_{-K^{\prime}}^{K^{\prime}}\frac{dz}{2\pi}f_{0}(x,z)f_{0}(y,z^{*})^{\dagger}\sigma+\int_{-K^{\prime}-\mathrm{i}K}^{K^{\prime}+\mathrm{i}K}\frac{dz}{2\pi}f_{0}(x,z)f_{0}(y,z^{*})^{\dagger}\sigma$
$\displaystyle\qquad-\int_{-K^{\prime}+\frac{\mathrm{i}K}{2}}^{K^{\prime}+\frac{\mathrm{i}K}{2}}\frac{dz}{2\pi}f_{0}(x,z)f_{0}(y,z^{*})^{\dagger}\sigma-\int_{-K^{\prime}-\frac{\mathrm{i}K}{2}}^{K^{\prime}-\frac{\mathrm{i}K}{2}}\frac{dz}{2\pi}f_{0}(x,z)f_{0}(y,z^{*})^{\dagger}\sigma.$
(C.3)

Adding vertical paths which cancel out by periodicity, we obtain the desired
form

$\displaystyle\delta(x-y)I_{4}=\oint_{C_{1}+C_{2}}\frac{dz}{2\pi}f(x,z)f(y,z^{*})^{\dagger}\sigma,$
(C.4)

where the definitions of the contours $C_{1}$ and $C_{2}$ are already
described in the sentences after Eq. (4.15).

##  Appendix D Derivation of the GLM equation

In this appendix we drive the GLM equation (6.9). By the relation between
right and left Jost solutions (5.4) and the integral representation (6.1),

$\displaystyle
f_{+}(x,z)\frac{1}{a(z)}-f_{0}(x,z)=\int_{-\infty}^{x}\mathrm{d}y\Gamma(x,y)f_{0}(y,z)+\left(f_{0}(x,-z-\mathrm{i}K)+\int_{-\infty}^{x}\mathrm{d}y\Gamma(x,y)f_{0}(y,-z-\mathrm{i}K)\right)\frac{b(z)}{a(z)}.$
(D.1)

Multiplying $\frac{f(w,z^{*})^{\dagger}\sigma}{2\pi},\ w\leq x$ to both sides
of the above equation from right, we integrate the expression by $z$ along the
contour $C_{1}+C_{2}$. First, let us evaluate the right-hand side. Defining

$\displaystyle\Omega_{\text{sc}}(x,w):=\int_{C_{1}+C_{2}}\frac{\mathrm{d}z}{2\pi}f_{0}(x,-z-\mathrm{i}K)\frac{b(z)}{a(z)}f_{0}(w,z^{*})^{\dagger}\sigma,$
(D.2)

and using the completeness relation (4.15), we find

$\displaystyle\text{R.H.S.}=\Gamma(x,w)+\Omega_{\text{sc}}(x,w)+\int_{-\infty}^{x}\mathrm{d}y\Gamma(x,y)\Omega_{\text{sc}}(y,w).$
(D.3)

In the contour integral of $\Omega_{\text{sc}}$, actually the integrals for
scattering states only contribute, and hence it can be rewritten as Eq.
(6.12).
Next, let us evaluate the left-hand side. In the SG equation, generally $a(z)$
can have higher-order zeros [24], but here we only consider the case where all
zeros are of the first order. Let zeros of $a(z)$ be $z_{1},\dots,z_{n}$. As
already mentioned before, by Eq. (5.11), the zeros of $a(z)$ simultaneously
emerge at $z,-z^{*},z-\mathrm{i}K,-\mathrm{z}^{*}-\mathrm{i}K$. Therefore, not
all $z_{j}$’s can be freely chosen. Here, however, we formally assign
independent labels to all zeros. Then, the left-hand side can be calculated
using the residue theorem:

L.H.S.
$\displaystyle=\sum_{j=1}^{n}\frac{\mathrm{i}f_{+}(x,z_{j})f_{0}(w,z_{j}^{*})^{\dagger}}{\dot{a}(z_{j})}$
$\displaystyle=\sum_{j=1}^{n}\frac{\mathrm{i}f_{-}(x,-z_{j}-\mathrm{i}K)b(z_{j})f_{0}(w,z_{j}^{*})^{\dagger}}{\dot{a}(z_{j})},$
(D.4)

where the relation $f_{+}(x,z_{j})=b(z_{j})f_{-}(x,-z_{j}-\mathrm{i}K)$
arising from the condition $a(z_{j})=0$ has been used. If we define

$\displaystyle\Omega_{\text{bd}}(x,w):=\sum_{j=1}^{n}f_{0}(x,-z_{j}-\mathrm{i}K)c_{j}^{2}f_{0}(w,z_{j}^{*})^{\dagger}\sigma,\quad
c_{j}^{2}:=\frac{b(z_{j})}{\mathrm{i}\dot{a}(z_{j})},$ (D.5)

then

$\displaystyle\text{L.H.S.}=-\Omega_{\text{bd}}(x,w)-\int_{-\infty}^{x}\mathrm{d}y\Gamma(x,y)\Omega_{\text{bd}}(y,w).$
(D.6)

To summarize, writing $\Omega=\Omega_{\text{sc}}+\Omega_{\text{bd}}$, we
obtain the final result

$\displaystyle\Gamma(x,w)+\Omega(x,w)+\int_{-\infty}^{x}\mathrm{d}y\Gamma(x,y)\Omega(y,w)=0,\quad
w\leq x.$ (D.7)

If we simplify Eq. (D.5) using Eqs. (4.13) and (4.14), it is rewritten as Eq.
(6.10) with (6.11).

##  Appendix E Calculation of determinant

Here we calculate the asymptotic form of the multi-soliton solution for
$x\to+\infty$. In particular, we provide the shift of the lattice $x_{0}$ in
Eq. (7.19). First, we write $\mathcal{P}=S_{1}\tilde{\mathcal{P}}S_{2}$ and
$\mathcal{Q}=S_{3}\tilde{\mathcal{Q}}S_{4}$, where

$\displaystyle S_{1}$
$\displaystyle=\operatorname{diag}\left(\frac{\Theta_{4}(\mathrm{i}z_{1})}{\Theta_{1}(\mathrm{i}z_{1})},\dots,\frac{\Theta_{4}(\mathrm{i}z_{n})}{\Theta_{1}(\mathrm{i}z_{n})}\right),\quad$
$\displaystyle S_{2}$
$\displaystyle=\operatorname{diag}\left(\frac{\Theta_{2}(\mathrm{i}z_{1})}{\Theta_{3}(\mathrm{i}z_{1})},\dots,\frac{\Theta_{2}(\mathrm{i}z_{n})}{\Theta_{3}(\mathrm{i}z_{n})}\right),$
(E.1) $\displaystyle S_{3}$
$\displaystyle=\operatorname{diag}\left(\frac{\Theta_{4}(\mathrm{i}z_{1})}{\Theta_{2}(\mathrm{i}z_{1})},\dots,\frac{\Theta_{4}(\mathrm{i}z_{n})}{\Theta_{2}(\mathrm{i}z_{n})}\right),\quad$
$\displaystyle S_{4}$
$\displaystyle=\operatorname{diag}\left(\frac{\Theta_{1}(\mathrm{i}z_{1})}{\Theta_{3}(\mathrm{i}z_{1})},\dots,\frac{\Theta_{1}(\mathrm{i}z_{n})}{\Theta_{3}(\mathrm{i}z_{n})}\right),$
(E.2)

and

$\displaystyle\tilde{\mathcal{P}}_{ij}$
$\displaystyle=\frac{-\Theta_{4}\Theta_{2}(\frac{x}{\sqrt{m}}+\mathrm{i}(z_{i}+z_{j}))}{\Theta_{2}(\mathrm{i}(z_{i}+z_{j}))\Theta_{1}(\frac{x}{\sqrt{m}})},$
(E.3) $\displaystyle\tilde{\mathcal{Q}}_{ij}$
$\displaystyle=\frac{\Theta_{4}\Theta_{1}(\frac{x}{\sqrt{m}}+\mathrm{i}(z_{i}+z_{j}))}{\Theta_{2}(\mathrm{i}(z_{i}+z_{j}))\Theta_{2}(\frac{x}{\sqrt{m}})}.$
(E.4)

Then, using the formula in Ref. [27], appendix D, we obtain

$\displaystyle\frac{\det\tilde{\mathcal{P}}}{\det\mathcal{G}}$
$\displaystyle=(-1)^{n}\frac{\Theta_{4}(\frac{x}{\sqrt{m}})\Theta_{1}(\frac{x}{\sqrt{m}}+nK+2\mathrm{i}\sum_{j=1}^{n}z_{j})}{\Theta_{1}(\frac{x}{\sqrt{m}})\Theta_{4}(\frac{x}{\sqrt{m}}+nK+2\mathrm{i}\sum_{j=1}^{n}z_{j})},$
(E.5) $\displaystyle\frac{\det\tilde{\mathcal{Q}}}{\det\mathcal{G}}$
$\displaystyle=(-1)^{n}\frac{\Theta_{4}(\frac{x}{\sqrt{m}})\Theta_{2}(\frac{x}{\sqrt{m}}+nK+2\mathrm{i}\sum_{j=1}^{n}z_{j})}{\Theta_{2}(\frac{x}{\sqrt{m}})\Theta_{4}(\frac{x}{\sqrt{m}}+nK+2\mathrm{i}\sum_{j=1}^{n}z_{j})}.$
(E.6)

This relation itself is generally valid even when all $z_{j}$’s are
independent and $n$ is not even.
As discussed in Sec. 9, for the real and bounded solution of the SG equation,
$n$ is even and $z_{j}$ and $z_{j}-\mathrm{i}K$ simultaneously appear. In this
case, we label the zeros, e.g., as $z_{n/2+j}=z_{j}-\mathrm{i}K,$
$j=\frac{n}{2}+1,\dots,n$. Then, $\det
S_{1}=\prod_{j=1}^{n/2}\frac{\Theta_{4}(\mathrm{i}z_{j})\Theta_{3}(\mathrm{i}z_{j})}{\Theta_{1}(\mathrm{i}z_{j})\Theta_{2}(\mathrm{i}z_{j})},\quad\det
S_{2}=(-1)^{n/2}\det S_{1}^{-1},\quad\det
S_{4}=\prod_{j=1}^{n/2}\frac{\Theta_{1}(\mathrm{i}z_{j})\Theta_{2}(\mathrm{i}z_{j})}{\Theta_{3}(\mathrm{i}z_{j})\Theta_{4}(\mathrm{i}z_{j})},\quad\det
S_{3}=(-1)^{n/2}\det S_{4}^{-1}$, and hence

$\displaystyle\frac{\det\mathcal{P}}{\det\mathcal{G}}$
$\displaystyle=\frac{\Theta_{4}(\frac{x}{\sqrt{m}})\Theta_{1}(\frac{x}{\sqrt{m}}+2\mathrm{i}\sum_{j=1}^{n}z_{j})}{\Theta_{1}(\frac{x}{\sqrt{m}})\Theta_{4}(\frac{x}{\sqrt{m}}+2\mathrm{i}\sum_{j=1}^{n}z_{j})}=\frac{\cos\frac{\phi_{0}(x-x_{0})}{2}}{\cos\frac{\phi_{0}(x)}{2}},$
(E.7) $\displaystyle\frac{\det\mathcal{Q}}{\det\mathcal{G}}$
$\displaystyle=\frac{\Theta_{4}(\frac{x}{\sqrt{m}})\Theta_{2}(\frac{x}{\sqrt{m}}+2\mathrm{i}\sum_{j=1}^{n}z_{j})}{\Theta_{2}(\frac{x}{\sqrt{m}})\Theta_{4}(\frac{x}{\sqrt{m}}+2\mathrm{i}\sum_{j=1}^{n}z_{j})}=\frac{\sin\frac{\phi_{0}(x-x_{0})}{2}}{\sin\frac{\phi_{0}(x)}{2}}.$
(E.8)

## References

  * [1] G. L. Lamb, Jr.,  Elements of soliton theory,  (Wiley, Hoboken, New Jersey, 1980).
  * [2] P. M. Chaikin and T. C. Lubensky,  Principles of Condensed Matter Physics,  (Cambridge University Press, Cambridge, 2000).
  * [3] Y. Togawa, T. Koyama, K. Takayanagi, S. Mori, Y. Kousaka, J. Akimitsu, S. Nishihara, K. Inoue, A. S. Ovchinnikov, and J. Kishine, Phys. Rev. Lett., 108, 107202 (2012).
  * [4] A. B. Borisov, J. Kishine, I. G. Bostrem, and A. S. Ovchinnikov, Phys. Rev. B, 79, 134436 (2009).
  * [5] T. Fulton and R. Dynes, Solid State Communications, 12(1), 57–61 (1973).
  * [6] P. S. Lomdahl, O. H. Soerensen, and P. L. Christiansen, Phys. Rev. B, 25, 5737–5748 (1982).
  * [7] T. Kato and M. Imada, J. Phys. Soc. Jpn., 66(5), 1445–1449 (1997).
  * [8] A. B. Borisov and V. V. Kiseliev, Solid State Communications, 59(7), 445–448 (1986).
  * [9] A. B. Borisov and V. V. Kiseliev, Physica D: Nonlinear Phenomena, 31(1), 49–64 (1988).
  * [10] A. B. Borisov and V. V. Kiseliev, Inverse Problems, 5(6), 959–982 (1989).
  * [11] D. E. Pelinovsky and R. E. White, Proc. R. Soc. A: Math. Phys. Eng. Sci., 476(2242) (2020).
  * [12] M. J. Ablowitz, D. J. Kaup, A. C. Newell, and H. Segur, Phys. Rev. Lett., 30, 1262–1264 (1973).
  * [13] M. J. Ablowitz, D. J. Kaup, A. C. Newell, and H. Segur, Stud. Appl. Math., 53, 249 (1974).
  * [14] A. S. Budagov and L. A. Takhtadzhyan, Dokl. Akad. Nauk SSSR, 235, 805–808 (1977).
  * [15] L. D. Faddeev and L. A. Takhtajan,  Hamiltonian Methods in the Theory of Solitons,  (Springer, Berlin, 1987).
  * [16] E. D. Belokolos, A. I. Bobenko, V. Z. Enol’skii, A. R. Its, and V. B. Matveev,  Algebro-Geometric Approach to Nonlinear Integrable Equations,  (Springer, Berlin, 1994).
  * [17] L. A. Takhtadzhyan and L. D. Faddeev, Theor. Math. Phys., 21(2), 1046–1057 (1974).
  * [18] L. A. Takhtadzhyan, Sov. Phys. JETP, 39, 228 (1974).
  * [19] I. M. Krichever, Funct. Anal. Appl., 11(1), 12–26 (1977).
  * [20] S. Kharchev and A. Zabrodin, J. Geom. Phys., 94, 19–31 (2015).
  * [21] D. A. Takahashi and M. Nitta, Ann. Phys., 354, 101–156 (2015).
  * [22] M. Abramowitz and I. A. Stegun,  Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables,  (Dover, Mineola, New York, 9th edition, 1965).
  * [23] D. A. Takahashi, Prog. Theor. Exp. Phys., 2016(4), 043I01 (2016).
  * [24] H. Tsuru and M. Wadati, J. Phys. Soc. Jpn., 53(9), 2908–2921 (1984).
  * [25] T. Yajima and M. Wadati, J. Phys. Soc. Jpn., 56(9), 3069–3081 (1987).
  * [26] V. E. Zakharov and A. B. Shabat, Funct. Anal. Its Appl., 8, 226–235 (1974).
  * [27] D. A. Takahashi, Phys. Rev. E, 93, 062224 (2016).
  * [28] A. Hurwitz and R. Courant,  Vorlesungen Über Allgemeine Funktionentheorie Und Elliptische Funktionen,  (Springer-Verlag, Berlin, Heidelberg, 4th edition, 1964).
  * [29] D. A. Takahashi, M. Kobayashi, and M. Nitta, Phys. Rev. B, 91, 184501 (2015).
  * [30] D. A. Takahashi, K. Ohashi, T. Fujimori, and M. Nitta, Phys. Rev. A, 96, 023626 (2017).
  * [31] D. A. Takahashi and M. Nitta, J. Low Temp. Phys., 175(1-2), 250–257 (2014).