Gabriel Cirio

Maud Marchal

Aurélien Le Gentil

Anatole Lécuyer

Élien Aur

Le Gentil

Anatole L Écuyer

Tap, squeeze and stir the virtual world: Touching the different states of matter through 6DoF haptic interaction

Keywords: I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling-Physically based modeling, H.5.1 [Information Interfaces and Presentation]: User Interfaces-Haptics I/O, H.5.1 [Information Interfaces and Presentation]: Multimedia Information Systems-Artificial, Augmented, and Virtual Realities

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

The computer graphics community has developed many different physically-based models to recreate materials with which we interact on a daily basis: smoke, fire, water, clay, plastic and metal, to name a few. We are used to see these simulations in movies or video games, but what about touching them? Haptic rendering algorithms compute force feedback in real-time, and allow to feel the objects in a virtual world through a haptic device. Due to the high complexity of physical simulations and the high update rate required for haptics, these algorithms often focus on enabling the haptic interaction with a single type of matter: either rigid, deformable or fluid media. However, even the most common activity in real life often involves interacting with more than one type of media: walking on pavement, parquet and puddles, having a meal on a plate with meat and a beverage, washing our hands by touching the faucet, the soap and the water. In order to bodily interact with such virtual worlds, we currently need to use different algorithms for each type of media, and set up coupling interfaces between these heterogeneous models. This supposes to increase the complexity and the computational cost of an already highly complex and time consuming simulation.

In this work, we propose a single model with a single haptic rendering technique to simulate and interact simultaneously with rigid, deformable and fluid media in real-time with 6DoF haptic feedback. This novel approach is the first to tackle the problem of physicallybased multistate haptic rendering. Based on a Smoothed-Particle Hydrodynamics (SPH) model [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF] brought to high update rates * e-mail: {gcirio, mmarchal, alegentil, alecuyer}@inria.fr through a GPU implementation, we are able to simulate fluid volumes, deformable and rigid bodies in real-time. Using a novel coupling scheme that takes advantage of the SPH model, we show how to render interaction forces and torques through a 6DoF haptic device. We measure the performance of our technique in terms of computation time and haptic feedback.

PREVIOUS WORK

Previous work on haptic interaction with deformable bodies [START_REF] Burdea | Virtual reality training for the diagnosis of prostate cancer[END_REF] or fluids [START_REF] Cirio | Six Degrees-of-Freedom Haptic Interaction with Fluids[END_REF] only allow the interaction with rigid bodies on top of another state. The CORDIS-ANIMA [START_REF] Cadoz | CORDIS-ANIMA: Modeling and Simulation System for Sound and image Synthesisthe General Formalism[END_REF] simulation framework was the first to simulate different states with haptic feedback, although using a mass-spring formalism and not the physically-based approach of continuum mechanics. To the best of our knowledge, there is no haptic rendering technique that provides haptic feedback for fluid, deformable and rigid states of matter in the same simulation. Recent work on physically-based simulations has used the SPH model [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF] for the simulation of multiple states of matter in a unified framework. However, update rates are usually not real-time, or fall short for haptic interaction. In [START_REF] Muller | Interaction of fluids with deformable solids[END_REF], deformable bodies are modeled in real-time with SPH and the Moving Least Squares (MLS) algorithm to compute the elastic forces. The approach is improved in [START_REF] Keiser | A unified lagrangian approach to solid-fluid animation[END_REF], allowing solid, deformable and fluid animation and interaction in a unified approach. A fully SPH-based approach is presented in [START_REF] Solenthaler | A unified particle model for fluid-solid interactions[END_REF], later corrected in [START_REF] Becker | Corotated SPH for deformable solids[END_REF] for a rotationally invariant formulation. Although these approaches allow a certain degree of interaction, they are not fast enough and do not solve the haptic rendering issue.

Our approach is the first to achieve multistate haptic rendering. Previous haptic rendering techniques were focused on a single state, e.g. fluid, deformable or solid, making it difficult to have different types of media coexist in the same simulation while providing convincing haptic feedback. Inspired by our previous work presented in [START_REF] Cirio | Six Degrees-of-Freedom Haptic Interaction with Fluids[END_REF], and bringing [START_REF] Solenthaler | A unified particle model for fluid-solid interactions[END_REF] and [START_REF] Becker | Corotated SPH for deformable solids[END_REF] for deformable bodies to the required speeds, we are able to simulate simultaneously fluid volumes, deformable bodies and rigid objects in a physically-based manner and at high update rates. Through the SPH model, we enable the haptic coupling with 6DoF haptic devices and the force feedback interaction with matter in its differrent states.

PHYSICALLY-BASED SIMULATION FOR HAPTIC RENDER-ING OF MULTISTATE MEDIA

The seamless cohabitation of different theoretical models for different states of matter within the same virtual world is made possible through the use of a single common computational approach: the SPH model. With all three states of matter simulated with the same data structures and the same algorithmic, the coupling between the different media becomes a straightforward task, and is included in the simulation. Figure 1 illustrates this principle. Any virtual object of the virtual world can be either fluid, deformable or rigid, and the state is hidden inside the model. The interactions between the different states are also computed seamlessly based on the same model. In the end, a virtual world is populated by defining the state of each virtual object together with its corresponding physical properties (stiffness for deformable bodies, viscosity for fluids). Users can interact with the virtual world through a haptic device using the interaction forces, without ever dealing with the models behind. In this section, we describe the foundations of our physicallybased simulation algorithms.

SPH-based Physical Simulation

The SPH model [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF] computes the smoothed quantity A i of a particle i at any position x in space through the general formula:

A i = ∑ j A j V j W (x -x j , h) (1)
where A j is the discrete quantity A sampled for neighboring particle j at position x j , V j is the volume of j, and W is the smoothing kernel of support h, where particles beyond h are not taken into account.

Fluid and Rigid Body Model

The motion of fluids is driven by the Navier-Stokes equations. Using the SPH model as presented in [START_REF] Muller | Particle-based fluid simulation for interactive applications[END_REF] and used in [START_REF] Cirio | Six Degrees-of-Freedom Haptic Interaction with Fluids[END_REF], pressure and viscosity forces are computed at each time step:

f pressure i = -V i ∑ j V j P i + P j 2 ∇W (x i -x j , h) (2)
f viscosity i = µV i ∑ j V j (v j -v i )∇ 2 W (x i -x j , h) (3)
where P is the pressure, v the velocity, µ the viscosity coefficient and g the gravity field. ∇ and ∇ 2 are respectively the gradient and Laplacian of the physical quantities. Parameters such as the viscosity can be changed to obtain different fluid behaviors, from smoke (no viscosity) to honey (very high viscosity).

Rigid bodies are modeled with the same SPH particles used in the fluid simulation, creating a unified particle model as described in [START_REF] Solenthaler | A unified particle model for fluid-solid interactions[END_REF]. This simple rigid body model has the main advantage of being meshless: since the overall shape of the rigid body is not important, it allows the seamless use of arbitrary-shaped rigid bodies, including concave objects. Moreover, it provides a gain in efficiency by computing a single model, since it makes the approach unified with fluid and deformable bodies. Additional collision detection techniques are avoided, since collision computations are included in the neighbor search.

Deformable Body Model

In order to simulate a deformable body with SPH governed by the elasticity theory, we compute elastic forces as in [START_REF] Solenthaler | A unified particle model for fluid-solid interactions[END_REF] using the improvements proposed in [START_REF] Becker | Corotated SPH for deformable solids[END_REF] to make the computations rotationally invariant by using a corotational approach. For more details, we refer the reader to the respective articles, as well as to [START_REF] Muller | Interaction of fluids with deformable solids[END_REF].

The elastic force exerted on a particle i is computed from the gradient of the locally rotated displacement field, ∇ ūi . Using an SPH formulation, this gradient is defined as:

∇ ūi = ∑ j V 0 i ū ji ∇W (x 0 i -x 0 j , h) T (4)
where x 0 i and x 0 j are the initial positions (in the undeformed state) of particles i and j respectively, and ū ji is the locally rotated displacement difference between neighboring particles j and i. Then, the strain ε can be computed from ∇ ūi , using a Green-Saint-Venant strain tensor formulation.

The elastic force f ji exerted on a neighboring particle j by particle i is defined as the negative gradient of strain energy with respect to displacement, and can be computed as:

f ji = -2V 0 i (I + ∇ ūT i )σ i V 0 j ∇W (x 0 i -x 0 j , h) (5)
where V 0 i is the initial volume of particle i and σ is the stress. σ is computed from the strain ε through the linear relation σ = Cε, where C is defined by only two parameters, the Young modulus and the Poisson ratio.

The elastic force f elasticity i exerted on particle i is symmetrized and computed as:

f elasticity i = ∑ j -R i f ji + R j f i j 2 ( 6
)
where R is the rotation matrix computed through the corotational approach.

Interaction Forces

In a unified, parallel and time-critical framework, unified interaction forces fit well for computation time reasons. Using the same interaction forces as in a fluid-fluid case, namely pressure and viscosity forces, improves the gain of parallel computation, while providing a reasonable amount of control over the forces through density and viscosity values. Hence, as in [START_REF] Cirio | Six Degrees-of-Freedom Haptic Interaction with Fluids[END_REF], interaction forces f interaction i are computed between particles belonging to different bodies as:

f interaction i = f pressure i + f viscosity i (7)
Moreover, by using the SPH model for the computation of interaction forces between fluid, deformable and rigid body particles, it removes the need of additional collision detection algorithms.

Integration and Simulation Loop

The different models previously described are run simultaneously under the same SPH framework, improving the efficiency of the simulation compared to having three heterogeneous models. Hence, the simulation loop runs in four consecutive steps:

1. Compute corotational matrices R i for each particle belonging to a deformable body; 2. Compute new density and elastic properties; 3. Compute new internal forces (between particles of the same body) and new interaction forces (between particles of different bodies); 4. Integrate (compute new velocity and position) through Leap-Frog integration for fluids and deformable bodies. For rigid bodies, integrate by applying rigid body dynamics.

The entire physical model is implemented on GPU using the CUDA framework, extending our previous implementation [START_REF] Cirio | Six Degrees-of-Freedom Haptic Interaction with Fluids[END_REF] to support deformable bodies. Space is subdivided in a 3D uniform grid, so that search for neighbors is limited to the surrounding cells. We compute the different physical quantities and forces related to fluids and elastic theories, as well as interaction forces. In order to further accelerate computations, we use texture lookups to retrieve the different data arrays to benefit from texture data caching. We also use an OpenGL Vertex Array Object for the particle position array, allowing the graphic rendering algorithm to access position data without copying it back to CPU memory.

MULTIMODAL RENDERING 4.1 6DoF Haptic Rendering

We base our force feedback computation on interaction forces and on an SPH rigid body as proxy. By using an SPH proxy, we avoid additional computations for haptic coupling: the proxy is an object of the virtual world, and it interacts with all other media around it through the unified interaction forces described in section 3.1.3. When the proxy interacts with different media on the virtual world, the simulation computes the corresponding interaction forces as any other body from the virtual world. In this way, the user can feel the virtual world through an arbitrary-shaped rigid body.

Once the interaction forces are computed on the proxy, they are fed to the rigid body dynamics algorithm to perform the haptic coupling between the proxy and the 6DoF haptic device. In the rigid body dynamics step, all the interaction forces exerted on the proxy rigid body body are summed to obtain a total force F haptic body and a total torque τ haptic body . The force and torque feedbacks coming from the haptic device are then added as external forces. Hence:

F haptic body = ∑ i body f interaction i body + f device (8) τ haptic body = ∑ i body τ interaction i body + τ device (9)
Then, the new position and velocity are computed by integrating the forces over the simulation time step. They are sent to the haptic device, closing the haptic loop in admittance mode. A Virtual Coupling mechanism [START_REF] Colgate | Issues in the haptic display of tool use[END_REF] is introduced between the haptic device and the rigid body, reconciling a high update rate haptic device with a lower rate simulation, leading to an increase of stability.

The resulting 6DoF haptic coupling scheme is unified, in the sense that it allows the interaction with different media (fluid, deformable and solid) without distinction, since the haptic forces are computed in a unified way.

Graphic Rendering

As for the physical simulation, we use a unified visual rendering algorithm for fluid and deformable bodies. Previous work [START_REF] Cirio | Six Degrees-of-Freedom Haptic Interaction with Fluids[END_REF] has shown how to render fluids with a screen-space technique based on [START_REF] Van Der Laan | Screen space fluid rendering with curvature flow[END_REF], with a trade-off between quality and performance. A threestep approach first computes per pixel fluid data (front and back fluid volume depth from view), then smoothes the front surface depth in screen space using a fast bilateral filter, and finally composes the frame.

We follow the same procedure to render deformable bodies. However, since most bodies are not transparent, as opposed to fluids, we do not handle transparency for deformable bodies, which leads to an increase in performance. In order to achieve a correct occlusion between the different media, we perform a depth test in the fragment shader.

PERFORMANCE

In this section we measure the performance of our haptic rendering technique, in terms of computation time and haptic feedback.

The simulations were carried out using a Virtuose 6DoF forcefeedback device from Haption Company, and a laptop computer with a Core 2 Extreme X7900 processor at 2.8GHz, 4GB of RAM, and a Nvidia Quadro FX 3600M GPU with 512MB of memory.

Computation Time

We measured the computation time of three scenes with an increasing number of particles. Each scene has either a fluid volume, a deformable body or a rigid body, all of the same size and shape (a cuboid with a square base of 0.5m side and a variable height according to the number of particles). We only measured the physical simulation computation time, without taking into account graphic rendering. The simulation ranges from 2,000 to 30,000 particles. Figure 2 plots the measurements for the different types of media. The higher complexity of the elastic theory implementation compared to the fluid theory and the rigid body dynamics is clearly visible: when the number of particles increases, the computation time for deformable bodies increases roughly 2 times faster than for fluids, and 6 times faster than for rigid bodies.

A determinant factor of our approach is its computation time, since haptic feedback requires update rates higher than the usual 30Hz of real-time visual rendering. It has been shown that, for fluids, update rates of 70Hz are satisfactory for haptic rendering, since forces due to fluids do not change rapidly [START_REF] Cirio | Six Degrees-of-Freedom Haptic Interaction with Fluids[END_REF]. Our implementation allows to simulate large volumes of fluid (around 20,000 particles) above 70Hz, as shown in Figure 2. In the case of deformable bodies, low update rate simulations (30-50Hz) use interpolation [START_REF] Saupin | Contact model for haptic medical simulations[END_REF], and faster models rely on pre-computed data and small deformation approaches [START_REF] James | A unified treatment of elastostatic contact simulation for real time haptics[END_REF], or limit the complexity of the virtual objects [START_REF] Burdea | Virtual reality training for the diagnosis of prostate cancer[END_REF]. Our model is an efficient trade-off between these approaches, with a fast simulation with complex deformable bodies (150Hz for 5,000 particles) and a virtual coupling mechanism for a smooth haptic rendering.

Haptic Feedback

We provide a visual feedback of the forces exerted on a coupled rigid body interacting with the three types of media, in order to show the correspondence of the simulation behavior with the expected physical behavior.

A 20cm radius sphere, coupled to a 6DoF Virtuose, is dropped onto a 50x50x20cm block of particles, representing a volume of fluid, a deformable body or a rigid body. Figure 3 shows the forces involved in the interaction between the coupled sphere and the media. The three interaction patterns are computed and recorded individually, but are shown on the same graph for comparison purposes.

When the media is solid, the sphere bounces quickly, generating a force peak. When the media is deformable, the sphere compresses the media's surface until its elastic behavior pushes the sphere back into the air. When the media is fluid, the sphere penetrates the fluid until pressure makes the sphere rise to the surface and float.

APPLICATION SCENARIOS

Our multistate haptic rendering technique can be applied to different fields, with haptic feedback enhancing the interaction experience. The medical field can substantially benefit from our approach, since surgeons are in constant interaction with fluids (blood), deformable bodies (organs, tissues) and rigid bodies (bones). We designed a medical simulation scene as an illustration of potential application, as shown in Figure 4, with these three states. The user can interact with the different body parts through a rigid probe, and experience haptic feedback through a 6DoF haptic device. Where interactive medical simulations of this kind were restricted to deformable bodies only, we enable haptic feedback on virtual worlds with complex rigid and deformable objects, together with volumes of fluid, without dealing with different models and their coupling.

Other potential application areas include the entertainment field (exploration of virtual worlds), training (manipulations of dangerous materials) and industrial simulations (virtual prototyping).

CONCLUSION

This paper introduces the first multistate haptic rendering approach, allowing haptic interaction simultaneously with media in fluid, deformable and rigid states. Based on the SPH physical model for all three types of media, our method avoids the complexity of dealing with different algorithms and their coupling. We achieve high update rates while simulating a physically-based virtual world governed by continuum mechanics theories, and show the multimodal rendering of the interaction, with forces and torques displayed through a 6DoF haptic device and with an optimized visual rendering algorithm. Preliminary tests with novice users were promising, with users giving positive feedback. Future work will focus on decoupling graphic and haptic rendering loops for an increased simulation performance. We will then conduct a perceptual evaluation in order to assess the capability of users to recognize the different states of matter.

Figure 1 :

 1 Figure1: Global view of our approach: a user seamlessly interacts with a virtual world populated with objects of different states of matter.

Figure 2 :

 2 Figure 2: Computation time for scenes made of either a rigid body (dashed curve), a deformable body (continuous curve) or a fluid volume (dotted curve), according to the number of particles.

Figure 3 :

 3 Figure 3: Forces exerted on a sphere successively dropped on rigid, deformable and fluid media.

Figure 4 :

 4 Figure 4: 6DoF haptic interaction in a medical scenario. Left: rendering of the particles composing the multistate virtual world. Middle: the user examines a deformable muscle. Right: fluid blood pours from the deformable intestine when the user penetrates it with the rigid probe.

ACKNOWLEDGEMENTS

This work was supported by the European Community under FP7 FET-Open grant agreement n • 222107 NIW -Natural Interactive Walking.


Yves Grandati

Multistep DBT and regular rational extensions of the isotonic oscillator

Keywords:

In some recent articles we developed a new systematic approach to generate solvable rational extensions of primary translationally shape invariant potentials. In this generalized SUSY QM partnership, the DBT are built on the excited states Riccati-Schrödinger (RS) functions regularized via specific discrete symmetries of the considered potential. In the present paper, we prove that this scheme can be extended in a multistep formulation. Applying this scheme to the isotonic oscillator, we obtain new towers of regular rational extensions of this potential which are strictly isospectral to it. We give explicit expressions for their eigenstates which are associated to the recently discovered exceptional Laguerre polynomials and show explicitely that these extensions inherit of the shape invariance properties of the original potential.

I. INTRODUCTION

Since the pionnering work of Gomez-Ullate et al [START_REF] Gómez-Ullate | The Darboux transformation and algebraic deformations of shape invariant potentials[END_REF], the exceptional orthogonal polynomials (EOP) and their connexion with rational extensions of solvable quantum potentials have been an active research subject [START_REF] Gómez-Ullate | Supersymmetry and algebraic Darboux transformations[END_REF][START_REF] Gómez-Ullate | An extended class of orthogonal polynomials defined by a Sturm-Liouville problem[END_REF][START_REF] Gómez-Ullate | An extension of Bochner's problem: exceptional invariant subspaces[END_REF][START_REF] Gómez-Ullate | Exceptional orthogonal polynomials and the Darboux transformation[END_REF][START_REF] Gómez-Ullate | On orthogonal polynomials spanning a non-standard flag[END_REF][START_REF] Quesne | Exceptional orthogonal polynomials, exactly solvable potentials and supersymmetry[END_REF][START_REF] Quesne | Solvable rational potentials and exceptional orthogonal polynomials in supersymmetric quantum mechanics[END_REF][START_REF] Bagchi | Isospectrality of conventional and new extended potentials, second-order supersymmetry and role of PT symmetry[END_REF][START_REF] Bagchi | An update on PT -symmetric complexified Scarf II potential, spectral singularities and some remarks on the rationally-extended supersymmetric partners[END_REF][START_REF] Odake | Infinitely many shape invariant potentials and new orthogonal polynomials[END_REF][START_REF] Odake | Another set of infinitely many exceptional (X l ) Laguerre polynomials[END_REF][START_REF] Ho | Properties of the exceptional (X l ) Laguerre and Jacobi polynomials[END_REF][START_REF] Odake | Infinitely many shape invariant potentials and cubic identities of the Laguerre and Jacobi polynomials[END_REF][START_REF] Sasaki | Exceptional Laguerre and Jacobi polynomials and the corresponding potentials through Darboux-Crum transformations[END_REF][START_REF] Dutta | Conditionally exactly solvable potentials and exceptional orthogonal polynomials[END_REF][START_REF] Ho | Zeros of the exceptional Laguerre and Jacobi polynomials[END_REF][START_REF] Grandati | Solvable rational extension of translationally shape invariant potentials[END_REF][START_REF] Grandati | Solvable rational extensions of the isotonic oscillator[END_REF][START_REF] Grandati | Solvable rational extensions of the Morse and Kepler-Coulomb potentials[END_REF]. The EOP appear to be the constitutive elements of the eigenstates of solvable rational extensions of the second category primary translationally shape-invariant potentials (TSIP) [START_REF] Grandati | Rational solutions for the Riccati-Schrödinger equations associated to translationally shape invariant potentials[END_REF][START_REF] Cooper | Supersymmetry in Quantum Mechanics[END_REF][START_REF] Dutt | Supersymmetry, shape invariance and exactly solvable potentials[END_REF][START_REF] Gendenshtein | Derivation of exact spectra of the Schrodinger equation by means of supersymmetry[END_REF] In a series of recent papers [START_REF] Grandati | Solvable rational extension of translationally shape invariant potentials[END_REF][START_REF] Grandati | Solvable rational extensions of the isotonic oscillator[END_REF][START_REF] Grandati | Solvable rational extensions of the Morse and Kepler-Coulomb potentials[END_REF] we have developped a new approach which allows to generate infinite sets of such regular solvable extensions starting from every TSIP in a very direct and systematic way without taking recourse to any ansatz or ad hoc deforming functions. This approach is based on a generalization of the usual SUSY partnership built from excited states. The corresponding Darboux-Bäcklund Transformations (DBT), which are covariance transformations for the class of Riccati-Schrödinger (RS) equations [START_REF] Grandati | Rational solutions for the Riccati-Schrödinger equations associated to translationally shape invariant potentials[END_REF], are based on regularized RS functions corresponding to unphysical (negative energy) eigenfunctions of the initial hamiltonian. They are obtained by using discrete symmetries acting on the parameters of the considered family of potentials. If the use of negative energy states has been already proposed in early years of SUSY MQ development [START_REF] Sukumar | Supersymmetric quantum mechanics of one-dimensional systems[END_REF], such a systematic scheme has never been envisaged. In the case of the isotonic oscillator, this gives the three infinite sets L1, L2 and L3 of regular rationally solvable extensions of this potential, the L1 and L2 series potentials being strictly isospectral to the isotonic one and inherits of its shape-invariance property. Very recently Gomez-Ullate, Kamran and Milson [START_REF] Gómez-Ullate | Two-step Darboux transformations and exceptional Laguerre polynomials[END_REF] introduce a new family of EOP built from the classical Laguerre polynomials via two step Darboux transformations. In the present article we show that the scheme developed in [START_REF] Grandati | Solvable rational extension of translationally shape invariant potentials[END_REF][START_REF] Grandati | Solvable rational extensions of the isotonic oscillator[END_REF][START_REF] Grandati | Solvable rational extensions of the Morse and Kepler-Coulomb potentials[END_REF] admits a very natural n-step extension which allows to recover and generalize the result of [START_REF] Gómez-Ullate | Two-step Darboux transformations and exceptional Laguerre polynomials[END_REF]. We obtain new infinite towers of rational extensions of the isotonic potential strictly isospectral to this last. Their eigenfunctions are associated to the new EOP families, the one introduced by Gomez-Ullate et al corresponding to the n = 2 particular case.

