# Dirac equation on a catenoid bridge: a supersymmetric approach

Ö. Yeşiltaş  Department of Physics, Faculty of Science, Gazi University, 06500
Ankara, Turkey<EMAIL_ADDRESS>J. Furtado  Universidade Federal do
Cariri(UFCA), Av. Tenente Raimundo Rocha,
Cidade Universitária, Juazeiro do Norte, Ceará, CEP 63048-080, Brasil  J.E.G.
Silva  Universidade Federal do Cariri(UFCA), Av. Tenente Raimundo Rocha,
Cidade Universitária, Juazeiro do Norte, Ceará, CEP 63048-080, Brasil

###### Abstract

In this paper we study the Dirac equation for an electron constrained to move
on a catenoid surface. We decoupled the two components of the spinor and
obtained two Klein-Gordon-like equations. Analytical solutions were obtained
using supersymetric quantum mechanics for two cases, namely, the constant
Fermi velocity and the position dependent Fermi velocity cases.

##  I Introduction

The geometry of the graphene layer geim ; novoselov ; katsnelson  plays an
important role in the electronic structure. The possibility of developing new
electronic devices due to the effect of graphene’s curvature is promoting a
growing interest in quantum mechanics at curved surfaces. Electronic
properties of a graphene nanotorus were studied considering the dynamics
governed by the Schrödinger equation GomesSilva:2020fxo , dainamle  as well as
by Dirac equation Yesiltas:2018zoy , dainam2 . Geometry fluctuations produce
pseudomagnetic fields ribbons  whose effects can be seen at ripples contijo
and corrugated layers corrugated . The curvature at the tip of a conical layer
produces a topological phase furtado , whereas helical strips induce chiral
properties dandoloff1 ; atanasovhelicoid ; atanasov .

Another interesting geometry studied in the last years is the catenoid
geometry. As a minimal surface, two dimensional wormhole geometry is
equivalent to a catenoid which are possible solutions for traversable
wormholes Wo . In Ref. wormhole ; picak  a bridge connecting a bilayer
graphene was proposed using a nanotube. In order to obtain a smooth bridge,
Ref.Dandoloff , dandoloff  suggested a catenoid surface to describe the
bilayer and the bridge using only one surface. This can be achieved due to the
catenoid curvature which is concentrated around the bridge and vanishes
asymptotically spivak . For a non-relativistic electron, the surface curvature
induces a geometric potential in the Schrödinger equation. The effects of the
geometry and external electric and magnetic fields upon the graphene catenoid
bridge was explored in Ref.euclides , where a single electron is governed by
the Schrödinger equation on the surface. Incidentally, the influence of a
position-dependent mass problem upon the electron on a catenoid bridge was
studied in Ref. euclides2 , where it was proposed an isotropic position-
dependent mass as a function of the Gaussian and mean curvatures.

The supersymmetry (SUSY) arose as an effort to obtain an unified description
of all interactions of nature wess , but it was immediately realized that it
provides powerful techniques to obtain analytical solutions for problems in
several branches of physics s1 . The factorization method was first introduced
by Dirac in order to obtain the spectrum of the quantum harmonic oscillator
dirac . Later, Schrödinger applied the same method for the radial part of the
Coulomb problem schrodinger . The supersymmetric quantum mechanics (SUSY QM)
is an intense field of research and nowadays there is a wide set of
Hamiltonians which are analytically solvable through the factorization method
sukhatme . On the other hand, there are more representations of SUSY QM such
as shape invariance R , Darboux transformations F  and Hamilton hierarchies DJ
. SUSY QM has received lot of interest for beautiful mathematical insight in
both relativistic and non-relativistic quantum mechanics 1 ; 2 ; 3 .

In this paper we study some properties of a relativistic electron constrained
to move on a catenoid surface. Unlike the non-relativistic electron, the Dirac
electron does not couple directly to the surface curvature. Instead, the
momentum operator (covariant derivative) is modified by the surface tetrads
and connection 1-forms. In section II, we study the massless Dirac equation on
the catenoid and we decoupled the two components of the spinor. As a result,
we obtain two second-order Klein-Gordon-like equations with an effective
potential depending on the surface curvature and on the effective Fermi
velocity. In section III we present briefly the SUSY QM approach which was
used to find analytical solutions for the Dirac equation on the catenoid in
two particular cases, namely, the constant Fermi velocity case and the
position-dependent Fermi velocity case. In section IV we present our final
remarks.

##  II Massless Dirac equation on a catenoid

The catenoid can be described by the meridian $u\in(-\infty,\infty)$ and the
parallel $\phi\in[0,2\pi)$, as

$\vec{r}=\sqrt{R^{2}+u^{2}}(\cos\phi\hat{i}+\sin\phi\hat{j})+R\sinh^{-1}\left(u/R\right)\hat{k},$
(1)

where $R$ is the radius of the catenoid bridge and $u=u(z)=R\sinh(z/R)$ , as
depicted in Fig.(1). Using this coordinate system, the $2+1$ spacetime
interval reads

$ds^{2}=dt^{2}-du^{2}-(R^{2}+u^{2})d\phi^{2}.$ (2)

Figure 1: Catenoid

Consider an orthogonal moving frame defined by the tangent vectors
$\hat{e}_{1}=\frac{\partial\vec{r}}{\partial u}$,
$\hat{e}_{2}=\frac{1}{\sqrt{R^{2}+u^{2}}}\frac{\partial\vec{r}}{\partial\phi}$
and the normal unit vector $\hat{e}_{3}=\hat{e}_{1}\times\hat{e}_{2}$. The
respective co-frame is given by $\omega^{0}=dt$, $\omega^{1}=du$,
$\omega^{2}=\sqrt{R^{2}+u^{2}}d\phi$ and $\omega^{3}=0$, where
$\omega^{a}=e^{a}_{\mu}dx^{\mu}$. Thus,
$d\vec{r}=\omega^{1}\otimes\hat{e}_{1}+\omega^{2}\otimes\hat{e}_{2}$ and the
metric can be rewritten as $ds^{2}=\eta_{ab}\omega^{a}\otimes\omega^{b}$,
being $\eta_{ab}=diag(+1,-1,-1)$. Using the torsion-free Cartan equation
$T^{a}=d\omega^{a}+\omega^{a}_{b}\wedge\omega^{b}=0$ leads to the one-form
connections $\omega^{a}_{b}=\Gamma_{cb}^{a}\omega^{c}$ of form
$\omega^{1}_{3}=\frac{R}{R^{2}+u^{2}}\omega^{1}$,
$\omega^{2}_{3}=-\frac{R}{R^{2}+u^{2}}\omega^{2}$ and
$\omega^{2}_{1}=\frac{u}{R^{2}+u^{2}}\omega^{2}$. Since
$d\hat{e}_{3}=\omega^{1}_{3}\vec{e}_{1}+\omega^{2}_{3}\vec{e}_{2}=K_{b}^{a}\omega^{a}\otimes\vec{e}_{b}$,
then the second fundamental form $K^{a}_{b}$ has components

