Flat-band engineering of quasi-one-dimensional systems via supersymmetric transformations

We introduce a systematic method to spectrally design quasi-one-dimensional crystal models described by the Dirac equation in the low-energy regime. The method is based on the supersymmetric transformation applied to an initially known pseudo-spin-1/2 model. This allows extending the corresponding susy partner so that the new model describes a pseudo-spin-1 system. The spectral design allows the introduction of a flat-band and discrete energies at will into the new model. The results are illustrated in two examples where the Su-Schriefer-Heeger chain is locally converted into a stub lattice

Both experimental [10] and theoretical [26,27] efforts have been made to provide useful tools and methods for engineering flat-band systems.For instance, via repetition of microarrays [26], utilizing polynomials of tight-binding Hamiltonian [27], through compact localized states classification [28], and using graph theory [29,30].The latter frequently rely on the tight-binding approach.In this article, we resolve to the Dirac approximation valid for low-energy systems, where the dynamics is described by Dirac equation with pseudo-spin one.Our approach is based on supersymmetric quantum mechanics so that the existence of a flat-band is granted by construction.Supersymmetric transformations, equivalently known in the literature as Darboux transformations, are specific non-unitary mappings between evolution equations of two quantum systems.
The latter can be used for the construction of the new solvable models, where the potential term of the initial system gets deformed, yet, the knowledge of the solutions is preserved.These transformations have been broadly explored in non-relativistic quantum mechanics [31], and during the last years in the construction of solvable models described by one-or two-dimensional Dirac equations.Supersymmetric transformations for low-dimensional Dirac operators was discussed in [32,33], and employed in a series of works, see e.g.[34][35][36][37][38][39][40][41][42][43][44].Most of these works focus on the analysis of pseudo-spin−1/2 quantum systems, but it was recently applied in the context of pseudo-spin-1 flat-band systems in [45].
The Su-Schriefer-Heeger (SSH) model is a one-dimensional chain of dimerized atoms, used originally for the analysis of solitonic effects in macromolecules [46,47], and known for possessing non-trivial topological properties [48].The low-energy approximation of its tight-binding Hamiltonian corresponds to the one-dimensional Dirac operator [47].Interestingly, solitonic states emerge in SSH ladders on domain walls, where the dimerization of atoms gets inverted.Domain walls on SSH-type chains of coupled dimers were experimentally realized on chlorine vacancies in the c(2×2) adsorption layer on Cu(100) in [13].The existence of topological domain wall states was discussed in [49,50], whereas the supersymmetric transformation has been applied to induce a topological gapped state in the SSH chain [51,52].Furthermore, the transmission properties of pseudo-spin-1 Dirac equations described through decorated, "bearded," SSH chains have been discussed [53].Spectral and symmetry properties of trimer SSH chain with next-nearest-interaction where considered in [54].
In this article, the supersymmetric transformation is exploited to connect known pseudo-spin-1/2 quantum models with new unknown pseudo-spin-1 models.Particularly, the transformation allows tuning the emerging flat band of the new model while adding new bound state energies, assuming that the proper boundary conditions are met.To this end, a pseudo-spin-1/2 model is trivially extended into a pseudo-spin-1 system by adding an isolated coupling term, such that the dispersion relations are kept invariant.The supersymmetric transformation of such an energy operator is then matched with a new and non-trivially extended pseudo-spin-1 Hamiltonian, where the added coupled term is no longer isolated and now describes an interaction with the rest of the elements in the system.Particularly, a graphene-like system is used as the initial model in the transformation, leading to explicit models that behave asymptotically as an SSH chain with altered dimerization patterns resembling the domain wall.Such an SSH chain gets locally decorated by additional atoms, forming a stub lattice in the localized region.This allows for a systematic mechanism to spectrally manufacture pseudo-spin-1 models based on relatively simple pseudo-spin-1/2 counterparts.
The manuscript is structured as follows.Section 2 summarizes the periodic structure and dispersion bands of the generalized stub lattice, where the special case where a flat-band emerges is considered.In Section 3, the general framework of the Darboux transform (susy transform) for arbitrary pseudo-spin systems is briefly introduced.Here, the transformation is implemented for an extended pseudo-spin-1/2 so that the susy partner renders a non-trivial pseudo-spin-1 model.Applications of the latter are further exemplified in Section 4 and Section 5, where explicit cases of quasi-one-dimensional pseudo-spin-1 systems are derived and discussed.Further discussions and future perspectives for future applications of the present results are detailed in Section 6.

