Photonic bandgap engineering using second-order supersymmetry

First-order supersymmetry (SUSY) adapted from quantum physics to optics manipulates the transverse refractive index of guided-wave structures using a nodeless ground state to obtain intended modal content. Second-order SUSY can be implemented using excited states as a seed function, even with the presence of nodes. We apply second-order SUSY to the coupled-mode equations by recasting them as the Dirac equation. This enables the engineering of non-uniform surface corrugation of waveguide gratings and coupling potential, which encapsulates the Bragg interaction between counterpropagating modes. We show that the added bound states appear as transmission resonances inside the bandgap of the finite grating. The probability density of each state provides the longitudinal modal energy distribution in the waveguide grating. The smooth modal energy distribution of the states obtained by SUSY can mitigate longitudinal spatial hole burning in high power laser operation. We demonstrate that degenerate second-order SUSY allows the insertion of two states, which can coalesce into Friedrich-Wintgen type bound states in the continuum (BIC) for one-dimensional grating. We show that the eigenfunctions of BIC states are doubly degenerate with opposite parity, and the corresponding transmission resonances have phase changes of 2π across these states. One-dimensional BIC states can find application as robust high-speed all-optical temporal integrators by lifting restrictions on the length of various sections in the phase-shifted grating. Supersymmetric (SUSY) transformations, which originated in quantum physics offer a robust, physics-based approach to design photonic structures. The authors apply second-order SUSY to engineer non-uniform surface corrugation of waveguide gratings and coupling potential, and its application in photonics.

P hotonic integrated circuits (PICs) provide a robust platform for data communications 1,2 , quantum technologies [3][4][5] , and various sensing applications 6,7 . Optical waveguides coupled to a resonator are vital components in devices such as filters, adddrop (de-)multiplexers, and distributed feedback (DFB) lasers for PICs. These devices rely on engineering the geometry or arrangement of materials within the waveguide to achieve intended functionalities.
Inverse design 8,9 and machine learning 10 are at the forefront of the design techniques for modern photonic structures. These methods achieve the targeted optical response by optimizing parameters in the design space, but often fail to provide the insight behind physical processes in the optimized structure. On the contrary, supersymmetric (SUSY) transformations, which originated in quantum physics and were recently adapted to photonics, offer a robust, physics-based approach to design photonic structures 11 . In quantum mechanics, SUSY allows the design of a potential and its eigenfunctions simultaneously from a known potential 12 . The equivalence of Helmholtz and Schrödinger equations has allowed the application of first-order SUSY (1-SUSY) in optics to design the transverse refractive index profile of waveguides and their modal content 11,13 . Engineering the transverse refractive index using SUSY has enabled design of digital multimode devices 14 , mode sorter 15 , and lasers 16 .
The behavior of guided modes is tailored by engineering the refractive index or geometry of the waveguide in the direction of propagation. Recently, metasurfaces on top of waveguides have emerged as an attractive method to manipulate guided waves and tailor the coupling of guided waves to free-propagating waves 17,18 . The advancement in fabrication techniques for such hybrid structures has made it possible to realize complex geometries experimentally 19 . However, an intuitive physics-based approach to design such structures is currently missing. We demonstrate that the second-order SUSY (2-SUSY) offers an efficient approach to design waveguides with engineered surface corrugation to obtain desired functionalities.
Commonly used 1-SUSY transforms the potential V 0 of the Hamiltonian H 0 ¼ À∂ 2 z þ V 0 in its nonsingular partner potential V 1 whose spectra can differ at most in the ground state energy level. This type of transformation is known as unbroken SUSY and is defined by first-order intertwining operator A ¼ ∂ z þ WðzÞ, where W(z) is the superpotential. For an arbitrary solution u(z) to the initial Hamiltonian at the energy value ϵ such that ½À∂ 2 z þ V 0 uðzÞ ¼ ϵuðzÞ, the superpotential becomes W (z) = ∂ z [u(z)]. The corresponding partner potential then becomes V 1 ¼ V 0 À 2∂ 2 z ½ln u. It is evident that V 1 is nonsingular in a given region if the generating function u(z) does not vanish in this region. Thus, generating function either equal to ground states energy E 0 or lower ϵ ≤ E 0 can be used in the 1-SUSY transformation. In optics, each SUSY operation transforms the refractive index distribution (potential) and the propagating modes (eigenfunctions) using a nodeless fundamental mode (ground state). Broken SUSY produces a potential isospectral to the initial potential, which has been utilized to produce photonic configurations with identical reflection and transmission characteristics 20 and complex potentials with real eigenvalue spectrum 21 . Isospectrality of 1-SUSY has also been utilized to preserve bandgaps while transforming ordered potential to potentials analogous to Brownian motion 22 .
The main limitation of 1-SUSY is the inability to modify the excited part of the eigenspectrum for a nonsingular partner potential. 2-SUSY is implemented through second-order intertwining operator B, which is obtained by two generating functions u 1 and u 2 . The generating functions are solutions of the initial Hamiltonian ½À∂ 2 z þ V 0 u 1;2 ðzÞ ¼ ϵ 1;2 u 1;2 , which are not required to satisfy the boundary conditions. The relationship for ϵ 1 and ϵ 2 is used to classify the types of 2-SUSY. If the two energies are unequal and real ϵ 1 ≠ ϵ 2 2 R, then the 2-SUSY transformed potential is 23,24 . When the energies are equal and real ϵ 1 ¼ ϵ 2 2 R, then the partner potential is given by where w 0 is an arbitrary constant 25 . If ϵ 1 is complex, then for a nonsingular partner potential ϵ 2 ¼ ϵ Ã 1 and the partner potential becomes We have provided an extended analysis of 2-SUSY for potential in the Schrödinger equation in Supplementary Note 1. It is evident that the partner of 2-SUSY is obtained by the Wronskian of the two generating functions. Thus, the partner potential is nonsingular when the Wronskian is nodeless rather than the generating functions. We discuss the critical differences in the partner potentials generated by 1-SUSY and 2-SUSY for infinite square well in the section "Results".
Recently, the equivalence between Maxwell's and Dirac equation has been applied to understand the electromagnetic spin and orbital angular momentum 27 , and examine the relationship between interface states and topology 28 . The Dirac equation is a set of coupled first-order differential equations, and a matrix intertwining operator relates two Dirac equations. The matrix form of the intertwining operators form a second-order polynomial in Hamiltonian, which is a hallmark of second-order or nonlinear SUSY 29 . Here, we show that the 2-SUSY can design the optical response of any system described by the coupled-mode equations by rewriting them as the Dirac equation. We study the amplitude and phase of the transmission spectrum for the states added by 2-SUSY at prescribed detuning (eigenvalue) by modifying the coupling parameter (potential). We use the extra degree of freedom provided by 2-SUSY to insert two states at the same eigenvalue in one dimension with opposite parity. We show that the degenerate states in waveguide gratings coalesce to form bound states in the continuum (BIC) strictly in one dimension. We discuss the practical applications of gratings obtained using 2-SUSY transformation.