$\displaystyle K^{a}_{b}$ $\displaystyle=$
$\displaystyle\frac{R}{R^{2}+u^{2}}\left(\begin{array}[]{cc}1&0\\\
0&-1\end{array}\right).$ (5)

Therefore, the catenoid has gaussian and mean curvatures

$\displaystyle K=det(K^{a}_{b})=-\frac{R^{2}}{(R^{2}+u^{2})^{2}}$ ,
$\displaystyle M=\frac{1}{2}K^{a}_{a}=0.$ (6)

Accordingly, the catenoid has negative gaussian curvature vanishing far from
the the bridge throat.

Intrinsically, only the $\omega^{2}_{1}$ is non zero. Using the metric
compactibility condition $\omega_{ab}+\omega_{ba}=0$ and using
$\omega^{12}=\omega^{12}_{\mu}dx^{\mu}$, then

$\displaystyle\omega^{12}_{\phi}=\frac{u}{\sqrt{R^{2}+u^{2}}}.$ (7)

The massless Dirac equation on a surface reads

$iv_{F}\hbar D_{\mu}\Psi=0,$ (8)

where $v_{F}$ is the Fermi velocity. The covariant spinor derivative is
$D_{\mu}=\partial_{\mu}-\Gamma_{\mu}$, where
$\Gamma_{\mu}=\frac{1}{4}\Gamma^{ab}_{\mu}\gamma_{a}\gamma_{b}$ is the spinor
connection. From Eq.(7) the only non-vanishing spinor connection component is
$\Gamma_{\phi}=\frac{1}{4}\Gamma^{12}_{\phi}\gamma_{1}\gamma_{2}$. In (2+1)-D
we can adopt the following representation for the flat $\gamma^{a}$ matrices
$\gamma_{0}=\sigma_{3}$, $\gamma_{1}=-i\sigma_{2}$ and
$\gamma_{2}=-i\sigma_{1}$. The curved $\gamma^{\mu}$ matrices are related to
the flat $\gamma^{a}$ matrices by $\gamma^{\mu}=e^{\mu}_{a}\gamma^{a}$, so
that,

$\displaystyle\gamma^{t}$ $\displaystyle=$ $\displaystyle
e^{t}_{0}\gamma^{0}=\gamma_{0}$ (9) $\displaystyle\gamma^{u}$ $\displaystyle=$
$\displaystyle\gamma^{1}$ (10) $\displaystyle\gamma^{\phi}$ $\displaystyle=$
$\displaystyle e_{2}^{\phi}\gamma^{2}=\frac{1}{\sqrt{R^{2}+u^{2}}}\gamma^{2},$
(11)

and the spinor connection is

$\Gamma_{\phi}=\frac{i}{2}\frac{u}{\sqrt{R^{2}+u^{2}}}\gamma_{0}.$ (12)

The massless Dirac equation takes the form

$\displaystyle i\hbar
v_{F}\left(\gamma^{0}\partial_{0}+\gamma^{1}\partial_{u}+\frac{1}{\sqrt{R^{2}+u^{2}}}\gamma^{2}\partial_{\phi}-i\frac{u}{2(R^{2}+u^{2})}\gamma^{2}\sigma_{3}\right)\Psi=0.$

Therefore, the catenoid geometry leads to a geometric potential as
$V_{1}(u)=\frac{u}{4(R^{2}+u^{2})}\gamma_{1}$, which vanishes asymptotically
and exhibits a parity odd behaviour near the throat. The Dirac equation can be
written as

$\displaystyle i\hbar\frac{\partial\Psi}{\partial t}=H_{D}\Psi,$ (14)

where the Hamiltonian is given by

$H_{D}=-i\hbar
v_{F}\left(\sigma_{1}\left(\partial_{u}-\frac{u}{2(R^{2}+u^{2})}\right)-\sigma_{2}\frac{1}{\sqrt{R^{2}+u^{2}}}\partial_{\phi}\right).$
(15)

The independence of the catenoid metric with respect to time and the angular
coordinate $\phi$ leads to the conservation of energy and of the angular
momentum along the direction $z$. Thus, the wave function can be written as
$\Psi(t,u,\phi)=e^{i\frac{E}{\hbar}t}e^{im\phi}\psi(u)$ what leads to the
Dirac hamiltonian

$\displaystyle H_{D}$ $\displaystyle=$ $\displaystyle-i\hbar
v_{F}\left(\begin{array}[]{cc}0&\partial_{u}-\frac{u}{2(R^{2}+u^{2})}-\frac{m}{\sqrt{R^{2}+u^{2}}}\\\
\partial_{u}-\frac{u}{2(R^{2}+u^{2})}+\frac{m}{\sqrt{R^{2}+u^{2}}}&0\end{array}\right).$
(18)

Using the Hamiltonian in eq.(18), the Dirac equation leads to coupled
equations for the two components of the spinor
$\psi=\left(\begin{array}[]{cc}\psi_{1}\\\ \psi_{2}\end{array}\right)$, as

$\displaystyle\left(\partial_{u}-\frac{u}{2(R^{2}+u^{2})}+\frac{m}{\sqrt{R^{2}+u^{2}}}\right)\psi_{1}$
$\displaystyle=$ $\displaystyle-i\frac{E}{\hbar v_{F}}\psi_{2}$ (20)
$\displaystyle\left(\partial_{u}-\frac{u}{2(R^{2}+u^{2})}-\frac{m}{\sqrt{R^{2}+u^{2}}}\right)\psi_{2}$
$\displaystyle=$ $\displaystyle-i\frac{E}{\hbar v_{F}}\psi_{1}.$ (21)

Decoupling the system of first-order equations in eq. (21) leads to

$\displaystyle-\psi_{1}^{\prime\prime}+\frac{u}{R^{2}+u^{2}}\psi_{1}^{\prime}+U_{1}\psi_{1}$
$\displaystyle=$ $\displaystyle\frac{E^{2}}{\hbar^{2}v_{F}^{2}}\psi_{1},$ (22)

where
$U_{1}=\left(\frac{u}{R^{2}+u^{2}}\right)^{\prime}-\left(\frac{u}{R^{2}+u^{2}}\right)^{2}+\left(\frac{m}{\sqrt{R^{2}+u^{2}}}\right)^{2}-\left(\frac{m}{\sqrt{R^{2}+u^{2}}}\right)^{\prime}$
and the prime stands for the derivative with respect to $u$. Employing the
change on the wave function of the form $\psi_{1}=(R^{2}+u^{2})^{1/4}\chi_{1}$
we obtain

