# The quantum fluctuation in an inhomogeneous background and its influence to
the phase transition in a finite volume system

Xiaogang Li1 Song Shu1<EMAIL_ADDRESS>Jia-Rong Li2 1\. Department of Physics
and Electronic Science, Hubei University, Wuhan, Hubei 430062, China  2\. Key
Laboratory of Quark and Lepton Physics (MOE) and Institute of Particle
Physics, Central China Normal University, Wuhan, Hubei 430079, China

###### Abstract

We have studied the grand potential and phase transitions of an inhomogeneous
finite volume spherical quark system. First the finite volume effect are
considered by applying the multiple reflection expansion method which is an
approximation for the density of states of the momentum in the grand potential
of the finite size system. Then the density of states of momentum is further
modified by the thermal fluctuations in the thermal system with the
inhomogeneous field background. The modification of the density of states is
calculated by the scattering phase shift from the quark Dirac equation in the
inhomogeneous field background. By the numerical calculation we have discussed
how the phase transition is changed by varying the configuration of the
inhomogeneous field background and show that the strong first order phase
transition could be changed to a crossover as a result.

##  I Introduction

The physics of the phase transition in small quantum system is very intricate,
especially in a tiny droplet of quark/nuclear matter created by the high
energy heavy ion collisionsLuo:2020pef ; Stephanov:1998dy . QCD itself and its
nonperturbative nature make it very hard to solve the problem of the phase
transition from the first principleShuryak:2021vnj . Through many years of
efforts made by QCD lattice calculation and effective model studies it is
known that at high temperatures and low baryon densities the transition from
nuclear matter to quark matter is a smooth crossover rather than a first order
phase transition. At relative high densities the phase diagram is still not
clear or not confirmedFischer:2018sdj ; Ding:2015ona ; Ding:2020rtq . The fire
ball created in the heavy ion collisions is very small with an estimated
radius about $2\sim 10\rm{fm}$. In theoretical study one should consider
finite size effectGraef:2012sh ; Klein:2017shl ; Fister:2015eca ;
Palhares:2009tf . Besides that the system may rather be inhomogeneous. There
are many discussions about the inhomogeneous phases in chiral phase
transitions or deconfinement phase transitionsBuballa:2014tba ;
Carignano:2014jla ; Nickel:2009wj ; Lakaschus:2020caq ; Fukushima:2020cmk ;
Fukushima:2013rx ; Kojo:2011cn , like chiral density wave(CDW) phases and
quarkyonic phases. However the quantum effect is nontrivial and involved to
deal with in the inhomogeneous systems. As for the CDW phase, a simple plane
wave chiral condensate is assumed to simplify the problem of solving the
complicated quark Dirac equations. The quantum effect of quark Dirac sea in
inhomogeneous background is important in studying the phase transition in
quark/nuclear matter. In this work we want to discuss the phase transition in
an inhomogeneous spherical quark system which has a finite volume. The finite
volume effect can be introduced in different ways. It could be added through
MIT boundary conditions or in multiple reflection expansion(MRE)
methodXu:2020loz ; Liu:2020elq ; Zhao:2018nqa ; Kiriyama:2002xy . The original
MRE method was derived for massless particlesBalian:1970ap . However it has
been further extended to calculate the thermodynamic quantities for quark
system with nonzero massMadsen:1993prl ; Gilson:1993prl . It is suitable for
studying spherical system with a finite volume. In this work as we discuss a
finite spherical system, the MRE scheme will be used. It will have a
correction to the density of states of the momentum in the momentum
integration in the grand potential. It should be emphasized that in previous
studies of the grand potential of the thermal system the mean field
approximation(MFA) in the homogeneous case has been often used to study the
phase transition. However in this study we argue that the quantum fluctuation
over the spatially inhomogeneous background will have very nontrivial
influence to the phase transition. The inhomogeneous background means that at
MFA the expectation value of the field has a spatial variation. As for a
thermal quark system in the inhomogeneous background the usually density
states of the momentum in the stand MFA grand potential will receive a
correction which can be derived by the scattering phase shift by solving the
quark Dirac equation in the inhomogeneous field background. The method of
evaluating the quantum fluctuations in inhomogeneous field background through
the scattering phase shift can be traced back to the work of Schwinger in the
study of the energy of the electron’s quantum fluctuations in certain
nontrivial configuration of electromagnetic fieldsSchwinger:1954pr . In later
years Jaffe and his collaborators further completed the renormalization
program and extended the method for calculating the quantum fluctuations in
inhomogeneous field background in electroweak systemsJaffe:1998plb ;
Jaffe:2002npb . The scattering phase shift method is practical to calculate
quantum fluctuations in quantum field system with spherically inhomogeneous
field background.

In this paper we will use Friedberg-Lee(FL) model to address these problems
and to illustrate the calculation scheme and its nontrivial influence to the
phase transition. In previous study of the FL model at MFA it indicates a
first order phase transition at finite temperatures and densitiesMao:2008prc ;
Gao:1992prd . The first order phase transition is relevant to the neutron star
mergers and cosmological electroweak phase transitionsPaschalidis:2017qmb ;
Ferreira:2020evu ; Trodden:1998ym . In a recent study the bubble dynamics is
studied in the strong first order phase transition in FL modelZhou:2020bzk .
In relativistic heavy ion collisions the crossover is the main feature at high
temperatures and low densities. In our previous study by a qualitative
analysis of the tunnelling amplitude of the two vacuums in FL model we
indicate that the first order phase transition may be turned into a crossover
by the quantum tunnellingShu:2011ij . In many chiral effective models the
study shows that the quantum fluctuations may weaken the first order phase
transition and strengthen the feature of the crossoverSkokov:2010sf ;
Herbst:2013ail ; Brandes:2021pti . However in these studies the quantum
fluctuations are treated in a homogeneous background. In this work we will
show in FL model how the quantum fluctuations in an inhomogeneous background
influence the phase transition in a finite volume. The organization of this
paper is as follows: in section II the FL model is introduced and the grand
potential is obtained at the MFA level. In section III it is illustrated how
the grand potential is modified in a finite size quark system under the
inhomogeneous background. And the order parameter of the phase transition is
derived accordingly. In section V, we present the numerical results and
discuss how the phase transition is changed by varying the spatial
configuration of the inhomogeneous field background. The last section is the
summary and outlook.

##  II The model and grand potential at MFA

FL model was originally introduced to describe the properties of hadrons. It
incorporates different bag models, like MIT and SLAC models, into one
effective model of quantum field theory which is renormalizableLee:1977prd1 ;
Lee:1977prd2 . We start from the Lagrangian of the FL model as

$\displaystyle{\cal
L}=\bar{\psi}(i\gamma_{\mu}\partial^{\mu}-g\sigma)\psi+\frac{1}{2}\partial_{\mu}\sigma\partial^{\mu}\sigma-U(\sigma),$
(1)

where

$\displaystyle
U(\sigma)=\frac{1}{2!}a\sigma^{2}+\frac{1}{3!}b\sigma^{3}+\frac{1}{4!}c\sigma^{4}+B.$
(2)

$\psi$ is the quark field and $\sigma$ is the phenomenological scalar field
which can roughly describe the nonperturbative long range effect of the gluon
field. $a,b,c$, and $g$ are the model parameters which are fixed by fitting
the general properties of nucleons. $B$ is the bag constant. In the MFA the
sigma field is treated as the classical mean field as
$\sigma(\vec{r},t)\rightarrow\bar{\sigma}$. The classical potential
$U(\bar{\sigma})$ generally has two minima which are corresponding to two
vacuums. One is the perturbative vacuum at $\bar{\sigma}=0$ which is a local
minimum. The other is the physical vacuum at $\bar{\sigma}=\sigma_{v}$ which
is a global minimum. At the background of the physical vacuum if there are
three valence quarks, they will rather be localized than propagate freely in
space with heavy masses. This is the soliton solution in this model
corresponding to a nucleon in a physical interpretation. In a simple
understanding one can regard this bound state solution as a confinement state
of three valence quarks through a bag. As long as the vacuum is kept in the
physical vacuum, the physical ground state of the three valence quarks will
keep confined as a soliton. So in this model in order to produce deconfinement
one has to change the vacuum. The possible way is to heat the vacuum. In FL
model one can introduce thermal grand potential through standard method of
finite temperature field theory. At finite temperature the partition function
of the system is

$Z=\int[id\psi^{\dagger}][d\psi][d\sigma]\exp\left[\int_{0}^{1/T}d\tau\int
d^{3}{\bf r}{\cal L}\right].$ (3)

At MFA we can write the grand potential as

$\Omega(\bar{\sigma},T)=-\frac{T\log
Z}{V}=\Omega_{\sigma}+\Omega_{\textmd{kin}},$ (4)

where the sigma field part of the grand potential has the following form

$\Omega_{\sigma}=\frac{1}{V}\int
d^{3}\textbf{r}\left[\frac{1}{2}(\nabla\bar{\sigma})^{2}+U(\bar{\sigma})\right],$
(5)

and $\Omega_{\textmd{kin}}$ is the kinetic part of the potential which can be
written as

$\Omega_{\textmd{kin}}=-\frac{T}{V}\textmd{Trlog}\left[\frac{1}{T}(i\partial_{0}+i\gamma^{0}\vec{\gamma}\cdot\vec{\nabla}-\gamma^{0}g\bar{\sigma})\right].$
(6)

In the homogeneous background field of sigma the grand potential can be
derived as

$\Omega_{\textmd{MFA}}(\bar{\sigma},T)=U(\bar{\sigma})-4N_{f}N_{c}T\int
dk\rho_{\textmd{hom}}(k)\log(1+e^{-E_{k}/T}),$ (7)

where $T$ is temperature, $V$ is the volume and
$E_{k}=\sqrt{k^{2}+g^{2}\bar{\sigma}^{2}}$. $\rho_{\textmd{hom}}$ is density
of states of momentum for homogeneous case which is

$\rho_{\textmd{hom}}(k)=\frac{1}{V}\frac{dN}{dk}=\frac{k^{2}}{2\pi^{2}}.$ (8)

One should notice that at $T=0$ the grand potential will be
$\Omega(\bar{\sigma},T=0)=U(\bar{\sigma})$. The FL model at finite
temperatures has been extensively studied in the literature. The general
results is that when the temperature is increased, the grand potential will
change its configuration. As a result the physical vacuum at
$\bar{\sigma}=\sigma_{v}$ is lifted up and becomes the local minimum, while
the perturbative vacuum at $\bar{\sigma}=0$ turns into a global minimum which
is the true vacuum. In the new thermal vacuum the valence quarks will
propagate freely in space instead of being bounded, which means the system is
deconfined. The deconfinement phase transition in FL model is a strong first
order phase transition.

##  III The finite size thermal quark system under the inhomogeneous
background

In the above discussion the thermal quark system is in the homogeneous case
without boundaries. However the thermal quark/nuclear matter created in the
real heavy ion collisions is always finite size system and rather
inhomogeneous. Theoretically one has to consider the corrections to the grand
potential of MFA. Though it is not easy, there are still different ways to
work out the corrections. In this section we will make use of MRE method to
give the finite size corrections to the grand potential, while the corrections
from inhomogeneous background will be obtained through scattering phase shift
method.

First we discuss calculation of of the $logZ$ in an inhomogeneous background.
For a system in a spherical box with a radius $R$ the integration over the
momentum in $logZ$ will be changed into a summation over the discrete momentum
modes. The general form of $logZ$ by the summation over the discretized
momentum modes can be written as

$logZ=\sum_{nlm}f({\bf k}_{nlm},T),$ (9)

where $n,l,m$ are the radial, angular momentum and projection of angular
momentum quantum numbers. Furthermore we consider a spherically symmetric
inhomogeneous background and asymptotically the eigenfunctions of the
Hamiltonian of the system will be given by trigonometric functions with
argument $kr+\delta_{l}(k)$. If we impose the condition that the
eigenfunctions vanish at the boundary, it is found that
$kR+\delta_{l}(k)=n\pi$ where $\delta_{l}(k)$ is the phase shift in the wave
function at the boundary generated by the scattering of the inhomogeneous
background which is in contrast to that of the free case in the homogeneous
background. When the radius $R$ of the system is large enough, the states
become infinitesimally close and we haveSchwinger:1954pr ; Jaffe:1998plb

$\frac{dn}{dk}=\frac{R}{\pi}+\frac{1}{\pi}\frac{d\delta_{l}(k)}{dk}.$ (10)

The density of states in momentum space can be defined as $dn/dk$ divided by
the volume $V$. The first term in the right hand side of the equation (10) is
related to the density of states at homogeneous case. When the system volume
goes to infinity this homogeneous part will be identical to the density of
states $\rho_{hom}(k)$ as in equation (8). Now if we assume that the volume is
sufficiently large at a microscopic level but finite at a macroscopic level,
then by applying the MRE method the density states of the large but finite
volume spherical system in homogeneous background will have the following
formKiriyama:2002xy ; Madsen:1993prl