Results
Key differences between 1-SUSY and 2-SUSY. We illustrate the key differences between 1-SUSY and 2-SUSY through the example of the infinite square well potential. In Fig. 1a, we show the 1-SUSY partner of the infinite square well potential shown in Fig. 1b. In Fig. 1c, we present a 2-SUSY partner of the infinite square well where the ground state at the energy E 0 and the first excited state at the energy E 1 are deleted by using unequal and real factorization energy. 2-SUSY allows manipulation of two adjacent eigenstates, and new states can be inserted at any position between two states. In Fig. 1d, we have deleted the first excited state at the energy E 1 and inserted a new state at the energy value E = 3.2 (a.u.). In Fig. 1e, we have deleted the ground state at the energy E 0 and inserted a new state at the energy value E = 2 (a.u.). Thus, in both cases, we have effectively moved one state of the potential to an arbitrary position. This type of transformation is not achievable by 1-SUSY as it would produce an isospectral potential. In Fig. 1f, we have presented an isospectral potential obtained by 2-SUSY when the factorization energies are complex. In Supplementary Fig. 1, we have presented the 2-SUSY partner potential for harmonic oscillator potential.
Optical-quantum analogy. Light propagation and interaction in the waveguide with non-uniform surface corrugation or index variation is described by a set of first-order coupled differential equations. The physical mechanism of wave propagation in such media is called DFB or Bragg reflection, which results from the cumulative reflection from each surface of the structure. The physical interaction in a non-uniform structure is characterized by the coupling parameter varying longitudinally along the waveguide.
We design a one-dimensional DFB cavity formed by weak corrugation of the height of the waveguide, as shown in Fig. 2a. In the cavity, uniform periodic perturbation results in constant interaction between counter-propagating modes, which is characterized by constant coupling potential given by Fig. 2b. The coupling leads to formation of a bandgap in the transmission spectrum and the dispersion of the cavity as shown in Fig. 2c. The cavity can be modified to sustain n discrete modes by introducing nonuniformity in the surface perturbation, see Fig. 2e. Figure 2d shows the dispersion of unperturbed waveguides, forming the scattering (transmission) continuum for the localized modes. The field intensity in the nth mode is confined in the longitudinal direction along the cavity. Each mode decays independently into both waveguides 1 and 2 with the decay rate given by the sum of the two processes. The transmitted field intensity consists of the superposition of the input electromagnetic energy and the field originating from the decay of the localized states. The transmitted field intensity in waveguide 2 originates entirely from the decay of the localized states. Thus, for frequencies at which localized modes are supported, we observe a peak in the transmission spectrum.
The average height of the waveguide is h 0 . The corrugated region has a finite length L, and the profile is hðzÞ ¼ h 0 þ ΔhαðzÞ cos½2πz=Λ þ θðzÞ for 0 < z < L, where Λ is the grating period. The slow variation of the amplitude and the phase of the grating is described by α(z) and θ(z). The perturbation in the height of the waveguide introduces a coupling between the two counter-propagating modes at the same optical frequency ω, near the Bragg frequency ω B = c 0 π/n 0 Λ. The electric field of the propagating modes around a Bragg frequency (ω B ) is given by Eðz; tÞ ¼ ½uðzÞe ÀiðωtÀk B zÞ þ vðzÞe Àiðω B tþk B zÞ þ c:c:, where u(z) and v(z) are the envelopes of the two counter-propagating waves, k B is the Bragg wave number, and c.c. stands for complex conjugate. The slowly varying envelopes for a weak grating strength (Δh/h 0 < < 1) satisfy the following coupled-mode equations 30 : where δ = k − k B is detuning from Bragg frequency, the coupling is described using coupling constant of uniform grating q 0 by parameter as q(z) = q 0 α(z)e −iβ(z)31 . The coupling constant for a uniform grating is given by where h eff is the effective height of the corrugated region, n w is the refractive index of the waveguide, and n eff is the effective index of the planar waveguide. The coupling parameter q(z), which depends on the waveguide refractive index, and variation in the height of the grating forms the potential in the eigenvalue equation and detuning from the Bragg frequency is the eigenvalue. The size of the bandgap for a uniform coupling potential is given by (|δ| < q 0 ), which is identical to the gap between electron and positron states of Dirac particle.
Thus, SUSY formalism for the Dirac equation is equivalent to application of SUSY in the coupled-mode equation. SUSY acts as an inverse design method to obtain transformed coupling potential required for the desired transmission spectrum and modal energy distribution simultaneously. As the Dirac equation describes a relativistic particle, the potential is classified as a scalar, vector, and pseudoscalar according to the behavior under Lorentz transformation. We can define another unitary operator U ¼ 1 ffiffi 2 p ð1 þ σ 2 Þ to obtain a new Dirac Hamiltonian where the coupling parameter is pseudoscalar potential. Next, we describe the differences in the spectrum for the two types of potentials.
SUSY design of cavity. Improvements in experimental methods in realizing Dirac materials have led to enormous theoretical activity regarding exact solutions of Dirac-like equations. The bound state spectrum of the Dirac equation with scalar potential has received attention with interest in the topological phase of matter 35 . The position-dependent mass plays the role of a system parameter to obtain a topologically protected edge state, which are the bound states of a scalar potential at the zero crossings 36 . 2-SUSY modifies the coupling potential to support bound states and desired frequency in the bandgap by engineering a single defect region. The inhomogeneity introduced by the transformation makes the effective index real at that frequency, allowing a sustained cavity mode. In Fig. 3a, we show the coupling parameter as scalar potential, which is transformed to support two states in the positive and negative detuning. SUSY transformation of scalar potentials produces a symmetric transmission spectrum where the positive states are mirrored in the negative spectrum. The pseudoscalar interpretation of the coupling potential provides greater freedom for the insertion of states. Figure 3b shows the coupling potential and modal energy distribution for two mirrored bound states in addition to a state at zero detuning. The initial function can be modified to produce an asymmetric spectrum for the pseudoscalar case. In Fig. 3c, d, we present the grating design corresponding to the coupling potentials obtained in the scalar and pseudoscalar case. If the cavity is lossless, then the energy stored by the modes couples evanescently into the scattering (transmission) channel. The cavity loss for a corrugated waveguide is controlled by tuning the strength of perturbation or changing the length of the perturbed region. As field intensity distribution decays exponentially away from the defect region, the finite length of the corrugated region leads to coupling between confined energy to the transmission continuum. The leakage of energy combined with any distributed loss in the cavity creates resonances with finite spectral linewidth. We observe a phase change of π across each state, which is a typical characteristic of a resonance. The SUSY method's flexibility enables the design of frequency combs and the energy distribution for each cavity mode in one shot. Dirac-type equations have two associated Schrödinger equation. Thus, the coupled-mode-equations can be recast into the Helmholtz equation to apply first-order SUSY 37 , which restricts the variety of spectrum obtained.
For a cavity with gain, the modes inserted can be used for lasing. For a state inserted at zero detuning, the coupling obtained by SUSY reverses the sign at the center of the corrugated region, similar to a phase-shifted grating. However, the smooth change in the sign of the coupling contrasts the phase-shift layer, which acts as a point defect in the periodic structure that results in abruptness in the intra-cavity field distribution, as shown in Fig. 3e, f. The advantage of smooth modal energy distribution is evident for high-power operation where irregular field distribution in phase-shifted DFB grating leads to longitudinal spatial hole burning. Various approaches, such as multiple phase shifts and chirped grating periods, lower the abruptness, reducing the fabrication tolerances. The smooth variation of fields near the defect region allows the devices to have easier fabrication, especially additive manufacturing techniques. Figure 3g, h illustrates the difference between a waveguide grating obtained using SUSY to insert state at zero detuning and phase-shifted grating. SUSY requires one defect region for multiple states, which is less cumbersome to fabricate than multiple defects at precise locations for traditional methods.
Bound states in the continuum. The inverse method to obtain a potential, which supports BICs from the wavefunction result in weakly localized potentials that oscillate at infinity 38 . These potentials are unrealistic for application in photonic systems as they are highly sensitive to errors introduced during fabrication 39  Coupling potential and field distribution obtained by SUSY transformation. a The coupling potential (shaded gray) derived by the SUSY method supports four modes at frequencies in the stopband. The modal energy distribution, which is the probability density of eigenfunction, is shown in red. The inserted states appear as transmission peaks and the corresponding phase change of π. Phase wrapping from [−π, π] leads to the sharp features at some resonant states. The finite length of the grating truncates the potential, which leads to the finite linewidth of each state. b An equivalent representation of the coupled-mode equation, where the coupling is pseudoscalar potential. It enables the design of a grating that supports states at the center of the stopband. c, d The design of grating with an envelope was obtained using SUSY transformation in a and b, respectively. The grating period has been magnified for visualization. e Pseudoscalar potentials produce grating with adiabatic taper (g), leading to a phase shift for the zero detuning state. The modal energy distribution is smooth in the defect region. f In contrast, modal energy distribution in phase-shifted grating has an abrupt change (h), which leads to longitudinal hole burning in high-power applications.  44 .
Here, we show the inverse method based on SUSY enables the construction of potentials with two eigenvalues separated by an infinitesimal parameter ϵ. The analytical nature of the formalism prevents any numerical error that might occur in the analysis of these eigenvalues using a computational intensive inverse method. As this parameter is set to zero, the eigenvalues coalesce, creating degenerate resonances in the scattering spectrum of the optical structure. Similar method is used to create potential with degenerate eigenvalues in Schrödinger equation on half-axis 45 and full axis 46 . However, these potentials are oscillatory and have singularity on the full axis, creating hurdles for realizing the photonic system.
Similarly, degenerate SUSY is applied to the coupled-mode equation, which creates overlapping resonances. Under the degenerate limit, SUSY adds two states at δ 1 and δ 1 + ϵ, where ϵ → 0. The coupling potential is obtained using Taylor expansions of second seed spinor at δ 1 + ϵ to avoid indeterminate solution as shown in the section "Methods". For a waveguide, finite-length corrugation connected to input-output waveguides enables interference of the resonances coupled to the same continuum. The interference of two degenerate states as shown in Fig. 4a produces a BIC, where the transmission coefficient is real, and the phase change across the state is zero. Figure 4a shows the coupling potential obtained after degenerate SUSY transformation in coupled-mode equations. The amplitude of the transmission coefficient at the BIC frequency is unity, and the corresponding phase change is 2π. Effectively the phase does not change across the BIC state in the DFB grating cavity. The finite length grating shown in Fig. 4b leads to the apodization of the potential. Thus, the transmission peak corresponding to the BIC has a finite width, allowing the coupling of the light in these states. In the ideal case, without the grating apodization, the lifetime of the BIC state will be infinite and zero phase difference across the state.
In practice, BIC states in photonic structures appear extremely narrow resonances which have found applications for lossless propagating PICs 47 and laser design 48 . BICs in strictly onedimensional structures have a promising application as photonic temporal integrators 49 . According to signal processing theory, a temporal integrator is implemented using a linear filtering device with a unit step temporal impulse response. In electronics, a capacitor is used for temporal integration. The electric charge accumulates at the capacitor, and the integrated signal is proportional to the voltage measured. In practice, the voltage at the capacitor decays exponentially. This property is used to design an optical temporal integrator using phase-shifted grating 49 . However, the length of both sections of phase-shifted grating must be equal, which restricts the central frequency of the integrator to the zero detuning. As the BIC state is a result of interference, the confinement of light does not depend on the location of the defect but the asymptotic behavior. Figure 4c compares the temporal response of the integrator formed by zero detunings and the BIC state. The BIC integrator provides a longer integration time window, which leads to higher processing speed 50 for an all-optical temporal integrator.

