Unconventional edge states in a two-leg ladder

Some popular mechanisms for restricting the diffusion of waves include introducing disorder (to provoke Anderson localization) and engineering topologically non-trivial phases (to allow for topological edge states to form). However, other methods for inducing somewhat localized states in elementary lattice models have been historically much less studied. Here we show how edge states can emerge within a simple two-leg ladder of coupled harmonic oscillators, where it is important to include interactions beyond those at the nearest neighbour range. Remarkably, depending upon the interplay between the coupling strength along the rungs of the ladder and the next-nearest neighbour coupling strength along one side of the ladder, edge states can indeed appear at particular energies. In a wonderful manifestation of a type of bulk-edge correspondence, these edge state energies correspond to the quantum number for which additional stationary points appear in the continuum bandstructure of the equivalent problem studied with periodic boundary conditions. Our theoretical results are relevant to a swathe of classical or quantum lattice model simulators, such that the proposed edge states may be useful for applications including waveguiding in metamaterials and quantum transport.


INTRODUCTION
The theory of the localization of electrons due to disorder, as created by randomness in the atomic lattice, is a towering achievement in condensed matter physics [1,2].This phenomena of an absence of diffusion can be observed throughout waves physics, including in acoustics, electromagnetics and photonics [3,4].Localization in perfectly ordered systems is also possible, most famously thanks to Wannier-Stark localization whereby an applied electric field breaks the periodicity of the lattice [5,6].Latterly, the rise of topological matter has celebrated the presence of highly localized states hosted by certain lattices, including topologically-protected edge states which exist in the gaps formed between bands filled with extended states [7][8][9].Furthermore and perhaps counterintuitively, the reality of so-called bound states in the continuum -in the broadest sense, localized states existing within a continuous spectrum of delocalized waves -have also been demonstrated in a series of photonic architectures [10][11][12][13] nearly a century after the original idea of von Neumann and Wigner [14].
Within condensed matter theory, the consideration of couplings going beyond the nearest-neighbour approximation is a perennially popular activity within standard treatments of lattice models [15][16][17][18][19]. Laterly, spiritually similar theoretical studies have been performed for acoustical, optical, and phononic platforms amongst others, with interesting implications for topological phases and for the creation of rather exotic dispersion relations [20][21][22][23][24][25][26].Experimentally, physical systems with non-negligible interactions beyond those of the nearest-neighbour type have been reliably produced, in some cases with tunability, using atomic spins inside a cavity [27]; acoustic metamaterials [28]; elastic metamaterials [29][30][31]; coupled optical cavities [32,33], and exciton-polariton condensates in microcavities [34,35].However, thus far the creation of edge states due to the inclusion of longer-range coupling has been mostly ignored, perhaps because of the seemingly contradictory concepts at play.
Here we propose exploiting longer-range couplings in order to generate edge states -without recourse to disorder, aperiodicity or topological non-trivialities.As a prototypical system, we employ a bosonic two-leg ladder model [36][37][38][39][40][41], which has already been explored experimentally to great effect using ultracold atoms for example [42][43][44][45][46][47].While previous studies of so-called 'biased' two-leg ladders have considered either interchain energy detunings across the two legs of the ladder or different twobody onsite interaction terms [48][49][50][51], here we consider a bias in the range of the couplings themselves.We suppose that the twoleg ladder is composed of an upper chain A [marked in pink in Fig. 1 (a)] which exhibits both nearest-nearest and next-nearestneighbour couplings, and a lower chain B [depicted in yellow in panel (a)] which only sustains nearest-neighbour interactions.Perhaps surprisingly, this coupling imbalance raises the possibility for edge states to emerge amongst the expected plethora of extended states which are spread throughout the ordered lattice, as implied in the probability density plots already previewed in Fig. 1 (b, c).Notably, the next-nearest-neighbour coupling introduces additional stationary points in the bandstructure (found in the continuum limit, after employing periodic boundary conditions) which occur for a certain quantum number q, which is associated with a definite eigenfrequency [52].This particular bandstructure eigenfrequency is in turn exactly the energy of the edge state (found after diagonalization of the analogous finite-sized lattice problem with the implied open boundary conditions), which suggests a sort of bulk-edge correspondence between the generation of auxiliary stationary points and the appearance of edge states.Furthermore, the discovered edge states appear inside the notional band of extended states, such they may be thought of as a kind of bound state in the continuum, a flavour of state which is of a high current interest within modern condensed matter physics [10][11][12][13].
The Hamiltonian operator Ĥ of the bosonic two-leg ladder sketched in Fig. 1 (a) is composed of just three parts where ĤA describes the N coupled harmonic oscillators comprising the chain A subsystem [depicted as pink spheres in Fig. 1 (a)] via the following Hamiltonian (we have set ℏ = 1 here and throughout) The nearest-neighbour coupling strength J 1 > 0 [red arrow in Fig. 1 (a)], and the next-nearest-neighbour coupling strength J 2 [cyan arrow] satisfies the inequality 0 ≤ J 2 ≤ J 1 .The N oscillators in the chain B subsystem [yellow spheres in Fig. 1 (a)], on the opposite side of the rung to their chain A counterparts, are governed by a standard nearest-neighbour tight-binding Hamiltonian ĤB , where The creation and annihilation operators c † n and c n create and destroy a bosonic excitation on the n-th rung of the ladder, where the two possible flavours of operator c n = {a n , b n } refer to the site at the top or bottom of the n-th rung, and where the bosonic commutation relation [c n , c † n ] = 1 is always observed.The common resonance frequency of all 2N oscillators throughout the two-leg ladder is ω 0 , although all of our results can easily be generalized to the case when the upper chain and lower chain have different resonance frequencies (this situation is briefly discussed later on in an appendix).Each oscillator is coupled to its adjacent partner along the rung (of the opposite flavour, A or B) via the coupling Hamiltonian The inter-chain coupling strength g > 0 [purple arrow in Fig. 1 (a)], such that the two dimensionless parameters of primary interest for the full model of Eq. ( 1) are the next-nearest neighbour coupling ratio J 2 /J 1 , with a size between 0 ≤ J 2 /J 1 ≤ 1, and the inter-chain coupling ratio g/J 1 ≥ 0. Together these coupling ratios control two important aspects of the overall model.Increasing J 2 /J 1 allows for more stationary points in the continuum Hamiltonian bandstructure to emerge [52], which raises the possibility for interference effects.For example, the group velocity and phase velocity may be in the same direction for one mode and in the opposite direction for another mode somewhere elsewhere in the Brillouin zone [22].Meanwhile, increasing g/J 1 allows for more mixing between the chain A and chain B subsystems which, due to their different boundary conditions coming from their different coupling ranges, raises the possibility for novel states to be forged in the combined ladder system.
In what follows, we discuss the band theory of the total Hamiltonian Ĥ of Eq. (1) in the infinitely long limit (and with periodic boundary conditions) in Sec. 2, which hints at novel states near to additional stationary points which arise due to interactions beyond nearest-neighbour when the coupling strength J 1 /4 ≤ J 2 ≤ J 1 .We consider the proper finite system defined by Eq. (1) in Sec. 3 (with its implied open boundary conditions), where we quantify the degree of localization of the eigenstates using the participation ratio.We also judge the character of the eigenvalue statistics using the mean value of the consecutive energy level spacings, which allows us to paint the extended state-to-edge state phase diagram of the model in Fig. 3. Some conclusions from our theoretical study are drawn in Sec. 4. We relegate to the appendices a complete treatment of the chain A subsystem model in the App.A, including a derivation of the analytic form of the density of states beyond nearest neighbour interactions, and an analysis of the chain B subsystem model in App.B, including a derivation of the analytic probability density function associated with the adjacent energy level spacings.Finally, some supplementary results for the full two-leg ladder are housed within App.C, while App.D briefly considers a small extension of that model to allow for interchain energy detunings.