$\rho_{\textmd{MRE}}(k)=\rho_{\textmd{hom}}(k)+\Delta\rho_{\textmd{MRE}}(k),$
(11)

with

$\Delta\rho_{\textmd{MRE}}(k)=\frac{3k}{R}f_{S}(\frac{k}{m})+\frac{6}{R^{2}}f_{C}(\frac{k}{m}),$
(12)

in which $R$ is the radius of the sphere and $m$ is the mass. In the right
hand side the first and second terms are the surface and the curvature
contributions with the functions $f_{S}(k/m)$ and $f_{C}(k/m)$ as

$f_{S}(\frac{k}{m})=-\frac{1}{8\pi}\left(1-\frac{2}{\pi}\arctan\frac{k}{m}\right),$
(13)
$f_{C}(\frac{k}{m})=\frac{1}{12\pi^{2}}\left[1-\frac{3k}{2m}\left(\frac{\pi}{2}-\arctan\frac{k}{m}\right)\right].$
(14)

For the homogeneous case the second term in the density of sates of the
equation (10) has been neglected. Therefore within the framework of MRE the
grand potential of a finite size spherical system in FL model can be written
as

$\Omega_{\textmd{MRE}}(\bar{\sigma},T)=U(\bar{\sigma})-4N_{f}N_{c}T\int
dk\rho_{\textmd{MRE}}(k)\log(1+e^{-E_{k}/T}),$ (15)

In the next we take into account the inhomogeneous field background. Then the
second term in equation (10) could not be neglected. If we consider a spatial
spherically symmetric sigma field background
$\bar{\sigma}(\vec{r})=\bar{\sigma}(r)$ the density of states of momentum will
receive a further correction derived from the second term in equation (10) of
which the result is

$\Delta\rho_{\textmd{inho}}(k)=\frac{3}{4\pi^{2}R^{3}}\sum\limits_{l}(2l+1)\frac{d\delta_{l}(k)}{dk},$
(16)

where $l$ is the quantum number of angular momentum. In FL model the density
of states of momentum is changed due to the thermal quarks scattering in the
inhomogeneous field background. $\delta_{l}(k)$ is the scattering phase shift
of quarks in the inhomogeneous background field of sigma $\bar{\sigma}(r)$.
The scattering phase shift can be determined by the partial wave scattering
method from the following stationary Dirac equation

$\left[-i\gamma_{0}\vec{\gamma}\cdot\vec{\nabla}+\gamma_{0}g\bar{\sigma}(r)\right]\psi=E\psi.$
(17)

The background field $\bar{\sigma}(r)$ can be decomposed into two parts as the
homogeneous vacuum part and the inhomogeneous spatial part as in the following

$\bar{\sigma}(r)=\bar{\sigma}+\tilde{\sigma}(r),$ (18)

By this decomposition we need to ensure that when $r\rightarrow\infty$, the
background field $\bar{\sigma}(r)$ goes to its vacuum value $\bar{\sigma}$,
meanwhile $\tilde{\sigma}(r)\rightarrow 0$. In this work the
$\tilde{\sigma}(r)$ is set to be a Woods-Saxon form like spatial configuration
which can be written as

$\tilde{\sigma}(r)=-\frac{\alpha}{1+e^{\beta(r-r_{0})}},$ (19)

where $\alpha$, $\beta$ and $r_{0}$ are the adjustable parameters to change
the spatial configuration of the background field $\tilde{\sigma}(r)$. In this
assumption the calculation of the scattering phase shift from equation (17)
has been thoroughly discussed in our previous work Shu:2016rrr ; Shu:2020hdz
where the scattering phase shift has been derived and regularized. The whole
calculation of the renormalization for the phase shift is very nontrivial,
however at MFA level the zero-point energy will be neglected and we only
preserve the regularized phase shift in this work. In order to make the
calculation self-contained in this work, a brief derivation of the phase shift
is presented in the Appendix A. The final result is the regularized phase
shift $\bar{\delta}_{l}(k)$ which can be numerically evaluated. By
substituting the regularized scattering phase shift into equation (16) one can
further determine the correction of the density of states due to the
inhomogeneous background.

Compared to that of the homogeneous system the density of states for a finite
volume inhomogeneous system will have the correction term of
$\Delta\rho_{\textmd{inho}}$. So the final result of the density of states for
the finite size inhomogeneous system is

$\rho(k)=\rho_{\textmd{hom}}(k)+\Delta\rho_{\textmd{MRE}}(k)+\Delta\rho_{\textmd{inho}}(k)=\frac{k^{2}}{2\pi^{2}}+\frac{3k}{R}f_{S}(\frac{k}{m})+\frac{6}{R^{2}}f_{C}(\frac{k}{m})+\frac{3}{4\pi^{2}R^{3}}\sum\limits_{l}(2l+1)\frac{d\delta_{l}(k)}{dk},$
(20)

Then the kinetic part of the grand potential for an inhomogeneous sphere can
be written as

$\Omega_{\textmd{kin}}=-4N_{f}N_{c}T\int dk\rho(k)\log(1+e^{-E_{k}/T}),$ (21)

where $E_{k}=\sqrt{k^{2}+g^{2}\bar{\sigma}^{2}}$. Besides the kinetic part the
sigma field part of the grand potential for the inhomogeneous sphere has the
following form

$\Omega_{\sigma}=\frac{4\pi}{V}\int
drr^{2}\left[\frac{1}{2}\left(\frac{d\bar{\sigma}(r)}{dr}\right)^{2}+U(\bar{\sigma}(r))\right].$
(22)

According to the decomposition of the sigma field the homogeneous part can be
separated out from the integration and the result is

$\Omega_{\sigma}=U(\bar{\sigma})+\frac{4\pi}{V}\int
drr^{2}\left[\frac{1}{2}\left(\frac{d\tilde{\sigma}(r)}{dr}\right)^{2}+\tilde{U}(\tilde{\sigma}(r))\right],$
(23)

where the potential $\tilde{U}(\tilde{\sigma})$ satisfies that
$\tilde{U}(\tilde{\sigma})\rightarrow 0$ when $r\rightarrow\infty$ and
$\tilde{\sigma}\rightarrow 0$ and its form is

$\tilde{U}(\tilde{\sigma}(r))=\frac{1}{2!}(a+b\sigma_{v0}+\frac{1}{2}c\sigma_{v0}^{2})\tilde{\sigma}(r)^{2}+\frac{1}{3!}(b+c\sigma_{v0})\tilde{\sigma}(r)^{3}+\frac{1}{4!}c\tilde{\sigma}(r)^{4},$
(24)

in which $\sigma_{v0}$ is fixed by the absolute minimum of the
$U(\bar{\sigma})$ as

$\sigma_{v0}=\frac{3|b|}{2c}\left[1+\left(1-\frac{8ac}{3b^{2}}\right)^{1/2}\right].$
(25)

The bag constant $B$ in $U(\bar{\sigma})$ is thus fixed as

$B=-(\frac{a}{2!}\sigma_{v0}^{2}+\frac{b}{3!}\sigma_{v0}^{3}+\frac{c}{4!}\sigma_{v0}^{4}).$
(26)

Finally we obtain the grand potential of finite volume spherical inhomogeneous
system as

$\Omega_{\textmd{inhom}}(\bar{\sigma},\tilde{\sigma},T)=U(\bar{\sigma})+\frac{4\pi}{V}\int
drr^{2}\left[\frac{1}{2}\left(\frac{d\tilde{\sigma}(r)}{dr}\right)^{2}+\tilde{U}(\tilde{\sigma}(r))\right]-4N_{f}N_{c}T\int
dk\rho(k)\log(1+e^{-E_{k}/T}).$ (27)

One should notice that the calculation of $\rho(k)$ is very nontrivial. In
principle it should vary with the mass $m=g\bar{\sigma}$, in which case the
calculation will be very involved. As an approximation we assume that
$\rho(k)$ is not sensitive to the variation of the quark mass. As a result the
quark mass in $\rho(k)$ is taken as a fixed value at $m=g\sigma_{v0}$ and the
form of $\rho(k)$ is only dependent on the spatial configuration of the sigma
background field $\tilde{\sigma}(r)$ and the system size $R$. Then the
homogeneous mean field $\bar{\sigma}$ in the grand potential could be treated
as a variable. At a certain temperature it can be determined by the
thermodynamic equilibrium condition as

$\left.\frac{\partial\Omega_{\textmd{inhom}}(\bar{\sigma},\tilde{\sigma},T)}{\partial\bar{\sigma}}\right|_{\bar{\sigma}=\sigma_{v}}=0.$
(28)

The sigma mean field value $\bar{\sigma}=\sigma_{v}$ at the global minimum of
the grand potential $\Omega(\bar{\sigma},\tilde{\sigma},T)$ is an order
parameter of the phase transition in the finite volume inhomogeneous system at
finite temperatures.

##  IV Numerical results and discussions

Now we will present the numerical results of the systems with the MFA and the
finite-volume effect by MRE, and the system with inhomogeneous background in a
finite-volume sphere respectively. In order to illustrate the phase
transition, the grand potential and the order parameter will be focused on. We
will give the variation curves of grand potential with respect to
$\bar{\sigma}$ at different temperatures in different conditions, and the
corresponding curves of its global minimum or namely the order parameter
$\sigma_{v}$ with respect to $T$ will be also illustrated. The model
parameters are set as
$a=17.7\rm{fm}^{-2},b=-1457.4\rm{fm}^{-1},c=20000,g=12.16$ which are widely
used in the literatureGao:1992prd ; Wilets:1985prd . Thus the quark mass
$m=g\sigma_{v0}$ and the bag constant $B$ are determined according to (25) and
(26). The quark flavors and colors are set to be $N_{f}=2$ and $N_{c}=3$.

As the deconfinement phase transition in FL model at MFA has been extensively
studied in the literatures, here we just discuss the modification to the phase
transition of MFA due to the effect of MRE for the given system size $R$. From
the result of the grand potential (15) with the density of states
$\rho_{\rm{MRE}}$ as shown in (11), the variation curves of the grand
potentials $\Omega_{\rm{MRE}}(\bar{\sigma})$ as functions of $\bar{\sigma}$ at
different temperatures at $R=3\rm{fm}$ are shown in FIG.1LABEL:sub@omg_M+M. In
order to compare the potential curves of different temperatures the curves
have been shifted to have the same fixed zero value of the grand potential at
$\bar{\sigma}=0$. One can see there are two minima corresponding to the two
vacua. At relatively low temperatures the physical vacuum is at the minimum
with nonzero value of $\bar{\sigma}$. With the temperature increasing the
physical vacuum is lifted up. At a certain temperature $T_{c}=136.5MeV$ the
two minima degenerate. The temperature $T_{c}$ is the critical temperature of
the deconfinement phase transition at the system size of $R=3\rm{fm}$. When
the temperature is further increased the physical vacuum is changed to be at
the minimum with zero value of $\bar{\sigma}$ which becomes the global minimum
of the grand potential. As the sigma field value $\sigma_{v}$ at the physical
vacuum is the order parameter, in FIG.1LABEL:sub@sgm_M+M we show how the order
parameter varies with the temperature. The solid line is corresponding to the
case of the $R=3\rm{fm}$. One can see that the order parameter has a sudden
jump at the critical temperature $T_{c}$ which clearly represents a first
order phase transition. If the system size $R$ is increased as $R=4\rm{fm}$
and $R=5\rm{fm}$, one can see that the critical temperature $T_{c}$ becomes
lower. If it approaches the limiting case as $R\to\infty$, the MFA result
which has a critical temperature $T_{c}=120.1MeV$ will be recovered. So one
can see that the description of the first order phase transition at MFA is not
much changed by the finite volume effect via MRE except for a higher critical
temperature at a smaller system size.

(a)

(b)

Figure 1: (a)The variation curves of the grand potentials as functions of
$\bar{\sigma}$ at $R=3\rm{fm}$ for different temperatures. (b)The variation
curve of the order parameter $\sigma_{v}$ as a function of temperature $T$ for
$R=3\rm{fm},4\rm{fm},5\rm{fm}$ and MFA respectively.