$\displaystyle-\chi_{1}^{\prime\prime}+V_{eff,1}\chi_{1}$ $\displaystyle=$
$\displaystyle\frac{E^{2}}{\hbar^{2}v_{F}^{2}}\chi_{1},$ (23)

where the effective potential $V_{eff,1}$ is given by

$\displaystyle
V_{eff,1}(u)=\Bigg{[}\left(\frac{m}{\sqrt{R^{2}+u^{2}}}\right)^{2}-\left(\frac{m}{\sqrt{R^{2}+u^{2}}}\right)^{\prime}\Bigg{]}$
(24)

Analogously for the second component of the spinor we have:

$\displaystyle-\chi_{2}^{\prime\prime}(u)+V_{eff,2}(u)\chi_{2}=\frac{E^{2}}{\hbar^{2}v_{F}^{2}}\chi_{2}(u),$
(25)

where the potential $V_{eff,2}(u)$ is exactly like $V_{eff,1}(u)$ but under
the interchanging $m\rightarrow-m$. The above equations were obtained under
the consideration of a constant Fermi velocity $v_{F}$. The position dependent
Fermi velocity case will be addressed in the following subsection.

###  II.1 Position-dependent Fermi velocity

It is interesting to consider a position-dependent Fermi velocity, i.e.,
$v_{F}=v_{F}(u)$. The dependence of the Fermi velocity as a function only of
$u$ lies in the symmetry of the catenoid on the angular variable, so that no
dependence of $\phi$ is expected. Starting from (20) and (21), the process of
decoupling the two equations renders for $\psi_{1}$

$\displaystyle\psi_{1}^{\prime\prime}+\left[\Sigma(u)-\left(\frac{1}{v_{F}(u)}\right)^{\prime}v_{F}(u)+\Lambda(u)\right]\psi_{1}^{\prime}+\left[\Sigma(u)^{\prime}-\left(\frac{1}{v_{F}(u)}\right)^{\prime}v_{F}(u)\Sigma(u)+\Lambda(u)\Sigma(u)\right]\psi_{1}=-\frac{E^{2}}{\hbar^{2}v_{F}(u)^{2}}\psi_{1}$

where,

$\displaystyle\Sigma(u)$ $\displaystyle=$
$\displaystyle\frac{m}{\sqrt{R^{2}+u^{2}}}-\frac{u}{2(R^{2}+u^{2})},$ (27)
$\displaystyle\Lambda(u)$ $\displaystyle=$
$\displaystyle-\frac{m}{\sqrt{R^{2}+u^{2}}}-\frac{u}{2(R^{2}+u^{2})}.$ (28)

In order to obtain a Klein-Gordon-like equation we must perform the following
change of variables $\psi_{1}(u)=\kappa(u)\omega_{1}(u)$, imposing that
$\omega_{1}(u)$ obeys a Sturm-Liouville equation. Such imposition leads to a
condition on the function $\kappa(u)$ which states that

$\kappa(u)=\exp\left[\frac{1}{4}\ln(R^{2}+u^{2})-\frac{1}{2}\ln(v_{F}(u))\right].$
(29)

Hence, we obtain

$\displaystyle-\omega_{1}(u)^{\prime\prime}+U_{eff,1}(u)\omega_{1}(u)=\frac{E^{2}}{\hbar^{2}v_{F}(u)^{2}}\omega_{1}(u),$
(30)

where,

$\displaystyle U_{eff,1}(u)=V_{eff,1}(u)+\bar{V}_{eff,1}(u),$ (31)

being

$\bar{V}_{eff,1}(u)=-\frac{\left(v_{F}(u)^{\prime}\right)^{2}-2v_{F}(u)\left[\frac{2mv_{F}(u)^{\prime}}{\sqrt{R^{2}+u^{2}}}+v_{F}(u)^{\prime\prime}\right]}{4v_{F}(u)^{2}}.$
(32)

Note that the potential $\bar{V}_{eff,1}(u)$ contains only the contributions
from the position-dependent Fermi velocity, so that if we consider a constant
Fermi velocity, $\bar{V}_{eff,1}(u)=0$ and we recover eq.(24). Analogously,
for the second component of the spinor we have an equation similar to eq.(30)
but with a potential $U_{eff,2}(u)$ equals to $U_{eff,1}(u)$ under the
interchanging $m\rightarrow-m$.

##  III Supersymmetric Quantum Mechanics Approach

Let us give a brief review of some aspects of supersymmetric quantum mechanics
(SUSY QM) junker  that we use in our search of exact solutions for the
effective potential models given in the previous Sections. The techniques
based on the factorization method which takes place in SUSY QM helps to
identify the Hamiltonians with solvable potentials. First, considering the one
dimensional eigenvalue equation below

$\emph{H}~{}\psi(z)=\left(-\frac{d^{2}}{dz^{2}}+U(z)\right)\psi(z)=E\psi(z)$
(33)

where $U(z)$ is a real scalar function. If we find two operators which are
adjoint of each other and first order differental operators as

$\mathcal{A}=\frac{d}{dz}+W(z),~{}~{}~{}~{}\mathcal{A^{{\dagger}}}=-\frac{d}{dz}+W(z),$
(34)

where the superpotential $W(z)$ is a real function and element of the first
order operators. Using the product of these operators, we can introduce two
partner Hamiltonians defined by

$\displaystyle\emph{H}=\emph{H}_{1}$ $\displaystyle=$
$\displaystyle\mathcal{A}^{{\dagger}}\mathcal{A}=-\frac{d^{2}}{dz^{2}}+U_{1}(z)$
(35) $\displaystyle\emph{H}_{2}$ $\displaystyle=$
$\displaystyle\mathcal{A}\mathcal{A}^{{\dagger}}=-\frac{d^{2}}{dz^{2}}+U_{2}(z).$
(36)

We can call the eigenfunctions of the Hamiltonians $\emph{H}_{1}$,
$\emph{H}_{2}$ as $\psi^{(1)}_{n}$ and $\psi^{(2)}_{n}$ while eigenvalues as
$E^{(1)}_{n}$ and $E^{(2)}_{n}$ respectively. Then, one can introduce the
partner potentials in terms of superpotentials as

$U_{1}(z)=W(z)^{2}-W^{\prime}(z),~{}~{}~{}~{}U_{2}(z)=W(z)^{2}+W^{\prime}(z).$
(37)

These partner Hamiltonians are sharing the same energy eigenvalues except the
groundstate. In case of unbroken SUSY, the ground state is not degenerate with
zero energy $E^{(1)}_{0}$ which also leads to