Generalized stub lattice
Let us consider the tight-biding model of the generalized stub lattice, where the hopping amplitudes are considered real but general otherwise.The model here is such that it converges to the SSH lattice or stub lattice for specific choices of the hopping parameters.The tight-binding Hamiltonian H gs of a generalized stub lattice (see Fig. 1a) is where A † n , B † n , and C † n are fermionic creation operators for electrons on site A n B n or C n .The lattice is supposed to be infinite, i.e. n acquires all integer values.The quantities t 1 , t 2 and t 3 are the real hopping amplitudes, and µ A,B,C correspond to real on-site energies.The index n runs over all elementary cells through an infinite periodic lattice.The primitive translation vector is (a, 0) so that the first Brillouin zone becomes k ∈ 0, 2π a .The Fourier transform of H1 provides us with the following operator: together with the secular equation det(H gs (k) − λ) = 0, which reads as Although the latter can be solved for any λ (for more details, see [55]), we are interested in the configuration where the flat-band is present.Indeed, this occurs for which, henceforth, is the case under consideration.This leads to the dispersion relations of the form From the latter, it is worth remarking that dispersion bands never touch; i.e., the band structure is always gapped.This is a stark difference with respect to other flat band systems such as the two-dimensional Lieb lattice [55], where vanishing next-nearest neighbors close the band gap.Furthermore, the trimer SSH model [16] does not hold a flat band in its spectrum.
Our model can be reduced to a couple of special cases when the parameters are fixed correspondingly.That is, for t 3 = 0, the C atoms are effectively isolated from the linear chain formed by A and B atoms.The C atoms host the flat-band states with energy µ C .These states are strongly localized at the C atoms as there is no interaction with other atoms.In turn, for t 1 ̸ = t 2 , the linear chain of A − B atoms coincides with an infinite SSH model.The case t 1 = t 2 = t 3 reproduces the stub lattice.The band structure of the system is illustrated in Fig. 1b.
The energy band E + (E − ) has its minimum (maximum) at K = π a .Expanding E + around this point, k = K + δk, we get When the momentum k is considered in the range where the quartic and higher terms in the expansion are negligible a 2n δk 2n ∼ 0, n ≥ 2, the dispersion relation turns into the expression known for massive one-dimensional Dirac fermions.By expanding the Hamiltonian (2) at K up to the first order in δk, we get the Dirac-type operator for which acts, in general, on three-component wave functions Ψ ∈ C 3 .For the sake of simplicity, it is more convenient to make the additional transformation (ψ 1 , ψ 2 , ψ 3 ) → (ψ 1 , i ψ 2 , i ψ 3 ).The resulting operator in coordinate representation reads as There are situations where Dirac equation does not provide reasonable approximation of lowenergy dynamics.It can happen when the energy gap is so large that the admissible energies are already out of the range where linear approximation of dispersion would be faithful.In our models, we assume that there can be reached a reasonable control over the parameters, e.g by performing the experiments on the optical lattices, that makes it possible to stay within the range of energies where Dirac approximation works well.
In the next section, we will present the method that allows to construct (7) with possibly inhomogeneous hopping t 1 and t 3 and facilite calculation of the associated eigenstates.

