DARBOUX TRANSFORMATIONS FOR DUNKL-SCHRÖDINGER EQUATIONS WITH ENERGY-DEPENDENT
POTENTIAL AND POSITION-DEPENDENT MASS
Axel Schulze-Halberg† and Pinaki Roy‡,∗

$\dagger~{}$Department of Mathematics and Actuarial Science and Department of
Physics, Indiana University Northwest, 3400 Broadway, Gary IN 46408, USA,
E-mail<EMAIL_ADDRESS>

$\ddagger~{}$Atomic Molecular and Optical Physics Research Group, Advanced
Institute of Materials Science, Ton Duc Thang University, Ho Chi Minh City,
Vietnam

$\ast~{}$Faculty of Applied Sciences, Ton Duc Thang University, Ho Chi Minh
City, Vietnam,
E-mail<EMAIL_ADDRESS>
###### Abstract

We construct arbitrary-order Darboux transformations for Schrödinger equations
with energy-dependent potential and position-dependent mass within the Dunkl
formalism. Our construction is based on a point transformation that
interrelates our equations with the standard Schrödinger case. We apply our
method to generate several solvable Dunkl-Schrödinger equations.

Keywords: Dunkl operator, Darboux transformation, Schrödinger equation,
position-dependent mass, energy-dependent potential

##  1 Introduction

Dunkl operators [9] are commuting differential-difference operators associated
with finite reflection groups. As generalizations of partial derivatives,
these operators appear in a wide range of mathematical applications, such as
Fourier analysis related to root systems [2], intertwining operator and
angular momentum algebra [12] [11], integrable Calogero-Moser-Sutherland (CMS)
models [24] [10], harmonic analysis [5], the generation of orthogonal
polynomials [17], and nonlinear wave equations [18]. Another field of
application arises in quantum physics when conventional derivatives are
replaced by Dunkl operators. This replacement can be done in the momentum
operators, such that the associated Hamiltonian generates governing equations
that contain modified derivatives and reflection operators. In the simplest
case of a one-dimensional nonrelativistic system [6], the reflection operator
entering in the Schrödinger equation coincides with the parity operator. If
this equation admits solutions that have parity, then for these solutions the
parity operator can be replaced by a numerical parameter, which in some cases
allows for closed-form construction of the latter solutions. Based on this
concept, results have been obtained for a variety of systems, for example the
three-dimensional Coulomb problem [15], three-dimensional relativistic and
nonrelativistic systems with radial symmetry [19] [20], and the isotropic two-
dimensional Dunkl oscillator in the plane [13] [14], just to mention a few.
While the latter applications encompass several types of governing equations,
in the present work we will focus on the one-dimensional Schrödinger equation
that is equipped with a position-dependent mass and an energy-dependent
potential. Implementation of these concepts within the Dunkl formalism lead to
a straightforward modification of the weighted Hilbert space associated with
the Hamiltonian [21], being very similar to the standard formulation without
Dunkl operators [23]. In order to find cases of our Dunkl-Schrödinger equation
that admits closed-form solutions, we adapt point transformations and Darboux
transformations to the present scenario. Both of these techniques have been
widely applied to nonrelativistic quantum systems with position-dependent
mass, for example in the study of finite gap systems [4], the problem of
operator ordering in Hamiltonians [16], and the construction of coherent
states [8]. In contrast to these and related applications, to the best of our
knowledge point transformations and Darboux transformations have not been
applied yet to Dunkl-Schrödinger equations equipped with both a position-
dependent mass and an energy-dependent potential. After collecting some basic
facts about the Darboux transformation in section 2, we outline basic quantum
theory in section 3 as it applies to systems with position-dependent mass and
energy-dependent potentials. Section 4 is devoted to point transformations and
their application to a particular system that admits bound states. In section
5 we develop Darboux transformations for the systems under consideration here,
and we apply both the standard and the confluent Darboux algorithm to generate
solvable Dunkl-Schrödinger systems that support bound states.

##  2 Preliminaries: Darboux transformations

In order to make this work self-contained, we will now present a brief summary
of basic facts about Darboux transformations. Our starting point is the pair
of Schrödinger equations

$\displaystyle\Phi^{\prime\prime}(y)+\left[\epsilon-U(y)\right]\Phi(y)$
$\displaystyle=$ $\displaystyle 0$ (1)
$\displaystyle\hat{\Phi}^{\prime\prime}(y)+\left[\epsilon-\hat{U}(y)\right]\hat{\Phi}(y)$
$\displaystyle=$ $\displaystyle 0,$ (2)

where $\epsilon$ denotes the real-valued stationary energy, the functions $U$,
$\hat{U}$ stand for the potentials, and $\Phi$, $\hat{\Phi}$ are solutions of
their respective equation. In order to establish an interrelation between the
two equations (1) and (2) by means of a Darboux transformation, we distinguish
between the standard algorithm and the confluent algorithm.

  * •

Standard algorithm: We choose solutions (called transformation functions)
$u_{1},u_{2},...,u_{n}$ of our initial equation (1), that pertain to the
stationary energies (called transformation energies)
$\epsilon_{1},\epsilon_{2},...,\epsilon_{n}$, respectively. This means that we
have

$\displaystyle u_{j}^{\prime\prime}(y)+\left[\epsilon_{j}-U(y)\right]u_{j}(y)$
$\displaystyle=$ $\displaystyle 0,~{}~{}~{}j=1,2,...,n.$

Then, the function

$\displaystyle\hat{\Phi}(y)$ $\displaystyle=$
$\displaystyle\frac{W_{u_{1},u_{2},...,u_{n},\Phi}(y)}{W_{u_{1},u_{2},...,u_{n}}(y)},$
(3)

is a solution of equation (2), provided the transformed potential $\hat{U}$
satisfies the constraint

$\displaystyle\hat{U}(y)$ $\displaystyle=$ $\displaystyle
U(y)-2~{}\frac{d^{2}}{dy^{2}}~{}\log\left[W_{u_{1},u_{2},...,u_{n}}(y)\right].$
(4)

We refer to (3) as Darboux transformation of order $n$.

  * •

Confluent algorithm: In contrast to the previous case, we have a single
transformation energy $\epsilon_{1}$ only, and the transformation functions
obey the following system

$\displaystyle u_{1}^{\prime\prime}(y)+\left[\epsilon_{1}-U(y)\right]u_{1}(y)$
$\displaystyle=$ $\displaystyle 0$ (5) $\displaystyle
u_{j}^{\prime\prime}(y)+\left[\epsilon_{1}-U(y)\right]~{}u_{j}(y)$
$\displaystyle=$ $\displaystyle-u_{j-1}(y),~{}~{}~{}j=2,...,n.$ (6)

The expressions for the transformed solution (3) and potential (4) remain the
same as before.

##  3 The Dunkl-Schrödinger quantum system

Our goal for this section is to set up a Hamiltonian for position-dependent
mass within the Dunkl formalism. In addition, we will define a domain for our
Hamiltonian, such that it renders hermitian within a suitable Hilbert space.
To this end, let us first introduce the simplest case of a Dunkl operator [9]
in the form

$\displaystyle D_{x}$ $\displaystyle=$
$\displaystyle\frac{d}{dx}+\frac{\nu}{x}-\frac{\nu}{x}~{}{\cal R},$ (7)

where $\nu$ is a real number, and ${\cal R}$ stands for the reflection or
parity operator. While in the scenario of constant mass we have $\nu>-1/2$, in
the presence of a position-dependent mass admissible values for the parameter
$\nu$ are determined through the Hilbert space that we will define below.
Next, we use our operator (7) to build the Hamiltonian in the form

$\displaystyle H$ $\displaystyle=$ $\displaystyle-
D_{x}~{}\frac{1}{2~{}m(x)}~{}D_{x}+V_{E}(x),$ (8)

introducing a position-dependent mass $m$, and an energy-dependent potential
$V_{E}$, both of which are assumed to be sufficiently smooth functions that
are defined on the whole real line except possibly at the origin. Let us
simplify subsequent calculations by requiring the mass to have parity. More
precisely, we set

$\displaystyle{\cal R}~{}m(x)$ $\displaystyle=$ $\displaystyle\mu~{}m(x),$ (9)

such that an odd mass is represented by $\mu=-1$, while an even mass
corresponds to $\mu=1$. In the same way we describe the action of the
reflection operator on a function $\Psi$ that the Hamiltonian (8) acts on. We
set

$\displaystyle{\cal R}~{}\Psi(x)$ $\displaystyle=$
$\displaystyle\delta~{}\Psi(x),$ (10)

where in the odd case we have $\delta=-1$ and in the even case $\delta=1$
holds. Now, the Hamiltonian (8) is hermitian in the weighted Hilbert space
$L_{w}^{2}(\mathbb{R})$, where the weight function $w$ is given by

$\displaystyle w(x)$ $\displaystyle=$
$\displaystyle|x|^{2\nu-\delta\nu+\frac{\delta\nu}{\mu}}.$ (11)

The Dunkl-Schrödinger equation is obtained from our Hamiltonian (8) in the
form

$\displaystyle
D_{x}~{}\frac{1}{2~{}m(x)}~{}D_{x}~{}\Psi(x)+\left[E-V_{E}(x)\right]\Psi(x)$
$\displaystyle=$ $\displaystyle 0.$ (12)

where $\Psi$ represents a solution to the equation. Upon taking into account
our weight function (11) and the energy dependence of our potential, the
modified probability density $P$ for a solution $\Psi$ of (12) in
$L_{w}^{2}(\mathbb{R})$ is given by

$\displaystyle P_{\Psi}(x)$ $\displaystyle=$
$\displaystyle\Psi(x)^{\ast}\Psi(x)~{}|x|^{2\nu-\delta\nu+\frac{\delta\nu}{\mu}}\left[1-\frac{\partial
V_{E}(x)}{\partial E}\right].$ (13)

