Confluent Supersymmetric Partners of Quantum Systems Emerging from the Spheroidal Equation

Abstract

We construct confluent supersymmetric partners of quantum systems that emerge from the spheroidal equation. Properties of the systems and of their transformed counterparts are discussed.

1. Introduction

The supersymmetry (SUSY) formalism is a method for generating solvable models in quantum mechanics. Mathematically based on the Darboux transformation [1–3], it converts a known quantum system into a transformed counterpart, the so-called SUSY partner. As this formalism is compatible with basically any exactly-solvable quantum model governed by the Schrödinger equation and many systems described by the Dirac equation, a large amount of applications can be found in the literature. An overview and details on the general topic can be found in one of the reviews [4–6] and references therein. The formalism of SUSY can be split into two different versions, commonly referred to as standard and confluent algorithm, respectively. While the vast majority of applications uses the standard algorithm, its confluent counterpart has received considerably less attention. As a consequence, many of its properties are not yet well understood. Recent works on the confluent algorithm include its application to the inverted oscillator [7] and to the Dirac equation for pseudoscalar potentials [8]. Further results concern its equivalence with the zero-mode supersymmetry scheme [9], a representation using a differential formula [10,11] and higher-order confluent SUSY transformations [12]. The purpose of the present note is to apply the confluent SUSY algorithm to models that are governed by the spheroidal equation [13], a particular case of the confluent Heun class [14]. The spheroidal equation has many applications in physics, including a sextic oscillator system with centrifugal barrier [15], a quantum model featuring two Coulombian centers [16], the dynamical theory of zonal tides in a global ocean [17], just to mention a few. While solutions of the spheroidal equation cannot be given in closed form, for certain parameter settings one can construct them through a series expansion in terms of Gegenbauer polynomials [18]. In this work we will construct such series solutions for two systems that are associated with the spheroidal equation. The first system is a generalization of the Mathews-Lakshmanan oscillator [19–21], while the second features a potential of Pöschl-Teller type that contains an aditional cosine term. After obtaining the solutions, we will generate confluent SUSY partners of the two systems. The remainder of this paper is organized as follows. In Section 2 we give a summary of the confluent SUSY algorithm. Section 3 is devoted to the spheroidal equation, the algorithm for the construction of its solutions and a discussion of their properties. In Section 4 the findings are applied to the two particular quantum systems mentioned above, resulting in the generation of their confluent SUSY partner systems.

Implementation notes: All results on spectral values and solutions that are presented subsequently were obtained by implementing the algorithm described in Section 2.2 using MATLAB (The MathWorks, Inc., Natick, MA, USA). Calculations were done at 40-digit precision. All digits of numerical values stated throughout this work are significant.

2. The Confluent SUSY Algorithm

We will now give a brief review of the confluent SUSY formalism. For further details, particularly on spectral design, we refer the reader to [11,22,23] and references therein. We start out by considering an initial equation that has the general form

$$f{\mathrm{\Psi }}^{″}+g{\mathrm{\Psi }}^{\prime }+(E-V_0)\mathrm{\Psi }=0,$$

(1)

where the coefficient functions f, g are required to be sufficiently smooth and the potential V₀ must be continuous. Next, for a positive integer N ≥ 2, assume that N solutions u₀, v_j, j = 0,…, N − 2, of the equation

$$f{u}^{″}+g{u}^{\prime }+(\mathrm{\lambda }-V_0)u=0,$$

(2)

are given, introducing a constant λ, commonly referred to as factorization energy. Note that the solutions u₀, v_j, j = 0,…, N − 2 do not have to be normalizable in the L²-sense, that is, they are not required to represent bound states. We will return to this point at the end of this paragraph. Now we define functions u_j, j = 1,…, N − 1, as follows

$$u_j={\displaystyle \sum _{k=0}^{j-1}\frac{{\partial }^k{v}_k}{\partial {\mathrm{\lambda }}^k}+\frac{1}{j!}\frac{{\partial }^j{u}_0}{\partial {\mathrm{\lambda }}^j},j=1,\dots ,N-1.}$$

