Majorana zero mode-soliton duality and in-gap and BIC bound states in modified
Toda model coupled to fermion

H. Blas(a), J.J. Monsalve(a) , R. Quicaño(b) and J. R. V. Pereira(a)

(a)Instituto de Física

Universidade Federal de Mato Grosso

Av. Fernando Correa, $N^{0}$ 2367

Bairro Boa Esperança, Cep 78060-900, Cuiabá - MT - Brazil

(b)Instituto de Matemática y Ciencias Afines-IMCA

Facultad de Ciencias

Universidad Nacional de Ingeniería

Av. Tupac Amaru, $N^{0}$ 210, Rimac, Lima-Perú

<EMAIL_ADDRESS>
A two-dimensional field theory of a fermion chirally coupled to Toda field
plus a scalar self-coupling potential is considered. Using techniques of
integrable systems we obtain analytical zero modes, in-gap states and bound
states in the continuum (BIC) for topological configurations of the scalar
field. Fermion-soliton duality mappings are uncovered for the bound state
spectrum, which interpolates the weak and strong coupling sectors of the model
and give rise to novel Thirring-like and multi-frequency sine-Gordon models,
respectively. The non-perturbative effects of the back-reaction of the fermion
bound states on the kink are studied and it is shown that the zero mode would
catalyze the emergence of a new kink with lower topological charge and greater
slope at the center, in the strong coupling limit of the model. For special
topological charges and certain relative phases of the fermion components it
is shown that the kinks can host Majorana zero modes. Our results may find
applications in several branches of non-linear physics, such as confinement in
QCD2, braneworld models, high $T_{c}$ superconductivity and topological
quantum computation. We back up our results with numerical simulations for
continuous families of topological sectors.

##  1 Introduction

The study of fermion-soliton systems have attracted significant interest in
the context of the localization and bound states of the fermion, since the
interaction of the fermions and the topological solitons are related to
interesting physical effects, such as the fermion number fractionization and
polarization of the fermionic vacuum [1, 2, 3], string superconductivity in
cosmology [4], and monopole catalysis of the baryon decay [5, 6]. Moreover, in
braneworld models, following the original work [7], some properties have been
studied, such as the localization of bulk chiral fermions on a brane [8]. Some
brane Lagrangians consider the kink soliton of the $\phi^{4}$ or sine-Gordon
type systems in one dimension (see e.g. [9]). In this context the particle
transfer and fermion spectrum in braneworld collisions have also been
considered [10, 11]. In addition, the kink field deforms the energy levels of
the fermionic vacuum leading to bound states and modified continuum states as
compared to the free fermion spectrum. These change the zero-point fermion
energy, giving rise to the Casimir effect in the presence of the soliton [12].

In most of the treatments the back-reaction of the fermion on the kink has
been neglected and the kink was taken as prescribed. Since the models of
physical interest are in general non-integrable they rarely possess analytical
solutions. So, some recent studies resorted to numerical simulations in order
to address the back-reaction of the fermion on the kink by computing self-
consistent solutions of the system [13, 14].

Many aspects of soliton-fermion systems in two dimensions, such as the quantum
corrections, have been studied in the literature; however, some important
features at the zero order are yet to be investigated. For instance, in the
earlier studies there have been assumed that only the bosonic self-interaction
part of the potential to be responsible for the creation of the kink
topological number. However, in those models the dynamics of the fermion-kink
coupling sector must be taken into account, since this term might contribute
to deform the asymptotic values of the kink in the strong coupling limit. This
requires the study of the self-consistent solutions uncovering the transition
from weak to strong coupling sectors of the model.

In this paper we consider a modification of the so-called $sl(2)$ affine Toda
model coupled to matter field (ATM) by introducing an additional potential of
self-interacting Toda field. The integrability and the soliton solutions
associated to some particular vacua of the ATM model have been studied in [15,
16]. One of its main properties is that for certain solutions there is an
equivalence between the $U(1)$ Noether current, involving only the Dirac
field, and the topological current associated to the pseudo-scalar soliton.
This property was established in [15] at the classical level and in [16] at
the quantum mechanical level. It implies that the density of the $U(1)$ charge
has to be concentrated in the regions where the pseudo-scalar field has non-
vanishing slope. The one-soliton of the theory is a kink type soliton with
fractional topological charge, and therefore the charge density is
concentrated inside the soliton. The general $sl(n)$ ATM model constitutes an
excellent laboratory to test ideas about confinement [16, 17, 18, 19], the
role of solitons in quantum field theories [16] and duality transformations
interchanging solitons and particles [17, 20, 21].

The ATM model has also been considered in many physical situations related to
fermion-soliton systems. It appears in the study of the dynamical Peierls’
energy gap generation in one-dimensional charge-density wave systems [22],
quantum field theory description of tunneling in the integer quantum Hall
effect [23], discussions of fractional charges induced in the ground state of
a Fermi system by its coupling to a chiral field [2, 24, 25, 26, 3, 27], the
spectrum and string tension of the low-energy effective Lagrangian of QCD2
[18] and a two dimensional model of high $T_{c}$ superconductivity [28]. The
model has also been considered in order to discuss the realization of chiral
symmetries in $(1+1)$ dimensions [29]. Numerical self-consistent solutions
have been discussed for the model with additional $\phi^{4}$ potential [30]
and the Casimir energy of the system has been computed for a prescribed kink
[12].

Here, using some techniques of integrable systems, we perform the study of the
soliton and fermion bound state spectra of a modified ATM model. Although this
model is non-integrable in general, its analytic soliton type solutions
associated to specific vacua, can be obtained using the Hirota’s tau function
formalism. The static version of the model is considered and the soliton and
fermion bound state spectra is uncovered. In this context, the in-gap and
continuum bound states (BIC), as well as the fermion zero-modes, are examined.
These soliton solutions allow us to perform a mapping to either modified sine-
Gordon or novel Thirring-like models, in the strong and weak coupling sectors,
respectively. We perform this correspondence for each of the sectors, i.e. in-
gap, BIC and zero-mode bound states corresponding to a family of soliton
profiles of the scalar field. So, one gets extensions and modifications of the
duality mappings of the usual ATM model [17, 20]. As a by product of the
process, one finds novel scalar and fermionic models, describing the dual
sectors of the modified ATM model, respectively. The equivalence between the
$U(1)$ and topological charges of the usual ATM model gets an additional
contribution as a non-local charge density due to the deformation potential
which depends only on the Toda field.

For special topological charges and certain relative phases of the fermion
components it is shown that the kinks would host Majorana zero modes. The
mapping between the kink and the Majorana zero mode is constructed and shown
to describe either a double sine-Gordon or a novel real fermion Thirring-like
model in the strong and weak coupling sector of the model, respectively. The
back-reaction of the Majorana bound state on the kink is shown to catalyze the
emergence of a new kink with lower topological charge and greater slope at the
center, in the strong coupling limit of the model. However, for the model with
vanishing self-coupling potential and in the intermediate coupling sector the
Majorana fermion-soliton system exhibits an equivalence between the
topological charge density and the Majorana state density. In the presence of
the scalar potential this equivalence gets an additional non-local
contribution depending on the Toda field.

The choice of complicated background configurations associated to the scalar
field renders the problem analytically unsolvable. So, one could resort to
certain numerical simulations. However, in the numerical approach some
important physical properties and concepts might not be clear as in the
analytical approach described above. Nevertheless, in the numerical treatment
the scalar field can represent topological kinks with many different shapes
and topological charge taking any real value hosting fermion bound states. So,
in order to check our results we have performed various numerical simulations
and they support strongly our analytical approach. In order to perform the
numerical simulations we follow the relaxation method which is discussed in
the appendix F.

The paper is organized as follows: The section 2 summarizes the main
properties of a modified ATM model. In section 3 the properties of the fermion
coupled to topological kinks are examined. The conditions for the threshold,
BIC and Majorana bound states are established. A special scalar self-coupling
potential is introduced and its vacuum structure is studied in some detail. In
section 4 self-consistent kink and fermion bound state solutions are obtained
through the tau function approach. A finite set of topological charges and the
associated fermion bound states, as well as the energy spectrum are uncovered.
It is established a second order algebraic equation for the bound state energy
eigenvalue $E$. In section 5 we construct a fermion-soliton mapping and the
back-reaction of the fermion on the kink is studied, as well as the strong
coupling sector of the model. In section 6 we find the zero modes of the model
and the parameters of the Majorana zero modes are obtained. The results of our
numerical simulations are discussed in section 7. In section 8 we discuss our
main results. The appendices A, B, …,E present useful equations and identities
used in the tau function approach. The appendix F presents the algorithm
suitable for solving the model using the relaxation method.

##  2 The model

We consider a fermion chirally coupled to a scalar field defined by the
following Lagrangian111Our notation: $x_{\pm}=t\pm x$, and so,
$\partial_{\pm}=\frac{1}{2}(\partial_{t}\pm\partial_{x})$, and
$\partial^{2}=\partial^{2}_{t}-\partial^{2}_{x}=4\partial_{-}\partial_{+}$. We
use $\gamma_{0}=\left(\begin{array}[]{cc}0&i\\\ -i&0\end{array}\right)$,
$\gamma_{1}=\left(\begin{array}[]{cc}0&-i\\\ -i&0\end{array}\right)$,
$\gamma_{5}=\gamma_{0}\gamma_{1}=\left(\begin{array}[]{cc}1&0\\\
0&-1\end{array}\right)$, and $\psi=\left(\begin{array}[]{c}\psi_{R}\\\
\psi_{L}\end{array}\right),\,\,\bar{\psi}=\psi^{\dagger}\gamma_{0}$.

$\displaystyle{\cal L}={1\over
2}\partial_{\mu}\Phi\,\partial^{\mu}\Phi+i{\bar{\psi}}\gamma^{\mu}\partial_{\mu}\psi-M\,{\bar{\psi}}\,e^{i\beta\,\Phi\,\gamma_{5}}\,\psi-V(\Phi),$
(2.1)

with the potential

$\displaystyle V(\Phi)\equiv
A_{1}\cos{(\beta_{1}\Phi)}+A_{2}\cos{(\beta_{2}\Phi)}+A_{3}\cos{(\beta_{3}\Phi)}.$
(2.2)

The parameter $M$ defines a free fermion mass and $\beta$ the coupling
strength between the fermion and scalar. The function $V(\Phi)$ defines a
multi-frequency scalar potential and the real parameters
$\beta_{j}\,(j=1,2,3)$ define the self-coupling parameters, respectively.

The model (2.1) for vanishing potential $V=0$ belongs to a wide class of
integrable soliton theories, the so-called $\hat{sl}(2)$ affine Toda systems
coupled to matter fields (ATM) [15]. The integrability properties, the
construction of the general solution including the solitonic ones, the
soliton-fermion duality properties, as well as its symplectic structures were
discussed in [16, 17, 20]. This model has been shown to describe the low-
energy effective Lagrangian of QCD2 with one flavor and $N$ colors [18] and
the BCS coupling in spinless fermions in a two dimensional model of high T
superconductivity in which the solitons play the role of the Cooper pairs
[28]. Moreover, the model has earlier been studied as a model for fermion
confinement in a chiral invariant theory [31] and the mechanism of fermion
mass generation without spontaneously chiral symmetry breaking in two-
dimensions [29].

For an arbitrary potential $V$ and $\beta\neq 0$ the dynamical model would
lose the charge conjugation and energy-reflection symmetries, resulting in a
non-symmetrical energy spectrum with respect to the $E=0$ line. We will
consider the model (2.1) with potential (2.2) and show that the dynamics of
the system exhibits kink type and associated spinor bound states such that the
charge conjugation and energy-reflection symmetries hold.

The eqs. of motion become

$\displaystyle
i\gamma^{\mu}\partial_{\mu}\psi-M\,e^{i\beta\Phi\,\gamma_{5}}\,\psi$
$\displaystyle=$ $\displaystyle 0,$ (2.3)
$\displaystyle\partial^{2}\Phi-M\beta\bar{\psi}\psi\sin{(\beta\Phi)}+iM\beta\bar{\psi}\gamma_{5}\psi\cos{(\beta\Phi)}+V^{\prime}[\Phi]$
$\displaystyle=$ $\displaystyle
0,\,\,\,\,\,\,\,V^{\prime}[\Phi]\equiv\frac{dV(\Phi)}{d\Phi}.$ (2.4)

The global $U(1)$ transformation

$\displaystyle\Phi\rightarrow\Phi\;;\qquad\psi\rightarrow
e^{i\delta_{1}}\,\psi\;;\qquad{\psi^{*}}\rightarrow
e^{-i\delta_{1}}\,{\psi^{*}},$ (2.5)

gives rise to the next Noether current and the conservation law

$\displaystyle
J^{\mu}={\bar{\psi}}\,\gamma^{\mu}\,\psi,\,\,\,\partial_{\mu}\,J^{\mu}=0.$
(2.6)

Next, taking into account the equations of motion (2.3)-(2.4) one can write
the anomalous conservation law

$\displaystyle\partial_{\mu}J_{5}^{\mu}$ $\displaystyle=$
$\displaystyle-\frac{2}{\beta}V^{\prime}(\Phi).$ (2.7) $\displaystyle
J_{5}^{\mu}$ $\displaystyle\equiv$
$\displaystyle\bar{\psi}\gamma^{\mu}\gamma_{5}\psi+\frac{2}{\beta}\partial^{\mu}\Phi$
(2.8)

In fact, for vanishing potential, i.e. $V(\Phi)=0$, the Lagrangian (2.1) is
invariant under the global chiral symmetry

$\Phi\rightarrow\Phi-\frac{2\lambda_{1}}{\beta},\,\,\,\,\,\,\psi\rightarrow
e^{i\gamma_{5}\lambda_{1}}\,\psi\;;\qquad{\psi^{*}}\rightarrow
e^{-i\gamma_{5}\lambda_{1}}\,{\psi^{*}},$ (2.9)

which is associated to the conservation law

$\partial_{\mu}J_{5}^{\mu}=0,$ (2.10)

with the current $J_{5}^{\mu}$ defined in (2.8). So, a non-vanishing potential
breaks the global chiral symmetry above and the expression
$-\frac{2}{\beta}V^{\prime}(\Phi)$ in the r.h.s. of the quasi-conservation law
(2.7) can be dubbed as an anomaly.

Moreover, the Lagrangian (2.1) for vanishing potential, $V(\Phi)=0$, is
invariant under the transformation

$\displaystyle\Phi\rightarrow\Phi+\frac{2\pi
n}{\beta},\,\,\,n\in\leavevmode\hbox{\sf Z\kern-3.99994ptZ},$ (2.11)

with the fermion field unchanged. The vacua are infinitely degenerate and one
can define the topological current and its associated charge

$\displaystyle
j^{\mu}_{top}\equiv{\beta\over{2\pi}}\epsilon^{\mu\nu}\partial_{\nu}\,\Phi$
(2.12) $\displaystyle Q_{\rm topol.}\equiv\int\,dx\,j_{top}^{0}$ (2.13)

which takes non-zero values depending only on the asymptotic values of $\Phi$,
at $x=\pm\infty$. However, for non-vanishing potential the Lagrangian (2.1)
will be invariant under (2.11) provided that

$\displaystyle\beta_{j}=n_{j}\beta,\,\,\,n_{j}\in\leavevmode\hbox{\sf
Z\kern-3.99994ptZ},\,\,\,j=1,2,3.$ (2.14)

Then, the model (2.1) with $V\neq 0$, defined for some set of special coupling
parameters $\beta_{j}$ satisfying (2.14), will inherit the degenerate vacua
(2.11) of the integrable ATM model with $V=0$.

The model (2.1) with additional mass term and potential $V=\lambda\Phi^{4}$
has been considered in [30] through a numerical method. It was one of the
first attempts to solve the equations self-consistently and non-perturbatively
in which the authors obtained the bound states of the fermion and the shape of
the soliton dynamically. In several treatments of this type of models the
pseudoscalar field was a priori prescribed and the back reaction of the spinor
on the soliton has been neglected. For several aspects of the interaction of
fermions with kinks, kink-antikink and domain walls, see e.g. [32, 33, 34].
The recent studies including the back-reaction were mainly performed using
numerical methods [35, 36]. In a recent work the analytical scattering states
of the fermion on a prescribed sine-Gordon kink has been considered [37]. So,
the analytical treatment of the back-reaction of the fermion bound state on
the kink is lacking in the literature. Here, we consider (2.1) with $V\neq 0$
and provide the first analytical results through the tau function approach of
soliton theory.

We will see below that the model (2.1) with potential $V\neq 0$ would exhibit
scalar field solitons approaching wide range of asymptotic values, apart from
the values $\beta\Phi(\pm\infty)\equiv 2\pi n,\,n\in\leavevmode\hbox{\sf
Z\kern-3.99994ptZ}$ as in (2.11). In fact, one can impose relevant boundary
conditions compatible with the dynamics and the vacuum points defined by the
zeroes of the potential $V$. So, as the potential parameters are continuously
varied, one can thought that the topologically nontrivial configurations of
the scalar field $\Phi$ of the model (2.1) evolve continuously and slowly from
a topologically trivial configuration to soliton configurations (with
$\Phi(\pm\infty)\equiv\pm\phi_{0}$, $\phi_{0}\in\leavevmode\hbox{\rm
I\kern-1.79993ptR}$), such that the corresponding spinor bound states assume
continuously varying configurations. Below, we will obtain analytical kink and
fermion bound state solutions for a finite set of values of $\phi_{0}$ and
resort to numerical simulations to get those solutions for
$\phi_{0}\in\leavevmode\hbox{\rm I\kern-1.79993ptR}$ belonging to some
continuous intervals.

##  3 Topological configurations coupled to the fermion

Our aim is to investigate exact analytical solutions of the system of eqs.
(2.3)-(2.4) such that nontrivial topological configurations can be achieved by
the field $\Phi$, interacting with the bound states formed by the fermion
field $\psi$. In the next steps we rewrite the equations of motion of the
model in terms of the scalar field $\Phi$ and a set of real fermion field
components $\xi_{a}\,\,(a=1,...,4)$. So, let us consider the Ansatz

$\displaystyle\psi_{{}_{E}}(x,t)$ $\displaystyle=$ $\displaystyle
e^{-iEt}\Big{[}\begin{array}[]{c}\xi_{3}(x)+i\,\xi_{4}(x)\\\
\xi_{1}(x)+i\,\xi_{2}(x)\end{array}\Big{]}.$ (3.3)

So, the system (2.3)-(2.4) becomes

$\displaystyle\xi^{{}^{\prime}}_{1}+E\,\xi_{2}-M\xi_{4}\sin{\beta\Phi}+M\xi_{3}\cos{\beta\Phi}$
$\displaystyle=$ $\displaystyle 0,$ (3.4)
$\displaystyle\xi^{{}^{\prime}}_{2}-E\,\xi_{1}+M\xi_{3}\sin{\beta\Phi}+M\xi_{4}\cos{\beta\Phi}$
$\displaystyle=$ $\displaystyle 0,$ (3.5)
$\displaystyle\xi^{{}^{\prime}}_{3}-E\,\xi_{4}+M\xi_{2}\sin{\beta\Phi}+M\xi_{1}\cos{\beta\Phi}$
$\displaystyle=$ $\displaystyle 0,$ (3.6)
$\displaystyle\xi^{{}^{\prime}}_{4}+E\,\xi_{3}-M\xi_{1}\sin{\beta\Phi}+M\xi_{2}\cos{\beta\Phi}$
$\displaystyle=$ $\displaystyle 0,$ (3.7)
$\displaystyle\partial_{t}^{2}\Phi-\partial_{x}^{2}\Phi+2M\beta\Big{[}(\xi_{1}\xi_{3}+\xi_{2}\xi_{4})\cos{\beta\Phi}-(\xi_{1}\xi_{4}-\xi_{2}\xi_{3})\sin{\beta\Phi}\Big{]}+V^{\prime}[\Phi]$
$\displaystyle=$ $\displaystyle 0,$ (3.8)