This defines the associated modified norm as

$\displaystyle N(\Psi)$ $\displaystyle=$
$\displaystyle\int\limits_{\mathbb{R}}P_{\Psi}(x)~{}dx~{}=~{}\int\limits_{\mathbb{R}}\Psi(x)^{\ast}\Psi(x)~{}|x|^{2\nu-\delta\nu+\frac{\delta\nu}{\mu}}\left[1-\frac{\partial
V_{E}(x)}{\partial E}\right]dx.$ (14)

If the potential $V_{E}$ does not depend on the energy, the term in square
brackets equals one. In this case we must require

$\displaystyle 2\nu-\delta\nu+\frac{\delta\nu}{\mu}$ $\displaystyle>$
$\displaystyle-1,$ (15)

such that the term $\sim|x|$ does not contribute singularities to the
integrand. For constant mass we have $\mu=1$, implying the known restriction
$\nu>-1/2$. If the potential $V_{E}$ is energy-dependent, then the restriction
(15) can be dropped, since the term in square brackets can contribute
singularities or remove them, depending on the behavior of the potential
$V_{E}$. In such a situation we must find restrictions for the parameter $\nu$
on a case-by-case basis. This is illustrated in our application section 4.2.
Let us now write (12) in expanded form by substituting (7) and by
incorporation of (10). This gives

$\displaystyle\frac{1}{2~{}m(x)}~{}\Psi^{\prime\prime}(x)+\Bigg{[}-\frac{m^{\prime}(x)}{2~{}m(x)^{2}}+\frac{\nu}{m(x)~{}x}-\frac{\nu~{}\delta}{2~{}m(x)~{}x}+\frac{\nu~{}\delta}{2~{}\mu~{}m(x)~{}x}\Bigg{]}~{}\Psi^{\prime}(x)+$
(16) $\displaystyle+$
$\displaystyle\Bigg{[}-\frac{\nu}{2~{}m(x)~{}x^{2}}-\frac{\nu~{}m^{\prime}(x)}{2~{}m(x)^{2}~{}x}+\frac{\nu~{}\delta}{2~{}m(x)~{}x^{2}}+\frac{\nu~{}\delta~{}m^{\prime}(x)}{2~{}m(x)^{2}~{}x}+\frac{\nu^{2}}{2~{}m(x)~{}x^{2}}-$
$\displaystyle-$
$\displaystyle\frac{\nu^{2}~{}\delta}{2~{}m(x)~{}x^{2}}+\frac{\nu^{2}~{}\delta}{2~{}\mu~{}m(x)~{}x^{2}}-\frac{\nu^{2}}{2~{}\mu~{}m(x)~{}x^{2}}+E-V_{E}(x)\Bigg{]}~{}\Psi(x)~{}=~{}0.$

If we distinguish the parameter values for $\delta$ and $\mu$, we obtain four
special cases of this equation. Let us show the two cases that arise from odd
and even position-dependent mass functions. We will study the two remaining
cases in detail below. For an odd-parity mass ($\mu=-1$) our equation takes
the form

$\displaystyle\frac{1}{2~{}m(x)}~{}\Psi^{\prime\prime}(x)+\Bigg{[}-\frac{m^{\prime}(x)}{2~{}m(x)^{2}}+\frac{\nu}{m(x)~{}x}-\frac{\nu~{}\delta}{m(x)~{}x}\Bigg{]}\Psi^{\prime}(x)+\Bigg{[}-\frac{\nu}{2~{}m(x)~{}x^{2}}-\frac{\nu~{}m^{\prime}(x)}{2~{}m(x)^{2}~{}x}+$
$\displaystyle+~{}\frac{\nu~{}\delta}{2~{}m(x)~{}x^{2}}+\frac{\nu~{}\delta~{}m^{\prime}(x)}{2~{}m(x)^{2}~{}x}+\frac{\nu^{2}}{m(x)~{}x^{2}}-\frac{\nu^{2}~{}\delta}{m(x)~{}x^{2}}+E-V_{E}(x)\Bigg{]}\Psi(x)~{}=~{}0.$
(17)

If our mass function has even parity ($\mu=1)$, then the Dunkl-Schrödinger
equation (16) reads

$\displaystyle\frac{1}{2~{}m(x)}~{}\Psi^{\prime\prime}(x)+\Bigg{[}-\frac{m^{\prime}(x)}{2~{}m(x)^{2}}+\frac{\nu}{m(x)~{}x}\Bigg{]}\Psi^{\prime}(x)+\Bigg{[}-\frac{\nu}{2~{}m(x)~{}x^{2}}-\frac{\nu~{}m^{\prime}(x)}{2~{}m(x)^{2}~{}x}+\frac{\nu~{}\delta}{2~{}m(x)~{}x^{2}}+$
$\displaystyle+~{}\frac{\nu~{}\delta~{}m^{\prime}(x)}{2~{}m(x)^{2}~{}x}+E-V_{E}(x)\Bigg{]}\Psi(x)~{}=~{}0.$
(18)

A particular case of this equation is the scenario of constant mass. We obtain
for $m=1/2$:

$\displaystyle\Psi^{\prime\prime}(x)+\frac{2~{}\nu}{x}~{}\Psi^{\prime}(x)+\Bigg{[}\frac{\nu~{}\delta-\nu}{x^{2}}+E-V_{E}(x)\Bigg{]}\Psi(x)~{}=~{}0.$
(19)

Note that the mass being constant implies $\mu=1$.

##  4 Point transformations

The first method for solving the Dunkl-Schrödinger equation (16) consists in a
point transformation that takes it into conventional Schrödinger form by
gauging away the first-derivative term $\sim\Psi^{\prime}$. Appropriate
choices for the free parameter functions then allow for the construction of a
solvable case.

###  4.1 Transformation to standard form

We apply the following point transformation

$\displaystyle\Phi(y)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{1}{m[x(y)]~{}x^{\prime}(y)}}~{}x(y)^{\nu-\frac{\delta\nu}{2}+\frac{\delta\nu}{2\mu}}~{}\Psi[x(y)].$
(20)

The coordinate change $x=x(y)$ remains arbitrary for now, except that it must
have an inverse. It is sufficient if the inverse exists on either the negative
or the positive axis, since we are considering solutions with parity only. The
specfic form of (20) ensures that the coefficient of the term
$\sim\Psi^{\prime}$ vanishes. Substitution of (20) renders (16) in the form

$\displaystyle\Phi^{\prime\prime}(y)$ $\displaystyle+$
$\displaystyle\Bigg{\\{}2~{}E~{}m[x(y)]~{}[x^{\prime}(y)]^{2}-2~{}m[x(y)]~{}V_{E}[x(y)]~{}[x^{\prime}(y)]^{2}+\frac{\delta~{}\nu~{}[x^{\prime}(y)]^{2}}{2~{}x(y)^{2}}+\frac{\delta~{}\nu~{}[x^{\prime}(y)]^{2}}{2~{}\mu~{}x(y)^{2}}-$
(21) $\displaystyle-$
$\displaystyle\frac{\nu^{2}~{}[x^{\prime}(y)]^{2}}{2~{}x(y)^{2}}-\frac{\nu^{2}~{}[x^{\prime}(y)]^{2}}{2~{}\mu~{}x(y)^{2}}+\frac{\delta~{}\nu~{}m^{\prime}[x(y)]~{}[x^{\prime}(y)]^{2}}{2~{}m[x(y)]~{}x(y)}+\frac{\delta~{}\nu~{}m^{\prime}[x(y)]~{}[x^{\prime}(y)]^{2}}{2~{}\mu~{}m[x(y)]~{}x(y)}-$
$\displaystyle-$
$\displaystyle\frac{3~{}\\{m^{\prime}[x(y)]\\}^{2}~{}[x^{\prime}(y)]^{2}}{4~{}m[x(y)]}+\frac{[x^{\prime}(y)]^{2}~{}m^{\prime\prime}[x(y)]}{2~{}m[x(y)]}-\frac{3~{}[x^{\prime\prime}(y)]^{2}}{4~{}[x^{\prime}(y)]^{2}}+\frac{x^{\prime\prime\prime}(y)}{2~{}x^{\prime}(y)}\Bigg{\\}}~{}\Phi(y)~{}=~{}0.$

Let us now recall that the conventional Schrödinger equation with energy-
dependent potential $U_{E}$ reads

$\displaystyle\Phi^{\prime\prime}(y)+\left[E-U_{E}(y)\right]\Phi(y)$
$\displaystyle=$ $\displaystyle 0.$ (22)

Hence, for the transformed Dunkl-Schrödinger equation (21) to take the form
(22), the following condition must be imposed on the potential $U_{E}$:

$\displaystyle U_{E}(y)$ $\displaystyle=$ $\displaystyle
E-2~{}E~{}m[x(y)]~{}[x^{\prime}(y)]^{2}+2~{}m[x(y)]~{}V_{E}[x(y)]~{}[x^{\prime}(y)]^{2}-\frac{\delta~{}\nu~{}[x^{\prime}(y)]^{2}}{2~{}x(y)^{2}}-$
(23)
$\displaystyle-~{}\frac{\delta~{}\nu~{}[x^{\prime}(y)]^{2}}{2~{}\mu~{}x(y)^{2}}+\frac{\nu^{2}~{}[x^{\prime}(y)]^{2}}{2~{}x(y)^{2}}+\frac{\nu^{2}~{}[x^{\prime}(y)]^{2}}{2~{}\mu~{}x(y)^{2}}-\frac{\delta~{}\nu~{}m^{\prime}[x(y)]~{}[x^{\prime}(y)]^{2}}{2~{}m[x(y)]~{}x(y)}-\frac{\delta~{}\nu~{}m^{\prime}[x(y)]~{}[x^{\prime}(y)]^{2}}{2~{}\mu~{}m[x(y)]~{}x(y)}+$
$\displaystyle+~{}\frac{3~{}\\{m^{\prime}[x(y)]\\}^{2}~{}[x^{\prime}(y)]^{2}}{4~{}m[x(y)]}-\frac{[x^{\prime}(y)]^{2}~{}m^{\prime\prime}[x(y)]}{2~{}m[x(y)]}+\frac{3~{}[x^{\prime\prime}(y)]^{2}}{4~{}[x^{\prime}(y)]^{2}}+\frac{x^{\prime\prime\prime}(y)}{2~{}x^{\prime}(y)}.$