BAND THEORY
Let us start by considering the full two-leg ladder model in the limit of large system size (N → ∞) and after employing periodic boundary conditions, which allows for an analytic treatment of the continuum theory.Our hope is for the continuum model to have some descriptive power over the finite-sized system, in a similar manner to how the so-called bulk-edge correspondence in the physics of topological matter suggests that a topological invariant can predict the presence or absence of topological edge states [8,9].We introduce the pair of exponential Fourier transforms for the bosonic operators a n and b n appearing in the ladder model of Eq. ( 1) as follows with the quantum number q = 2πm/N d, where the integer m ∈ [−N/2, +N/2] implies that the first Brillouin zone exists for −π/d < q ≤ π/d.The length scale d enters Eq. ( 5) due to the periodicity of the lattice, which repeats itself with identical rungs existing after every spacing d along the ladder [cf.Fig. 1 (a)].Substitution of Eq. ( 5) into the total Hamiltonian operator Ĥ of Eq. ( 1) leads to the q-space form of the Hamiltonian where the two eigenfrequencies ω A q and ω B q , which are associated with modes belonging to the individual chains A and B respectively, read which immediately highlights an important coupling asymmetry amongst the two subsystems.Rearranging the transformed Eq. ( 6) into the neater form where the 2 × 2 Bloch Hamiltonian H q and the column vector Ψ q (which contains the relevant operators) are together defined by reveals some important properties of the model.Namely, the model respects only one of three important Bloch Hamiltonian symmetries (time-reversal, inversion and chiral) as follows H q = H * −q , (time-reversal symmetry is observed) (11) (chiral symmetry is not observed) (13) where the the inversion operator I = σ x and the chiral operator C = σ z .Both operators are given in terms of one of the three Pauli spin matrices σ = (σ x , σ y , σ z ), while I is the identity matrix.As usual, time-reversal symmetry is respected since no magnetic (or even pseudo-magnetic) field is present [cf.Eq. ( 11)].Inversion symmetry is in general broken because as soon as the the next-nearest-neighbour coupling J 2 ̸ = 0, then the diagonal terms in the Bloch Hamiltonian H q are unequal (ω A q ̸ = ω B q ), such that the system is indeed altered upon interchanging type A oscillators with type B oscillators throughout the ladder [cf.Eq. ( 12)].Chiral (or sublattice) symmetry is not obeyed due to the nonzero couplings terms J 1 ̸ = 0 and J 2 ̸ = 0 in the on-diagonal terms of the Bloch Hamiltonian H q , which ensure that the spectrum is not symmetric about ω 0 (the trivial breaking of chiral symmetry by the constant energy term ω 0 , appearing in both ω A q and ω B q , is not consequential as it is just an arbitrary energy shift) [cf.Eq. ( 13)].This brief symmetry analysis suggests that the proposed model likely cannot be classified with a traditional topological index as per the Altland-Zirnbauer table [53,54], such that any edge states arising in the two-leg ladder should probably be deemed to be topologically trivial.Indeed, the winding number suggested by the curve formed in a parametric plot of the Bloch Hamiltonian H q = H I I + σ • H over the first Brillouin zone is zero.Here the threedimensional Hamiltonian vector H = (H x , H y , H z ), and the employed Pauli matrix decomposition implies that the on-diagonal term H I = ω 0 + 2J 1 cos (qd) + J 2 cos (2qd), the off-diagonal term H x = g is a constant, the leading diagonal imbalance term H z = J 2 cos (2qd), and finally H y = 0.
The eigenvalues of the Bloch Hamiltonian H q [cf.Eq. ( 10)] suggest the pair of eigenfrequencies ω q,τ of the two-leg ladder, which are grouped into two bands (as codified by the index τ = ±1) as follows The expression of Eq. ( 14) makes use of ωq , the average eigenfrequency of the uncoupled chains A and B encompassing the ladder, and the effective coupling constant Ω q , which are together defined by where ∆ q is the detuning across the two subsystem chains A and B. The eigendecomposition of Eq. (10), aided by the eigenfrequencies ω q,τ of Eq. ( 14), directly leads to the diagonal form of the original two-leg ladder Hamiltonian operator Ĥ of Eq. ( 1) like so where the twin Bogoliubov operators β q,τ introduced above are defined by the operator superpositions β q,+ = sin θ q a q − cos θ q b q , β q,− = cos θ q a q + sin θ q b q , with the two trigonometric Bogoliubov coefficients which have emerged are given by which are determinable from the two quantities Ω q and ∆ q , as defined in Eq. ( 16) and Eq. ( 17) respectively.The knowledge of the eigenfrequencies and eigenstates of the two-leg ladder allows for some general properties of the system to be inferred within standard band theory, and we restrict our analysis to the first Brillouin zone −π/d < q ≤ π/d only due to the periodicity of the considered lattice geometry.
The group velocity v q,τ = ∂ q ω q,τ for a mode in the band τ = ±1 and with the quantum number q, follows directly from Eq. ( 14) as Clearly there are always, that is for both bands τ = ±1 and for any values of J 2 /J 1 and g/J 1 , roots of the group velocity v q,τ for modes at the centre and edge of the Brillouin zone with the quantum numbers q = 0 and q = π/d respectively.For larger coupling ratios J 2 /J 1 , these two guaranteed roots are joined by additional stationary points.For example, in the extreme limit of g ≫ J 2 the group velocity v q,τ reduces to the band index independent quantity v q,τ = −2d [J 1 sin (qd) + J 2 sin (2qd)].Upon transforming the variable from qd to z = e iqd , the roots of the group velocity can be then be found by solving the quartic equation The already known roots of q = 0 and q = π/d correspond to the twin solutions of z = ±1, essentially reducing the problem to the effectively quadratic equation (z + 1)(z − 1)(J 2 z 2 + J 1 z + J 2 ) = 0.The final two roots are then found to be located at , or equivalently at qd = ±(π − arctan (2J 2 /J 1 ) 2 − 1).This brief analysis suggests additional stationary points only occur above a critical coupling ratio J 2 /J 1 > 1/2, at least in the considered regime g ≫ J 2 .A similar calculation may be performed in the opposing limit of g ≪ J 2 , when the two subsystems effectively decouple, which reveals a different threshold coupling ratio of J 2 /J 1 > 1/4 for additional stationary points to emerge (see App.A for the derivation).These additional stationary points in the bandstructure are found out to be highly consequential for the emergence of edge states, as is discussed later on.
The bandgap δω, defined as the energy difference between the minimum of the upper band ω q,+ and the maximum of the lower band ω q,− for any value of the quantum number q, is given by the intuitive formula which may be positive, negative or zero in the case of the considered ladder model.In particular, a negative bandgap δω suggests band overlap since the lower band will necessarily protrude into the territory of the nominally upper band.In the simplest case of very weak next-nearest neighbour coupling J 2 ≪ J 1 , a negative bandgap only exists when the inequality g < 2J 1 holds, since min{ω q,+ } = ω 0 − 2J 1 + g and max{ω q,− } = ω 0 + 2J 1 − g in this limit, such that the bandgap of Eq. ( 22) becomes δω = 2(g − 2J 1 ) within this simple regime.The influence of the bandgap δω on the appearance of edge states is discussed in detail in the next section.
Stationary points, bandstructure and degrees of localization in the two-leg ladder model.We consider five increasingly large values of the next-nearest neighbour coupling J2 upon descending the columns of the figure (as marked on the left-hand side), and we set the inter-chain coupling g = 3J1/2.First column: the group velocity vq,τ as a function of the quantum number q, in the first Brillouin zone −π/d < q ≤ π/d.The results for the τ = +1 (cyan lines) and τ = −1 (pink lines) bands are shown.The calculation is carried out with periodic boundary conditions for a two-leg ladder in the limit of 2N → ∞ oscillators [cf.Eq. ( 21)].In panels (g, j, m), the significant J2 couplings have given rise to extra roots in the group velocity (marked by pink and cyan circles), corresponding to additional stationary points in the bandstructure ωq,τ , which are similarly marked in panels (h, k, n).Second column: the eigenfrequencies ωq,τ forming the periodic bandstructure [cf.Eq. ( 14)] and corresponding to the results of the first column.Dashed pink and cyan horizontal lines: the extrema of each (τ = ±1) band, which are similarly carried through to the third column, and which imply a series of negative bandgaps [cf.Eq. ( 22)].Solid blue horizontal lines in panels (k) and (n): the modes associated with additional stationary points in the lower τ = −1 band, which are carried through to panels (l) and (o) in the third column.Third column: the scaled participation ratio PR (m) /2N of each state m [cf.Eq. ( 23)], which reside at the eigenfrequency ωm.Solid grey verticle lines: guide for the eye at 2/3.The numerical calculation in the third column is carried out with open boundary conditions for a two-leg ladder with 2N = 200 oscillators.
We present graphically some visualizations of a few key results coming from the bandstructure of the two-leg ladder in the first two columns of Fig. 2, for the example case with the inter-chain coupling g = 3J 1 /2.The first column of panels (a, d, g, j, m) in Fig. 2 display the group velocity v q,τ in the first Brillouin zone using Eq. ( 21), where descending the column increases the next-nearest neighbour coupling strength J 2 (as marked on the left-hand side of Fig. 2).The results for the upper band τ = +1 are denoted with magenta lines, while the lower τ = −1 band are associated with cyan lines.As discussed after Eq. ( 21), the guaranteed stationary points at q = 0 and q = π/d are joined by additional stationary points for the cases of larger next-nearest neighbour couplings J 2 , as marked by pink circles and cyan circles in the relevant Fig. 2 (g, j, m), where J 2 ≥ J 1 /2.The second column of Fig. 2 shows the eigenfrequencies ω q,τ corresponding to the cases of the first column and forming the exact bandstructure of the two-leg ladder via Eq.( 14).The dashed pink and cyan lines mark the extrema of each (τ = ±1) band, which imply a series of negative bandgaps [cf.Eq. ( 22)] within the chosen parameter regime.Importantly, the final three panels (h, k, n) in Fig. 2 contain pink circles and cyan circles in the same manner as the group velocity panels (g, j, m), which correspond to stationary points in the bandstructure.In particular, the solid blue lines in panels (k) and (n) mark the frequencies ω q,τ of the modes which are associated with additional stationary points in the lower τ = −1 band, which lies inside the energetic sector created by the upper and lower bounds of the lower τ = −1 band (as denoted by the two dashed cyan lines).As such, if these modes are edge states, and if they are surrounded by a collection of extended states comprising the rest of the effective band, the highlighted modes will in essence be a type of bound state in the continuum [10][11][12][13].In the next Sec.3, we finally consider the proper finite-sized two-leg ladder system as defined by Eq. ( 1), which allows us to probe the existence of edge states within the lattice model and their relation to certain bulk properties.

