text:
\section{Introduction}

 At present more than 200
pulsars\footnote{\color{blue}\detokenize{
https://confluence.slac.stanford.edu/display/GLAMCOG/Public+List+of+LAT-Detected+Gamma-Ray+Pulsars}} are
known to emit $\gamma$-rays in the high-energy
(HE; E $\gtrsim$ 100\,MeV) band, and this number is
steadily growing~\citep{2013ApJS..208...17A}.
Understanding the emission mechanism of $\gamma$-ray pulsars and its relationship to radio
pulsar emission is one of the most pressing questions in high-energy astrophysics.

Since the launch of the {\it Fermi}\ satellite, the $\gamma$-ray data unambiguously point to an
outer-magnetospheric origin of the high-energy emission.  The observed spectra lack the
super-exponential cut-off expected in ``polar cap'' models~\citep{2008Sci...322.1221A}, and there are far
more radio-quiet $\gamma$-ray pulsars than are predicted in the polar cap scenario. There
are a few independent {\it geometrical} models of pulsar high-energy
emission that  assume enhanced particle acceleration and radiation along the last open
field
lines of a magnetic dipole~\citep{1995ApJ...438..314R, 2010ApJ...715.1282B}.
Though {\it geometrical} models of the $\gamma$-ray emission from pulsars have been quite
successful, the underlying physics of particle acceleration and emission remains an open
question.  In modeling  the electrodynamics of the pulsar gap, it is
typically assumed that the primary emission mechanism  is curvature
radiation~\citep[\eg][]{1996ApJ...470..469R,Baring99}.

Our current theoretical models of pulsar radiation are not yet able
 to produce solid  predictions.  Indeed, there are several unanswered questions:
the global kinetic structure of the magnetosphere, the location of the gaps,
the current closure, and the location of pair production \citep{2014ApJ...795L..22C,2015MNRAS.448..606C,2015MNRAS.449.2759B}.  Also, the process
is non-local---photons emitted at one place can produce an avalanche on different field
lines, so the gap may  not self-close even though it produces a cascade.

The VERITAS Collaboration  detected  pulsed emission from the Crab pulsar at energies above
100\,GeV~\citep{2011Sci...334...69V}. This  discovery was confirmed by MAGIC,
who have recently claimed the detection of the pulsed emission up to
1.5\,TeV~\citep{2016A&A...585A.133A}. Most importantly, the
photon spectrum of the pulsed emission between 100\,MeV and 400\,GeV is best described by
a broken power law and is statistically  preferred over a power law with an
exponential cut-off.

 Conventionally,   $\gamma$-ray emission from pulsars is attributed to incoherent
curvature radiation.
The detection of the Crab pulsar by VERITAS~\citep{2011Sci...334...69V}  implies the
importance of inverse-Compton (IC) scattering for the production of $\gamma$-ray photons.
This follows from a very general argument for an upper limit
on the possible energies  of curvature
photons~\citep[following a similar approach applicable to synchrotron
emission;][]{1996ApJ...457..253D,2010MNRAS.405.1809L}.
Assuming  that the  radius of curvature $R_c $ of the  {magnetic field}\ lines  is a fraction  $\xi$ of the light cylinder radius $R_c = \xi R_{LC}$ and
balancing radiative losses with acceleration by the {electric  field}\ (which is a fraction of the
{magnetic field}, $E=\eta B$), one can find {\it the maximum energy of curvature emission
within the Crab pulsar magnetosphere}:
$
\epsilon_{br}
 =
(150\, {\rm GeV}) \, \eta^{3/4} \sqrt{\xi} .
$ \citep{2012ApJ...754...33L}.
The  possibility  that the emission above the break energy is produced  as a tail to
the curvature emission  is excluded by the fact that
the spectrum of the Crab pulsar reported by  MAGIC and VERITAS {\it does not show the
exponential cut-off} indicative of radiation reaction-limited curvature emission. Phase-dependent variations of the break energy  at few GeVs should not affect the shape of the spectrum at
sufficiently high energies, \eg above 10\,GeV.

Importantly, it is the high-energy tails that carry information about the emission mechanism and thus deserve
special consideration. The two competing models of pulsar $\gamma$-ray emission predict
qualitatively different spectra well above the break energy (curvature emission predicts
an exponential cut-off, while inverse-Compton scattering predicts a power-law).

One of the problems of the broadband double power-law fits,  commonly used for a typical {\it Fermi}-LAT
spectrum of a bright $\gamma$-ray pulsar, is that  most of the errors accumulate due to the  {\it
arbitrary} parametrization of the spectral near the peak, where the two power laws connect.  A different parametrization of this spectral part
 can produce substantially different fits~\citep{2012ApJ...757...88L}.

\cite{2013MNRAS.431.2580L} showed
that  the broadband  spectrum of the Crab pulsar from UV to very high-energy
$\gamma$-rays---nearly ten decades in energy---can be reproduced  within the framework of
the
cyclotron self-Compton model \citep[see also][]{2015ApJ...811...63H}. Emission is produced  by two counter-streaming beams
within
the outer gaps, at distances above $\sim$20 NS radii. The outward moving beam  produces
UV--X-ray photons via Doppler-boosted cyclotron emission, and GeV photons by
inverse-Compton
scattering the cyclotron photons produced by the inward going beam. The scattering occurs
deep in the Klein-Nishina regime, whereby the  IC component provides a direct
measurement
of the particle distribution within the magnetosphere.   The required plasma multiplicity
is
high, on the order of $10^6$ to $10^7$, but is consistent with the average particle flux injected into
the
pulsar wind nebula. Klein-Nishina reduction in the scattering cross-section (and the
corresponding  reduction
of the electron energy loss rate)  allows the primary  leptons to be accelerated to very
high energies with hard spectra.  The secondary plasma is less energetic but more dense
and has approximately the same energy content as the primary beam, producing
cyclotron self-Compton radiation. The  cyclotron component has a broad peak in the
UV--X-ray range, while the IC component extends to hundreds of GeV.  Below $\sim$1\,GeV
the  curvature emission from the primary beam can contribute a substantial flux fraction.
The IC emission from the primary beam extends well into the TeV regime but has proven
difficult to detect by Cherenkov telescopes.

How special is the Crab in terms of the high-energy power-law tail?
\cite{2012ApJ...757...88L} analyzed  the high energy emission from Geminga (using
post-processed data)  and concluded that the high energy part of the spectrum is better
fit with a power-law shape than with an exponential cut-off.  The non-detection by
VERITAS~\citep{2015ApJ...800...61A} is consistent with the fairly steep power-law slope
inferred
by \cite{2012ApJ...757...88L}.  \cite{2014ApJ...797L..13L} have reported the
detection of pulsed emission from Vela $> 50$\,GeV, noting that the observed
spectrum above 10\,GeV of the
Vela pulsar is harder than a simple exponential function predicted in the exponential
cut-off spectral shape of the tail.
\cite{2015ApJ...804...86M} has performed a statistical analysis of high
energy tails using a stacked analysis of 115 {\it Fermi}-LAT-detected
pulsars, finding no significant stacked excess at energies above 50\,GeV.

In this paper we present the results of a complementary analysis of high-energy pulsar
tails, defined as emission above 10\,GeV.  We chose  four bright $\gamma$-ray pulsars: PSR
J0007+7303, PSR J0633+1746 (Geminga), PSR J0835--4510 (Vela) and  PSR J2021+3651,
and we perform detailed spectral fits for each pulsar.

\section{Fermi-LAT Data Analysis \& Results}

The {\it Fermi}-LAT is an electron-positron pair-conversion telescope sensitive to gamma-ray photons with energies
between 20\,MeV and 300\,GeV.  It has a FoV of $\sim2.5$\,sr and attains full-sky coverage approximately every
three hours.  For a complete description of the instrument, see~\cite{2009ApJ...697.1071A,2012ApJS..203....4A}.

To obtain the spectral parameters for each of the four pulsars, the {\it Fermi}-LAT \texttt{Science Tools
v10r0p5} with Pass 8 reprocessed instrument response functions and the standard quality cuts described
in~\cite{2012ApJS..199...31N} are used. Roughly 7\,yr of ``source''-class events with energies between 100\,MeV
and 100\,GeV collected between 4 August 2008 and 12 August 2015 within a $20^\circ$ region-of-interest (ROI) of
the 2PC location of each pulsar are processed with the maximum likelihood fitting routine.

The spectral analysis for each pulsar except Vela
presented here follows the same steps outlined in the second {\it Fermi}-LAT pulsar
catalog~\citep[2PC][]{2013ApJS..208...17A}.
To generate spectral energy distributions of the {\it Fermi}-LAT data for each pulsar, a
binned maximum likelihood
analysis is done in logarithmically spaced energy bins spanning
the range 100\,MeV to 100\,GeV. Spectral models
for all sources in the 3FGL catalog~\citep{2015ApJS..218...23A} in the ROI in addition to the galactic and
isotropic diffuse backgrounds (gll\_iem\_v06.fits, iso\_P8R2\_SOURCE\_V6\_v06.txt) are included in the likelihood
fitting.  Sources within a circle $4^{\circ}$ in radius are modeled as power laws of the form
\begin{equation}
 {dN \over dE} = N_0 {\left(E \over E_0\right)}^{\gamma},
\end{equation}
where $N_0$ is the flux normalization, $E_0$ is fixed to 300\,GeV, and $\gamma$ is the spectral index.
These sources have only their normalization parameter left free in the fitting
routine.  All other parameters are fixed except for the normalization on the galactic and isotropic diffuse
background models.  In each bin, the pulsar is modeled as a point source with a simple power-law spectral shape
with a fixed spectral index of $-2$.  An iterative spectral binning is performed starting
at three bins per decade in energy and increasing to four bins per decade, five, and so
on. The spectrum containing the greatest number of bins above 10\,GeV is kept for
subsequent fitting, while the others are discarded.

The {\it Fermi}-LAT data analysis for the Vela pulsar is done differently than described above due to three nearby
bright diffuse sources (Vela X, Vela Jr., and Puppis A) that must be taken into account in order to achieve the
greatest sensitivity.  To extract spectral parameters for the diffuse sources, the data are first phase-folded
with the the \texttt{Tempo2} pulsar timing package~\citep{2006MNRAS.369..655H} using a publicly available timing
solution from a webpage maintained by M.
Kerr\footnote{\texttt{\detokenize{www.slac.stanford.edu/~kerrm/fermi_pulsar_timing/}}}.
Phase-folding the data allows
selection of the off-pulse region to remove the contribution from the Vela pulsar, which allows more robust
likelihood fitting for the nearby diffuse sources.   Since the timing solution is not valid for the full 7\,yr
data set, only data recorded by the {\it Fermi}-LAT during the period of validity of the timing solution are used
in the following analysis.  First, the off-pulse phase region between 0.8 -- 1.0 is selected to remove the Vela
pulsar.  Subsequently, all normalization parameters in the model are scaled to 20\% of their nominal values, and
the Vela pulsar model is removed.  All sources within $4^{\circ}$ of the Vela pulsar 2PC location are modeled as
power laws with free normalization parameters, including the three diffuse sources and the galactic and isotropic
diffuse components.  Likelihood analyses are done in logarithmically spaced energy bins to allow computation of
the normalization parameters for each source in each bin.  The normalization parameters are then scaled by a
factor of 5, and the spectral model for the Vela pulsar is added.  The normalization parameters for the three
diffuse sources are then fixed to the previously computed values in each energy bin, and the Vela pulsar SED is
generated in the same way as described in the previous paragraph.

Phase-averaged spectral energy distributions for each pulsar are shown in
Figure~\ref{seds}.  Each SED is fitted
with both a power-law and an exponential cut-off above 10\,GeV to test the preferred
spectral shape.  The choice of 10\,GeV is an arbitrary starting point for probing
 the highest energies of the spectra.  However, 10\,GeV is above a few times the
break energy for each pulsar reported in the 2PC, which helps to reduce contamination by
the harder portion of the spectra at lower energies.  Spectral
indices for
the power-law fits and fit probabilities are summarized in Table~\ref{tab:fit_summary}.
\newline
\newcolumntype{C}{>{\centering\arraybackslash}X}
\begin{table*}[h]
\centering
\begin{tabularx}{1.0\textwidth}{CCCCC}
  \hline
  Pulsar & Power-law spectral index above 10\,GeV & Power-law fit chi-square / ndf (probability) & Exponential
cut-off fit chi-square / ndf (probability) & Notes
\\ \hline
  PSR J0007+7303 & $-3.69 \pm 0.016$ & 5.26 / 3 (0.15) & 0.95 / 3 (0.81) & \\
  PSR J0633+1746 & $-5.12 \pm 0.16$ & 5.54 / 2 (0.06) & 4.87 / 2 (0.09) & Geminga \\
  PSR J0835--4510 & $-4.54 \pm 0.08$ & 14.1 / 5 (0.02) & 8.06 / 5 (0.15) & Vela \\
  PSR J2021+3651 & $-4.73 \pm 0.54$ & 1.10 / 1 (0.29) & 2.40 / 1 (0.12) & \\
  \hline
\end{tabularx}
\caption{ Spectral indices for the power-law fits, reduced chi-squares, and fit probabilities for the SED of each pulsar. }
\label{tab:fit_summary}
\end{table*}

\begin{figure*}[h]
\center
\epsscale{1.0}
\includegraphics[width=6.6in]{combined_4panel_SED_v4.pdf}
\caption{ Spectral energy distributions derived from the {\it Fermi}-LAT data for the four
pulsars: PSR J0007+7303, PSR J0633+1746 (Geminga), PSR J0835--4510 (Vela), and PSR
J2021+3651. In each panel, the solid line shows the power-law fit, while the dashed line
shows the exponential cut-off fit. Each fit is applied above 10\,GeV. \label{seds}}
\end{figure*}

The fit probabilities for either the power-law fit or the fit with an exponential cut-off are comparable. Therefore, one has to conclude that both models describe the data equally well and,
therefore, both models are equally preferred. The probability of a fit to the Vela data with a power law is the worst of all four power law fits. In the case of Vela, the data points at
the lowest energies included in the fit contribute the most to the $\chi^2$. In that energy range the spectrum shows clear curvature, highlighting the problem of arbitrarily
choosing 10\,GeV as the energy threshold for the fit. It also emphasizes that the sensitivity of the method strongly depends on how far away from the break in the spectrum the fit
begins.

In order to alleviate the subjective choice of where to start the fit, the spectra are also fitted with a power law plus sub-exponential cutoff (PLSEC) of the form
\begin{equation}
 {dN \over dE} = N_0 {\left(E \over E_0\right)}^{\gamma_1} \textrm{exp}\left(-{\left(E \over E_c\right)}^{\gamma_2}\right),
\end{equation}
where $E_0$ is fixed to 300\,GeV, and $N_0$, $E_c$, $\gamma_1 $ and $\gamma_2$ are free parameters, and a smooth, broken power
law (SBPL) of the form
\begin{equation}
 {dN \over dE} = N_0 {\left(E \over E_0\right)}^{\gamma_1}\left( 1 + {\left(E \over E_b	\right)}^{\gamma_1 - \gamma_2 \over \beta} \right)^{-\beta},
\end{equation}
where $E_0$ is again fixed to 300\,GeV, and $N_0$, $\gamma_1$, $\gamma_2$, $E_b$, and $\beta$ are all left free in the fit.
The parameters of the fit functions are constrained by fitting them to
data points below 10\,GeV. For both the PLSEC and the SBPL, including the minimum number of data points below 10\,GeV to constrain the parameters of the fit yields acceptable fit
probabilities. However, including more data points at lower energies in the fits, the fits fail with probabilities quickly sinking below $10^{-5}$.
For all energy ranges for the fits, the PLSEC clearly under-predicts the flux above 10\,GeV, and the SBPL law clearly over-predicts the
data---neither of the fitted functions match the data above 10\,GeV.

We note that the upper limit at 48.4\,GeV for PSR J0007+7303 is a factor of $\sim$3.5 below the power-law fit, but it still lies within 1.2\,$\sigma$ of the fit taking the error on the fit
into account.
Therefore, the upper limit at 48.4\,GeV is not inconsistent with the power-law scenario.

\section{Discussion}

We have performed an analysis of the high energy tails, defined as emission above 10\,GeV,
of four of the brightest {\it Fermi}-LAT-detected $\gamma$-ray pulsars. Inn our
approach, the binning was adjusted iteratively starting at three bins per decade
and increasing, keeping the spectrum with the most bins above 10\,GeV for subsequent
fitting. In the second {\it Fermi}-LAT pulsar
catalog~\citep[2PC][]{2013ApJS..208...17A} the binning was arbitrary, mostly tailored to
weaker sources. As a result, for the bright sources considered in this work, the binning
was not optimal.

The fit results for all four of the pulsars analyzed in this study do not allow any
firm claim  to be made about the shape of the spectra above 10\,GeV.  In no case is the
power-law fit or exponential cut-off fit significantly preferred over the other, though
the exponential cut-off is marginally preferred in the case of PSR J0007+7303 (which is
the dimmest pulsar of all considered in this work).
\iffalse
Previously, in the case of Geminga,
\cite{2012ApJ...757...88L} used post-processed data, with the standard binning, while in
the present work binning was adjusted to produce maximum likelihood of a fit.
\fi

One of the unquantified biases in our analysis is the choice of the lower energy cut-off,
taken here to be 10\,GeV. This choice  is motivated by our attempt to avoid fitting the
roll-off region near the break energies.

\iffalse
Alternatively we chose the lower cut at three
times the peak energy. In the sample, PSR J0007+7303 has the highest peak energy of  $4.7$
GeV and is the weakest of all, so a choice of higher energy cut reduces the statistical
significance.   The break energies in Vela and PSR J2021+3651  are  both  $3.0$ GeV.  The
only exception is Geminga, with the break energy of $2.2$ GeV.
\fi

Our work leaves open the question of the shape of the high-energy gamma-ray spectral
tails of these four pulsars.

\acknowledgments

\bibliographystyle{apj}


text:
\section{Introduction}

In addition of being a descriptive model for some physical phenomena [1-36]
(including but not limited to, energy density many-body problem, electronic
properties of semiconductors), the position-dependent-mass (PDM) von Roos
Hamiltonian [33] (in $\hbar =2m_{\circ }=1$ units)%
\begin{equation}
H=-\frac{1}{2}\left[ m\left( x\right) ^{\alpha }\partial _{x}m\left(
x\right) ^{\beta }\partial _{x}m\left( x\right) ^{\gamma }+m\left( x\right)
^{\gamma }\partial _{x}m\left( x\right) ^{\beta }\partial _{x}m\left(
x\right) ^{\alpha }\right] +V\left( x\right)
\end{equation}%
is shown to be a mathematically useful model as it enriches the class of
exactly-solvable quantum mechanical problems. Whilst the parameters $\alpha
,\beta ,$ and $\gamma $ are subjected to the so-called von Roos constraint%
\begin{equation}
\alpha +\beta +\gamma =-1\text{ };\text{ \ }\alpha ,\beta ,\gamma \in
\mathbb{R}
,
\end{equation}%
the only feasibly admissible case that ensures the continuity conditions at
the heterojunction boundaries between two crystals is for $\alpha =\gamma $.
This, in effect, reduces the domain of the acceptable values of the
so-called \emph{"ambiguity parameters"} $\alpha ,\beta ,$ and $\gamma $ and
suggests that the PDM-Hamiltonian (1) be rewritten as%
\begin{equation}
H=-m\left( x\right) ^{\alpha }\partial _{x}m\left( x\right) ^{\beta
}\partial _{x}m\left( x\right) ^{\alpha }+V\left( x\right) \text{ };\text{ \
}2\alpha +\beta =-1.
\end{equation}

Recently, moreover, Dutra and Almeida [14] have carried out a reliability
test on the orderings available in the literature. They have used an exactly
solvable Morse model and concluded that the orderings of Gora and Williams
[15] ($a=\beta =\gamma =0,$ $\alpha =-1$), and Ben Danial and Duke [16] ($%
a=\alpha =\gamma =0,$ $\beta =-1$) should be discarded for they result in
complex energies. Nevertheless, they have classified the ordering of Zhu and
Kroemer [17] ($a=0,$ $\alpha =\gamma =-1/2,$ $\beta =0$), and that of Li and
Kuhn [18] ($a=\alpha =0,\beta =\gamma =-1/2$) as good orderings. Therefore,
the continuity conditions at the heterojunction boundaries and Dutra's and
Almeida's [14] reliability test would ultimately single out Zhu and Kroemer (%
$a=0,$ $\alpha =\gamma =-1/2,$ $\beta =0$) as a \emph{"reliable good
ordering"} . Fitting into this category/classification, Mustafa and
Mazharimousavi [19] have used a PDM-pseudo-momentum operator and suggested a
new \emph{"reliable good ordering"} ($\beta =-1/2,$ $\alpha =\gamma =-1/4$).
Yet, in their study of classical and quantum PDM harmonic oscillator, Cruz
et al [20] have considered different ambiguity parameters settings and
argued that only one of these orderings, $\beta =-1/2,$ $\alpha =\gamma
=-1/4 $ (i.e., that of Mustafa and Mazharimousavi [19]), gives rise to a
potential term that is the same as the classical PDM oscillator. For more
details on the ambiguity associated with the uniqueness of the kinetic
energy operator, the reader my refer to [8, 14, 19, 20 and references cited
therein].

In this letter, only for mathematical and/or quantum mechanical curiosity on
the associated exact-solvability, we consider the PDM von Roos Hamiltonian
(1) for a quantum particle moving within the domain mandated by its own
PDM-function $m\left( x\right) $ accompanied by an infinite potential well $%
V\left( x\right) $. In section 2, we consider an unusual PDM-function of the
form [20]%
\begin{equation}
m\left( x\right) =\frac{m_{\circ }}{\left[ 1-\left( \frac{x}{a}\right) ^{2}%
\right] ^{2}}\text{ };\text{ }V\left( x\right) =\left\{
\begin{tabular}{l}
$0$ \ for $\left\vert x\right\vert <a$ \\
$\infty $ for $\left\vert x\right\vert \geq a$%
\end{tabular}%
\right. .
\end{equation}%
Under such settings, it is obvious that the classical motion of such a
particle is confined to the domain $\mathcal{D}\left( x\right) =\left(
-a,a\right) $. and the PDM-Hamiltonian in (1) would result a PDM-Schr\"{o}%
dinger equation of the form%
\begin{equation}
\left[ -\partial _{x}\frac{1}{m\left( x\right) }\partial _{x}+\tilde{V}%
\left( x\right) \right] \psi \left( x\right) =E\text{ }\psi \left( x\right)
\text{ };\text{ \ }\left\vert x\right\vert <a,
\end{equation}%
with%
\begin{equation}
\tilde{V}\left( x\right) =\frac{1}{2}\left( 1+\beta \right) \frac{m^{\prime
\prime }\left( x\right) }{m\left( x\right) ^{2}}-\left[ \alpha \left( \alpha
+\beta +1\right) +\beta +1\right] \frac{m^{\prime }\left( x\right) ^{2}}{%
m\left( x\right) ^{3}}+V\left( x\right) .
\end{equation}%
Where primes denote derivatives with respect to $x$ and the PDM-Schr\"{o}%
dinger equation in (5) is known as the target equation. Obviously, as $%
\alpha ,\beta ,$ and $\gamma $ change within constraint (2), a profile
change of $\tilde{V}\left( x\right) $ in (6) is unavoidable (hence, an
ordering ambiguity conflict erupts in the process). Consequences of a
point-canonical-transformation (PCT) mapping (often mediates a transition
between two different effective potentials) on such a PDM particle are also
discussed in section 2. Therein, we shall witness that such a mass setting
results in an effective reflectionless modified P\"{o}schl-Teller type
potential well (cf. e.g., [20, 37]). We conclude in section 3.

\section{Consequences of PCT-mapping}

Following the well known PCT recipe (cf., e.g., [19,24;27]) would, with a
substitution of the form $\psi \left( x\right) =m\left( x\right)
^{1/4}\,\phi \left( q\left( x\right) \right) $ in (5), result in $q^{\prime
}\left( x\right) =\sqrt{m\left( x\right) }$ and suggests the following point
canonical transformation%
\begin{equation}
q\left( x\right) =\int\nolimits^{x}\sqrt{m\left( t\right) }dt=a\tanh
^{-1}\left( \frac{x}{a}\right) \text{ ; \ \ }\left\vert x\right\vert <a,
\end{equation}%
with $m\left( x\right) $ given in (4). Consequently, the PDM-Schr\"{o}dinger
equation (5) is mapped into%
\begin{equation}
\left( -\frac{d^{2}}{dq^{2}}+V_{eff}\left( q\right) \right) \phi \left(
q\right) =E\phi \left( q\right) \text{ ;\ \ \ }q\in \left( -\infty ,\infty
\right) ,
\end{equation}%
where%
\begin{equation}
V_{eff}\left( q\right) =g_{1}\frac{m^{\prime \prime }\left( x\right) }{%
m\left( x\right) ^{2}}-g_{2}\frac{m^{\prime }\left( x\right) ^{2}}{m\left(
x\right) ^{3}},
\end{equation}%
with
\begin{equation}
\text{\ }g_{1}=\frac{1}{4}\left( 1+2\beta \right) \text{ },\text{ \ \ }%
g_{2}=\alpha \left( \alpha +\beta +1\right) +\beta +\frac{9}{16}.
\end{equation}%
One may then use (7) and substitute%
\begin{equation}
\frac{x}{a}=\tanh \left( \frac{q}{a}\right) \text{ ;\ \ \ }\left\vert
x\right\vert <a,
\end{equation}%
in (9) to imply%
\begin{equation}
V_{eff}\left( q\right) =-\frac{\lambda \left( \lambda -1\right) }{a^{2}\cosh
^{2}\left( q/a\right) }+\frac{8}{a^{2}}\left( 3g_{1}-2g_{2}\right) \text{ ; }%
\lambda \left( \lambda -1\right) =4\left( 5g_{1}-4g_{2}\right) .
\end{equation}%
Which, in turn, would lead to an over simplified reference Schr\"{o}dinger
equation of the form%
\begin{equation}
\frac{d^{2}\phi \left( q\right) }{dq^{2}}+\left( k^{2}+\mu ^{2}\frac{\lambda
\left( \lambda -1\right) }{\cosh ^{2}\left( \mu q\right) }\right) \phi
\left( q\right) =0\text{,}
\end{equation}%
where%
\begin{equation}
k^{2}=E-8\mu ^{2}\left( 3g_{1}-2g_{2}\right) =\mathcal{E}\text{ ; }\mu =1/a%
\text{\ .}
\end{equation}%
Remarkably, the "shifted-by-a-constant" potential well in (12) is known to
be reflectionless/transparent (at any energy) with a reflection coefficient
equals zero if $\lambda >1$ is a positive integer (cf., e.g., Lekner [21]
and/or Diaz et al. [25] for a comprehensive and detailed study on this
potential). Therefore, a quick recollection to the essentials of the
Modified P\"{o}schl-Teller potential by Diaz et al. [25] would (with a bound
state solution when $\mathcal{E}<0$ ) provide a clear correspondence to the
following recycled solution%
\begin{eqnarray}
\phi \left( q\right) &=&\left( \cosh \mu q\right) ^{\lambda }\left[
A\,_{2}F_{1}\left( a,b;\frac{1}{2};-\sinh ^{2}\mu q\right) +\right.  \notag
\\
&&\left. B\,\left( \sinh \mu q\right) \,_{2}F_{1}\left( a+\frac{1}{2},b+%
\frac{1}{2};\frac{3}{2};-\sinh ^{2}\mu q\right) \right]
\end{eqnarray}%
where
\begin{equation}
a=\frac{1}{2}\left( \lambda -\frac{\sqrt{\left\vert \mathcal{E}\right\vert }%
}{\mu }\right) \,;\text{ \ \ \ }b=\frac{1}{2}\left( \lambda +\frac{\sqrt{%
\left\vert \mathcal{E}\right\vert }}{\mu }\right) ,
\end{equation}%
and%
\begin{equation}
\mathcal{E}_{n}\mathcal{=-}\mu ^{2}\left( \lambda -1-n\right) ^{2};\text{ \
\ }n\in
\mathbb{N}
,\text{ }0\leq n<\lambda -1.
\end{equation}%
Hereby, as long as the reference Schr\"{o}dinger equation (13) is in point,
the corresponding solutions hold true. However, a mapping to the original
target problem defined through equations (4)-(6) would suggest that%
\begin{eqnarray}
\psi _{n}\left( x\right) &=&\left( \frac{1}{1-\left( \frac{x}{a}\right) ^{2}}%
\right) ^{\frac{\lambda +1}{2}}\left[ A\,_{2}F_{1}\left( a,b;\frac{1}{2};%
\frac{\left( \frac{x}{a}\right) ^{2}}{\left( \frac{x}{a}\right) ^{2}-1}%
\right) \right.  \notag \\
&&\left. +B\left( \sqrt{\frac{\left( \frac{x}{a}\right) ^{2}}{1-\left( \frac{%
x}{a}\right) ^{2}}}\right) \,_{2}F_{1}\left( a+\frac{1}{2},b+\frac{1}{2};%
\frac{3}{2};\frac{\left( \frac{x}{a}\right) ^{2}}{\left( \frac{x}{a}\right)
^{2}-1}\right) \right] ,
\end{eqnarray}%
and%
\begin{equation}
E_{n}=\mu ^{2}\left( 8\left( 3g_{1}-2g_{2}\right) -\left( \lambda
-1-n\right) ^{2}\right) ,\text{ \ \ }n\in
\mathbb{N}
,\text{ }0\leq n<\lambda -1,
\end{equation}%
with%
\begin{equation}
a=\frac{1}{2}\left( 1+n\right) \text{ };\text{ \ \ }b=\lambda -\frac{1}{2}%
\left( 1+n\right) .
\end{equation}

The ground state solution of which would (with $n=0$) read%
\begin{eqnarray}
\psi _{0}\left( x\right) &=&A_{0}\left( \frac{1}{1-\left( \frac{x}{a}\right)
^{2}}\right) ^{\left( \lambda +1\right) /2}F\left( \frac{1}{2},\lambda -%
\frac{1}{2};\frac{1}{2};\frac{\left( \frac{x}{a}\right) ^{2}}{\left( \frac{x%
}{a}\right) ^{2}-1}\right)  \notag \\
&=&A_{0}\left( \frac{1}{1-\left( \frac{x}{a}\right) ^{2}}\right) ^{\left(
2-\lambda \right) /2},
\end{eqnarray}%
with an eigenvalue%
\begin{equation}
E_{0}=\mu ^{2}\left( 8\left( 3g_{1}-2g_{2}\right) -\left( \lambda -1\right)
^{2}\right) .
\end{equation}%
Nonetheless, the boundary condition on the wave function%
\begin{equation*}
\psi \left( x\right) =m\left( x\right) ^{1/4}\phi \left( q\left( x\right)
\right)
\end{equation*}%
on the $x$-axis would mandate%
\begin{equation}
\underset{x\rightarrow \pm a}{\lim }\psi \left( x\right) =0.
\end{equation}%
Feasibly, the satisfaction of this condition is obviously achievable when $%
\lambda >2$. This would imply that%
\begin{equation}
\lambda =\frac{1}{2}\left( 1+\sqrt{1+80g_{1}-64g_{2}}\right) >2.
\end{equation}%
At this point, we choose to stick with the continuity conditions at the
heterojunction boundaries between two crystals and take $\alpha =\gamma
\Rightarrow \beta =-1-2\alpha $ to imply%
\begin{equation}
\lambda =\frac{1}{2}\left( 1+\sqrt{\left( 3+8\alpha \right) ^{2}}\right)
>2\implies \left\vert 3+8\alpha \right\vert >3.
\end{equation}%
This suggests, in addition to von Roos constraint (2), new constraints $%
\alpha =\gamma <-3/4$ or $\alpha =\gamma >0$ on the ambiguity parameters. In
table 1, we list the parametric values associated with the available
orderings in the literature, along with their classifications as to being
admissible (i.e., $\lambda >2$ ) or non-admissible (i.e., $\lambda \leq 2$
). In the same table, we list some new orderings that are feasibly
admissible. \ Obviously, non of the orderings known in the literature can be
labeled as admissible, i.e., $\lambda >2$, (within our current methodical
proposal, of course).

\section{Concluding Remarks}

We have considered the PDM-von Roos Hamiltonian for a quantum particle
endowed with an unusual position-dependent mass function with two
singularities in an infinite potential well. A
point-canonical-transformation is used and an obvious correspondence between
two effective (reference and target) Hamiltonians is obtained. Apparently,
the singular PDM settings led to (through a PCT) the reflectionless modified
P\"{o}schl-Teller potential well. The exact solution of which is known in
different perspectives. We have followed Diaz et al. [25] to come out with
the result that non of the available known orderings in the literature is
labeled admissible. Hereby, it should be noted that, our classifications on
the "admissibility" and "non-admissibility" of the ordering ambiguity
parameters (reported in table 1) are consequences of the boundary condition
in (23). The validity of which resides only within the setting of the
current methodical proposal associated with our model in (4).

Nevertheless, one could re-scale $k^{2}$ of equation (2) of Lekner [21] for $%
\lambda =\nu +1$, for example, and deduce all related reflectionless
positive energy states, non-reflecting wave packets, etc. Yet, by re-scaling
$k^{2}$ of equation (39.2) of Fl\"{u}gge [37] for $\lambda =\nu +1$, one
would recycle Fl\"{u}gge's results. Consequently, we may conclude that the
ordering ambiguity conflict associated with the uniqueness of the kinetic
energy operator does not only depend on the heterojunction boundaries
between two crystals and the Dutra's and Almeida's [14] reliability test.
Even the form of the potential and/or the form of the position dependent
mass may have their say in the process, so to speak. Such ambiguity in the
non-uniqueness representation of the position-dependent mass Schr\"{o}dinger
Hamiltonian should be attributed to the lack of the Galilean invariance
(cf., e.g., ref.[33] for more details on the issue).

To the best of our knowledge, such a study has not been considered
elsewhere, not only within the Hermitian PDM-Hamiltonians' settings, but
neither within the complex non-Hermitian settings. A gap that remains
interesting and merits exploration/attention.

.\newpage

redpajama_set_name: RedPajamaArXiv