Note that this condition arises simply by comparing the coefficients of $\Phi$
in the two equations (21) and (22). In the special case of constant mass
$m=1/2$ and $\mu=1$, the condition (23) reads

$\displaystyle U_{E}(y)$ $\displaystyle=$ $\displaystyle
E-E~{}[x^{\prime}(y)]^{2}+V_{E}[x(y)]~{}[x^{\prime}(y)]^{2}-\frac{\delta~{}\nu~{}[x^{\prime}(y)]^{2}}{x(y)^{2}}+\frac{\nu^{2}~{}[x^{\prime}(y)]^{2}}{x(y)^{2}}+\frac{3~{}[x^{\prime\prime}(y)]^{2}}{4~{}[x^{\prime}(y)]^{2}}-\frac{x^{\prime\prime\prime}(y)}{2~{}x^{\prime}(y)}.$

Since our conditions (23) and (LABEL:con0) contain three and two free
parameter functions, respectively, the same function $U_{E}$ can be obtained
by different choices for the latter parameter functions. As such, each $U_{E}$
gives rise to a class of Dunkl-Schrödinger systems with energy-dependent
potential and position-dependent mass. We will comment further on this
property of (23) and (LABEL:con0) in our application section 5.3. Let us now
present examples for the construction of a Dunkl-Schrödinger system with
energy-dependent potential by means of point transformations.

###  4.2 Application

Our starting point is the conventional Schrödinger equation (1). If we can
generate a potential $U_{E}$ that satisfies the condition (23), such that (1)
is solvable, then our point transformation (20) allows for the construction of
an associated Dunkl-Schrödinger equation that is also solvable. Let us make
the following choices for the position-dependent mass and the potential in the
Dunkl-Schrödinger system that we want to construct:

$\displaystyle m(x)$ $\displaystyle=$ $\displaystyle
p~{}\exp\left(-q~{}x^{2}\right)$ (25) $\displaystyle V_{E}(x)$
$\displaystyle=$ $\displaystyle
E-2~{}p~{}E~{}\exp\left(q~{}x^{2}\right)-\frac{p}{2}~{}\exp\left(q~{}x^{2}\right)$
(26) $\displaystyle x(y)$ $\displaystyle=$ $\displaystyle\sqrt{y},$ (27)

where $p$ and $q$ stand for arbitrary real numbers. We will use particular
cases of our specific mass profile (25) in the subsequent sections. Note that
our change of coordinate is defined and invertible on the positive half-axis,
which is sufficient for our purposes. Furthermore, we chose the particular
change of coordinate (27), such that the conventional Schrödinger equation (1)
solvable. Substitution of (27) into the potential condition (23) yields

$\displaystyle U_{E}(y)$ $\displaystyle=$ $\displaystyle
E+\frac{q^{2}}{4}-\left(-\frac{q}{4}+\frac{p^{2}}{4}+p^{2}~{}E-\frac{q~{}\delta~{}\nu}{2}\right)\frac{1}{y}-\left(\frac{3}{16}+\frac{\delta~{}\nu}{4}-\frac{\nu^{2}}{4}\right)\frac{1}{y^{2}}.$

In order to keep our calculations simple and transparent, let us take a
special case of the mass function (25) by applying the parameter setting
$p=q=1$. Then (LABEL:uegen) becomes

$\displaystyle U_{E}(y)$ $\displaystyle=$ $\displaystyle
E+\frac{1}{4}-\left(E-\frac{\delta~{}\nu}{2}\right)\frac{1}{y}-\left(\frac{3}{16}+\frac{\delta~{}\nu}{4}-\frac{\nu^{2}}{4}\right)\frac{1}{y^{2}}.$

The corresponding equation (1) then takes the form

$\displaystyle\Phi^{\prime\prime}(y)+\left[-\frac{1}{4}+\left(E-\frac{\delta~{}\nu}{2}\right)\frac{1}{y}+\left(\frac{3}{16}+\frac{\delta~{}\nu}{4}-\frac{\nu^{2}}{4}\right)\frac{1}{y^{2}}\right]\Phi(y)$
$\displaystyle=$ $\displaystyle 0.$ (29)

The general solution to this equation can be given in terms of hypergeometric
functions. We have

$\displaystyle\Phi_{\tiny{\mbox{gen}}}(y)$ $\displaystyle=$
$\displaystyle\exp\left(-\frac{y}{2}\right)y^{\frac{1}{2}+\frac{1}{4}\sqrt{1-4\delta\nu+4\nu^{2}}}\times$
(30) $\displaystyle\times$
$\displaystyle\Bigg{[}c_{1}~{}{}_{1}F_{1}\left(\frac{1}{2}-E+\frac{\delta~{}\nu}{2}+\frac{1}{4}~{}\sqrt{1-4~{}\delta~{}\nu+4~{}\nu^{2}},1+\frac{1}{2}~{}\sqrt{1-4~{}\delta~{}\nu+4~{}\nu^{2}},y\right)+$
$\displaystyle+$ $\displaystyle
c_{2}~{}U\left(\frac{1}{2}-E+\frac{\delta~{}\nu}{2}+\frac{1}{4}~{}\sqrt{1-4~{}\delta~{}\nu+4~{}\nu^{2}},1+\frac{1}{2}~{}\sqrt{1-4~{}\delta~{}\nu+4~{}\nu^{2}},y\right)\Bigg{]},$

where $c_{1}$ and $c_{2}$ represent arbitrary constants. Furthermore,
${}_{1}F_{1}$ and $U$ stand for the confluent hypergeometric functions of
first and second kind, respectively [1]. Since we are aiming for the
construction of bound states, we must work with a solution that is bounded on
the whole positive half-axis. This condition is not satisfied by the function
$U$, such that we can set $c_{1}=1$ and $c_{2}=0$ in (30). As a result, we
obtain the particular solution

$\displaystyle\Phi(y)$ $\displaystyle=$
$\displaystyle\exp\left(-\frac{y}{2}\right)y^{\frac{1}{2}+\frac{1}{4}\sqrt{1-4\delta\nu+4\nu^{2}}}\times$
$\displaystyle\times$
$\displaystyle{}_{1}F_{1}\left(\frac{1}{2}-E+\frac{\delta~{}\nu}{2}+\frac{1}{4}~{}\sqrt{1-4~{}\delta~{}\nu+4~{}\nu^{2}},1+\frac{1}{2}~{}\sqrt{1-4~{}\delta~{}\nu+4~{}\nu^{2}},y\right).$

Upon reversing our point transformation (20), we obtain the function $\Psi$ in
the form

$\displaystyle\Psi(x)$ $\displaystyle=$ $\displaystyle
x^{\frac{1}{2}-\nu+\frac{1}{2}\sqrt{1-4\delta\nu+4\nu^{2}}}\times$ (31)
$\displaystyle\times$
$\displaystyle{}_{1}F_{1}\left(\frac{1}{2}+E-\frac{\delta~{}\nu}{2}+\frac{1}{4}~{}\sqrt{1-4~{}\delta~{}\nu+4~{}\nu^{2}},1+\frac{1}{2}~{}\sqrt{1-4~{}\delta~{}\nu+4~{}\nu^{2}},-x^{2}\right).$

It is convenient to rewrite the confluent hypergeometric function by means of
an identity [25] as follows

$\displaystyle\Psi(x)$ $\displaystyle=$
$\displaystyle\exp\left(-x^{2}\right)x^{\frac{1}{2}-\nu+\frac{1}{2}\sqrt{1-4\delta\nu+4\nu^{2}}}\times$
(32) $\displaystyle\times$
$\displaystyle{}_{1}F_{1}\left(\frac{1}{2}-E+\frac{\delta~{}\nu}{2}+\frac{1}{4}~{}\sqrt{1-4~{}\delta~{}\nu+4~{}\nu^{2}},1+\frac{1}{2}~{}\sqrt{1-4~{}\delta~{}\nu+4~{}\nu^{2}},x^{2}\right).$

Since the point transformation (20) interrelates the Dunkl-Schrödinger
equation (16) and its conventional counterpart (1), the function (32) must
solve (16). Upon plugging the settings (25) and (26) for $p=q=1$ into (16), we
obtain the latter equation as

$\displaystyle\frac{1}{2}~{}\Psi^{\prime\prime}(x)+\left(x+\frac{\nu}{x}\right)\Psi^{\prime}(x)+\left(\frac{1}{2}+2~{}E+\nu-\frac{\nu}{2~{}x^{2}}-\delta~{}\nu+\frac{\delta~{}\nu}{2~{}x^{2}}\right)\Psi(x)~{}=~{}0.$

Note that even if (32) is a solution to the latter equation, it does not
automatically solve the initial Dunkl-Schrödinger equation (12), unless it has
the correct parity behavior given by (10). In order to investigate this
behavior, we note that both the exponential and the confluent hypergeometric
function are even, such that the overall behavior depends on the exponent of
the monomial term. Let us distinguish the two cases of $\delta=-1$ and
$\delta=1$. In the first of these cases the exponent reads

$\displaystyle\frac{1}{2}-\nu+\frac{1}{2}~{}\sqrt{1+4~{}\nu+4~{}\nu^{2}}$
$\displaystyle=$ $\displaystyle\frac{1}{2}-\nu+\frac{1}{2}~{}|1+2~{}\nu|.$
(33)