Then we will focus on numerical results of the cases with the inhomogeneous
background field $\tilde{\sigma}(r)$ which spatial configuration is determined
by the parameters $\alpha$, $\beta$ and $r_{0}$ as shown in (19). As we know
at a semiclassical level there is a soliton solution in the FL model. The
sigma field in the soliton solution has also a spatial configuration of the
Woods-Saxon type in which the values of the parameters of $\alpha$, $\beta$
and $r_{0}$ can be fitted approximately as $\alpha\approx 0.26\rm{fm}^{-1}$,
$\beta\approx 8\rm{fm}^{-1}$ and $r_{0}\approx 1\rm{fm}$. In this study we
take the configuration of the sigma field in the soliton solution as an
initial configuration. The parameters $\alpha$, $\beta$ and $r_{0}$ can be
varied. Our purpose is to study how the phase transition is influenced by
varying the configuration of the inhomogeneous background field of the sigma.
The parameter $\alpha$ describes the depth of the potential well of sigma
field. The parameters $\beta$ and $r_{0}$ determine the spatial varying rate
of the potential wall and the width of the potential well of the
$\tilde{\sigma}(r)$ respectively. In our study the $\tilde{\sigma}(r)$ is not
the meaning of the hadron but to describe the spatially inhomogeneous fireball
in the heavy ion collision. In the following study the parameter $\alpha$ will
be fixed at $\alpha\approx 0.26\rm{fm}^{-1}$ which corresponding to a fixed
depth of the potential well in the $\tilde{\sigma}(r)$, while $\beta$ and
$r_{0}$ will be varied by hand. The reason for this treating is that we assume
in the heavy ion collisions the difference between the vacuum values of the
sigma field inside the fireball and outside the fireball will change
relatively slowly at the phase transition while the size and shape of the
fireball will change more remarkably.

In the calculation of the grand potential one may notice that there is a
summation over angular momentum number $l$ in the evaluation of the
inhomogeneous part of the density of the states (20). In principle the
summation index $l$ goes to infinity. However the contribution to the density
of states from the scattering of the partial wave will decrease with the
angular momentum number $l$ increasing. As a approximation in the real
numerical calculation we only do the summation over $l=0,1$ which means only
the s and p partial waves have been included in the summation.

Now we are in a position to study how the phase transition in a finite volume
is influenced by varying the spatial configuration of the background field
$\tilde{\sigma}(r)$ through changing the parameters $\beta$ and $r_{0}$. In
the following calculation the system size is taken as $R=3\rm{fm}$. Firstly
let $\beta$ fixed at $\beta=8\rm{fm}^{-1}$ and only $r_{0}$ varied. From FIG.
2LABEL:sub@sgm-r_r0v_1 one can see with $r_{0}$ increasing the width of the
potential well becomes larger. By the equation (17) the scattering phase shift
can be numerically calculated in the potential background with different
$r_{0}$. Then the order parameter $\sigma_{v}$ can be determined at different
temperatures through equation (28) at certain $r_{0}$. In
FIG.2LABEL:sub@sgm_MI_3fm_b=8_set1 the temperature dependence of the
$\sigma_{v}$ for different $r_{0}$ has been plotted. It could be seen that
when $r_{0}$ is very small ($r_{0}\lesssim 0.3\rm{fm}$) which means the
spatial distribution of the inhomogeneous area is very small, the quantum
fluctuation over this tiny inhomogeneous background is also very small. As a
result the phase transition is still a strong first order phase transition
which case is very similar to that of the MRE result. When $r_{0}$ gets larger
it is clear that quantum correction from the inhomogeneous background becomes
notable. In FIG.3LABEL:sub@pt_MI_3fm_(8,1d5) at $r_{0}=1.5\rm{fm}$ it is shown
that besides the global minimum of the physical vacuum a new local minimum
emerges near the zero value of the sigma field in the grand potential at
relatively low temperatures before the phase transition. When the temperature
is greater than the critical temperature this local minimum becomes the global
minimum which is corresponding to the physical vacuum. The corresponding
variation of the order parameter $\sigma_{v}$ with the temperature could be
seen in FIG.2LABEL:sub@sgm_MI_3fm_b=8_set1. It is clear that the first order
phase transition has been obviously weakened. When $r_{0}$ is further
increased ($r_{0}\gtrsim 1.9\rm{fm}$) it is found that the first order phase
transition has been changed into a crossover. One can see that the order
parameter continuously decreasing with temperature increasing in
FIG.2LABEL:sub@sgm_MI_3fm_b=8_set1. From FIG.3LABEL:sub@sgm_MI_3fm_b=8_set2 at
$r_{0}=2.1\rm{fm}$ one can clearly see that there is only one minimum in the
grand potential which is corresponding to unique vacuum in this system and
with temperature increasing it continuously moves towards a small value of the
sigma field which obvious displays a crossover.

(a)

(b)

Figure 2: (a)The spatial configurations of $\tilde{\sigma}(r)$ for the fixed
$\beta=8\rm{fm}^{-1}$ and different $r_{0}$. (b) The $\sigma_{v}$ as functions
of temperature for different inhomogeneous backgrounds represented by
different $r_{0}$ with fixed $\beta$ and $R=3\rm{fm}$.

(a)

(b)

Figure 3: (a)The grand potentials as functions of $\bar{\sigma}$ for
$\beta=8\rm{fm}^{-1}$ and $r_{0}=1.5\rm{fm}$ at $R=3\rm{fm}$ for different
temperatures. (b)The grand potentials as functions of $\bar{\sigma}$ for
$\beta=8\rm{fm}^{-1}$ and $r_{0}=2.1\rm{fm}$ at $R=3\rm{fm}$ for different
temperatures.

In the above analysis one can see that with the spatial size of the
inhomogeneous background increasing the quantum fluctuations become more and
more remarkable which gradually change the first order phase transition into
the crossover. However if one makes a deeper study, it is found that the
critical temperature of the phase transition is not always increasing with
$r_{0}$ increasing. There is a subtle competition and balance between the
contribution from the s wave scattering and p wave scattering. Let us explain
in the following . When $r_{0}\lesssim 0.3\rm{fm}$ the quantum corrections
from both s wave and p wave are small. The critical temperature of the phase
transition almost stays unchanged with $r_{0}$ increasing. When $r_{0}\sim
0.4\rm{fm}$ there emerges a s wave bound state energy level in the energy
spectrum of the Dirac equation (17) which makes a relatively notable
enhancement of the critical temperature of the phase transition as shown in
FIG.2LABEL:sub@sgm_MI_3fm_b=8_set1. At relative small $r_{0}$ the quantum
correction from the s wave is larger than that of the p wave. When s wave
bound orbital appears the first order phase transition is modified remarkably
which means the critical temperature of the phase transition is increased and
the first order phase transition is weakened. However when $r_{0}$ is further
increased from $0.4\rm{fm}$ to $0.9\rm{fm}$ one will find that the critical
temperature of the phase transition is decreasing and the first order phase
transition is getting stronger as shown in FIG.4LABEL:sub@sgm_MI_3fm_b=8_set2.
This abnormal varying trend in the above interval of varying $r_{0}$ is
because the s wave contribution is decreasing while the p wave contribution is
increasing. And the contribution from the p wave counteracts that from the s
wave as long as there is no bound state energy level in the p wave state.
According to the Levinson index theorem the number of the bound stats $n_{l}$
with the angular momentum $l$ in our case is given by
$n_{l}\pi=\bar{\delta}_{l}(k=0)$Ma:1984xv ; Ma:2006zzc . From
FIG.4LABEL:sub@ft1_(8,0d4to1d0) one can see that the p wave phase shift at
zero momentum keeps zero when $r_{0}$ changes from $0.4\rm{fm}$ to
$0.9\rm{fm}$. As a result there is no bound state level in p wave state and
the p wave contribution plays the role of strengthening the first order phase
transition in the above varying interval of $r_{0}$. However the situation
will change until $r_{0}\gtrsim 1.0\rm{fm}$ when there emerges a bound state
energy level in p wave state which can be seen from
FIG.4LABEL:sub@ft1_(8,0d4to1d0) that $\bar{\delta}_{1}(k=0)=\pi$. There is a
sign flipping of the p wave contribution meanwhile its influence becomes
greater than that of the s wave. As a result the critical temperature of the
phase transition is suddenly increased dramatically and the first order phase
transition is remarkably weakened which can be seen from
FIG.4LABEL:sub@sgm_MI_3fm_b=8_set2. After that the p wave contribution
increases with $r_{0}$ increasing and plays the dominant role in modifying the
phase transition. Therefore the first order phase transition is further
weakened with increasing $r_{0}$ and gradually turned into a crossover at
$r_{0}\gtrsim 1.9\rm{fm}$ as shown in FIG.2LABEL:sub@sgm_MI_3fm_b=8_set1.

(a)

(b)

Figure 4: The regulated phase shift for the p wave as functions of momentum
$k$ for fixed $\beta$ and different $r_{0}$. (b)The $\sigma_{v}$ as functions
of temperature for different inhomogeneous backgrounds represented by
different $r_{0}$ with fixed $\beta$ and $R=3\rm{fm}$.

Secondly let us fix $r_{0}=1.5\rm{fm}$ and make $\beta$ varied. In
FIG.5LABEL:sub@sgm-r_bv_1 it is shown that with $\beta$ decreased the
variation of the sigma field over radial distance near the potential wall
become milder. So the parameter $\beta$ determines the steep level of the wall
of the potential well of the sigma field. It is found that with $\beta$
decreased the first order phase transition is weakened and the critical
temperature of the phase transition has been increased as shown in
FIG.5LABEL:sub@sgm_MI_3fm_r0=1d5. When $\beta\gtrsim 4\rm{fm}^{-1}$ there is
one bound state energy level in each of the s wave and p wave state. And the
quantum correction from the p wave scattering is more remarkable than that of
s wave. With $\beta$ decreased the first order phase transition is weakened
and the critical temperature of the phase transition is increased. However it
should be noted that there is a nontrivial variation of the modification to
the phase transition between $\beta=4\rm{fm}^{-1}$ and $\beta=3\rm{fm}^{-1}$.
At $\beta\sim 3\rm{fm}^{-1}$ there emerges a second bound state level in the s
wave state which can be observed from the s wave phase shift. From our
numerical calculation it is found that the s wave phase shift value at zero
momentum jumps from $\bar{\delta}_{0}(k=0)=\pi$ to
$\bar{\delta}_{0}(k=0)=2\pi$. The correction from the s wave scattering is
once again greater than that of the p wave and becomes the dominant part in
the contribution to the density of the states. As a result the first order
phase transition is turned into a crossover. After that the pseudo critical
temperature of the crossover will slowly increases with $\beta$ decreasing as
long as there is no additional bound state level appearing in the s wave or p
wave state.

(a)

(b)

Figure 5: (a)The spatial configurations of $\tilde{\sigma}(r)$ for fixed
$r_{0}=1.5\rm{fm}$ and different $\beta$. (b) The $\sigma_{v}$ as functions of
temperature for different inhomogeneous backgrounds represented by different
$\beta$ with fixed $r_{0}$ and $R=3\rm{fm}$.

Here we want to make some comments on the bound state energy level and how it
influences the phase transition. In the inhomogeneous background field of
sigma the bound state may be formed in the thermal quark system. In principle
these discrete bound state energies may directly contribute additional
isolated terms in the grand potential. However as these energy terms are
constant energies and averaged by the space volume, they just give a small
constant background energy which only shift the whole grand potential upward
or downward. Therefore these isolated terms are neglected in our calculation.
However the discrete bound state can also influence the grand potential by the
density of states. Through the scattering mode of the continuum over the bound
state does the bound state indirectly give the modification to the phase
transition. As the phase shift of the different scattering partial wave will
be changed dramatically at the time when the bound state level emerges, it
will give remarkable modification to the density of states. As a result the
grand potential has been substantially modified accordingly. Our study shows
that the possible bound state in the special partial wave state in the thermal
quark system will largely enhance the quantum scattering effect of the
corresponding partial wave which results in weakening the first order phase
transition or even turning the first order phase transition into a crossover.

Finally we will give some brief discussion about the volume dependence of the
quantum scattering effect in the inhomogeneous background. From the expression
of the density of states (20) one can see that the size $R$ of the system
appears in the denominator. In the surface and curvature term it is order of
$1/R$ and $1/R^{2}$ respectively, and in the inhomogeneous term it is order of
$1/R^{3}$. Therefore by a qualitative analysis it can be deduced that when the
system size is large enough the quantum fluctuations of the inhomogeneous
background can be neglected and the MRE result will be recovered, and if the
system size goes to infinity the MFA result will be recovered. In our case of
study if the system size $R$ is greater than about $10\rm{fm}$, with $r_{0}$
increasing or $\beta$ decreasing the first order phase transition may be
weakened in some degree while the crossover will not appear. Therefore the MRE
result is nearly recovered. When the system size $R$ is greater than about
$20\rm{fm}$ the quantum scattering effect over the inhomogeneous background
could be almost neglected and the MRE effect also becomes very small. Hence
the result is very close to the case of the MFA.

##  V Summary and outlook

In summary we have studied the grand potential and the phase transition in the
thermal quark system with a finite volume under the spatially inhomogeneous
background field. By varying the spatial configuration of the inhomogeneous
background the scattering phase shift for the different partial wave has been
numerically calculated in the corresponding inhomogeneous background at a
certain small spatial volume of the system. It is found that when the bound
state is formed in the thermal system the density of states of momentum will
be substantially modified by the scattering phase shift of the corresponding
partial wave. As a result the first order phase transition in the system with
the given volume has been weakened remarkably and even turned into the
crossover. However with the system volume increasing the quantum corrections
in the inhomogeneous background will decrease. At sufficient large volume of
the system the quantum corrections from the inhomogeneous background can be
neglected.