where $E$ can be regarded as an eigenvalue of the Hamiltonian in the spinor
sector. The self-consistent solutions for a fermion coupled to static scalar
field in the form of a kink (domain wall) have been considered mainly using
numerical methods in the literature.The self-consistent approach to find
fermionic bound states must be contrasted to the results in which the scalar
field kink is considered as an external field or a priori prescribed, as
presented for example in [33] and [34] for of the $\Phi^{4}$ model and for the
case of the kink-antikink of the sine Gordon model, respectively. However, as
we will show below, the presence of the fermion bound state will exert a back
reaction on the kink by changing its profile and the relevant spectra of the
system. Moreover, we will show that the phases of the fermion components will
depend on the parameters of the potential $V$.

We will search for analytical solutions to the of the systems of eqs.
(3.4)-(3.8), such that the static version of (3.8) is considered. In order to
tackle this problem we follow the tau function approach developed for the
related integrable ATM model in [15, 16]. For vanishing potential $V(\Phi)=0$,
the integrable ATM model defined in (2.1) has been studied in [15, 16, 20,
18]. The tau function approach can be adapted conveniently to the non-
integrable system (3.4)-(3.8) provided that the potential $V(\Phi)$ (2.2)
displays the relationships

$\displaystyle\frac{\beta_{j}}{\beta}=\nu_{j},\,\,\,\,\,\,\,\,\,\nu_{j}\in\leavevmode\hbox{\sf
Z\kern-3.99994ptZ},$ (3.9)

between the frequency parameters $\beta_{j}\,(j=1,2,3)$, such that the set
$\\{\nu_{j}\\}$ defines a relatively co-prime integers. This process will
provide a variety of topologically nontrivial configurations for the scalar
field $\Phi$ and the relevant fermionic bound states associated to the field
$\xi_{a}$.

We will look for solutions to the system (3.4)-(3.8) such that they are
invariant under the parity transformation

$\displaystyle{\cal P}_{x}:x\rightarrow-x.$ (3.10)

So, let us consider certain field configurations satisfying

$\displaystyle{\cal P}(\Phi)$ $\displaystyle=$ $\displaystyle-\Phi$ (3.11)
$\displaystyle{\cal P}(V)$ $\displaystyle=$ $\displaystyle V.$ (3.12)

We will pursue for some static solutions to the field $\Phi$ representing a
kink type solution such that

$\displaystyle\Phi(x)$ $\displaystyle=$ $\displaystyle-\Phi(-x),$ (3.13)
$\displaystyle\Phi(\pm\infty)$ $\displaystyle=$
$\displaystyle\pm\phi_{o},\,\,\,\,\Phi(0)=0,\,\,\,\,\phi_{o}=const..$ (3.14)

Under the above assumptions the solutions to the system (3.4)-(3.8) will be
eigenfunctions of the parity operator. Then one must have

$\displaystyle\xi_{1}(x)=-\sigma\,\xi_{4}(-x),\,\,\,\,\xi_{2}(x)=\sigma\,\xi_{3}(-x),\,\,\,\,\,\sigma=\pm
1,$ (3.15)

where $\sigma$ are the eigenvalues of the parity operator ${\cal P}_{x}$. Let
us impose the following boundary conditions to the spinors

$\displaystyle\xi_{a}(\pm\infty)\rightarrow
0,\,\,\,a=1,...,4;\,\,\,\,\xi_{1}(0)=\xi_{0}=const.$ (3.16)

Additional consistency conditions must be satisfied by the solutions; in fact,
from the equations (3.4)-(3.7) one gets the relationships

$\displaystyle\xi_{1}^{2}+\xi_{2}^{2}-\xi_{3}^{2}-\xi_{4}^{2}=c_{1}=const.$
(3.17)

and

$\displaystyle\partial_{x}[\xi_{1}^{2}+\xi_{2}^{2}+\xi_{3}^{2}+\xi_{4}^{2}]+4M\Big{[}(\xi_{1}\xi_{3}+\xi_{2}\xi_{4})\cos{(\beta\Phi)}-(\xi_{1}\xi_{4}-\xi_{2}\xi_{3})\sin{(\beta\Phi)}\Big{]}=0.$
(3.18)

Since all the spinor components vanish at $x\rightarrow\pm\infty$, according
to the boundary condition (3.16), from now on we will set $c_{1}=0$ in (3.17).
In addition to the above set of conditions we require a finite energy and
localized energy density associated to the field configurations.

Useful identities related to (3.18), which do not involve the eingenvalue $E$,
can be obtained by combining some subset of equations of the system
(3.4)-(3.7). So, multiplying by $\xi_{1}$ (3.4) and (3.5) by $\xi_{2}$ and
adding the two expressions one has

$\displaystyle\frac{1}{2}\partial_{x}[\xi_{1}^{2}+\xi_{2}^{2}]+M\left\\{[\xi_{1}\xi_{3}+\xi_{2}\xi_{4}]\cos{(\beta\Phi)}-[\xi_{1}\xi_{4}-\xi_{2}\xi_{3}]\sin{(\beta\Phi)}\right\\}=0.$
(3.19)

Similarly, from (3.6)-(3.7) one can get the next identity

$\displaystyle\frac{1}{2}\partial_{x}[\xi_{3}^{2}+\xi_{4}^{2}]+M\left\\{[\xi_{1}\xi_{3}(x)+\xi_{2}\xi_{4}(x)]\cos{(\beta\Phi)}-[\xi_{1}\xi_{4}(x)-\xi_{2}\xi_{3}(x)]\sin{(\beta\Phi)}\right\\}=0.$
(3.20)

In fact, the identities (3.19) and (3.20) are related to each other due to the
parity symmetries (3.13) and (3.15). Moreover, adding the last eqs. (3.19) and
(3.20) one can get the relationship (3.18).

In addition, the system of equations (3.4)-(3.8) exhibits the following
discrete symmetry

$\displaystyle\Phi$ $\displaystyle\rightarrow$
$\displaystyle-\Phi,\,\,\,E\rightarrow-E$ (3.21) $\displaystyle\xi_{1}$
$\displaystyle\rightarrow$ $\displaystyle-i\xi_{4},\,\,\,\,\xi_{2}\rightarrow
i\xi_{3}$ (3.22) $\displaystyle\xi_{4}$ $\displaystyle\rightarrow$
$\displaystyle\,\,\,\,i\xi_{1},\,\,\,\,\xi_{3}\rightarrow-i\xi_{2}.$ (3.23)

This symmetry implies that one can construct two sets of solutions related to
each other. In fact, once we know the set $\\{\Phi,\psi_{E}\\}$ as a solution
of the system, with $\psi_{E}$ defined in (3.3), one can construct another
set, $\\{-\Phi,\psi_{-E}\\}$, provided that the new components of the spinor
$\psi_{-E}$ are constructed according to (3.22)-(3.23). Below we will discuss
this symmetry in the context of the charge conjugation symmetry and the
related formulation for the Majorana fermions in a system of this type.

Some comments are in order here. First, the static version of (3.8) can be
written as a suggestive equation involving some current densities. So, taking
into account the identity (3.18) into the equation (3.8) one can write

$\displaystyle\partial_{x}\Big{[}J^{0}+\frac{2}{\beta}\partial_{x}\Phi\Big{]}$
$\displaystyle=$ $\displaystyle\frac{2}{\beta}V^{\prime}(\Phi)$ (3.24)
$\displaystyle J^{0}$ $\displaystyle\equiv$
$\displaystyle\xi_{1}^{2}+\xi_{2}^{2}+\xi_{3}^{2}+\xi_{4}^{2}.$ (3.25)

Moreover, by integrating (3.24) in the whole line and taking into account the
relevant b.c.’s one gets

$\displaystyle\int_{-\infty}^{+\infty}dx\,V^{\prime}(\Phi)=0.$ (3.26)

Second, taking into account the parity property (3.13) and the eq. (3.24) one
can show that $J^{0}(x)$ is an even function, i.e. $J^{0}(-x)=J^{0}(x)$.

Third, the equation (3.24) is invariant under the parity symmetry
(3.11)-(3.15) provided that $V$ is an even function, i.e. $V(-\Phi)=V(\Phi)$.
As a realization of this symmetry we will show below the transformation of the
kink into an anti-kink, and the related changes in the signs of the density
charges of the fermion. Moreover, this symmetry implies certain relationship
between the topological charge of the kink (anti-kink) and the fermion (anti-
fermion) bound state Noether charge in the model.

Fourth, as mentioned above, for vanishing potential one has the integrable ATM
model [16, 17]. It has been shown that in this case the model exhibits the
equivalence of the both topological and Noether currents at the classical and
quantum levels, respectively. So, in order to write an analogous equivalence
equation one integrates the b.h.s. of (3.24) in the interval $[-\infty,x]$ in
order to get

$\displaystyle J^{0}(x)$ $\displaystyle=$
$\displaystyle-\frac{2}{\beta}\partial_{x}\Phi+\frac{2}{\beta}\partial^{-1}[\frac{dV(\Phi)}{d\Phi}].$
(3.27)

This defines a first order integro-differential equation relating the Noether
(2.6), topological (2.12) and a new non-local charge densities, respectively

$\displaystyle J^{0}(x)$ $\displaystyle=$
$\displaystyle(-\frac{4\pi}{\beta^{2}})\,j^{0}_{top}+{\cal J}^{0}_{nonloc}$
(3.28) $\displaystyle{\cal J}^{0}_{nonloc}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{\beta}\int^{x}_{-\infty}\frac{dV(\Phi)}{d\Phi}dx^{\prime}-\frac{1}{\beta}\int^{+\infty}_{x}\frac{dV(\Phi)}{d\Phi}dx^{\prime}.$
(3.29)

Fifth, remarkably, one can show that the system of first order equations
comprising the eqs. (3.4)-(3.7) and the first order integro-differential eq.
(3.27) imply the static version of the second order differential eq. (3.8).
So, a solution of the first-order system of eqs. (3.4)-(3.7) and (3.27) will
be solution of the static version of the second order system of eqs.
(3.4)-(3.8).

Sixth, the strength of the back-reaction of the fermion on the kink can be
controled by the coupling constant $\beta$. In fact, by taking the limit
$\beta\rightarrow 0$ of the integro-differential equation in (3.27) and upon
taking the $x-$derivative of the remaining eq. one is left with an effective
equation for the scalar field $\Phi$ alone. In fact, in this limit the fermion
decouples from the scalar as it is clear from the Lagrangian (2.1).

Next, we write some identities which will be useful in the construction of the
analytic solutions and in order to write a second order algebraic equation for
the eigenvalue $E$. A useful relationship arises for the charge density
$J^{0}$ once the equation (3.18) is integrated in $x$. In fact, taking into
account (3.25) one can rewrite (3.18) as

$\displaystyle
J^{0}=-4M\int_{-\infty}^{x}\,dx\,\Big{[}(\xi_{1}\xi_{3}+\xi_{2}\xi_{4})\cos{(\beta\Phi)}-(\xi_{1}\xi_{4}-\xi_{2}\xi_{3})\sin{(\beta\Phi)}\Big{]}.$
(3.30)

Taking the derivatives of the system (3.4)-(3.7) one can write the identity

$\displaystyle(E^{2}-M^{2})[\xi^{2}_{1}+\xi^{2}_{2}+\xi^{2}_{3}+\xi^{2}_{4}]$
$\displaystyle=$
$\displaystyle-[\xi_{1}\xi^{\prime\prime}_{1}+\xi_{2}\xi^{\prime\prime}_{2}+\xi_{3}\xi^{\prime\prime}_{3}+\xi_{4}\xi^{\prime\prime}_{4}]+$
(3.31) $\displaystyle
2M\beta\Phi^{\prime}\\{[\xi_{1}\xi_{4}-\xi_{2}\xi_{3}]\cos{\beta\Phi}+[\xi_{1}\xi_{3}+\xi_{2}\xi_{4}]\sin{\beta\Phi}\\}$

Let us consider the normalization

$\displaystyle\int_{-\infty}^{+\infty}dx\,[\xi^{2}_{1}+\xi^{2}_{2}+\xi^{2}_{3}+\xi^{2}_{4}]=1.$
(3.32)

Integrating in the whole line the expression (3.31) and taking into account
the normalization (3.32), as well as the relevant boundary conditions one can
write

$\displaystyle(E^{2}-M^{2})$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}dx\,[(\xi^{\prime}_{1})^{2}+(\xi^{\prime}_{2})^{2}+(\xi^{\prime}_{3})^{2}+(\xi^{\prime}_{4})^{2}]-$
(3.33) $\displaystyle
2M\beta\int_{-\infty}^{\infty}dx\,\Phi\frac{d}{dx}\\{[\xi_{1}\xi_{4}-\xi_{2}\xi_{3}]\cos{\beta\Phi}+[\xi_{1}\xi_{3}+\xi_{2}\xi_{4}]\sin{\beta\Phi}\\}.$

In addition, using the eqs. of motion one can write the identity

$\displaystyle\frac{d}{dx}\\{[\xi_{1}\xi_{4}-\xi_{2}\xi_{3}]\cos{\beta\Phi}+[\xi_{1}\xi_{3}+\xi_{2}\xi_{4}]\sin{\beta\Phi}\\}$
$\displaystyle=$ $\displaystyle-[2E-\beta\Phi^{\prime}]\times$
$\displaystyle\Big{[}(\xi_{1}\xi_{3}+\xi_{2}\xi_{4})\cos{(\beta\Phi)}-(\xi_{1}\xi_{4}-\xi_{2}\xi_{3})\sin{(\beta\Phi)}\Big{]}.$
$\displaystyle=$
$\displaystyle\frac{1}{4M}[2E-\beta\Phi^{\prime}]J^{\prime}_{0},$ (3.35)

where in order to get the last line we have used (3.18) and the expression of
$J^{0}$ in (3.25). Then, taking into account the last relationship one can
write (3.33) as

$\displaystyle E^{2}$ $\displaystyle=$ $\displaystyle
M^{2}+\int_{-\infty}^{\infty}dx\,\Big{\\{}(\xi^{\prime}_{1})^{2}+(\xi^{\prime}_{2})^{2}+(\xi^{\prime}_{3})^{2}+(\xi^{\prime}_{4})^{2}\Big{\\}}-\beta
E\int_{-\infty}^{\infty}dx\,\Phi
J^{\prime}_{0}+\frac{\beta^{2}}{2}\int_{-\infty}^{\infty}dx\,\Phi\Phi^{\prime}J^{\prime}_{0}$
(3.36)

This equation provides a second order algebraic equation for the eingenvalue
$E$. We will see below that analogous second order algebraic equation for the
eingenvalue $E$ emerges in the context of the tau function formalism applied
to solve the system (3.4)-(3.8) in order to find kink type and bound state
solutions.

###  3.1 Threshold states

The threshold or half-bound states are defined as states where the fermion
field approaches a constant value at spatial infinity. The solutions $\xi_{a}$
of this type of state are finite but they do not decay fast enough at infinity
to be square integrable [38, 39, 40]. In order to examine these type of states
we will write an equation involving only the bilinears of the
$\xi_{a}\,^{\prime}s$, their derivatives and the eigenvalue $E$. In fact, from
the system (3.4)-(3.7) and taking into account the trigonometric identity
$\sin^{2}{\beta\Phi}+\cos^{2}{\beta\Phi}=1$ one can write the following
relationship

$\displaystyle{\cal A}_{2}(E^{2}-M^{2})+2{\cal A}_{1}E+{\cal A}_{0}=0,$ (3.37)

with

$\displaystyle{\cal A}_{2}$ $\displaystyle\equiv$
$\displaystyle\xi_{1}^{2}+\xi_{2}^{2}+\xi_{3}^{2}+\xi_{4}^{2},$ (3.38)
$\displaystyle{\cal A}_{1}$ $\displaystyle\equiv$
$\displaystyle\xi^{\prime}_{1}\xi_{2}-\xi_{1}\xi^{\prime}_{2}-\xi^{\prime}_{3}\xi_{4}+\xi_{3}\xi^{\prime}_{4},$
(3.39) $\displaystyle{\cal A}_{0}$ $\displaystyle\equiv$
$\displaystyle(\xi^{\prime}_{1})^{2}+(\xi^{\prime}_{2})^{2}+(\xi^{\prime}_{3})^{2}+(\xi^{\prime}_{4})^{2}.$
(3.40)

So, one can analyze (3.37) for these type of states in order to inspect the
values of $E$ between the bound and continuum energy spectrum. So, let us
consider

$\displaystyle\xi_{a}(\pm\infty)=c^{\pm}_{a},\,\,\,\,\xi^{\prime}_{a}(\pm\infty)=0,\,\,\,\,\,\,c_{a}^{\pm}=\mbox{const.},\,\,\,a=1,2,3,4.$
(3.41)

Applying the condition (3.41) one has that ${\cal A}_{0}^{\pm}={\cal
A}_{1}^{\pm}=0$ in (3.39)-(3.40) and then the equation (3.37) becomes

$\displaystyle{\cal A}_{2}^{\pm}[E^{2}_{thr}-M^{2}]=0.$ (3.42)

Since by definition one has ${\cal A}_{2}^{\pm}\neq 0$, the algebraic equation
(3.42) provides the energies of the threshold states

$\displaystyle E_{thr}^{\pm}=\pm M.$ (3.43)

###  3.2 Bound states in the continuum (BIC)

The equation (3.36) can be used to analyze the behavior of the spectra around
the threshold states, $E=E_{th}^{\pm}$ as defined in (3.43). In fact, for kink
type solitons in the limit of small soliton width one can consider the
function $\Phi=\frac{\phi_{o}}{2}(H(x)-H(-x)),\,(\phi_{o}>0)$, where $H(x)$ is
the Heaviside step function. Then, substituting this form of kink for $\Phi$
into the equation (3.36) one can write the next second order algebraic
equation for $E$

$\displaystyle E^{2}-\lambda_{o}\,E-(M^{2}+I)=0,$ (3.44)

where $\lambda_{o}\equiv\phi_{o}J_{0}(0)$ and
$I\equiv\int_{-\infty}^{\infty}dx\,[(\xi^{\prime}_{1})^{2}+(\xi^{\prime}_{2})^{2}+(\xi^{\prime}_{3})^{2}+(\xi^{\prime}_{4})^{2}]$.
Therefore, solving the second order equation (3.44) for the unknown $E$ one
has

$\displaystyle E$ $\displaystyle=$ $\displaystyle\pm
M\Delta+\frac{1}{2}\lambda_{o},\,\,\,\,\,\,\,\,\lambda_{o}\equiv\phi_{o}J_{0}(0)$
(3.45) $\displaystyle\Delta$ $\displaystyle\equiv$
$\displaystyle\Big{\\{}1+\frac{\lambda_{0}^{2}+4I}{4M^{2}}\Big{\\}}^{1/2},\,\,\,\,\lambda_{o}>0,\,\,\,\,I>0.$
(3.46)

Notice that, since $I>0$ then $\Delta>1$. Therefore, from the equation (3.45)
one can see the appearance of bound states in the continuum (BIC), i.e. states
with $E>E^{+}_{th}$, with $E^{+}_{th}$ being the positive threshold state in
(3.43) . In fact, one has a state which takes a value
$E=+M\Delta+\frac{1}{2}\lambda_{0}$ (positive BIC) which is located above the
positive threshold $E^{+}_{th}=+M$ in (3.43). Similarly, one has a state with
energy $E=-M\Delta+\frac{1}{2}\lambda_{0}$. This will be an in-gap state for
$E>E^{-}_{th}$ provided that $M>\frac{I}{\lambda_{o}}$.

Moreover, in order to examine the negative BIC states one must rederive an
analogous equation to (3.44) specialized for an anti-kink soliton. So, from
(3.36) for anti-kink $\Phi=-\frac{\phi_{o}}{2}(H(x)-H(-x))$ one has

$\displaystyle E^{2}+\lambda_{o}\,E-(M^{2}+I)=0.$ (3.47)

Therefore, the solutions become

$\displaystyle E$ $\displaystyle=$ $\displaystyle\pm
M\Delta-\frac{1}{2}\lambda_{o},$ (3.48)