Coupling via Darboux transformation
Let us start the section with a brief review of the Darboux transformation for Dirac-type operators of the form H = −iγ∂ x + V , where γ and V can be generic N × N matrices.The Darboux transformation for N = 2 was discussed in [32], while the general case was considered in [33].In general, the Darboux transformation relates the initially known stationary equation (H − ϵ)Ψ = 0 with the new unknown equation ( H − ϵ) Ψ = 0, where H = −iγ∂ x + V is also a Dirac-type operator with an altered potential term.Furthermore, the transformation maps the solutions of the first equation into the solutions of the second equation.The latter is not necessarily a one-to-one mapping, and is achieved through a first-order and non-unitary differential operator L, the exact form of which is shown below.
The transformation is based on N eigenstates Ψ a , a = 1, 2, . . ., N , of H, (H − λ a )Ψ a = 0.The eigenstates are used to compose an N × N matrix U = (Ψ 1 Ψ 2 . . .Ψ N ).There holds so that we can define the new Dirac-type operator with U x ≡ ∂ x U .The latter is related with H through the intertwining relation LH = HL, where L is the first-order differential operator with P x ≡ −i∂ x the momentum operator.Here, the operator L effectively maps the eigenstates of H into the eigenstates of H, with the exception of the states Ψ a , a = 1, . . ., N that belong to the kernel of L. Indeed, there holds The Hamiltonian H can have the new bound states of energies The columns of (U † ) −1 correspond to formal eigenstates of H.When j − th column of the later matrix is square integrable, then it forms the bound state of H with energy λ j , see [32].
Let us consider a generic, one-dimensional pseudo-spin-1/2 Dirac system described by the following stationary equation where m = m(x), v = v(x) and A = A(x) are real functions so that H 1/2 is hermitian.We assume that it is possible to find formal solutions of the equation for any real ϵ.
We trivially extend H 1/2 by an additional degree of freedom so that the new operator have the form This represents a system where two subsystems coexist without any mutual interaction.In one of them, dynamics is driven by H 1/2 while in the second one, dynamics is frozen as the energy operator is constant.
It is straighforward to find the eigenvectors of the extended operator H 1 from the eigenvectors of H 1/2 .We shall use them to perform the supersymmetric (susy) transformation of H 1 .In order to do so, we fix the matrix U , see (8), in the following manner The columns of U are formed by the eigenvectors corresponding to the eigenvalues ϵ or λ, respectively.The components ψ a , ϕ a , a = 1, 2, 3 and ξ 1 and ξ 2 can be fixed as real-valued functions.The functions ξ 1 and ξ 2 can be arbitrary, but they should not be zero identically as the transformed Hamiltonian with coupled subsystems could not be hermitian in that case, see Appendix.
With the matrix U fixed, we can construct L and H 1 through ( 10) and ( 13), such that the intertwining relation is satisfied.The new potential V 1 is not hermitian in general.Nevertheless, we can exploit the freedom in the choice of the functions ξ 1 and ξ 2 in order to recover the hermiticity of V 1 in (13).To this end, it is sufficient to fix ξ 1 in the following manner: with c 1 a real integration constant.It is worth noticing that W 0 ≡ ϕ 2 ψ 1 −ϕ 1 ψ 2 is a real constant as well.Indeed, the relation ∂ x W 0 = 0 can be derived when taking into account that (ψ 1 , ϕ 1 ) t and (ψ 2 , ϕ 2 ) t are eigenvectors of H 1/2 corresponding to the same eigenvalue.
The new Hamitonian H 1 defined in ( 13) has the following form where All the nonvanishing components (17 Such a situation would be undesirable as it would be necessary to introduce additional boundary conditions at the singularities.The additional boundary conditions could compromise the calculation of physically relevant eigenstates of H 1 .Indeed, physical eigenstates of H 1 could be mapped into the formal eigenstates of H 1 that would not belong to its domain.Therefore, the elements of the matrix U should be set such that d is a node-less function. The components U are not independent.Indeed, ξ 1 is given in terms of ξ 2 , see (15).The functions ψ a can be expressed in terms of ϕ a , a = 0, 1, 2, respectivelly, Additionally, ϕ 1 can be expressed via ϕ 2 as they are two linearly independent solutions of Therefore, d ≡ d(x) is determined by three functions only, d = d(ϕ 0 , ϕ 2 , ξ 2 ) where ξ 2 is arbitrary in principle.The freedom in its choice can be exploited to keep V 1 free of any additional singularities.We discuss the explicit choice of ξ 2 in the models presented in the next section.
The formulas (17) suggest that a major simplification of the potential V 1 occurs when either ψ 0 = 0 or ϕ 0 = 0.Then, there holds V 13 = 0 or V 23 = 0, respectively.We are interested in the latter as V 1 acquires the form of the Dirac operator (7) in this case.Although it is not possible to set ϕ 0 = 0 for a generic Hamiltonian H 1 , it is possible for cases where H 1 acquires the specific form where m is a real constant.When we fix ϵ = m, then we can find the corresponding eigenstate (ψ 0 , ϕ 0 , 0) where ϕ 0 = 0 and ψ 0 = exp A(x)dx .The matrix U then reads as The components of the simplified potential term are as follows The equation ( 19) reduces to Schrödinger-type equation.Assuming that ϕ 2 is fixed, we get ϕ 1 as follows, After substituting ( 15) and ( 24) into ( 23), we obtain the potential components Comparing the potential terms in H 1 with h gs in (7), we find that the two operators coincide provided that The hopping amplitudes t 1 and t 3 in the effective Hamiltonian H 1 of the quasi-one dimensional chain would be inhomogeneous.In this context, the operator H 1 described two systems without any mutual interaction.In contrast, the operator H 1 corresponds to a qualitatively different physical reality; the two subsystems are coupled by V 13 .In the next section, we will apply the presented framework for construction of two explicit models that can be matched with a decorated SSH model.We will discuss three explicit models where the inhomogeneity makes it possible to convert SSH chain to stub lattice locally.

Tunable flat-band in the gap
Let us fix the Hamiltonian H 1/2 in (12) as the energy operator of a massive particle with pseudospin-1/2 under the influence of a null external magnetic field with the gauge rule A(x) = A 0 ∈ R.
The trivially extended operator H 1 then reads as where m and λ are real constants, and the corresponding eigenvectors can be found for any ϵ.
In accordance with the results of the previous section, we fix ϵ = m and As mentioned in ( 18) and ( 19), the components ψ 1,2 and ϕ 1 can be obtained in terms of ϕ 2 .We will assume that |λ| ̸ = |m|.Then we can write where c 0 is a real constant.We used the fact that ϕ 1 and ϕ 2 have to solve the same differential equation ( 19) of the second order.We assume that they are linearly independent, i.e. the Wronskian W 0 of the two solutions is nonvanishing, W 0 ̸ = 0.The component ξ 1 is fixed as in (15).
The functions ϕ 2 and ξ 2 are to be selected such that the components (17) of V 1 are free of singularities.We make the following choice, Here we assume that A, m, and λ are fixed such that κ is real.The matrix U than satisfies (21).
The Hamiltonian H 1 and the intertwining operator L have the following explicit forms where Here we combined c 0 , c 1 and W 0 into a single parameter ω, It can be concluded from (35) that both V 12 and V 13 are regular provided that we fix ω such that For |ω| > ω crit , the interaction V 13 vanishes asymptotically.The potential term V 12 is asymptotically constant but it changes its sign, For most of the eligible values of either ω or λ, the term V 13 (x) represents a rather narrow well or a bump, dependently on the sign of ω, whereas V 12 (x) forms a smoothed potential step.When |ω| approaches ω crit , the magnitude of V 13 increases.Simultaneously, it gets wider so that it resembles a smoothed rectangular well (for ω > 0) or barrier (for ω < 0).The potential V 12 turns into a smoothed two-step barrier with an intermediate plateau.The width of the plateau is very sensitive to the proximity of |ω| to ω crit , see Fig. 2 for illustration.
When ω = ω crit , V 13 simplifies considerably.We have Asymptotic behavior of the interactions is different in this case.We have We can see that V 12 is changing its sign asymptotically again.The interaction V 13 acquires non-zero constant value for large negative x.
The potential V 1 can be matched with the interaction term of the quasi-one-dimensional chain (27) where the generalized stub lattice gets converted into SSH chain with a parallel chain of non-interacting atoms.The interaction between the two chains is localized in case of (35) while in case of (39), it gets extended over the half-axis.The dimerization pattern on the SSH chain changes as the ratio of t 1 /t 2 (where t 1 = V 12 + t 2 ) is inverted along the x axis.This way, it resembles the SSH chains with a domain wall.Fig. 3 illustrates the case ω = ω crit , where a semi-infinite generalized stub lattice decomposes at the origin into an SSH chain with parallel non-interacting atoms.
As noted below (11), supersymmetric transformation can generate bound states with discrete energies in the new system.The candidates for the new bound states are formed by the columns of the matrix (U † ) −1 that satisfies The eigenstate corresponding to the eigenvalue m is not normalizable, while the other two columns are eigenvectors with the eigenvalue λ.Their explicit form and square integrability is not of our interest.The reason is that λ corresponds to the flat-band energy and, therefore, it is infinitely degenerated anyway.Indeed, we can find infinite number of independent normalizable eigenvectors of the form L(0, 0, ξ(x)) T where ||ξ|| < ∞, with T the transposition operation.
By construction, the spectrum of H 1 is composed by two energy bands of negative and positive energies.The energy λ of the flat band can take any value within the energy gap, The barrier represented by V 1 in ( 35) is perfectly transparent.The quasi-particles tunnel through it without being back scattered.It it reminiscent to the Klein tunneling of quasiparticles in graphene through electrostatic barriers.Here, the particles can pass through the  Figure 3: Inhomogeneous hopping parameters t 1 = t 2 + V 12 (red), t 3 = V 13 (orange) and t 3 (gray-dotted) for ω = ω c computed from (37), (39), and (28).We fixed , and ρ = 1 (left) and ρ = 0.5 (right).In the insets, we depict the quasi-onedimensional chain with the corresponding interactions (the thicker the line, the stronger the coupling), with t 1 red, t 2 black, and t 3 the vertical line.barrier without reflection independently on their energy.It can be understood with the use of the intertwining operator L that makes it possible to map the eigenstates of H 1 into those of H 1 .The Hamiltonian H 1 corresponds to the free-particle energy operator with a constant potential.Let us suppose that its physical eigenstate ψ k (x) = e ikx (a, b, c) t corresponds to the plane waves with a fixed momentum k ∈ R. The intertwining operator L converts these states into the scattering states of H 1 .We can write where a, b, c are complex-valued constant components of the three-component wave function.
The matrix U x U −1 converges to a constant matrix for large |x|.Therefore, Lψ k (x), |x| → ∞, acquires the form of the plane wave whose momentum is not altered by the potential barrier, i.e. there is no back-scattering.

Coexistence of discrete and flat-band energy levels
The Hamiltonian H 1 can inherit spectral properties of the initial, uncoupled operator H 1 .The latter one, by construction, shares the spectrum of H 1/2 except the flat-band energy.Therefore, spectral design of H 1 starts with the proper choice of H 1/2 .In this section, there will be presented the models where the flat-band coexists with discrete energies.The models with one and two discrete energies will be introduced.We will illustrate how Darboux transformation can be used in two steps.In the first one, 2 × 2 Darboux transformation can generate H 1/2 with requested structure of discrete energies.In the second step, Darboux transformation applied on 3 × 3 operator H 1 produces the Hamiltonian H 1 .

Two bound states and a flat-band
In this subsection, we design a solvable model of a pseudo-spin-1 system with two discrete energies and a flat band in the energy spectrum.We will discuss in detail the role of the Darboux transformation at different stages of the construction.
First, we shall construct a solvable model described by pseudo-spin−1/2 Hamiltonian that has two discrete energies in its spectrum.Let us set the initial Hamiltonian as the energy operator of a free particle system as The spectrum of H 1/2 consists of two bands divided by energy gap that stretches between ±m 0 .Darboux transformation can be used to convert H 1/2 into the new pseudo-spin-1/2 Hamiltonian that would possess two discrete energies.We can rely here on the existing results.Solvable systems constructed from the free-particle model via Darboux transformation were discussed in [41,43] where models with diverse configurations of discrete energies within the energy gap were presented.
We demand that there are two discrete real energies λ 0 and −λ 0 in the new system.In order to construct such a Hamiltonian, we fix the matrix U 1/2 in the following manner, see [41] for more details, where u 11 = cosh k 0 x, u 21 = cosh(k 0 x + a 0 ), and The intertwining operator L 1/2 and the new Hamiltonian H 1/2 can be written as They satisfy The two linearly independent eigenstates of H 1/2 corresponding to the eigenvalue λ are They satisfy The eigenstates F (λ) and G(λ) of H 1/2 for an eigenenergy λ ̸ = λ 0 can be found with help of the intertwining operator L 1/2 , The Hamiltonian H 1/2 has two square integrable bound states v ± with energy ±λ 0 that form the columns of the matrix (U † 1/2 ) −1 , see ( 42), Therefore, we constructed the operator H 1/2 with the two discrete energies in the spectrum.Now, we use the operator H 1/2 to define the extended operator H 1 , see (20), Darboux transformation L of the extended system is defined in terms of 3 × 3 matrix U , see (10) and (21).We fix the components of U in terms of solutions F (λ) and G(λ) for a given λ, |λ| < λ 0 , in the following manner, Then it satisfies The Hamiltonian H 1 can be constructed as in (9).It inherits the discrete energies ±λ 0 of H 1 .The corresponding bound states w ± can be found with help of the intertwining operator L, The parameters ρ and κ of the flat-band solution can be arbitrary in principle.Substituting (58) into ( 23), the elements V 12 and V 13 of the potential of the new Hamiltonian H 1 are represented by rather extensive formulas that we will not present here explicitly.The numerical tests reveal the range of ρ and κ where V 12 and V 13 are regular, see Fig. 4a).
Nevertheless, the model gets remarkably simplified when ρ and κ are fixed as follows, Then the intertwining operators H 1 and L acquire particularly simple form where and A(x) is defined in (49).We can see that the dimerization pattern of the SSH chain undergoes the change at x = −a 0 /k 0 .The SSH chain gets coupled with the parallel chain of atoms by V 13 that acquires a constant value.The components V 12 and V 13 are plotted in Fig. 4b).The atomic chain with the corresponding hopping parameters is illustrated in Fig. 5.
With the current choice (60) of κ and ρ, the bound states (59) are where The other eigenstates of H 1 for eigenvalues λ can be found from those of H 1/2 , see (55), Discussion of the scattering properties of the model with the potential (75) can be conducted in close analogy with the previous model.The matrix U x U −1 in the operator L = ∂ x − U x U −1 tends asymptotically to a constant matrix, and, therefore, the operator L cannot change the momentum of the plane wave that corresponds to the scattering state of H 1 .Therefore, the current setting described by H 1 is also free of back-scattering.
In the construction, we used Darboux transformation at two different occasions.First, it was used to derive pseudo-spin-1/2 Hamiltonian H 1/2 with the requested discrete energies.In that case, the intertwining operator L 1/2 was represented by 2 × 2 matrix operator (49).Then we applied another Darboux transformation given in terms of 3 × 3 operator L. In the specific case (60), it acquired a compact form (61).It provided us with the Hamiltonian H 1 .It is worth noticing in this context that the two intertwining relations mediated by L 1/2 and L can brought into a compatible form by an extension of the operator L 1/2 .Indeed, the intertwining relation (50) can be written as Then it is possible to write down a single intertwining relation that connects the trivially extended initial Hamiltonian H 1/2 0 0 λ and the target Hamiltonian H 1 ,