It should be reminded that the result here may be model dependent. The FL
model is only an effective model and lacks chiral symmetry. The chiral phase
transition could not be discussed in this model. So the interesting problems
are how the scattering phase shift in the chiral models can be calculated in a
spherical system with an inhomogeneous background and how the chiral phase
transition will be modified by the inhomogeneous background with the possible
bound states emerged. In the future study we will use chiral models like the
linear sigma model or the NJL model to address these issues.

## Appendix A

For the spherical symmetrical system one can assume the form of the quark
field as

$\psi(\vec{r})=\frac{1}{r}\left(\begin{array}[]{c}F(r)\\\
i\vec{\sigma}\cdot\vec{e}_{r}G(r)\end{array}\right)y_{jm}^{l},$ (29)

where $y_{jm}^{l}$ is the two-component Pauli spinor harmonic and
$\vec{e}_{r}$ is the spatial unit vector. $F(r)$ and $G(r)$ are the upper and
lower components of the Dirac spinor. From the Dirac equation (17) one can
further derive the first order coupled differential equations of $F(r)$ and
$G(r)$ as

$\left(\frac{d}{dr}-\frac{l+1}{r}\right)F(r)+(g\bar{\sigma}(r)+E)G(r)=0,$ (30)
$\left(\frac{d}{dr}+\frac{l+1}{r}\right)G(r)+(g\bar{\sigma}(r)-E)F(r)=0,$ (31)

Then the coupled equations (30) and (31) can be rewritten into the following
decoupled second order differential equations as

$F^{\prime\prime}-\frac{g\bar{\sigma}(r)^{\prime}}{E+g\bar{\sigma}(r)}F^{\prime}+\left[\frac{l+1}{r}\frac{g\bar{\sigma}(r)^{\prime}}{E+g\bar{\sigma}(r)}-\frac{l(l+1)}{r^{2}}+(E^{2}-g^{2}\bar{\sigma}(r)^{2})\right]F=0,$
(32)
$G^{\prime\prime}+\frac{g\bar{\sigma}(r)^{\prime}}{E-g\bar{\sigma}(r)}G^{\prime}+\left[\frac{l+1}{r}\frac{g\bar{\sigma}(r)^{\prime}}{E-g\bar{\sigma}(r)}-\frac{(l+1)(l+2)}{r^{2}}+(E^{2}-g^{2}\bar{\sigma}(r)^{2})\right]G=0,$
(33)

where the prime denotes the differentiation with respect to $r$. For the
evaluation of the phase shift both equations will give the identical result.
In the following discussion Eq.(32) will be used for the calculation of the
phase shift. When $r>>R$ the asymptotic form of Eq.(32) is

$F^{\prime\prime}-\left[\frac{l(l+1)}{r^{2}}-k^{2}\right]F=0,$ (34)

where $k^{2}=E^{2}-g^{2}\sigma^{2}_{v0}$. The solutions will be spherical
Hankel functions. Meanwhile it should satisfy that $F(r)\to 0$ with $r\to 0$.
There are two linearly independent solutions

$F^{(1)}_{l}(r)=e^{i\beta_{l}(k,r)}rh^{(1)}_{l}(kr),$ (35)
$F^{(2)}_{l}(r)=e^{-i\beta^{*}_{l}(k,r)}rh^{(2)}_{l}(kr),$ (36)

where $h^{(1)}_{l}(kr)$ and $h^{(2)}_{l}(kr)$ are the Hankel functions of the
first and second kinds and $h^{(2)}_{l}(kr)=h^{(1)*}_{l}(kr)$. The function
$\beta_{l}(k,r)$ should satisfy $\beta_{l}(k,r)\to 0$ as $r\to\infty$. Then
the scattering solution is

$F_{l}(r)=F^{(2)}_{l}(r)+e^{i2\delta_{l}(k)}F^{(1)}_{l}(r),$ (37)

and obeys $F_{l}(0)=0$, which leads to the result of the scattering phase
shift

$\delta_{l}(k)=-\textmd{Re}\beta_{l}(k,0),$ (38)

where Re means the real part. It is obvious that the phase shift could be
evaluated from $\beta_{l}$. By substituting $F^{(1)}_{l}$ into Eq.(32) one
could obtain the equation of $\beta_{l}$

$i\beta^{\prime\prime}_{l}rh_{l}+2i\beta^{\prime}_{l}(h_{l}+rh^{\prime}_{l})-\beta^{\prime
2}_{l}rh_{l}-\frac{g\bar{\sigma}(r)^{\prime}}{E+g\bar{\sigma}(r)}(i\beta^{\prime}_{l}rh_{l}+h_{l}+rh^{\prime}_{l})+\left[\frac{l+1}{r}\frac{g\bar{\sigma}(r)^{\prime}}{E+g\bar{\sigma}(r)}-g^{2}(\bar{\sigma}(r)^{2}-\sigma^{2}_{v0})\right]rh_{l}=0.$
(39)

In the fixed background soliton field of $\bar{\sigma}(r)$ this equation could
be numerically solved to obtain the phase shift $\delta_{l}(k)$.

To regularize the phase shift one has to make Born expansion which means
$\beta_{l}$ should be expanded in powers of $g$ as

$\beta_{l}=g\beta_{l1}+g^{2}\beta_{l2}+\cdots.$ (40)

Substituting the expansion (40) into Eq.(39) and neglecting the higher order
terms $O(g^{3})$ one can obtain a set of coupled differential equations about
$\beta_{l1}$ and $\beta_{l2}$ as

$i\beta^{\prime\prime}_{l1}rh_{l}+(2i\beta^{\prime}_{l1}-\frac{\bar{\sigma}(r)^{\prime}}{E})(h_{l}+rh^{\prime}_{l})+(l+1)\frac{\bar{\sigma}(r)^{\prime}}{E}h_{l}=0,$
(41) $i\beta^{\prime\prime}_{l2}rh_{l}-\beta^{\prime
2}_{l1}rh_{l}-\frac{i\bar{\sigma}(r)^{\prime}}{E}\beta^{\prime}_{l1}rh_{l}+(2i\beta^{\prime}_{l2}+\frac{\bar{\sigma}(r)^{\prime}\bar{\sigma}(r)}{E^{2}})(h_{l}+rh^{\prime}_{l})-\left[\frac{(l+1)}{r}\frac{\bar{\sigma}(r)^{\prime}\bar{\sigma}(r)}{E^{2}}+(\bar{\sigma}(r)^{2}-\sigma^{2}_{v0})\right]rh_{l}=0.$
(42)

These equations could be numerically solved to obtain the first and second
Born approximations of the phase shifts namely $\delta^{(1)}_{l}$ and
$\delta^{(2)}_{l}$ as

$\delta^{(1)}_{l}=-g\textmd{Re}\beta_{l1}(k,r=0),\ \ \
\delta^{(2)}_{l}=-g^{2}\textmd{Re}\beta_{l2}(k,r=0).$ (43)

Finally the regularized phase shift can be defined as

$\bar{\delta}_{l}(k)\equiv\delta_{l}(k)-\delta_{l}^{(1)}(k)-\delta_{l}^{(2)}(k),$
(44)

which can be numerically calculated.

## References

  * (1) X. Luo, S. Shi, N. Xu and Y. Zhang, Particles 3, no.2, 278-307 (2020).
  * (2) M. A. Stephanov, K. Rajagopal and E. V. Shuryak, Phys. Rev. Lett. 81, 4816-4819 (1998).
  * (3) E. Shuryak, Lect. Notes Phys. 977, 1-520 (2021).
  * (4) C. S. Fischer, Prog. Part. Nucl. Phys. 105, 1-60 (2019).
  * (5) H. T. Ding, F. Karsch and S. Mukherjee, Int. J. Mod. Phys. E 24, no.10, 1530007 (2015).
  * (6) H. T. Ding, Nucl. Phys. A 1005, 121940 (2021).
  * (7) G. Graef, M. Bleicher and Q. Li, Phys. Rev. C 85, 044901 (2012).
  * (8) B. Klein, Phys. Rept. 707-708, 1-51 (2017).
  * (9) L. Fister and J. M. Pawlowski, Phys. Rev. D 92, no.7, 076009 (2015).
  * (10) L. F. Palhares, E. S. Fraga and T. Kodama, J. Phys. G 38, 085101 (2011).
  * (11) M. Buballa and S. Carignano, Prog. Part. Nucl. Phys. 81, 39-96 (2015).
  * (12) S. Carignano, M. Buballa and B. J. Schaefer, Phys. Rev. D 90, no.1, 014033 (2014).
  * (13) D. Nickel, Phys. Rev. D 80, 074025 (2009).
  * (14) P. Lakaschus, M. Buballa and D. H. Rischke, Phys. Rev. D 103, no.3, 034030 (2021).
  * (15) K. Fukushima, T. Kojo and W. Weise, Phys. Rev. D 102, no.9, 096017 (2020).
  * (16) K. Fukushima and C. Sasaki, Prog. Part. Nucl. Phys. 72, 99-154 (2013).
  * (17) T. Kojo, Y. Hidaka, K. Fukushima, L. D. McLerran and R. D. Pisarski, Nucl. Phys. A 875, 94-138 (2012).
  * (18) Y. Z. Xu, C. Shi, X. T. He and H. S. Zong, Phys. Rev. D 102, 114011 (2020).
  * (19) R. L. Liu, M. Y. Lai, C. Shi and H. S. Zong, Phys. Rev. D 102, no.1, 014014 (2020).
  * (20) Y. P. Zhao, P. L. Yin, Z. H. Yu and H. S. Zong, Nucl. Phys. B 952, 114919 (2020).
  * (21) O. Kiriyama and A. Hosaka, Phys. Rev. D 67, 085010 (2003)
  * (22) R. Balian and C. Bloch, Ann. Phys. (N.Y.) 60, 401 (1970).
  * (23) J. Madsen, Phys. Rev. Lett. 70, 391 (1993); Phys. Rev. D 47, 5156 (1993); 50, 3328 (1994).
  * (24) E.P. Gilson and R.L. Jaffe, Phys. Rev. Lett. 71, 332 (1993).
  * (25) J. Schwinger, Phys. Rev. 94, 1362(1954).
  * (26) E. Farhi, N. Graham, P.E. Haagensen and R.L. Jaffe, Phys.Lett. B 427, 334-342(1998).
  * (27) E. Farhi, N. Graham, R.L. Jaffe and H. Weigel, Nucl.Phys. B 630, 241-268(2002).
  * (28) H. Mao, M. Yao and W.Q. Zhao, Phys. Rev. C 77, 065205 (2008).
  * (29) S. Gao, E.K. Wang and J.R. Li, Phys. Rev. D 46, 3211(1992).
  * (30) V. Paschalidis, K. Yagi, D. Alvarez-Castillo, D. B. Blaschke and A. Sedrakian, Phys. Rev. D 97, no.8, 084038 (2018).
  * (31) M. Ferreira, R. Câmara Pereira and C. Providência, Phys. Rev. D 101, no.12, 123030 (2020).
  * (32) M. Trodden, Rev. Mod. Phys. 71, 1463-1500 (1999).
  * (33) S. Zhou, S. Shu and H. Mao, Chin. Phys. C 45, no.4, 043104 (2021).
  * (34) S. Shu and J. R. Li, Nucl. Phys. A 901, 1-13 (2013).
  * (35) V. Skokov, B. Friman, E. Nakano, K. Redlich and B. J. Schaefer, Phys. Rev. D 82, 034029 (2010).
  * (36) T. K. Herbst, J. M. Pawlowski and B. J. Schaefer, Phys. Rev. D 88, no.1, 014007 (2013).
  * (37) L. Brandes, N. Kaiser and W. Weise, Eur. Phys. J. A 57, no.7, 243 (2021).
  * (38) R. Friedberg, T.D. Lee, Phys. Rev. D 16, 1096 (1977).
  * (39) R. Friedberg, T.D. Lee, Phys. Rev. D 15, 1694 (1977).
  * (40) S. Shu, Sci. China Phys. Mech. Astron. 60, no.4, 041021 (2017).
  * (41) S. Shu, X. Li and J. R. Li, Nucl. Phys. A 1014, 122256 (2021).
  * (42) R. Goldflam and L. Wilets, Phys. Rev. D 25, 1951(1982); M. Bickeboller, M.C. Birse, H. Marschall and L. Wilets, Phys.Rev.D 31, 2892(1985).
  * (43) Z. Q. Ma and G. J. Ni, Phys. Rev. D 31, 1482 (1985).
  * (44) Z. Q. Ma, J. Phys. A 39, R625-R659 (2006).


# Linearly stable and unstable complex soliton solutions with real energies in
the Bullough-Dodd model

Francisco Correa∘, Andreas Fring∙ and Takanobu Taira∙
$\bullet$ Department of Mathematics, City, University of London,
Northampton Square, London EC1V 0HB, UK
$\circ$ Instituto de Ciencias Físicas y Matemáticas, Universidad Austral de
Chile,
Casilla 567, Valdivia, Chile
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract:

