{'fraction_non_alphanumeric': 0.04340874096884993, 'fraction_numerical': 0.024458130996091437, 'mean_word_length': 4.322773837667454, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The confluent supersymmetry algorithm for Dirac equations with pseudoscalar potentials

We introduce the confluent version of the quantum-mechanical supersymmetry (SUSY) formalism for the Dirac equation with a pseudoscalar potential. Application of the formalism to spectral problems is discussed, regularity conditions for the transformed potentials are derived, and normalizability of the transformed solutions is established. Our findings extend and complement former results.


Introduction
The formalism of supersymmetry (SUSY) has become popular in quantum mechanics as a method to construct new solvable models. The principal scheme of SUSY goes back to the Darboux transformation that was first introduced in [8] as an algorithm to map solutions of linear differential equations onto each other. As such, it applies to quantum-mechanical systems governed by equations of Schrödinger type. While the concept of Darboux transformations has been developed further and became applicable to a variety of linear and nonlinear equations [10] [19], it has been most successfully used within the SUSY framework in quantum-mechanical applications that involve spectral problems, see e.g. the reviews [7] [12] and references therein. The reason for this lies in the fact that besides generating new systems admitting closed-form solutions, the SUSY scheme can be used for the modification of associated spectra. In particular, discrete spectral values corresponding to bound state solutions of a quantum system, can be created or removed. There is a large amount of literature dedicated to the study of such spectral problems, newer examples include applications to particular systems, involving parametric double-well potentials [27], confined harmonic oscillator intercations [11], or the Swanson model [29]. More examples of such applications can be found from the references in [12]. Another, related branch of research that involves the SUSY scheme concerns the construction of rational extensions to solvable potentials, such that the extended system admits solutions in terms of exceptional orthogonal polynomials [16] [17], see for example the recent studies [22] [26] [25]. Roughly speaking, one can distinguish two different types of SUSY algorithms that apply to models governed by the Schrödinger equation: the standard algorithm that is used in the vast majority of applications, and the less known confluent algorithm. The difference between them manifests in the so-called auxiliary or seed solutions that are necessary to apply the respective algorithm. In the standard case, the latter functions are required to be solutions to the initial Schrödinger equation, while in the confluent case they must solve an inhomogeneous system of equations, often referred to as a Jordan chain. There is considerably less literature on the confluent SUSY scheme than on its standard counterpart. An introduction, discussion of mathematical properties and some applications can be found in [30], [20], [13]- [15], as well as in [4] and [28]. Despite the latter references, properties of the confluent SUSY algorithm remain unknown, for example the link between the auxiliary solutions and regularity of the transformed potential in the general case. Furthermore, generalizations of the standard SUSY scheme to models not governed by Schrödinger equations, do not include the confluent case. An example for such a scenario is the Dirac equation. The SUSY algorithm was introduced [21] and applied [6], but only for the standard context. The purpose of the present work is to complement the latter results by setting up the confluent SUSY scheme for the Dirac equation, applicable to systems governed by pseudoscalar potentials. For applications of such potentials, the reader may refer to [24] [9] and references therein. The remainder of this paper is organized as follows. Section 2 gives a brief review of the confluent SUSY formalism for the Schrödinger equation, while section 3 is devoted to the construction of the confluent algorithm for the Dirac equation. A discussion regarding spectral problems and regular transformed potentials can be found in section 4. The final section 5 presents several examples, demonstrating the generation of solutions to spectral problems as well as the manipulation of the corresponding spectra.

The confluent SUSY scheme
We start out by considering the following stationary Schrödinger equation in atomic units where the energy E is a real-valued constant and U 0 = U 0 (x) denotes the potential. In order to apply a N -th order confluent SUSY transformation (or N -SUSY transformation) to a solution Ψ of (1), we must first determine N ≥ 2 functions u 0 , u 1 , ..., u N −1 , that solve the following system of equations, introducing a real constant λ. The system (2), (3) is commonly referred to as Jordan chain of order N , while the functions u j , j = 0, ..., N − 1 are often called auxiliary solutions or transformation functions. In the conventional SUSY scheme, the auxiliary solutions must satisfy the initial Schrödinger equation at pairwise different energies [3], which precisely constitutes the difference to the confluent algorithm that we are focusing on: here, all auxiliary solutions are associated with the same energy value λ. Now, once the system (2), (3) has been solved, we can construct the following function Φ: where the symbol W stands for the Wronskian of the auxiliary solutions in its index. The function Φ is a solution to the Schrödinger equation for the transformed potential U 1 , given by the expression Here, C stands for an arbitrary constant. While the expressions for the transformed solution (4) and its associated potential (6) look formally the same as in the conventional SUSY scheme [7], they are profoundly different due to the system (2), (3) that determines the auxiliary solutions in the confluent case. By means of the variation-of-constants formula, the following integral representation of the auxiliary solutions can be constructed [28]: whereû stands for any solution of the first equation (2). An alternative representation for the second-order case N = 2 in (2), (3) is given by [5] whereû stands for any solution of (2) and we understand u 0 as a function depending on both the variable x, as well as the parameter λ that enters through equation (2) and it is a sufficiently smooth function of both arguments. It is important to point out that the particular solutions, given by the second terms on the right hand sides of (7) for N = 2 and (8), respectively, in general are not the same. The second representation (8) is very convenient if the integrals in (7) cannot be carried out, as will be demonstrated in the application section 5.2.

The confluent SUSY scheme for the Dirac equation
In this work we focus on the following stationary, two-component Dirac equation in one spatial dimension [21] [23] where σ 2 stands for the second Pauli matrix, E denotes a real-valued number, and the spinor Ψ = (Ψ 1 , Ψ 2 ) T represents the solution. Furthermore, we assume the potential V 0 to be of pseudoscalar form, that is, for a constant, positive mass m, Pauli matrices σ 1 , σ 3 , and a function q 0 . In components, the Dirac equation (9) for the potential (10) can be written as Solving the first of these equations for Ψ 2 gives the relation Substitution of this expression into (12) leads to the equation Hence, any spinor that solves the initial Dirac equation (9), is also a solution to the equations (13), (14) and vice versa. Let us now rewrite (14) as follows where ǫ = E 2 − m 2 and U 0 = q 2 0 + q ′ 0 . Equation (15) is of Schrödinger form and as such admits the SUSY formalism. While the conventional SUSY scheme has already been discussed and applied, we will now focus on the confluent case and make it available to the Dirac equation. To this end, assume that we have found N ≥ 2 solutions u 0 , ..., u N −1 of the Jordan chain (2), (3). Then the function Φ 1 , given by (4) is a solution of the Schrödinger equation The potential U 1 is given by the expression (6), that is, for an arbitrary constant C. Now, in order to convert the transformed Schrödinger potential (17) to a pseudoscalar Dirac potential, we must identify a function q 1 that is the counterpart of q 0 in (10). Such a function q 1 can be obtained once the transformed Schrödinger equation (16) is rewritten in the form Comparison of the two equations (16) and (18) leads to the following condition for the sought function q 1 : where U 1 is given in (17). Equation (19) is of Riccati type, the general solution of which can only be obtained from a known particular solution that will be determined now. First, we apply the linearizing transformation q 1 =q ′ /q to (19), leading to the result After comparison with equation (16) we conclude that a particular solution of (20) is provided by a solution Φ 1 of (16), evaluated at the energy ǫ = 0, that is, As a direct consequence we obtain a particular solution q 1,p of the Riccati equation (19) as It is well-known that the general solution of our Riccati equation can be constructed by means of a particular solution q 1,p via the formula where c denotes an arbitrary constant. Since we know that our particular solution is given in the form (22), we can incorporate the latter function into (23), leading to the simplified version Since the constant c is arbitrary, expression (24) provides a one-parameter family of functions that result from the confluent SUSY transformation. In the last step it remains to convert (24) into a transformed pseudoscalar Dirac potential V 1 , given by The explicit form of this potential is obtained after insertion of (24). We get recall that the functionq is defined in (21). The one-parameter family of pseudoscalar potentials (26) enters in the transformed Dirac equation as follows The first component Φ 1 of the solution spinor Φ = (Φ 1 , Φ 2 ) T to our equation (27) is given in (4), while the second component Φ 2 is determined through a relation analogous to (13): where the explicit form of q 1 can be found in (24). Due to the length of the expressions involved here, we omit to show the expanded form of Φ 2 .

Regularity conditions and spectral problems
While in the previous section we introduced the confluent SUSY algorithm as an entirely computational scheme for the Dirac equation, we will now focus on spectral problems, as they appear in quantum mechanics. In particular, we are interested in finding constraints, under which the transformed potential is free of singularities and the associated transformed solution is L 2 -normalizable. Before we go into more detail, let us assume the Dirac equation (9) to be defined on a real interval D = (a, b), equipped with boundary conditions of Dirichlet type at the endpoints of D, that is, Ψ(a) = Ψ(b) = (0, 0) T , which are allowed to be understood in the sense of a limit, and the norm of a spinor given by We further assume that this spectral problem admits a discrete spectrum (E n ) of energies, where n belongs to some index set, such that ǫ n = E 2 n − m 2 is uniformly bounded from below by its minimum ǫ 0 (ground state energy).
The transformed problem, governed by equation (27) for the boundary conditions Φ(a) = Φ(b) = (0, 0) T , must have a nonsingular potential V 1 in order for the solutions to be physically meaningful. If we do not want singularities in the domain of V 1 , we have to take care of two issues: to obtain a regular Schrödinger potential U 1 from the confluent SUSY transformation and to generate q 1 without additional singularities.
The condition imposed by the confluent algorithm [12] to obtain a Schrödinger potential U 1 without new singularities is that the Wronskian W u 0 ,u 1 ,...,u N−1 in (6) must be nodeless in (a, b). In other words, the Wronskian is allowed to vanish only at the points a or b. This condition has different characteristics for each order N of the SUSY transformation. For example, in the second order case N = 2, the derivative of the Wronskian W u 0 ,u 1 can be expressed, using the Jordan chain (2), (3), as which implies that for a real potential U 0 in (15), the Wronskian is a non increasing monotone function and can be written as where w 0 is an integration constant. If we find a transformation function such that u 0 (a) = 0 or u 0 (b) = 0, then it can be always found a domain for w 0 to avoid zeros in the Wronskian. For the third order case N = 3 in addition of the condition for the value in one of the boundaries of the domain D it is also required that λ ≤ ǫ 0 in (2), for more details of this case see [15]. Now, the regularity of the pseudoscalar Dirac potential (26) also depends on the function q 1 , as given in (22) and (24). Let us consider the particular case (22), our condition is to find a solutionq = Φ 1|ǫ=0 that is free of zeros in (a, b). It generates a nonsingular parametrizing function q 1 if we choose the latter function as the particular solutionq in (22). In case we use the general solution (24) for q 1 , rather than (22), we can avoid singularities in (a, b) by additionally imposing that x 1/q 2 dt is bounded from above or below. This choice implies that there is a constant c in (24), such that the denominator in the second term does not vanish. When the confluent SUSY algorithm is applied to the Schrödinger equation some properties are conserved, in particular the spectral properties which reminds almost equal, the spectral manipulation that can be done are: to add an eigenvalue, to delete one or to obtain a new equation with the same spectrum. The position of this possible extra eigenvalue is in general arbitrary but every order N of the transformation has its own peculiarities. To observe if one extra eigenvalue is inserted or suppressed in the spectrum of the transformed equation special attention have to be paid when ǫ = λ. Notice that u 0 is a solution of (1) and the corresponding transformed solution (4) is the trivial. A linear independent solution v 0 can be obtained asking a direct substitution demonstrate that v 0 is solution of (2). Now, if we transform this solution using the rule (4), then To obtain the last equality first the columns of the Wronskian determinant were interchanged, then all the second and higher derivatives of the transformation functions were replaced by the transformation functions themselves and their first derivative using the Jordan chain (2), (3). With elementary row operations an upper triangular block determinant can be obtained, and since by construction W u 0 ,v 0 = 1, the equation (32) has been demonstrated. If the solution (32) is L 2 −normalizable, then λ is an element of the spectrum of the problem (5) with Dirichlet boundary conditions Φ 1 (a) = Φ 1 (b) = 0. Note again the importance of study the zeros of W u 0 ,...,u N−1 , in this case to know if a element of the spectrum is added or deleted. The confluent SUSY algorithm for the Dirac equation has an extra liberty in the spectral design point of view. The presence of the constant C in (17) has no importance for the Schrödinger equation because it can be interpreted as a shift of the zero energy point, which is arbitrary. If ǫ n is an eigenvalue of (15) then ǫ n + C will be of (16), in the situation that it is not suppressed by the confluent SUSY algorithm. As a result The corresponding eigenvalue of the transformed Dirac equation given in (27) will be |E n | = √ ǫ n + C + m 2 . The constant C will shift every element of the spectrum in a non-linear fashion, and also the gap between the positive and negative elements of the spectrum can be modulated.

Applications
We will now present three applications of the confluent SUSY algorithm for the Dirac equation, involving potentials of physical interest. In the first application we study a Coulomb-like potential, defined on the positive semiaxis. The second application is devoted to a linear potential plus a rational extension that is defined on the whole real line and features a link to the harmonic oscillator. In the last application we consider a trigonometric interaction on a bounded interval.

Coulomb-type system
Let us consider the following boundary-value problem of Dirichlet type for the Dirac equation where the pseudoscalar potential V 0 is given by for a natural number ℓ. Comparison of (35) with the general form (10) shows that the parametrizing function q 0 of the Dirac potential reads The boundary-value problem (33), (34) admits a discrete spectrum of energy values E n,ℓ , where n is a nonnegative integer: The corresponding solutions Ψ = (Ψ 1 , Ψ 2 ) T of our problem (33), (34) have the following explicit form where L stands for a generalized Laguerre function [1] and E n,ℓ is defined in (37). Before we continue with the construction of our SUSY transformation, let us briefly recall that the first component (38) of the Dirac solution also solves the Schrödinger equation (15) for the potential U 0 , which in the present case is given by This is a shifted Coulomb potential with centrifugal barrier, which is well-known to admit a discrete spectrum, confirming (37). We will now apply the confluent 2-SUSY algorithm to the boundary-value problem (33), (34), such that two values are removed from the discrete spectrum of the transformed problem. More precisely, we are aiming at removing a pair of spectral values ±|E n 0 ,ℓ |, where n 0 is a nonnegative integer. In order to perform our transformation, we need two transformation functions u 0 and u 1 , together with an associated constant λ, that solve the Jordan chain (2), (3) for N = 2. We choose the first of these functions as u 0 = Ψ 1|n=n 0 , which solves (2) for the constant λ = E 2 n 0 ,ℓ − m 2 , see (37). We will comment below on the consequences of choosing λ as we did. Now, the second transformation function u 1 we generate by means of the integral representation (7) forû = 0, that is, This integral can be evaluated in closed form, but due to the length of the resulting expression we omit to show the explicit form of u 1 here. Now that we have the two required transformation functions, let us use the conditions for regularity and normalizability of the transformed potential and solution, respectively. These conditions, developed in section 4, state that the Wronskian (30) of the transformation functions must be nonzero inside the domain D of our boundary-value problem, given by D = (0, ∞). According to (30), we have the following Wronskian Since the integrand is positive, its integral is a strictly increasing function, which due to its continuity can have at most one zero. If we choose the integration constant w 0 as then we know that (41) has exactly one zero, which is located at x = 0. We are therefore guaranteed that the function (17) is free of singularities inside its domain D = (0, ∞). Next, let us comment on the choice of λ = E 2 n 0 ,ℓ − m 2 , which determines integrability of the transformed solution Φ 1 to equation (16). Since by the election λ matches one of the spectral values admitted by the boundary-value problem (5), for the potential (39) and Ψ 1 (0) = Ψ 1 (∞) = 0, and also because the Wronskian W u 0 ,u 1 vanishes at the origin, the corresponding solution (32) will not be L 2 −nomalizable and as a consequence the element will be removed from the transformed Schrödinger problem with the same boundary condition and using C = 0. At the Dirac level, this corresponds to the pair of numbers ±|E n 0 ,ℓ | that is missing in the discrete spectrum of the transformed problem. In order to set up the potential associated with the latter problem, we must construct the parametrizing function q 1 , using either (22) or (24). For the sake of simplicity, we will take the first option, that is, Unfortunately, the explicit form of this function is very long, such that we do not show it here. The same is true for the transformed solution Φ = (Φ 1 , Φ 2 ) T , which is constructed according to (4) and (28). Recall that the transformation functions are given by u 0 = Ψ 1|n=n 0 and (40), where q 1 is determined by (43). The functions (44), (45) provide a solution to the transformed boundary-value problem, given by for a pseudoscalar potential V 1 that is obtained by inserting the parametrizing function (43) into the general form (25): The discrete spectrum of the transformed boundary-value problem (46), (47) is given by (37) except for the two values that correspond to n = n 0 . Before we complete this example, let us visualize special cases of the initial and transformed parametrizing functions q 0 and q 1 , as given in (36) and (43), respectively, in figure 1. The figure reveals that the functions q 1 diverge stronger at zero than their counterparts q 0 , which can be understood by analyzing the transformation and by determining the behaviour of the transformed parametrizing function q 1 close to x = 0. To this end, let us recall that the transformation function u 0 is defined by u 0 = Ψ 1|n=n 0 , where Ψ 1 can be found in (38). Since close to x = 0 we have u 0 ∝ x ℓ+1 , the Wronskian (41), (42) shows the following asymptotic behaviour Let us express this Wronskian in the convenient way where the function W catches the behaviour of the Wronskian away from x = 0, and as such is bounded there. Substitution of the expression (48) into the general form (17) for the transformed potential function U 1 , the latter function renders in the form Since we would like to determine the behavior of the transformed function q 1 in (43), it is important to study the possible behavior of the solutions of (16) for the potential (49). This equation allows two types of solutions according to their behavior close to the origin: the first vanishes in the origin as Φ 1|n=−1 ∝ x ℓ+3 and the second diverges as Φ 1|n=−1 ∝ x −(ℓ+2) , in this example we used the latter and as a consequence Φ 1|n=−1 can then be written as where the function ϕ is bounded at x = 0. The latter expression for Φ 1|n=−1 can now be plugged into (43), yielding the constraint note that the term ϕ ′ /ϕ is bounded at x = 0. On comparing (50) to (36), we observe that the coefficient of the −x −1 term is greater for q 1 , which explains why the transformed function q 1 diverges faster than q 0 at x = 0. Figure 1: the functions q 0 (dotted curve) and q 1 (solid curve). In order to generate q 1 , we removed the spectral values ±|E 2,1 | in the left figure, and ±|E 2,2 | in the right figure.

Harmonic oscillator system
In contrast to the previous example, this time we will start out from a boundary-value problem for the Schrödinger equation (14), and generate the Dirac SUSY partners arising from the confluent algorithm. Let us consider the following problem where the real constant A is a free parameter that will prove useful below. By inspection of (14) we observe that in the present case which is a shifted harmonic oscillator interaction. The boundary-value problem (51), (52) admits an infinite discrete spectrum (ǫ n ), where n is a nonnegative integer: The associated, L 2 -normalizable solutions are given by where H denotes a Hermite polynomial [1]. Before we continue, let us state the general solution u of our Schrödinger equation (51), which will prove useful when we construct the Dirac problem.
for two free constants c 1 and c 2 . Furthermore, the symbol 1 F 1 stands for the confluent hypergeometric function [1]. We will now set up a boundary-value problem for the Dirac equation (9) that can be related to (51), (52). To this end, we determine the parametrizing function q 0 of the pseudoscalar Dirac potential to be constructed. According to U 0 = q 2 0 + q ′ 0 and (53), we must solve the equation , with respect to q 0 . Since this equation has the same form as its counterpart (19), we can use the same solution method. To keep calculations simple, we will restrict ourselves to the particular solution inspired by (22). We therefore need to evaluate a solution of the Schrödinger equation (51) at zero energy, that is, ǫ = 0. After setting ǫ = 0 in the general solution (56) we obtain While in general this expression consists of intricate special functions, it can be made rational and real-valued by appropriate choices of the constants c 1 , c 2 and A. In particular, one can obtain rational functions q 0 if the constant A is chosen as A = −4k − 1 for a nonnegative integer k and if c 2 = 0, as the following examples show: Let us pick the case (57) for our parametrizing function q 0 . This choice determines a pseudoscalar potential (10) of the form which enters in the corresponding boundary-value problem for the Dirac equation that reads The discrete spectrum admitted by this problem can be found by plugging (54) into the equation ǫ = E 2 − m 2 and solving for E: The first component Ψ 1 of the corresponding solution Ψ = (Ψ 1 , Ψ 2 ) to this problem is given in (55), while the second component Ψ 2 can be generated by means of (13), that is, We are now ready to apply a confluent 2-SUSY transformation to the Schrödinger problem (51), (52). To this end, we need two transformation functions u 0 , u 1 that solve the system (2), (3). Let us choose the function u 0 as a special case of the general solution (56) for the settings A particular effect of this choice for c 1 and c 2 is that the function u 0 vanishes at negative infinity, while it becomes unbounded at positive infinity. This behavior is one of the condition needed to obtain a Wronskian free of zeros. It remains to determine a second transformation function u 1 in order to solve the system (2), (3). In the present case the integral representation (7) for j = 1 does not yield the sought function u 1 in closed form, because the integrals cannot be evaluated, given that u 0 is given by (56) and (59) for a general value of λ. Therefore, in order to avoid integration, we will resort to the differential formula (8), which can be written in the form for an arbitrary real constant B and a function v 0 defined by (31). Remind that u 0 and v 0 are linearly independent solutions of equation (2) with Wronskian W u 0 ,v 0 = 1. The useful expression for the derivative with respect the first parameter of the hypergeometric function can be found in [2]. Let us point out that we will not have to resolve the integral that was introduced in (60) due to our knowledge about the Wronskian of u 0 and v 0 . Now that we have determined the latter two transformation functions through (56), (59), and (60), we can perform a confluent 2-SUSY transformation on the boundary-value problem (51), (52), recall that we chose A = −5.
According to (17) for N = 2, the potential U 1 in the transformed Schrödinger equation (16) is given by where the constant C was set to zero. Observe that the Wronskian W u 0 ,u 1 in the logarithm's argument can be simplified by substituting the explicit form (60) of u 1 and then expanding the determinant [5]. As a result, the integral contained in (60) does not appear in the Wronskian anymore. For the sake of brevity we do not show the explicit form of the transformed potential (61). The transformed solution Φ 1 of the equation (16) is constructed by means of (4) for N = 2, that is, Note that the Wronskian in the numerator was obtained by Laplace expansion, see [5] for details. As in the case of (61), the integral in (60) does not appear, such that we do not need to resolve it. The L 2 -normalizable functions (62) satisfy the spectral problem governed by (16) and (61), subject to the boundary conditions implying that Φ 1 must vanish at the infinities. We omit to substitute the explicit forms of the functions involved in (62) due to their length. Since we chose the parameters (59), the discrete spectrum of the transformed problem is given by (54), and contains an additional spectral value λ. The associated L 2 -normalizable solution can be found substituting u 0 with the parameters (59) and (60) in (32). This yields where the Wronskian in the numerator was simplified similarly to its counterpart in (62). Again, for reasons of brevity, we do not include the full form of (63). In figure 2 we show a plot of typical potential functions U 0 and U 1 , as given in (53) for A = −5 and (61), respectively. We are now ready to generate a Dirac boundary-value problem from the transformed Schrödinger equation (16) with potential (61). To this end, we must generate the parametrizing function q 1 that will enter in the pseudoscalar potential of the Dirac equation. According to (22), in the simplest case we can construct q 1 through a solution of (16), (61) that is taken at ǫ = 0. If this solution does not have any nodes, then the function q 1 will be regular on the whole real line. Let us construct q 1 as follows, using our transformation functions u 0 and u 1 as in (62): where u is the general solution (56) that is taken here at ǫ = c 2 = 0. Now that we have determined the parametrizing function q 1 for the pseudoscalar potential of the transformed Dirac equation, it remains to construct the second component Φ 2 of its solution Φ = (Φ 1 , Φ 2 ) T . According to (28), we have In this expression q 1 must be substituted from (64), while the function Φ 1 can be either taken from the solution set (62) or it is given by the single solution (63). The constant E are elements of the discrete spectrum that is determined from (58) as follows (E n ) = ± 2 n + 6 + m 2 : n = 0, 1, 2, 3...
where the last contribution contains the newly created spectral value λ associated with (63). The functions (62), (65) solve the transformed boundary-value problem the pseudoscalar potential of which is given by (26), that is,  observe that the constant C enters in all discrete spectral values.

A third-order confluent SUSY transformation
While the previous examples referred to the confluent SUSY algorithm of second order, let us now present a transformation of third order. To this end, we consider the folllowing boundaryvalue problem of Dirichlet type where the pseudoscalar potential V 0 is given by Comparison of (70) with the general form (10) shows that the parametrizing function q 0 reads q 0 = 2 tan(x) + 2 tan(x) .
The boundary-value problem (68)-(70) admits an infinite discrete spectrum (E n ) and a corresponding orthogonal set of solutions (Ψ n ), where n is a nonnegative integer. The spectral values are given by while the associated solution spinor components Ψ 1 , Ψ 2 can be expressed through hypergeometric functions as follows Ψ 1 = cos 3 (x) sin 2 (x) 2 F 1 −n, n + 5, 5 2 , sin 2 (x) (73) Ψ 2 = 4 n (n + 5) cos 4 (x) sin 3 (x) 5 m + 5 (2 n + 5) 2 + m 2 2 F 1 −n + 1, n + 6, m + (2 n + 5) 2 + m 2 2 F 1 −n, n + 5, Observe that the first argument of the hypergeometric functions is a nonpositive integer, such that they can be represented through Jacobi polynomials. For the sake of brevity, we omit to state the latter representation explicitly, and instead refer the reader to [1]. Figure 4 shows normalized probability densities associated with the solutions (73), (74) for the first three values of n. Now, in order to perform a third-order confluent SUSY transformation, we must find three auxiliary solutions u 0 , u 1 , u 2 that solve the Jordan chain (2), (3) for N = 3. As a first step, we choose the function u 0 as follows: This choice determines the energy λ in (2), (3) as λ = 25, which we obtained from (72) for n = 0. Next, we use the integral representation (7) to find the remaining two auxiliary solutions u 1 and u 2 . Substitution of (75) yields the result note that we choseû = 0. Furthermore, we set all integration constants to zero except for the inner integral in u 2 , where we picked the integration constant as −1/20. Now that the three auxiliary solutions have been determined, we are ready to perform the confluent SUSY transformation. To this end, we plug (73), (75)-(77) into (4): This function constitutes the first component of our transformed solution spinor Φ = (Φ 1 , Φ 2 ) T . Since the resulting expression for Φ 1 is very long, we omit to state it here in explicit form. Before we can determine the second component, the parametrizing function q 1 of the transformed potential in (27) must be constructed using (22) or (24). We choose the simplest case (22), given by Observe that choosing C = 0 in (17) and n = −5/2 results in ǫ = E 2 − m 2 = 0, that is, the function Φ 1 is evaluated at zero energy ǫ. Similar to Φ 1 , the function q 1 has a very long and involved form, such that we do not display it here. The same holds for the remaining component Φ 2 of the transformed solution spinor Φ that is generated by substituting (78) and (79) into (28). Figure 5 shows the parametrizing functions q 0 and q 1 of our initial and transformed potentials, respectively. Furthermore, normalized probability densities for our transformed solution spinor are visualized in figure 6 for the first three values of n. Note that there is no solution for n = 0 anymore, as the corresponding bound state was was deleted by our choice of u 0 in (75).

Concluding remarks
We have introduced the confluent SUSY formalism for the Dirac equation with a pseudoscalar potential. The present discussion can be easily extended to scalar potentials, as this requires only a reparametrization [21]. The formalism was applied to three different interactions each one with different domain of definition, the first a Coulomb-like potential, then to a linear potential with a rational extension and finally to a trigonometric interaction. Some of the possible spectral manipulation that the formalism permits were illustrated. While in principle confluent transformations of arbitrary order can be performed using the formalism presented in this work, the calculations become very complicated, even when resorting to a symbolic calculator.