At this point we cannot proceed until we know the admissible values for the
parameter $\nu$. According to section 3, due to the energy dependence of our
potential $V_{E}$, we must find these admissible values from the behavior of
the term $\partial V_{E}/\partial E$ that enters in (14). Substitution of (27)
yields

$\displaystyle 1-\frac{\partial V_{E}(x)}{\partial E}$ $\displaystyle=$
$\displaystyle 2~{}\exp\left(x^{2}\right),$ (34)

which gives the standard restriction $\nu>-1/2$ because the confluent
hypergeometric function is regular at the origin in the cases of bound states
that we are interested in. Upon using our restriction in (33), we obtain

$\displaystyle\frac{1}{2}-\nu+\frac{1}{2}~{}|1+2~{}\nu|~{}=~{}\frac{1}{2}-\nu+\frac{1}{2}~{}(1+2~{}\nu)~{}=~{}1.$

Since the exponent is an odd number, the function $(\ref{psi0})$ has odd
parity, and as such solves the initial Dunkl-Schrödinger equation. Next we
consider the case $\delta=1$, the exponent of the middle term on the right
side of (32) takes the form

$\displaystyle\frac{1}{2}-\nu+\frac{1}{2}~{}\sqrt{1-4~{}\nu+4~{}\nu^{2}}$
$\displaystyle=$ $\displaystyle\frac{1}{2}-\nu+\frac{1}{2}~{}|1-2~{}\nu|$ (37)
$\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{lll}1-2~{}\nu&\mbox{if}&\nu<\frac{1}{2}\\\\[4.30554pt]
0&\mbox{if}&\nu\geq\frac{1}{2}\end{array}\right\\}.$

Hence, if $\delta=1$ and $\nu\geq 1/2$, the function (32) has even parity,
such that it solves our initial Dunkl-Schrödinger equation. If $\nu<1/2$, then
the exponent depends on $\nu$, resulting in the condition

$\displaystyle 1-2~{}\nu$ $\displaystyle=$ $\displaystyle 2~{}k,$ (38)

for an arbitrary integer $k$. If we take into account that $\nu>-1/2$ due to
(34), there is no solution of (38). Hence, in the present case we have the
overall constraint

$\displaystyle\nu$ $\displaystyle\geq$ $\displaystyle\frac{1}{2}.$

Now, since we want to construct solutions of bound state type, we restrict the
stationary energy, such that the confluent hypergeometric function in (32)
degenerates to a polynomial. This restriction is given by equating the first
argument of the latter function to a nonpositive integer, leading to

$\displaystyle E$ $\displaystyle=$ $\displaystyle
n+\frac{1+\delta~{}\nu}{2}+\frac{1}{4}~{}\sqrt{1-4~{}\delta~{}\nu+4~{}\nu^{2}},$
(39)

where $n$ is a nonnegative integer. Let us state examples of our solutions
(32) for the parameter setting $\nu=1/2$, $\delta=-1$ (odd parity), and the
parameter $n$ as given in (39). We find

$\displaystyle\Psi(x)_{n=0}$ $\displaystyle=$
$\displaystyle\exp\left(-x^{2}\right)x$ (40) $\displaystyle\Psi(x)_{n=1}$
$\displaystyle=$
$\displaystyle\exp\left(-x^{2}\right)x\left(1-\frac{x^{2}}{2}\right)$ (41)
$\displaystyle\Psi(x)_{n=2}$ $\displaystyle=$
$\displaystyle\exp\left(-x^{2}\right)x\left(1-x^{2}+\frac{x^{4}}{2}\right).$
(42)

The graphs of these solutions are shown in figure 1. Now we state the
corresponding solutions for the case of even parity ($\delta=1$), leaving the
remaining parameters unchanged. This gives

$\displaystyle\Psi(x)_{n=0}$ $\displaystyle=$
$\displaystyle\exp\left(-x^{2}\right)$ (43) $\displaystyle\Psi(x)_{n=1}$
$\displaystyle=$ $\displaystyle\exp\left(-x^{2}\right)\left(1-x^{2}\right)$
(44) $\displaystyle\Psi(x)_{n=2}$ $\displaystyle=$
$\displaystyle\exp\left(-x^{2}\right)\left(1-2~{}x^{2}+\frac{x^{4}}{2}\right).$
(45)

Graphs of these solutions are displayed in figure 2. Next, we verify
normalizability of our solution (32) with respect to the modified norm (14).
Taking into account that (34) is a positive function, and implementing the
present settings (27) and $\mu=1$, the probability density and the associated
norm read

$\displaystyle P_{\Psi}(x)$ $\displaystyle=$
$\displaystyle\Psi(x)^{\ast}\Psi(x)~{}|x|^{2\nu}\exp\left(x^{2}\right)$
$\displaystyle N(\Psi)$ $\displaystyle=$ $\displaystyle
2~{}\int\limits_{\mathbb{R}}P_{\Psi}(x)~{}dx,$ (46)

observe that $\delta$ does not enter in the norm due to the mass having even
parity. The integral (46) exists for (32) under the energy restriction (39),
such that the latter solution is normalizable, and thus represents bound
states. Figure 3 shows normalized probability densities associated with the
solutions (40)-(42) and (43)-(45), respectively.

Figure 1: Graphs of the odd solutions (40)-(42). Figure 2: Graphs of the even
solutions (43)-(45). Figure 3: Graphs of normalized probability densities
pertaining to the odd solutions (40)-(42).

##  5 Darboux transformations

Besides point transformations that we studied in the previous section, a
further method for generating solvable equations of Schrödinger form is given
by Darboux transformations. These cannot be applied to the Dunkl-Schrödinger
equation (12) directly due to its form, such that we first need to convert it
suitably. More precisely, we will use the point transformation (20) for this
task. In addition, we must distinguish between the cases of odd-parity mass
and even-parity mass.

###  5.1 Odd-parity mass functions

Before we start the construction of our Darboux transformation, let us point
out that odd-parity mass functions must be negative either on the positive or
on the negative axis. Even though there are some applications for negative
mass functions, to the best of our knowledge they are generally considered
non-physical. For this reason we will develop the Darboux transformation, but
omit to state an example for it. Now, if the mass has odd parity, we have
$\mu=-1$ according to (9), such that the Dunkl-Schrödinger equation takes the
form (17) after introduction of the parameter $\delta$ in (10). Next, we apply
our point transformation (20), such that the resulting equation is given by
(21) for $\mu=-1$. In explicit form we have

$\displaystyle\Phi^{\prime\prime}(y)$ $\displaystyle+$
$\displaystyle\Bigg{\\{}2~{}E~{}m[x(y)]~{}[x^{\prime}(y)]^{2}-2~{}m[x(y)]~{}V_{E}[x(y)]~{}[x^{\prime}(y)]^{2}-\frac{3~{}\\{m^{\prime}[x(y)]\\}^{2}~{}[x^{\prime}(y)]^{2}}{4~{}m[x(y)]}+$
(47) $\displaystyle+$
$\displaystyle\frac{[x^{\prime}(y)]^{2}~{}m^{\prime\prime}[x(y)]}{2~{}m[x(y)]}-\frac{3~{}[x^{\prime\prime}(y)]^{2}}{4~{}[x^{\prime}(y)]^{2}}+\frac{x^{\prime\prime\prime}(y)}{2~{}x^{\prime}(y)}\Bigg{\\}}~{}\Phi(y)~{}=~{}0.$

Inspection of this equation shows that the parameters $\nu$ and $\delta$ from
the Dunkl operator (7) are not present anymore. In other words, the Dunkl
scenario enters in our system entirely by means of the point transformation
(20), provided the mass function has odd parity. Now, in order to apply a
Darboux transformation to (47), we must determine a transformation parameter
that plays the role of the stationary energy in the conventional Schrödinger
case. We choose this parameter to be $E$, which implies that its coefficient
must be equal to one. The condition reads

$\displaystyle 2~{}m[x(y)]~{}[x^{\prime}(y)]^{2}$ $\displaystyle=$
$\displaystyle 1.$ (48)

Implementation of this constraint in (47) renders our equation in the form

$\displaystyle\Phi^{\prime\prime}(y)+\left\\{E-V_{E}[x(y)]-\frac{\\{m^{\prime}[x(y)]\\}^{2}}{8~{}m[x(y)]^{3}}\right\\}\Phi(y)$
$\displaystyle=$ $\displaystyle 0.$ (49)

We can now apply a Darboux transformation to (49). To this end, we use our
formulas (3) and (4), where the potential $U$ must be replaced by $V_{E}$. We
then have for the transformed solution and the associated potential

$\displaystyle\hat{\Phi}(y)~{}=~{}\frac{W_{u_{1},...,u_{n},\Phi}(y)}{W_{u_{1},...,u_{n}}(y)}\qquad\qquad\qquad\hat{V}_{E}[x(y)]~{}=~{}V_{E}[x(y)]-2~{}\frac{d^{2}}{dy^{2}}~{}\log\left[W_{u_{1},...,u_{n}}(y)\right],$

which enter in the transformed Schrödinger equation (2) that is the partner of
(49), and reads in the present case

$\displaystyle\hat{\Phi}^{\prime\prime}(y)+\left\\{E-\hat{V}_{E}[x(y)]-\frac{\\{m^{\prime}[x(y)]\\}^{2}}{8~{}m[x(y)]^{3}}\right\\}\hat{\Phi}(y)$
$\displaystyle=$ $\displaystyle 0.$

The next step consists in reverting the point transformation (20) for
reinstating the Dunkl-Schrödinger form. This reversion is achieved by setting

$\displaystyle\hat{\Psi}(x)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{m(x)}{y^{\prime}(x)}}~{}x^{-\nu+\delta\nu}~{}\hat{\Phi}[y(x)],$