$\mathcal{A}\psi^{(1)}_{0}(z)=0.$ (38)

One can express the groundstate wavefunction in terms of the superpotential
using (38)

$\psi^{(n)}_{0}(z)=C_{n}\exp\left(\int^{z}W(t)dt\right).$ (39)

When it comes to the discrete spectrum of $H_{2}$, the relationship between
the energies of these Hamiltonians are given by

$E^{(2)}_{n}=E^{(1)}_{n+1},~{}~{}~{}~{}E^{(1)}_{0}=0.$ (40)

The eigenfunctions of the partner Hamiltonians are linked as below

$\psi^{(2)}_{n}=\frac{1}{\sqrt{E^{(1)}_{n+1}}}\mathcal{A}\psi^{(1)}_{n+1},$
(41)
$\psi^{(1)}_{n+1}=\frac{1}{\sqrt{E^{(2)}_{n}}}\mathcal{A}^{{\dagger}}\psi^{(2)}_{n}.$
(42)

We remind that if $\psi^{(1)}_{n+1}$ is normalized, then, $\psi^{(2)}_{n}$ is
also normalized(non-relativistic case). The annihilation of the groundstate of
the operator $H_{1}$ means that there will be no SUSY partner for this state
only sukhatme .

###  III.1 Constant Fermi velocity case

Let us consider the effective potential model for the constant Fermi velocity
given below

$\displaystyle\emph{H}_{1}$ $\displaystyle=$
$\displaystyle-\frac{d^{2}}{du^{2}}+V_{1,eff}(u)$ (43) $\displaystyle
V_{1,eff}(u)$ $\displaystyle=$
$\displaystyle\frac{m^{2}}{R^{2}+u^{2}}+\frac{mu}{(R^{2}+u^{2})^{3/2}},$ (44)

and

$\displaystyle\emph{H}_{1}\chi^{(1)}_{n}=\epsilon^{2}\chi^{(1)}_{n}(u),$ (45)

where $\epsilon=\frac{E}{v_{F}}$. The corresponding superpotential of the
effective potential of the Hamiltonian $\emph{H}_{1}$ is

$W(u)=\frac{m}{\sqrt{R^{2}+u^{2}}}.$ (46)

Then, one can obtain the groundstate wavefunction as

$\chi^{(1)}_{0}=2^{-m}(u+\sqrt{R^{2}+u^{2}})^{-m}.$ (47)

On the other hand, the partner Hamiltonian is given by

$\displaystyle\emph{H}_{2}$ $\displaystyle=$
$\displaystyle-\frac{d^{2}}{du^{2}}+V_{2,eff}(u)$ (48) $\displaystyle
V_{2,eff}(u)$ $\displaystyle=$
$\displaystyle\frac{m^{2}}{R^{2}+u^{2}}-\frac{mu}{(R^{2}+u^{2})^{3/2}},$ (49)

where we can write the first order operators
$\mathcal{A},\mathcal{A}^{{\dagger}}$ as

$\displaystyle\mathcal{A}$ $\displaystyle=$
$\displaystyle\frac{d}{du}+\frac{m}{\sqrt{R^{2}+u^{2}}}$ (50)
$\displaystyle\mathcal{A}^{{\dagger}}$ $\displaystyle=$
$\displaystyle-\frac{d}{du}+\frac{m}{\sqrt{R^{2}+u^{2}}}.$ (51)

Next, we use a point transformation which is $u=R\tan x$ so that the equations
(50) and (51) become

$\displaystyle\mathcal{A}$ $\displaystyle=$
$\displaystyle\frac{1}{R}\cos^{2}x\frac{d}{dx}+\frac{m}{R}\cos x$ (52)
$\displaystyle\mathcal{A}^{{\dagger}}$ $\displaystyle=$
$\displaystyle-\frac{1}{R}\cos^{2}x\frac{d}{dx}+\frac{m}{R}\cos x+\frac{\sin
2x}{R}.$ (53)

We can obtain the transformed Hamiltonian $\emph{H}_{1}$ as

$\emph{H}_{1}=-\frac{\cos^{4}x}{R^{2}}\frac{d^{2}}{dx^{2}}+\frac{2\cos^{2}x\sin
2x}{R^{2}}\frac{d}{dx}+\frac{m\cos^{2}x(m+3\sin x)}{R^{2}}.$ (54)

Using $\chi^{(1)}_{n}(x)=\sec^{2}x\bar{\chi}^{(1)}_{n}(x)$ in (54), we obtain

$-\bar{\chi}^{\prime\prime}_{1}(x)+\xi\bar{\chi}_{1}(x)=0,$ (55)

with

$\xi=3m\sec x\tan
x+(m^{2}+2)\sec^{2}x-4-\frac{E^{2}R^{2}}{v^{2}_{F}}\sec^{4}x,$ (56)

where we use $\chi^{(1)}_{n}(x)\rightarrow\chi_{1}(x)$ for the sake of
simplicity. Let us analyze now two important cases, namely, the zero energy
case and the case around the throat of the catenoid.

####  III.1.1 Case I

In case of zero energy, i.e. $\epsilon=\frac{ER}{v_{F}}=0$, we are lead to

$-\bar{\chi}^{\prime\prime}_{1}(x)+(3m\sec x\tan
x+(m^{2}+2)\sec^{2}x-4)\bar{\chi}_{1}(x)=0,$ (57)

where we can relate it to the system known as trigonometric Scarf-I potential
in the literature sukhatme  whose eigenvalues and eigenfunctions are well
known. Thus, one can get the solutions for zero energy as:

$\displaystyle\chi^{(1)}_{n}(x)$ $\displaystyle=$ $\displaystyle
N_{1}~{}[1+\exp(2ix)]^{2}\exp\left\\{-2i[x-m\arctan\exp(ix)]\right\\}.$ (58)

Let us turn back to (55) and apply $r=\sin x$ and
$\bar{\chi}_{1}(r)=(r^{2}-1)^{1/4}\bar{\bar{\chi}}_{1}(r)$ to get

$\displaystyle-(1-r^{2})\bar{\bar{\chi}}_{1}^{\prime\prime}(r)+2r\bar{\bar{\chi}}^{\prime}_{1}(r)+\left(\frac{r^{2}-2}{4(1-r^{2})}+3m\frac{r}{(1-r^{2})}+\frac{m^{2}+2}{(1-r^{2})}-4-\frac{\epsilon^{2}}{(1-r^{2})^{2}}\right)\bar{\bar{\chi}}_{1}(r)=0.$
(59)

Because (59) is not the type of hypergeometric differential equation, we shall
look at the approximation methods.
Singularity analysis:
Let us consider (59) expressed below:

$\bar{\chi}^{\prime\prime}_{1}(r)-\frac{2r}{1-r^{2}}\bar{\chi}^{\prime}(r)-\frac{1}{(1-r^{2})^{2}}\left(\frac{r^{2}-2}{4}+3mr+(m^{2}+2)-4(1-r^{2})-\frac{\epsilon^{2}}{1-r^{2}}\right)\bar{\chi}=0$
(60)

where $p(r)=-\frac{2r}{1-r^{2}}$ and
$q(r)=-\frac{1}{(1-r^{2})^{2}}\left(\frac{r^{2}-2}{4}+3mr+m^{2}+2-4(1-r^{2})-\frac{\epsilon^{2}}{1-r^{2}}\right)$.
The singular points of (60) are $r_{0}=\pm 1$. If $(r-r_{0})p(r)$ and
$(r-r_{0})^{2}q(r)$ are both analytic at $r=r_{0}$, then, $r_{0}$ is a
’regular singular point’, if not, it becomes ’irregular singular point’.
Hence, $r_{0}=1$ and $r_{0}=-1$ are irregular singular points because $q(r)$
is analytic at $r_{0}=\pm 1$, but $q(r)$ is not analytic due to the term
$\frac{\epsilon^{2}}{1-r^{2}}$. If $r_{0}$ is an irregular singular point,
$(r-r_{0})^{k}p(x)$ and $(r-r_{0})^{2k}q(r)$ are analytic where $k$ is the
least integer number, then, the irregular singular point at $r=r_{0}$ has a
rank of $k-1$. In our problem, $k=2$ and the rank is $1$.
Next, we can continue to look for the solutions of (59).

####  III.1.2 Case II: $r\rightarrow 0$ limit:

Let us analyze the behavior of the solution of Eq.59 near the origin. By
expanding the coefficients of Eq.59 in power-series up to first order in $r$,
we obtain

$\chi^{\prime\prime}_{1}-2r\chi^{\prime}_{1}+\left(\frac{10-4m^{2}+4\epsilon^{2}}{4}-3mr\right)\chi_{1}=0,$
(61)

whose solution can be written as

$\chi_{1}(r)=c_{1}e^{-\frac{3mr}{2}}H_{\alpha}\left(\frac{3m}{2}+r\right)+c_{2}e^{-\frac{3mr}{2}}M\left(-\alpha/2,1/2,\left(\frac{3m}{2}+r\right)^{2}\right),$
(62)

where $\alpha=\frac{10+5m^{2}+4\epsilon^{2}}{8}$, $H_{\alpha}$ is the Hermite
polynomial of degree $\alpha$ and
$M\left(-\alpha/2,1/2,\left(\frac{3m}{2}+r\right)^{2}\right)$ is the Kummer
confluent hypergeometric function. The polynomial degree $\alpha$ is a
positive integer $\alpha=n$ provided that,

$\epsilon=\frac{\sqrt{8n-5(2+m^{2})}}{2}.$ (63)

A deep analysis on the spectrum near the catenoid throat can be accomplished
as the following. Let us look at the power series of
$\frac{\epsilon^{2}}{(1-r^{2})^{2}}$ about the point $r=0$. Up to first-order
in $r$, we have

$\frac{\epsilon^{2}}{(1-r^{2})^{2}}\approx\epsilon^{2}.$ (64)

Using the first term of (64) in (59), it becomes

$-(1-r^{2})\bar{\bar{\chi}}_{1}^{\prime\prime}(r)+2r\bar{\bar{\chi}}^{\prime}_{1}(r)+\left(\frac{r^{2}-2}{4(1-r^{2})}+3m\frac{r}{(1-r^{2})}+\frac{m^{2}+2}{(1-r^{2})}-4-\epsilon^{2}\right)\bar{\bar{\chi}}_{1}(r)=0.$
(65)

If we consider the functions of Sturm-Lioville type equation given in (59) and
(65),

$-(1-r^{2})\bar{\bar{\chi}}_{1}^{\prime\prime}(r)+2r\bar{\bar{\chi}}^{\prime}_{1}(r)+V_{i}(r)\bar{\bar{\chi}}_{1}=0,~{}~{}~{}~{}i=1,2$
(66)

where

$\displaystyle V_{1}(r)$ $\displaystyle=$
$\displaystyle\frac{r^{2}-2}{4(1-r^{2})}+3m\frac{r}{(1-r^{2})}+\frac{m^{2}+2}{(1-r^{2})}-4-\frac{\epsilon^{2}}{(1-r^{2})^{2}}$
(67) $\displaystyle V_{2}(r)$ $\displaystyle=$
$\displaystyle\frac{r^{2}-2}{4(1-r^{2})}+3m\frac{r}{(1-r^{2})}+\frac{m^{2}+2}{(1-r^{2})}-4-\epsilon^{2}$
(68)

In the low energies values, $V_{1}(r)$ and $V_{2}(r)$ have similar behaviour
while they act different in the higher energy values as is shown in the figure
below.

Figure 2: Graph of (67) and (68). $\epsilon=0.5$, $m=2$ are used for each
curve. The red curve corresponds to $V_{1}(r)$, green one is the curve of
$V_{2}(r)$.

.

Figure 3: Graph of (67) and (68). $\epsilon=15$, $m=2$ are used for each
curve. The red curve corresponds to $V_{1}(r)$, blue one is the curve of
$V_{2}(r)$.

.

On the other hand, (67) and (68) can be expanded in series when $r\rightarrow
0$ as

$\displaystyle V_{1}(r)$ $\displaystyle=$
$\displaystyle-\frac{5}{2}+m^{2}-\epsilon^{2}+3mr+(\frac{7}{4}+m^{2}-2\epsilon^{2})r^{2}+O(r^{3})$
(69) $\displaystyle V_{2}(r)$ $\displaystyle=$
$\displaystyle-\frac{5}{2}+m^{2}-\epsilon^{2}+3mr+(\frac{7}{4}+m^{2})r^{2}+O(r^{3}),$
(70)

the system reduced into the linear potential plus harmonic oscillator system.
Let us consider (66) with $V_{1}(r)$. A mapping
$\bar{\bar{\chi}}_{1}=\frac{Z(r)}{\sqrt{1-r^{2}}}$ can be applied to (66) and
dividing the whole equation by $1-r^{2}$ leads to

$-Z^{\prime\prime}(r)+U(r)Z(r)=0,$ (71)

where
$U(r)=\frac{14-4m^{2}-12mr+(-31+4m^{2})r^{2}+12mr^{3}+17r^{4}+4\epsilon^{2}}{4(-1+r^{2})^{3}}$.
We can expand $U(r)$ in series when $r\rightarrow 0$, then,