We investigate different types of complex soliton solutions with regard to
their stability against linear pertubations. In the Bullough-Dodd scalar field
theory we find linearly stable complex ${\cal{PT}}$-symmetric solutions and
linearly unstable solutions for which the ${\cal{PT}}$-symmetry is broken.
Both types of solutions have real energies. The auxiliary Sturm-Liouville
eigenvalue equation in the stability analysis for the ${\cal{PT}}$-symmetric
solutions can be solved exactly by supersymmetrically mapping it to an
isospectral partner system involving a shifted and scaled inverse
$\cosh$-squared potential. We identify exactly one shape mode in form of a
bound state solution and scattering states which when used as linear
perturbations leave the solutions stable. The auxiliary problem for the
solutions with broken ${\cal{PT}}$-symmetry involves a complex shifted and
scaled inverse $\sin$-squared potential. The corresponding bound and
scattering state solutions have complex eigenvalues, such that when used as
linear perturbations for the corresponding soliton solutions lead to their
decay or blow up as time evolves.

††conference: Stability of complex soliton solutions in the Bullough-Dodd
model

##  1 Introduction

The success of ${\cal{PT}}$-symmetric quantum mechanics [1] and the
applications of the concepts in optical analogues [2] have inspired research
in many different areas of physics. One of them is the study of complex
${\cal{PT}}$-symmetric soliton solutions in nonlinear integrable systems in
1+1 dimensions [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13] or models that admit
complex Skyrmion type solutions in 2+1 dimensions [14]. Despite being complex,
these sort of solutions possess nonetheless real energies, a property that can
be attributed to general arguments that make use of the underlying antilinear
symmetry, generally referred to as ${\cal{PT}}$-symmetry [3, 7, 13, 14].

The other aspect, that makes such type of solutions potentially interesting
physical objects, is their stability. So far this vital property has been
largely ignored in their studies and here the main purpose is to investigate
the stability properties of particular complex one-soliton solutions for a
specific integrable scalar field theory, the Bullough-Dodd model [15, 16].
Besides shedding light on the important question of whether a complex soliton
solution will collapse, i.e. grow or decay, when slightly perturbed, it is
also well-known that as a by-product of a stability analysis one may obtain
the so-called shape mode as a bound state solution for the auxiliary Sturm-
Liouville eigenvalue equation. These solutions are instrumental in the
understanding of energy transfer effects in the dynamics of multi-soliton
solutions that can be unravelled using a collective coordinate or moduli space
approach [17, 18]. For the model studied here we will investigate these moduli
space properties in upcoming work [19].

Our manuscript is organised as follows: In section 2 we assemble some
generalities about the stability analysis to be carried out. In addition, we
briefly recall the general set up for a supersymmetric formulation of quantum
mechanics. While this is not essential for a stability analysis, it greatly
facilitates the construction of exact solutions to the auxiliary Sturm-
Liouville eigenvalue equations. In section 3 we construct complex one-soliton
solutions for the Bullough-Dodd model and analyse in detail their stability.

##  2 Stability analysis, zero modes, shape modes and supersymmetry

The systems considered here are complex scalar field theories with a
Lagrangian density of the general form

$\mathcal{L=}\frac{1}{2}\partial_{\mu}\varphi\partial^{\mu}\varphi-V(\varphi),$
(1)

with $\varphi(x,t)\in\mathbb{C}$ and potential $V(\varphi)\in\mathbb{C}$. We
denote the Lagrangian, i.e. the volume integral over this density, as
$L=\int\nolimits_{-\infty}^{\infty}$ $\mathcal{L}dx$. Our space-time metric is
taken to be Lorentzian of the form $\mathop{\mathrm{d}iag}(1,-1)$, so that the
Euler-Lagrange equation resulting from (1) reads

$\ddot{\varphi}-\varphi^{\prime\prime}+\frac{\partial
V(\varphi)}{\partial\varphi}=0.$ (2)

We use standard conventions and denote partial derivatives with respect to
time $t$ and space $x$ by overdots and dashes, respectively. The static
solutions to these equations, i.e. the solutions to (2) with
$\dot{\varphi}=0$, are fields $\phi(x,\mathbf{u})$ depending on the space
coordinate $x$ and a set of constants that are identified as moduli or
collective coordinates $\mathbf{u}=(u^{1},\ldots,u^{n})$. As is well-known,
one may construct the time-dependent solutions $\varphi$ from the time-
independent solutions $\phi$ simply by a Lorentz boost
$\phi[(x-vt)/\sqrt{1-v^{2}}]=\varphi(x,t)$ with $v$ denoting the velocity. The
energies of particular solutions are computed as

$E[\varphi]=\int\nolimits_{-\infty}^{\infty}dx\varepsilon(\varphi),\qquad\varepsilon(\varphi)=\left(\frac{1}{2}\dot{\varphi}^{2}+\frac{1}{2}(\varphi^{\prime})^{2}+V(\varphi)\right),$
(3)

where $\varepsilon(\varphi)$ denotes the energy density. When the Hamiltonian
associated to (1) admits a modified ${\cal{CPT}}$-symmetry and $\varphi$ is
mapped by this symmetry to itself or possibly to a second degenerate solution,
the energy $E[\varphi]$ is guaranteed to be real, irrespective of whether
$\varphi\in\mathbb{R}$, $\varphi\in\mathbb{C}$ or
$V^{\dagger}(\varphi)=V(\varphi)$ [14]. This reasoning may also be extended to
ensure the reality of other conserved quantities.

With regard to the main purpose of this paper, we now briefly recall the main
argument of a standard stability analysis beginning with the linearisation of
the Euler-Lagrange equation (2) by replacing

$\varphi\rightarrow\varphi_{s}+\varepsilon\delta\varphi\qquad\varepsilon\ll
1,$ (4)

where $\varphi_{s}$ solves (2) and $\delta\varphi=:\chi$ is a small
perturbation. This converts the Euler-Lagrange equation into

$\ddot{\varphi_{s}}-\varphi_{s}^{\prime\prime}+\left.\frac{\partial
V(\varphi)}{\partial\varphi}\right|_{\varphi_{s}}+\varepsilon\left(\ddot{\chi}-\chi^{\prime\prime}+\left.\chi\frac{\partial^{2}V(\varphi)}{\partial\varphi^{2}}\right|_{\varphi_{s}}\right)+\mathcal{O}(\varepsilon^{2})=0.$
(5)

Evidently the first three terms in (5) vanish by construction and when
assuming the perturbation to be separable of the form $\chi(x,t)=e^{i\lambda
t}\Phi(x)$, the term of first order in $\varepsilon$ reduces to a Sturm-
Liouville eigenvalue problem of the same form as the time-independent
Schrödinger equation

$-\Phi_{xx}+V_{1}\Phi=\lambda^{2}\Phi,~{}~{}~{}~{}\ \ \ \ \ \text{with \
}V_{1}(x):=\left.\frac{\partial^{2}V(\varphi)}{\partial\varphi^{2}}\right|_{\varphi_{s}}.$
(6)

For the eigenfunction solutions to this equation with $\lambda\in\mathbb{R}$
the linear perturbation will simply introduce an oscillation in time around
the solution $\varphi_{s}$, whereas when $\lambda\in\mathbb{C}$ the solutions
will grow or decay with time and are therefore unstable.

Let us next see how the energy of the perturbed solution behaves. Expanding up
to second order in $\varepsilon$ and integrating two terms by parts, we easily
derive

$\displaystyle E[\varphi_{s}+\chi]$ $\displaystyle=$ $\displaystyle
E[\varphi_{s}]+\\!\\!\\!\int\limits_{-\infty}^{\infty}dx\left[\left(\left.\frac{\partial
V(\varphi)}{\partial\varphi}\right|_{\varphi_{s}}-\varphi_{s}^{\prime\prime}\right)\chi+\frac{1}{2}\chi\left(\frac{\dot{\chi}^{2}}{\chi}-\chi^{\prime\prime}+\left.\chi\frac{\partial^{2}V(\varphi)}{\partial\varphi^{2}}\right|_{\varphi_{s}}\right)\right]\,\,\,\,\,\,\,\,$
$\displaystyle+\left.\chi(\chi^{\prime}+\varphi_{s}^{\prime})\right|_{-\infty}^{\infty}+\mathcal{O}(\varepsilon^{3}).$

We notice that first bracket of the integrand is simply the Euler-Lagrange
equation for the static solution. The second bracket of the integrand vanishes
when $\chi$ is separable as stated above and $\Phi(x)$ satisfies the Sturm-
Liouville equation (6). The surface term in the second line of (2) is zero
when $\Phi(x)$ vanishes asymptotically
$\lim_{x\rightarrow\pm\infty}\Phi(x)=0$, so that the energy of the perturbed
solution becomes identical to the one of the static solution $\varphi_{s}$.

The _zero mode_ , or translation mode, defined as the solution $\Phi_{0}$ to
(6) with $\lambda=0$, plays a special role for a number of reasons. First of
all, it is well-known [20, 21, 22] that the dimensionality of the
corresponding moduli spaces may be obtained from these modes. Moreover it
allows to directly compute the superpotential and to set up a supersymmetric
scheme that will be important for the concrete system we consider below.

The zero mode is easily computed: Assuming that the static solution $\phi(x)$
belongs to a class of solutions characterised by a continuous parameter, say
$a$, such that $\phi(x)=\psi(x,a)|_{a=0}$, we may think of
$\phi_{s}+\varepsilon\delta\phi$ in (4) as

$\psi(x,a)=\phi(x)+a\left.\frac{\partial\psi(x,a)}{\partial
a}\right|_{a=0}+\mathcal{O}(a^{2}),$ (8)

so that we can identify the zero mode as

$\Phi_{0}(x):=\left.\frac{\partial\psi(x,a)}{\partial a}\right|_{a=0}.$ (9)

We may now repeat the argument for different continuous parameters
$a\rightarrow a_{1},\ldots a_{n}$, so that the number $n$ of independent zero
modes constitutes the dimension of the moduli space.

The excited bound states or scattering states for the eigenvalue equation in
(6) may be constructed in a number of ways. Here it turns out that the
potentials we are interested in possess superpartner potentials that have been
solved in the literature, so that we may employ the scheme and simply map the
known solutions to our problem. Let us therefore briefly recall the general
set up of supersymmetric quantum mechanics by following the ideas of Witten
and others [23, 24, 25, 26]. Assuming that the eigenvalue equation (6) has a
discrete eigenvalue spectrum, $E_{n}^{(1)}=\lambda_{n}^{2}$, we re-write it as

$H_{1}\Phi_{n}^{(1)}=\left(-\frac{d^{2}}{dx^{2}}+V_{1}\right)\Phi_{n}^{(1)}=L_{+}L_{-}\Phi_{n}^{(1)}=E_{n}^{(1)}\Phi_{n}^{(1)},$
(10)

with $\Phi_{0}=:\Phi_{0}^{(1)}$ and $E_{0}=0$. We assumed here that the
Hamiltonian factorises into the product of two first order differential
operators defined as

$L_{\pm}:=\mp\frac{d}{dx}+W(x),~{}~{}\ \ \ \text{with \
}W(x):=-\frac{\Phi_{0}^{\prime}}{\Phi_{0}},$ (11)

where $W(x)$ is referred to as the superpotential related to the potential in
(10) as $V_{1}=W^{2}-W^{\prime}$. Defining next a new Hamiltonian $H_{2}$ with
the product of the operators $L_{\pm}$ in reverse order, we have

$H_{2}\Phi_{m}^{(2)}=\left(-\frac{d^{2}}{dx^{2}}+V_{2}\right)\Phi_{m}^{(2)}=L_{-}L_{+}\Phi_{m}^{(2)}=E_{m}^{(2)}\Phi_{m}^{(2)},$
(12)

with a second potential $V_{2}=W^{2}+W^{\prime}$. Considering now

$\displaystyle H_{2}\left(L_{-}\Phi_{n}^{(1)}\right)$ $\displaystyle=$
$\displaystyle
L_{-}L_{+}L_{-}\Phi_{n}^{(1)}=E_{n}^{(1)}\left(L_{-}\Phi_{n}^{(1)}\right),$
(13) $\displaystyle H_{1}\left(L_{+}\Phi_{m}^{(2)}\right)$ $\displaystyle=$
$\displaystyle
L_{+}L_{-}L_{+}\Phi_{m}^{(2)}=E_{m}^{(2)}\left(L_{+}\Phi_{m}^{(2)}\right),$
(14)

and comparing these equations with (10), (12) we conclude

$\Phi_{n+1}^{(1)}=N_{n+1}^{(1)}L_{+}\Phi_{n}^{(2)},~{}~{}\
~{}\Phi_{n}^{(2)}=N_{n}^{(2)}L_{-}\Phi_{n}^{(1)},~{}~{}\
~{}E_{n}^{(2)}=E_{n+1}^{(1)},~{}~{}\
~{}E_{0}^{(1)}=0,~{}~{}~{}n\in\mathbb{N}_{0}.$ (15)