(3)

As an example, let us mention that in the simplest case N = 2 we obtain from Equation (3) a single function u₁, given by

$$u_1=v_0+\frac{\partial u_0}{\partial \mathrm{\lambda }},$$

(4)

where v₀ and u₀ are solutions of Equation (2). We will now use the functions u_j, j = 0,…, N − 1, called transformation functions, to set up a SUSY transformation of order N, applied to the solution Ψ of Equation (1). It reads

$$\mathrm{\Phi }=f^{\frac{N}{2}}\frac{W_{u_0,\dots ,u_{N-1,}\mathrm{\Psi }}}{W_{u_0,\dots ,u_{N-1}}}.$$

(5)

Here, the symbol W stands for the Wronskian of the functions in its index. The Equation (5) is a solution of the supersymmetric partner equation to Equation (1), given by

$$f{\mathrm{\Phi }}^{″}+g{\mathrm{\Phi }}^{\prime }+(E-V_1)\mathrm{\Phi }=0.$$

(6)

The new potential V₁ is linked to its initial counterpart V₀ in Equation (1) by means of the following relation [24]

$$V_1=V_0+N(N-2)\left[{\left(\frac{f^{\prime }}{2f}\right)}^2-\frac{f^{″}}{2}\right]+\frac{Ng{f}^{\prime }}{2f}-N{g}^{\prime }-\frac{f^{\prime }{W}_{u_0,u_1,\dots ,u_{N-1}}^{\prime }}{W_{u_0,u_1,\dots ,u_{N-1}}}-2f\frac{d}{dx}\frac{W_{u_0,u_1,\dots ,u_{N-1}}^{\prime }}{W_{u_0,u_1,\dots ,u_{N-1}}}.$$

(7)

The functions V₀ and V₁ are usually called supersymmetric partner potentials. There are two principal differences between the confluent SUSY algorithm and its standard counterpart, the first of which contains the transformation order. The confluent algorithm uses a SUSY transformation of order two or higher, while in the conventional case first-order transformations are possible. In other words, there is no first-order confluent SUSY transformation. In applications, such as the recent works in optics [25,26], one must consider SUSY transformations of at least order two, before the confluent SUSY algorithm becomes applicable. The second principal difference between the confluent SUSY algorithm and its standard counterpart lies in the definition of the transformation functions u_j, j = 0,…, N − 1, affecting both transformed solution (5) and potential (7). In the standard framework of SUSY each transformation function is a solution of the initial Equation (1), associated to its own factorization energy, while the confluent algorithm requires a single factorization energy and defines the transformation functions through the differential formula Equation (3). Note that the latter functions can also be expressed through an integral representation. Since this representation is equivalent to Equation (3) and for our purposes less convenient, we omit to state it here. Before we conclude this paragraph, let us comment once more on the transformation functions u_j, j = 0,…, N − 1. The main purpose of these functions is to render the transformed potential (7) real-valued and free of singularities inside its domain. This condition can be fulfilled if the Wronskian that appears in the denominator of the terms in Equation (7), is free of zeros, which essentially relates to linear independence of the transformation functions. However, in the present models that will be considered in the subsequent sections, we only have access to solutions representing bound states. As a consequence, we do not need to consider transformation functions that are of more general type.

3. The Spheroidal Equation

The spheroidal equation, equipped with boundary conditions of Dirichlet type, can be written in the following way

$$(1-y^2){\chi }^{″}-2y{\chi }^{\prime }+\left[A-h^2{y}^2+\frac{M^2}{y^2-1}\right]\chi =0,y\in (-1,1)$$

(8)

$$\chi (-1)=\chi (1)=0,$$

(9)