Flat-band and a single bound state
In this case, we shall consider the model with a single discrete energy and a flat-band energy.
Following the strategy explained in the previous subsection, we shall select the Hamiltonian H 1/2 such that it has one discrete energy in the energy gap.We shall use here the results of [43,56] where such an operator was constructed via Darboux transformation and possessed the vector potential A(x) = m tanh mx.The trivially extended operator H 1 then reads as where M , m and λ are real constants.The stationary equation H 1 Ψ = ϵ Ψ is exactly solvable for any ϵ.Indeed, fixing the wave function Ψ = (ψ, ϕ, 0) and |ϵ| ̸ = M , the stationary equation decouples into ϕ = −i ∂x−A E+M ψ and (−∂2 x − ϵ 2 + M 2 + m 2 )ψ = 0. Hence, the upper component is the eigenstate of the Schrödinger Hamiltonian of the free particle.This is due to the fact that the potential term in (69) is related to the reflectionless Pöschl-Teller model 2 , see [43,56].The operator H 1 has a square integrable bound state Ψ −M with energy −M , The elements of the matrix U are fixed in the following manner: The parameter µ controls reality of these functions, i.e. µ = 1 for k 2 < m 2 and µ = −i for k 2 > m 2 .The components ϕ 1 and ϕ 2 can be calculated from ϕ a = −i ∂x−A E+M ψ a , with a = 1, 2. Furthermore, we fix the flat-band energy as The matrix U satisfies the relation so that Eq. ( 17) renders the potential components where We shall fix the parameters such that d(x) is non-vanishing for x ∈ R in order to keep the potential regular.The function d(x) is linear in δc.The terms that do not depend on δc are bounded.Let us fix k > m > 0, µ = −i.
Then the coefficient of δc is strictly positive and we can always fix δc such that the first term of d(x) is greater than the sum of the remaining terms.This way, we can keep d(x) > 0.
We are interested in the critical value of δc.We have We find that w(x) is an odd and strictly decreasing function, We define the critical value δc crit as follows Then both V 12 and V 13 are regular for |δc| ≥ |δc crit |.The potential term V 12 changes its sign asymptotically.In the limit of large |x|, it has the following behavior, The explicit form of V 12 and V 13 for |δc| > |δc crit | is in (75).
When δc = −δc crit , we have where The components V 12 and V 13 are depicted explicitly in Fig. 6.The behavior of both V 12 and V 13 is sensitive to the proximity of |δc| to |δc crit |.The width of the plateau in the two-step function and the width of the well increases as |δc| tends to |δc crit |.This time, the plateau corresponds to a non-vanishing energy.In Fig. 6, we present the plots of V 12 and V 13 for δc that is very close to δc crit but with varying m and M such that M 2 + m 2 is kept constant.For small values of m, V 12 forms a two-step function that resembles the potential from the previous model, see Fig. 6a), Fig. 6b).As m increases, there is formed a potential well in V 12 , see Fig. 6c) and 6d).In Fig. 7, there is the plot of the hopping amplitudes t 1 = t 2 + V 12 and t 3 = V 13 for δc = δc crit of the generalized stub lattice.In this figure, there is also the generalized stub lattice with hopping amplitudes t 1 , t 2 and t 3 in correspondence with (27).
Density of probability of the bound state Ψ−M is plotted in Fig. 6.The spectrum of this model has the form where |M | > |k| > |m|, which is further illustrated in Fig. 8.