with $\Delta$ and $\lambda_{o}$ defined in (3.45)-(3.46). Then, one has a
negative BIC state with energy value $E=-M\Delta-\frac{1}{2}\lambda_{0}$,
which is located below the negative threshold $E^{-}_{th}=-M$ in (3.43). In
addition, one has a state with energy $E=+M\Delta-\frac{1}{2}\lambda_{0}$.
This is an in-gap state ($E<+M$) provided that $M>\frac{I}{\lambda_{o}}$.
Below we will find some positive and negative BIC bound states in the context
of the tau function approach.

###  3.3 Majorana zero-mode states

The system of equations (3.4)-(3.7) can be written in matrix form as

$\displaystyle H\xi$ $\displaystyle=$ $\displaystyle E\,\xi,$ (3.49)
$\displaystyle\xi$ $\displaystyle\equiv$
$\displaystyle\left(\begin{array}[]{c}\xi_{R}\\\
\xi_{L}\end{array}\right)=\left(\begin{array}[]{c}\xi_{3}+i\xi_{4}\\\
\xi_{1}+i\xi_{2}\end{array}\right),$ (3.54)

or in components

$\displaystyle\left(\begin{array}[]{cc}i\frac{d}{dx}&iMe^{-i\beta\Phi}\\\
-iMe^{i\beta\Phi}&-i\frac{d}{dx}\end{array}\right)\left(\begin{array}[]{c}\xi_{R}\\\
\xi_{L}\end{array}\right)$ $\displaystyle=$ $\displaystyle
E\left(\begin{array}[]{c}\xi_{R}\\\ \xi_{L}\end{array}\right).$ (3.61)

We shall assume the operator $H$ to be a hermitian operator with eigenfuntions
$\xi$ and a spectrum of real eingenvalues. One can verify that $H$ satisfies

$\displaystyle\Gamma^{-1}H\Gamma=H^{\star},\,\,\,\,\,\,\Gamma\equiv\pm
i\gamma_{1}.$ (3.62)

So, taking into account (3.62) and the complex conjugate of (3.49) one can
write the identity

$\displaystyle\xi=\Gamma\xi^{\star}.$ (3.63)

In terms of the Dirac field $\psi_{{}_{E}}$ defined in (3.3), and taking into
account $\psi_{{}_{-E}}^{\star}\equiv e^{iEt}\xi^{\star}$, the last
relationship can be written as

$\displaystyle\psi_{{}_{E}}=\Gamma\psi_{{}_{-E}}^{\star}.$ (3.64)

This relationship exhibits the charge conjugation symmetry of the system and
then, one expects that the particles and holes will have identical spectra.
So, one can formulate Majorana fermions for a system of this kind since there
exists a particle-hole symmetry, or, in a relativistic model, a charge
conjugation symmetry.

A Majorana fermion is defined by considering the particle and hole with the
same energy as a single particle. Then, from the relationship (3.64) and for
zero-mode solutions ($E=0$) one can write

$\displaystyle\psi_{{}_{0}}=\Gamma\psi_{{}_{0}}^{\star}.$ (3.65)

In fact, the fermions obeying the reality condition (3.65) are defined as
Majorana fermions. These states are well known to lead to interesting
phenomena in condensed matter and we will search for them as fermion bound
states. Recently, the model (3.61) has been considered as the continuum limit
of a one-dimensional system of spin-orbit coupled Dirac Hamiltonian and a
s-wave superconducting pairing, such that the s-wave pairing is represented by
$-iMe^{i\beta\Phi}$ [41]. The bosonization techniques have been applied to
related fermionic systems with finite length which exhibit Majorana zero modes
[42]. So, the complementary descriptions provided by the fermionic and bosonic
formulations of the superconducting phase have been considered. Below, in the
context of the tau-function approach, we will perform analogous Majorana zero
mode-soliton duality description of the model (3.61) such that the static
scalar field is a self-consistent solution of (3.8).

So, setting $E=0$ into the system of eqs. (3.4)-(3.7) one has that the zero-
mode states $\xi_{a}$ must satisfy

$\displaystyle(\xi^{\prime}_{1})^{2}+(\xi^{\prime}_{2})^{2}+(\xi^{\prime}_{3})^{2}+(\xi^{\prime}_{4})^{2}-M^{2}[(\xi_{1})^{2}+(\xi_{2})^{2}+(\xi_{3})^{2}+(\xi_{4})^{2}]=0.$
(3.66)

Taking into account the normalization (3.32) one can rewrite (3.66) as

$\displaystyle\int_{-\infty}^{+\infty}dx\,[(\xi^{\prime}_{1})^{2}+(\xi^{\prime}_{2})^{2}+(\xi^{\prime}_{3})^{2}+(\xi^{\prime}_{4})^{2}]=M^{2}.$
(3.67)

Next, imposing the reality condition for $E=0$ (3.65) (or in the form 3.63) on
the spinor $\xi$ one has the following two cases:
Case I ($\Gamma=-i\gamma_{1}$)

$\displaystyle\xi_{1}=-\xi_{3},\,\,\,\,\xi_{2}=\xi_{4}.$ (3.68)

Case II ($\Gamma=i\gamma_{1}$)

$\displaystyle\,\,\,\,\,\,\,\,\,\,\xi_{1}=\xi_{3},\,\,\,\,\,\,\,\xi_{2}=-\xi_{4}.$
(3.69)

Below, in the tau function approach, we will find zero-mode bound states
satisfying the above conditions (3.67) and (3.68)-(3.69).

###  3.4 A special potential $V(\Phi)$ with
$\frac{\beta_{j}}{\beta}=j;\,j=1,2,3$

We will consider a particular potential (2.2) with the next relationships
imposed to the rational fractions between the coupling constants $\beta_{j}$
in (3.9)

$\displaystyle\beta_{1}$ $\displaystyle=$
$\displaystyle\beta,\,\,\,\beta_{2}=2\beta,\,\,\,\beta_{3}=3\beta.$ (3.70)

Notice that a particular set for the relatively prime integers
$\nu_{1}=1,\nu_{2}=2$ and $\nu_{3}=3$ has been chosen in (3.9). So, the
potential has the extrema

$\displaystyle\Phi_{c}=\frac{\pi}{\beta}m,\,\,m\in\leavevmode\hbox{\sf
Z\kern-3.99994ptZ};\,\,\,\,\,\,\,\,\pm\frac{\theta_{0}}{\beta}+\frac{2\pi}{\beta}n,\,\,\,\,n\in\leavevmode\hbox{\sf
Z\kern-3.99994ptZ},$ (3.71)

where

$\displaystyle\theta_{0}\equiv\frac{1}{\beta}\arccos{\Big{[}\frac{1}{6}\left(-\frac{A_{2}}{A_{3}}+\epsilon_{o}\sqrt{(\frac{A_{2}}{A_{3}})^{2}-3\frac{A_{1}}{A_{3}}+9}\right)}\Big{]};\,\,\,\,\,\,\,\,\,\,\epsilon_{o}\equiv\pm
1,\,\,\,\,\theta_{0}\in[-\pi,+\pi].$ (3.72)

The second extrema in (3.71) is defined for
$A_{2}^{2}-3A_{1}A_{3}+9A_{3}^{2}\geq 0$. From (3.72) one can get the next
condition for the parameters

$\displaystyle
A_{1}=-4A_{2}\cos{(\beta\theta_{0})}-3[1+2\cos{(2\beta\theta_{0})}]A_{3},\,\,\,\,\,\,\,\,\,\,\,A_{3}\neq
0.$ (3.73)

Consider the case $A_{3}\neq 0$ and define

$\displaystyle
z\equiv\frac{1}{6}\frac{A_{2}}{A_{3}}+\epsilon_{o}\sqrt{(\frac{A_{2}}{A_{3}})^{2}-3\frac{A_{1}}{A_{3}}+9},\,\,\,\,\,\epsilon_{o}\equiv\pm
1.$ (3.74)

So, one can parametrize the quantities $(A_{1}/A_{3})$ and $(A_{2}/A_{3})$ as

$\displaystyle\frac{A_{1}}{A_{3}}=3-12z\cos{(\beta\theta_{0})},\,\,\,\,\frac{A_{2}}{A_{3}}=3z-3\cos{(\beta\theta_{0})}.$
(3.75)

It is useful to analyze separately the cases $|z|>1$ and $|z|\leq 1$. So, for
the case $|z|>1$ and taking into account the above parametrization the
potential becomes

$\displaystyle V_{1}(\Phi)$ $\displaystyle=$ $\displaystyle
A_{3}\Big{\\{}[3-12z\cos{(\beta\theta_{0})}]\cos{(\beta\Phi)}+[3z-3\cos{(\beta\theta_{0})}]\cos{(2\beta\Phi)}+\cos{(3\beta\Phi)}\Big{\\}}-$
$\displaystyle\Big{\\{}-\frac{1}{3}A_{3}\Big{[}12z-3\cos{(\beta\theta_{0})}+6z\cos{(\beta\theta_{0})}+\cos{(3\beta\theta_{0})}\Big{]}\Big{\\}},$
$\displaystyle|z|>1,\,\,\,\mbox{sign}(A_{3})=\mbox{sign}(z),$

where the constant expression has been subtracted in order to define the zero
level of the potential at the vacua $\pm\frac{\theta_{0}}{\beta}$. So, (3.4)
describes a family of potentials with kink-type solutions interpolating the
vacua $\pm\frac{\theta_{0}}{\beta}$. The parameters $\\{z,A_{3}\\}$ must be
chosen from either the set $\\{z>1,\,A_{3}>0\\}$ or the set
$\\{z<-1,\,A_{3}<0\\}$, for each value of $\theta_{0}$ in the interval
$\theta_{0}\in[-\pi,+\pi]$ .

For the particular values of $\theta_{0}=\pm\frac{\pi}{3},\,\pm\frac{2\pi}{3}$
we show the potentials in the Figs. 1. These special values for $\theta_{0}$,
as well as the relationship (3.73), will be obtained below for kink-type and
fermion bound state solutions in the tau function approach. Notice that the
values $\Phi_{c}=(3n\theta_{0}),\,n\in\leavevmode\hbox{\sf
Z\kern-3.99994ptZ}$, for these particular $\theta_{0}^{\prime}s$ correspond to
the first set of the extrema in (3.71) and in this case represents a false
vacua, as shown in the Figs. 1.

For arbitrary real values of $\theta_{0}$ in the interval
$\theta_{0}\in[-\pi,+\pi]$ we will use a numerical simulation in order to find
the relevant spectra of kinks and fermion bound states.

Next, let us define the potential in the region $|z|\leq 1$ such that some
points of the extrema in (3.71) can be parametrized by
$\frac{1}{\beta}\arccos{(-z)}$. So, consider the potential

$\displaystyle V_{2}(\Phi)$ $\displaystyle=$ $\displaystyle
A_{3}\Big{\\{}[3-12z\cos{(\beta\theta_{0})}]\cos{(\beta\Phi)}+[3z-3\cos{(\beta\theta_{0})}]\cos{(2\beta\Phi)}+\cos{(3\beta\Phi)}\Big{\\}}-c_{0},\,\,\,\,\,\,|z|\leq
1,$

where a constant value $c_{0}$ has been subtracted in order to define the zero
level of the potential at the true vacua. The extrema in this parametrization
become

$\displaystyle\Big{\\{}\frac{\pi k}{\beta},\,\,\,\,k\in\leavevmode\hbox{\sf
Z\kern-3.99994ptZ}\Big{\\}};\,\,\,\,\,\Big{\\{}\pm\frac{\theta_{0}}{\beta}+\frac{2\pi
n}{\beta},\,\,\,\,n\in\leavevmode\hbox{\sf
Z\kern-3.99994ptZ}\Big{\\}};\,\,\,\,\,\,\,\,\Big{\\{}\pm\frac{\alpha_{0}}{\beta}+\frac{2\pi
m}{\beta},\,\,\,\,m\in\leavevmode\hbox{\sf Z\kern-3.99994ptZ}\Big{\\}}.$
(3.77)

with

$\displaystyle\alpha_{0}\equiv\arccos{(-z)}.$ (3.78)

Some points of the extrema in (3.77) become a true vacua while the other
subset of points become a false vacua. Let us consider the potential (3.4)
possessing the true vacua at $\pm\frac{\alpha_{0}}{\beta}$

$\displaystyle V_{2}^{(1)}(\Phi)$ $\displaystyle=$ $\displaystyle
A_{3}\Big{\\{}[3-12z\cos{(\beta\theta_{0})}]\cos{(\beta\Phi)}+[3z-3\cos{(\beta\theta_{0})}]\cos{(2\beta\Phi)}+\cos{(3\beta\Phi)}\Big{\\}}-$
$\displaystyle
A_{3}\Big{\\{}-3z+2z^{3}+(3+6z^{2})\cos{(\beta\theta_{0})}\Big{\\}},\,\,\,\,\,\,\,|z|\leq
1,$

with

$\displaystyle\\{A_{3}<0,\,z>-\cos{\theta_{0}}\\}\,\,\,\,\,\mbox{or}\,\,\,\,\,\,\\{A_{3}>0,\,z<-\cos{\theta_{0}}\\}$
(3.80)

such that the constant $c_{0}$ in (3.4) has been assumed to take a particular
form. In the Fig. 2 (left) we show an example such that the first two sets of
(3.77) become a false vacua while the remaining subset becomes a true vacua
with minima $\pm\alpha_{0}+2\pi m$.

Similarly, the potential with vacua as the first set of points in (3.77) for
even $k$ can be written as

$\displaystyle
V_{2}^{(2)}(\Phi)=A_{3}\Big{\\{}[3-12z\cos{(\beta\theta_{0})}]\cos{(\beta\Phi)}+[3z-3\cos{(\beta\theta_{0})}]\cos{(2\beta\Phi)}+\cos{(3\beta\Phi)}\Big{\\}}-$
$\displaystyle
A_{3}\Big{\\{}4+3z-3(1+4z)\cos{(\beta\theta_{0})}\Big{\\}},\,\,\,\,\,\,|z|\leq
1,$ (3.81)

such that

$\displaystyle A_{3}<0,\,\,\,\,z>\cos{\theta_{0}}-2,$ (3.82)

where a particular form of $c_{0}$ has been set to define the zero point
energy. In the Fig. 2 (right) we show another example such that the last two
sets of (3.77) become a false vacua while the first subset becomes a true
vacua with minima $\frac{2n\pi}{\beta},\,n\in\leavevmode\hbox{\sf
Z\kern-3.99994ptZ}$.

Finally, the potential with vacua as the first set of points in (3.77) for odd
$k$ can be written as

$\displaystyle
V_{2}^{(3)}(\Phi)=A_{3}\Big{\\{}[3-12z\cos{(\beta\theta_{0})}]\cos{(\beta\Phi)}+[3z-3\cos{(\beta\theta_{0})}]\cos{(2\beta\Phi)}+\cos{(3\beta\Phi)}\Big{\\}}-$
$\displaystyle
A_{3}\Big{\\{}-4+3z-3(1-4z)\cos{(\beta\theta_{0})}\Big{\\}},\,\,\,\,\,\,|z|\leq
1,$ (3.83)

such that

$\displaystyle A_{3}>0,\,\,\,\,z<\frac{1}{3}(2-\cos{\theta_{0}}),$ (3.84)

where a convenient $c_{0}$ has been fixed in order to define the zero point
energy. In the Fig. 3 we show an example such that the last two sets of (3.77)
become a false vacua while the first subset becomes a true vacua with minima
$\frac{(2n+1)\pi}{\beta},\,n\in\leavevmode\hbox{\sf Z\kern-3.99994ptZ}$.

We will concentrate below, analytically (secs. 4 and 5) and numerically (see
sec. 7), to the various types of fermion bound states induced by the kink-type
solutions associated to the degenerate vacua described by the potential (3.4).
For another set of parameters such as $A_{1}\neq
0,\,\beta_{1}=\beta,\,A_{2}=A_{3}=0$ in (2.2) we will show that Majorana bound
states are hosted by soliton type solutions interpolating the first type of
vacua in (3.71) for $m=\pm 1$ (see sections 6 and 6.2).

Figure 1: The plots of the potential $V_{1}(\Phi)$ in (3.4) for
$A_{3}=0.1,z=4.5,\,\beta=1$. Left Fig. for $\theta_{0}=\frac{\pi}{3}$ and
vacua $(2n+1)\theta_{o},\,n\in\leavevmode\hbox{\sf Z\kern-3.99994ptZ}$. Right
Fig. for $\theta_{0}=\frac{2\pi}{3}$ and vacua
$n\theta_{o},\,n\in\leavevmode\hbox{\sf Z\kern-3.99994ptZ}$. Figure 2: Left,
the potential $V_{2}^{(1)}(\Phi)$ in (3.4) for
$A_{3}=-0.1,z=0.2,\,\beta=1,\,\theta_{0}=0.8$ and vacua $\alpha_{o}=\pm 1.77$.

Figure 3: Left, the potential $V_{2}^{(2)}(\Phi)$ in (3.4) for
$A_{3}=-0.1,z=0.3,\,\beta=1,\,\theta_{0}=1.15,\,\alpha_{o}=1.88$ and vacua
$\pm 2\pi n,\,n\in\leavevmode\hbox{\sf Z\kern-3.99994ptZ}$. Right, the
potential $V_{2}^{(3)}(\Phi)$ in (3.4) for
$A_{3}=0.1,z=0.3,\,\beta=1,\,\theta_{0}=0.8,\,\alpha_{o}=1.88$ and vacua
$\pm(2n+1)\pi,\,n\in\leavevmode\hbox{\sf Z\kern-3.99994ptZ}$.

##  4 Soliton solutions and the tau function approach

In this section we solve the system of equations (3.4)-(3.8) analytically and
determine a class of kink and fermion bound states with eigenvalue $E$. In
order to perform this task we will follow the tau function approach and write
the fields $\xi_{a}$ and $\Phi$ in terms of the so-called tau functions. Even
though the tau function and Hirota methods have been applied to integrable
systems, those methods can also be applied to some non-integrable models in
order to find 1-soliton type solutions as in the generalized sine-Gordon
models [19] and, sometimes, in order to construct the two-soliton solutions,
as in the non-integrable modified regularized long wave equation (mRLW) [43,
44].

So, let us consider the scalar and the fermion fields of the model to depend
on the tau functions
$\\{\tau_{1},\,\tau_{0},\,\tau_{\xi}^{a},\,\widetilde{\tau}_{\xi}^{a}\\}$
through the Ansatz

$\displaystyle e^{i\beta\Phi}$ $\displaystyle=$ $\displaystyle
e^{-i\theta_{0}}\,\frac{\tau_{1}}{\tau_{0}},\,\,\,\,$ (4.1)
$\displaystyle\xi_{1}$ $\displaystyle=$
$\displaystyle\,[\frac{\tau_{\xi}^{1}}{\tau_{0}}]\,e^{-i\theta_{1}}+\,[\frac{\widetilde{\tau}_{\xi}^{1}}{\tau_{1}}]\,e^{i\theta_{1}}$
(4.2) $\displaystyle\xi_{2}$ $\displaystyle=$
$\displaystyle-i\\{\,[\frac{\tau_{\xi}^{2}}{\tau_{0}}]\,e^{-i\theta_{2}}-\,[\frac{\widetilde{\tau}_{\xi}^{2}}{\tau_{1}}]\,e^{i\theta_{2}}\\}$
(4.3) $\displaystyle\xi_{3}$ $\displaystyle=$
$\displaystyle\,[\frac{\tau_{\xi}^{3}}{\tau_{0}}]\,e^{-i\theta_{3}}+\,[\frac{\widetilde{\tau}_{\xi}^{3}}{\tau_{1}}]\,e^{i\theta_{3}}$
(4.4) $\displaystyle\xi_{4}$ $\displaystyle=$ $\displaystyle
i\\{\,[\frac{\tau_{\xi}^{4}}{\tau_{0}}]\,e^{-i\theta_{4}}-\,[\frac{\widetilde{\tau}_{\xi}^{4}}{\tau_{1}}]\,e^{i\theta_{4}}\\},$
(4.5)