where A, h and M are constants. Comparison with Equation (1) shows that Equation (8) is a special case. While Equation (8) cannot be solved in closed form, under certain conditions it admits series solutions in terms of Gegenbauer polynomials. For the sake of convenience, we will now briefly summarize how to construct such series solutions. For further details, the reader may refer to [18]. It is known that there exists an infinite, discrete and ascending sequence of characteristic values A = A_j, j ∈ ℕ ∪ {0}, for each of which the problem Equation (8), Equation (9) admits a solution that has exactly j zeros inside (−1, 1). The solutions form an orthogonal set with respect to the usual inner product of the Hilbert space L²(−1, 1). All solutions have parity, but do not feature a closed-form representation. Instead, they can be written as a series expansion in terms of Gegenbauer polynomials T [27] via one of the formulas

$${\chi }_{2j}={(1-y^2)}^{\frac{M}{2}}{\displaystyle \sum _{n=0}^{\infty }{(d_{2n})}_{|A=A_{2j}}}{T}_{2n}^M(y)$$

(10)

$${\chi }_{2j+1}={(1-y^2)}^{\frac{M}{2}}{\displaystyle \sum _{n=0}^{\infty }{(d_{2n+1})}_{|A=A_{2j+1}}}{T}_{2n+1}^M(y),$$

(11)

for j ∈ ℕ ∪ {0}, where Equation (10) and Equation (11) represent even and odd parity functions, respectively. The expansion coefficients in the series are determined through a three-term recursion relation that reads

$$a_n{d}_{n-2}+b_n{d}_n+c_n{d}_{n+2}=0,n\in \mathbb{N}\cup \left\{0\right\},d_{-2}=d_{-1}=0,$$

(12)

where the following abbreviations are involved:

$$a_n=\frac{n(n-1)h^2}{(2n+2M-1)(2n+2M-3)}$$

(13)

$$b_n=\frac{h^2\{2[(n+M)(n+M+1)-M^2]-1\}}{(2n+2M+3)(2n+2M-1))}+(n+M)(n+M+1)-A$$

(14)

$$c_n=\frac{h^2(n+2M+1)(n+2M+2)}{(2n+2M+3)(2n+2M+5)}.$$

(15)

It is important to point out that the series Equations (10) and (11) are in general divergent. Convergence is guaranteed only if the parameter A in Equation (14) takes a characteristic value A = A_j, j ∈ ℕ ∪ {0}. The algorithm for determining these values, as well as the expansion coefficients in the series Equations (10) and (11), consists of two steps.

Step 1. We solve the recursion (12) with respect to d₂/d₀, first substituting n = 0 and the second time n = 2. This gives the following two relations

$$\begin{array}{cc}\frac{d_2}{d_0}=-\frac{b_0}{c_0}& \frac{d_2}{d_0}=\frac{-a_2}{b_2+c_2\frac{d_4}{d_2}}.\end{array}$$

(16)

Combining these two relations, we obtain

$$\begin{array}{cc}-\frac{b_0}{c_0}& =\frac{-a_2}{b_2+c_2\frac{d_4}{d_2}}.\end{array}$$

(17)

Now we will modify the right side of this equation. We shift all indices n in Equation (12) up by two, set n = 4, solve for d₄/d₂, and substitute the result into the right-hand side of Equation (16). An iteration of this process gives rise to the following continued fraction expansion:

$$\begin{array}{cc}-\frac{b_0}{c_0}& =\frac{-a_2}{b_2+c_2\frac{-a_4}{b_4+c_4\frac{-a_6}{b_6+c_6\dots }}}.\end{array}$$

(18)

We stop this iteration after N > 0 steps by reasoning that d_N+2/d_N is very small, so we set d_N+2/d_N = 0. This turns Equation (18) into an algebraic equation that can be solved for the parameter A, yielding approximations for the characteristic values belonging to solutions of even parity. Observe that Equation (18) generally has several solutions, each of which is associated with a characteristic value. Also, the number of solutions increases with the number N of steps that are done. In a similar way we obtain a continued fraction expansion for odd indices:

$$\begin{array}{cc}-\frac{b_1}{c_1}& =\frac{-a_3}{b_3+c_3\frac{-a_5}{b_5+c_5\frac{-a_7}{b_7+c_7\dots }}},\end{array}$$