Discussion
SUSY transformation to design of quasi-isospectral refractive index landscape allows for global phase matching critical to wave devices. First-order SUSY is limited to the manipulation of the ground state of the potential. Second-order SUSY enables the design of the excited part of the spectrum and arbitrary control of insertion of new states. The extra degree of freedom provided by second-order SUSY can pave the path for a new class of SUSY  Fig. 4 Bound states in the continuum. a Coupling potential obtained by degenerate SUSY for two degenerate states shown in black and blue. Both states decay into the same channel, two states at the same frequency can interfere, leading to Friedrich-Wintgen BIC. For a finite grating, the transmission amplitude for BIC states is one, and the phase shift across the state sweeps 2π. In the ideal case for an infinite grating, the phase does not change. b The structure of grating for coupling potential supporting degenerate states. The increase in height is still within the weak grating limit. The grating period has been magnified for visualization. c Comparison of the temporal response of zero detuning state and the BIC state. The ideal temporal integrator has a unit step impulse response. However, the leakage of energy reduces the temporal bandwidth and imposes a strict restriction on the different sections of the grating. As the localization of fields in the BIC state occurs due to interference, the BIC integrators do not have length restrictions. designed optical structures. We have applied second-order SUSY to coupled-mode equations by rewriting them into Dirac-form to design non-uniform corrugation based on the desired spectral response. Waveguide-resonator systems underpin PICs that provide a robust platform for fast optical information processing. The grating-based device is engineered by extra degrees of freedom provided by slowly varying nonuniform gratings, where the local period or depth of modulation varies longitudinally. We have investigated the amplitude and phase characteristics of the transmission spectrum for various cavity engineered from its coupling parameter. The SUSY method provides an elegant single-step procedure to simultaneously derive the structure of the grating and the modal energy distribution as compared to the synthesis of waveguide grating from the coupling parameter using the layer peeling method 51 . We observe that the bound states of the coupling potential appear as the resonant state for a finite photonic cavity. The phase of transmission resonances changes by π over the small frequency range. The lifetime can be tuned by engineering the strength of the waveguide for optical filters.
Second-order SUSY allows manipulation of excited states and insertion of a new state at an arbitrary position between two adjacent energy levels. This transformation cannot be achieved by first-order SUSY as the insertion of new states lead to isospectral potentials. This property of the second-order SUSY method can be utilized to add overlapping resonances in the grating. We observe that the transmission coefficient in the finite cavity is real, and the phase changes abruptly by 2π across the BIC state, which contrasts the infinite case where the phase does not change. The high transmission over a narrow spectrum allows the realization of all-optical temporal integrator.
Coupled-mode equations have been generalized to include higher-order effects, which become essential in deeper gratings 52 . We have demonstrated the method by using shallow corrugation on the grating surface, but the formalism used to obtain the profile is very generic. Hence the method described can be adapted for designing a strongly interacting photonic crystal. Additionally, coupled first-order differential equations describe the evolution of the modal amplitudes of two states, which evolve either in time or with propagation distance 53 . Thus, the SUSY method described can be used to obtain the new design for a large number of photonic systems. SUSY procedure of BIC can be repeated to design any number of states with degenerate eigenvalues to create resonance with exceptionally low spectral linewidth. The (1 + 1) Dirac equation also appears in the Ablowitz-Kaup-Newell-Segur (AKNS) hierarchy for the Lax eigenvalue equation of the nonlinear evolution equations 54 . The formalism of degenerate Darboux transformation can be also used to create FW BICs in one dimension is also used to study the generation of higher-order rogue waves 55 .