where $\theta_{0},\,\theta_{a}\,(a=1,...,4)$ are real parameters.

We will concentrate below on the analytical kink and fermionic bound state
profiles for the case $\beta_{n}=n\beta\,(n=1,2,3)$ in (3.70) and the
topological sector defined by the charge $Q_{topol}=\frac{\theta_{0}}{\pi}$.

Next, we substitute the Ansatz (4.1)-(4.5) into the system (3.4)-(3.8) for the
fermion components coupled to the field $\Phi$ and write the differential
equations in terms of the relevant tau functions. The four eqs. (3.4)-(3.7) in
polynomial form in terms of the tau functions are provided in the appendix A,
see eqs.(A.1)-(A.4), respectively.

Notice that for the chosen particular potential with the coupling constants
defined in (3.70) the derivative of the potential $V^{\prime}(\Phi)$ terms in
the eq. (3.8) can be written as integer powers of the tau functions. So, in
the appendix B we present the scalar field equation (3.8) written in the form
of a polynomial equation in integer powers of the tau functions (see eq.
(B.1)).

In addition, the consistency conditions (3.17) (for $c_{1}=0$) and (3.18) in
terms of the tau functions are provided in the appendix C, see eqs.
(C.1)-(C.2), respectively.

###  4.1 Tau functions for kink-type solitons and fermion bound states

A self-consistent solution to the system of equations (3.4)-(3.8) for the
topological sectors defined by the parameter $\theta_{0}$ will be found for
the following tau functions

$\displaystyle\tau_{0}$ $\displaystyle=$ $\displaystyle
1+e^{-i\theta_{0}}e^{2\kappa x},\,\,\,\tau_{1}=1+e^{i\theta_{0}}e^{2\kappa
x},$ (4.6) $\displaystyle\tau_{\xi}^{1}$ $\displaystyle=$
$\displaystyle\rho_{1}e^{\kappa x},\,\,\tau_{\xi}^{2}=\rho_{2}e^{\kappa
x},\,\,\tau_{\xi}^{3}=\rho_{3}e^{\kappa
x},\,\,\tau_{\xi}^{4}=\rho_{4}e^{\kappa x},$ (4.7)
$\displaystyle\widetilde{\tau}_{\xi}^{1}$ $\displaystyle=$
$\displaystyle\rho_{1}e^{\kappa
x},\,\,\widetilde{\tau}_{\xi}^{2}=\rho_{2}e^{\kappa
x},\,\,\widetilde{\tau}_{\xi}^{3}=\rho_{3}e^{\kappa
x},\,\,\widetilde{\tau}_{\xi}^{4}=\rho_{4}e^{\kappa x},$ (4.8)

where $\kappa,\theta_{0},\rho_{a}\,(a=1,...,4)$ are real parameters.

In general, the tau functions become linear combinations of the monomials of
the types $1,e^{\kappa x},e^{2\kappa x},...(\kappa\in\leavevmode\hbox{\rm
I\kern-1.79993ptR})$. So, substituting these tau functions into the relevant
equations in the appendices A, B and C one gets some polynomial equations of
the form $c_{0}^{l}+c_{1}^{l}e^{\kappa x}+c_{2}^{l}e^{2\kappa
x}+....+c_{n}^{l}e^{n_{k}\kappa x}=0,\,(n_{l}\in\leavevmode\hbox{\rm
I\kern-1.79993ptN},\,\,\,l=1,...,7)$. So, by setting to zero the coefficients
(i.e. $c_{n}^{l}\equiv 0$) we will get some parameter relationships and
determine the relevant eigenvalues $E$ depending on the set of parameters
$\\{M,\kappa,\theta_{0},\sigma\\}$. Next, we will implement this process
below.

Firstly, it is convenient to get the simplest forms of the parameter
relationships through the auxiliary equations such as the parity eq. (3.15)
and the consistency equations (3.17) and (3.18), respectively. So, imposing
the parity conditions (3.15) to the $\xi_{a}$ fields defined in (4.2)-(4.5)
for the tau functions (4.6)-(4.8) provide the next parameter relationships

$\displaystyle\rho_{4}$ $\displaystyle=$
$\displaystyle\sigma\rho_{1},\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\theta_{4}=-\theta_{1}+\theta_{0}-\frac{1}{2}\pi$
(4.9) $\displaystyle\rho_{3}$ $\displaystyle=$
$\displaystyle-\sigma\rho_{2},\,\,\,\,\,\,\,\,\,\,\,\theta_{3}=-\theta_{2}+\theta_{0}+\frac{1}{2}\pi.$
(4.10)

Additional parameter relationships can be obtained from the equations (3.17)
and (3.18) written in terms of the tau functions (C.1) and (C.2),
respectively. So, from (C.1) upon substituting the tau functions (4.6)-(4.8)
and taking into account the relationships (4.9)-(4.10) one gets

$\displaystyle\rho_{2}^{2}$ $\displaystyle=$
$\displaystyle\rho_{1}^{2}\csc{(2\theta_{2}-\theta_{0})}\sin{(2\theta_{1}-\theta_{0})}.$
(4.11)

Next, instead of replacing the tau functions (4.6)-(4.8) directly into (C.2)
(which corresponds to (3.18)) it is more convenient to use its simplified
version as in (3.19). In order to get other simple expressions for the
parameters relationships it is convenient to use the integrated form of
(3.19). So, upon substituting the tau function expressions (4.6)-(4.8) into
(3.19) one can integrate it in the interval $[-\infty,x]$ and taking into
account the boundary conditions (3.14) and (3.16) the outcome is a first order
polynomial in $\mbox{sech}(2\kappa x)$. So, setting to zero its zero’th and
first order term coefficients and after some algebraic manipulations taking
into account (4.11) one gets, respectively

$\displaystyle\kappa=\sigma
M\frac{\sin{(\theta_{1}-\theta_{2})}}{\cos{(\theta_{1}+\theta_{2})}-2\sin{\theta_{1}}\cos{\theta_{2}}\cot{\theta_{0}}}$
(4.12)

and

$\displaystyle\rho_{2}$ $\displaystyle=$
$\displaystyle\frac{\rho_{1}}{4}\csc{(2\theta_{2}-\theta_{0})}\times$
$\displaystyle\Big{[}\frac{\cos{(2\theta_{2})}-\cos{(2\theta_{1})}+4\cos{[2(\theta_{1}-\theta_{2})]}-\cos{[2(\theta_{1}-\theta_{0})]}+\cos{[2(\theta_{2}-\theta_{0})]}-4\cos^{2}{(\theta_{1}+\theta_{2}-\theta_{0})}}{2\cos{\theta_{2}}\sin{\theta_{1}}\cos{\theta_{0}}-\cos{(\theta_{1}+\theta_{2})}\sin{\theta_{0}}}\Big{]}.$

Next, taking into account (4.11) and (4.1) one can get an equation relating
only the angles $\\{\theta_{0},\theta_{1},\theta_{2}\\}$

$\displaystyle\cos{[2(\theta_{1}-\theta_{2}+\theta_{0})]}+\cos{[2(\theta_{1}-\theta_{2}-\theta_{0})]}-4\cos{(2\theta_{1})}+4\cos{(2\theta_{2})}+18\cos{[2(\theta_{1}-\theta_{2})]}-$
$\displaystyle
4\cos{[2(\theta_{1}-\theta_{0})]}+4\cos{[2(\theta_{2}-\theta_{0})]}-6\\{1+2\cos{[2(\theta_{1}+\theta_{2}-\theta_{0})]}\\}-2\cos{(2\theta_{0})}=0.$
(4.14)

Now, let us consider the equations of appendices A and B. Substituting the
expressions (4.1)-(4.5), given in terms of the tau functions (4.6)-(4.8), into
the system of equations (A.1)-(A.2) one gets a $8$th order polynomial in
powers of the exponential $e^{\kappa x}$. The even order terms vanish
identically and the odd order terms, i.e. the first, third, fifth and seventh
orders, provide non-vanishing relationships. So, from (A.1) one has

$\displaystyle\kappa\rho_{1}\cos{\theta_{1}}-E\rho_{2}\sin{\theta_{2}}-M\sigma[\rho_{2}\sin{(\theta_{2}-\theta_{0})}\cos{\theta_{0}}+\rho_{1}\sin{\theta_{0}}\cos{(\theta_{1}-\theta_{0})}]=0,$
(4.15)
$\displaystyle\cos{\theta_{0}}[\kappa\rho_{1}\cos{\theta_{1}}-3\rho_{2}(E+M\sigma)\sin{\theta_{2}}]+\sin{\theta_{0}}[3\kappa\rho_{1}\sin{\theta_{1}}-M\rho_{1}\sigma\cos{\theta_{1}}+\rho_{2}(E+2M\sigma)\cos{\theta_{2}}]=0,$
(4.16) $\displaystyle
4\kappa\rho_{1}\cos{\theta_{1}}+\rho_{2}(4E+5M\sigma)\sin{\theta_{2}}+\rho_{2}(2E+M\sigma)\sin{(\theta_{2}-2\theta_{0})}-$
$\displaystyle
2\rho_{1}[\kappa\cos{(\theta_{1}-2\theta_{0})}+M\sigma\cos{(\theta_{1}-\theta_{0})}\sin{\theta_{0}}]=0,$
(4.17)
$\displaystyle\cos{\theta_{0}}[\kappa\rho_{1}\cos{\theta_{1}}+\rho_{2}(E+M\sigma)\sin{\theta_{2}}]-\sin{\theta_{0}}[M\rho_{1}\sigma\cos{\theta_{1}}+E\rho_{2}\cos{\theta_{2}}-\kappa\rho_{1}\sin{\theta_{1}}]=0.$
(4.18)

Similarly, from (A.2) one has

$\displaystyle
2\kappa\rho_{2}\sin{\theta_{2}}+\rho_{1}(2E+M\sigma)\cos{\theta_{1}}+M\sigma\rho_{1}\cos{(\theta_{1}-2\theta_{0})}-2M\rho_{2}\sigma\sin{\theta_{0}}\sin{(\theta_{2}-\theta_{0})}=0,$
(4.19)
$\displaystyle\cos{\theta_{0}}[\kappa\rho_{2}\sin{\theta_{2}}+3\rho_{1}(E+M\sigma)\cos{\theta_{1}}]+\sin{\theta_{0}}[\rho_{1}(E+2M\sigma)\sin{\theta_{1}}-3\kappa\rho_{2}\cos{\theta_{2}}-M\rho_{2}\sigma\sin{\theta_{2}}]=0,$
(4.20)
$\displaystyle-4\kappa\rho_{2}\sin{\theta_{2}}+\rho_{1}(4E+5M\sigma)\cos{\theta_{1}}+\rho_{1}(2E+M\sigma)\cos{(\theta_{1}-2\theta_{0})}+$
$\displaystyle\rho_{2}[2\kappa\sin{(\theta_{2}-2\theta_{0})}-M\sigma\cos{\theta_{2}}+M\sigma\cos{(\theta_{2}-2\theta_{0})}]=0,$
(4.21)
$\displaystyle\cos{\theta_{0}}[\rho_{1}(E+M\sigma)\cos{\theta_{1}}-\kappa\rho_{2}\sin{\theta_{2}}]+\sin{\theta_{0}}[M\rho_{2}\sigma\sin{\theta_{2}}+E\rho_{1}\sin{\theta_{1}}+\kappa\rho_{2}\cos{\theta_{2}}]=0.$
(4.22)

Next, we consider the problem of finding the bound state energy $E$. The four
eqs. (4.15)-(4.18) can be written as an homogeneous linear system of equations
for the unknown independent variables
$\\{\cos{\theta_{1}},\sin{\theta_{1}},\cos{\theta_{2}},\sin{\theta_{2}}\\}$.
So, in order to find a non-trivial solution, one imposes the condition of
vanishing determinant for the $4\times 4$ matrix formed by the coefficients in
that linear system of equations, and one gets a quadratic algebraic equation
for the bound state energy $E$

$\displaystyle[4\kappa^{2}-M^{2}\sin^{2}{(\theta_{o})}]E^{2}+M\sigma[4\kappa^{2}-2M^{2}\sin^{2}{\theta_{o}}-\kappa
M\sigma\sin{(2\theta_{o})}]E-$ $\displaystyle M^{2}\sin{\theta_{o}}[2\kappa
M\sigma\cos{\theta_{o}}-(\kappa^{2}-M^{2})\sin{\theta_{o}}]=0$ (4.23)

Following similar reasoning as above for the system of equations (4.19)-(4.22)
one gets the same algebraic equation (4.23). So, the equation (4.23) is a
second order polynomial in the variable $E$, and its exact solutions become

$\displaystyle E^{\sigma}=-\sigma
M+\frac{M}{4}\Big{\\{}\frac{2[\sigma\pm\cos{\theta_{o}}]+\frac{M}{2\kappa}(1\pm
2\sigma)\sin{\theta_{o}}}{1-(\frac{M}{2\kappa})^{2}\sin^{2}{\theta_{0}}}\Big{\\}}$
(4.24)

Notice that the eigenvalue $E$ depends on the parameters
$\sigma,\,M,\theta_{0}$ and $\kappa$. The fraction $\frac{M}{2\kappa}$, which
appear in the expression for $E^{\sigma}$ in (4.24), can be determined from
(4.12). So, one has

$\displaystyle\frac{M}{2\kappa}=-\frac{\sigma}{2}\cot{\theta_{0}}+\frac{1}{2}\sigma\csc{\theta_{0}}\,[\frac{\sin{(\theta_{0}-2\theta_{+})}}{\sin{(2\theta_{-})}}],$
(4.25)

where

$\displaystyle\theta^{\pm}\equiv\frac{\theta_{1}\pm\theta_{2}}{2}.$ (4.26)

Since the set of angles $\\{\theta_{1},\theta_{2},\theta_{0}\\}$ satisfy the
relationship (4.14) one can write this eq. as

$\displaystyle
6-18\cos{(4\theta^{-})}+4\cos{(2\theta_{0})}\sin^{2}{(2\theta^{-})}+12[1-2\sin^{2}{(\theta_{0}-2\theta^{+})}]+16\cos{\theta_{0}}\sin{(2\theta^{-})}\sin{(\theta_{0}-2\theta^{+})}=0.$

The last equation is a second order polynomial in the ‘variable’
$\sin{(\theta_{0}-2\theta^{+})}$. So, one can write the two solutions in the
form

$\displaystyle\sin{(\theta_{0}-2\theta^{+})}=\frac{1}{3}\Big{[}\cos{\theta_{0}}+\epsilon_{0}\sqrt{2}\sqrt{7+\cos{(2\theta_{0})}}\Big{]}\,\sin{(2\theta^{-})},\,\,\,\,\,\,\,\epsilon_{0}\equiv\pm
1,$ (4.27)

where $\theta^{\pm}$ are defined in (4.26).

Therefore, substituting the expression for
$[\frac{\sin{(\theta_{0}-\theta^{+})}}{\sin{\theta^{-}}}]$ from (4.27) into
the equation (4.25) one can get

$\displaystyle\frac{M}{2\kappa}=-\frac{\sigma}{3}\Big{[}\cot{\theta_{0}}+\epsilon_{2}\frac{\sqrt{2}}{2}\csc{\theta_{0}}\sqrt{7+\cos{(2\theta_{0})}}\Big{]},\,\,\,\epsilon_{2}=\pm
1.$ (4.28)

Remarkably, the last equation shows that $\frac{M}{2\kappa}$ depends only on
the parameters $\\{\sigma,\theta_{0},\epsilon_{2}\\}$. Substituting the
expression for $\frac{M}{2\kappa}$ provided in (4.28) into (4.24) one can get
the eigenvalues $E^{(\sigma)}$

$\displaystyle E^{\sigma}(\theta_{0})=-\sigma
M\left\\{1-\frac{1}{8}\Big{[}\frac{36-6(1-4\sigma\epsilon_{1})\cos{(\theta_{o})}-3\sqrt{2}(1+2\sigma\epsilon_{1})\epsilon_{2}\sqrt{7+\cos{(2\theta_{o})}}}{5-\cos{(2\theta_{o})}-\epsilon_{2}\sqrt{2}\cos{\theta_{o}}\sqrt{7+\cos{(2\theta_{o})}}}\Big{]}\right\\},\,\,\,\,\,\,\,\,\,\epsilon_{1}=\pm
1.$ (4.29)

This is our final expression for the eigenvalue $E$ depending on the set of
parameters $\\{M,\sigma,\theta_{0}\\}$ and the signs $\epsilon_{a}=\pm
1\,(a=1,2)$. Notice that $E^{\sigma}(-\theta_{0})=E^{\sigma}(\theta_{0})$ as
$\kappa$ changes sign under $\theta_{0}\rightarrow-\theta_{0}$. In fact, the
change of sign for $\theta_{0}$ implies different asymptotic values for $\Phi$
at $x=\pm\infty$ and the interchange soliton/anti-soliton.

In addition, combining the relationship for $\frac{M}{2\kappa}$ in (4.12) and
the above equation (4.28) one gets again the equation (4.27). In fact, once a
particular value for $\theta_{0}$ is known one has that the equation (4.27)
defines the allowed set of angles $\\{\theta_{1},\theta_{2}\\}$. It is known
that one-dimensional bound-states are non-degenerate [45]; so, in our case the
bound-states with equal set of $\\{\theta_{1},\theta_{2}\\}$ must be labeled
by their relevant parities. The solutions of (4.27) for each $\theta_{0}$ will
define a continuous family of parameters $\\{\theta_{1},\theta_{2}\\}$ and the
associated parameter $z$ (3.74) of the potential $V_{1}(\Phi)$ in (3.4) will
depend on these parameters through $A_{j}(j=1,2,3)$, as we will show below.

Next, let us analyze the second order equation (3.8). So, substituting the
expressions (4.1)-(4.5), given in terms of the tau functions (4.6)-(4.8), into
the equation (3.8) one gets a $12^{th}$ order polynomial in powers of the
exponential $e^{\kappa x}$. The terms with odd order power and the terms of
order $6,8,10$ and $12$ vanish identically, provided that the previous
parameters relationships are assumed. The terms with the even power $0,2$ and
$4$ provide certain non-trivial relationships between the parameters. So, the
zero’th, second and fourth order power terms furnish the following
relationships, respectively

$\displaystyle A_{1}$ $\displaystyle=$
$\displaystyle-4A_{2}\cos{\theta_{0}}-3[1+2\cos{(2\theta_{0})}]A_{3},$ (4.30)
$\displaystyle A_{2}$ $\displaystyle=$
$\displaystyle\frac{\kappa^{2}}{\beta^{2}\sin^{2}{(\theta_{0})}}-$
$\displaystyle
M\sigma\csc^{2}{\theta_{0}}\Big{\\{}\rho_{1}^{2}\sin{\theta_{1}}\sin{(\theta_{1}-\theta_{0})}+\rho_{2}^{2}\cos{\theta_{2}}\cos{(\theta_{2}-\theta_{0})}-\rho_{1}\rho_{2}\cos{(\theta_{1}-\theta_{2})}\cos{\theta_{0}}\Big{\\}},$
$\displaystyle A_{3}$ $\displaystyle=$
$\displaystyle\frac{1}{6}M\sigma\Big{[}\rho_{1}^{2}+\rho_{2}^{2}-2\rho_{1}\rho_{2}\cos{(\theta_{1}-\theta_{2})}\Big{]}\csc^{2}{(\theta_{0})}.$
(4.32)

In fact, the equation (4.30) is precisely the relationship (3.73) which has
earlier been defined in order to introduce a particular vacua (3.71) of the
model.

####  4.1.1 The tau function approach and the $\theta_{o}$ values

Below, we will find the values of the parameter $\theta_{o}$. We proceed by
establishing additional parameter relationships from (3.28)-(3.29) and (3.30)
for the above solutions (4.1)-(4.5). So, let us write the relevant expressions
for the charge densities of the topological and non-local currents as