(19)

that leads to the characteristic values associated with the odd solutions.

Step 2. We set d₀ = 1. Relation (16) then gives d₂, which can be plugged into the recursion relation (12) for n = 2 in order to obtain d₄. The continuation of this process yields the following expansion coefficients d₆, d₈, d₁₀,…. The coefficients for odd indices are constructed in a similar manner, starting out with d₁ = 1. The recursion (12) for n = 1 generates d₃, which can be substituted into (12) for n = 3 in order to give d₅. This iteration produces the subsequent expansion coefficients d₇, d₉, d₁₁,…. The solutions of the boundary-value problem (8), (9) can now be approximated by plugging the expansion coefficients and the characteristic values into the series representations (10), (11).

4. Applications

We will now use point transformations for converting the spheroidal equation, such that it describes different physically interesting models. In particular, we focus on two systems, featuring a double oscillator potential of Mathews-Lakshmanan type and a generalized trigonometric Pöschl-Teller interaction, respectively.

4.1. A Double Oscillator System

The first model that will be considered here can be generated by means of the following point transformation.

$$\chi [y(x)]={(\mu x^2-1)}^{-\frac{1}{4}}\mathrm{\Psi }(x),y(x)=-\sqrt{\mu }x,$$

(20)

where µ > 0 is a real constant. In addition, we make the following definitions in Equation (8)

$$A=\frac{2E}{\mu }+\frac{1}{{\mu }^2}-\frac{1}{\mu }-\frac{1}{4}h=-\frac{\sqrt{q}}{\mu }M=\frac{1}{2\sqrt{\mu }}\sqrt{\mu +\frac{4}{\mu }-4},$$

(21)

for a real constant q. Note that although h can become imaginary, it enters in the spheroidal equation and in the recursion Equations (12)–(15) for the series coefficients merely as the real-valued h², therefore preventing the spectral values from becoming complex. After substituting both Equation (20) and the settings (21) into the boundary-value problem (8), (9), we obtain after some elementary manipulations

$$(1-\mu x^2){\mathrm{\Psi }}^{″}-\mu x{\mathrm{\Psi }}^{\prime }+\left[2E-q{x}^2-\frac{(1-\mu )x^2}{1-\mu x^2}\right]\mathrm{\Psi }=0,x\in \left(-\frac{1}{\sqrt{\mu }},\frac{1}{\sqrt{\mu }}\right)$$

(22)

$$\mathrm{\Psi }\left(-\frac{1}{\sqrt{\mu }}\right)=\mathrm{\Psi }\left(\frac{1}{\sqrt{\mu }}\right)=0.$$

(23)

Note that we applied the point transformation (20) to the initial boundary conditions (9) as well, yielding the result Equation (23). The model described by Equations (22) and (23) can be interpreted as the quantum counterpart of a classical system that is governed by the following Lagrangian

$$ℒ=\frac{{\dot{x}}^2}{2-2\mu x^2}-\frac{q}{2}{x}^2-\frac{(1-\mu )x^2}{2-2\mu x^2},$$

(24)

The specific form of the kinetic term in Equation (24) can be explained by introducing a position-dependent mass function. This function m and the potential V featured in Equation (24) have the explicit form

$$m=\frac{1}{1-\mu x^2}{V}_0=\frac{q}{2}{x}^2+\frac{(1-\mu )x^2}{2-2\mu x^2}.$$

(25)

The interaction V represents a nonlinear, two-parameter double oscillator, containing a harmonic term and a contribution of Mathews-Lakshmanan type. By virtue of its two parameters µ and q, the potential (25) has several well-known limiting cases, a few of which are listed in

Table 1. Limiting cases of the potential V in Equation (25).

Parameter Setting	Potential	Name of Special Case
q = −1, µ → 0	V0 = 0	Free particle
q = 1, µ → 0	V0 = x2	Harmonic oscillator
q = −3, µ → 0	V0 = − x2	Inverted harmonic oscillator
q = 0	$V_0=\frac{(1-\mu )x^2}{2-2\mu x^2}$	Mathews-Lakshmanan oscillator