$-Z^{\prime\prime}(r)+\left(3mr+(-\frac{11}{4}+2m^{2}-3\epsilon^{2})r^{2}-\frac{7}{2}+m^{2}-\epsilon^{2}\right)Z(r)=0.$
(72)

(72) can also be written as

$\displaystyle-Z^{\prime\prime}(r)$ $\displaystyle+$ $\displaystyle
U_{1}(r)Z(r)=\epsilon^{2}Z(r),$ (73) $\displaystyle U_{1}$ $\displaystyle=$
$\displaystyle
3mr+\left(-\frac{11}{4}+2m^{2}-3\epsilon^{2}\right)r^{2}-\frac{7}{2}+m^{2}$
(74)

We note that (74) is known as energy dependent potentials in the literature
becasue of it’s coefficient including energy parameter edp . Moreover, it is
important that $2m^{2}>\frac{11}{4}+3\epsilon^{2}$ which makes the coefficient
of $r^{2}$ term as positive due to the physical solutions. The solutions of
(74) can be found as

$Z_{n}(r)=D_{n}\left(\frac{6m+8m^{2}r-(11+12\epsilon^{2})r}{(-11+8m^{2}-12\epsilon^{2})^{3/4}}\right)$
(75)

where $D_{n}(ar)$ are the parabolic cylinder functions and $\epsilon^{2}$
satisfies the relationship

$\epsilon^{2}=\frac{1}{2f(\epsilon)^{3}/4}\left(-16m^{4}+(11+12\epsilon^{2})(-7+f(\epsilon))-8m^{2}(-12-3\epsilon^{2}+f(\epsilon))\right)=n,$
(76)

here, $f(\epsilon)=\sqrt{-11+8m^{2}-12\epsilon^{2}}$. We are not giving the
solutions of (76) which are very huge, one can compute the numerical energy
eigenvalues. In the recent paper, our interest is fixing the $\epsilon^{2}$
term in (59) without sacrificing more originality. Hence, let us make another
transformation for the function $\bar{\bar{\chi}}_{1}(r)$ which leads to a
transformation of (65) into a Jacobi type differential equation

$\bar{\bar{\chi}}_{1}(r)=(r+1)^{b}(1-r)^{a}Y(r)$ (77)

then, (65) turns into

$\displaystyle(1-r^{2})Y^{\prime\prime}(r)-2\left(a-b+(1+a+b)r\right)Y^{\prime}(r)+[\frac{1}{4}(17-4a^{2}-4b^{2})+\epsilon^{2}$
$\displaystyle+$ (78)
$\displaystyle\frac{7-16a^{2}+12m+4m^{2}}{8(r-1)}+\frac{-7+16b^{2}+12m-4m^{2}}{8(r+1)}-b-a(1+2b)]Y(r)=0.$

(78) is in the form of Jacobi differential equation which is given by

$\displaystyle(1-r^{2})y^{\prime\prime}(r)+\left(\beta-\alpha-(\alpha+\beta+2)r\right)y^{\prime}(r)+n(n+\alpha+\beta+1)y(r)=0.$
(79)

Matching (78) and (79) gives

$\displaystyle a$ $\displaystyle=$
$\displaystyle\frac{1}{4}\sqrt{7+12m+4m^{2}},$ (80) $\displaystyle b$
$\displaystyle=$ $\displaystyle\frac{1}{4}\sqrt{7-12m+4m^{2}}.$ (81)

The parametrs $a,b$ makes the rational terms in (78) equal to zero. and we
also obtain,

$\displaystyle\alpha$ $\displaystyle=$ $\displaystyle 2a$ (82)
$\displaystyle\beta$ $\displaystyle=$ $\displaystyle 2b.$ (83)

$\bar{\bar{\chi}}_{1}(r)$ can be obtained as

$\bar{\bar{\chi}}_{1}(r)=(r+1)^{b}(1-r)^{a}P_{n}^{(2a,2b)}(r),$ (84)

where $P_{n}^{(2a,2b)}(r)$ are the Jacobi polynomials. Hence, we can give the
whole solutions of our system for $n=\nu$ (45)

$\displaystyle\chi^{(1)}_{n}(u)=N_{1}\left(1-\frac{u}{\sqrt{u^{2}+R^{2}}}\right)^{\frac{1}{4}\sqrt{7+12m+4m^{2}}-1}\left(\frac{u}{\sqrt{u^{2}+R^{2}}}+1\right)^{\frac{1}{4}\sqrt{7-12m+4m^{2}}-1}P^{(2a,2b)}_{n}\left(\frac{u}{\sqrt{u^{2}+R^{2}}}\right),$
(85)

$N_{1}$ is the corresponding normalization constant and the energy eigenvalues
can be found as

$E^{(1)}_{n}=\pm\frac{V_{F}}{2\sqrt{2}R}\left[-27+4m^{2}+2(M_{1}+M_{2})+M_{1}M_{2}+8n^{2}+4n(M_{1}+M_{2})\right]^{1/2},$
(86)

where $M_{1}=\sqrt{7+4m(m-3)}$ and $M_{2}=\sqrt{7+4m(m+3)}$. The results agree
with the solutions already obtained in the literature sukhatme . It is
important to highlight here that since the energy must be real, it imposes a
constraint condition between the principal quantum number $n$ and the orbital
angular momentum $m$, such that the inside of the square root should be
positive for physical solutions. In fig.(4a) we can see the probability
density and in fig.(4b) the probability ring for $n=1$ and $m=-2$. These
values for $n$ and $m$ were chosen in order to obey the constraint imposed by
the realness of the energy. From the plots we can see a region of confinement
of the electron around the throat of the catenoid. Similarly, for $n=3$ and
$m=-2$ we can see in fig.(5) that besides the confinement region at the throat
of the catenoid, we have also two other regions of confinement, one in the
upper sheet and another in the lower sheet.

We can also calculate the $\bar{\chi}^{(2)}_{n}(u)$ which is the solution of
(48) using (41):

$\displaystyle\chi^{(2)}_{n}$ $\displaystyle=$ $\displaystyle
N_{2}\frac{1}{\sqrt{E_{n+1}}}\left(\frac{d}{u}+\frac{m}{\sqrt{R^{2}+u^{2}}}\right)\chi^{(1)}_{n+1}(u)$

Considering the partner Hamiltonian $\mathcal{H}_{2}$, let us introduce the
corresponding transformed partner potential $V_{2,eff}(x)$ corresponding to
(48).

$V_{2,eff}(x)=m^{2}\sec^{2}x+m\tan x\sec x$ (87)

which is sharing the same energy spectrum with $\emph{H}_{1}$.

(a)

(b)