FINITE-SIZED SYSTEMS
Let us now consider the two-leg ladder Hamiltonian Ĥ as originally introduced in Eq. ( 1) for a ladder of finite length N , such that there are 2N coupled oscillators in the overall system.Upon staying in real space [instead of moving into q-space using the transformation of Eq. ( 5)] the Hamiltonian matrix H will be of size 2N × 2N [instead of the 2 × 2 Bloch Hamiltonian matrix H q of Eq. ( 10), as found within q-space].The diagonalization of the finite square matrix H (with implied open boundary conditions) is readily achieved, leading to 2N eigenfrequencies ω m , where the integer index m = {1, 2, ..., 2N }.For each state labelled by m, the n-th element of the associated eigenvector can be written as c n (m), and there will be 2N elements in total, corresponding to each site in the two-leg ladder.The knowledge of the eigenvectors allows for the spread of the state across the whole lattice to be quantified using a popular localization measure: the participation ratio PR (m) of a state m.This quantity may be defined via the formula [55,56] where the numerator is simply a normalization, while the denominator discriminates between states with different spatial extents.
In the case of the two-leg ladder considered here, the sums in Eq. ( 23) are taken over the whole array n = {1, 2, ..., 2N }.A fully localized state is trapped on just one site n, so that PR (m f. localized ) = 1.Conversely, a fully extended state will feel all 2N sites equally, such that its participation ratio will scale exactly with the system size as PR (m f. extended ) = 2N .Otherwise, typical extended states will observe a linear with system size relationship like PR (m extended ) ≃ β(2N ), where the prefactor β is some dimensionless number which is known exactly for certain simple cases (see Refs. [57,58] for example, while App.B contains a derivation of the common prefactor β = 2/3 for a single regular chain).Away from these common scenarios, that is those with participation ratios scaling like (2N ) 0 or (2N ) 1 , states exhibiting the more peculiar expression PR (m) ≃ β(2N ) α are of great interest.Here the exponent α satisfies the inequality 0 < α < 1, such that the state has a distinctly sublinear dependence on the overall system size 2N .These unusual states are then arguably localized to some degree, depending upon the exact size of the key exponent α.In some contexts (including quasiperiodic models), these states are sometimes known as multifractal states, where α is related to the fractal dimension of the state [59][60][61].
We explore the various degrees of localization of the states supported in the two-leg ladder in the third column of Fig. 2. As in the first and second columns of Fig. 2, the inter-chain coupling remains at g = 3J 1 /2 and going down the column of panels (c, f, i, l, o) increases the next-nearest neighbour coupling J 2 , but now we consider a finite-sized array of 2N = 200 oscillators as described by Eq. ( 1).We plot in Fig. 2 (c, f, i, l, o) the eigenfrequencies ω m , where each state is labelled by the index m, as a function of its participation ratio PR (m) [cf.Eq. ( 23)].The simplest case of zero next-nearest neighbour coupling (J 2 = 0) is displayed in Fig. 2 (c), where the common participation ratio of PR (m) ≃ (2/3)(2N ) for all states m is most apparent, signalling rather standard extended states (indeed, the prefactor of β ≃ 2/3 is obtained analytically in App.B).In Fig. 2 (f, i), where the next-nearest neighbour coupling J 2 is nonzero but still relatively small, there are fluctuations around the value of 2/3 (a number marked with a solid grey vertical line throughout the third column) suggesting typical extended states but without the complete uniformity of panel (c).Notably especially in panel (i), where J 2 = J 1 /2, two collections of states have started to accumulate around the band edges of the upper and lower band (as found from the continuum calculation, and as marked with the dashed magenta and cyan horizontal lines in each panel).The cases with the largest next-nearest neighbour coupling J 2 in the lowest two panels finally reveal the formation of significantly less extended states, with the smallest scaled participation ratio being min{PR (m)}/2N ≃ 0.171 in Fig. 2 (o), where J 2 = J 1 .This state of minimum participation ratio in panel (o) resides at a certain eigenfrequency ω m , as marked by a solid blue line, which corresponds to an eigenfrequency ω q,τ coming from the continuum model.This particular eigenfrequency ω q,τ is shown in the adjacent panel (n) by the same solid blue line.14), Eq. ( 21) and Eq.(22).In this panel, we consider a two-leg ladder with 2N = 200 oscillators.Panel (b): we take a horizontal cut through panel (a) at g = 3J1/2, and then we make vertical cuts at five different values of J2 (as noted in the plot legend), before extending the results as a function of the system size 2N in a log-log plot (up to 2N = 700).Panel (c): a map of the mean consecutive energy level spacings rm, as a function of both g and J2 [cf.Eq. ( 27)].In this panel, we consider a two-leg ladder with 2N = 2000 oscillators and the indicative black lines are the same as the lines provided in panel (a).
The specific quantum number q associated with this particular eigenfrequency ω q,τ is distinguished by a blue circle in panel (n), since it is a stationary point [all additional stationary points are marked with circles in the first two columns of Fig. 2, as discussed after Eq. ( 21) and as confirmed in the group velocity plots given in panels (g, j, m)].Notably, the most localized state in panel (o) is also seen to be housed inside the notional spectral band of extended states (limited by the dashed horizontal lines), a defining characteristic of a bound state in the continuum [10][11][12][13].
We visualize two characteristic states in the two-leg ladder in Fig. 1 (b, c), where the couplings g = 3J 1 /2 and J 2 = J 1 , and we consider a system with 2N = 200 oscillators [which matches the parameter situation in Fig. 2 (o)].We plot the probability density over the 2N sites of the array, where the results for the A and B sites n are shown with pink and yellow lines respectively.Notably, the example extended state shown in in Fig. 1 (b) feels around two-thirds of the array, its scaled participation ratio PR (m) /2N ≃ 0.675 and its oscillatory structure is noticeably balanced across both subsystems chain A and chain B. In Fig. 1 (c) we finally reveal the spatial structure of the more localized state discussed in Fig. 2 (o), of scaled participation ratio PR (m) /2N ≃ 0.171, which is seen to be a type of edge state due to its significant probability density at the ends of the ladder (across both chain A and chain B).The correspondence between the appearance of an edge state in the finite system and the emergence of additional stationary points in the continuum bandstructure is somewhat reminiscent of the celebrated bulk-edge correspondence in topological matter (where a calculation of a topological invariant from a continuum theory predicts the presence or absence of topological edge states in the finite version of the model [8,9]).
We chart the degrees of localization of states supported by the two-leg ladder system more generally in Fig. 3.As a function of both the next-nearest neighbour coupling ratio J 2 /J 1 and the inter-chain coupling ratio g/J 1 , we plot the minimum of the scaled participation ratio min{PR (m)}/2N in Fig. 3 (a), where we consider a finite system again composed of 2N = 200 oscillators.Blue and green regions in panel (a) represent cases where the ladder only supports highly extended states, which display the typical behaviour of PR (m) /2N ∼ 2/3, while more interestingly the orange-red regions in the map suggest that the system may host somewhat localized states, since the localization measure PR (m) /2N ∼ 1/10.The map of Fig. 3 (a) demonstrates the interplay between the key ratios J 2 /J 1 and g/J 1 : the next-nearest neighbour coupling ratio J 2 /J 1 needs to be sufficiently large (in this case, J 2 > J 1 /4) for a red phase to first appear, while if the inter-chain coupling ratio g/J 1 is too large (roughly, g > 2J 1 ) no red phase is supported.The general character of the diagram of Fig. 3 (a) can be understood using some key results from the band theory of the continuum version of the model.Explicitly, as was discussed after Eq. ( 21) and also after Eq. ( 22), some limiting cases of the properties of the two-leg ladder with periodic boundary conditions suggest that additional stationary points arise when additional stationary points arise when a negative bandgap arises when The inequality of Eq. ( 24) explains the critical behaviour along the bottom horizontal axis of Fig. 3 (a), as marked by the emergence of the red-orange colours, corresponding to the first appearance of stationary points due to next-nearest neighbour interactions.Meanwhile Eq. ( 25) is suggested at the top horizontal axis of Fig. 3 (a), due to the transition from blue to green regions.The band overlap arising from a negative bandgap δω arises when Eq. ( 26) is true, which actually approximately holds rather nicely for all of the range 0 ≤ J 2 ≤ J 1 .The solid black line in Fig. 3 (a) demarcates the exact line in parameter space below which the bandgap δω is negative [cf.Eq. ( 22)].Similarly, the dashed black line marks areas with different numbers of stationary points for the upper band (τ = +1), while the dotted black line does the same job for the lower band (τ = −1).These three guides for the eye, arising from purely bulk calculations, certainly seem to coincide with interesting behaviour of the corresponding finite system with edges.As an aside, when tending towards zero inter-chain coupling (g → 0) we find a singular limit: just above the horizontal axis of Fig. 3 (a) the red-orange colours when J 2 > J 1 /4 refer to significantly smaller participation ratios, but this behaviour is not found exactly at g = 0 (see for example the analyses of the individual chain A and chain B models in App.A and App.B respectively).The occurrence of this singular limit lends some support to the notion that it is important for greater degrees of localization that the system contains subsystems with different boundary conditions, as here the upper and lower chain subsystems have contrasting boundary conditions due to their different coupling ranges (which is only not relevant exactly when g = 0, as then the subsystems are completely decoupled).
A more accurate categorization of the nature of the extended states and edge states supported by the two-leg ladder requires an analysis which tracks changes due to size-dependencies.In Fig. 3 (b) we effectively take a horizontal cut through panel (a) at the specific inter-chain coupling g = 3J 1 /2, and then we make vertical cuts at five different values of J 2 (as noted in the plot legend).We then extending the result, the minimum of the scaled participation ratio min{PR (m)}/2N , as a function of the total number of oscillators in the ladder 2N .This process results in the log-log plot presented as Fig. 3 (b).Four of the lines in Fig. 3 (b) effectively tend towards a constant (or oscillate between two constants in the case of the dark green line, where J 2 = 0) which suggests that the participation ratio scales linearly with 2N : this data truly describes an extended phase of the system.However, the red line in Fig. 3 (b), representing the case of J 2 = J 1 , does not seem to tend towards a constant for the considered system sizes such that the participation ratio scales sublinearly with 2N (at least within this regime up to 2N = 700).This suggests an unconventional phase containing edge states [cf. the edge state depicted in Fig. 1 (c) for the case of 2N = 200 oscillators], albeit not fully localized states since the relationship is not completely independent of 2N .Similar plots to Fig. 3 (b) may be found in App.C for other values of the specific inter-chain coupling g, providing further information about the overall phase diagram of Fig. 3 (a).Notably, we draw our conclusions from these aforementioned finite size effect plots based upon system sizes from ∼ 10 1 to ∼ 10 3 , since this regime can plausibly be explored experimentally [27][28][29].We leave to a later study the numerical crunching of two-leg ladders containing larger numbers of oscillators, where the already-discovered non-monotonic nature of the data may revise some of the scalings between some ranges of exceedingly large 2N (in which case, the potential consequences will only effect the experimental study of extremely large systems).We also briefly discuss in App.D the impact of a nonzero energy detuning between the resonance frequencies in the upper chain and lower chain as a complement to the data presented in Fig. 3, an analysis which also may be useful to guide certain future experiments and which hints at the possibility for inducing stronger degrees of localization thanks to energy detuning.
We further characterize the behaviour of the two-leg ladder by considering the energy level statistics of the model [62,63].For a generic 2N × 2N matrix Hamiltonian, we may list the resulting eigenvalues E m in ascending order, where the index m = {1, 2, ..., 2N }.Therefore, there are 2N − 1 consecutive energy level spacings S m , where S m = E m+1 − E m .The ratio of two consecutive gaps r m may then be defined as [64,65] which observes the bounds 0 ≤ r m ≤ 1, where the index runs for m = {2, 3, ..., 2N − 1}.The dimensionless quantity r m measures the correlations between the adjacent energy level gaps in a given spectrum.The mean value r m , found after averaging over all 2N − 1 ratios r m , acts as a kind of judge of the chaoticity appearing in the system.Some typical values of r m for different categories of energy-level statistics are provided in Refs.[57,65]: Poissonian level statistics are associated with r m = 2 ln 2 − 1 ≃ 0.39, while Gaussian orthogonal ensembles are characterized with r m = 4 − 2 √ 3 ≃ 0.54.In ordered (picket fence-like) models r m should essentially be unity.Localization is implied with Poissonian-like values of r m , since localized states can reside in different single-particle basis states of nevertheless similar energies, such that they will not interact or display any significant level repulsion [64].In the case of the considered two-leg ladder model, we associate the eigenfrequencies ω m with the energy levels E m alluded to above, and we apply averaging over Eq. ( 27) in order to generate the results presented in Fig. 3 (c).In the map of panel (c), we plot r m as function of both the next-nearest neighbour coupling ratio J 2 /J 1 and the inter-chain coupling ratio g/J 1 , while the black lines acting as guides for the eye are exactly those of Fig. 3 (a).Notably, the yellow region of Fig. 3 (c) corresponds to an ordered system with r m ≃ 1, implying no collections of localized states.Meanwhile, different flavours of extended phases are suggested by the orange and red regions at the top of panel (c), which are approximately encased by the indicative black lines.The lower half of Fig. 3 (c) sees red areas gradually merge into purple islands, in particular the enclosed region ending at J 2 = J 1 , which are associated with r m ≃ 0.4.These kinds of averaged energy level statistics instead suggest a (relatively) more localized phase, as is consistent with the participation ratio map of Fig. 3 (a).Similar conclusions based upon the averaging of the quantity defined in Eq. ( 27) have recently been drawn in both theoretical and experimental studies of a variety of other physical systems hosting collections of relatively localized states [66][67][68].