The paper is organized as follows. First we recall the basic elements of the method in its one step form and the results obtained for the isotonic system. The regularity of the RS functions is here obtained in a self consistent way, using disconjugacy properties of the Schrödinger equation [START_REF] Hartman | Ordinary differential equations[END_REF][START_REF] Coppel | Disconjugacy[END_REF][START_REF] Bôcher | Leçons sur les méthodes de Sturm[END_REF] without recourse of the Kienast-Lawton-Hahn theorem [START_REF] Szegö | Orthogonal polynomials[END_REF][START_REF] Erdélyi | Higher transcendental functions[END_REF]. Then we present on a formal level the general m-step scheme and in the specific case of the isotonic oscillator, we give a criterion of concrete application in terms of asymptotic behaviour near the origin. Next we consider in a detailed manner the particular m = 2 case, initially envisaged by Gomez-Ullate et al [START_REF] Gómez-Ullate | Two-step Darboux transformations and exceptional Laguerre polynomials[END_REF]. We then prove by induction a sufficient condition to be able to build a m-step chain of regular and strictly isospectral extended potentials. We finally show that the shape invariance properties of the isotonic potential are hereditary and reached by all these extended potentials.

During the writing of this article, two papers on the same subject appeared. The first one [START_REF] Odake | Exactly Solvable Quantum Mechanics and Infinite Families of Multi-indexed Orthogonal Polynomials[END_REF] contains results parallels to those developed here. The authors adopt a point of view similar to the one initiated in [START_REF] Grandati | Solvable rational extension of translationally shape invariant potentials[END_REF][START_REF] Grandati | Solvable rational extensions of the isotonic oscillator[END_REF][START_REF] Grandati | Solvable rational extensions of the Morse and Kepler-Coulomb potentials[END_REF] and extended in the present paper. The so-called "virtual states deletion method" corresponding exactly to the use of DBT based on unphysical eigenfunctions associated to regularized RS functions. The second one [START_REF] Quesne | Higher-order SUSY, exactly solvable potentials, and exceptional orthogonal polynomials[END_REF] discusses quite extensively the 2 and 3 steps cases under a slightly different angle and proposes an interesting conjecture on possible "degeneracies" in the set of the new generated ELP. The content of the present article has been established independently of these two works.

II. DARBOUX-B ÄCKLUND TRANSFORMATIONS (DBT) AND DISCRETE SYMMETRIES

If ψ λ (x; a) is an eigenstate of H(a) = -d 2 /dx 2 + V (x; a), a ∈ R m , x ∈ I ⊂ R, associated to the eigenvalue E λ (a) (E 0 (a) = 0) ψ ′′ λ (x; a) + (E λ (a) -V (x; a)) ψ λ (x; a) = 0, (1)
then the Riccati-Schrödinger (RS) function w λ (x; a) = -ψ ′ λ (x; a)/ψ λ (x; a) satisfies the corresponding Riccati-Schrödinger (RS) equation [START_REF] Grandati | Rational solutions for the Riccati-Schrödinger equations associated to translationally shape invariant potentials[END_REF]

-w ′ λ (x; a) + w 2 λ (x; a) = V (x; a) -E λ (a). (2)
It is a well-known fact that the set of general Riccati equations is invariant under the group G of smooth SL(2, R)-valued curves M ap(R, SL(2, R)) [START_REF] Cariñena | Integrability of Riccati equation from a group theoretical viewpoint[END_REF][START_REF] Cariñena | Group theoretical approach to the intertwined hamiltonians[END_REF]. The particular subclass of Riccati-Schrödinger equations is, as for it, preserved by a specific subset of G. These transformations, called Darboux-Bäcklund Transformations (DBT), are built from any solution w ν (x; a) of the initial RS equation Eq [START_REF] Gómez-Ullate | Supersymmetry and algebraic Darboux transformations[END_REF] as [START_REF] Grandati | Rational solutions for the Riccati-Schrödinger equations associated to translationally shape invariant potentials[END_REF][START_REF] Cariñena | Integrability of Riccati equation from a group theoretical viewpoint[END_REF][START_REF] Cariñena | Group theoretical approach to the intertwined hamiltonians[END_REF] w λ (x; a)

A(wν ) → w (ν) λ (x; a) = -w ν (x; a) + E λ (a) -E ν (a) w ν (x; a) -w λ (x; a) , (3)
where

E λ (a) > E ν (a). w (ν)
λ is then a solution of the RS equation:

-w

(ν)′ λ (x; a) + w (ν) λ (x; a) 2 = V (ν) (x; a) -E λ (a), (4)
with the same energy E λ (a) as in Eq( 2) but with a modified potential

V (ν) (x; a) = V (x; a) + 2w ′ ν (x; a). (5)
This can be schematically resumed as

   w λ A(wν ) w (ν) λ V A(wν ) V (ν) . (6)
The corresponding eigenstate of

H (ν) (a) = -d 2 /dx 2 + V (ν) (x; a) can be written ψ (ν) λ (x; a) = exp -dxw (ν) λ (x; a) ∼ 1 E λ (a) -E ν (a) A (w ν ) ψ λ (x; a), (7)
where A (a) is a first order operator given by

A (w ν ) = d/dx + w ν (x; a). (8)
Eq [START_REF] Quesne | Exceptional orthogonal polynomials, exactly solvable potentials and supersymmetry[END_REF] can still be written as

ψ (ν) λ (x; a) ∼ W (ψ ν , ψ λ | x) ψ ν (x; a) , (9)
where W (y 1 , y 2 | x) is the wronskian of the functions y 1 , y 2 .

From V , the DBT generates a new potential V (ν) (quasi) isospectral to the original one and its eigenfunctions are directly obtained from those of V via Eq [START_REF] Quesne | Exceptional orthogonal polynomials, exactly solvable potentials and supersymmetry[END_REF]. Nevertheless, in general, w ν (x; a) and then the transformed potential V (ν) (x; a) are singular at the nodes of ψ ν (x; a). For instance, if ψ n (x; a) (ν = n) is a bound state of H(a), V (n) is regular only when n = 0, that is when ψ n=0 is the ground state of H, and we recover the usual SUSY partnership in quantum mechanics [START_REF] Cooper | Supersymmetry in Quantum Mechanics[END_REF][START_REF] Dutt | Supersymmetry, shape invariance and exactly solvable potentials[END_REF].

We can however envisage to use any other regular solution of Eq(2) as long as it has no zero on the considered real interval I, even if it does not correspond to a physical state. For some systems, it is possible to obtain such solutions by using specific discrete symmetries Γ i which are covariance transformations for the considered family of potentials

a Γi → a i V (x; a) Γi → V (x; a i ) = V (x; a) + U (a) . (10)
Γ i acts on the parameters of the potential and transforms the RS function of a physical excited eigenstate

w n into a unphysical RS function v n,i (x; a) = Γ i (w n (x; a)) = w n (x; a i ) associated to the negative eigenvalue E n,i (a) = Γ i (E n (a)) = U (a) -E n (a i ) < 0. -v ′ n,i (x; a) + v 2 n,i (x; a) = V (x; a) -E n,i (a). ( 11
)
To v n,i corresponds an unphysical eigenfunction of H(a)

φ n,i (x; a) = exp -dxv n,i (x; a) (12)
associated to the eigenvalue E n,i (a).

If the transformed RS function v n,i (x; a) of Eq( 11) is regular on I, it can be used to build a regular extended potential (see Eq [START_REF] Gómez-Ullate | Exceptional orthogonal polynomials and the Darboux transformation[END_REF] and Eq(7))

V (n,i) (x; a) = V (x; a) + 2v ′ n,i (x; a) (13)
(quasi)isospectral to V (x; a). The eigenstates of V (n,i) are given by (see Eq(3))

   w (n,i) k (x; a) = -v n,i (x; a) + E k (a)-En,i(a) vn,i(x;a)-w k (x;a) ψ (n,i) k (x; a) = exp -dxw (n,i) k (x; a) ∼ 1 √ E k (a)-En,i(a) A (v n,i ) ψ k (x; a), (14)
for the respective energies E k (a). The nature of the isospectrality depends if 1/φ n,i (x; a) satisfies or not the appropriate boundary conditions. If it is the case, then 1/φ n,i (x; a) is a physical eigenstate of H (n,i) (a) = -d 2 /dx 2 + V (n,i) (x; a) for the eigenvalue E n,i (a) and we only have quasi-isospectrality between V (x; a) and V (n,i) (x; a). If it is not the case, the isospectrality between V (n,i) (x; a) and V (x; a) is strict.

The above construction can be summarized by the following diagram

E k > 0 : w k A(vn,i) w (n,i) k }Physical RS functions Γ j ↓ E k,j < 0 : v k,j A(vn,i) v (n,i) k,j } Regularized (unphysical) RS functions V A(vn,i) V (n,i) Potentials (15) with ψ (n,i) k = A (v n,i ) ψ k φ (n,i) k,j = A (v n,i ) φ k,j . (16)
This procedure can be viewed as a "generalized SUSY QM partnership" where the DBT can be based on excited states RS functions properly regularized by the symmetry Γ j .

III. DISCONJUGACY AND REGULAR EXTENSIONS

To control the regularity of v n,i we can make use of the disconjugacy properties of the Schrödinger equation for negative eigenvalues.