Figure 4: (a) Probability density for $\chi_{n}^{(1)}(u)$ when $n=1$ and
$m=-2$. (b) Probability ring for $\chi_{n}^{(1)}(u)$ when $n=1$ and $m=-2$

(a)

(b)

Figure 5: (a) Probability density for $\chi_{n}^{(1)}(u)$ when $n=3$ and
$m=-2$. (b) Probability rings for $\chi_{n}^{(1)}(u)$ when $n=3$ and $m=-2$.

###  III.2 position-dependent Fermi velocity case

We consider now the transformed equation for the position-dependent Fermi
velocity case presented in the previous Section

$-\omega^{\prime\prime}(u)+\left[\frac{m^{2}}{R^{2}+u^{2}}+\frac{mu}{(R^{2}+u^{2})^{3/2}}+\frac{((v_{F}(u))^{\prime})^{2}}{4v_{F}(u)^{2}}-\frac{1}{2v_{F}(u)}\left(v^{\prime\prime}_{F}(u)+\frac{2mv^{\prime}_{F}(u)}{\sqrt{R^{2}+u^{2}}}\right)\right]\omega(u)=\frac{E^{2}}{v_{F}(u)^{2}}\omega(u).$
(88)

Using a point transformation which is $u=R\tan x$ in (88) leads to

$-\omega^{\prime\prime}(x)+2\tan
x\omega^{\prime}(x)+\Delta(x)\omega(x)=\frac{E^{2}R^{2}}{v_{F}(x)^{2}}\omega(x),$
(89)

where

$\Delta(x)=m^{2}\sec^{2}x+m\tan x\sec x+\frac{R\sec
x}{4}\frac{v^{\prime}_{F}(x)}{v_{F}(x)}-\frac{v^{\prime\prime}_{F}(x)}{2v_{F}(x)}+\frac{v^{\prime}_{F}(x)}{v_{F}(x)}\tan
x-\frac{mR\sec x}{v_{F}(x)}v^{\prime}_{F}(x)$ (90)

Using $\omega(x)=\sec x\omega_{1}(x)$ and a suggestion on the Fermi velocity
which is

$v_{F}(x)=\lambda\sec^{2}x$ (91)

makes (88) as

$-\omega^{\prime\prime}_{1}(x)+\left[-1+(m^{2}-1)\sec^{2}x+\frac{1}{2}(2m+R-4mR)\sec
x\tan x\right]\omega_{1}(x)=\frac{E^{2}R^{2}}{\lambda^{2}}\omega_{1}(x),$ (92)

where $\lambda$ is a real parameter. The trigonometric ansatz for the Fermi
velocity could be related to periodic inhomogeneities in the graphene lattice,
such as the presence of heteroatoms or topological defects. The potential
function in (92) is known as trigonometric Scarf-I potential in the literature
whose solutions are very well-known sukhatme . Then, the solutions of (92) can
be written as

$\omega_{1}(x)=(1-\sin x)^{\frac{A-B}{2}}(1+\sin
x)^{\frac{A+B}{2}}P^{(A-B-1/2,A+B-1/2)}_{n}(\sin x),$ (93)

where

$\displaystyle A$ $\displaystyle=$
$\displaystyle\frac{1}{8(-R+m(4R-2))}\times$
$\displaystyle\times[-4R+8m(2R-1)+4\sqrt{2}m^{2}c-\sqrt{2}(3+c\sqrt{-3+4m^{2}(1+m-2R)(-3-2R+4m(-1+m+2R))})],$
$\displaystyle B$ $\displaystyle=$ $\displaystyle\frac{c}{2\sqrt{2}},$ (95)
$\displaystyle c$ $\displaystyle=$
$\displaystyle\sqrt{-3+4m^{{}^{2}}+\sqrt{(-3+4m(1+m-2R)+2R)(-3-2R+4m(-1+m+2R))}},$
(96)

and energy spectrum is given as

$E^{(1)}_{n}=\pm\frac{\lambda}{R}\sqrt{(A+n)^{2}-1},~{}~{}n=0,1,2,...$ (97)

Also in this case, the quantity $(A+n)^{2}-1$ must be positive definite in
order to assure that the energy is real. Such imposition leads to a constraint
condition between the $n$ and $m$. The solutions of $\omega$ can be written in
terms of $u$ as

$\displaystyle\omega_{n}(u)=\omega^{(1)}_{n}(u)$ $\displaystyle=$
$\displaystyle
N_{1}\sqrt{1+\frac{u^{2}}{R^{2}}}\left(1-\frac{u}{\sqrt{R^{2}+u^{2}}}\right)^{\frac{A-B}{2}}\left(1+\frac{u}{\sqrt{R^{2}+u^{2}}}\right)^{\frac{A+B}{2}}\times$
(98) $\displaystyle\times
P^{(A-B-1/2,A+B-1/2)}_{n}\left(\frac{u}{\sqrt{R^{2}+u^{2}}}\right).$

We also remind that the corresponding superpotential of the system in (92) is:

$W(x)=-A\tan x+B\sec x.$ (99)

We can introduce the partner Hamiltonian of the system given in (87) as

$-\omega_{2}^{\prime\prime}(u)+\left[\frac{m^{2}}{R^{2}+u^{2}}-\frac{mu}{(R^{2}+u^{2})^{3/2}}-\frac{((v_{F}(u))^{\prime})^{2}}{4v_{F}(u)^{2}}+\frac{1}{2v_{F}(u)}\left(v^{\prime\prime}_{F}(u)+\frac{2mv^{\prime}_{F}(u)}{\sqrt{R^{2}+u^{2}}}\right)\right]\omega_{2}(u)=\frac{E^{2}}{v_{F}(u)^{2}}\omega_{2}(u).$
(100)

Using (99) and (93), we can calculate the solutions of (100) which shares the
same energy states (97) with (87) except the ground state:

$\displaystyle\omega_{2}(x)$ $\displaystyle\sim$ $\displaystyle\sec
x\frac{(\cos\frac{x}{2}-\sin\frac{x}{2})(1-\sin
x)^{A-1}}{2(\cos\frac{x}{2}+\sin\frac{x}{2})}\left\\{(1+2A+n)\cos^{2}xP\left(n,\frac{1}{2}+A+B,\frac{1}{2}+A-B,\sin
x\right)\right.$ (101)
$\displaystyle\left.-2P\left(n+1,-\frac{1}{2}+A+B,-\frac{1}{2}+A-B,\sin
x\right)(A-B+2A\sin x)\right\\},$

And the whole solutions read as