where $y=y(x)$ is the inverse function of the initial coordinate change
$x=x(y)$. Now, the function $\hat{\Psi}$ solves the transformed counterpart of
(17) that is given by

$\displaystyle\frac{1}{2~{}m(x)}~{}\hat{\Psi}^{\prime\prime}(x)+\Bigg{[}-\frac{m^{\prime}(x)}{2~{}m(x)^{2}}+\frac{\nu}{m(x)~{}x}-\frac{\nu~{}\delta}{m(x)~{}x}\Bigg{]}\hat{\Psi}^{\prime}(x)+\Bigg{[}-\frac{\nu}{2~{}m(x)~{}x^{2}}-\frac{\nu~{}m^{\prime}(x)}{2~{}m(x)^{2}~{}x}+$
$\displaystyle+~{}\frac{\nu~{}\delta}{2~{}m(x)~{}x^{2}}+\frac{\nu~{}\delta~{}m^{\prime}(x)}{2~{}m(x)^{2}~{}x}+\frac{\nu^{2}}{m(x)~{}x^{2}}-\frac{\nu^{2}~{}\delta}{m(x)~{}x^{2}}+E-\hat{V}_{E}(x)\Bigg{]}\hat{\Psi}(x)~{}=~{}0.$
(50)

Let us point out that any solution $\hat{\Psi}$ of this equation only solves
the transformed Dunkl-Schrödinger equation

$\displaystyle
D_{x}~{}\frac{1}{2~{}m(x)}~{}D_{x}~{}\hat{\Psi}(x)+\left[E-\hat{V}_{E}(x)\right]\hat{\Psi}(x)$
$\displaystyle=$ $\displaystyle 0,$ (51)

if $\hat{\Psi}$ has the correct parity given by the parameter $\delta$ in
(50). In summary, we have established Darboux transformations between the
Dunkl-Schrödinger partner equations (17) and (50).

###  5.2 Even-parity mass functions

As in the previous case, our starting point is the Dunkl-Schrödinger equation
(12) that we consider in the form (18) for the parameter setting $\mu=1$.
Application of our point transformation (20) yields

$\displaystyle\Phi^{\prime\prime}(y)$ $\displaystyle+$
$\displaystyle\Bigg{\\{}2~{}E~{}m[x(y)]~{}[x^{\prime}(y)]^{2}-2~{}m[x(y)]~{}V_{E}[x(y)]~{}[x^{\prime}(y)]^{2}+$
$\displaystyle+$
$\displaystyle\frac{\delta~{}\nu~{}[x^{\prime}(y)]^{2}}{x(y)^{2}}-\frac{\nu^{2}~{}[x^{\prime}(y)]^{2}}{x(y)^{2}}+\frac{\delta~{}\nu~{}m^{\prime}[x(y)]~{}[x^{\prime}(y)]^{2}}{m[x(y)]~{}x(y)}-$
$\displaystyle-$
$\displaystyle\frac{3~{}\\{m^{\prime}[x(y)]\\}^{2}~{}[x^{\prime}(y)]^{2}}{4~{}m[x(y)]}+\frac{[x^{\prime}(y)]^{2}~{}m^{\prime\prime}[x(y)]}{2~{}m[x(y)]}-\frac{3~{}[x^{\prime\prime}(y)]^{2}}{4~{}[x^{\prime}(y)]^{2}}+\frac{x^{\prime\prime\prime}(y)}{2~{}x^{\prime}(y)}\Bigg{\\}}~{}\Phi(y)~{}=~{}0.$

In contrast to the previously considered case of odd-parity mass, this time
our equation contains both parameters $\nu$ and $\delta$ from the Dunkl
operator (7). While so far we have followed the same approach as for the odd-
parity mass function, at this point we need to proceed differently. The reason
lies in the presence of the Dunkl parameters $\nu$ and $\delta$ in our
equation (LABEL:dsse1even). More precisely, if we used the setting (48) and
took the stationary energy $E$ as our parameter for the Darboux
transformation, then the transformed potential would depend on $\nu$ and
$\delta$, which is not desirable. Therefore, we would have to assign numerical
values to these Dunkl parameters, which would constrain parity of the solution
to (LABEL:dsse1even). Furthermore, we would need to introduce new constants
that play the role of the Dunkl parameters, such that the Darboux transformed
counterpart of (LABEL:dsse1even) can be mapped back onto Dunkl-Schrödinger
form. While in general these issues must be solved on a case-by-case basis,
they can be avoided entirely by restricting the mass function and the
coordinate change as follows

$\displaystyle m(x)~{}=~{}p~{}x^{q}\qquad\qquad\qquad x(y)~{}=~{}\exp(y),$
(53)

where $p$ is an arbitrary constant, and $q$ is an even integer. Note that our
coordinate change is invertible on the whole real axis. Before we continue,
let us point out that the mass function in (53) for $q\neq 0$ is either
singular at the origin or vanishes there. While in most applications this is
not a desired behavior, here we introduce the mass function in (53) merely due
to a specific mathematical property. More precisely, the effect of the
settings (53) on our equation (LABEL:dsse1even) is that the Dunkl parameters
now appear as a constant term in the coefficient of the solution $\Phi$:

$\displaystyle\Phi^{\prime\prime}(y)$ $\displaystyle+$
$\displaystyle\Bigg{(}-\frac{(q+1)^{2}}{4}+(q+1)~{}\delta~{}\nu-\nu^{2}+2~{}p~{}\exp\left[(q+2)~{}y\right]\left\\{E-V_{E}\left[\exp(y)\right]\right\\}\Bigg{)}~{}\Phi(y)~{}=~{}0.$

We can rewrite this equation in the form

$\displaystyle\Phi^{\prime\prime}(y)+\left[\epsilon-
U_{E}(y)\right]~{}\Phi(y)~{}=~{}0,$ (55)

introducing the abbreviations $\epsilon$ and $U_{E}$ that are given by

$\displaystyle\epsilon$ $\displaystyle=$
$\displaystyle(q+1)~{}\delta~{}\nu-\nu^{2}$ (56) $\displaystyle U_{E}(y)$
$\displaystyle=$
$\displaystyle\frac{(q+1)^{2}}{4}-2~{}p~{}\exp\left[(q+2)~{}y\right]\left\\{E-V_{E}\left[\exp(y)\right]\right\\}.$
(57)

Therefore, we can now choose $\epsilon$ as the Darboux transformation
parameter instead of the stationary energy. As a result, the transformed
potential will not depend on the Dunkl parameters, and the transformed
equation can be converted back to Dunkl-Schrödinger form in a straightforward
manner. Application of the Darboux transformation (3) and (4) while replacing
$U$ by $U_{E}$, we obtain the transformed solution and the associated
potential as

$\displaystyle\hat{\Phi}(y)~{}=~{}\frac{W_{u_{1},...,u_{n},\Phi}(y)}{W_{u_{1},...,u_{n}}(y)}\qquad\qquad\qquad\hat{U}_{E}(y)~{}=~{}U_{E}(y)-2~{}\frac{d^{2}}{dy^{2}}~{}\log\left[W_{u_{1},...,u_{n}}(y)\right],$

These functions enter in the transformed counterpart of (55) that reads

$\displaystyle\hat{\Phi}^{\prime\prime}(y)+\left[\epsilon-\hat{U}_{E}(y)\right]~{}\hat{\Phi}(y)~{}=~{}0.$
(58)

Upon rewriting this equation in explicit form, we obtain the transformed
partner of (LABEL:sse0x) as

$\displaystyle\hat{\Phi}^{\prime\prime}(y)$ $\displaystyle+$
$\displaystyle\Bigg{(}-\frac{(q+1)^{2}}{4}+(q+1)~{}\delta~{}\nu-\nu^{2}+2~{}p~{}\exp\left[(q+2)~{}y\right]\left\\{E-\hat{V}_{E}\left[\exp(y)\right]\right\\}\Bigg{)}\hat{\Phi}(y)~{}=~{}0.$

Observe that the potential $\hat{V}_{E}$ can be obtained from the transformed
version of (57) that is given by

$\displaystyle\hat{U}_{E}(y)$ $\displaystyle=$
$\displaystyle\frac{(q+1)^{2}}{4}-2~{}p~{}\exp\left[(q+2)~{}y\right]\left\\{E-\hat{V}_{E}\left[\exp(y)\right]\right\\}.$

Solving with respect to $\hat{V}_{E}$ yields the result

$\displaystyle\hat{V}_{E}\left[\exp(y)\right]$ $\displaystyle=$ $\displaystyle
E-\exp\left[-(q+2)~{}y\right]\left[\frac{(q+1)^{2}}{8~{}p}-\frac{\hat{U}_{E}(y)}{2~{}p}\right].$
(60)

Now that the Darboux transformation is complete, the remaining task consists
in converting our equation (LABEL:sse2xt) back to Dunkl-Schrödinger form by
reversing the point transformation (20). Recall that in the present case we
made the settings (53), such that the reversed point transformation reads

$\displaystyle\hat{\Psi}(x)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{m(x)}{y^{\prime}(x)}}~{}x^{-\nu}~{}\hat{\Phi}[y(x)]~{}=~{}\sqrt{p}~{}x^{\frac{1}{2}+\frac{q}{2}-\nu}~{}\hat{\Phi}\left[\log(x)\right],$
(61)

note that the function $y=y(x)=\log(x)$ stands for the inverse of
$x=x(y)=\exp(y)$. Next, we rewrite the transformed potential (60) in terms of
the variable $x$. This gives

$\displaystyle\hat{V}_{E}(x)$ $\displaystyle=$ $\displaystyle
E-x^{-q-2}\left\\{\frac{(q+1)^{2}}{8~{}p}-\frac{\hat{U}_{E}\left[\log(x)\right]}{2~{}p}\right\\}.$
(62)