Thus apart from $E_{0}^{(1)}$ the two Hamiltonians $H_{1}$ and $H_{2}$ are
isospectral, hence the referral as supersymmetric. When
$L_{-}^{\dagger}=L_{+}$ the normalisation constants are simply
$N_{n+1}^{(1)}=1/\sqrt{E_{n}^{(2)}}$, $N_{n}^{(2)}=1/\sqrt{E_{n+1}^{(1)}}$,
but we will not assume this as we allow here for complex superpotentials $W$.
The scheme is directly extended to scattering solutions with continuous
eigenvalues spectra.

The crucial feature for our purposes is here the observation from (15) that
when one of the eigenvalue problems (10) or (12) has been solved, one may
directly deduce the solutions for the other. It is in this sense that we will
apply this scheme as the potentials $V_{1}$ we construct from the zero modes
are not obvious to solve, whereas the solutions for the supersymmetric partner
potentials $V_{2}$ can be found in the literature.

##  3 The Bullough-Dodd model

As an example, we are interested here in the stability properties of the
complex soliton solutions of the Bullough-Dodd model [15, 16], that is a one
scalar field integrable field theory already quite well studied in its
classical [27] and quantum [28] aspects. It is described by the Lagrangian
density of the form

$\mathcal{L}_{\text{BD}}=\frac{1}{2}\partial_{\mu}\varphi\partial^{\mu}\varphi-e^{\varphi}-\frac{1}{2}e^{-2\varphi}+\frac{3}{2}\qquad\text{with}\;\;\varphi\in{\mathbb{C}}.$
(16)

The resulting classical nonlinear equation of motion

$\ddot{\varphi}-\varphi^{\prime\prime}+e^{\varphi}-e^{-2\varphi}=0,$ (17)

can be solved by standard techniques from classical theory of integrable
systems. Here we exploit the fact that the equation may be solved using
Hirota’s direct method [29] when parameterizing the fields as
$\varphi=\ln(\tau_{0}/\tau_{1})$ with $\tau_{i}$ denoting the so-called
$\tau$-functions. To identify various one-soliton solutions we can therefore
make the general Ansatz

$\varphi(x,t)=\ln\left(\frac{\alpha_{0}+\alpha_{1}e^{kx+lt}+\alpha_{2}e^{2kx+2lt}}{\beta_{0}+\beta_{1}e^{kx+lt}+\beta_{2}e^{2kx+2lt}}\right),$
(18)

with unknown constants $\alpha_{i},\beta_{i}$, $i=0,1,2$, and $k,l$.
Substituting this Ansatz into equation (17) and subsequently reading off the
coefficients of $e^{jkx+jlt}$ for $j=0,\ldots,8$ leads to 9 equations, that we
do not report here, which may solved for the unknown constants in (18). In
this manner we find various types of solutions.

###  3.1 Complex one-soliton solutions of type I

The first type of solutions we obtain are one-soliton solutions of the form

$\varphi_{I}^{\pm}(x,t)=\ln\left[\frac{\cosh\left(\beta+\sqrt{k^{2}-3}t+kx\right)\pm
2}{\cosh\left(\beta+\sqrt{k^{2}-3}t+kx\right)\mp
1}\right],\quad\beta\in\mathbb{C}.$ (19)

When $\beta$ is real and $|k|>\sqrt{3}$ the solution $\varphi_{I}^{+}$ is real
but becomes singular for $x_{0}=-[\beta+\sqrt{k^{2}-3}t]/k$, whereas
$\varphi_{I}^{-}$ is only real for
$x<x_{l}=-[\beta+\operatorname{arccosh}(2)+\sqrt{k^{2}-3}t]/k$ and
$x>x_{r}=[-\beta+\operatorname{arccosh}(2)-\sqrt{k^{2}-3}t]/k$ as the argument
of the logarithm becomes negative in the complementary regime. Samples of
these types of solutions in the different regimes are depicted in figure 1. In
the case $|k|<\sqrt{3}$ both solutions $\varphi_{I}^{\pm}$ are always complex,
see figure 3. The type I solutions found here formally coincide with the
solutions constructed in [30].

Despite the fact that these solutions are complex, their energies are real
governed by an antilinear ($\cal{PT}$-symmetry), and for N-soliton solutions
together with the fact that these theories are integrable [7, 13]. The
$\cal{PT}$-symmetry is easily identified as

${\cal{PT}}:x\rightarrow-x,\quad t\rightarrow-t,\quad
i\rightarrow-i,\quad\varphi\rightarrow\varphi.$ (20)

The Bullough-Dodd Lagrangian (16) is trivially invariant under this symmetry
and the solutions $\varphi_{I}^{\pm}(x,t)$ evidently respect it for purely
imaginary constants $\beta$, i.e.
${\cal{PT}}:\varphi_{I}^{\pm}\rightarrow\varphi_{I}^{\pm}$. Thus to take
$\beta\in i\mathbb{R}$ seems to be very suggestive as this choice not only
converts the solutions into ${\cal{PT}}$-symmetric ones, but also regularises
the singularities of the solutions at the cost of a discontinuity resulting
from the branch cut of the logarithm as depicted in figures 1 and 2. The
behaviour of these complex shifted solutions is identical to those for the
$|k|<\sqrt{3}$ and $\beta\in\mathbb{R}$ already depicted in figure 3, when we
identify in these solutions $\beta\rightarrow\sqrt{k^{2}-3}t$ and
$\sqrt{k^{2}-3}t\rightarrow i\beta$ for fixed time.

Figure 1: Scattering solutions (19) of the Bullough-Dodd classical equation of
motion (17) for $|k|>\sqrt{3}$. Panel (a): Travelling real one-cusp solution
$\varphi_{I}^{+}$ for $k=2$ with $\beta=3$ and panel (b) ${\cal{PT}}$
regularized travelling complex solution $\varphi_{I}^{+}$ for $k=2$ with
$\beta=i3/5$. Solid lines correspond to real and dotted lines to imaginary
parts.

Figure 2: Scattering solutions (19) of the Bullough-Dodd classical equation of
motion (17) for $|k|>\sqrt{3}$. Panel (a): Complex solution $\varphi_{I}^{-}$
for $k=2$ with $\beta=3$ and panel (b) ${\cal{PT}}$ regularized complex
solution $\varphi_{I}^{-}$ for $k=2$ with $\beta=i3/5$. Solid lines correspond
to real and dotted lines to imaginary parts.

Figure 3: Complex breather solutions (19) of the Bullough-Dodd classical
equation of motion (17) for $|k|<\sqrt{3}$. Panel (a): $\varphi_{I}^{+}$ for
$k=1<\sqrt{3}$ with $\beta=2$ and panel (b) $\varphi_{I}^{-}$ for
$k=1<\sqrt{3}$ with $\beta=2$. Solid lines correspond to real and dotted lines
to imaginary parts.

Indeed the energies for the solutions (19), evaluated according to (3), are
always real

$E[\varphi_{I}^{\pm}]=-6|k|,$ (21)

and independent of the complex parameter $\beta$.

####  3.1.1 Zero modes for type I complex one-soliton solutions

The Sturm-Liouville eigenvalue problem (6) of the stability analysis for the
Bullough-Dodd model reads

$-\Phi_{xx}+V(x)\Phi=\lambda^{2}\Phi\qquad\text{with}\quad
V(x)=\left(e^{\phi(x)}+2e^{-2\phi(x)}\right).$ (22)

The static solution for the cusp and oscillatory solutions
$\varphi^{\pm}_{I}(x,t)$ acquire the form

$\phi^{\pm}_{I}(x)=\ln\left[\frac{\cosh\left(\beta+\sqrt{3}x\right)\pm
2}{\cosh\left(\beta+\sqrt{3}x\right)\mp 1}\right],$ (23)

which when substituted into (22) yield the potentials

$\displaystyle V_{1}^{+}(x)$ $\displaystyle=$ $\displaystyle
1-\frac{3}{1-\cosh\left(\beta+\sqrt{3}x\right)}+\frac{8\sinh^{4}\left[\frac{1}{2}\left(\beta+\sqrt{3}x\right)\right]}{\left[2+\cosh\left(\beta+\sqrt{3}x\right)\right]^{2}},$
(24) $\displaystyle V_{1}^{-}(x)$ $\displaystyle=$ $\displaystyle
1-\frac{3}{1+\cosh\left(\beta+\sqrt{3}x\right)}+\frac{8\cosh^{4}\left[\frac{1}{2}\left(\beta+\sqrt{3}x\right)\right]}{\left[2-\cosh\left(\beta+\sqrt{3}x\right)\right]^{2}}.$
(25)

We think now of the larger class of functions depending on the continuous
parameter $a$ as $\psi^{\pm}(x,a)=\phi^{\pm}(x+a/\sqrt{3})$. By means of (9)
the zero mode is then easily computed to

$\Phi_{0}^{\pm}(x)=-\frac{3\left\\{\tanh\left[\frac{1}{2}\left(\beta+\sqrt{3}x\right)\right]\right\\}^{\mp
1}}{2\pm\cosh\left(\beta+\sqrt{3}x\right)}.$ (26)

One verifies that the $\Phi_{0}^{\pm}(x)$ indeed satisfy (22) for the
potentials (25) with eigenvalue $\lambda=0$. Notice that when taking
$a\rightarrow ia$, as seems to be natural when $k<\sqrt{3}$, so that we
identify $\psi^{\pm}(x,a)=\phi^{\pm}(x+ia/\sqrt{3})$, the resulting zero mode
is the same as in (26) only multiplied by $i$, which simply corresponds to a
different normalization constant.

####  3.1.2 Bound state solutions for type I complex one-soliton solutions

It is not straightforward or obvious how to solve the eigenvalue equation for
the potentials $V_{1}^{\pm}(x)$ directly. However, we notice now that the
supersymmetric partner Hamiltonians are inverse $\cosh^{2}$-potentials that
have been treated extensively in the literature. Since
$V_{1}^{+}(x,\beta)=V_{1}^{-}(x,\beta+i\pi)$, we just focus on the “+”-case
from now on, as all quantities related to the “-”-system are computed simply
by a shift in $\beta$. Using the relation (11) we compute the superpotential
from the zero modes to

$W^{+}(x)=\frac{\sqrt{3}}{2}\frac{3+2\cosh\left(\beta+\sqrt{3}x\right)+\cosh\left[2\left(\beta+\sqrt{3}x\right)\right]}{\left[\cosh\left(\beta+\sqrt{3}x\right)+2\right]\sinh\left(\beta+\sqrt{3}x\right)},$
(27)

which is sufficient to define the corresponding intertwining operators
$L_{\pm}$ in (11). Using next the defining relation for the second Hamiltonian
$H_{2}$ we compute the associated second potential from (12) to

$V_{2}=3-\frac{3}{2}\text{sech}^{2}\left(\frac{\beta}{2}+\frac{\sqrt{3}x}{2}\right).$
(28)

When shifting the overall energy by $3$, this is an inverse
$\cosh^{2}$-potential discussed in the literature. It possess bound states as
treated in section 23 in problem 5 [31] and scattering states for which the
potential is entirely reflectionless [32], see also [33]. From [31] we read
off the eigenfunction for the energy $E$ with potential
$U=-U_{0}\operatorname{sech}^{2}(\alpha x)$, $\alpha,U_{0}\in\mathbb{R}$ as

$\Psi(x)=\left[1-\tanh^{2}(\alpha
x)\right]^{\epsilon/2}\,_{2}F_{1}\left[\epsilon-s,s+\epsilon+1;\epsilon+1;\frac{1}{1+e^{2x\alpha}}\right],$
(29)

where $\,{}_{2}F_{1}(a,b;c;z)$ denotes the hypergeometric function and

$\epsilon:=\frac{\sqrt{-2Em}}{\alpha\hbar},\qquad
s:=\frac{1}{2}\left(\sqrt{1+\frac{8mU_{0}}{\alpha^{2}\hbar^{2}}}-1\right).$
(30)

The quantization condition for the energy emerges from the requirement
$\lim_{x\rightarrow\pm\infty}\Psi(x)=0$ as $\epsilon-s=-n$ for
$n\in\mathbb{N}_{0}$, which when solved yields

$E_{n}=-\frac{\alpha^{2}\hbar^{2}}{8m}\left(\sqrt{\frac{8mU_{0}}{\alpha^{2}\hbar^{2}}+1}-2n-1\right)^{2}.$
(31)

As argued in [31], the number of bound states is limited by the constraint
$n<s$. When shifting the overall energy by $3$ we obtain $V_{2}=U$ for the
parameter identifications $\hbar=1$, $m=1/2$, $U_{0}=3/2$,
$\alpha=\sqrt{3}/2$. Since for these values $s=1$, there is only one bound
state for $n=0$. Thus we obtain the eigensystem solutions