$\displaystyle\omega_{2}(u)$ $\displaystyle=$
$\displaystyle\omega^{(2)}_{n}(u)=\frac{N_{2}}{E^{(1)}_{n+1}}\sqrt{1+\frac{u^{2}}{R^{2}}}\frac{\sqrt{1+\frac{R}{\sqrt{R^{2}+u^{2}}}}\left(1-\frac{u}{\sqrt{R^{2}+u^{2}}}\right)^{A-1}}{2\sqrt{1-\frac{R}{\sqrt{R^{2}+u^{2}}}}}\left\\{(1+2A+n)\frac{R}{\sqrt{R^{2}+u^{2}}}\right.\times$
(102) $\displaystyle\times
P\left(n,\frac{1}{2}+A+B,\frac{1}{2}+A-B,\frac{u}{\sqrt{R^{2}+u^{2}}}\right)-2P\left(n+1,-\frac{1}{2}+A+B,-\frac{1}{2}+A-B,\frac{u}{\sqrt{R^{2}+u^{2}}}\right)\times$
$\displaystyle\left.\times(A-B+2A\frac{u}{\sqrt{R^{2}+u^{2}}})\right\\}$

##  IV final remarks

We studied the Dirac equation for an electron constrained to move on a
catenoid surface realized as a connection bridge between two layers of
graphene. We decoupled the two components of the spinor and in order to obtain
two Klein-Gordon-like equations. Analytical solutions were obtained using
supersymmetric quantum mechanics for two cases, namely, the constant Fermi
velocity case and the position dependent Fermi velocity case. In case of
constant Fermi velocity, both zero energy modes and approximated bound state
energies are obtained. The bound state solutions for the constant Fermi
velocity are obtained in the limit of $r\rightarrow 0$.

For the constant Fermi velocity case, the supersymmetric quantum mechanical
approach allows us to identify the potential as a very well known
trigonometric Scarf-I potential. The eigenvalues of the Hamiltonian operator
and the eigenfunctions were found exactly. It is important to highlight that
the solutions of the system imposes a relation between the principal quantum
number $n$ and the orbital angular momentum $m$ in order to guarantee the
positively of the energy.

For the position-dependent Fermi velocity case, similarly, the supersymmetric
quantum mechanical approach allows us to identify the potential as a very well
known trigonometric Scarf-I potential. In this case exact solutions were also
found under a trigonometric ansatz for the Fermi velocity.

The paper not only presents important properties about the dynamics of an
electron constrained to move on a catenoid bridge but also opens up new
possibilities of investigation. The electron-phonon interaction as well as the
thermodynamic properties will be addressed in a future work.

Data Availability Statement
The authors confirm that the data supporting the findings of this study are
available within the article and its supplementary materials.

## Acknowledgments

J.E.G.Silva thanks the Conselho Nacional de Desenvolvimento Científico e
Tecnológico (CNPq), grants no 312356/2017-0 for financial support.

## References

  * (1) A. K. Geim, K. S. Novoselov, Nature Materials 6, 183 (2007).
  * (2) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  * (3) M. Katsnelson, Graphene: Carbon in two dimensions, Cambridge University Press, Cambridge, (2012).
  * (4) J. E. Gomes Silva, J. Furtado and A. C. A. Ramos, Eur. Phys. J. B 93, no.12, 225 (2020)
  * (5) Dai-Nam Le et al, Eur. Phys. Lett. 10005 (2019).
  * (6) Ö. Yeşiltaş, Adv. High Energy Phys. 2018, 6891402 (2018).
  * (7) Dai-Nam Le et al, Physica E 107 60 (2019).
  * (8) F. Guinea, A. K. Geim, M. I. Katsnelson, and K. S. Novoselov, Phys. Rev. B 81, 035408 (2010).
  * (9) F. de Juan, A. Cortijo, M. A. H. Vozmediano, Phys. Rev. B 76, 165409, (2007).
  * (10) V. Atanasov, A. Saxena, Phys. Rev. B 81, 205409 (2010).
  * (11) C. Furtado , F. Moraes , A.M. de M. Carvalho, Phys. Lett. A 372, 5368, (2008).
  * (12) R. Dandoloff, T. T. Truong, Phys. Lett. A 325, 233 (2004).
  * (13) V. Atanasov, R. Dandoloff and A. Saxena Phys. Rev. B 79, 033404 (2009).
  * (14) V. Atanasov, A. Saxena, Phys. Rev. B 92, 035440 (2015).
  * (15) R. Dandoloff, A. Saxena and B. Jensen, Phys. Rev. A 81 014102 (2010).
  * (16) J. González, J. Herrero, Nucl. Phys. B 825, 426 (2010).
  * (17) R. Pincak, J. Smotlacha, Quantum Matter 5, 114 (2016).
  * (18) R. Dandoloff, Phys. Lett. A, 373 (2009).
  * (19) R. Dandoloff, A. Saxena, B. Jensen, Phys. Rev. A 81, 014102 (2010).
  * (20) M. Spivak, A comprehensive introduction to differential geometry, Publish or Perish, Houston, (1999).
  * (21) J.E.G. Silva et al., Phys. Lett. A 384, 126458 (2020).
  * (22) J.E.G. Silva et al., Eur. Phys. J. B 94, 127 (2021).
  * (23) J. Wess and J. Bagger, _Supersymmetry and Supergravity_ , Princeton University Press, 1992 ( ISBN:0691085560 ,0691025304).
  * (24) A. Sinha and P. Roy, Int. J. Mod. Phys. A 21, 5807–5822 (2006).
  * (25) P.A.M. Dirac, The Principles of Quantum Mechanics, 2nd edn, Clarendon, Oxford (1935)
  * (26) E. Schrödinger, Proc. R. Irish Acad. 46 A 9-16 (1940)
  * (27) F. Cooper, A. Khare and U. Sukhatme, _Supersymmetry in Quantum Mechanics_ , World Scientific, 2001( ISBN:9810246129, 9789810246129).
  * (28) R. L. Hall and O. Yesiltas, Int. J. of Mod. Phys. E, 1910 1923 (2010).
  * (29) F. Correa and V. Jakubsky, Phys. Rev. A, 95(3) 033807 (2017).
  * (30) D. Bermudez and D. J. Fernandez, Symm. Int. and Geo.-Meth. and Appl., 7 025 (2011).
  * (31) G. Junker, Symmetry-Basel, 12(10) 1590 (2020).
  * (32) G. Levai, Czech. J. of Phys. 54(10) 1121 (2004).
  * (33) M. Znojil, Symmetry-Basel, 12(6) 892 (2020).
  * (34) G. Junker, Supersymmetric Methods in Quantum, Statistical and Solid State Physics, IOP Publishing Ltd; Enlarged and Revised edition (April 20, 2019), ISBN-13: 978-0750320245.
  * (35) J. Formanek, R. J. Lombard, J. Mares, Czech. J. of Physics, 54 289 (2004).