.

Next, we will show that the point transformation (20) preserves orthogonality and normalizability of the bound-state solutions stemming from the initial system (8), (9), as long as we require the solutions of Equations (22) and (23) to be located in the weighted Hilbert space $L_w^2(-\sqrt{1/\mu },\sqrt{1/\mu })$ weight function

$$w=\sqrt{\frac{1}{1-\mu x^2}}.$$

(26)

Let k, l be two nonnegative integers ands let N be a constant. We obtain by combining orthogonality of the Equations (10) and (11) with the point transformation (20) and (26) the following result

$$N{\delta }_{kl}={\displaystyle \underset{-1}{\overset{1}{\int }}{\chi }_k^{*}{\chi }_{l}dy}={\displaystyle \underset{-\sqrt{\frac{1}{\mu }}}{\overset{\sqrt{\frac{1}{\mu }}}{\int }}{\mathrm{\Psi }}_k^{*}{\mathrm{\Psi }}_l\sqrt{\mu }}\sqrt{\frac{1}{1-\mu x^2}}dx=\sqrt{\mu }{\displaystyle \underset{-\sqrt{\frac{1}{\mu }}}{\overset{\sqrt{\frac{1}{\mu }}}{\int }}{\mathrm{\Psi }}_k^{*}{\mathrm{\Psi }}_l\sqrt{\frac{1}{1-\mu x^2}}}dx.$$

(27)

Combining the very left and the very right term in this chain of equations, we obtain

$${\displaystyle \underset{-\sqrt{\frac{1}{\mu }}}{\overset{\sqrt{\frac{1}{\mu }}}{\int }}{\mathrm{\Psi }}_k^{*}{\mathrm{\Psi }}_l}\sqrt{\frac{1}{1-\mu x^2}}dx=\frac{N}{\sqrt{\mu }}{\delta }_{kl},$$

(28)

implying orthogonality of the set (Ψ_j) with respect to the weighted Hilbert space $L_w^2(-\sqrt{1/\mu },\sqrt{1/\mu })$. Normalizability in the same Hilbert space follows then from

$$\Vert {\mathrm{\Psi }}_k\Vert ={\displaystyle \underset{-\sqrt{\frac{1}{\mu }}}{\overset{\sqrt{\frac{1}{\mu }}}{\int }}{\mathrm{\Psi }}_k^{*}{\mathrm{\Psi }}_k}\sqrt{\frac{1}{1-\mu x^2}}dx=\frac{N}{\sqrt{\mu }}{\displaystyle \underset{-1}{\overset{1}{\int }}{\chi }_k^{*}{\chi }_{k}dy}=\frac{N}{\sqrt{\mu }}\Vert {\chi }_k\Vert .$$

(29)

Now that we have established orthogonality and normalizability of the solutions, we are able to construct an explicit form for them. First, the point transformation (20) gives

$$\mathrm{\Psi }(x)={(\mu x^2-1)}^{\frac{1}{4}}\chi (-\sqrt{\mu x})$$

(30)

In the next step we substitute the series form (10), (11) for the solutions of the spheroidal Equation (8) into Equation (30). Note that due to the settings (21), the characteristic values of A will translate into discrete spectral values of E, denoted by E_j, j a nonnegative integer. After adding an index for taking into account the bound-state character of the solutions, we obtain

$${\mathrm{\Psi }}_{2j}={(1-\mu x^2)}^{\frac{1}{4}+\frac{1}{4\sqrt{\mu }}\sqrt{\mu +\frac{4}{\mu }-4}}{\displaystyle \sum _{n=0}^{\infty }{(d_{2n})}_{|E=E_{2j}}}{T}_{2n}^{\frac{1}{2\sqrt{\mu }}\sqrt{\mu +\frac{4}{\mu }-4}}(-\sqrt{\mu }x)$$

(31)