$\displaystyle\Phi_{0}^{(2)}(x)$ $\displaystyle=$ $\displaystyle
N_{0}\text{sech}\left[\frac{1}{2}\left(\beta+\sqrt{3}x\right)\right],\qquad
E_{0}^{(2)}=\frac{9}{4},$ (32) $\displaystyle\Phi_{1}^{(1)}(x)$
$\displaystyle=$ $\displaystyle
L_{+}\Phi_{0}^{(2)}=N_{0}\frac{3\sqrt{3}\cosh\left[\frac{1}{2}\left(\beta+\sqrt{3}x\right)\right]\coth\left(\beta+\sqrt{3}x\right)}{\cosh\left(\beta+\sqrt{3}x\right)+2},\qquad
E_{1}^{(1)}=\frac{9}{4}.\qquad\,\,$ (33)

Notice that the singularity in $\Phi_{1}^{(1)}(x)$ for $x_{0}=-\beta\sqrt{3}$
when $\beta\in\mathbb{R}$, inherited form the superpotential $W(x)$, is
eliminated in the $\cal{PT}$-symmetric solution when $\beta\in i\mathbb{R}$.
Since we have obtained a well-defined asymptotically vanishing solution,
$\lim_{x\rightarrow\pm\infty}\Phi_{1}^{(1)}(x)=0$ with real eigenvalues, we
conclude that the complex type I soliton solution remains stable when linearly
perturbed by the shape mode $\Phi_{1}^{(1)}$. Let us next consider the case
when $\lambda$ is taken to be a continuous parameter.

####  3.1.3 Scattering state solutions for type I complex one-soliton
solutions

The scattering states for the inverse $\cosh$-squared potential are discussed
in detail in [32]. When adjusting the parameters therein to our equation (22)
for $V_{2}$, we obtain the two scattering solutions

$\displaystyle\Psi_{e}^{(2)}$ $\displaystyle=$
$\displaystyle\cosh^{2}\left(\frac{\hat{x}}{2}\right)\,_{2}F_{1}\left[\gamma_{-},\gamma_{+};\frac{1}{2};-\zeta^{2}(\hat{x})\right]=\cos\left(\tilde{x}\right)-\frac{\sqrt{3}}{2\kappa}\tanh\left(\frac{\hat{x}}{2}\right)\sin\left(\tilde{x}\right),\qquad$
(34) $\displaystyle\Psi_{o}^{(2)}$ $\displaystyle=$
$\displaystyle\zeta(\hat{x})\cosh^{2}\left(\frac{\hat{x}}{2}\right)\,_{2}F_{1}\left[\frac{1}{2}+\gamma_{-},\frac{1}{2}+\gamma_{+};\frac{3}{2};-\zeta^{2}(\hat{x})\right]$
$\displaystyle=$
$\displaystyle\frac{2\sqrt{3}\kappa\sin\left(\tilde{x}\right)+3\tanh\left(\frac{\hat{x}}{2}\right)\cos\left(\tilde{x}\right)}{4\kappa^{2}+3},$

with $\gamma_{\pm}:=1\pm i\kappa/\sqrt{3}$, $\zeta(x):=\sinh\left(x/2\right)$,
$\tilde{x}:=\kappa\left(\frac{\beta}{\sqrt{3}}+x\right)$,
$\hat{x}:=\beta+\sqrt{3}x$ and continuous eigenvalues
$\lambda^{2}=3+\kappa^{2}$ for $\kappa$ being a free parameter. As we expect
from Abel’s identity, the Wronskian for these two solutions is constant,
$W(\Psi_{e}^{(2)},\Psi_{o}^{(2)})=\Psi_{e}^{(2)}d\Psi_{o}^{(2)}/dx-d\Psi_{e}^{(2)}/dx\Psi_{o}^{(2)}=\sqrt{3}/2$,
indicating that they are indeed linearly independent. Thus according to (15),
the two scattering states for the potential $V_{1}$ are obtained as

$\displaystyle\Psi_{e}^{(1)}$ $\displaystyle=$ $\displaystyle
L_{+}\Psi_{e}^{(2)}=\frac{1}{\kappa}\left[f_{-}(\hat{x})\cos\left(\tilde{x}\right)-f_{+}(\hat{x})\sin\left(\tilde{x}\right)\right],$
(36) $\displaystyle\Psi_{o}^{(1)}$ $\displaystyle=$ $\displaystyle
L_{+}\Psi_{o}^{(2)}=\frac{2\sqrt{3}}{3+4\kappa^{2}}\left[f_{+}(\hat{x})\cos\left(\tilde{x}\right)+f_{-}(\hat{x})\sin\left(\tilde{x}\right)\right],$
(37)

where

$f_{+}(x):=\frac{3}{2}-\kappa^{2}-\frac{9}{2\cosh\left(x\right)+4},\qquad
f_{-}(x):=\frac{3\sqrt{3}\kappa\cosh(x)\coth\left(\frac{x}{2}\right)}{2\cosh(x)+4}.$
(38)

Once again the singularity at $x_{0}=-\beta\sqrt{3}$ is eliminated in the
$\cal{PT}$-symmetric solution by taking $\beta\in i\mathbb{R}$. Now the
Wronskian for these two solutions is computed to
$W(\Psi_{e}^{(1)},\Psi_{o}^{(1)})=(3+\kappa^{2})\sqrt{3}/2$. When
$\kappa=\kappa_{0}:=\pm i\sqrt{3}$ we recover the zero mode $\Phi_{0}^{+}$
from $\Psi_{e}^{(1)}$, whereas
$\lim_{x\rightarrow\kappa_{0}}\Psi_{o}^{(1)}(x)=0$.

Since the scattering solutions are well-defined with real continuous
eigenvalues one may be tempted to conclude that the complex type I soliton
solution remains stable when linearly perturbed by either of the two
scattering solutions $\Psi_{e}^{(1)}$ or $\Psi_{o}^{(1)}$. However, as the
limits for the asymptotic values of $x$ are not well defined, the energies of
the perturbed solution (2) do not give definite values so that the
perturbation by these scattering solutions leads to ill-defined objects. A
natural way to overcome this issue would be to restrict the domain of the
theory to finite interval.

An interesting solution arises for $\kappa=0$

$\lim_{\kappa\rightarrow
0}\Psi_{o}^{(1)}=\sqrt{3}\left(1-\frac{3}{\cosh\left(\beta+\sqrt{3}x\right)+2}\right)=:\Psi_{0}^{(1)},$
(39)

which has a well-defined finite limit $\sqrt{3}$ for $x\rightarrow\pm\infty$.
Since the derivative of this solution is asymptotically vanishing as well as
the derivative of the static solution $\phi_{I}$, the second line in (2) is
vanishing up to order $\varepsilon^{2}$ so that the energy of the perturbed
solution is the same as the original solution. Since we have also a positive
eigenvalue for the auxiliary eigenvalue problem we conclude that the solution
remains stable when perturbed by $\Psi_{0}^{(1)}$.

###  3.2 Complex one-soliton solutions of type II

The second type of one-soliton solutions obtained from (18) by fixing the
unknown constants from direct substitution into the Bullough-Dodd equation
(17) are always complex and do not possess a real regime in their range

$\varphi_{II}^{\pm}(x,t)=\ln\left\\{\omega^{\mp
2}\left[1-\frac{6e^{t\sqrt{k^{2}+3\omega^{\pm
1}}+kx\pm\beta}}{\left(1+e^{t\sqrt{k^{2}+3\omega^{\pm
1}}+kx\pm\beta}\right)^{2}}\right]\right\\}.$ (40)

Here $\omega=e^{i\pi/3}$ denotes the third root of unity. A sample solution is
depicted in figure 4, from which we observe that the real part of this
solution is an oscillation between a regular shaped soliton solution and a
double peakon solution.

The energies when directly computed with (3) for the Lagrangian densities
$\mathcal{L}_{\text{BD}}$ are infinite. However, when shifting the vacuum from
$\varphi_{0}=0$ to the complex plane $\varphi_{0}^{\pm}=\mp i2\pi/3$ amounts
to defining the new Lagrangians densities

$\mathcal{L}_{\text{BD}}^{\pm}=\mathcal{L}_{\text{BD}}+\frac{3}{2}\left(\omega^{\mp
2}-1\right).$ (41)

Computing the energies of the solution $\varphi_{II}^{\pm}$ with these
Lagrangians leads to the same real values as those computed for the type I
solution with $\mathcal{L}_{\text{BD}}$

$E^{+}[\varphi_{II}^{+}]=E^{-}[\varphi_{II}^{-}]=-6|k|.$ (42)

The ${\cal{PT}}$-symmetry for the individual Lagrangians is now broken, but
instead we have ${\cal{PT}}$:
$\mathcal{L}_{\text{BD}}^{+}\rightarrow\mathcal{L}_{\text{BD}}^{-}$. Since the
energies are real this suggests that the criteria provided in [14] may still
be slightly enlarged. One may easily convince oneself that the energies of two
${\cal{PT}}$-symmetrically related degenerate solutions are real even when
they originate from two different ${\cal{PT}}$-symmetrically related
Hamiltonians. Having obtained real energies suggests that these solutions are
potentially well-defined physical objects, but let us investigate whether the
solutions are also stable.

Figure 4: Complex one-soliton solutions of type II (40) with $\beta=2$,
$k=1/2$ for the Bullough-Dodd classical equation of motion (17). Solid lines
correspond to real and dotted lines to imaginary parts.

####  3.2.1 Zero modes for complex type II one-soliton solutions

The static solutions corresponding to $\varphi_{II}^{\pm}(x,t)$ are now
obtained as

$\phi_{II}^{\pm}(x)=\ln\left[\omega^{\mp 2}\frac{\cos(\sqrt{3}\omega^{\pm
1/2}x\pm\beta)-2}{\cos(\sqrt{3}\omega^{\pm 1/2}x\pm\beta)+1}\right].$ (43)

Taking $\beta\in\mathbb{R}$ the two solutions are related to each other by the
${\cal{PT}}$-symmetry in (20) as
${\cal{PT}}:\phi_{II}^{+}\rightarrow\phi_{II}^{-}$. Alternatively, we may also
change the sign in $\beta$ in $\phi_{II}^{-}$ and take $\beta\in i\mathbb{R}$
in which case the two solutions are still related by the same symmetry. The
potentials for the eigenvalue equation as introduced in (22) result to

$V_{1}^{\pm}=\omega^{\mp 2}\left\\{1-\frac{3}{\cos\left(\sqrt{3}\omega^{\pm
1/2}x\pm\beta\right)+1}+2\frac{[\cos(\sqrt{3}\omega^{\pm
1/2}x\pm\beta)+1]^{2}}{[\cos(\sqrt{3}\omega^{\pm
1/2}x\pm\beta)-2]^{2}}\right\\}.$ (44)

Since the two potentials are not ${\cal{PT}}$-symmetric by themselves one may
be tempted to conclude already at this stage that there are no real eigenvalue
solutions to (22). However, all we can deduce with certainty at this point is
that the symmetry as defined in (20) is broken and will not ensure the reality
of $\lambda$. We can not guarantee that there are no other types antilinear
symmetries that would explain a real spectrum.

Thus let us show explicitly how the eigensystem may be solved in this case.
The function $\psi$ in (9) is now taken to
$\psi^{\pm}(x,a)=\phi^{\pm}(x+a\omega^{\mp 1/2})$, from which we easily
calculate the zero mode to

$\Phi_{0}^{\pm}(x)=\mp\frac{3^{3/2}\tan\left[\frac{1}{2}(\sqrt{3}\omega^{\pm
1/2}x\pm\beta)\right]}{\cos\left(\sqrt{3}\omega^{\pm 1/2}x\pm\beta\right)-2}.$
(45)

We notice that in these solutions there is no need to take $\beta$ to be
purely imaginary as the singularities in the zero mode are already regularized
because $x$ is multiplied by a complex constant.

####  3.2.2 Bound state solutions for type II complex one-soliton solutions

Next we employ relation (11) to compute the superpotentials from the zero
modes

$W^{\pm}(x)=\pm\frac{\sqrt{3}\omega^{\pm
1/2}}{2}\frac{[3+\cos\left[2\left(\sqrt{3}\omega^{\pm
1/2}x\pm\beta\right)\right]-2\cos\left(\sqrt{3}\omega^{\pm
1/2}x\pm\beta\right)]}{[\cos\left(\sqrt{3}\omega^{\pm
1/2}x\pm\beta\right)-2]\sin\left(\sqrt{3}\omega^{\pm 1/2}x\pm\beta\right)},$
(46)

from which we calculate the corresponding intertwining operators $L_{\pm}$ in
(11). The partner potentials are then obtained by a direct calculation from
(12) to

$V_{2}^{\pm}=\frac{3\omega^{\mp 2}}{1-\sec\left(\sqrt{3}\omega^{\pm
1/2}x\pm\beta\right)}=\frac{3}{2}\frac{\omega^{\pm
1}}{\sin^{2}\left[\frac{1}{2}(\pm\beta+\sqrt{3}\omega^{\pm
1/2}x)\right]}-3\omega^{\pm 1}.$ (47)