$\displaystyle{\cal J}^{0}_{nonl}(x)=\frac{a_{0}+a_{1}\cosh{(2\kappa
x)}+a_{2}\cosh{(4\kappa x)}}{4\kappa[\cos{\theta_{0}}+\cosh{(2\kappa
x)}]^{2}}-\frac{a_{2}}{2\kappa},$ (4.33)

with

$\displaystyle a_{0}$ $\displaystyle\equiv$ $\displaystyle
16A_{2}\cos{\theta_{0}}\sin^{3}{(\theta_{0})}+A_{3}\\{(45+\log{8})\sin{\theta_{0}}-36\sin{(2\theta_{0})}+3[7+\log{16}]\sin{(3\theta_{0})}-$
(4.35) $\displaystyle 9\sin{(4\theta_{0})}+\log{8}\sin{(5\theta_{0})}\\},$
$\displaystyle a_{1}$ $\displaystyle\equiv$ $\displaystyle
12(A_{2}-3A_{3})\sin{\theta_{0}}+4A_{3}(15+\log{8})\sin{(2\theta_{0})}-4(A_{2}+9A_{3})\sin{(3\theta_{0})}+$
(4.37) $\displaystyle 6A_{3}(1+\log{4})\sin{(4\theta_{0})},$ $\displaystyle
a_{2}$ $\displaystyle\equiv$ $\displaystyle
9A_{3}[\sin{\theta_{0}}-2\sin{(2\theta_{0})}].$ (4.38)

Notice that the above parameters $a_{2}$ and $a_{3}$ satisfy the identity:
$2\kappa a_{3}+a_{2}=0$. So, the non-local current (4.33) satisfies the
boundary condition

$\displaystyle{\cal J}^{0}_{nonl}(\pm\infty)=0.$ (4.39)

Next, the topological charge density $j^{0}_{top}$ can be written as

$\displaystyle-\frac{4\pi}{\beta^{2}}\,j^{0}_{top}$ $\displaystyle=$
$\displaystyle-\frac{2}{\beta}\,\partial_{x}\Phi$ (4.40) $\displaystyle=$
$\displaystyle-\frac{4\kappa}{\beta^{2}}\,[\frac{\sin{\theta_{0}}}{\cos{\theta_{0}}+\cosh{(2\kappa
x)}}].$

By integrating the r.h.s. of (3.30) one can get the charge density $J^{0}$
written as

$\displaystyle J^{0}$ $\displaystyle=$
$\displaystyle\frac{4}{[\cos{\theta_{0}}+\cosh{(2\kappa
x)}]^{2}}[b_{0}+b_{1}\cosh{(2\kappa x)}],$ (4.41)

where the expressions for $b_{0}$ and $b_{1}$ are provided in the appendix D.

Next, let us consider the current charge densities (4.33), (4.40) and (4.41)
and substitute them into the relationship (3.28). So, one gets the next first
order polynomial in $\cosh{(2\kappa x)}$

$\displaystyle c_{0}+c_{1}\cosh{(2\kappa x)}$ $\displaystyle\equiv$
$\displaystyle 0,$ (4.42)

where an overall factor $[\cos{\theta_{0}}+\cosh{(2\kappa x)}]^{-2}$ has been
omitted and the coefficients $c_{0},c_{1}$ are provided in the appendix E.

So, the coefficients in the identity (4.42) must satisfy $c_{0}=c_{1}=0$.
Next, take into account the relationships between $\rho_{1}^{2}$ and
$\rho_{2}^{2}$ in (4.11), $\kappa$ and $M$ in (4.12) and $\rho_{1}$ and
$\rho_{2}$ in (4.1). So, taking into account the form of $c_{0}$ in (E.1) and
setting $c_{0}=0$ according to the identity (4.42), one can get the following
equation

$\displaystyle[3+5\cos{(2\theta_{0})}+\cos{(4\theta_{0})}]\,\csc{(2\theta_{2}-\theta_{0})}\,\csc^{2}{\theta_{0}}\,\sin{[2(\theta_{1}-\theta_{2})]}=0.$
(4.43)

Similarly, setting $c_{1}=0$ in (4.42) and taking into account (E.2) one can
get

$\displaystyle[1+2\cos{(2\theta_{0})}]\,\cos{(\theta_{1}-\theta_{2})}\cot{\theta_{0}}\csc{(2\theta_{2}-\theta_{0})}\sin{(\theta_{1}-\theta_{2})}^{2}=0.$
(4.44)

Notice that the equation for $\kappa$ (4.12) implies

$\displaystyle[\sin{\theta_{0}}\cos{(\theta_{1}+\theta_{2})}-2\cos{\theta_{0}}\sin{\theta_{1}}\cos{\theta_{2}}]\neq
0.$ (4.45)

Moreover, in order to define real values for $\rho_{1}$ and $\rho_{2}$ in
(4.11) one must have

$\displaystyle\theta_{1}-\theta_{2}\neq\frac{n\pi}{2}.$ (4.46)

Then the equations (4.43)-(4.44) must hold provided that (4.45) and (4.46) are
satisfied. So, the factors involving the angle $\theta_{0}$ in (4.43),
$[3+5\cos{(2\theta_{0})}+\cos{(4\theta_{0})}]$ and $[1+2\cos{(2\theta_{0})}]$
in (4.44) provide, upon equating them to zero, the common solutions

$\displaystyle\theta_{0}=\pm\frac{\pi}{3}+\pi
n,\,\,\,\,\,n\in\leavevmode\hbox{\sf Z\kern-3.99994ptZ}.$ (4.47)

Notice that $|\theta_{0}|\leq\pi$ in the kink(anti-kink) solution (4.51) and
this parameter determines the asymptotic value of the kink (anti-kink) at
$x=\pm\infty$. Therefore, one can choose the set

$\displaystyle\theta_{0}=\\{\pm\frac{\pi}{3},\,\,\pm\frac{2\pi}{3}\\}.$ (4.48)

These are the possible values of the parameter $\theta_{0}$ defined in (4.1)
and the tau functions $\tau_{0}$ and $\tau_{1}$ in (4.6) . So, for the
$\theta_{0}^{\prime}s$ in (4.48) and taking into account the expressions
(4.28) and 4.29) one can compute $\kappa$ and the energy eigenvalues $E$,
respectively. In addition, the $\theta_{0}$ values and the eq. (4.14) will
determine a relationship between the phases $\\{\theta_{1},\theta_{2}\\}$ of
the spinor components in (4.2)-(4.3). The other phases
$\\{\theta_{3},\theta_{4}\\}$ in (4.4)-(4.5) can be determined from
(4.9)-(4.10).

In the Figs. 4-7 we present the relevant values of the set
$\\{\frac{\kappa}{M},\,\frac{E}{M}\\}$ for the angles provided in (4.48) and a
given parity $\sigma$. A careful observation of the tables and points for the
both parities allows one to deduce that the discrete points with coordinates
$(\kappa\,,\,E)$ (for a fixed value of $M>0$) satisfy the symmetry

$\displaystyle\kappa\rightarrow-\kappa\,\,\mbox{and}\,\,\sigma\rightarrow-\sigma\,\mbox{implies}\,\,\,\,E\rightarrow-E.$
(4.49)

In the case $\theta_{0}=\pi/3$ (tables and points in the Figs. 4 and 5) a
careful examination allow us to say that all of the states become in-gap
states such that $|\frac{E}{M}|<1$ and satisfy the symmetry (4.49).

In the case $\theta_{0}=2\pi/3$ (see the tables in Figs. 6 and 7) one has the
both type of states, in-gap and BIC states. Remarkably, one notices the
appearance of some bound states in the continuum (BIC) with energies
$|\frac{E}{M}|>1$; i.e. $E=\pm 1.01M,\pm 1.1M$. These bound states remain
localized even though they coexist with the continuous spectrum of radiating
waves that can carry energy away, i.e. for energies above or below the
threshold ($|E|>M$) as defined in (3.43). These types of bound states were
first predicted by von Neumann and Wigner [46] in 3D linear Schrödinger
equation with a localized potential. For an account of the developments in the
BIC theory and the explanation of the physical mechanisms responsible for
their appearance in different physical areas see e.g. [47]. In the context of
soliton theory the BIC states appeared as embedded solitons which can be
robust despite being in resonance with the linear spectrum. For different
mechanisms underlying the existence of these solitons see e.g. [48, 49]).

Figure 4: The left table shows the $E/M$ and $\kappa/M$ values for
$\theta_{0}=\pi/3,\sigma=1$.The right Fig. shows the points with coordinates
$(\kappa/M\,,\,E/M)$.

Figure 5: The left table shows the $E/M$ and $\kappa/M$ values for
$\theta_{0}=\pi/3,\sigma=-1$.The right Fig. shows the points with coordinates
$(\kappa/M\,,\,E/M)$.

Figure 6: The left table shows the $E/M$ and $\kappa/M$ values for
$\theta_{0}=2\pi/3,\sigma=1$.The right Fig. shows the points with coordinates
$(\kappa/M\,,\,E/M)$. Notice the BIC states for $E=+1.01M,-1.1M$.

Figure 7: The left table shows the $E/M$ and $\kappa/M$ values for
$\theta_{0}=2\pi/3,\sigma=-1$.The right Fig. shows the points with coordinates
$(\kappa/M\,,\,E/M)$. Notice the BIC states for $E=-1.01M,+1.1M$.

Let us write the explicit forms of the fermion bound states and the kink
solutions. The fermion components become

$\displaystyle\xi_{a}=(-i)^{a-1}\,\rho_{a}\,e^{\kappa
x-i\theta_{a}}\Big{[}\frac{1}{1+e^{2\kappa
x-i\theta_{0}}}+(-1)^{a-1}\frac{e^{2i\theta_{a}}}{1+e^{2\kappa
x+i\theta_{0}}}\Big{]},\,\,\,\,a=1,2.$ (4.50)

The remaining components $\xi_{a}\,(a=3,4)$ can be obtained from the above
expressions (4.50) taking into account the parity condition relationships
(3.15). Notice that the expressions (4.50) define real $\xi_{a}$ functions.

Next, taking into account the relationship (4.1) and the tau functions (4.6),
as wel as the conditions (3.13)-(3.14), one can write the kink as

$\displaystyle\Phi(x)$ $\displaystyle=$
$\displaystyle\frac{2}{\beta}\arctan{\Big{[}\tan{(\frac{\theta_{0}}{2})}\,\,\tanh{(\kappa
x)}\Big{]}},\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\theta_{0}\in\\{\pm\frac{\pi}{3},\,\,\pm\frac{2\pi}{3}\\}.$
(4.51)

The asymptotic values of the field are
$\Phi(\pm\infty)\equiv\pm\frac{\theta_{0}}{\beta}$.

In addition, a condition for the spinor will be given by setting the spinor
component $\xi_{1}$ to a constant at $x=0$, i.e.

$\displaystyle\xi_{1}(0)\equiv\xi_{0}=\rho_{1}\cos{\theta_{1}}[1+\tan{\theta_{1}}\tan{(\frac{\theta_{0}}{2})}],$
(4.52)

where $\xi_{0}$ is a real constant.

Figure 8: The kinks $\Phi(x)$ for the set of parameters $\kappa=1,\beta=1$,
$\theta_{0}=\\{2\pi/3,\pi/3\\},\,\phi_{0}=\pi$. The kink with $\phi_{0}$ value
will host a Majorana fermion.

Figure 9: The bound state of the fermion as a function of $x$ for the positive
parity $\sigma=1$. It is plotted for
$\rho_{1}=0.41,\rho_{2}=-0.65,\,\theta_{1}=0.138,\,\theta_{2}=0.85,\,\kappa=1,\theta_{0}=2\pi/3$.
Notice the spinor components relationships $\xi_{1}(-x)=-\xi_{4}(x)$ and
$\xi_{2}(-x)=+\xi_{3}(x)$

Figure 10:  The bound state of the fermion as a function of $x$ for the
negative parity $\sigma=-1$. It is plotted for
$\rho_{1}=0.41,\rho_{2}=-0.64,\,\theta_{1}=0.14,\,\theta_{2}=0.85,\,\kappa=1,\theta_{0}=2\pi/3$.
Notice the spinor components relationships $\xi_{1}(-x)=+\xi_{4}(x)$ and
$\xi_{2}(-x)=-\xi_{3}(x)$

The above (4.51) kink-type solutions for the values
$\\{\theta_{0}=\,\pi/3,\,2\pi/3\\}$ are plotted in the Fig. 8. For comparison
we also plotted the sine-Gordon type kink,
$\frac{4}{\beta}\arctan{\Big{[}e^{2\kappa x}\Big{]}}-\frac{\pi}{\beta}$, with
asymptotic values $\phi_{o}=\pm\pi/\beta$. The sine-Gordon type kink also
appear as a solution of the model for another set of tau functions and this
kink will host a Majorana fermion, as we will show below in the section 6 .

One can define the topological charge as

$\displaystyle Q_{topol}$ $\displaystyle\equiv$
$\displaystyle\frac{\beta}{2\pi}\,\int_{-\infty}^{+\infty}\partial_{x}\Phi$
(4.53) $\displaystyle=$
$\displaystyle\frac{\beta}{2\pi}[\Phi(+\infty)-\Phi(-\infty)].$ (4.54)

In the case at hand it becomes $Q_{topol}=\frac{\theta_{0}}{\pi}$. So, our
analytical solutions are associated to the non-integer topological charge
values

$\displaystyle Q_{topol}=\pm\frac{1}{3},\,\pm\frac{2}{3},$ (4.55)

where the $+$ sign corresponds to solitons and the $-$ sign for antisolitons.

The fermion bound state corresponding to the soliton ($\theta_{0}>0$) can be
normalized as follows

$\displaystyle\frac{1}{{\cal N}}\int_{-\infty}^{+\infty}dx\,J^{0}$
$\displaystyle=$ $\displaystyle 1$ (4.56) $\displaystyle{\cal N}$
$\displaystyle\equiv$
$\displaystyle\frac{2}{\kappa}\Big{[}\rho_{1}^{2}\cos{(2\theta_{1}-\theta_{0})}-\rho_{2}^{2}\cos{(2\theta_{2}-\theta_{0})}+\theta_{0}(\rho_{1}^{2}+\rho_{2}^{2})\csc{\theta_{0}}\Big{]}.$
(4.57)

where the charge density $J^{0}=\sum_{a=1}^{4}\xi_{a}^{2}$ as defined in
(3.25) becomes

$\displaystyle J^{0}$ $\displaystyle=$
$\displaystyle\frac{4}{[\cos{\theta_{0}}+\cosh{(2\kappa
x)}]^{2}}\Big{\\{}\rho_{1}^{2}\cos{(\theta_{1}-\theta_{0})}\cos{\theta_{1}}+\rho_{2}^{2}\sin{(\theta_{2}-\theta_{0})}\sin{\theta_{2}}+$
(4.58)
$\displaystyle\frac{1}{2}[\rho_{1}^{2}+\rho_{2}^{2}+\rho_{1}^{2}\cos{(2\theta_{1}-\theta_{0})}\cos{\theta_{0}}-\rho_{2}^{2}\cos{(2\theta_{2}-\theta_{0})}\cos{\theta_{0}}]\cosh{(2\kappa
x)}\Big{\\}}.$

The normalization condition ${\cal N}=1$ in (4.56) and the eq. (4.11) can be
used to determine the particular values of $\\{\rho_{1},\rho_{2}\\}$, once the
angles $\\{\theta_{1},\theta_{2}\\}$ are fixed from the curves defined by
(4.27) for each $\theta_{0}$.

Regarding the full set of parameters of the system of equations (3.4)-(3.8)
and the solution (4.50)-(4.51) one has $17$ parameters, i.e.
$\\{\rho_{a},\,\theta_{a},\,A_{j},\,\theta_{0},\,\kappa,E,\,M,\,\beta,\sigma\\}\,\mbox{for}\,a=1,2,3,4\,\mbox{and}\,j=1,2,3$.
One can consider the parameters $\beta,M$ (which appear in the Lagrangian) and
$\rho_{1}$ as free parameters. In fact, $\rho_{1}$ can be fixed from the
condition (4.52) once a particular value for $\xi_{0}$ is chosen (this value
can be set according to the normalization condition (4.56)-(4.57)). In
addition, once a particular value for $\theta_{0}$ is chosen from the set
(4.48) and a parity eigenvalue $\sigma$ is determined for a particular
solution one can get the parameter $\kappa$ and the eigenvalue $E$, as
functions of $M$, from the eqs.(4.28) and 4.29), respectively. The $\rho_{2}$
can be determined from (4.11), and the set $\\{\theta_{3},\theta_{4}\\}$ and
$\\{\rho_{3},\rho_{4}\\}$ can be obtained from (4.9)-(4.10). So, the values
for $A_{j}$ can be determined from (4.30)-(4.32).

Moreover, the equation (4.27) defines a family of curves in the plane
$\theta_{1}\times\theta_{2}$ for each $\theta_{0}$ of the set (4.48). Notice
that a choice of a particular point $(\theta_{1},\theta_{2})$ in a given curve
defines a parameter value $z$ in (3.74) corresponding to the potential of type
$V_{1}(\Phi)$ in (3.4). In fact, the parameter $z$ depends, through the
parameters $A_{j}\,(j=1,2,3)$ in (4.30)-(4.32), on the angles
$\theta_{1},\theta_{2},\theta_{0}$ and the parameters $M,\sigma,\beta,\kappa$.

Figure 11: The Figs. show the family of curves in the plane
$\theta_{1}\times\theta_{2}$ defined by (4.27) for $\theta_{0}=\frac{2\pi}{3}$
and the cases $\epsilon_{0}=+1$(left figure) and $\epsilon_{0}=-1$(right
figure).

In the Figs. 9 and 10 we present, for some set of parameters, the plots of the
solutions $\xi_{a},\,\,a=1,...4$ defined in (4.2)-(4.5) considering the tau
functions (4.6)-(4.8) for the both signs of the parity eigenvalues $\sigma=\pm
1$, respectively. Their relevant parameters satisfy the relationships
discussed above, therefore the properties of the solitons in Fig. 8, such as
their asymptotic values, width and slope, will always depend on the relevant
spinor parameters, and vice versa.

In the Fig. 11 we show a family of curves in the plane
$\theta_{1}\times\theta_{2}$ as defined through the eq. (4.27) for
$\theta_{0}=\frac{2\pi}{3}$ and the cases $\epsilon_{0}=+1$(left figure) and
$\epsilon_{0}=-1$(right figure).

##  5 Back-reaction of spinor on the kink and strong coupling sector

Fermions bounded by kinks were previously considered by some authors in the
literature. However, in most of the studies the effect of the back-reaction of
the fermions on the soliton were considered to be nearly zero or completely
neglected. Some attempts to take into account the back-reaction of the fermion
on the kink have been performed mainly by numerical simulations [13, 30],
although self-consistent analytical solutions, to our knowledge, are still
missing. The main reason for that is the lack of integrability of the models
considered and the non-existence of general exact methods to solve the highly
nonlinear and non-integrable models, as well as the integral-differential
equations to solve the self-consistent system of equations.

Next, we examine the back-reaction of the Dirac fermion on the kink using the
exact solutions obtained above. So, from (4.50) and (4.51) one can construct
the following mappings between the spinor components and scalar field