Methods
Theoretical background. One-dimensional (1 + 1) Dirac equation with a generic (2 × 2) matrix potential (V 0 ) with well-defined eigenvalues (E) and spinors (Φ) for particle of mass (m) where S(z) is the Lorentz scalar potential, and P(z) is the pseudoscalar potential. In coupled-mode equations the potential V 0 is replaced by coupling potential q 0 . We define an intertwining operator ðLÞ relating two Dirac Hamiltonians (h 0 = −iσ 2 ∂ z + q 0 (z)) and (h 1 = −iσ 2 ∂ z + q 1 (z)), such that Lh 0 ¼ h 1 L. In the simplest form, the intertwining operator can be be assumed to of form L ¼ A∂ z þ B, where A and B are two (2 × 2) z-dependent matrices. Upon substitution, in the intertwining relation, we obtain three relationships. First is a restriction on the matrix elements A, where only two elements can be arbitrarily chosen, and others need to be derived using Aσ 2 − σ 2 A = 0. Second, a difference of matrix potentials Similar to Schrödinger case, Riccati equation is linearized upon substituting matrix G = −∂ z UU −1 , and we get an equation in U, This equation is similar to the Dirac equation where the eigenpair is a spinor and a number. However, in this case the eigenpair U and Λ are matrices. The solution of Eq. (4) is obtained by choosing U is composed of two spinors (Φ 1 , Φ 2 ) which are solution of the Dirac equation corresponding to eigenvalues, δ 1 and δ 2 , respectively, h 0 Φ 1,2 = δ 1,2 Φ 1,2 .
In principle, the eigenstates must be formal solutions regardless of boundary and normalization conditions.
The Darboux transformation is then defined with matrix A is identity as L ¼ 1∂ z À UU À1 and the eigenfunction is transformed to ψ ¼ LΦ.
Repeated transformations. In order to illustrate the steps required to repeat the SUSY transformation, the eigenfunction is rewritten as . We can substitute the where W(Φ 1 , Φ 2 ) = Φ 11 Φ 22 − Φ 12 Φ 21 is the Wronskian. Using the identities of matrix algebra, we can obtain the transformed eigenfunction Similarly, the potential is derived by Thus, multiple (n) SUSY transformation produces the potential q 1 = q 0 + D n σ 2 − σ 2 D n , where the four determinants are constructed from W(Φ 11 , Φ 21 , ⋯, Φ n1 ) by replacing the 2nth line with the nth derivatives of the first element of the spinors and (2n − 1)th line with the nth derivatives of the second element of the spinors.
Second-order SUSY in Dirac equation. Dirac-type Hamiltonian is by definition a first-order differential equation and the intertwinning operator and its adjoint form a second-order matrix differential operator LL y and L y L. Upon substitution of the form of operator L and L y , it is shown that Thus, the superalgebra is given by 4 × 4 matrix Hamiltonian and supercharges, respectively, fQ; Q y g ¼ ðH À ϵ 1 ÞðH À ϵ 2 Þ; fQ; Qg ¼ fQ y ; Q y g ¼ 0: Thus, the SUSY formalism is associated with the Dirac equation of the second order or nonlinear form.
Degenerate SUSY transformation. Application of SUSY twice for the same eigenvalue does not produce two states. In addition, the potential obtained after two SUSY transformations for the same eigenvalue has an indeterminate form. However, the issue can be resolved by applying the SUSY of two eigenvlaues for δ 1 and δ 1 + ϵ. As the parameter ϵ is infinitesimal, we can use Taylor expansion of seed eigenfunctions and elements of the determinant for DT. For an infinitesimal ϵ, the eigenfunction Φ 2 is approximated by the Taylor expansion Φ 2 (δ 1 + ϵ) = Φ 2 (δ 1 ) + [∂ ϵ Φ 2 (δ 1 + ϵ)]ϵ + O(ϵ). Similarly, potential requires an element-wise Taylor expansion of the 2n × 2n matrix D n give by ∂ n i ϵ D n ij ðδ 1 þ ϵÞ, where n i is the floor of the function [(i + 1)/2], and [i]. After transformation has been applied to insert states at δ 1 and δ 1 + ϵ, we let ϵ → 0 in the expression for potential and eigenfunctions.

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.