DISCUSSION
In conclusion, we have presented a basic theory for the generation of unconventional edge states in bipartite coupled mode systems due to the presence of asymmetric interactions.We have concentrated on the simplest case of a two-leg ladder model, where we allowed for different coupling ranges in the upper chain as compared to the lower chain.The long-range coupling leads to additional stationary points developing in the bandstructure which, in a remarkable display of a seemingly non-topological bulk-edge correspondence, predicts the existence of states residing at the edges of the analogous finite-sized lattice.Since these edge states are essentially surrounded by a band of extended states they are in some sense a kind of bound state in the continuum [10][11][12][13], albeit generated by the unusual mechanism of asymmetric coupling ranges.Our theoretical predictions could be realized in various physical platforms modelled by lattices with significant next-nearest neighbour interactions, for example with atomic, photonic, phononic and polaritonic setups [27][28][29][30][31][32][33][34][35].Perhaps the acoustic metamaterial platform of Ref. [28] is the most promising avenue for immediate experimental realization due to its tunability.Finally, although we have concentrated on a toy model of a two-leg ladder for simplicity, our results should generalize to any bipartite systems with asymmetries in their respective couplings, for example matter modes interacting with radiative modes [69][70][71][72].

Appendix A: The chain A model
In this appendix we consider the chain A model only, as defined by the Hamiltonian operator ĤA of Eq. ( 2).It describes a linear chain of quantum harmonic oscillators (each of resonance frequency ω 0 ) interacting with both nearest-neighbour couplings J 1 and -most importantly -also next-nearest-neighbour couplings J 2 [see the sketch of chain A in Fig. 1 (a)].We are especially interested in the continuum limit of a very long chain (N → ∞), and we commonly employ periodic boundary conditions in order to obtain an analytic understanding of both the band theory and the energy level statistics of this subsystem model.In particular, we analyze the additional stationary points in the bandstructure which appear for stronger next-nearest neighbour couplings, we present a detailed derivation of the analytic form of the density of states beyond nearest-neighbour interactions, and we explore in detail the probability density function of adjacent energy levels.
The exponential Fourier transform of Eq. ( 5) immediately diagonalizes the chain A Hamiltonian of Eq. ( 2), which reveals the eigenfrequencies ω A q associated with the quantum number q as follows ĤA = q ω A q a † q a q , (A 1) The assumed periodic boundary conditions ensure the 2π-periodicity of the eigensystem, allowing us to focus on the first Brillouin zone −π/d < q ≤ π/d only.Essentially, the nonzero inter-site couplings J 1 and J 2 in the chain A model allow for excitations to move throughout the one-dimensional lattice, which gives rise to kinetic energy terms [cf. the second and third terms in Eq (A 2)] in addition to the onsite energy term of each individual oscillator [cf. the first term in Eq (A 2)].The energy bands formed from ω A q are plotted in the first column of Fig. 4 using Eq (A 2), for four increasingly large values of the next-nearest neighbour coupling J 2 upon descending the panels (a, e, i, m).Most notably, with larger next-nearest couplings J 2 in Fig. 4 (i, m) the bandstructure exhibits additional stationary points, as marked by coloured circles in the figure .The extrema ω A ± of the bandstructure defined by Eq (A 2) are also noticeably affected by sufficiently strong J 2 , since the maximum and minimum eigenfrequencies are given by ) The maximum eigenfrequency ω A + always occurs at the quantum number q = 0, irrespective of the value of J 2 .However, the minimum eigenfrequency ω A − can take on two distinct values, either ω A I or ω A II , depending upon the crucial coupling ratio J 2 /J 1 .