$${\mathrm{\Psi }}_{2j+1}={(1-\mu x^2)}^{\frac{1}{4}+\frac{1}{4\sqrt{\mu }}\sqrt{\mu +\frac{4}{\mu }-4}}{\displaystyle \sum _{n=0}^{\infty }{(d_{2n+1})}_{|E=E_{2j+1}}}{T}_{2n+1}^{\frac{1}{2\sqrt{\mu }}\sqrt{\mu +\frac{4}{\mu }-4}}(-\sqrt{\mu }x),$$

(32)

We are now ready to compute confluent SUSY partners of our model (22), (23), the governing equation of which is clearly a special case of Equation (8) for the settings

$$f=1-\mu x^{2}g=-\mu x.$$

(33)

For the sake of simplicity, here we will restrict ourselves to transformations of second order, given by the relation (5) for N = 2. The principal reason for this restriction lies in the computational effort that is required for the calculation of the transformed potential (7), in particular for the Wronskians and their derivatives. Each function u_j, j = 0,…, N − 1, is constructed by means of the algorithm shown in Section 4, such that it is represented by a series of the form (31) or (32). Consequently, Wronskians for values of N ≥ 3 contain at least three of such series including at least their second derivatives. This means a large amount of computational work, especially since we must maintain our prescribed numerical precision. For these reasons, here we prefer to use only the case N = 2 (second-order transformations), but we will get back to this issue below. Now we need two transformation functions u₀ and u₁, the first of which must be a solution of Equation (22). After determining u₀ we can generate u₁ by means of Equation (4). While in general we can choose a transformation function u₀ for any factorization energy λ, in the present case only bound-state solutions of the series form (31), (32) are available. Next we focus on the function u₁ as given in Equation (4). Since v₀ must be a solution of Equation (8) for the same energy λ as u₀, we can only pick u₀ itself or a multiple of it. Note that applying the reduction-of-order method for generating another solution of Equation (8) would lead to integrals of rational functions that feature high-degree polynomials in the denominator. Such integrals can in general not be resolved in closed form. Now we observe that the transformation functions u₀ and u₁ enter in the potential (7) only through their Wronskian. Consequently, if v₀ in Equation (4) and u₀ were not linearly independent, we would obtain the same transformed potential as for v₀ = 0. For a nonnegative integer j we make the choice

$$u_1=\frac{\partial u_0}{\partial {\mathrm{\lambda }}_{|\mathrm{\lambda }=E_j}},$$

(34)

while u₀ will be one of the solutions given in Equations (31) and (32), pertaining to the spectral value E_j. We can now perform a confluent SUSY transformation of second order by substituting the settings Equation (34) into Equations (5) and (7). The representation (31) for u₀ and u₁ requires truncation of the series. Despite its convergence being very fast [18], we still have to take into account a large number of terms. For this reason we do not show the explicit calculations that lead to the subsequently shown results.

In the Figures 1–4 we display examples of potentials V₁ that result from a confluent SUSY transformation of second order. These potentials are obtained from Equation (7) by substituting N = 2 and the Wronskians of the transformation functions u₀, u₁. Each of the four plots contained in the latter figures shows the initial potential and two of its confluent SUSY partners at the factorization energies λ = E₀ and λ = E₁, respectively. We chose a variety of parameter settings for µ and q in order to start out from four different shapes of the initial potential.

4.2. A generalized trigonometric Pöschl-Teller system

We will now generate another quantum model from the spheroidal Equation (8). Since the process of construction follows similar steps as in the previous section, we will not give details for all our calculations. This time we employ the following point transformation

$$\chi [y(x)]=\sqrt{\frac{1}{\cos \left(\sqrt{\frac{\mu }{2}x}\right)}}\mathrm{\Psi }(x),y(x)=-\sin \left(\sqrt{\frac{\mu }{2}x}\right),$$

(35)

while we keep the parameter settings (21). The transformation (35) converts the boundary-value problem (8), (9) into the new form

$${\mathrm{\Psi }}^{″}+(E-V_0)\mathrm{\Psi }=0,x\in \left(-\frac{\pi }{\sqrt{2\pi }},\frac{\pi }{\sqrt{2\pi }}\right)$$