Again, this is a well studied potential corresponding to a special version of
the Pöschl-Teller potential [35, 33, 36, 37]. Adapting the shift and scaling
to the solutions that can be found in the literature, we obtain the
eigensystem solutions

$\displaystyle\Phi_{0}^{(2)\pm}$ $\displaystyle=$ $\displaystyle
N_{0}\frac{1}{\sin\left[\frac{1}{2}\left(\pm\beta+\sqrt{3}\omega^{\pm
1/2}x\right)\right]},\qquad E_{0}^{(2)}=-\frac{9}{4}\omega^{\pm 1},$ (48)
$\displaystyle\Phi_{1}^{(1)\pm}$ $\displaystyle=$ $\displaystyle
L_{+}\Phi_{0}^{(2)}=\frac{N_{0}}{2}\frac{3^{3/2}\omega^{\pm
1/2}\cos\left(\pm\beta+\sqrt{3}\omega^{\pm
1/2}x\right)}{\left[2-\cos\left(\pm\beta+\sqrt{3}\omega^{\pm
1/2}x\right)\right]\sin\left[\frac{1}{2}\left(\pm\beta+\sqrt{3}\omega^{\pm
1/2}x\right)\right]},\,\,E_{1}^{(1)}=E_{0}^{(2)}.$

The solutions are asymptotically vanishing,
$\lim_{x\rightarrow\pm\infty}\Phi_{1}^{(1)\pm}(x)=0$ and come in complex
conjugate pairs, which is the signature property of a spontaneously broken
${\cal{PT}}$-symmetry. Instead we have
${\cal{PT}}:\Phi_{1}^{(1)+}(x)\rightarrow\Phi_{1}^{(1)-}(x)$. We conclude that
the complex soliton solutions of type II become unstable when perturbed by the
bound state solutions $\Phi_{1}^{(1)\pm}$.

####  3.2.3 Scattering solutions for type II complex one-soliton solutions

Finally we discuss the scattering solutions for the auxiliary Sturm-Liouville
eigenvalue equation obtained from the type II complex one-soliton solutions.
By similar means as in the previous section we exploit the fact that the
potential becomes a shifted and scaled inverse $\sin$-squared potential. For
the $V_{1}^{\pm}$-potential we then find the two linearly independent
scattering solutions

$\displaystyle\Psi_{e}^{(2)\pm}$ $\displaystyle=$ $\displaystyle
g^{\pm}(x)\sin(\kappa x)+h^{\pm}(x)\cos(\kappa x),$ (49)
$\displaystyle\Psi_{o}^{(2)\pm}$ $\displaystyle=$ $\displaystyle
g^{\pm}(x)\cos(\kappa x)+h^{\pm}(x)\sin(\kappa x),$ (50)

with

$\displaystyle g^{\pm}(x)$ $\displaystyle=$
$\displaystyle\pm\frac{3^{3/2}\omega^{\pm
1/2}\kappa}{2}\frac{\cos\left(\beta+\sqrt{3}\omega^{\pm
1/2}x\right)\tan\left[\frac{1}{2}\left(\beta+\sqrt{3}\omega^{\pm
1/2}x\right)\right]}{\cos\left(\beta+\sqrt{3}\omega^{\pm 1/2}x\right)-2}$ (51)
$\displaystyle h^{\pm}(x)$ $\displaystyle=$
$\displaystyle\frac{3}{2}\omega^{\pm
1}\left[\frac{3}{2-\cos\left(\beta+\sqrt{3}\omega^{\pm
1/2}x\right)}-1\right]-\kappa^{2}$ (52)

for the complex continuous eigenvalues
$(\lambda^{\pm})^{2}=\kappa^{2}-3\omega^{\pm}$. One may be tempted to choose
the free parameter $\kappa$ in such a way that $(\lambda^{\pm})^{2}$ becomes
real. However, when $\kappa$ acquired a complex part the solutions are easily
seen to diverge at $x\rightarrow\pm\infty$. Thus overall we conclude that the
complex type II one-solitons solutions become unstable when linearly perturbed
with these scattering solutions $\Psi_{e,o}^{(2)\pm}$.

##  4 Conclusions

We investigated several types of complex one-soliton solutions in the
Bullough-Dodd model with regard to their stability when linearly perturbed.
The type I solutions were previously known, whereas the type II solutions have
been newly constructed. The overall energy for both types of solutions were
found to be real when shifting the vacuum of the potential appropriately.

We found that the auxiliary Sturm-Liouville eigenvalue equation of the
stability analysis for the type I solutions can be supersymmetrically mapped
to an isospectral partner system involving a shifted and scaled inverse
$\cosh$-squared potential. Remarkably these potentials also emerge in the
stability analysis of the $\phi^{4}$-theory and sine-Gordon model, but simply
with different scalings and overall shifts [34]. We concluded that the type I
solutions remain stable when linearly perturbed with bound state and a
scattering solutions with vanishing continuous parameter.

For the type II solutions the original auxiliary Sturm-Liouville eigenvalue
equation can be supersymmetrically mapped to a complex shifted and scaled
inverse $\sin$-squared potential. The corresponding bound state and scattering
solutions for these equations have complex eigenvalues or are asymptotically
divergent, such when used to perturb the type II solutions they become
unstable or the energies become infinite.

The underlying ${\cal{PT}}$-symmetry serves for several purposes: i) When
introducing complex rather than real shifts the singularities in the solutions
are regularized. ii) As previously pointed out, it governs the reality of the
energy of the solutions, but as we found here the set criteria provided in
[14] may still be enlarged. The energies of two ${\cal{PT}}$-symmetrically
related degenerate solutions are guaranteed to be real even when they solve
the Euler-Lagrange equations of two ${\cal{PT}}$-symmetrically related
Hamiltonians that might even be different. iii) As we have seen here the
${\cal{PT}}$-symmetry also governs the stability of the solutions. However, in
order to obtain stable solutions, that is real eigenvalues in the auxiliary
Sturm-Liouville equation, the perturbation must be ${\cal{PT}}$-symmetric by
themselves.

Here we observed that the ${\cal{PT}}$-symmetry of the soliton solution
$\varphi$ was inherited by the perturbing field $\chi$. The type I solutions
and their perturbing fields were found to be ${\cal{PT}}$-symmetry by
themselves, which guaranteed the reality of the energy and their stability. On
the other hand the type II solutions were mapped into different degenerate
solutions, so that the reality of their energy was still ensured. However, on
the level auxiliary of the Sturm-Liouville equation the ${\cal{PT}}$ is broken
for a specific solution so that the eigenvalues become complex and the
solutions unstable.

Acknowledgments: FC was partially supported by Fondecyt grant 1211356. TT is
supported by EPSRC grant EP/W522351/1.

## References

  * [1] C. M. Bender, P. E. Dorey, C. Dunning, A. Fring, D. W. Hook, H. F. Jones, S. Kuzhel, G. Levai, and R. Tateo,  PT Symmetry: In Quantum and Classical Physics,  (World Scientific, Singapore) (2019).
  * [2] R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter, and D. N. Christodoulides,  Non-Hermitian physics and ${\cal{PT}}$ symmetry,  Nature Phys. 14, 11 (2018).
  * [3] A. Fring,  $\cal{PT}$-Symmetric deformations of the Korteweg-de Vries equation,  J. Phys. A: Math. Theor. 40, 4215–4224 (2007).
  * [4] M. Wadati,  Construction of parity-time symmetric potential through the soliton theory,  J. Phys. Soc. Japan 70, 074005 (2008).
  * [5] P. E. G. Assis and A. Fring,  Compactons versus solitons,  Pramana J. of Phys. 74, 857–865 (2010).
  * [6] A. Cavaglia, A. Fring, and B. Bagchi,  $\cal{PT}$-symmetry breaking in complex nonlinear wave equations and their deformations,  J. Phys. A: Math. Theor. 44, 325201(42) (2011).
  * [7] J. Cen and A. Fring,  Complex solitons with real energies,  J. Phys. A: Math. Theor. 49, 365202 (2016).
  * [8] F. Correa and A. Fring,  Regularized degenerate multi-solitons,  J. of High Energy Phys. 2016(9), 1–16 (2016).
  * [9] A. Khare and A. Saxena,  Novel ${\cal{PT}}$-invariant solutions for a large number of real nonlinear equations,  Phys. Lett. A 380, 856–862 (2016).
  * [10] M. J. Ablowitz and Z. H. Musslimani,  Integrable nonlocal nonlinear Schrödinger equation,  Phys. Rev. Lett. 110, 064105 (2013).
  * [11] J. Cen, F. Correa, and A. Fring,  Integrable nonlocal Hirota equations,  J. of Math. Phys. 60, 081508 (2019).
  * [12] J. Cen,  Nonlinear classical and quantum integrable systems with PT-symmetries,  PhD thesis, City, University of London, City Research Online digital archive (2019).
  * [13] A. Fring and T. Taira,  Complex BPS solitons with real energies from duality,  J. of Phys. A: Math. Theor. 53, 455701 (2020).
  * [14] F. Correa, A. Fring, and T. Taira,  Complex BPS Skyrmions with real energy,  Nucl. Phys. B 971, 115516 (2021).
  * [15] R. Dodd and R. Bullough,  Polynomial conserved densities for the sine-Gordon equations,  Proc. of the Royal Soc. of London. A. Math. and Phys. Sciences 352(1671), 481–503 (1977).
  * [16] A. V. Zhiber and A. B. Shabat,  The Klein–Gordon equation with nontrivial group,  in Doklady Akademii Nauk, 247, 1103–1107, Russian Academy of Sciences, 1979.
  * [17] I. Takyi and H. Weigel,  Collective coordinates in one-dimensional soliton models revisited,  Phys. Rev. D 94, 085008 (2016).
  * [18] N. S. Manton, K. Oleś, T. Romańczukiewicz, and A. Wereszczyński,  Kink moduli spaces: Collective coordinates reconsidered,  Phys. Rev. D 103, 025024 (2021).
  * [19] F. Correa, A. Fring, and T. Taira,  in preparation.
  * [20] R. Jackiw,  Quantum meaning of classical field theory,  Rev. of Mod. Phys. 49, 681 (1977).
  * [21] E. J. Weinberg,  Parameter counting for multimonopole solutions,  Phys. Rev. D 20, 936 (1979).
  * [22] S. K. Kim, K. S. Soh, and J. H. Yee,  Index theory for the nonrelativistic Chern-Simons solitons,  Phys. Rev. D 42, 4139 (1990).
  * [23] E. Witten,  Dynamical breaking of supersymmetry,  Nucl. Phys. B188, 513 (1981).
  * [24] F. Cooper and B. Freedman,  Aspects of supersymmetric quantum mechanics,  Annals Phys. 146, 262 (1983).
  * [25] E. Witten,  Constraints on supersymmetry breaking,  Nucl. Phys. B202, 253 (1982).
  * [26] F. Cooper, A. Khare, and U. Sukhatme,  Supersymmetry and quantum mechanics,  Phys. Rept. 251, 267–385 (1995).
  * [27] V. A. Andreev,  Bäcklund transformations of the Bullough-Dodd-Zhiber-Shabat equation and symmetries of integrable equations,  Theor. and Math. Phys. 79, 448–450 (1989).
  * [28] A. Fring, G. Mussardo, and P. Simonetti,  Form-factors of the elementary field in the Bullough-Dodd model,  Phys. Lett. B 307, 83–90 (1993).
  * [29] R. Hirota,  The direct method in soliton theory, volume 155,  Cambridge University Press, 2004.
  * [30] P. E. G. Assis and L. A. Ferreira,  The Bullough–Dodd model coupled to matter fields,  Nucl. Phys. B 800, 409–449 (2008).
  * [31] L. D. Landau and E. M. Lifschitz,  Quantum mechanics non-relativistic theory, second ed.,  (Pergamon Press, Oxford) (1965).
  * [32] J. Lekner,  Reflectionless eigenstates of the $\operatorname{sech}^{2}$ potential,  Am. J. of Phys. 75, 1151–1157 (2007).
  * [33] S. Flügge,  Practical quantum mechanics,  (Springer Science & Business Media) (2012).
  * [34] R. Jackiw,  Quantum meaning of classical field theory,  Rev. of Mod. Phys. 49, 681 (1977).
  * [35] G. Pöschl and E. Teller,  Bemerkungen zur Quantenmechanik des anharmonischen Oszillators,  Phys. Rev. 85, 166-179 (1952).
  * [36] F. Correa and M.S. Plyushchay,  Spectral singularities in ${\cal{PT}}$-symmetric periodic finite-gap systems,  Phys. Rev. D 86, 085028 (2012).
  * [37] S. Dey and A. Fring,  Bohmian quantum trajectories from coherent states,  Phys. Rev. A 88, 022116 (2013).