A second order differential equation like Eq(1) is said to be disconjugated on I ⊂ R (V (x; a) is supposed to be continuous on I) if every solution of this equation has at most one zero on I [START_REF] Hartman | Ordinary differential equations[END_REF][START_REF] Coppel | Disconjugacy[END_REF]. As it is well known, this zero is necessarily simple and at this value the considered solution changes its sign. For a closed or open interval I, the disconjugacy of Eq(1) is equivalent to the existence of solutions of this equation which are everywhere non zero on I [START_REF] Hartman | Ordinary differential equations[END_REF][START_REF] Coppel | Disconjugacy[END_REF]. In the following we will consider I = ]0, +∞[.

We have also the following result:

Theorem [START_REF] Hartman | Ordinary differential equations[END_REF][START_REF] Coppel | Disconjugacy[END_REF] If there exists a continuously differentiable solution on I of the Riccati inequation

-w ′ (x) + w 2 (x) + G(x) ≤ 0, (17)
then the equation

ψ ′′ (x) + G(x)ψ(x) = 0 ( 18
)
is disconjugated on I.

In our case, since E n,i (a) ≤ 0, we have

-w ′ 0 (x; a) + w 2 0 (x; a) = V (x; a) ≤ V (x; a) -E n,i (a), (19)
w 0 (x; a) being continuously differentiable on I. The above theorem ensure the existence of nodeless solutions φ(x; a) of Eq [START_REF] Gómez-Ullate | The Darboux transformation and algebraic deformations of shape invariant potentials[END_REF] with E λ (a) = E n,i (a), that is, of regular RS functions v(x; a) solutions of Eq [START_REF] Odake | Infinitely many shape invariant potentials and new orthogonal polynomials[END_REF]. To prove that a given solution φ(x; a) belongs to this category, it is sufficient to determine the signs of the limit values φ(0 + ; a) and φ(+∞; a). If they are identical then φ is nodeless and if they are opposite, then φ presents a unique zero on I.

In the first case V (x; a) + 2v ′ (x; a) constitutes a regular (quasi)isospectral extension of V (x; a). In fact, Bôcher [START_REF] Bôcher | Leçons sur les méthodes de Sturm[END_REF] has even established a more precise theorem:

Bôcher's disconjugacy theorem If we can find u continuously differentiable on [a, b] such that Consider first the case where we have G(x) > 0, ∀x ∈ ]c, d[. We can then take u(x) = 0 which implies that any solution ψ(x) of Eq [START_REF] Grandati | Solvable rational extension of translationally shape invariant potentials[END_REF] as well as its derivative ψ ′ (x) have at most one zero in ]c, d[. If ξ 0 is the corresponding zero of ψ ′ (x) then ψ(x) is strictly monotonous on ]c, ξ 0 [ and ]ξ 0 , d[.

u ′ (x) + u 2 (x) < G(x) (20)
If we consider now the case of Eq(1) on I = ]0, +∞[ for a strictly negative eigenvalue E n,i (a) < 0. From Eq [START_REF] Grandati | Solvable rational extensions of the isotonic oscillator[END_REF] we can take u(x) = -w 0 (x; a) to apply Bôcher's theorem from which we deduce that not only Eq(1) is disconjugate on ]0, +∞[ but also the SUSY partner ψ (0) (x; a) of any solution of Eq(1) can have at most one zero on ]0, +∞[. Indeed (cf Eq(9))

ψ(x; a) = W (ψ 0 , ψ | x) ψ 0 (x; a) = ψ (0) (x; a). (21)
IV. ISOTONIC OSCILLATOR Consider the isotonic oscillator potential (ie the radial effective potential for a three dimensional isotropic harmonic oscillator with zero ground-state energy)

V (x; ω, a) = ω 2 4 x 2 + a(a -1) x 2 + V 0 (ω, a), x > 0, (22)
with a = l + 1 ≥ 1 and V 0 (ω, a) = -ω a + 1 2 . It is the unique exceptional primary translationally shape invariant potential of the second category [START_REF] Grandati | Rational solutions for the Riccati-Schrödinger equations associated to translationally shape invariant potentials[END_REF]. The corresponding Schrödinger equation is the Liouville form of the Laguerre equation on the positive half-line and its physical spectrum, associated to the asymptotic Dirichlet boundary conditions

ψ 0 + ; ω, a = 0 = ψ (+∞; ω, a) (23)
is given by

(z = ωx 2 /2, α = a -1/2) E n (a) = 2nω, ψ n (x; ω, a) = x a e -z/2 L α n (z) . ( 24
)
To ψ n corresponds the RS function [START_REF] Grandati | Rational solutions for the Riccati-Schrödinger equations associated to translationally shape invariant potentials[END_REF]

w n (x; ω, a) = w 0 (x; ω, a) + R n (x; ω, a), (25)
with

w 0 (x; ω, a) = ω 2 x - a x (26)
and

R n (x; ω, a) = -2nω ωx -(2a + 1) /x- ... 2 (n -j + 1) ω ωx -(2 (a + j) -1) /x- ... 2ω ωx -(2 (a + n) -1) /x = ωxL α+1 n-1 (z) /L α n (z) . (27)
The shape invariance property satisfied by V is [START_REF] Grandati | Rational solutions for the Riccati-Schrödinger equations associated to translationally shape invariant potentials[END_REF][START_REF] Cooper | Supersymmetry in Quantum Mechanics[END_REF][START_REF] Dutt | Supersymmetry, shape invariance and exactly solvable potentials[END_REF][START_REF] Gendenshtein | Derivation of exact spectra of the Schrodinger equation by means of supersymmetry[END_REF] V (x; ω, a) = V (x; ω, a) + 2w ′ 0 (x; ω, a) = V (x; ω, a + 1) + 2ω [START_REF] Coppel | Disconjugacy[END_REF] We have three possible discrete symmetries for V (see Eq [START_REF] Bagchi | An update on PT -symmetric complexified Scarf II potential, spectral singularities and some remarks on the rationally-extended supersymmetric partners[END_REF]) which are given by 1)

ω Γ+ → (-ω) , V (x; ω, a) Γ+ → V (x; ω, a) + ω(2a + 1) w n (x; ω, a) Γ+ → v n,1 (x; ω, a) = w n (x; -ω, a), (29)
2)

a Γ- → 1 -a, V (x; ω, a) Γ- → V (x; ω, a) + ω(2a -1) w n (x; ω, a) Γ- → v n,2 (x; ω, a) = w n (x; ω, 1 -a), (30) 3)
(ω, a)

Γ3=Γ+•Γ- → (-ω, 1 -a) V (x; ω, a) Γ3 → V (x; ω, a) + 2ω w n (x; ω, a) Γ3 → v n,3 (x; ω, a) = w n (x; -ω, 1 -a). (31)
In the (ω, α) parameters plane, Γ + and Γ -correspond respectively to the reflections with respect to the axes ω = 0 and α = 0. The RS functions v n,i , i = +, -, 3, satisfy the respective RS equations

-v ′ n,i (x; ω, a) + v 2 n,i (x; ω, a) = V (x; ω, a) -E n,i (ω, a) , (32) with 
  E n,+ (ω, a) = E -(n+a+1/2) (ω) < 0 E n,-(ω, a) = E n+1/2-a (ω) E n,3 (ω, a) = E -(n+1) (ω) < 0. ( 33
)
These eigenvalues are always negative in the i = + and i = 3 cases and the inequality E n,-(ω, a) ≤ 0 necessitates to satisfy the constraint α = a -1/2 > n. When E n,i (ω, a) ≤ 0 the Schrödinger equation

φ ′′ (x; ω, a) + (E n,i (ω, a) -V (x; ω, a)) φ(x; ω, a) = 0 ( 34
)
is disconjugated on ]0, +∞[. We are then sure that everywhere non zero solutions of Eq [START_REF] Cariñena | Integrability of Riccati equation from a group theoretical viewpoint[END_REF] exist, the question being to if φ n,i (x; ω, a) = exp -dxv n,i (x; ω, a) is such a solution.

V. L1 AND L2 SERIES OF EXTENSIONS AND CORRESPONDING EXCEPTIONAL LAGUERRE POLYNOMIALS

The L1 and L2 series of extensions is obtained using the Γ + and Γ -symmetries respectively in which case we have synthetically

E n,i (ω, a) = -2ω (a + i (n + 1/2)) φ n,i (x; ω, a) = φ 0,i (x; ω, a)L iα n (-iz) , i = ±, (35)
with

φ 0,i (x; ω, a) ∼ x iα+1/2 exp (iz/2) , ( 36
)
where α is supposed to statisfy the constraint α > -in. The corresponding RS function is given by

v n,i (x; ω, a) = v 0,i (x; ω, a) + Q n,i (x; ω, a), (37)
with

v 0,i (x; ω, a) = -i ω 2 x - iα + 1/2 x = - 1 x (i (z + α) + 1/2) ( 38
)
and

Q n,i (x; ω, a) = -iωxL iα+1 n-1 (-iz) /L iα n (-iz) . (39)
Since [START_REF] Szegö | Orthogonal polynomials[END_REF][START_REF] Erdélyi | Higher transcendental functions[END_REF]

   L α n (x) → x→0 + (α+1) n n! = n+α n L α n (x) ∼ x→+∞ (-1) n n! x n , (40)
where (X) n = (X) ... (X + n -1) is the usual Pochhammer symbol [START_REF] Erdélyi | Higher transcendental functions[END_REF], we have for α > 1/2

φ n,+ (+∞; ω, a) = +∞ φ n,+ (0 + ; ω, a) = 0 + , (41)
and for α > n

φ n,+ (+∞; ω, a) = 0 ± φ n,+ (0 + ; ω, a) = ±∞, (42)
with ± = (-1) n . Because of the constraint α > -in, Eq(34) is disconjugated both for i = + and i = -. Combined to Eq [START_REF] Sukumar | Supersymmetry, factorization of the Schrödinger equation and an hamiltonian hierarchy[END_REF] and Eq [START_REF] Andrianov | The factorization method and quantum systems with equivalent energy spectra[END_REF], this ensures that φ n,i keeps a constant strictly positive sign on ]0, +∞[ which means that v n,i (x; ω, a) is regular on this interval. Note that this result implies in particular that is L α n (x) nodeless on the negative half line and that L -α n (x) has no zero on the positive half line when α > n, which is in agreement with the Kienast-Lawton-Hahn theorem [START_REF] Szegö | Orthogonal polynomials[END_REF][START_REF] Erdélyi | Higher transcendental functions[END_REF].

V (n,i) (x; ω, a) = V (x; ω, a) + 2v ′ n,i (x; ω, a) is then a regular extension of V (x; ω, a). Since 1/φ n,i (x; ω, a) diverges at the origin, it cannot be a physical eigenstate of H (n,i) . Consequently, V (n,i) and V are strictly isospectral.

The (unnormalized) physical eigenstates of V (n,i) , ψ

(n,i) k (x; ω, a) = exp -dxw (n,i) k (x; ω, a) , satisfy    H (n,i) (ω, a)ψ (n,i) k (x; ω, a) = E k (ω) ψ (n,i) k (x; ω, a) -w (n,i)′ k (x; ω, a) + w (n,i) k (x; ω, a) 2 = V (n,i) (x; ω, a) -E k (ω) , (43)
where, using the shape invariance property of V (cf Eq(28)), we can write

V (n,i) (x; ω, a) = V (x; ω, a -1 ) + 2Q ′ n,i (x; ω, a). (44)
From Eq [START_REF] Odake | Infinitely many shape invariant potentials and cubic identities of the Laguerre and Jacobi polynomials[END_REF] we have also

ψ (n,i) k (x; ω, a) = A (v n,i ) ψ k (x; ω, a), k ≥ 0. ( 45
)
More precisely, using the identities [START_REF] Szegö | Orthogonal polynomials[END_REF][START_REF] Erdélyi | Higher transcendental functions[END_REF]]

L α n (x) + L α+1 n-1 (x) = L α+1 n (x) xL α+1 n-1 (x) -αL α n-1 (x) -nL α-1 n (x) = 0, (46)
we obtain