$\displaystyle\xi_{1}$ $\displaystyle=$
$\displaystyle\frac{2\rho_{1}}{\sin{\theta_{0}}}\sqrt{\sin{(\frac{\theta_{0}-\beta\Phi}{2})}\sin{(\frac{\beta\Phi+\theta_{0}}{2})}}\,\,\,\,\cos{(\theta_{1}-\frac{\theta_{0}+\beta\Phi}{2})},$
(5.1) $\displaystyle\xi_{2}$ $\displaystyle=$
$\displaystyle-\frac{2\rho_{2}}{\sin{\theta_{0}}}\sqrt{\sin{(\frac{\theta_{0}-\beta\Phi}{2})}\sin{(\frac{\beta\Phi+\theta_{0}}{2})}}\,\sin{(\theta_{2}-\frac{\theta_{0}+\beta\Phi}{2})},$
(5.2) $\displaystyle\xi_{3}$ $\displaystyle=$
$\displaystyle-\frac{2\sigma\rho_{2}}{\sin{\theta_{0}}}\sqrt{\sin{(\frac{\theta_{0}-\beta\Phi}{2})}\sin{(\frac{\beta\Phi+\theta_{0}}{2})}}\,\sin{(\theta_{2}-\frac{\theta_{0}-\beta\Phi}{2})},$
(5.3) $\displaystyle\xi_{4}$ $\displaystyle=$
$\displaystyle-\frac{2\sigma\rho_{1}}{\sin{\theta_{0}}}\sqrt{\sin{(\frac{\theta_{0}-\beta\Phi}{2})}\sin{(\frac{\beta\Phi+\theta_{0}}{2})}}\,\,\cos{(\theta_{1}-\frac{\theta_{0}-\beta\Phi}{2})}.$
(5.4)

Notice that these spinor components satisfy the parity symmetry defined in
(3.15) provided that the kink satisfies (3.13). In fact, the kink presented in
(4.51) satisfies (3.13).

Moreover, using the above mappings one can get the following identity

$\displaystyle
2M\beta\Big{[}(\xi_{1}\xi_{3}+\xi_{2}\xi_{4})\cos{\beta\Phi}-(\xi_{1}\xi_{4}-\xi_{2}\xi_{3})\sin{\beta\Phi}\Big{]}=\beta
B_{1}\sin{(\beta\Phi)}+2\beta B_{2}\sin{(2\beta\Phi)}+3\beta
B_{3}\sin{(3\beta\Phi)},$ (5.5)

with

$\displaystyle B_{1}$ $\displaystyle=$
$\displaystyle\frac{M\sigma}{2\sin^{2}{\theta_{0}}}\Big{\\{}\rho_{1}^{2}+\rho_{2}^{2}-2\rho_{1}\rho_{2}\cos{(\theta_{1}-\theta_{2})}-4\rho_{1}^{2}\cos{(2\theta_{1}-\theta_{0})}\cos{\theta_{0}}+4\rho_{2}^{2}\cos{(2\theta_{2}-\theta_{0})}\cos{\theta_{0}}\Big{\\}}$
$\displaystyle B_{2}$ $\displaystyle=$
$\displaystyle-\frac{M\sigma}{2\sin^{2}{\theta_{0}}}\Big{\\{}[\rho_{1}^{2}+\rho_{2}^{2}-2\rho_{1}\rho_{2}\cos{(\theta_{1}-\theta_{2})}]\cos{\theta_{0}}-\rho_{1}^{2}\cos{(2\theta_{1}-\theta_{0})}+\rho_{2}^{2}\cos{(2\theta_{2}-\theta_{0})}\Big{\\}}$
(5.6) $\displaystyle B_{3}$ $\displaystyle=$
$\displaystyle\frac{M\sigma}{6\sin^{2}{\theta_{0}}}\Big{\\{}\rho_{1}^{2}+\rho_{2}^{2}-2\rho_{1}\rho_{2}\cos{(\theta_{1}-\theta_{2})}\Big{\\}}.$
(5.7)

One can verify the next relationship

$\displaystyle
B_{1}=-4B_{2}\cos{\theta_{0}}-3[1+2\cos{(2\theta_{0})}]\,B_{3}.$ (5.8)

Notice that the r.h.s. of the identity (5.5) depends only on the scalar field.
We proceed by constructing a static equation of motion involving only the
scalar field $\Phi$. This follows by making use of the identity (5.5) into the
static version of the equation (3.8). In fact, the identity (5.5) defines an
equivalence of the spinor-scalar interactions terms (l.h.s of (5.5)) to a sum
of multi-frequency sine terms (r.h.s. of (5.5)). So, substituting the l.h.s of
(5.5) into the relevant terms of the static version of the equation of motion
(3.8) one can get

$\displaystyle\Phi^{\prime\prime}+\beta C_{1}\sin{(\beta\Phi)}+2\beta
C_{2}\sin{(2\beta\Phi)}+3\beta C_{3}\sin{(3\beta\Phi)}=0,$ (5.9)

where

$\displaystyle C_{j}=A_{j}-B_{j},\,\,\,\,\,j=1,2,3.$ (5.10)

Inspecting the form of the parameters $A_{3}$ in (4.32) and $B_{3}$ in (5.7)
one concludes that $B_{3}=A_{3}$, then

$\displaystyle C_{3}=0.$ (5.11)

So, (5.9) defines a static equation of motion for a double SG model, such that
an effective potential can be defined as

$\displaystyle V_{DSG}=C_{1}\cos{(\beta\Phi)}+C_{2}\cos{(2\beta\Phi)}.$ (5.12)

The equation (5.9) can be regarded as an effective model for the scalar field
$\Phi_{DSG}$ which incorporates the back-reaction of the fermion on the kink
from the original system (3.4)-(3.8). The effective model (5.9) carries the
effect of the spinor on the kink through the parameters $C_{j}$ which depend
on $B_{j}$ (5.10). In fact, the $B_{j}$ depend on the spinor parameters
$\\{\rho_{1,2},\theta_{1,2},\sigma,\theta_{0},\kappa\\}$. The parameters
$C_{j}$ also incorporate, through the terms $A_{i}$, the kink parameters
$\\{\kappa,\theta_{0}\\}$, as well as the parameters $\\{\beta,M\\}$ which
appear in the Lagrangian.

Taking into account (4.30) and (5.8) one can show the next relationship

$\displaystyle C_{1}=-4C_{2}\cos{\theta_{0}}.$ (5.13)

Next, let us examine the vacua of the potential $V_{DSG}$ in (5.12). Since the
parameters satisfy (5.13) one gets the critical points of the potential (5.12)
to be

$\displaystyle\Phi_{DSG}^{vac}=\frac{\pi}{\beta}m,\,\,\pm\frac{\theta_{0}}{\beta}+\frac{2\pi}{\beta}n;\,\,\,\,\,m,n\in\leavevmode\hbox{\sf
Z\kern-3.99994ptZ}.$ (5.14)

Remarkably, these vacua take the same values as the ones we have defined
corresponding to the special potential in (3.71).

The model (5.9) with potential (5.12) and the parameters relationships (5.11)
and (5.13) possesses the following two types of kink solutions.

Type A

By direct substitution in the static equation of motion (5.9) for $C_{3}=0$
and $C_{1}=-4C_{2}\cos{\theta_{0}}$, one can verify the solution

$\displaystyle\Phi_{DSG}^{A}=\frac{2}{\beta}\arctan{\Big{[}\tan{(\frac{\theta_{0}}{2})}\tanh{(\gamma
x)}\Big{]}},$ (5.15)

with

$\displaystyle\gamma=\sqrt{\beta C_{2}}\sin{\theta_{0}},\,\,\,C_{2}>0.$ (5.16)

Notice that the solution (5.15) resembles to the one in (4.51). However, we
will show below that the $\theta_{0}$ dependence of their relevant widths
$\frac{1}{\kappa}$ and $\frac{1}{\gamma}$, respectively for (4.51) and (5.15),
are completely different in the both cases. In fact, the solution (4.51) has
been self-consistently obtained for the soliton-fermion system (3.4)-(3.8)
with potential $V$ defined for $\beta_{j}=j\,\beta\,(j=1,2,3)$) such that the
parameter $\theta_{0}$ takes the particular values of
$\theta_{0}=\pm\frac{\pi}{3},\,\pm\frac{2\pi}{3}$.

So, in the model (3.4)-(3.8) one can examine the degree of the back-reaction
of the spinor components on the soliton by comparing the widths of the both
solitons $\Phi$ in (4.51) ($\sim\frac{1}{\kappa}$) and $\Phi^{A}_{DSG}$ in
(5.15) ($\sim\frac{1}{\gamma}$), respectively. In doing so, one can argue that
the soliton (5.15) of the model (5.9)-(5.11) carries the spinor back-reaction
of the soliton-fermion system (3.4)-(3.8) through the parameter $\gamma$,
which defines its with as $\sim\frac{1}{\gamma}$.

Figure 12: The kinks $\Phi$ (green) and $\Phi_{DSG}^{A}$ (red) as a function
of $x$ for weak and large coupling parameter $\beta$. Left: For
$\beta=0.015,\,\gamma=4.6$ . Right: For $\beta=1.85,\,\gamma=0.09$. In the
both cases
$\theta_{0}=2\pi/3,M=1,\kappa=0.56,\rho_{1}=0.41,\rho_{2}=-0.65,\theta_{1}=0.138,\theta_{2}=0.85,\sigma=1$.

The changes in the local features of the kink, once the asymptotic value
$\theta_{0}=|\Phi(\pm\infty)|$ is fixed, will be considered below in order to
study the back-reaction of the fermion on the kink; so, let us compute the
slope of the kink $\Phi^{A}_{DSG}$ at the origin $x=0$

$\displaystyle\mu^{A}=\frac{2}{\beta}\gamma\,\tan{(\theta_{0}/2)},$ (5.17)

where

$\displaystyle\gamma$ $\displaystyle=$
$\displaystyle\frac{2\sqrt{2}}{\beta^{3/2}}\Big{\\{}2\kappa^{2}+M\beta^{2}\sigma(2\beta-1)\times$
$\displaystyle\Big{[}\left(\rho_{1}^{2}+\rho_{2}^{2}-2\rho_{1}\rho_{2}\cos{(\theta_{1}-\theta_{2})}\right)\cos{\theta_{0}}-\rho_{1}^{2}\cos{(2\theta_{1}-\theta_{0})}+\rho_{2}^{2}\cos{(2\theta_{2}-\theta_{0})}\Big{]}\Big{\\}}^{1/2},$

such that $\kappa$, in terms of $M,\theta$ and $\sigma$, provided in (4.28)
must be taken into account.

On the other hand, the slope of the kink $\Phi$ (4.51) at the origin $x=0$
becomes

$\displaystyle\mu=\frac{2}{\beta}\kappa\,\tan{(\theta_{0}/2)}.$ (5.19)

In (4.28) one has $\kappa=f(\theta_{0},\sigma)M$ with $f(\theta_{0},\sigma)$
being a proportionality factor; so, the slope $\mu$ of the kink $\Phi$ defined
in (5.19) is directly proportional to the parameter $M$ which appears in the
Lagrangian (2.1) as the fermion mass in the free fermion limit of the model
($\beta\rightarrow 0$). Whereas, $\gamma$ in (5) carries the effect of all the
parameters of the model, including the spinor parameters. Therefore, one can
argue that the slope $\mu^{A}$ in (5.17) measures the degree of back-reaction
of the spinor components on the kink. Notice that the slopes satisfy

$\displaystyle\frac{\mu^{A}}{\mu}=\frac{\gamma}{\kappa}.$ (5.20)

Let us analyze the strong coupling sector ($\beta\rightarrow$ large) of the
model by considering the kink solutions. So, for finite $M$ (remember that
$\kappa\sim M$) and fixed $f(\theta_{0},\sigma)$ the slope of the kink $\Phi$
(5.19) in the limit $\beta\rightarrow$ large, tends to zero i.e.
$\mu\rightarrow 0$. However, it is possible to obtain a kink $\Phi^{A}_{DSG}$
with finite slope in (5.17) by maintaining $\frac{\gamma}{\beta}=$ finite in
the limit $\beta\rightarrow$ large. In fact, taking into account (5) and the
relationship between the parameters $\rho_{1}$ and $\rho_{2}$ in (4.11), in
the limit ($\beta\rightarrow$ large) one can consider $\rho_{1}\rightarrow$
large, in order to get

$\displaystyle\mu^{A}$ $\displaystyle=$
$\displaystyle[\frac{\gamma}{\beta}]2\tan{\frac{\theta_{0}}{2}}$ (5.21)
$\displaystyle\sim$
$\displaystyle\sqrt{M}\,\frac{\rho_{1}}{\beta}=\mbox{finite},$ (5.22)

where a set of values for $\\{\theta_{1},\theta_{2},\theta_{0}\\}$ and the
parity $\sigma$ have been fixed.

In view of the above discussions, one can argue that the double SG model (5.9)
describes the strong coupling sector of the model defined by the Lagrangian
(2.1) for kink solitons of the type (5.15). In order to see qualitatively the
back-reaction of the fermion $\psi$ on the soliton we plot the solitons of
types $\Phi$ and $\Phi_{DSG}^{A}$ in the Fig. 12 for weak and strong coupling
constant parameter $\beta$ and the particular value
$\theta_{0}=\frac{2\pi}{3}$.

Type B

The model (5.9)-(5.11) also exhibits the next kink and bounce type solutions
[19]. Consider

$\displaystyle\Phi_{DSG}^{B}=\frac{2}{\beta}\arctan{\Big{[}(\frac{1}{d})\frac{1+h\exp{(2kx)}}{\exp{(kx)}}\Big{]}}.$
(5.23)

The kink type solution arises for the relationships

$\displaystyle
C_{2}=-\frac{1}{4\beta}(\frac{k^{2}}{1+\cos{\theta_{0}}}),\,\,\,\,\,h=-\frac{1}{4}(\frac{d^{2}\cos{\theta_{0}}}{1+\cos{\theta_{0}}}),\,\,\,d=2\sqrt{1+\sec{\theta_{0}}}.$
(5.24)

For the solution $\theta_{0}=\frac{\pi}{3}$ in (4.48) one gets $d=2\sqrt{3}$.
Then, one has the kink-type solution

$\displaystyle\Phi_{DSG}^{B}=\frac{2}{\beta}\arctan{\Big{[}\frac{1}{\sqrt{3}}\sinh{kx}\Big{]}}.$
(5.25)

This solutions interpolates the values $\pm\pi$ of the set of points in the
vacua (5.14).

Moreover, for $\theta_{0}=\frac{2\pi}{3}$ one has $d=2$ and the next bounce-
like solution

$\displaystyle\Phi_{DSG}^{bounce}=\frac{2}{\beta}\cosh{(kx)}.$ (5.26)

##  6 Zero-mode bound states

Regarding the zero modes, notice that the both conditions (3.68) and (3.69) in
order to have Majorana fields satisfy trivially the relationship (3.17). Then,
we do not have at our disposal a relationship between the parameters analogous
to the equation (4.11), which has been used in the calculations to arrive at
the fermion spectra for $E\neq 0$. Consequently, the constructions above
performed to arrive at the equations (4.12)-(4.14) and (4.43)-(4.44) do not
hold anymore. So, one must search for new bound state solutions to the whole
system of equations (3.4)-(3.8) specialized to the zero-modes $E=0$.

So, from this point forward in this section we set $E=0$ in the system of
equations (3.4)-(3.7). So, let us set $E=0$ in the equivalent system of
equations (3.61). So, one has in components

$\displaystyle\frac{d}{dx}\xi_{0R}$ $\displaystyle=$ $\displaystyle-
Me^{-i\beta\Phi}\xi_{0L},$ (6.1) $\displaystyle\frac{d}{dx}\xi_{0L}$
$\displaystyle=$ $\displaystyle-Me^{i\beta\Phi}\xi_{0R},$ (6.2)

where $\xi_{0R}=\xi_{3}+i\xi_{4}$ and $\xi_{0L}=\xi_{1}+i\xi_{2}$. Notice that
the system (6.1)-(6.2) describes a Majorana fermion coupled to a scalar field,
provided that either the identification (3.68) or the one in (3.69) is taking
into account. In fact, the case with $\beta=0$ describes a free massive
Majorana system [50, 51]. Recently, the model (6.1)-(6.2) has been considered
as the simplest continuum model of a one-dimensional superconducting fermionic
symmetry-protected topological (SPT) phase in condensed matter such that the
mean-field superconducting (SC) pairing potential is represented by the field
$\Phi$ [42, 41].

Moreover, the static second order equation derived from (3.8) for the kink
$\Phi$ and fermion components $\xi_{0}^{T}=\left(\xi_{0R},\,\xi_{0L}\right)$
becomes

$\displaystyle-\frac{d^{2}\Phi}{dx^{2}}+M\beta[\xi^{\star}_{0R}\xi_{0L}e^{i\beta\Phi}+\xi^{\star}_{0L}\xi_{0R}e^{-i\beta\Phi}]+V^{\prime}(\Phi)=0.$
(6.3)

Then, in order to find the kinks and zero-mode fermion bound states one must
solve the system of equations (6.1)-(6.3). However, in the following we will
construct those solutions by solving a system of first order integro-
differential equations. Notice that the integro-diferential equation (3.27)
has been written for any bounded static configuration associated to an
arbitrary value of the parameter $E$, and in particular, it must hold for the
zero-modes with $E=0$. So, let us consider the identity (3.27) written for the
field $\xi_{0}$ and the kink $\Phi$

$\displaystyle\frac{1}{\beta}\frac{d\Phi}{dx}+\frac{1}{2}(\xi^{\star}_{0R}\xi_{0R}+\xi^{\star}_{0L}\xi_{0L})=\int_{-\infty}^{x}d\hat{x}\,V^{\prime}(\Phi).$
(6.4)

Remarkably, one can show by direct computation that the system of first order
equations (6.1)-(6.2) together with the first order integro-differential eq.
(6.4) imply the second order differential eq. (6.3). So, a solution of the
first-order system of eqs. comprising (6.1)-(6.2) and (6.4) will also be a
solution of the system of eqs. (6.1)-(6.3). This property has also been used
in the preceding sections in order to find the fermion-kink solitons for the
non-vanishing spectrum, i.e. for $E\neq 0$.

Moreover, the identity (3.66) for $E=0$ and upon integration in the whole line
and written for the $\xi_{0}$ spinor components becomes

$\displaystyle\int_{-\infty}^{+\infty}dx[(\xi^{\star}_{0R})^{{}^{\prime}}\xi^{\prime}_{0R}+(\xi^{\star}_{0L})^{{}^{\prime}}\xi^{\prime}_{0L}]-M^{2}\int_{-\infty}^{+\infty}dx\left(\xi^{\star}_{0R}\xi_{0R}+\xi^{\star}_{0L}\xi_{0L}\right)=0.$
(6.5)

The bound state solutions must satisfy the last identity, as we will show
below.

Next, in order to find the zero-mode fermion bound states coupled to kink-type
solitons we resort to the tau function approach applied to the system of first
order equations (6.1)-(6.2) and the first order integro-differential equation
(6.4). Let us assume the following Ansatz

$\displaystyle e^{i\beta\Phi/2}$ $\displaystyle=$ $\displaystyle
e^{i\alpha_{0}}\,\frac{\tau_{0}}{\tau_{1}},$ (6.6) $\displaystyle\xi_{0R}$
$\displaystyle=$ $\displaystyle\frac{\tau_{R}}{\tau_{0}},$ (6.7)
$\displaystyle\xi_{0L}$ $\displaystyle=$
$\displaystyle\frac{\tau_{L}}{\tau_{1}}.$ (6.8)

So, replacing the relationships (6.6)-(6.8) into the equation (6.1) and (6.2),
respectively, one can write the following equations in terms of the tau
functions

$\displaystyle
e^{-2i\alpha_{0}}M\tau_{1}\tau_{L}-\tau_{R}\tau^{\prime}_{0}+\tau_{0}\tau^{\prime}_{R}$
$\displaystyle=$ $\displaystyle 0,$ (6.9) $\displaystyle
e^{2i\alpha_{0}}M\tau_{0}\tau_{R}-\tau_{L}\tau^{\prime}_{1}+\tau_{1}\tau^{\prime}_{L}$
$\displaystyle=$ $\displaystyle 0.$ (6.10)

Let us choose the parameters $\beta_{j}$ in the potential (2.2) to be
$\beta_{j}=\beta\,j,(j=1,2,3)$. We will see that this choice of the
$\beta_{j}$ parameters provides a bound state solution. Likewise, substituting
the tau functions in (6.6)-(6.8) into (6.4) one gets