II
energy ω A q , follow directly from taking derivatives of the dispersion of Eq. (A 2) with respect to the quantum number q, like so These two fundamental quantities are plotted in the second and third columns of Fig. 4 respectively, and they correspond to the bandstructures ω A q plotted in the adjacent panels of the first column of Fig. 4. In particular, the critical values of q corresponding to additional stationary points in ω A q are similarly marked by coloured circles in the panels (j) and (n), which are exactly the additional zeroes in the group velocity v A q due to the non-negligible next-nearest coupling J 2 .Notably, these additional stationary points are also marked by coloured circles in the inverse effective mass 1/m A q plots of panels (k) and (o), which implies that these states have a rather small effective mass m A q and hence that they are more prone to moving around the lattice.Importantly, while the inverse effective mass 1/m A q is a sign-changing quantity as a function of the quantum number q throughout the third column in Fig. 4, the lowest two panels (k) and (o) have additional roots in the inverse effective mass at certain quantum numbers q = q * , which are relatively near to the Brillouin zone edges (that is, they are not the roots in the vicinity of qd = ±π/2).These q * modes indeed satisfy 1/m A q * = 0, where These special states with the quantum number q * correspond to excitations with infinite effective masses, which perhaps implies a greater tendency for them to localize.This is discussed quantitatively later on around Eq. (A 34) and in Fig. 6 (b) when a finite array is considered, and the underpinning mechanism is linked to the ideas of Sajeev John and co-workers [73].
The density of states g (ω) counts the number of states at a given frequency ω.It is defined for a system of N energy levels, each residing at a certain frequency ω n , via the informal summation formula where δ(x) is Dirac's delta function, and where the integral over the density of states ∞ −∞ g (ω) dω = N indeed recovers the total number of energy levels.With regard to the considered Hamiltonian operator ĤA of Eq. (A 1), where the quantum number q = 2πm/N d due to the employed periodic boundary conditions, one may rewrite Eq. (A 11) in the continuum limit of N → ∞ oscillators as the density of states integral where the integral is over the whole of the first Brillouin zone, and where the eigenfrequencies ω A q are provided by the dispersion relation of Eq. (A 2).The Dirac delta function of some function f (x) is conveniently expressable using the identity where x n are the roots of the function f (x) so that f (x n ) = 0, and where the first derivative f ′ (x) = d dx f (x).Working in the dimensionless variable x, and introducing the dimensionless quantity a and the dimensionless coupling strength parameter b, as we may rewrite the argument of the Dirac delta function appearing in the density of states of Eq. (A 12) as follows of states are described by the pleasingly short expressions which all appear graphically in the top row of Fig. 4 [cyan lines in panels (a, b, c, d)], where we consider J 2 = 0 as the baseline case from which the more complicated behaviours in the lower rows eventually emerge.In particular, the cosine band of Eq. (A 28) leads to standard Van Hove singularities in the density of states of Eq. (A 33) and Fig. 4 (d), which are typical in one-dimensional systems.Clearly from Eq. (A 31), there are only stationary points in the bandstructure at qd = 0 and qd = π, while zeroes of the inverse effective mass only occur at qd = ±π/2, as shown graphically in Fig. 4 (b, c).In a future study, it may be interesting to consider the influence of the counter-rotating terms [74,75] discarded from the coupling part of the Hamiltonian, as well as to examine the impact of including an interaction term [76,77] in the Hamiltonian in order to probe a nonlinear array of anharmonic oscillators.
Thus far, we have exploited periodic boundary conditions and considered the limit of an infinitely large system (N → ∞) in order to study some fundamental band properties of chain A analytically.We shall now consider open boundary conditions with the Hamiltonian operator ĤA for a finite system of N oscillators, which allows us to study the degrees of localization and energy level statistics of the model.The theoretical treatment requires the diagonalization of the N × N matrix form of Eq. ( 2), which results in a finite number of eigenfrequencies ω A m for the countable number of states, which are labelled from m = {1, 2, ..., N }.The values of the associated eigenvectors, on oscillator n and for a state m, are given by the amplitude c A n (m).We show some typical eigenfrequency results, as a function of the coupling ratio J 2 /J 1 , for a chain of N = 10 oscillators with the yellow lines in Fig. 5 (a).As a guide for the eye, we include the continuum band edge results of ω A ± as the thick purple lines, using Eq.(A 3) and Eq.(A 4), which act as the energetic bounds.We see from Fig. 5 (a) how the collection of ten eigenfrequencies ω A m evolve with increasing J 2 /J 1 , and in particular we notice the crucial point around J 2 = J 1 /4, where the minimum eigenfrequency finally starts to lower due to the appearance of a new stationary point in the analogous infinite system [cf.Eq. (A 4)], and where several energy levels start to cluster together, which is important for what follows.
We quantify the degrees of localization of the states m sustained by chain A using the participation ratio PR (m), as defined in Eq. (23).Typically, this measure takes the form of PR (m) = (2/3)(N + 1) or similar in conventional, one-dimensional tightbinding systems [as is shown later on in Eq. (B 6)].The linear dependence on N is characteristic of all extended states, while the prevalent prefactor of 2/3 arises from the geometry of the linear chain with open boundary conditions.In Fig. 6 (a) we are interested in the three measures, all of which are linked to the participation ratio and are scaled by the convenient factor N + 1, as a function of the coupling ratio J 2 /J 1 .The red line in Fig. 6 (a) denotes max{PR (m)}, the maximum of the participation ratios (found after consideration of all states m).With vanishing J 2 and for J 2 > J 1 /2, this quantity is essentially the constant We plot as a function of the next-nearest neighbour coupling J2 (in units of the nearest neighbour coupling J1), and here we consider a finite system of N = 300 oscillators.Red line: max{PR (m)}, the maximum of the participation ratios of all states m.Orange line: PR (m), the mean value of the participation ratios.Yellow line: min{PR (m)}, the minimum of the participation ratios.Panel (b): some chosen eigenfrequencies as a function of J2.Thin yellow line: the specific eigenfrequency ω A M corresponding to the state with the index M , that with the smallest participation ratio.Pink line: the eigenfrequency ω A q * , corresponding to the state with the quantum number q * , which is associated with the appearance of an additional root in the inverse effective mass in the continuum limit [cf.Eq. (A 34)].Purple lines: ω A ± , the extrema of the continuum model band [cf.Eq. (A 3) and Eq.(A 4)].Panel (c): a log-log plot of the minimum of the scaled participation ratio min{PR (m)}/(N + 1) as a function of the system size N .We consider the coupling ratios J2/J1 = {0, 1/4, 1/2, 1} with the cyan, pink, green and orange lines respectively.The result for J2 = 0 (cyan line) is analytic [cf.Eq. (B 7)].max{PR (m)}/(N + 1) ≃ 2/3 [this relation is exact for J 2 = 0].This implies a linear relationship between the participation ratio and the system size, while the 2/3 prefactor is typical of extended states.Interestingly, for 0 < J 2 < J 1 /2 a larger prefactor associated with max{PR (m)} is demonstrated, which perhaps indicates the presence of some unusually super-extended states due to the non-negligible next-nearest neighbour hoppings.Similarly, the yellow line in Fig. 6 (a) denotes min{PR (m)}, graphing the minimum of the participation ratio (over all states m) scaled by N + 1.The additional stationary points suddenly occurring around J 2 = J 1 /4 are seemingly correlated with an unusually low prefactor of ∼ 0.4 in this case, which may perhaps be termed sub-extended states since the linear scaling with N is maintained, but less than one half of the total number of sites are felt by the state.Finally, the mean of the participation ratio PR (m), calculated as the average over all N states, is represented by the orange line in Fig. 6 (a).The orange curve suggests that the average state does not behave all that differently below and above the critical coupling of J 2 /J 1 = 1/4.
The calculations leading to the results of Fig. 6 (a) were carried out for a finite chain of N = 300 oscillators, which also led to the results of panel Fig. 6 (b).There we plot some key eigenfrequencies as a function of the coupling ratio J 2 /J 1 .Similarly to Fig. 5 (a), the purple lines mark the continuum band edges ω A ± as a guide for the eye at the energetic bounds.More importantly, the yellow line in Fig. 6 (b) tracks the special eigenfrequency ω A M , coming from the particular state m = M with the smallest participation ratio min{PR (m)}.For the important coupling regime of J 1 /4 < J 2 < J 1 , where the participation ratio of the state M is appreciably smaller than for all other states [cf.Fig. 6 (a)], the key eigenfrequency ω A M is remarkably seen to be not too dissimilar to the analytic expression as plotted with the pink line in Fig. 6 (b).We arrived at Eq. (A 34) by substituting the specific the quantum number q = q * [cf.Eq. (A 10)], which is associated with a zero in the inverse effective mass 1/m A q from Eq. (A 9), into the bandstructure of Eq. (A 2).The approximate likeness of ω A M and Eq.(A 34), as suggested by Fig. 6 (b), can be interpreted as this particular excitation with q * having a large effective inertia, such that the state is much more resistant to extending over the entire onedimensional array of coupled oscillators.
We confirm the extended nature of all of the states sustained in the chain A model in Fig. 6 (c), which plots the minimum of the scaled participation ratio min{PR (m)}/(N + 1) as a function of the overall system size N .We consider four values of J 2 , corresponding to the four cases considered previously in Fig. 4. Clearly, all of the cases quickly tend towards constant values, revealing the linear dependence of this localization measure with the total number of oscillators N in the chain.In particular, the case of J 2 = 0 (cyan line) is describable analytically [as shown later on in Eq. (B 7), prompting the choice of scaling with N + 1 rather than simply N ] and it exhibits a triangular waveform due to an even-odd relationship with the system size N .
The energy level statistics of chain A are investigated in Fig. 7.In particular, we consider two closely related quantities describing energy level correlations in Fig. 7 (a) and (b) [62,63].Meanwhile Fig. 7 (c) is focussed on the statistics of consecutive energy level spacings [64,65], as was also investigated in Fig. 3 (c) for the two-leg ladder model of the main text.We previously described with yellow lines in Fig. 5 (a) the eigenfrequencies ω A m , as a function of the coupling ratio J 2 /J 1 , for a relatively small chain of N = 10 oscillators.Taking four cuts through Fig. 5 (a) at the coupling ratios J 2 /J 1 = {0, 1/4, 1/2, 1} results in the energy level diagram sketched in Fig. 5 (b), where some level clusterings become more obviously apparent than in Fig. 5 (a).With these energy levels E m = ω A m , ordered such that m = 1 corresponds to the smallest eigenfrequency and ascending to the largest eigenfrequencies when m = N , we may define the N − 1 energy level separations S m = E m+1 − E m , where the index m = {1, 2, ..., N − 1}.The mean energy level spacing is the average of these energetic distances D = S m , such that one may define the N − 1 normalized energy level spacings t m = S m /D.We plot these normalized and dimensionless energy level spacings t m in Fig. 5 (c), corresponding to the eigenfrequencies of Fig. 5 (b).In the thermodynamic limit of N → ∞, the distribution of t m can then be described within the statistical theory of energy levels [62,63] as we now discuss.
The cumulative distribution function C(t) of the adjacent energy levels, as a function of the continuum energy level spacing t, is plotted in Fig. 7 (a) for four representative values of J 2 within the chain A model.There are no energy level spacings with a spacing below a vanishing amount such that C(t → 0) = 0, while for some sufficiently large value of the dimensionless energy level spacing t = T all energy level spacings in the chain A model should be smaller than it (due to the collection of energy levels packing themselves into an effective band), such that C(t ≥ T ) = 1.The simplest case of J 2 = 0 is denoted by the cyan line in Fig. 7 (a), whose cumulative distribution curve is obtained analytically with the exact inverse trigonometric expression where the threshold dimensionless energy level spacing T = π/2, as is marked at the top of Fig. 7 (a).The derivation of Eq. (A 35) is provided later on in App.B, where the chain B model is analyzed.The results for three other cases with J 2 ̸ = 0, corresponding to the values considered in Fig. 5 (b, c), are shown in Fig. 7 (a) after a numerical calculation for a finite system composed of N = 15000 oscillators.This leads to successively larger values energy level spacings T for which C(T ) = 1 first occurs, due to the increasingly large bandwidths of the effective bands formed by the collections of energy levels.These approximate threshold T values are noted at the top of Fig. 7 (a).The probability density functions p(s) associated with the four curves of Fig. 7 (a), such that the integral t 0 p (s) ds = C (t), are plotted in Fig. 7 (b).The probability density function p(s) in the simplest case of J 2 = 0 follows from taking the derivative of Eq. (A 35) with respect to t and evaluating the result at the dimensionless value s, leading to the neat formula which indeed satisfies the normalization ∞ 0 p (s) ds = 1 in order to conserve the total probability.Furthermore, the first moment obeys ∞ 0 sp (s) ds = 1, since the mean level spacing was normalized to unity.The inverse-square root distribution of Eq. (A 36) is shown by the cyan line Fig. 7 (b), demonstrating a lack of energy level repulsion since p(s → 0) = (2/π) 2 ≃ 0.405.This is followed by an increasing likelihood of larger level separations until the critical value of s = π/2 is reached, after which there is a step function drop-off to zero due to the finite bandwidth of the chain A model.The rough numerical results for the three nonzero values of J 2 considered in Fig. 7 (a) are also shown in Fig. 7 (b), which display significantly enhanced clusterings of the energy levels due to the s → 0 behaviour, as was already hinted at in the level spacing diagram of Fig. 5 (c).The critical behaviour (abrupt dropping to zero) of the J 2 ̸ = 0 cases in Fig. 7 (b) appears at increasingly large values of s, as is expected from the analytic case of J 2 = 0 and Eq.(A 36), where the criticality lies at s = π/2.
For a given spectrum, the ratio of consecutive energy level spacings r m , as defined by Eq. ( 27), provides a measure of the correlations between adjacent energy level gaps [64].The mean value r m provides a certain sense of the degree of chaoticity, being essentially unity in ordered models and dropping to r m = 2 ln 2 − 1 ≃ 0.39 for models with Poissonian level statistics, which sometimes implies a more localized phase [57,65].The evolution of mean ratio r m for the chain A model with the coupling ratio J 2 /J 1 is depicted in Fig. 7 (c).We consider the results for three increasingly large magnitudes of the finite system size N with increasingly dark green lines.Strikingly, while for weak next-nearest couplings J 2 < J 1 /4 the thermodynamic limit of N → ∞ sees the mean ratio approach r m = 1, corresponding to an ordered system, the regime of strong next-nearest couplings J 1 /4 < J 2 ≤ J 1 is characterized by a decay from unity.This drop-off starts abruptly at the key coupling J 2 = J 1 /4, falling towards a value of r m ≃ 0.55 when J 2 = J 1 (in the thermodynamic N → ∞ limit).This mean value r m is somewhat similar to the case of a Gaussian orthogonal ensemble, where r m = 4 − 2 √ 3 ≃ 0.54 [81], suggesting a kind of extended phase in the chain A model which is consistent with the participation ratio results of Fig.   alongside the cumulative distribution functions of some other celebrated models.The picket fence model with equal energy level spacings necessitates a step function behaviour (red line), the uncorrelated energy levels model with random energy level spacings yields a smooth step function (green line), and a chaotic model described by a Gaussian orthogonal ensemble essentially exhibits a Gaussian function (yellow line) [62,63].The probability density function p(s) associated with Eq. (B 13) was given in Eq. (A 36), and it is represented by the cyan line in Fig. 8 (b) along with the other distributions corresponding to the cumulative curves of Fig. 8 (a).Namely, the picket fence model has a Dirac delta function distribution (red line) due to the completely equal level spacings, while the uncorrelated model has a Poisson distribution (green line) implying that small level spacings are highly likely since p(s → 0) = 1, while large level spacings are vanishingly unlikely with an exponential fall-off p(s) = e −s .The chaotic model has a Wigner distribution (yellow line) which suggests significant energy level repulsion due to the behaviour p(s → 0) = 0, along with vanishing chances for large level spacings due to the very fast Gaussian decay p(s → ∞) ≃ e −πs 2 /4 .Interestingly, the inverse-square root distribution of the chain B model as described by Eq. (A 36) (cyan line) demonstrates the absence of energy level repulsion since p(s → 0) = (2/π) 2 ≃ 0.405, after which increasingly large level separations are increasingly likely up to the critical value of s = π/2.Above this value there is a brutal step function drop-off to zero due to the effective formation of a finite cosine band of energy levels [cf.Eq. (B 4)].