ψ (n,i) k (x; ω, a) = (v n,i (x; ω, a) -w k (x; ω, a)) ψ k (x; ω, a) (47)
∼ L i n,k,α (z) x a+i exp(-z/2) L iα n (-iz)
, where

L + n,k,α (z) = L α n (-z) L α+1 k (z) + L α+1 n-1 (-z) L α k (z) L - n,k,α (z) = (k + n + α) L α k (z) L -α n (z) -(-n + α) L α k (z) L -α n-1 (z) -(k + α) L α k-1 (z) L -α n (z) , (48)
are polynomials of degree n + k, namely exceptional Laguerre polynomials (ELP) of the L1 and L2 series respectively. From the orthogonality conditions on the eigenstates for fixed values of n and a, we retrieve the fact that the ELP L i n,k,α (z) constitute orthogonal families with the corresponding weights

W i n (z) = z a+i/2 exp(-z) (L iα n (-iz)) 2 . ( 49
)
VI. FORMAL SCHEME FOR MULTI-STEP DBT

A. General elements

The question of successive iterations of DBT is very natural and is at the center of the construction of the hierarchy of hamiltonians in the usual SUSY QM scheme [START_REF] Sukumar | Supersymmetry, factorization of the Schrödinger equation and an hamiltonian hierarchy[END_REF]. Staying at the formal level, it can be simply described by the following straightforward generalization of Eq(6)

     w λ A(wν 1 ) w (ν1) λ A(w (ν 1 ) ν 2 ) w (ν1,ν2) λ ... A(w (ν1,,...,ν m-1 ) νm ) w (ν1,...,νm) λ V A(wν 1 ) V (ν1) A(w (ν 1 ) ν 2 )
V (ν1,ν2) ...

A(w (ν1,,...,ν m-1 ) νm ) V (ν1,...,νm) , (50)
where w (ν1,...,νm) λ is a RS function associated to the eigenvalue E λ of the potential

V (ν1,...,νm) (x; a) = V (x; a) + 2 m-1 j=1 w (ν1,,...,νj-1) νj (x; a) ′ . ( 51
)
The corresponding eigenfunction is given by (cf Eq [START_REF] Quesne | Exceptional orthogonal polynomials, exactly solvable potentials and supersymmetry[END_REF] and Eq(9))

ψ (ν1,...,νm) λ (x; a) = A w (ν1,...,νm-1) νm ψ (ν1,...,νm-1) λ (x; a) = A w (ν1,...,νm-1) νm ... A (w ν1 ) ψ λ (x; a), (52)
that is,

ψ (ν1,...,νm) λ (x; a) = w (ν1,...,νm-1) νm (x; a) -w (ν1,...,νm-1) λ (x; a) ψ (ν1,...,νm-1) λ (x; a) (53)
Other useful representations for the extended potentials and of their eigenfunctions are possible. From Sylvester's theorem we can deduce [START_REF] Muir | A treatise on the theory of determinants[END_REF] the following property for the Wronskians

W (y 1 , ..., y m , y | x) = W (W (y 1 , ..., y m ) , W (y 1 , ..., y m-1 , y) | x) W (y 1 , ..., y m-1 | x) . (55)
Using the identity [START_REF] Muir | A treatise on the theory of determinants[END_REF] W (uy 1 , ...,

uy m | x) = u m W (y 1 , ..., y m | x) , (56)
this gives

W (y 1 , ..., y m , y | x) W (y 1 , ..., y m | x) = W W (y1,...,ym) W (y1,...,ym-1) , W (y1,...,ym-1,y) W (y1,...,ym-1) | x W (y1,...,ym|x) W (y1,...,ym-1|x) (57)
and comparing to Eq [START_REF] Shnol | Equidistant spectra of anharmonic oscillators[END_REF], we obtain ψ (ν1,...,νm) λ

(x; a) = W (ψ ν1 , ..., ψ νm , ψ λ | x) W (ψ ν1 , ..., ψ νm | x) , (58)
which is the well known Crum formula [START_REF] Crum | Associated Sturm-Liouville systems[END_REF][START_REF] Matveev | Darboux Transformations and Solitons[END_REF] for the eigenfunctions. Inserting this result in Eq(51), we then deduce the Crum formula for the potential [START_REF] Crum | Associated Sturm-Liouville systems[END_REF][START_REF] Matveev | Darboux Transformations and Solitons[END_REF] V (ν1,...,νm) (x; a) = V (x; a) + 2 (log W (ψ ν1 , ...,

ψ νm | x)) ′′ . ( 59
)
Since all the functions ψ νj implied in the wronskians in Eq(58) and Eq(59) are eigenfunctions of the same hamiltonian H(a), the properties of the determinants allow us to replace in these wronskians the even derivatives ψ by -E νj l ψ ′ νj . For instance

W (ψ ν1 , ..., ψ ν 2l | x) = ψ ν1 (x; a) ... ψ ν 2l (x; a) ψ ′ ν1 (x; a) ... ψ ′ ν 2l (x; a) ... ... (-E ν1 ) l-1 ψ ν1 (x; a) ... (-E ν 2l ) l-1 ψ ν 2l (x; a) (-E ν1 ) l-1 ψ ′ ν1 (x; a) ... (-E ν 2l ) l-1 ψ ′ ν 2l (x; a) (60)
= ∆ (ν1,...,ν 2l ) (x; a)

2l j=1 ψ νj (x; a) , where
∆ (ν1,...,ν 2l ) (x; a) = (-1) l 1 ... 1 w ν1 (x; a) ... w ν 2l (x; a) ... ... (E ν1 ) l-1 ... (E ν 2l ) l-1 (E ν1 ) l-1 w ν1 (x; a) ... (E ν 2l ) l-1 w ν 2l (x; a) (61)
and

W ψ ν1 , ..., ψ ν 2l+1 | x = ψ ν1 (x; a) ... ψ ν 2l (x; a) ψ ′ ν1 (x; a) ... ψ ′ ν 2l (x; a) ... ... (-E ν1 ) l-1 ψ ′ ν1 (x; a) ... -E ν 2l+1 l-1 ψ ′ ν 2l+1 (x; a) (-E ν1 ) l ψ ν1 (x; a) ... -E ν 2l+1 l ψ ν 2l+1 (x; a) (62)
= ∆ (ν1,...,ν 2l+1 ) (x; a)

2l+1 j=1 ψ νj (x; a) , where ∆ (ν1,...,ν 2l+1 ) (x; a) = 1 ... 1 w ν1 (x; a) ... w ν 2l+1 (x; a) ... ... (E ν1 ) l-1 w ν1 (x; a) ... (E ν 2l ) l-1 w ν 2l+1 (x; a) (E ν1 ) l ... (E ν 2l +1 ) l . ( 63
)
Then we obtain the representation of the ψ (ν1,...,νm) λ in terms the so-called Crum-Krein determinants ∆ (ν1,...,νm) [START_REF] Crum | Associated Sturm-Liouville systems[END_REF][START_REF] Krein | On a continuous analogue of a Christoffel formula from the theory of orthogonal polynomials[END_REF][START_REF] Bagrov | Darboux transformation, factorization, and supersymmetry in one-dimensional quantum mechanics[END_REF] ψ (ν1,...,νm) λ (x; a) = ψ λ (x; a) ∆ (ν1,...,νm,λ) (x; a) ∆ (ν1,...,νm) (x; a) .