$\displaystyle-\frac{\beta^{2}\tau_{L}\tau_{R}+2i\tau_{1}\tau^{\prime}_{0}-2i\tau_{0}\tau^{\prime}_{1}}{\beta^{2}\tau_{0}\tau_{1}}+\frac{i\beta}{2}\int^{x}_{-\infty}d\hat{x}\,\left(\frac{e^{2i\alpha_{0}}\,\tau_{0}^{4}-e^{-2i\alpha_{0}}\tau_{1}^{4}}{\tau_{0}^{6}\tau_{1}^{6}}\right)\times$
$\displaystyle\Big{[}3A_{3}e^{4i\alpha_{0}}\tau_{0}^{8}+2A_{2}e^{2i\alpha_{0}}\tau_{0}^{6}\tau_{1}^{2}+(A_{1}+3A_{3})\tau_{0}^{4}\tau_{1}^{4}+2A_{2}e^{-2i\alpha_{0}}\tau_{0}^{2}\tau_{1}^{6}+3A_{3}e^{-4i\alpha_{0}}\tau_{1}^{8}\Big{]}=0.$
(6.11)

Notice that the non-local terms above will vanish for $A_{j}=0,\,j=1,2,3$,
which corresponds to $V=0$. Below, we will search for special tau functions
such that the whole non-local term above, once the $x-$integral is performed,
maintain similar tau function structure as the first local term in (6.11),
i.e. a rational expression of bi-linear terms in the tau function variables.

So, let us consider the following tau functions,

$\displaystyle\tau_{0}$ $\displaystyle=$ $\displaystyle
1+i\,\epsilon\,e^{2\kappa x},\,\,\,\,\,\tau_{1}=1-i\,\epsilon\,e^{2\kappa
x},\,\,\,\,\,\,\,\,\,\epsilon=\pm 1,$ (6.12) $\displaystyle\tau_{L}$
$\displaystyle=$ $\displaystyle\zeta_{1}e^{i\alpha_{1}}e^{\kappa
x},\,\,\,\,\,\,\,\tau_{R}=\zeta_{2}e^{i\alpha_{2}}e^{\kappa x},$ (6.13)

where the parameters $\kappa,\zeta_{a},\alpha_{a}\,(a=1,2)$ are real numbers.

Then, for the above tau functions the fields take the forms

$\displaystyle\Phi$ $\displaystyle=$
$\displaystyle\frac{4}{\beta}\arctan{\Big{[}\epsilon\,e^{2\kappa
x}\Big{]}}+\frac{2\alpha_{o}}{\beta},$ (6.14) $\displaystyle\xi_{0R}$
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\zeta_{2}\Big{[}\mbox{sech}(2\kappa
x)\Big{]}^{1/2}\,e^{i\Theta_{R}},\,\,\,\,\Theta_{R}(x)=-\arctan{\Big{[}\left(\frac{\sinh{[\kappa
x-\frac{1}{2}\epsilon\log{(\tan{\alpha_{2}})}]}}{\cosh{[\kappa
x+\frac{1}{2}\epsilon\log{(\tan{\alpha_{2}})}]}}\right)^{\epsilon}\Big{]}}$
(6.15) $\displaystyle\xi_{0L}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\zeta_{1}\Big{[}\mbox{sech}(2\kappa
x)\Big{]}^{1/2}\,e^{i\Theta_{L}},\,\,\,\,\Theta_{L}(x)=-\arctan{\Big{[}\left(\frac{\sinh{[\kappa
x-\frac{1}{2}\epsilon\log{(\tan{\alpha_{1}})}]}}{\cosh{[\kappa
x+\frac{1}{2}\epsilon\log{(\tan{\alpha_{1}})}]}}\right)^{\epsilon}\Big{]}}$
(6.16)

Notice that $\Phi$ in (6.14) represents a kink type static soliton
interpolating the asymptotic values $\frac{2}{\beta}(\alpha_{0}+\epsilon\pi)$
and $\frac{2\alpha_{0}}{\beta}$. Whereas, the fermion components (6.15)-(6.16)
define a zero-mode fermion bound state, i.e. a solution for $E=0$.

Regarding the parity properties of this solution, the parameters must satisfy
certain relationships implied by (3.13)-(3.15). So, from the kink solution
(6.14) one has that

$\displaystyle\alpha_{o}=\epsilon\frac{\pi}{2},\,\,\,\,\epsilon=\pm 1,$ (6.17)

in order to satisfy the parity condition (3.13). In fact, for
$\epsilon=-1\,(\epsilon=+1)$ one has a kink (antikink) solution. This kink
$\Phi$ resembles the sine-Gordon static soliton; however, the kink in (6.14)
interpolates two vacua defined by

$\displaystyle\beta\Phi(\pm\infty)=\pm\phi_{0},\,\,\,\,\,\,\,\,\phi_{0}\equiv\pi,$
(6.18)

in order to satisfy $\Phi(-x)=-\Phi(x)$ in (3.13). So, following the
definition in (4.54) one gets the topological charges

$\displaystyle Q_{topol}=\pm 1.$ (6.19)

Moreover, the parity conditions (3.15) for the fermion components furnish the
next parameter relationships

$\displaystyle\alpha_{2}$ $\displaystyle=$
$\displaystyle\alpha_{1}+n\pi,\,\,\,\ n\in\leavevmode\hbox{\sf
Z\kern-3.99994ptZ},$ (6.20) $\displaystyle\zeta_{2}$ $\displaystyle=$
$\displaystyle[\frac{\sigma}{\cos{(n\pi)}}]\,\zeta_{1}.$ (6.21)

The tau functions (6.12)-(6.13) will be substituted into the above system of
equations (6.9)-(6.11) and the condition (6.5) must be imposed in order to
find the parameters relationships and study the properties of the solutions
for $\Phi$ and $\xi_{0}$.

Firstly, let us substitute the tau functions (6.12)-(6.13) into the pair of
equations (6.9)-(6.10). So, one gets the next relationships

$\displaystyle\alpha_{o}$ $\displaystyle=$
$\displaystyle-\frac{\pi}{2}(m+n),\,\,\,\,m\in\leavevmode\hbox{\sf
Z\kern-3.99994ptZ}.$ (6.22) $\displaystyle\kappa$ $\displaystyle=$
$\displaystyle-\epsilon M.$ (6.23)

By comparison of the above result (6.22) with (6.17) one gets the next
condition for the integers $m$ and $n$

$\displaystyle m+n=-\epsilon.$ (6.24)

For later convenience let us write the slope of the kink (anti-kink) (6.14) at
$x=0$ for the parameter $\kappa$ defined in (6.23). So, one has

$\displaystyle\mu^{K}\equiv-\epsilon\frac{4M}{\beta}.$ (6.25)

Notice that the identity (6.5) is a condition for the zero-mode solutions and
it does not depend explicitly on the form of the potential $V$. So, one can
derive a general relationship for the parameters by substituting (6.7)-(6.8)
in the form (6.12)-(6.13) into (6.5). The relevant terms of (6.5) become

$\displaystyle\int_{-\infty}^{+\infty}dx[(\xi^{\star}_{0R})^{{}^{\prime}}\xi^{\prime}_{0R}+(\xi^{\star}_{0L})^{{}^{\prime}}\xi^{\prime}_{0L}]$
$\displaystyle=$
$\displaystyle\frac{\pi}{4}\,\kappa\,[\zeta_{1}^{2}+\zeta_{2}^{2}]$ (6.26)
$\displaystyle\int_{-\infty}^{+\infty}dx\left(\xi^{\star}_{0R}\xi_{0R}+\xi^{\star}_{0L}\xi_{0L}\right)$
$\displaystyle=$
$\displaystyle\frac{\pi}{4}\,\frac{1}{\kappa}\,[\zeta_{1}^{2}+\zeta_{2}^{2}]$
(6.27)

Substituting the above quantities into (6.5) one gets again (6.23).

So, the kink and anti-kink of type (6.14) host the zero-mode bound states
(6.15)-(6.16). In order to compare qualitatively the kink (6.14) to the ones
of the section 4.1.1 we have plotted the kink (6.14) in the Fig. 8. Notice
that its asymptotic values are $\pm\pi$ as in (6.18), whereas the kinks of
sec. 4.1.1 interpolate the asymptotic values
$\theta_{0}=\pm\frac{\pi}{3},\pm\frac{2\pi}{3}$, respectively.

Below, in order to verify the equation (6.11) let us consider the two cases:
For vanishing, $V=0$, and non-vanishing, $V\neq 0$, potentials, respectively.

###  6.1 Zero-mode bound states for $V=0$

In this case let us consider $A_{j}=0\rightarrow V=0$. So, by substituting
(6.7)-(6.8) in terms of the tau functions (6.12)-(6.13) into (6.11) with
$A_{j}=0$ one derives the following expression for the parameter $\kappa$

$\displaystyle\kappa=-\epsilon\,\beta^{2}\left(\frac{\zeta_{1}^{2}+\zeta_{2}^{2}}{16}\right).$
(6.28)

Next, let us examine the Majorana zero-modes for $V=0$.

####  6.1.1 Majorana zero-modes for $V=0$

In order to have a Majorana spinor $\xi$ it has been imposed the reality
condition (3.65). It has provided the both cases (3.68)-(3.69). So, we
consider below each case separately.

Case I. The conditions (3.68) provides the next relationships

$\displaystyle\zeta_{2}=\left(\frac{\sin{\alpha_{1}}}{\sin{\alpha_{2}}}\right)\zeta_{1}$
(6.29)
$\displaystyle\alpha_{1}+\alpha_{2}=q\pi,\,\,\,\,\,\,q\in\leavevmode\hbox{\sf
Z\kern-3.99994ptZ}.$ (6.30)

Then, from (6.20) and (6.30) one can get

$\displaystyle\alpha_{1}=\frac{\pi}{2}(q-n),\,\,\,\,\alpha_{2}=\frac{\pi}{2}(q+n).$
(6.31)

Inspecting the relationship (6.29) one must have $\sin{\alpha_{1}}\neq 0$ and
$\sin{\alpha_{2}}\neq 0$; therefore, in (6.31) the integers $n$ and $q$ must
satisfy

$\displaystyle q$ $\displaystyle\equiv$ $\displaystyle s+r+1$ (6.32)
$\displaystyle n$ $\displaystyle\equiv$ $\displaystyle
s-r,\,\,\,\,\,\,\,s,r\in\leavevmode\hbox{\sf Z\kern-3.99994ptZ}.$ (6.33)

In fact, these identities make the quantities $(q\pm n)$ to be odd integers
for any set of integers $s,r\in\leavevmode\hbox{\sf Z\kern-3.99994ptZ}$,
implying $\sin{(\alpha_{j})}\neq 0\,(j=1,2)$.

Case II. The conditions (3.69) provides the next relationships

$\displaystyle\zeta_{2}=\left(\frac{\cos{\alpha_{1}}}{\cos{\alpha_{2}}}\right)\zeta_{1}$
(6.34)
$\displaystyle\alpha_{1}+\alpha_{2}=p\pi,\,\,\,\,\,\,p\in\leavevmode\hbox{\sf
Z\kern-3.99994ptZ}.$ (6.35)

Then, from (6.20) and (6.35) one can get

$\displaystyle\alpha_{1}=\frac{\pi}{2}(p-n),\,\,\,\,\alpha_{2}=\frac{\pi}{2}(p+n).$
(6.36)

Likewise, inspecting the relationship (6.29) one must have
$\cos{\alpha_{1}}\neq 0$ and $\cos{\alpha_{2}}\neq 0$; therefore, in (6.36)
the integers $n$ and $p$ must satisfy

$\displaystyle p$ $\displaystyle\equiv$ $\displaystyle s+r$ (6.37)
$\displaystyle n$ $\displaystyle\equiv$ $\displaystyle
s-r,\,\,\,\,\,\,\,s,r\in\leavevmode\hbox{\sf Z\kern-3.99994ptZ}.$ (6.38)

In fact, these identities make the quantities $(p\pm n)$ to be even integers
for arbitrary set of integers $\\{s,r\\}$, implying $\cos{(\alpha_{j})}\neq
0\,(j=1,2)$.

In the Fig. 13 we plot some examples of the Majorana bound state inside the
kink for the both cases I and II, respectively.

Figure 13: The kink $\Phi$ (blue) and the Majorana zero-mode components
$\xi_{1}$ (green) and $\xi_{2}$ (brown) as functions of $x$. Left (Case I):
For
$\beta=1,\,\alpha_{0}=-\frac{\pi}{2},\,\zeta_{1}=-\zeta_{2}=-\sqrt{8},\kappa=M=1,\alpha_{1}=\pi/2$.
Right (Case II): For
$\beta=1,\,\alpha_{0}=-\frac{\pi}{2},\,\zeta_{1}=\zeta_{2}=-\sqrt{8},\kappa=M=1,\alpha_{1}=\pi$.
The red lines represent the Majorana density $(\xi_{1}^{2}+\xi_{2}^{2})$ vs
$x$, in each case.

####  6.1.2 Fermion back-reaction and soliton-particle(Majorana) duality

Here we follow an analogous approach considered in the sine-Gordon/massive
Thirring mapping in order to exhibit the soliton/particle duality in the
affine Toda model coupled to a Dirac fermion [17, 52]. So, in our case the
zero-mode fermion components can be written as functionals of the scalar field
for the fermion and kink solutions of the type (6.15)-(6.16) and (6.14),
respectively.

Then, from (6.15)-(6.16) and (6.14) one has the next two cases, either for the
soliton or the anti-soliton field.

1\. Majorana coupled to a soliton:
$\alpha_{0}=-\frac{\pi}{2},\,\,\,\epsilon=1$

$\displaystyle\xi_{0L}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\\{\tan{[\frac{1}{4}(\beta\Phi+\pi)]}\\}^{1/2}\zeta_{1}e^{i\alpha_{1}}[1+i\,e^{\frac{i}{2}\beta\Phi}]$
(6.39) $\displaystyle\xi_{0R}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\\{\tan{[\frac{1}{4}(\beta\Phi+\pi)]}\\}^{1/2}\zeta_{1}e^{i\alpha_{2}}[1-i\,e^{-\frac{i}{2}\beta\Phi}],\,\,\,\,\,\,\,\,\,\zeta_{1}^{2}=\frac{8M}{\beta^{2}}.$
(6.40)

2\. Majorana coupled to an anti-soliton:
$\alpha_{0}=\frac{\pi}{2},\,\,\,\epsilon=-1$

$\displaystyle\xi_{0L}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\\{\cot{[\frac{1}{4}(\beta\Phi+\pi)]}\\}^{1/2}\zeta_{1}e^{i\alpha_{1}}[1-i\,e^{\frac{i}{2}\beta\Phi}]$
(6.41) $\displaystyle\xi_{0R}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\\{\cot{[\frac{1}{4}(\beta\Phi+\pi)]}\\}^{1/2}\zeta_{1}e^{i\alpha_{2}}[1+i\,e^{-\frac{i}{2}\beta\Phi}],\,\,\,\,\,\,\,\,\,\,\zeta_{1}^{2}=\frac{8M}{\beta^{2}}.$
(6.42)

Notice that these mappings between the Majorana components and the scalar
field soliton are non-perturbative by construction. Even though they are
classical, they resemble the bosonization formulas at the quantum level [50].
In fact, following the ideas of [17, 52], this mapping can be regarded as the
strong coupling (DSG) $\leftrightarrow$ weak coupling (Majorana fermion model)
sectors of the affine Toda coupled to Majorana fermion.

In the case of the affine Toda model coupled to (charged) Dirac field (eq.
(2.1) with $V=0$) it has been obtained an analogous classical mapping relating
the particle (fermion) to a soliton (boson field) such that they describe the
strong (sine-Gordon) and weak (massive Thirring) coupling sectors of the
model, respectively [17]. At the quantum level, the bosonized Toda model
coupled to Dirac field has been discussed in relation to a confinement
mechanism in the context of the equivalence between the sine-Gordon and
massive Thirring models [16].

Remarkably, one can use the above mappings in order to decouple the strong and
weak coupling sectors of the model at the level of the equations of motion.
First, let us substitute the mappings (6.39)-(6.40) and (6.41)-(6.42) into the
soliton(anti-soliton)-fermion interaction terms of the model in (6.3). So, one
has

$\displaystyle\xi^{\star}_{0R}\xi_{0L}e^{i\beta\Phi}+\xi^{\star}_{0L}\xi_{0R}e^{-i\beta\Phi}$
$\displaystyle\equiv$
$\displaystyle-\epsilon\zeta_{1}\zeta_{2}\cos{(\alpha_{1}-\alpha_{2})}\Big{[}\sin{(\beta\Phi)}+\sin{(2\beta\Phi)}\Big{]},\,\,\,\,\,\epsilon=\mp
1.$ (6.43) $\displaystyle=$
$\displaystyle-\epsilon\zeta_{1}^{2}(-1)^{N}\Big{[}\sin{(\beta\Phi)}+\sin{(2\beta\Phi)}\Big{]},\,\,\,\,\,\,N\equiv
r+n-s,$ (6.44)

where the last relationship holds for the both cases I and II defining the
Majorana fermions in the section 6.1.1. Notice that the r.h.s. of (6.44) is
non-vanishing for any set of values of the integers $r,n,s$ in the both cases
I and II, (6.32)-(6.33) and (6.37)-(6.38), respectively.

Since the r.h.s. of (6.44) is non-vanishing one can conclude that the back-
reaction of the Majorana fermion on the soliton is non-vanishing and depends
on the special values of the fermion parameters
$\zeta_{1}^{2}=\zeta_{2}^{2}\equiv\frac{8M}{\beta^{2}}$ and the non-vanishing
factor $(-1)^{N},\,N\in\leavevmode\hbox{\sf Z\kern-3.99994ptZ}$, determined by
the relationship between the phases $\alpha_{1}$ and $\alpha_{2}$ of the
Majorana fermion components in the cases I and II presented in the sec. 6.1.1.
Recently, for the scattering of the iso-doublet femions on an external sine-
Gordon field and considering a small fermion-soliton coupling, it has been
reported a vanishing of the back-reaction of the Majorana fermion on the
soliton [37].

Therefore, substituting the identity (6.44) into (6.3) (with $V=0$) one has

$\displaystyle-\frac{d^{2}\Phi}{dx^{2}}-8(-1)^{N}\epsilon\frac{M^{2}}{\beta}\Big{[}\sin{(\beta\Phi)}+\sin{(2\beta\Phi)}\Big{]}=0.$
(6.45)

This is the static equation of the double sine-Gordon model (DSG) defined for
the scalar field $\Phi$ with the relevant potential

$\displaystyle
V_{DSG}(\Phi)=8(-1)^{N}\epsilon\frac{M^{2}}{\beta}\Big{[}\cos{(\beta\Phi)}+\frac{1}{2}\cos{(2\beta\Phi)}-\cos{(\frac{2\pi}{3})}-\frac{1}{2}\cos{(\frac{4\pi}{3})}\Big{]}.$
(6.46)

The Fig. 14 shows the plot of this potential with vacua
$\beta\Phi=\frac{2\pi}{3}\,n\,(n\in\leavevmode\hbox{\sf Z\kern-3.99994ptZ})$.
In fact, the kink type solutions of the DSG interpolate these vacuum points.
Notice that the points $\beta\Phi=n\pi,(n\in\leavevmode\hbox{\sf
Z\kern-3.99994ptZ})$ of the potential (6.46) are critical points which do not
support stable soliton type solutions. On the other hand, the soliton (6.14)
hosts the Majorana fermion and it interpolates the vacuum values
$\Phi=\pm\pi$. So, the vacua of the effective potential (6.46) do not support
kink-type solitons of the type (6.14).

In analogy to the mapping between the strong/weak coupling sectors of the
affine Toda model as the soliton/particle correspondence [16] and the
discussion in section 5 of the present paper on the strong coupling sector of
the model (2.1), one can think of the DSG model (6.45) as describing the
strong coupling sector of the model (2.1) with $V=0$ with equations of motion
(6.1)-(6.3) and (6.3) in the zero-mode sector.