FIG. 1 .
FIG. 1. Unconventional edge states in the bosonic two-leg ladder model.Panel (a): a sketch of the two-leg ladder system under consideration, composed of two linear chains (denoted A and B) of quantum harmonic oscillators.Each oscillator is coupled (with the strength g, purple arrows) to its partner across the rung.The upper chain (A, pink balls) exhibits both nearest-neighbour coupling (with the strength J1, red arrows) and next-nearest-neighbour coupling (with the strength J2, cyan arrows).The lower chain (B, yellow balls) displays nearest-neighbour coupling (with J1, red arrows) only.The rungs are equally spaced at intervals of the distance d along the ladder.Panel (b): the probability density of a typical extended state over the 2N sites of the two-leg ladder.The results for the A and B sites n are shown with pink and yellow lines respectively.Panel (c): the equivalent profile of a typical edge state.In the figure, we consider a two-leg ladder with a total of 2N = 200 oscillators, and the couplings g = 3J1/2 and J2 = J1.The scaled participation ratio PR (m) /2N of the states considered is the rather typical ≃ 0.675 in panel (b) and perhaps surprisingly just ≃ 0.171 in panel (c).

FIG. 3 .
FIG.3.Degrees of localization and energy level statistics in the two-leg ladder model.Panel (a): a map of the minimum of the scaled participation ratio min{PR (m)}/2N , as a function of both the inter-chain coupling strength g and the next-nearest neighbour coupling strength J2, both in units of the nearest-neighbour coupling J1 [cf.Eq. (23)].Dashed black line: boundary marking areas with different numbers of stationary points for the upper τ = +1 band.Dotted black line: boundary marking areas with different numbers of stationary points for the lower τ = −1 band.Solid black line: the borderline below which there is a negative bandgap.These indicative black lines are calculated from the continuum results of Eq. (14), Eq. (21) and Eq.(22).In this panel, we consider a two-leg ladder with 2N = 200 oscillators.Panel (b): we take a horizontal cut through panel (a) at g = 3J1/2, and then we make vertical cuts at five different values of J2 (as noted in the plot legend), before extending the results as a function of the system size 2N in a log-log plot (up to 2N = 700).Panel (c): a map of the mean consecutive energy level spacings rm, as a function of both g and J2 [cf.Eq. (27)].In this panel, we consider a two-leg ladder with 2N = 2000 oscillators and the indicative black lines are the same as the lines provided in panel (a).