Upon implementing the settings (61) and (62) in our equation (LABEL:sse2xt),
we obtain its explicit form as

$\displaystyle\frac{1}{2~{}m(x)}~{}\hat{\Psi}^{\prime\prime}(x)+\Bigg{[}-\frac{m^{\prime}(x)}{2~{}m(x)^{2}}+\frac{\nu}{m(x)~{}x}\Bigg{]}\hat{\Psi}^{\prime}(x)+\Bigg{[}-\frac{\nu}{2~{}m(x)~{}x^{2}}-\frac{\nu~{}m^{\prime}(x)}{2~{}m(x)^{2}~{}x}+\frac{\nu~{}\delta}{2~{}m(x)~{}x^{2}}+$
$\displaystyle+~{}\frac{\nu~{}\delta~{}m^{\prime}(x)}{2~{}m(x)^{2}~{}x}+E-\hat{V}_{E}(x)\Bigg{]}\hat{\Psi}(x)~{}=~{}0,$

recall that the mass is restricted to the form shown in (53). If the function
$\hat{\Psi}$ has the correct parity as indicated by the value of $\delta$,
then it is a solution of the transformed Dunkl-Schrödinger equation (51). In
summary, we have established Darboux transformations between the Dunkl-
Schrödinger partner equations (17) and (50), provided the position-dependent
mass function obeys the restriction (53).

###  5.3 Applications

In order to illustrate the formalism developed in the previous paragraphs, we
will now perform examples of second-order Darboux transformations with the
goal of generating solvable Dunkl-Schrödinger equations with an energy-
dependent potentials that admits solutions of bound state type. While in the
first application the standard algorithm is used, the second application is
devoted to the confluent algorithm. Our starting point for both applications
is the Dunkl-Schrödinger equation (12) that we will now equip with the
following settings

$\displaystyle m(x)~{}=~{}\frac{1}{2}\qquad\qquad\qquad
V_{E}(x)~{}=~{}\frac{x^{2}}{E}.$ (63)

There are several motivations for choosing a constant mass function. Besides
resulting in comparably simple calculations, our constant mass is a particular
case of both (25) and (53) for the settings $p=1/2$, $q=0$. We will comment on
different mass functions at the end of this paragraph. Now, upon substituting
(63), the Dunkl-Schrödinger equation (12) takes the specific form (19).
Furthermore, a constant mass implies even parity, such that we have $\mu=1$.
In the next step we observe that the above choice for our mass function
matches (53) if we set $p=1/2$ and $q=0$. After incorporation of these
settings into equation (19), we obtain

$\displaystyle\Psi^{\prime\prime}(x)+\frac{2~{}\nu}{x}~{}\Psi^{\prime}(x)+\Bigg{(}\frac{\delta~{}\nu-\nu}{x^{2}}+E-\frac{x^{2}}{E}\Bigg{)}\Psi(x)~{}=~{}0,$
(64)

which is a special case of (LABEL:sse0x). A particular solution to this
equation can be given in the form

$\displaystyle\Psi(x)$ $\displaystyle=$
$\displaystyle\exp\left(\frac{x^{2}}{2\sqrt{E}}\right)x^{\frac{1}{2}-\nu+\frac{1}{2}\sqrt{1-4\delta\nu+4\nu^{2}}}~{}L_{-\frac{1}{2}+\frac{E^{\frac{3}{2}}}{4}-\frac{1}{4}\sqrt{1-4\delta\nu+4\nu^{2}}}^{\frac{1}{2}\sqrt{1-4\delta\nu+4\nu^{2}}}\left(\frac{x^{2}}{\sqrt{E}}\right),$
(65)

note that $L$ stands for the associated Laguerre function [1]. Recall that the
function (65) is a solution of the Dunkl-Schrödinger equation (12) only if it
has parity that matches the correct value of $\delta$. Let us now apply our
point transformation (20) that takes (64) into a form that is needed for
applying our Darboux transformation. We choose

$\displaystyle\Phi(y)~{}=~{}\exp\left[\frac{1}{2}~{}y~{}(2~{}\nu-1)\right]\Psi\left[\exp(y)\right]\qquad\qquad\qquad
x(y)~{}=~{}\exp(y).$ (66)

Observe that the left part is obtained from (20) by inserting the present
settings (63) for $\mu=1$, $p=1/2$ and $q=0$. Upon inserting the explicit form
(65), the function $\Phi$ is given by

$\displaystyle\Phi(y)$ $\displaystyle=$
$\displaystyle\exp\left[-\frac{1}{2~{}\sqrt{E}}~{}\exp(2~{}y)+\frac{y}{2}\sqrt{1-4~{}\delta~{}\nu+4~{}\nu^{2}}\right]L_{-\frac{1}{2}+\frac{E^{\frac{3}{2}}}{4}-\frac{1}{4}\sqrt{1-4\delta\nu+4\nu^{2}}}^{\frac{1}{2}\sqrt{1-4\delta\nu+4\nu^{2}}}\left[\frac{\exp(2~{}y)}{\sqrt{E}}\right].$

As such, it is a particular solution of the equation

$\displaystyle\Phi^{\prime\prime}(y)+\left[\delta~{}\nu-\nu^{2}-\frac{1}{4}+E~{}\exp(2~{}y)-\frac{\exp(4~{}y)}{E}\right]\Phi(y)$
$\displaystyle=$ $\displaystyle 0,$ (68)

which is a special case of (55). As expected, the Dunkl parameters $\delta$
and $\nu$ appear as a constant term that is independent of the variable $y$.
Let us follow the definitions from (56) and (57) by setting

$\displaystyle\epsilon~{}=~{}\delta~{}\nu-\nu^{2}\qquad\qquad\qquad
U_{E}(y)~{}=~{}\frac{1}{4}-E~{}\exp(2~{}y)+\exp(4~{}y).$ (69)

These definitions can be implemented in equation (68) that takes the
Schrödinger-like form (55). Before we proceed with the Darboux transformation,
let us recall a remark we made in section 4 about conditions (23) and
(LABEL:con0). According to the remark, the same function $U_{E}$ can be
constructed from various choices of the parameters $m$, $V_{E}$ and the
coordinate change $x$. Let us now demonstrate this using the $U_{E}$ given in
(69), recall that we obtained the latter function through the settings (63).
If we replace these settings by

$\displaystyle m(x)~{}=~{}\frac{1}{2}~{}x^{2}\qquad\qquad\qquad
V_{E}(x)~{}=~{}E+\frac{1}{E}-\frac{E}{x^{2}}-\frac{2}{x^{4}},$ (70)

and perform our point transformation (66), we obtain the equation

$\displaystyle\Phi^{\prime\prime}(y)+\left[3~{}\delta~{}\nu-\nu^{2}-\frac{1}{4}+E~{}\exp(2~{}y)-\frac{\exp(4~{}y)}{E}\right]\Phi(y)$
$\displaystyle=$ $\displaystyle 0.$ (71)

We observe that up to a numerical factor, this equation has the same form as
its counterpart (68). We can match both forms by formally redefining the
parameter $\nu$ in (71) by means of new parameters $\bar{\delta}$ and
$\bar{\nu}$ as

$\displaystyle\nu$ $\displaystyle=$
$\displaystyle\frac{3~{}\delta}{2}+\frac{1}{2}~{}\sqrt{9-4~{}\bar{\delta}~{}\bar{\nu}+4~{}\bar{\nu}^{2}}.$

Substitution renders our equation (71) in the form

$\displaystyle\Phi^{\prime\prime}(y)+\left[\bar{\delta}~{}\bar{\nu}-\bar{\nu}^{2}-\frac{1}{4}+E~{}\exp(2~{}y)-\frac{\exp(4~{}y)}{E}\right]\Phi(y)$
$\displaystyle=$ $\displaystyle 0,$

which coincides with (68) up to parameter naming. Hence, the settings (63) for
a constant mass system and (70) for a position-dependent mass system yield the
same equation. We can therefore consider the results obtained by application
of the Darboux transformation in the subsequent paragraphs as valid for both
the constant mass case as well as certain position-dependent mass scenarios.
We are now ready to perform our Darboux transformation, such that we must
distinguish between the standard and the confluent algorithm.

####  5.3.1 Standard Darboux algorithm

This algorithm requires that we provide two transformation functions
$u_{1},u_{2}$ that solve (55) for the present settings. We choose these
functions as special cases of our soution (LABEL:phiini) as

$\displaystyle
u_{1}(y)~{}=~{}\Phi(y)_{\epsilon=\epsilon_{1}}\qquad\qquad\qquad
u_{2}(y)~{}=~{}\Phi(y)_{\epsilon=\epsilon_{2}}.$ (72)

In order to keep our calculations transparent, we want these functions to take
a relatively simple form. This can be achieved by selecting the parameters
$\epsilon_{1}$ and $\epsilon_{2}$ as

$\displaystyle\epsilon_{1}~{}=~{}\frac{1}{4}\qquad\qquad\qquad\epsilon_{2}~{}=~{}-\frac{3}{4}.$

Substitution of these values into our transformation functions (72) gives
after simplifying and collecting terms

$\displaystyle u_{1}(y)$ $\displaystyle=$
$\displaystyle\exp\left[-\frac{1}{2~{}\sqrt{E}}~{}\exp(2~{}y)\right]L_{\frac{E^{\frac{3}{2}}}{4}-\frac{1}{2}}\left[\frac{\exp(2~{}y)}{\sqrt{E}}\right]$
(73) $\displaystyle u_{2}(y)$ $\displaystyle=$
$\displaystyle\exp\left[y-\frac{1}{2~{}\sqrt{E}}~{}\exp(2~{}y)\right]L_{\frac{E^{\frac{3}{2}}}{4}-1}^{1}\left[\frac{\exp(2~{}y)}{\sqrt{E}}\right].$
(74)