(36)

$$\mathrm{\Psi }\left(-\frac{\pi }{\sqrt{2\pi }}\right)=\mathrm{\Psi }\left(\frac{\pi }{\sqrt{2\pi }}\right)=0$$

(37)

where the potential V₀ is given explicitly by

$$V_0=\frac{1-\mu }{2\mu {\cos }^2\left(\sqrt{\frac{\mu }{2}x}\right)}-\frac{q}{4\mu }\cos \left(\sqrt{2\mu x}\right)+\frac{q}{4\mu }-\frac{1}{2\mu }+\frac{1}{2}$$

(38)

We observe that in contrast to the prior system (22), (23), the present governing Equation (36) has standard Schrödinger form. By combining the point transformation (35) with the series (10) and (11), we obtain solutions Ψ_j, j a nonnegative integer, of bound-state type to our boundary-value problem (36), (37) in the form

$${\mathrm{\Psi }}_{2j}={\cos }^{\frac{1}{2}+\frac{1}{2\sqrt{\mu }}\sqrt{\mu +\frac{4}{\mu }-4}}\left(\sqrt{\frac{\mu }{2}z}\right){\displaystyle \sum _{n=0}^{\infty }{(d_{2n})}_{|E=E_{2j}}}{T}_{2n}^{\frac{1}{2\sqrt{\mu }}\sqrt{\mu +\frac{4}{\mu }-4}}\left[-\sin \left(\sqrt{\frac{\mu }{2}z}\right)\right]$$

(39)

$${\mathrm{\Psi }}_{2j+1}={\cos }^{\frac{1}{2}+\frac{1}{2\sqrt{\mu }}\sqrt{\mu +\frac{4}{\mu }-4}}\left(\sqrt{\frac{\mu }{2}z}\right){\displaystyle \sum _{n=0}^{\infty }{(d_{2n+1})}_{|E=E_{2j+1}}}{T}_{2n+1}^{\frac{1}{2\sqrt{\mu }}\sqrt{\mu +\frac{4}{\mu }-4}}\left[-\sin \left(\sqrt{\frac{\mu }{2}z}\right)\right]$$

(40)

Orthogonality and normalizabiliy of these solutions with respect to the Hilbert space $L^2\left(-\pi /\sqrt{2\mu },\pi /\sqrt{2\mu }\right)$ can be established similarly to Equations (27)–(29), such that we omit to give details. In order to perform our confluent SUSY transformations, we need to determine factorization functions. This time we will compute second and third-order transformations, such that we need three transformation functions u₀, u₁ and u₂. Note that the computational effort for determining the Wronskians and their derivatives in Equation (7) is significantly lower then in the previous application, because our governing Equation (36) takes a much simpler form when compared to its counterpart (22). Now, for a nonnegative integer j we pick

$$u_1=\frac{\partial u_0}{\partial {\mathrm{\lambda }}_{|\mathrm{\lambda }=E_j}}{u}_2=\frac{1}{2}\frac{{\partial }^2{u}_0}{\partial {\mathrm{\lambda }}^2{}_{|\mathrm{\lambda }=E_j}},$$

(41)

while u₀ will be taken from the solution set given by Equations (39) and (40), pertaining to the spectral value E_j. Figures 5–8 show potentials V₁ that result from confluent SUSY transformations of second and third order, respectively. In order to start out from different shapes of the initial potential, we picked a variety of parameter settings for q and µ. Due to the length of the expressions involved in the calculations of the SUSY partner potentials, we omit to give details.

5. Conclusions

In this work we have constructed confluent SUSY partners of two quantum systems that are related to the spheroidal equation by means of a point transformation. Although our models do not admit closed-form solutions, we can generate SUSY partners through a series representation for them. Since these representations are only available for bound-state type solutions, we are restricted in the choice of transformation functions and factorization energies. Our approach allows for an extension to more general models by varying the point transformation that connects the spheroidal equation and the targeted system.