FIG. 5 .
FIG. 5. Eigenfrequencies of the chain A for a finite system with open boundary conditions.Panel (a): the eigenfrequencies ω A m for a finite system of N = 10 oscillators, as a function of the next-nearest neighbour coupling J2 (in units of the nearest neighbour coupling J1), as marked by the yellow lines.Thick purples lines: ω A ± , the extrema of the equivalent band formed in the continuum limit N → ∞ [cf.Eq. (A 3) and Eq.(A 4)].Panel (b): the specific eigenfrequencies ω A m , taken from panel (a) at the chosen values of J2/J1 = {0, 1/4, 1/2, 1}, as associated with the cyan, pink, green and orange lines respectively.Purple lines: the extrema ω A ± in the continuum limit.Panel (c): the relative energy level spacings tm between the eigenfrequencies ω A m taken from panel (b).

FIG. 6 .
FIG.6.Degrees of localization in the chain A model.Panel (a): various localization measures, scaled by N + 1, related to the participation ratio PR (m) of a state m.We plot as a function of the next-nearest neighbour coupling J2 (in units of the nearest neighbour coupling J1), and here we consider a finite system of N = 300 oscillators.Red line: max{PR (m)}, the maximum of the participation ratios of all states m.Orange line: PR (m), the mean value of the participation ratios.Yellow line: min{PR (m)}, the minimum of the participation ratios.Panel (b): some chosen eigenfrequencies as a function of J2.Thin yellow line: the specific eigenfrequency ω A M corresponding to the state with the index M , that with the smallest participation ratio.Pink line: the eigenfrequency ω A q * , corresponding to the state with the quantum number q * , which is associated with the appearance of an additional root in the inverse effective mass in the continuum limit [cf.Eq. (A 34)].Purple lines: ω A ± , the extrema of the continuum model band [cf.Eq. (A 3) and Eq.(A 4)].Panel (c): a log-log plot of the minimum of the scaled participation ratio min{PR (m)}/(N + 1) as a function of the system size N .We consider the coupling ratios J2/J1 = {0, 1/4, 1/2, 1} with the cyan, pink, green and orange lines respectively.The result for J2 = 0 (cyan line) is analytic [cf.Eq. (B 7)].