The transformed solution $\hat{\Phi}$ is obtained from the rule (3) for the
case $n=2$, that is, we have

$\displaystyle\hat{\Phi}(y)$ $\displaystyle=$
$\displaystyle\frac{W_{u_{1},u_{2},\Phi}(y)}{W_{u_{1},u_{2}}(y)}.$ (75)

Furthermore, the transformed potential $\hat{U}_{E}$ reads

$\displaystyle\hat{U}_{E}(y)$ $\displaystyle=$ $\displaystyle
U_{E}(y)-2~{}\frac{d^{2}}{dx^{2}}\log\left[W_{u_{1},u_{2}}(y)\right],$ (76)

recall that $U_{E}$ is given in (69). We omit to state the explicit form of
the functions (75) and (76) due to their length, but we will give particular
cases below after having reverted the point transformation. The solution (75)
and the potential (76) enter in equation (58), which in the next step we map
back to Dunkl-Schrödinger form by the inverted version of (66). It is given by

$\displaystyle\hat{\Psi}(x)~{}=~{}x^{\frac{1}{2}-\nu}~{}\hat{\Phi}\left[\log(x)\right]\qquad\qquad\qquad
y(x)~{}=~{}\log(x).$ (77)

As mentioned above, the explicit form of this function is very large, such
that we restrict ourselves to stating a special case here for the settings
$|\nu=5/2$, $\delta=1$, and $E=6$. Substitution yields

$\displaystyle\hat{\Psi}(x)_{|\nu=5/2,\delta=-1,E=4}~{}=~{}\frac{\exp\left(-\frac{x^{2}}{2}\right)I_{0}\left(\frac{x^{2}}{2}\right)}{(24-6~{}x^{2}+x^{4})~{}I_{0}\left(\frac{x^{2}}{2}\right)-x^{2}~{}(x^{2}-4)~{}I_{1}\left(\frac{x^{2}}{2}\right)},$
(78)

where $I$ stands for the modified Bessel function of the first kind [1]. This
solution represents a bound state for our transformed system, as we will
demonstrate below. Next, we compute the transformed potential $\hat{V}_{E}$
from (62) by inserting the present settings $p=1/2$ and $q=0$, along with
(76). We obtain

$\displaystyle\hat{V}_{E}(x)$ $\displaystyle=$ $\displaystyle
E-\frac{1}{4~{}x^{2}}+\frac{1}{x^{2}}~{}\hat{U}_{E}\left[\log(x)\right].$ (79)

Note that $\hat{V}_{E}$ depends on the energy $E$. As in case of the
transformed solution, we will give a particular potential (79) for $E=4$.
Evaluation yields

$\displaystyle\hat{V}_{4}(x)$ $\displaystyle=$
$\displaystyle\Bigg{\\{}4~{}\Bigg{[}(24-6~{}x^{2}+x^{4})~{}I_{0}\left(\frac{x^{2}}{2}\right)-x^{2}~{}(x^{2}-4)~{}I_{1}\left(\frac{x^{2}}{2}\right)\Bigg{]}^{2}\Bigg{\\}}^{-1}\times$
$\displaystyle\times$
$\displaystyle\Bigg{[}(9216-7488~{}x^{2}+1920~{}x^{4}-180~{}x^{6}+4~{}x^{8}+x^{10})~{}I_{0}\left(\frac{x^{2}}{2}\right)^{2}-$
$\displaystyle-$ $\displaystyle
2~{}x^{2}~{}(x-2)~{}(x+2)~{}(x^{2}+16)~{}(24-6~{}x^{2}+x^{4})~{}I_{0}\left(\frac{x^{2}}{2}\right)I_{1}\left(\frac{x^{2}}{2}\right)+$
$\displaystyle+$ $\displaystyle
x^{2}~{}(x^{2}-4)^{2}~{}(72+16~{}x^{2}+x^{4})~{}I_{1}\left(\frac{x^{2}}{2}\right)^{2}\Bigg{]}.$

Figure 4 shows a graph of this function, along with the initial potential from
(63) and other transformed special cases for specific choices of the
stationary energy.

Figure 4: Graphs of the initial potential from (63) and transformed
counterparts (79) for the settings (63) and the transformation functions (73)
and (74).

The solution (77) and potential (79) enter in the transformed Dunkl-
Schrödinger equation that is the partner to (64). It reads

$\displaystyle\hat{\Psi}^{\prime\prime}(x)+\frac{2~{}\nu}{x}~{}\hat{\Psi}^{\prime}(x)+\Bigg{[}\frac{\delta~{}\nu-\nu}{x^{2}}+E-\hat{V}_{E}(x)\Bigg{]}\hat{\Psi}(x)~{}=~{}0.$
(80)

Let us now recall that the function (77) is an admissible solution for this
equation only if its parity matches the correct value of $\delta$. Inspection
of the explicit form of (77) reveals that parity is determined by a single
monomial term. We have

$\displaystyle\hat{\Psi}(x)$ $\displaystyle\stackrel{{\scriptstyle
x<<1}}{{=}}$ $\displaystyle
x^{\frac{1}{2}-\nu+\frac{1}{2}\sqrt{1-4\delta\nu+4\nu^{2}}},$ (81)

where irrelevant multiplicative constants were discarded. The exponent in (81)
must be an odd integer if $\delta=-1$, and it must be an even integer for
$\delta=1$. Let us evaluate the exponent for the first case. We find

$\displaystyle\frac{1}{2}-\nu+\frac{1}{2}~{}\sqrt{1+4~{}\nu+4~{}\nu^{2}}$
$\displaystyle=$
$\displaystyle\frac{1}{2}-\nu+\frac{1}{2}~{}|1+2~{}\nu|~{}=~{}1,$

note that in the last step we used the restriction $\nu>-1/2$. In the second
case $\delta=1$, the exponent (81) takes the form

$\displaystyle\frac{1}{2}-\nu+\frac{1}{2}~{}\sqrt{1-4~{}\nu+4~{}\nu^{2}}$
$\displaystyle=$ $\displaystyle\frac{1}{2}-\nu+\frac{1}{2}~{}|-1+2~{}\nu|$
(84) $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{lll}1-2~{}\nu&\mbox{if}&\nu<\frac{1}{2}\\\\[4.30554pt]
0&\mbox{if}&\nu\geq\frac{1}{2}\end{array}\right\\}.$

Since $\nu>-1/2$, the exponent cannot attain any even integer value for
$\nu<1/2$. Thus, we must have $\nu\geq 1/2$, in which case the exponent
vanishes. Consequently, the correct parity of our transformed solution (77) is
established if we choose $\nu\geq 1/2$. In the next step we will establish
normalizability of our solutions (77) for certain values of the energy. To
this end, we recall that the modified norm of a solution to the Dunkl-
Schrödinger equation (12) is given by (14). In order to be meaningful, the
integral this modified norm must exist, and the integrand must be nonnegative.
Existence of the integral can be achieved by choosing the stationary energy
such that the first argument of the Laguerre function in (65) becomes a
nonnegative integer. This condition reads

$\displaystyle-\frac{1}{2}+\frac{E^{\frac{3}{2}}}{4}-\frac{1}{4}\sqrt{1-4~{}\delta~{}\nu+4~{}\nu^{2}}$
$\displaystyle=$ $\displaystyle n,$

where $n$ is a nonnegative integer. Solving with respect to $E$ gives

$\displaystyle E$ $\displaystyle=$
$\displaystyle\left(4~{}n+2+\sqrt{1-4~{}\delta~{}\nu+4~{}\nu^{2}}\right)^{\frac{2}{3}}.$
(85)

As a consequence of this choice for the stationary energy, the Laguerre
function will degenerate to a polynomial, which renders (65) bounded on the
whole real line. This property persists under our Darboux transformation, such
that the transformed solution (77) is also bounded on the whole real line.
This extends to the whole integrand in (14), such that the integral exists. In
addition we have

$\displaystyle 1-\frac{\partial V_{E}(x)}{\partial E}$ $\displaystyle\geq$
$\displaystyle 0.$

We do not verify these conditions explicitly here due to the length of the
involved expressions. Instead we show graphs of normalized probability
densities associated with our transformed solution (77) for different
parameter values in figure 5.

Figure 5: Graphs of the normalized transformed probability density for the
parameter settings $\delta=-1$ and $\nu=5/2$.

Recall that the probability density for a solution $\Psi$ is given by (13).
Inspection of the figure indicates that the probability densities are
nonnegative functions. Note that the energy values chosen in the figures 4 and
5 are taken from (85) for $n=0,1,2$. Since we know now that our transformed
solution (77) represents bound states if the stationary energy is chosen from
(85), let us present graphs of particular cases. Figure 6 and figure 7 show
graphs of the solutions (77) for the odd-parity case and the even-parity case,
respectively.

Figure 6: Graphs of the normalized transformed probability density for the
parameter settings $\delta=-1$ and $\nu=5/2$. Figure 7: Graphs of the
normalized transformed probability density for the parameter settings
$\delta=1$ and $\nu=7/2$.

Observe that one of the graphs shown in figure 6 pertains to the function
(78).

####  5.3.2 Confluent Darboux algorithm

We go back to our initial equation (68) that we consider to be in the form
(55) due to our settings (69). In contrast to its standard counterpart, the
confluent algorithm of our Darboux transformations requires transformation
functions that solve the system (5), (6). In the present case of a second-
order transformation the latter system reads

$\displaystyle
u_{1}^{\prime\prime}(y)+\left[\epsilon_{1}-U_{E}(y)\right]u_{1}(y)$
$\displaystyle=$ $\displaystyle 0$ (86) $\displaystyle
u_{2}^{\prime\prime}(y)+\left[\epsilon_{1}-U_{E}(y)\right]~{}u_{2}(y)$
$\displaystyle=$ $\displaystyle-u_{1}(y),$ (87)