Until now, all these results are purely formal. The central problem is now to choose in an appropriate way the family of eigenfunctions (ψ ν1 , ..., ψ νm ) from which the successive DBT are built (see Eq(58) and Eq [START_REF] Fernández | Second order SUSY transformations with 'complex energies[END_REF]) in order to ensure the regularity of the successive extended potentials V (ν1,...,νm) .

A first answer has been given more than half a century ago by Crum [START_REF] Crum | Associated Sturm-Liouville systems[END_REF]. The Crum proposal [START_REF] Crum | Associated Sturm-Liouville systems[END_REF] is to take the successive eigenstates of the discrete bound states spectrum starting from the ground state (ν 1 , ..., ν m ) = (0, ..., m -1) (65)

Krein [START_REF] Krein | On a continuous analogue of a Christoffel formula from the theory of orthogonal polynomials[END_REF][START_REF] Bagrov | Darboux transformation, factorization, and supersymmetry in one-dimensional quantum mechanics[END_REF][START_REF] Bagrov | Darboux transformation and elementary exact solutions of the Schrödinger equation[END_REF] extended this result showing it is possible to choose sets of unnecessarily successive eigenstates if they satisfy a given condition, namely

(ν 1 , ..., ν m ) = (n 1 , ..., n m ) ∈ N m , (66)
with

(n -n 1 )(n -n 2 )...(n -n m ) ≥ 0, ∀n ∈ N. ( 67
)
This last is verified if the set (n 1 , ..., n m ) is constituted by "aggregates" of an even number of eigenstates. The possibility to employ sets of two juxtaposed eigenstates has been rediscovered later by Adler [START_REF] Adler | A modification of Crum method[END_REF]. These results have been used extensively in the context of higher order SUSY (see for instance [START_REF] Sukumar | Supersymmetric quantum mechanics of one-dimensional systems[END_REF][START_REF] Sukumar | Supersymmetry, factorization of the Schrödinger equation and an hamiltonian hierarchy[END_REF][START_REF] Andrianov | The factorization method and quantum systems with equivalent energy spectra[END_REF][START_REF] Andrianov | Higher-derivative supersymmetry and the Witten index[END_REF][START_REF] Bagrov | Darboux transformation, factorization, and supersymmetry in one-dimensional quantum mechanics[END_REF][START_REF] Bagrov | Darboux transformation and elementary exact solutions of the Schrödinger equation[END_REF][START_REF] Samsonov | Darboux transformation and exactly solvable potentials with quasi-equidistant spectrum[END_REF][START_REF] Bagrov | N-order Darboux transformation and a spectral problem on semiaxis[END_REF][START_REF] Samsonov | New possibilities for supersymmetry breakdown in quantum mechanics and second order irreducible Darboux transformations[END_REF][START_REF] Fernández | A simple generation of exactly solvable anharmonic oscillators[END_REF][START_REF] Fernández | Second-order supersymmetric periodic potentials[END_REF][START_REF] Fernández | Second order SUSY transformations with 'complex energies[END_REF][START_REF] Mielnik | The finite difference algorithm for higher order supersymmetry[END_REF]). In [START_REF] Samsonov | New possibilities for supersymmetry breakdown in quantum mechanics and second order irreducible Darboux transformations[END_REF], Samsonov has still extended the Krein-Adler result. Note finally that Fernandez et al [START_REF] Fernández | A simple generation of exactly solvable anharmonic oscillators[END_REF] employed the specific "Wick rotated" eigenfunctions of negative energies introduced by Shnol' [START_REF] Grandati | Solvable rational extension of translationally shape invariant potentials[END_REF][START_REF] Samsonov | Darboux transformation and exactly solvable potentials with quasi-equidistant spectrum[END_REF][START_REF] Shnol | Equidistant spectra of anharmonic oscillators[END_REF] to build successive extensions of the harmonic potential.

B. Application to the isotonic oscillator

Our aim is to consider the possibility to generalize our construction in a multi-step version on the basis of the regularized RS functions (or unphysical eigenfunctions) obtained above. Staying at a formal level, it would correspond to the following generalization of the diagram Eq [START_REF] Sasaki | Exceptional Laguerre and Jacobi polynomials and the corresponding potentials through Darboux-Crum transformations[END_REF] (see Eq [START_REF] Fernández | Second-order supersymmetric periodic potentials[END_REF])

E k > 0 : w k A(vn 1 ,i 1 ) w (n1,i1) k A(v (n 1 ,i 1 ) n 2 ,i 2 ) w (n1,i1,n2,i2) k ... A(v (n1,i1,,...,n m-1 ,i m-1 ) nm,im
) w (n1,i1,,...,nm,im) k

Γ i ↓ E k,i < 0 : v k,i A(vn 1 ,i 1 ) v (n1,i1) k,i A(v (n 1 ,i 1 ) n 2 ,i 2 ) v (n1,i1,n2,i2) k,i ... A(v (n1,i1,,...,n m-1 ,i m-1 ) nm,im ) v (n1,i1,,...,nm,im) k,i V A(vn 1 ,i 1 ) V (n1,i1) A(v (n 1 ,i 1 ) n 2 ,i 2 ) V (n1,i1,n2,i2) ... A(v (n1,i1,,...,n m-1 ,i m-1 ) nm,im ) V (n1,i1,,...,nm,im) (68) with v (n0) k = v k .
Under this way, we generate chains of isospectral extensions V (n1,i1,,...,nm,im) of the potential V with (see Eq [START_REF] Fernández | Second order SUSY transformations with 'complex energies[END_REF] and Eq(59))

V (n1,i1,...,nm,im) (x; a) = V (x; a) + 2 v (n1,i1,...,nm-1,im-1) nm,im (x; a) ′ (69) = V (x; a) + 2 m j=1 v (n1,i1,...,nj-1,ij-1) nj ,ij (x; a) ′ = V (x; a) + 2 (log W (φ n1,i1 , ..., φ nm,im | x)) ′′ w (n1,i1,,...,nm,im) k
and the v (n1,i1,,...,nm,im) k,i are RS functions of the extended potential V (n1,i1,...,nm,im) associated respectively to the eigenvalues E k and E k,i .

The corresponding eigenfunctions are given by (see Eq [START_REF] Shnol | Equidistant spectra of anharmonic oscillators[END_REF] and Eq(58)) = W (φn 1 ,i 1 ,...,φn m,im ,φ k,i |x) W (φn 1 ,i 1 ,...,φn m,im |x) .

           ψ (n1,
(70)

Suppose that at the step m-1 all the extended potentials in the chain are regular and strictly isospectral. The problem is now to control the regularity of the potential V (n1,i1,...,nm,im) and of its physical eigenstates ψ (n1,i1,,...,nm-1,im-1) k as well as its strict isospectrality with V (n1,i1,...,nm-1,im-1) . This is achieved if the unphysical eigenfunction φ (n1,i1,,...,nm-1,im-1) nm,im , associated to the DBT A(v (n1,i1,,...,nm-1,im-1) nm,im

), is nodeless and if 1/φ (n1,i1,,...,nm-1,im-1) nm,im cannot be into the set of physical eigenstates of V (n1,i1,...,nm-1,im-1) .

The disconjugacy of the Schrödinger equation φ ′′ (x; ω, a) + E nm,im (ω, a) -V (n1,i1,...,nm-1,im-1) (x; ω, a) φ(x; ω, a) = 0, (71) for a negative eigenvalue E nm,im guarantees that v (n1,i1,...,nm-1,im-1) nm,im is at most singular in one point. The regularity of the extended potential V (n1,i1,...,nm,im) and its strict isospectrality with the preceding one in the chain are then satisfied as soon as the unphysical eigenfunction φ (n1,i1,,...,nm-1,im-1) nm,im satisfies appropriate boundary conditions, which was already the argument used above in the one step case. Namely, the strict isospectrality and the regularity are satisfied when φ (n1,i1,...,nm-1,im-1) nm,im tends to ∞ at one extremity of ]0, +∞[ and to 0 at the other extremity, with the same sign for both limits.

In the case of the isotonic oscillator, there are in fact three type of regularization transformations Γ i , i = +, -, 3 associated respectively to the three series of EOP L1, L2 and L3 [START_REF] Quesne | Solvable rational potentials and exceptional orthogonal polynomials in supersymmetric quantum mechanics[END_REF][START_REF] Grandati | Solvable rational extensions of the isotonic oscillator[END_REF]. At each step we can a priori use any of the three type if regularized RS function. However, since the use of DBT associated to RS functions of the L3 series leads only to quasi isospectrality, the argument above is no more valid. Consequently, in the following we exclude this case and consider only DBT built from the RS functions of the L1 or L2 series.

From Eq(64), we see that in the present case φ (n1,i1,...,nm-1,im-1) nm,im (x; ω, a) = φ 0,im (x; ω, a) R (n1,i1,,...,nm,im) (x; ω, a) ,

where R (n1,i1,,...,nm,im) (x; ω, a) = L imα nm (-i m z) ∆ (n1,i1,,...,nm,im) (x; ω, a) ∆ (n1,i1,,...,nm-1,im-1) (x; ω, a)

is a rational function. Due to the presence of the exponential in φ 0,im , the behaviour at infinity of φ (n1,i1,...,nm-1,im-1) nm,im is the same as for φ 0,im , namely

   φ (n1,i1,...,nm-1,im-1) nm,+ (x; ω, a) ∼ x→+∞ exp (z/2) x α+1/2 R (n1,i1,,...,nm,+) (+∞; ω, a) → x→+∞ ±∞ φ (n1,i1,...,nm-1,im-1) nm,- (x; ω, a) ∼ x→+∞ exp (-z/2) x -α+1/2 R (n1,i1,,...,nm,-) (+∞; ω, a) → x→+∞ 0 ± , (74)
(the sign of the limit being the one of R (n1,i1,,...,nm,+) (+∞; ω, a) or R (n1,i1,,...,nm,-) (+∞; ω, a) respectively). The absence of node for φ (n1,i1,...,nm-1,im-1) nm,im and the strict isospectrality between V (n1,i1,...,nm-1,im-1) and V (n1,i1,...,nm,im) , are simultaneously ensured if and only if we have

   φ (n1,i1,...,nm-1,im-1) nm,+ (x; ω, a) ∼ x→0 + x α+1/2 R (n1,i1,,...,nm,+) (0 + ; ω, a) → x→0 + 0 ± φ (n1,i1,...,nm-1,im-1) nm,- (x; ω, a) ∼ x→0 + x -α+1/2 R (n1,i1,,...,nm,-) (0 + ; ω, a) → x→0 + ±∞, (75)
(the sign of the limit being the one of R (n1,i1,,...,nm,+) (0 + ; ω, a) or R (n1,i1,,...,nm,-) (0 + ; ω, a) respectively), with sign R (n1,i1,,...,nm,im) 0 + ; ω, a = sign R (n1,i1,,...,nm,im) (+∞; ω, a) .

VII. TWO-STEP CASE

Before to consider the general case, we first consider the two-step case. We then have to check that the unphysical eigenstates φ (n1,i1) n2,i2 satisfy the appropriate boundary conditions. Using Eq [START_REF] Cariñena | Group theoretical approach to the intertwined hamiltonians[END_REF] in Eq(70), we have for m = 2

φ (n1,i1) n2,i2 (x; ω, a) = W φ 0,i1 (x; ω, a)L i1α n1 (-i 1 z) , φ 0,i2 (x; ω, a)L i2α n2 (-i 2 z) | x φ 0,i1 (x; ω, a)L i1α n1 (-i 1 z) ,
where α > -i j n j , j = 1, 2, or with Eq(56)

φ (n1,i1) n2,i2 (x; ω, a) = φ 0,i1 (x; ω, a) W L i1α n1 (-i 1 z) , (x α exp (z/2)) i2-i1 L i2α n2 (-i 2 z) | x L i1α n1 (-i 1 z) . ( 76
) If i 1 = i 2 = i, then φ (n1,i) n2,i (x; ω, a) = φ 0,i (x; ω, a) W L iα n1 (-iz) , L iα n2 (-iz) | x L iα n1 (-iz) . ( 77
)
Using [START_REF] Muir | A treatise on the theory of determinants[END_REF] W (y 1 , ...,

y m | x) = dz dx m(m-1)/2 W (y 1 , ..., y m | z) (78)
and [START_REF] Szegö | Orthogonal polynomials[END_REF][START_REF] Erdélyi | Higher transcendental functions[END_REF] (

L α n (x)) ′ = -L α+1 n-1 (x) , (79)
this gives

φ (n1,i) n2,i (x; ω, a) = iωx iα+3/2 exp (iz/2) L iα n1 (-iz) L iα n1 (-iz) L iα n2 (-iz) L iα+1 n1-1 (-iz) L iα+1 n2-1 (-iz) . ( 80
)
We then have (cf Eq [START_REF] Matveev | Darboux Transformations and Solitons[END_REF])

       φ (n1,i) n2,i (x; ω, a) ∼ x→0 + iωx iα+3/2 ( n 1 +iα n 1 ) n1+iα n1 n2+iα n2 n1+iα n1-1 n2+iα n2-1 = iωx iα+3/2 n2+iα n2-1 n2-n1 n2 φ (n1,i) n2,i (x; ω, a) ∼ x→+∞ i n2 (n 2 -n 1 ) ω n2! x 2n2+iα-1/2 exp (iz/2) , (81)
that is,

   φ (n1,+) n2,+ (x; ω, a) ∼ x→0 + ωx (α+2)-1/2 n2+α n2-1 n2-n1 n2 → x→0 + 0 ± φ (n1,-) n2,-(x; ω, a) ∼ x→+∞ ω(n2-n1) n2! x 2n2+α-1/2 exp (z/2) → x→+∞ ±∞, , if a > 1, ( 82
)
where ± = sign(n 2 -n 1 ), and

   φ (n1,-) n2,-(x; ω, a) ∼ x→0 + (-1) n2 ωx -(α-2)-1/2 α-n2 n2-1 n2-n1 n2 → x→0 + ±∞ φ (n1,-) n2,-(x; ω, a) ∼ x→+∞ (-1) n2 ω(n2-n1) n2! x 2n2-α-1/2 exp (-z/2) → x→+∞ 0 ± , if a > sup (n 2 , n 1 ) , (83)
where ± = (-1) n2 sign(n 2 -n 1 ).

If i 1 = -i 2 = -i, then using the last of the following identities [START_REF] Andrews | Special functions[END_REF]   

(z α L α n (z)) ′ = (n + α) z α-1 L α-1 n (z) (e -z L α n (z)) ′ = -e -z L α+1 n (z) (z α e -z L α n (z)) ′ = (n + 1) z α-1 e -z L α-1 n+1 (z) , (84)
we obtain

φ (n1,-i) n2,i (x; ω, a) = ωx -iα+3/2 exp (-iz/2) 2 ω iα W L -iα n1 (iz) , z iα exp (iz) L iα n2 (-iz) | z L -iα n1 (iz) (85) = 2x iα-1/2 exp (iz/2) L -iα n1 (iz) L -iα n1 (iz) zL iα n2 (-iz) -iL -iα+1 n1-1 (iz) (n 2 + 1) L iα-1 n2 (-iz) .
We then have (cf Eq [START_REF] Matveev | Darboux Transformations and Solitons[END_REF])

   φ (n1,-i) n2,i (x; ω, a) ∼ x→0 + 2x iα-1/2 i n2+1 (α+in2)...α n2! φ (n1,-i) n2,i (x; ω, a) ∼ x→+∞ i n2+1 1 n2! ω 2 n2+1 x 2n2+iα+3/2 exp (iz/2) . (86)
Consequently (a > 1)

   φ (n1,-) n2,+ (x; ω, a) ∼ x→0 + 2x α-1/2 (α+n2)...α n2! → x→0 + 0 + φ (n1,-) n2,+ (x; ω, a) ∼ x→+∞ 1 n2! ω 2 n2+1 x 2n2+a+1 exp (z/2) → x→+∞ +∞ (87) and when α > n 2    φ (n1,+) n2,-(x; ω, a) ∼ x→0 + (-1) n2+1 2 (α-n2)...α n2! x -α-1/2 → x→0 + ±∞ φ (n1,+) n2,-(x; ω, a) ∼ x→+∞ (-1) n2+1 1 n2! ω 2 n2+1 x 2n2-a+2 exp (-z/2) → x→+∞ 0 ± , (88)
with ± = (-1) n2+1 . When a > sup(n 2 , n 1 ), we have the correct asymptotic behaviour which ensures the regularity of V (n1,i1,n2,i2) (x; ω, a) and V (n1,i1) (x; ω, a) as their strict isospectrality with V (x; ω, a) for every choice of (n 1 , i 1 , n 2 , i 2 ) with n 1 = n 2 .

We then obtain a chain with two strictly isospectral regular successive extensions of

V V A(vn 1 ,i 1 ) V (n1,i1) A(v (n 1 ,i 1 ) n 2 ,i 2 ) V (n1,i1,n2,i2) . (89)

VIII. GENERAL M-STEP CASE

Consider now the general case of a m-step DBT. We proceed by induction, supposing that we have built a chain of m -1 regular strictly isospectral potentials

V A(vn 1 ,i 1 ) V (n1,i1) A(v (n 1 ,i 1 ) n 2 ,i 2 )
...

A(v (n1,i1,,...,n m-2 ,i m-2 ) n m-1 ,i m-1 ) V (n1,i1,...,nm-1,im-1) . (90)
and that we have α > -i j n j , ∀j ≤ m. This implies in particular that the φ

(n1,i1,...,nj-1,ij-1) nj ,ij (x; ω, a), j ≤ m -1,
which are eigenfunctions of H (n1,i1,...,nm-1,im-1) associated to the respective negative eigenvalues E nj-1,ij-1 (ω, a) < 0, are all nodeless on ]0, +∞[ and that they keep a constant sign on this interval that we can always take as positive.

We have to determine what is the constraint to implement for the unphysical eigenfunction φ (n1,i1,...,nm-1,im-1) nm,im in order that it satisfies the good boundary conditions necessary to build a regular and strictly isospectral extension V (n1,i1,...,nm,im) at the next step.

From Eq(74) we know that

   φ (n1,i1,...,nm-2,im-2) nm-1,+ (x; ω, a) → x→+∞ +∞ φ (n1,i1,...,nm-2,im-2) nm-1,- (x; ω, a) → x→+∞ 0 + , (91)
and consequently we have also (cf Eq(75))

   φ (n1,i1,...,nm-2,im-2) nm-1,+ (x; ω, a) → x→0 + 0 + φ (n1,i1,...,nm-2,im-2) nm-1,- (x; ω, a) → x→0 + +∞. (92)
For x sufficiently large or in the vinicity of 0, we have E nm-1,± (ω, a) < V (n1,i1,...,nm-2,im-2) (x; ω, a). The Bôcher disconjugacy theorem then implies that in the neighbourhoods of 0 and of +∞, the functions φ (n1,i1,...,nm-2,im-2) nm-1,im-1 are strictly monotonous, that is strictly decreasing if i m-1 = -and strictly increasing if i m-1 = +.

We already know (cf Eq(74)) that φ (n1,i1,...,nm-1,im-1) nm,im

(+∞; ω, a) = ±∞, if i m = + 0 ± , if i m = - (93)
We prove now by induction the following result for the unphysical eigenfunctions of the potentials obtained at the (m -1)

th step (z = ωx 2 / 2) φ (n1,i1,...,nm-1,im-1) nm,im (x; ω, a) ∼ x→0 + φ (n1,i1,...,nm-1,im-1) nm,im (0)
x im(α+q(i1,...,im))-1/2 (1 + O(z)) , (94) where φ (n1,i1,...,nm-1,im-1) nm,im
(0) has the same sign as φ

(+∞; ω, a) and where q (i 1 , ..., i m ) = q + (i 1 , ..., i m ) -q -(i 1 , ..., i m ), q ± (i 1 , ..., i m ) is the number of i j equal to ± in the set (i 1 , ..., i m ). q ± (i 1 , ..., i m ) can be viewed as the number of state of "± charge" used in the chain of DBT at the m th step and the difference q (i 1 , ..., i m ) as the total charge associated to the chain.

We then have

   φ (n1,i1,...,nm-1,im-1) nm,+ (x; ω, a) → x→0 + 0 ± , if a > -q (i 1 , ..., i m-1 , +) + 1 = -q (i 1 , ..., i m-1 ) φ (n1,i1,...,nm-1,im-1) nm,- (x; ω, a) → x→0 + ±∞, if a > -q (i 1 , ..., i m-1 , -) = 1 -q (i 1 , ..., i m-1 ) , (95)
which are necessary conditions in order that the potentials V (n1,i1,...,nm,+) and V (n1,i1,...,nm,-) respectively are regular and strictly isospectral to V (n1,i1,,...,nm-1,im-1) . Note that these constraints are independent of the n j but depend only on the type of the states used to build the chain.

Refering to the results obtained in the 2-step case (see Eq(82), Eq(83) and Eq(88)), we can verify that this property is valid for m = 2 since q (+, +) = -q (-, -) = 2 q (+, -) = q (-, +) = 0.

(96)

Suppose that Eq(94) is verified till the (m -1) th step, that is,

φ (n1,i1,...,nm-2,im-2) n,± (x; ω, a) ∼ x→0 + x ±(α+q(i1,...,im-2))+1/2 (1 + O(z)) , (97)
with a > -q (i 1 , ..., i m-2 , +) = -1 -q (i 1 , ..., i m-2 ), if i m-1 = +, and a > 1 -q (i 1 , ..., i m-2 , -) = 2 -q (i 1 , ..., i m-2 ), if i m-1 = -.

Since φ (n1,i1,...,nm-2,im-2) nm-1,im-1 and φ (n1,i1,...,nm-2,im-2) k,i are both eigenfunctions of H (n1,i1,...,nm-2,im-2) , using the wronskian theorem [START_REF] Messiah | Mécanique quantique T[END_REF] in Eq(70) we can write the unphysical eigenfunctions at the (m -1) th step as φ (n1,i1,...,nm-1,im-1) k,i (x; ω, a) .

(x; ω, a) = E nm-
We have to consider several cases corresponding to the four possible couples (i m-1 , i). 1) Consider first the case i = i m-1 = + in which we take x 0 = 0. Using Eq(92), Eq(98) gives φ (n1,i1,...,nm-2,im- x (α+q(i1,...,im-2,+,+)

)-1/2 (1 + O(z)) , (100)
Consequently, in this case, Eq(94) is also verified at the m th step.

2) In the case i = i m-1 = -, if we take x 0 = +∞, we have similarly with Eq(91) and Eq(98)

φ (n1,i1,...,nm-2,im-2,nm-1,-) k,- (x; ω, a) ∼ +∞ x dξφ (n1,i1,...,nm-2,im-2) nm-1,- (ξ; ω, a)φ (n1,i1,...,nm-2,im-2) k,- (ξ; ω, a) φ (n1
,i1,...,nm-2,im-2) nm-1,-

(x; ω, a) . (101)
As in the preceding case, since φ x -(α+q(i1,...,im-2,-,-))-1/2 .

Eq(94) is again verified at the m th step.

3) It remains to consider the cases i = + = -i m-1 and i = -= -i m-1 . Since φ | x has the same positive sign at these two limits. Due to the disconjugacy of φ ′′ (x; ω, a) + E k,+ (ω, a) -V (n1,i1,...,nm-1,im-1) (x; ω, a) φ(x; ω, a) = 0,

(
we then deduce that φ (n1,i1,...,nm-2,im-

+ O(z). (0)
Consequently, we deduce

φ (n1,i1,...,nm-2,im-2,nm-1,-) k,+ (x; ω, a) ∼ x→0 +
x (α+q(i1,...,im-2,-,+))-

∼ x→0 + -2 (α + q (i 1 , ..., i m-2 )) φ (n1,i1,...,nm-2,im-2) k,- (0)x -(α+q(i1,...,im-2))-1/2 (1 + O(z)) ∼ x→0 + W φ n1,+ , ..., φ nq + ,+ , ..., φ nq + +1,-, ..., φ nm,-, ψ k | x (107) ∼ x (m+1)(m/2+a) e -(m+1)z/2 W e z L α n1 (-z) , ..., e z L α nq + (-z) , z -α L -α nq + +1 (z) , ..., z -α L -α nm (z), L α k (z) | x , (106)
that is, with Eq(79) and Eq(84)

W φ n1,+ , ..., φ nq + ,+ , ..., φ nq + +1,-, ..., φ nm,-, ψ k | x (108) ∼ x (m+1)(m/2+a) e (q+-(m+1)/2)z z -(m-q+)(α+m) det Ψ (n1,+,...,nq + ,+,nq

k i,j =    L α+i-1 nj (-z) , if j ≤ q + (n j -α -i + 2) i-1 z m-i+1 L -α-i+1 nj (z), if q + + 1 ≤ j ≤ m (-1) i-1 L α+i-1 k-i+1 (z) , if j = m + 1. (109)
In the same manner W φ n1,+ , ..., φ nq + ,+ , ..., φ nq + +1,-, ..., φ nm,-| x (110)

∼ x m((m-1)/2+a) e mz/2 W L α n1 (-z) , ..., L α nq + (-z) , z -α e -z L -α nq + +1 (z) , ..., z -α e -z L -α nm (z) | x , that is, W φ n1,+ , ..., φ nq + ,+ , ..., φ nq + +1,-, ..., φ nm,-| x (111)
∼ x m((m-1)/2+a) z (m-q+)(-α-m+1)) e -(m/2-q+)z det Φ (n1,+,...,nq + ,+,nq + +1,-,....,nm,-) , where (i, j ∈ {1, ..., m})

Φ (n1,+,...,nq + ,+,nq + +1,-,....,nm,-) i,j

= L α+i-1 nj -i+1 (-z) , if j ≤ q + (n j + 1) i-1 z m-i+1 L -α-i+1 nj+i-1 (z), if q + + 1 ≤ j ≤ m. (112)
Eq(70) gives then ψ (n1,+,...,nq + ,+,nq + +1,-,....,nm,-) k (x; ω, a) ∼ ψ 0 (x; ω, a + q (i 1 , ..., i m )) det Ψ (n1,+,...,nq + ,+,nq + +1,-,....,nm,-) k det Φ (n1,+,...,nq + ,+,nq + +1,-,....,nm,-)

.

The ψ (n1,i1,...,nm,im) k are then obtained as the product of the modified "gauge factor" (ie fundamental state of the initial potential) ψ 0 (x; ω, a + q (i 1 , ..., i m )), multiplied by the ratio of two polynomials, the polynomial denominator det Φ (n1,+,...,nq + ,+,nq + +1,-,....,nm,-) being common to all the eigenstates of the considered extension. The orthogonality condition on the ψ (n1,i1,...,nm,im) k ensure that the det Ψ (n1,+,...,nq + ,+,nq + +1,-,....,nm,-) k constitute an orthogonal family with respect to the weight W (n1,i1,...,nm,im) (z) = z α+q(i1,...,im) e -z det Φ (n1,+,...,nq + ,+,nq + +1,-,....,nm,-) 2 .

X. SHAPE INVARIANCE OF THE EXTENDED POTENTIALS

Consider a chain of regular strictly isospectral extensions

V A(vn 1 ,i 1 ) V (n1,i1) A(v (n 1 ,i 1 ) n 2 ,i 2 ) ... A(v (n1,i1,,...,n m-1 ,i m-1 ) nm,im ) V (n1,i1,,...,nm,im) , (114)
The superpartner of the potential V (n1,i1,...,nm,im) (x; ω, a) is given by [START_REF] Cooper | Supersymmetry in Quantum Mechanics[END_REF][START_REF] Dutt | Supersymmetry, shape invariance and exactly solvable potentials[END_REF] V (n1,i1,...,nm,im) (x; ω, a) = V (n1,i1,...,nm,im) (x; ω, a) + 2 w (n1,i1,...,nm,im) 0

(x; ω, a) ′ , (115)
w (n1,i1,...,nm,im) 0 (x; ω, a) being the RS function associated to the ground level of V (n1,i1,...,nm,im) (E 0 (ω) = 0). Since (see Eq(3)) w (n1,i1,...,nm,im) 0 (116) we have with Eq(69)

V (n1,i1,...,nm,im) (x; ω, a) = V (n1,i1,...,nm-1,im-1) (x; ω, a)

-2 E nm,im (ω, a) v (n1,i1,...,nm-1,im-1) nm,im (x; ω, a) -w (n1,i1,...,nm-1,im-1) 0

(x; ω, a) ′ . ( 118
)
We proceed by induction. Suppose that the potential V (n1,i1,...,nm-1,im-1) has the same shape invariance properties than V V (n1,i1,...,nm-1,im-1) (x; ω, a) = V (n1,i1,...,nm-1,im-1) (x; ω, a) + 2 w (n1,i1,...,nm-1,im-1) 0 (x; ω, a) ′ (119) = V (n1,i1,...,nm-1,im-1) (x; ω, a + 1) + 2ω.

As shown in [START_REF] Grandati | Solvable rational extensions of the isotonic oscillator[END_REF], this is effectively the case for m = 2. It results V (n1,i1,...,nm,im) (x; ω, a) = V (n1,i1,...,nm,im) (x; ω, a + 1) + 2ω -2 ∆ (n1,i1,...,nm,im) (x; ω, a)

′ , (120)
where ∆ (n1,i1,...,nm,im) (x; ω, a) = E nm,im (ω, a) v (n1,i1,...,nm-1,im-1) nm,im Suppose also that ∆ (n1,i1,...,nm-1,im-1) (x; ω, a) = 0, (122) which again is verified for m = 2 [START_REF] Grandati | Solvable rational extensions of the isotonic oscillator[END_REF]. We can then write v (n1,i1,...,nm-2,im-2) nm-1,im-1 (124) and (see Eq [START_REF] Quesne | Higher-order SUSY, exactly solvable potentials, and exceptional orthogonal polynomials[END_REF]) E nm,im (ω, a + 1) -E nm-1,im-1 (ω, a + 1) = E nm,im (ω, a) -E nm-1,im-1 (ω, a) .

(125)

Inserting Eq(124), Eq(123) and Eq(125) into Eq(121), we obtain after a little elementary algebra ∆ (n1,i1,...,nm,im) (x; ω, a) = 0.

(126)

By induction, this property is verified for every m. Then Eq(120) becomes V (n1,i1,...,nm,im) (x; ω, a) = V (n1,i1,...,nm,im) (x; ω, a + 1) + 2ω, (127) that is, the potential has also the same shape invariance properties than the isotonic potential. The translational shape invariance of the isotonic potential is then hereditary in such chain of extensions.

XI. CONCLUSION

We have shown that the "generalized SUSY QM partnership" that we have previously elaborated in a one-step scheme can be extended in a multi-step formulation. We have proven the necessary conditions to obtain chains of regular extensions of the isotonic potential of arbitrary length and have given explicit expressions for their eigenstates. We also established explicitely the hereditary character of the shape invariance properties of the isotonic potential which are common to all the potentials in a given chain.

The case of the generic potentials of the second category of primary TSIP [START_REF] Grandati | Rational solutions for the Riccati-Schrödinger equations associated to translationally shape invariant potentials[END_REF], namely Pöschl-Teller or Scarf potentials, can be considered in the same way. This work is in progress and a forthcoming paper is in preparation. For the first category exceptional TSIP, due to the strict isospectrality constraint, the only interesting case in which we can envisage such chain of extensions is the ERKC potential [START_REF] Grandati | Solvable rational extensions of the Morse and Kepler-Coulomb potentials[END_REF]. This is the object of further investigations.

  on all [a, b] , then nor ψ(x) = u(x)ψ(x) -ψ ′ (x) nor ψ(x) can have more than one zero in [a, b]. Moreover, if ψ(x) vanish in one point of [a, b] then ψ(x) is everywhere nonzero in this interval (and reciprocally). If this disconjugacy property is satisfied for every [a, b] ⊂ ]c, d[, then it is in fact satisfied in all the open interval ]c, d[.

  ...,νm-1) νm , ψ (ν1,...,νm-1) λ | x ψ (ν1,...,νm-1) νm (x; a) , where W (y 1 , ..., y m | x) is the wronskian of the functions y 1 , ..., y m W (y 1 , ..., y m | x) = y 1 (x) ... y m (x) ... ... y (m-1) 1 (x) ... y (m-1) m(x).

  l ψ νj and the odd derivatives ψ (2l+1) νj

(

  n1,i1,...,nm-2,im-2) nm-1,-and φ (n1,i1,...,nm-2,im-2) k,-are both strictly positive on ]0, +∞[, φ (n1,i1,...,nm-2,im-2,nm-1,-) k,-keeps a constant sign on this interval and we can always take it as positive. From Eq(74)we deduce that for any ε > 0, +∞ ε dξφ (n1,i1,...,nm-2,im-2) nm-1,-(ξ; ω, a)φ (n1,i1,...,nm-2,im-2) k,-(ξ; ω, a) is finite and positive. Since (a > 1 -q (i 1 , ..., i m-2 )) φ (n1,i1,...,nm-2,im-2) nm-1,-(x; ω, a)φ (n1,i1,...,nm-2,im-2) k,-(x; ω, a) ∼ x→0 + x -2(α+q(i1,...,im-2))+1 (1 + O(z)) , Eq(101) gives φ (n1,i1,...,nm-2,im-2,nm-1,-) k,-(x; ω, a) ∼ x→0 +

  (x; ω, a) = -v (n1,i1,...,nm-1,im-1) nm,im (x; ω, a) -E nm,im (ω, a) v (n1,i1,...,nm-1,im-1) nm,im(x; ω, a) -w (n1,i1,...,nm-1,im-1) 0 (x; ω, a) ,

  (x; ω, a) -w (n1,i1,...,nm-1,im-1) 0(x; ω, a)(121)+w (n1,i1,...,nm-1,im-1) 0(x; ω, a) + v (n1,i1,...,nm-1,im-1) nm,im(x; ω, a + 1).

  im-1 (ω, a) v (n1,i1,...,nm-2,im-2) nm-1,im-1(x; ω, a) -w (n1,i1,...,nm-2,im-2) 0 (x; ω, a). i1,...,nm-1,im-1) 0(x; ω, a) = -v (n1,i1,...,nm-2,im-2) nm-1,im-1 (x; ω, a) -En m-1 ,i m-1 (ω,a) v (n1,i1,...,n m-2 ,i m-2 ) n m-1 ,i m-1 (x;ω,a)-w (n1,i1,...,n m-2 ,i m-2 ) 0 (x;ω,a) v (n1,i1,...,nm-1,im-1) nm,im (x; ω, a) = -v (n1,i1,...,nm-2,im-2) nm-1,im-1 (x; ω, a) + En m,im (ω,a)-En m-1 ,i m-1 (ω,a) v (n1,i1,...,n m-2 ,i m-2 ) n m-1 ,i m-1(x;ω,a)-v (n1,i1,...,n m-2 ,i m-2 ) nm,im (x;ω,a)

  =W (φn 1 ,i 1 ,...,φn m,im ,ψ k |x) W (φn 1 ,i 1 ,...,φn m,im |x)

	k	i1,...,nm,im)	(x; a) =	W φ	(n1,i1,...,n m-1 ,i m-1 ) nm,im φ (n1,i1,...,n m-1 ,i m-1 ) ,ψ (n1,i1,...,n m-1 ,i m-1 ) k nm,im (x;a)	|x
	φ (n1,i1,...,nm,im) k,i	(x; a) =	W φ	(n1,i1,...,n m-1 ,i m-1 ) nm,im φ (n1,i1,...,n m-1 ,i m-1 ) ,φ (n1,i1,...,n m-1 ,i m-1 ) k,i nm,im (x;a)	|x

XII. ACKNOWLEDGMENTS

I would like to thank A. Bérard, R. Milson and C. Quesne for stimulating exchanges and very interesting suggestions.

x -(α+q(i1,...,im-2,+,-))-1/2 .

The result Eq(94) is then still valid at the m th step. We still have to verify that this condition is compatible with the regularity of the preceding potentials in the chain. For that, note that at the m th step the explicit expression of the resulting potential and its associated eigenfunctions are a priori independent of the ordering chosen for the sequence . This means that we can consider several chains leading to the same extended potential. Consider the set (n 1 , i 1 , ..., n m-2 , i m-2 , n m-1 , i m-1 ) in which we have regrouped first all the n j associated to a positive i j , that is, we have re-arranged the chain as

, -, where n ′ 1 , ..., n ′ m-1 is a permutation of (n 1 , ..., n m-1 ) and where q + = q + (i 1 , ..., i m-1 ).

Suppose first that q + = m -1, that is, i ′ j = +, ∀j ∈ {1, ..., m -1}. ∀j ∈ {1, ..., m -1} we have 1 -q + i ′ 1 , ..., i ′ j ≤ 0 and since a > 1, the φ

are all nodeless implying the regularity of all the potentials V (n ′

and consequently the regularity

) . We can iterate this reasoning to go backward along the chain and we see clearly that at each step, the opposite of the "total charge" decreases, ensuring the regularity of the associated potential. In this procedure, once we attain the level q + , we have eliminated all the negative i ′ j and retrieve then the first case. Consequently, if we satisfy the condition a > 1 -q (i 1 , ..., i m-1 ) at the step m, we are ensured of the existence of a chain of regular extensions which leads to a regular V (n1,i1,,...,nm,im) .

IX. EIGENSTATES OF THE EXTENDED POTENTIALS

We can generalize the calculations made in the two-step case to determine explicit expressions for the eigenstates of the successive extensions. For that we adopt the ordering mentioned above and consider the spectrum of the potential V (n1,+,...,nq + ,+,nq + +1,-,....,nm,-) where q + = q + (i 1 , ..., i m ). Using Eq(56) and Eq(78), we have