Discussion
We presented the method for the spectral design of quasi-one-dimensional crystals with flatband that can be described effectively by the Dirac equation.Our approach is based on the susy transformation of the trivially extended pseudo-spin-one operator H 1 , see (13).The later operator is block diagonal with pseudo-spin-1/2 Hamiltonian H 1/2 and a constant on the diagonal.The operator H 1/2 governs the dynamics of a dimerized chain of atoms.The constant term represents an additional, parallel chain of atoms that are not interacting with their neighbors and host the flat band states.Susy transformation of H 1 provides us with the operator H 1 that already possesses non-trivial interaction between the two atomic chains, see (17) or (23).
The susy transformation of H 1 is partially defined in terms of flat-band solutions.The latter functions can be selected such that H 1 is hermitian, see (15).
The method was explained and explicitly illustrated in three explicit models where SSH chain interacts locally with a parallel chain of otherwise non-interacting atoms.In the models, the form of the interaction was tunable via free parameters, see Fig. 2, Fig. 4 and Fig. 6.We demonstrated that in the limit where the parameters approach the critical values, the interaction approximates piecewise constant potential, see Fig. 3, Fig. 5 and Fig. 7.
The trivial extension of H 1/2 to H 1 by diagonal constant λ fixes the energy of the flat-band, see (13).The remaining spectral characteristics of H 1 (and H 1 ) are inherited from those of H 1/2 .The spectrum of H 1/2 itself can also be adjusted by an appropriate susy transformation.The model discussed in Section 5.1 exemplifies such a situation.Indeed, the block-diagonal operator in (57) can be obtained from the free-particle Hamiltonian via susy transformation that generates two discrete energies E = ±λ in the spectrum of (49).It is worth noticing that Darboux transformation applied on pseudo-spin-1/2 system can generate up to two new bound states.When greater number of bound states is needed, it is possible to make a sequence of Darboux transformations to get the target Hamiltonian H 1/2 with the requested structure of discrete energies.Pseudo-spin-1/2 Dirac operators obtained via chains of Darboux transformations were discussed in the context of non-hermitian optics in [42].
The models described by H 1 also inherit scattering characteristics of H 1 .Both in Section 4 and in Section 5, H 1 was fixed as the reflectionless operator.The susy partners H 1 shared this property as they did not support any backscattering on the potential barriers.This behavior resembles Klein tunneling that occurs in electrostatic field, see the recent analysis in [53].
The work was inspired by [45] where spin-one free particle Dirac operator was transformed by supersymmetric transformation into the new ones with non-trivial potentials.In one specific case, the susy transformation resulted in decoupling of the Hamiltonian, i.e. it acquired blockdiagonal form with 2 × 2, spin-1/2 operator and a constant on the diagonal, see [45] for more details.The inverse approach was followed in the current article and exploited in the context of quasi-one-dimensional atomic chains.The presented construction based on susy transformation of the block-diagonal operators shares the philosophy of [57] where the interaction between uncoupled systems was induced via unitary transformations.
The presented approach to spectral design of quasi-one-dimensional systems with flat-bands is very flexible.It is straightforward to adjust it for the construction of quasi-one-dimensional systems with a higher number of flat-bands and/or higher number of atoms in the elementary cell of the dimerized chain.The number of the flat-bands as well as their energies can be easily controlled by addition of non-interacting parallel atomic chains to the initial Hamiltonian H 1/2 , i.e. via its trivial extension by corresponding number of diagonal constant terms.The susy transformation would generate the coupling between the original SSH-type chain(s) and the other, initially non-interacting, chains.It would be also possible to consider effects of magnetic field that would alter the phase of the hopping parameters.It is worth noticing in this context that generation of the synthetic magnetic field for optical lattices was proposed in [17], [23].Analysis of boundary effects on a finite lattice with the use of the supersymmetric transformation represents another interesting research direction as the boundary effects can play a major role e.g. in topological properties of the atomic chains [48].It would be also interesting to apply this approach in the analysis of the systems with quasi-bound states.Nevertheless, analysis of these topics goes beyond the scope of the present work and should be reported elsewhere.
A On the structure the matrix U Let us define the following set of 3 × 3 matrices The set M is closed with respect to matrix multiplication and inverse operation, i.e. it forms a group, which is manifestly non-hermitian.The vector b can be nullified by a specific choice of the seed solutions.Nevertheless, the new potential δV becomes again block-diagonal.
The latter reveals that to avoid non-hermitian potential terms, we must work with at least two flat-band states associated with the same flat-band level, rendering a transformation matrix with the structure This ensures that the new Hamiltonians constructed using the supersymmetric transformation showcase all the desired properties.

Figure 1 :
Figure 1: (a) Sketch of the periodic pseudo-one-dimensional stub lattice with three atoms A, B, and C per unitary cell (dashed rectangle).The hopping parameters are denoted by t 1 , t 2 , and t 3 .(b) Dispersion relations (5) in terms of the hopping parameters and the on-site interactions µ A and µ B = µ C , where ∆± = (t 1 ± t 2 ) 2 + t 2 3 + (µ A −µ C ) 2

Figure 8 :
Figure 8: Spectrum of H 1 .For each fixed value k, the two shaded bands represent the continuum.The thick black line corresponds to flat-band energy (73) and the thin black line is the energy of the bound state LΨ −M .We fixed m = 0.5 and M = √ 1 − m 2 .


belongs to M. If we select the matrix U such that U ∈ M, then U x ∈ M, with U x ≡ ∂ x U , leads to the relation δV = i[S 1 , U x U −1