recall that $U_{E}$ is given in (69). In order to solve the system (86), (87),
we observe that the first equation is formally the same as (68), such that a
solution $u_{1}$ can be taken from (LABEL:phiini). After inserting the
definition of $\epsilon$ in (69) and replacing $\epsilon$ by $\epsilon_{1}$,
we find

$\displaystyle u_{1}(y)$ $\displaystyle=$
$\displaystyle\exp\left[-\frac{1}{2~{}\sqrt{E}}~{}\exp(2~{}y)+\frac{y}{2}\sqrt{1-4~{}\epsilon_{1}}\right]L_{-\frac{1}{2}+\frac{E^{\frac{3}{2}}}{4}-\frac{1}{4}\sqrt{1-4\epsilon_{1}}}^{\frac{1}{2}\sqrt{1-4\epsilon_{1}}}\left[\frac{\exp(2~{}y)}{\sqrt{E}}\right].$
(88)

In the present application we choose the parameter $\epsilon_{1}$ as

$\displaystyle\epsilon_{1}$ $\displaystyle=$ $\displaystyle-2.$

Upon substitution into (88), we obtain the explicit form

$\displaystyle u_{1}(y)$ $\displaystyle=$
$\displaystyle\exp\left[-\frac{1}{2~{}\sqrt{E}}~{}\exp(2~{}y)+\frac{3~{}y}{2}\right]L_{\frac{E^{\frac{3}{2}}}{4}-\frac{5}{4}}^{\frac{3}{2}}\left[\frac{\exp(2~{}y)}{\sqrt{E}}\right].$
(89)

In the next step we need to find the remaining transformation function $u_{2}$
by solving equation (87), which can be done by applying integral or
differential formulas [22]. According to these formulas, the simplest way of
constructing a particular solution to (87) is given by

$\displaystyle u_{2}(y)$ $\displaystyle=$
$\displaystyle\frac{\partial\Phi(y)}{\partial\epsilon}_{\mid\epsilon=\epsilon_{1}},$
(90)

where as before $\Phi$ was taken from (LABEL:phiini). We omit to state the
explicit form of the resulting transformation function $u_{2}$ after
application of the derivative. From this point on the confluent algorithm of
the Darboux transformation proceeds in the same way as the standard case.
Substitution of the functions (89), (89), (LABEL:phiini) into (75) and (76)
yields the transformed solution and potential, respectively. After reversing
our point transformation, we obtain a solution (77) and an associated
potential (79) for equation (80). If the latter solution has the correct
parity as indicated by the parameter $\delta$, then it solves the transformed
Dunkl-Schrödinger equation (51) for the present settings in (63). Due to the
parametric derivatives, the transformed solution and potential do not allow
for simplification, such that their explicit forms are long and involved. As a
consequence, it becomes difficult to analyze the probability density
associated with the solutions in explicit form. For this reason, we restrict
ourselves here to show graphs of the transformed potential only, see figure 8.

Figure 8: Graphs of the initial potential from (63) and transformed
counterparts (79) for the settings (63) and the transformation functions (89)
and (90).

##  6 Concluding remarks

When applying solution-generating methods like point or Darboux
transformations, the presence of an energy-dependent potential imposes an
additional constraint on the transformed system, in particular when looking
for bound states. This constraint arises from the norm that involves the
parametric derivative of the initial system’s potential. In general,
normalizability is not preserved by either transformation. However, if we
assume that the coordinate change $x=x(y)$ does not depend on the energy, then
we can derive the following relation from (23):

$\displaystyle 1-\frac{\partial U_{E}(y)}{\partial E}$ $\displaystyle=$
$\displaystyle 2~{}m[x(y)]~{}[x^{\prime}(y)]^{2}\left[1-\frac{\partial
V_{E}(y)}{\partial E}\right].$

We observe that if the mass is nonnegative and the coordinate change is real-
valued, then normalizability is preserved by the point transformation (20).
For the Darboux transformation an analogous statement cannot be obtained in
such a simple manner because the transformed potential depends on the
transformation functions that can exhibit a complicated dependence on the
energy. A further comment on the Darboux transformations and the restriction
(53) is in order. This restriction on the mass can be avoided by assigning the
parity parameter $\delta$ a fixed value of $-1$ or $1$. Since in this case the
transformed potential does not depend on this parameter anymore, the
restriction (53) is not necessary. On the other hand, the Darboux-transformed
equation must then be cast in a form that allows mapping it back to Dunkl-
Schrödinger form, and parity of the transformed solutions must be established.

##### Data availability statement.

The data that supports the findings of this study are available within the
article.

## References

  * [1] M. Abramowitz and I. Stegun, $"$Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables$"$, (Dover Publications, New York, 1964)
  * [2] J.P. Anker, $"$An introduction to Dunkl theory and its analytic aspects$"$, In: Filipuk, G., Haraoka, Y., Michalik, S. (eds) Analytic, Algebraic and Geometric Aspects of Differential Equations, Trends in Mathematics, Birkhäuser, Cham, 2017
  * [3] D. Barrios Rolania, J.C. Garcia-Ardila, and D. Manrique, $"$On the Darboux transformations and sequences of p-orthogonal polynomials$"$, Appl. Math. Comput. 382 (2020), 125337
  * [4] R. Bravo and M.S. Plyushchay, $"$Position-dependent mass, finite-gap systems, and supersymmetry$"$, Phys Rev. D 93 (2016), 105023
  * [5] F. Bouzeffour, A. Nemri, A. Fitouhi, and S. Ghazouani, $"$On harmonic analysis related with the generalized Dunkl operator$"$, 23 (2012), 609
  * [6] W.S. Chung and H. Hassanabadi, $"$One-dimensional quantum mechanics with Dunkl derivative$"$, Mod. Phys. Lett. A 34 (2019), 1950190
  * [7] F. Cooper, A. Khare, and U. Sukhatme, $"$Supersymmetry and Quantum Mechanics$"$, Phys.Rept. 251 (1995), 267
  * [8] B.G. da Costa, G.A.C. da Silva, and I.S. Gomez, $"$Supersymmetric quantum mechanics and coherent states for a deformed oscillator with position-dependent effective mass$"$, J. Math. Phys. 62 (2021), 092101
  * [9] C.F. Dunkl, $"$Differential-difference operators associated to reflection groups”, Transactions of the American Mathematical Society, 311 (1989), 167
  * [10] P. Etingof, $"$Calogero-Moser systems and representation theory$"$, Zurich Lect. Adv. Math. 4, Eur. Math. Soc., Zurich, 2007
  * [11] M. Feigin and M. Vrabec, $"$Intertwining operator for AG2 Calogero-Moser-Sutherland system$"$, J. Math. Phys. 60 (2019), 073503
  * [12] M. Feigin and T. Hakobyan, $"$On Dunkl angular momenta algebra$"$, JHEP 11 (2015), 107
  * [13] V.X. Genest, M.E.H. Ismail, L. Vinet, and A. Zhedanov, $"$The Dunkl oscillator in the plane II : representations of the symmetry algebra$"$, Commun. Math. Phys. 329 (2014), 999
  * [14] V.X. Genest, M.E.H. Ismail, L. Vinet, and A. Zhedanov, $"$The Dunkl oscillator in the plane I : superintegrability, separated wavefunctions and overlap coefficients$"$, J. Phys. A 46 (2013), 145201
  * [15] S. Ghazouani, I. Sboui, M.A. Amdouni, and M. Ben El Hadj Rhouma, $"$The Dunkl-Coulomb problem in three-dimensions: energy spectrum, wave functions and h-spherical harmonics$"$, J. Phys. A 52 (2019), 225202
  * [16] S. Karthiga, V. Chithiika Ruby, and R.M. Senthilvelan, $"$An inclusive SUSY approach to position dependent mass systems$"$, Ann. Phys. 382 (2018), 1645
  * [17] Y. Luo, S. Tsujimoto, L. Vinet, and A. Zhedanov, $"$Dunkl-supersymmetric orthogonal functions associated with classical orthogonal polynomials$"$, J. Phys. A 53 (2020), 085205
  * [18] H. Mejjaoli, $"$Nonlinear generalized Dunkl-wave equations and applications$"$, J. Math. Anal. Appl. 375 (2011), 118
  * [19] R.D. Mota, D. Ojeda-Guillen, M. Salazar-Ramirez, and V.D. Granados, $"$Exact Solutions of the 2D Dunkl–Klein–Gordon Equation: The Coulomb Potential and the Klein–Gordon Oscillator$"$, Mod. Phys. Lett. A 36 (2021), 2150171
  * [20] R.D. Mota, D. Ojeda-Guillen, M. Salazar-Ramirez, and V.D. Granados, $"$Exact Solution of the Relativistic Dunkl Oscillator in (2+1) Dimensions$"$, Ann. Phys. 411 (2019), 167964
  * [21] A. Schulze-Halberg, $"$Generalized Dunkl-Schrodinger equations: solvable cases, point transformations, and position-dependent mass systems$"$, Phys. Scr. 97 (2022), 085213
  * [22] A. Schulze-Halberg, $"$Regularity conditions for transformed potentials in the confluent supersymmetry algorithm$"$, Int. J. Mod. Phys. A 33 (2018), 1850214
  * [23] A. Schulze-Halberg and O. Yesiltas, $"$Generalized Schrödinger equations with energy-dependent potentials: formalism and applications$"$, J. Math. Phys. 59 (2018), 113503
  * [24] J.F. van Diejen and L. Vinet (eds.), $"$Calogero-Moser-Sutherland Models$"$, in CRM Series in Mathematical Physics, (Springer, New York, 2000)
  * [25] https://functions.wolfram.com/HypergeometricFunctions/Hypergeometric1F1/16/01/01/0001/