It is interesting to analyze the kink-like solutions of the DSG model. The
model (6.45) exhibits a kink solution of the type

$\displaystyle\Phi_{DSG}(x)$ $\displaystyle=$
$\displaystyle\frac{2}{\beta}\arctan{\Big{[}\tan{(\frac{\pi}{3})}\,\,\tanh{(\sqrt{3}Mx)}\Big{]}},$
(6.47)

provided that $\epsilon=1$ and $N$ is an even integer in the potential (6.46).
In fact, this kink interpolates the vacua $\pm 2\pi/3$ and it resembles to one
of the kinks plotted in the Fig. 8 with slope at $x=0$ provided by

$\displaystyle\mu^{DSG}=6\,\frac{M}{\beta}.$ (6.48)

Notice that the kink (6.14) which hosts the Majorana zero-mode interpolates a
different vacua, i.e. $\pm\pi$. However, at the origin $x=0$ one can compute
the ratio between the slopes of the DSG kink (6.48) and the slope of the SG-
type kink in (6.25). So, one has

$\displaystyle\frac{\mu^{DSG}}{\mu^{K}}=\frac{3}{2}.$ (6.49)

Remarkably, in the strong coupling sector of the model the back-reaction
effect of the zero mode (6.15)-(6.16) on the kink (6.14) reflects by shifting
its vacuum points and rising the kink slope according to (6.49).

Moreover, the DSG model (6.45) possesses a bounce-like solution which
interpolates between the vacuum value $\pi$ and the point
$\Phi_{DSG}=\frac{\pi}{2}$ and then it comes back. Since $\pi$ is a false
vacuum position this type of solution is not related to any stable particle in
the quantum theory (see [19] and references therein).

The DSG model has been mapped to various spin models at the quantum level, and
in particular, its relation with the deformed quantum Ashkin-Teller model has
been considered [50]. In these developments the bosonized Majorana fermion
field plays a fundamental role.

Moreover, the sine-Gordon type models play an important role in
superconductors, especially multi-band superconductors including layered high-
temperature superconductors. The vortex fluxes and the associated kink-like
solutions of the double SG appear in the study of multi-band superconductors
in the Ginzburg-Landau theory approach to superconductivity (see e.g. [53, 54]
and references therein).

Next, from the expressions (6.6)-(6.8) and (6.12)-(6.13) in order to achieve
the Majorana fermion sector corresponding to the weak coupling sector of the
model one can get the relationship

$\displaystyle e^{i\beta\Phi}=-\left(\frac{\xi_{0L}}{\xi_{0R}}\right)^{2}.$
(6.50)

So, substituting (6.50) into (6.1)-(6.2) one gets the non-linear field
equations for the Majorana field components

$\displaystyle\frac{d}{dx}\xi_{0R}$ $\displaystyle=$ $\displaystyle\mp
M\frac{\xi_{0R}^{3}}{|\xi_{0L}|^{2}},$ (6.51)
$\displaystyle\frac{d}{dx}\xi_{0L}$ $\displaystyle=$ $\displaystyle\mp
M\frac{\xi_{0L}^{3}}{|\xi_{0R}|^{2}},$ (6.52)

where the $\mp$ signs stand for the case I and II, respectively. This system
of equations can be regarded as describing the weak coupling sector of the
model (6.1)-(6.3).

Moreover, the equation (6.4) for $V=0$ can be written as

$\displaystyle\frac{1}{\beta}\frac{d\Phi}{dx}=-\frac{1}{2}(\xi^{\star}_{0R}\xi_{0R}+\xi^{\star}_{0L}\xi_{0L})$
(6.53)

In fact, the mapping (6.41)-(6.42) [or (6.39)-(6.40)] satisfies the above
relationship. This is reminiscent of the topological and Noether currents
equivalence for the charged fermion [16]. However, in the real field Majorana
fermion case the equivalence (6.53) can be regarded as an equivalence of the
topological charge density to the Majorana state density.

Our result above is analogous to the mapping of the Liouville equation for a
real scalar plus an auxiliary field to a relativistic non-linear and non-local
model for real Majorana type spinor fields in 1+1 dimensions in the context of
the dual $\sigma$ model [55]. Moreover, it is known the equivalence of a
gapless Majorana fermion and the critical point of a real scalar field in
$(1+1)d$. The diccionary being ${\cal
L}_{f}^{(1+1)d}=\bar{\xi}\gamma_{\mu}\partial_{\mu}\xi\leftrightarrow{\cal
L}_{b}^{(1+1)d}=(\partial_{\mu}\phi)^{2}+\phi^{4}$, which holds in the
infrared region (see [56] and references therein). In fact, this relationship
can be understood as a bosonization of $(1+1)d$ free Majorana fermions. In
this context, the Lagrangian formulation and the bosonization of the self-
interacting Majorana fermion system (6.51)-(6.52) deserves a further study.

On the other hand, Majorana bound states in topological superconductors are
pursued due to their potential applications in topological quantum
computation. Since there are other entities as topological bound states in
such systems, a so-called Majorana polarization quantity $P_{M}$ has been
defined in order to discriminate between proper Majorana fermions from other
bound states [57]. Analogously, let us define the quantity

$\displaystyle P_{M}=\frac{{\cal N}}{{\cal D}},$ (6.54)

where

$\displaystyle{\cal
N}\equiv\int_{-\infty}^{+\infty}dx\,2\,\xi_{1}\xi_{2};\,\,\,\,\,\,{\cal
D}\equiv\int_{-\infty}^{+\infty}dx\,[\xi_{1}^{2}+\xi_{2}^{2}].$ (6.55)

For the above soliton (antisoliton) solutions and the host Majoranas one finds

$\displaystyle{\cal
N}=\mp\frac{1}{4}\zeta_{1}^{2}\int_{-\infty}^{+\infty}dx[1+\cos{(\beta\Phi)}],\,\,\,\,{\cal
D}=\frac{1}{4}\zeta_{1}^{2}\int_{-\infty}^{+\infty}dx\cos{(\frac{\beta\Phi}{2})},$
(6.56)

where the $\mp$ signs corresponds to the cases I and II, respectively. Since
$\Phi(\pm\infty)=\pm\pi/\beta$ one has

$\displaystyle P_{M}=\mp\frac{\pi}{4}\approx\mp 0.7854.$ (6.57)

Figure 14: The plot of the DSG potential $V(\Phi)$ in (6.46). Notice the
minima for $\beta\Phi=\frac{2\pi}{3}\,n\,(n\in\leavevmode\hbox{\sf
Z\kern-3.99994ptZ})$. The values $\Phi=\pm\pi$ correspond to unstable critical
points.

In order to compare these values of $P_{M}$ with the behavior of the Dirac
spinor obtained in the above sections let us define the quantity

$\displaystyle P_{D}=\frac{{\cal N}_{D}}{{\cal D}_{D}},$ (6.58)

where

$\displaystyle{\cal
N}_{D}\equiv\int_{-\infty}^{+\infty}dx\,\,(\xi_{1}\xi_{2}-\xi_{3}\xi_{4});\,\,\,\,\,\,{\cal
D}_{D}\equiv\frac{1}{2}\int_{-\infty}^{+\infty}dx\,[\xi_{1}^{2}+\xi_{2}^{2}+\xi_{3}^{2}+\xi_{4}^{2}],$
(6.59)

where the Dirac spinor components have been considered.

Notice that the definition (6.58)-(6.59) reduces to the relationship
(6.54)-(6.55) once the identification (3.69) (or (3.68) with ${\cal
N}_{D}\rightarrow-{\cal N}_{D}$ ) is taken into account in order to define a
Majorana fermion. Computing the above integrals for the explicit Dirac spinor
solutions and setting $\theta_{1}=\theta_{2}=\theta_{o}$ one has

$\displaystyle
P_{D}=-\frac{(\sin{\theta_{o}})^{2}}{\theta_{o}}\,\,\,\rightarrow\,\,\,\,|P_{D}|<0.725.$
(6.60)

By comparison of the relations (6.57) and (6.60) it can be shown, for any
value of the parameter $\theta_{o}$, that the quantity $|P_{D}|$ is less than
the Majorana polarization quantity $|P_{M}|$ (i.e. $|P_{D}|<|P_{M}|$).

###  6.2 Zero-mode bound states for $V\neq 0$

In this case one can derive parameter relationships by substituting
(6.7)-(6.8) in terms of the tau functions (6.12)-(6.13) into (6.11). So, one
has the relationships

$\displaystyle\kappa$ $\displaystyle=$
$\displaystyle-\frac{1}{32}\epsilon\,\Big{\\{}\beta^{2}\left(\zeta_{1}^{2}+\zeta_{2}^{2}\right)+\beta^{3/2}\sqrt{\beta\left(\zeta_{1}^{2}+\zeta_{2}^{2}\right)^{2}+256A_{1}}\Big{\\}},\,\,\,\beta>0,$
(6.61) $\displaystyle A_{1}$ $\displaystyle=$
$\displaystyle\epsilon\frac{\kappa}{4\beta^{3}}\Big{[}16\epsilon\kappa+\beta^{2}(\zeta_{1}^{2}+\zeta_{2}^{2})\Big{]},$
(6.62) $\displaystyle A_{2}$ $\displaystyle=$ $\displaystyle A_{3}=0.$ (6.63)

Notice that for $A_{1}=0$ the relationship (6.61) for $\kappa$ reduces to the
previous one in (6.28). Taking into account (6.23) and (6.61) one can get the
expression for $A_{1}$

$\displaystyle
A_{1}=\frac{M}{4\beta^{3}}\Big{[}16M-\beta^{2}(\zeta_{1}^{2}+\zeta_{2}^{2})\Big{]}.$
(6.64)

The equation (6.11) is satisfied for the values in (6.61)-(6.63) without any
additional restrictions for the parameters, except that the phase angle
$\alpha_{0}$ must be $\epsilon\frac{\pi}{2}$ as in (6.17). So, one has the
potential

$\displaystyle V=A_{1}\cos{(\beta\Phi)}.$ (6.65)

So, the kinks of type (6.14) host the zero-mode bound states (6.15)-(6.16) in
the model (2.1) for the sine-Gordon potential $V$ provided in (6.65).
Recently, for an external (without back-reaction) sine-Gordon kink which
interacts weakly with fermionic iso-dublets through a Yukawa interaction term
it has been discussed the appearance of Majorana zero-modes, which do not
exert back-reactions on the kink in this approximation, i.e. the fermion-kink
interaction term in the Majorana sector vanishes [37]. Moreover, the non-
vanishing back-reaction of the Majorana zero-modes on the kink of our model
with the sine-Gordon potential (6.65) would be discussed in analogy to the
preceding sections 6.1.1 and 6.1.2.

##  7 Numerical simulations

We have shown above that the system of equations (3.4)-(3.8) can be solved
analytically and obtained various types of kinks which interpolate the vacua
$\Phi(\pm\infty)=\pm\frac{\pi}{3},\pm\frac{2\pi}{3},\pm\pi$, respectively. The
analytic kinks and their relevant spinor bound states, as well as the
associated eigenvalues $E$ have been found. However, it is desirable to get
analogous solutions, as well as for the eigenvalues $E$, such that the
parameter $\theta_{0}$ ($\Phi(\pm\infty)\equiv\pm\theta_{0}$) belongs to a
continuous interval. A general analytical method is not known for arbitrary
values of that parameter. Below we will use the so-called relaxation method to
numerically solve the system of equations; this method proceeds to find the
solution by starting with a guess solution and searching for the true
solution. In order to achieve a fast convergence, it is needed a good initial
trial guess solution satisfying the same boundary conditions for the kink and
the bound states as the real solutions. Thus, as suitable guess initial
solutions we will choose the form of the analytic kink and fermionic bound
states in (4.50) and (4.51), provided that certain values of the parameters
are assumed.

Our aim is to obtain numerical solutions for the static version of the system
of eqs. (3.4)-(3.8) corresponding to the potential of type $V_{1}(\Phi)$ in
(3.4), such that nontrivial topological configurations can be achieved by the
field $\Phi$, interacting with the bound states formed by the fermion field
$\psi$. It is expected that the spinor bilinears $\bar{\psi}\psi$ and
$\bar{\psi}\gamma_{5}\psi$ will develop some localized bound states when the
scalar field $\Phi$ becomes a kink type soliton. In the appendix F we will
consider the static version of the eqs. of motion (3.4)-(3.8) in order to
implement their relevant discretizations .

We will search for static solutions for $\Phi$ representing a kink type
soliton satisfying the parity and boundary conditions (3.13)-(3.14) and
solving the system (3.4)-(3.8) such that the bound states are invariant under
the parity symmetry (3.10) and (3.15). Under these assumptions we will search
for solutions to the system to be eigen-functions of the parity operator
${\cal P}_{x}$, with $\sigma=\pm 1$ eigenvalues.

Figure 15:  $\frac{E}{M}$ vs $\theta_{0}$. The parities of the bound states
are indicated as $(\pm)$. Notice the BIC states in the left column such that
$|\frac{E}{M}|>1$ in some regions of $\theta_{0}$.

In the appendix F we present the relevant developments of the numerical
simulation in the context of the relaxation method. At the final step of each
iteration, following the algorithm in this method, one gets the relevant
eingenvalue $E$ defined as a constant independent of the variable $x$.

In the Figs.15-18 we present our numerical results. In the Fig.15 it is shown
the numerical plot $E$ vs $\theta_{0}$ for the interval $[-\pi,\pi]$. In this
case the initial trial kink and bound state functions have been defined for
the analytical solutions (4.51) and (4.50) for
$\kappa=1,M=1,\beta=1,\theta_{1}=1.548,\theta_{2}=1.552,A_{3}=10^{5},z=4.0,\rho_{1}=0.5,\rho_{2}=2.19$.
So, it is numerically simulated the behavior of $E$ as the angle $\theta_{0}$
varies. For a given $\theta_{0}$ one gets the two energies $\pm E$
corresponding to the conjugation symmetry and to the positive and negative
parities. Notice that the plots in the Figs.15 approximately approach the
values $E=\pm 0.33,\pm 0.83$, within numerical accuracy. In fact, these values
correspond to the angle $\theta_{0}=\pi/3$ and $\kappa/M=1$ of the analytical
results presented in the Figs.4 and 5. In addition, the both Figs.15 show
qualitatively that the eigenvalue $E$ as a function of $\theta_{0}$ is an even
function, i.e. $E(-\theta_{0})=E(\theta_{0})$. Notice that the analytical
result for $E$ (4.29) exhibits this symmetry.

Moreover, the two plots in the right column of Fig.15 show that the spectra
belongs entirely to the in-gap states, i.e. $|\frac{E}{M}|<1$. Remarkably, the
left column of Fig. 15 exhibits the appearance of the bound states in the
continuum (BIC) such that $|\frac{E}{M}|>1$ for certain continuous regions of
the parameter; i.e. $\frac{\pi}{3}<\theta_{0}<\frac{2\pi}{3}$ and
$-\frac{2\pi}{3}<\theta_{0}<-\frac{\pi}{3}$.

Next, in the Fig.16 the initial trial kink and bound state functions in (4.51)
and (4.50) for
$\theta_{0}=\frac{\pi}{3},M=1,\beta=1,\theta_{1}=1.548,\theta_{2}=1.552,A_{3}=10^{5},z=4.0,\rho_{1}=0.5,\rho_{2}=2.19$
are used to simulate the behavior of $E$ under the variation of $\kappa$ in
the interval $\frac{\kappa}{M}\in[-1.5\,,\,1.5]$ for the fixed value of
$\theta_{0}=\frac{\pi}{3}$. Notice that all the states are confined to the in-
gap region of the spectra. The plots approach the values $E=\pm 0.33$ for
$\kappa/M=1$, within numerical accuracy. In fact, these values correspond to
the angle $\theta_{0}=\pi/3$ and $\kappa/M=1$ of the analytical results
presented in the Figs.4 and 5. In addition, the Figs.16 show qualitatively
that the eigenvalue $E$ changes sign as $E\rightarrow-E$ under the
transformation: $\kappa\rightarrow-\kappa$ and $\sigma\rightarrow-\sigma$.
This symmetry also appears in the analytical results, as stated in the eq.
(4.49).

Next, in the Fig.17 we present $E$ vs $\theta_{0}$ for the interval
$[-\pi,\pi]$ for another set of initial trial parameters. In fact, the initial
trial kink and bound state functions have been assumed for the analytical
solutions (4.51) and (4.50) as
$\kappa=-1,M=1,\beta=1,\theta_{1}=1.1,\theta_{2}=0.9,A_{3}=10^{5},z=4.0,\rho_{1}=1.1,\rho_{2}=0.9$.
Notice that the plots in the right column of the Fig.17 approximately approach
the value $E=\pm 1.1$. In fact, this value corresponds to the angle
$\theta_{0}=-2\pi/3$ and $\kappa/M=-1$ of the analytical results presented in
the Fig.6.

Finally, in the Figs.18 the initial trial kink and bound state functions used
in the Fig.17 are used to simulate the behavior of $E$ under the variation of
$\kappa$ in the interval $\frac{\kappa}{M}\in[-1.5\,,\,1.5]$ for the fixed
value of $\theta_{0}=\frac{2\pi}{3}$. Notice that all the states are confined
to the in-gap region of the spectra in this case. The plots approach the
values $E=\mp 0.06$ for $\kappa/M=\pm 1$, within numerical accuracy. In fact,
these values correspond to the angle $\theta_{0}=2\pi/3$ and $\kappa/M=\mp 1$
of the analytical results presented in the Figs.6 and 7. In addition, the
Figs.18 show qualitatively that the eigenvalue $E$ changes sign as
$E\rightarrow-E$ under the transformation: $\kappa\rightarrow-\kappa$ and
$\sigma\rightarrow-\sigma$. This symmetry also appears in the analytical
results, as stated in the eq. (4.49).

So, the exact values of the set $\\{E,\theta_{0},\sigma,M\neq 0\\}$ of our
analytical results provided in the tables and figures 4, 5, 6, and 7 for the
set $\theta_{0}=\pm\frac{\pi}{3},\pm\frac{2\pi}{3}$, can be approached
numerically provided that a convenient set of initial trial parameters are
assumed for each simulation.

Figure 16:  $\frac{E}{M}$ vs $\frac{\kappa}{M}$ for $\theta_{0}=\pi/3$. The
parities of the bound states are indicated as $(\pm)$. These plots show in-gap
states ($|\frac{E}{M}|<1$). Figure 17:  $\frac{E}{M}$ vs $\theta_{0}$. The
parities of the bound states are indicated as $(\pm)$. Notice the BIC states
in the both Figs. such that $|\frac{E}{M}|>1$ for certain discrete and
continuous regions of $\theta_{0}$. Figure 18:  $\frac{E}{M}$ vs
$\frac{\kappa}{M}$ for $\theta_{0}=2\pi/3$. The parities of the bound states
are indicated as $(\pm)$. These plots show in-gap states ($|\frac{E}{M}|<1$).

Moreover, inspecting the Figs. 15-18 one can see the appearance of certain
zero-modes for some values of the set of parameters
$\\{\theta_{0},\kappa,\sigma,\,M\neq 0\\}$. In fact, the left columns of Figs.
15-16 and the both Figs. 17-18 exhibit the presence of zero-modes. They are
solitons with topological charge $\frac{\theta_{0}}{\pi}$ hosting exotic bound
state fermions. It would be interesting to examine numerically their
properties as compared to the ones of the zero-mode Majorana fermions; e.g.
the degree of divergence from the Majorana component relationships
(3.68)-(3.69) and their relevant polarization $P_{D}$ as compared to the ones
for the Majorana fermion $P_{M}$ as in (6.58) and (6.54), respectively. So,
the numerical solutions satisfying those properties deserve a further
examination and we postpone it for a future work.

##  8 Discussion

We have considered a modified Toda model coupled to matter providing the
analytical and numerical solutions as bound states of the fermion and solitons
of the scalar field. A deformation of the ATM model was performed by adding a
multi-frequency scalar potential. The solutions related to a special set of