FIG. 8 .
FIG. 8.A comparison of the energy level statistics of the chain B model with some celebrated models.Panel (a): the cumulative distribution function C(t) as a function of the dimensionless energy level spacing t for four different models.Panel (b): the probability density functions p(s) corresponding to the curves of panel (a).Red lines: the picket fence model, where the equal energy level spacing leads to a Dirac delta function distribution.Green lines: the uncorrelated energy levels model, where the random energy level spacing leads to a Poisson distribution.Yellow lines: a chaotic model, where the Gaussian orthogonal ensemble leads to the simplest type of Wigner distribution.Cyan lines: the chain B model, as described analytically by the cumulative and density functions of Eq. (A 35) and Eq.(A 36) respectively.

FIG. 10 .
FIG. 10.Degrees of localization in the two-leg ladder model with detuning.Maps of the minimum of the scaled participation ratio min{PR (m)}/2N , as a function of both the inter-chain coupling strength g and the next-nearest neighbour coupling strength J2, both in units of the nearest-neighbour coupling J1 [cf.Eq. (23)].We consider five different detunings D = ωA − ωB, defined via the resonance frequencies ωA and ωB of the oscillators in the upper chain and lower chain respectively, across the row of panels as marked as the top of the figure.The case of the central panel (c) recovers the model of the main text, where D = 0 (since ωA = ωB = ω0).In this figure, we consider a two-leg ladder with 2N = 200 oscillators.

2 FIG. 11 .
FIG.11.Finite size effects in the two-leg ladder model with detuning.Log-log plots of the minimum (taken over all eigenstates m) of the scaled participation ratio PR (m) /2N as a function of 2N , the total number of oscillators in the chain A-B system with detuning [cf.Eq. (23)].We consider five values of the inter-chain detuning D [as marked by the plot legend in panel (a)], corresponding to the values in Fig.10 (a-e).This figure supplements the results of Fig.10by effectively taking horizontal cuts through Fig.10 (a-e) at three values of g [as marked at the top of this figure] and a single vertical cut at J2 = J1, and then extending the outcome as a function of the system size 2N (up to 2N = 700).