Robust edge states induced by electron-phonon interaction in graphene nanoribbons

The search of new means of generating and controlling topological states of matter is at the front of many joint efforts, including bandgap engineering by doping and light-induced topological states. Most of our understading, however, is based on a single particle picture. Topological states in systems including interaction effects, such as electron-electron and electron-phonon, remain less explored. By exploiting a non-perturbative and non-adiabatic picture, here we show how the interaction between electrons and a coherent phonon mode can lead to a bandgap hosting edge states of topological origin. Further numerical simulations witness the robustness of these states against different types of disorder. Our results contribute to the search of topological states, in this case in a minimal Fock space.


I. INTRODUCTION
The search of topological states of matter is now reshaping condensed matter physics. 1,2 The pioneering works in the 1980s 3,4 bloomed about 20 years later with the prediction 5,6 and discovery of topological insulators in two 7 and three dimensions. 8 This field is today more active than ever, with new trends and discoveries expanding its frontiers. This includes, for example, the search for gapless but topological phases such as Weyl semimetals, 9 topological states induced by time-dependent fields (the so-called Floquet topological insulators), [10][11][12] time-dependent lattice distortions 13 and, more recently, topological states in non-Hermitian systems. 14,15 Today, topological states have a main role in the global search for means of achieving on-demand properties. 16 However, most of the current understanding of topological states remains at the level of a single-particle. The effect of interactions, both on topological phases predicted on the basis of a single-particle picture or as a mean of inducing new ones, stands out as a major problem. Previous studies along this direction have shown that electron-phonon interaction can either suppress 17 or even induce [18][19][20][21] non-trivial topological phases as the temperature increases. But the interaction between electrons and coherent phonons can also induce dressed states (even in the low-temperature limit), thereby requiring a careful analysis of the excitation spectrum of the composed system (electron and phonons). This type of interaction typically requires going beyond the adiabatic limit and has been predicted to lead to a phonon-induced bandgap opening in carbon nanotubes. 22,23 Experiments have also evidenced a breakdown of the Born-Oppenheimer approximation in graphene with the same type of high-symmetry optical phonons. 24 The recent observation of chiral phonons in 2D materials 25 also adds much interest in the context of possible ARPES experiments. 26 Furthermore, other authors have put forward the possibility of using optical means to control the electron-phonon interaction. 27 Here we examine a model for electron-phonon interaction in a quasi-one dimensional system and show that it may lead to robust topological edge states in a sample that otherwise lacked them. Specifically, we consider a graphene nanoribbon in the presence of a strong electron-phonon interaction with a single high-symmetry optical phonon mode. By exploiting a Fock space picture incorporating non-perturbative and nonadiabatic effects, 28,29 we find that at the center of the phonon induced bandgaps (located at half the phonon energy above the Dirac point) there are edge states of topological origin induced by the electron-phonon interaction. Furthermore, our numerical simulations show that these states remain robust to different types of disorder and ribbon geometries.

II. HAMILTONIAN MODEL AND FOCK SPACE SOLUTION SCHEME
To investigate the effects of the electron-phonon (e-ph) interaction, we use the framework introduced in Refs. [28] and [29]. The purpose of this section is to write the system's Hamiltonian in a basis for the electron-phonon Fock space corresponding to a single electron plus the excitations of the phonon mode. Since the description is coherent and as such the quantum phases are fully preserved, the e-ph interaction does not produce any phase randomization. This approach has been used for a variety of problems including vibration assisted tunneling in STM experiments, 30 transport through molecules 31,32 and resonant tunneling in double barrier heterostructures. 33 Let us consider a tight-binding description of graphene nanoribbons (GNRs) through the Hamiltonian where the sum runs over nearest neighbor carbon atoms and γ ij represents the hopping amplitude connecting them. The fermion operatorĉ † i (ĉ i ) creates (annihilates) an electron at site i of the lattice. Carbon displacements δr i from their respective equilibrium positions r 0 i are incorporated as a renormalization of the bare hopping amplitude γ 0 = 2.7 eV through 34 i − r 0 j + δr i − δr j | accounts for the distance between carbons i and j, a 0 = |r 0 i − r 0 j | 1.42Å is the equilibrium C-C distance, and b 3.37 is the rate of decay.
As depicted in Fig. 1(a), we consider a single phonon mode characterized by a rigid displacement δr = a 0 Qu between  sublattices A and B, where Q sets the strength of the displacement and u = (cos φ, sin φ) its direction. 35 The positions of the carbon atoms thus depend on which sublattice they belong, i.e. r i = r 0 i ± δr/2, for i ∈ {A, B}, respectively. Assuming small displacements, i.e. Q 1, we linearize the hoppings as γ ij = γ 0 (1 − bQ cos α ij ), where α ij is the angle subtended by the C-C bond and the displacement direction. Now we impose quantization on the mechanical coordinate, such that the above hoppings introduce the electron-phonon interaction.
The full Hamiltonian therefore readŝ whereâ † (â) creates (annihilates) one phonon excitation of frequency ω, and γ x sets the strength of the e-ph interaction. The pure electronic Hamiltonian is given by the C-C hoppings in equilibrium, and writes as in Eq. (1) but with the replacement γ ij → γ 0 . In Eq.
(2) we also included the phonon Hamiltonian, given byĤ We work within a Fock space spanned by |i, n = |i ⊗ |n states, where |i describes a single electron state (usually referred to the site basis), while |n sets the number of phonon excitations in the lattice. In this basis, the Hamiltonian of Eq.
(2) can be represented as the following matrix: Here, H 0 and H 1 are electronic matrices representing the γ 0 and γ x hoppings in the hexagonal graphene lattice, 1 ph and 1 el are the identity matrices in phonon and electron subspaces, respectively, while the remaining phonon matrices are defined as: Expressed in the Fock space basis, the Hamiltonian in Eq. (4) can be visualized as the original one without interactions together with the replicas corresponding to different phonon excitations and the interactions between them. The presence of these excitations is accounted for by the additional energies n ω. In this representation, the e-ph coupling enters through the last term, which enables the absorption and emission of a single phonon each time the electron 'hops' between two carbon atoms.

III. ELECTRON-PHONON INDUCED EDGE STATES
As introduced in the previous section, here we consider a fully quantized vibrational mode of frequency ω, and describe the vibrating nanoribbons through a Fock space spanned by the states |i, n , which accounts for both the electronic and vibrational degrees of freedom.
This might remind the reader of a similar picture used for time-periodic Hamiltonians: Floquet theory. 36 Indeed, there are several parallels stemming from a seeming isomorphism between the Floquet space and the Fock space, 37 but a few crucial differences must be noticed: (i) Unlike for the case of time-dependent potentials, for the case of phonons, temperature plays a natural role in defining the phonon population. (ii) In Floquet theory, the replica index is unbound while in the case of phonons described here it is bounded from below (n ≥ 0). (iii) The matrix elements for phonon emission and absorption change with the phonon population. Both descriptions (i.e. Fock space and Floquet space) match only when the system is in a highly excited state, a fact that is far from correct for optical phonons with typical energies exceeding k B T at room temperature. Aside from these differences, inspired from what we learned from Floquet topological states, 10,11,38 one might then search for similar physics induced by the electron-phonon interaction.
Another issue one might notice is the role of time-reversal symmetry (TRS): While a gap in irradiated graphene requires circularly polarized light so as to break TRS, 10,3940 the phonons considered here do not break such a symmetry (though they do open a bandgap 22,23,41 at ω/2 in the bulk material). However, one needs also to point out that in this work we are restricting ourselves to a ribbon geometry (i.e. a quasi one-dimensional system) rather than a two-dimensional system. Interestingly, chiral phonons, 25,42 lying at the corners of the Brillouin zone, could be used in two-dimensional hexagonal lattices to break TRS (at least locally in the valleys). Carrying on with the Fock-Floquet analogy, one could expect in the latter case similar physics as that of irradiated graphene with circularly polarized light.

A. Vibration induced bandgaps
The proposed vibration of the lattice consists in a single mode characterized by a rigid displacement between the two sublattices. We will work in the case where the displacement direction coincides with the longitudinal direction of the ribbon, i.e. φ = 0, motivated by the strong e-ph coupling observed in the optical mode A 1 (L) in CNTs leading to a Peierls-like mechanism [43][44][45] and Kohn anomalies; 46 and also in graphene samples. 24 To begin with, we consider a graphene nanoribbon with armchair edge geometry (ac-GNR). The reason of this particular choice rests in the possibility it offers to decompose the system into a series of decoupled eigenmodes. To do so, we start from the non-interacting case (γ x = 0) and we use the basis transformation proposed in Ref. [47]. This transformation takes the ac-GNR into a series of N y eigenmodes, each one consisting on a dimer chain with alternating hoppings γ (0) q,a = 2γ 0 cos[πq/(2N y + 1)] and γ (0) q,b = γ 0 , with q = 1, . . . , N y the mode number [see Fig. 1(b)]. In this sense, one can identify each eigenmode q as an independent Su-Schrieffer-Heeger (SSH) model, 48 for which the topological properties are well-known. [49][50][51] Furthermore, the SSH model for time-dependent hoppings was also investigated in the context of Floquet topological states, both theoretically 52-55 and experimentally. 56 When including the e-ph interaction, we can extend this mode decomposition in the Fock space, such that for each eigenmode we obtain a series of dimer chains (or replicas of the non-interacting case), each one belonging to a different number n of phonon excitations, see Fig. 1(c). This is easy to see regarding the structure of the Fock Hamiltonian in Eq. (4), where H 0 and H 1 commute with each other for this particular ribbon geometry and phonon mode. So, for a given mode q, the intrachain hoppings alternate between γ (0) q,a and γ (0) q,b , as in the non-interacting case. The phonon energy n ω in the nreplica enters as a site energy along the whole chain, according to the second term in the r.h.s. of Eq. (4); and the interchain hoppings connecting the n − 1 and n replicas are given by γ The degree of complexity imposed by the interaction clearly difficults the possibility of having analytic solutions for the system. However, the assumed weak coupling between the replicas (γ x γ 0 ) allows us to estimate the effects of the vibration on the electronic band structure by using perturbation theory around half the phonon energy. Other quantities like the local density of states (LDoS), eigenenergy spectrum and wavefunction amplitudes will be addressed numerically through standard techniques.
To begin with, we neglect the coupling γ x between the different phonon replicas, such that the (q, n)-dispersion relation at zeroth order can be written as: This equation represents the band dispersion one would obtain for a dimer chain with unit cell length a = 3a 0 /2, but shifted in an integer number of ω, which accounts for the phonons energy. Among the infinite number of replicas, we will take as reference the lowest energy bands belonging to the zerophonon subspace. This is justified for optical phonons with energies largely exceeding the thermal energy at room temperature (the stretching mode in graphene has a phonon energy of about 200 meV, this is about 8 times the thermal energy at 300K). Hence, the phonon population is zero for practical purposes. According to Eq. (6), the number of replicas crossing each other will depend on the relation between the phonon energy ω and the bandwidth associated with each replica, given by 6γ 0 for N y > 1.
For ω γ 0 , the distance between bands belonging to different replicas is much larger than their widths and, therefore, they do not cross each other. In this limit, the perturbation on the zero-phonon replica due to replicas with a higher number of phonons becomes negligible: The electron lying in a lattice without phonons can never reach enough energy as to spontaneously emit a phonon.
The symmetry between the conduction (+) and valence (−) bands with respect to n ω for n = 0 and n = 1 ensures that the band crossing takes place at half the phonon energy, ε = ω/2. This fact, together with the maximum allowed energy for the zero-phonon band, imposes the condition ω < 6γ 0 for the first band crossing between two different phonon replicas. In general, for smaller values of ω, the zero-phonon replica will cross with the n phonon replica once the condition n ω < 6γ 0 is fulfilled.
The eigenmode decomposition allows us to simplify the analysis around the band crossing processes. As in this particular geometry they are decoupled, the only relevant crossings are those with the same value of q. We can think, then, in the crossing between n = 0 and n = 1 bands belonging to the eigenmode q. To have such a crossing, it is necessary that the n = 0 conduction band and the n = 1 valence band overlap, which imposes the following range: If such inequality can not be met, then the bands do not cross each other, and the zero-phonon replica gets virtually unperturbed. Conversely, if such inequality is fulfilled, the bands will cross at the k-points determined by the condition ε q,− (k). This equation has two solutions, ±k * q , due to the symmetric dispersion of the bands in Eq. (6) around k = 0. When we include the e-ph interaction through γ x , the crossing between the bands gets avoided, yielding a gap induced by the vibration. For these k * q the group velocity goes to zero, meaning that a new backscattering process was introduced. From the point of view of the n = 0 replica, an electron traveling in a static ribbon 57 with energy ε (0) q,+ (k * q ) may suffer a reflection, together with the emission of a phonon of energy ω. Alternatively, from the point of view of the n = 1 replica, an electron traveling in a vibrating ribbon, such that the composite e-ph system has energy ε (1) q,− (k * q ), may also be reflected, after absorption of the phonon excitation present in the lattice.
Supposing the limit case N y → ∞, we can take the cosine argument in γ (0) q,a as a continuous variable x = qπ/(2N y + 1) in the range 0 < x < π/2. Under the perturbative regime, we estimate the bandgap size from a reduced Hamiltonian which only includes those bands which are expected to cross (see App. A). The size of the gap, as a function of x, writes: where η = ω/4γ 0 and η ± = η ± 1/2. In Fig. 2 we show this vibration induced bandgap for several values of the phonon frequency in the limit N y → ∞. To some extent, one can regard this limit as taking the nanoribbon into a twodimensional graphene layer. Here, the set of eigenmodes becomes dense, meaning that there is always a specific x * -value in which the vibration induced gap necessarily closes. This is a consequence of the preserved TRS by the phonon mode. Going back to the finite N y system, the set of x-values is no longer dense, and the overall gap will be given by the q-mode closest to x * . This somewhat relaxes the need to break TRS as to open a gap in quasi one-dimensional systems. We notice that Eq. (8) has physical meaning as long as the argument of the square root is positive. This sets the condition for those eigenmodes in which there is a band crossing, and it results to be 0 < x < acos(η − ) for η > 1/2 and acos(η + ) < x < π − acos(η − ) for η < 1/2, in agreement with Eq. (7). On the other hand, the bandgap size depends linearly on γ x , and there is also some dependence with the phonon frequency through the η parameter. Interestingly, while for semiconducting GNRs the bandgap around the central region ε = 0 can only be closed by increasing N y , the vibration induced bandgap could be controlled, to some extent, through the modulation of the e-ph coupling. 27 We can, in turn, determine the maximum value of the gap as a function of ω. This yields two regimes: (i) for η < √ 3/2, a frequency independent regime with maximum gap ∆ max = 3γ x , and (ii) for η > √ 3/2, a frequency dependent regime where the maximum gap decreases and it closes in the limit ω = 6γ 0 .

B. Spectral properties and characterization of the edge states
Throughout the following analysis we will consider a highfrequency regime by choosing ω = 5γ 0 . 58 For this value, the conduction band belonging to the zero-phonon replica only cross with the valence band of the one-phonon replica, since the condition n ω < 6γ 0 can only be fullfilled by n = 1. Although the chosen ω largely exceeds the typical phonon energy of the optical mode, this high-frequency regime, together with the weak e-ph coupling assumption (γ x γ 0 ), allows us to truncate the infinite Fock space in the first two phonon replicas, thereby simplifying the discussion of the vibration effects on the electronic properties of the ribbon. In particular, this ensures both the valence band (ε < 0) and the semiconducting gap (ε = 0) regions being unaffected by the vibration, at least in lowest order in the e-ph coupling. This allows us to establish a clear distinction between the well-known "native" edge-states 59 at ε = 0, and the expected e-ph induced edgestates, located at ε = ω/2. In any case, the same analysis can be carried out for ω 200 meV, but keeping in mind that a competition between the native topology and the e-ph induced edge states may occur near the charge neutrality point. 54 In Fig. 3(a) we show the LDoS weighted over the zerophonon replica (red shaded area), evaluated at the lateral border of a semi-infinite ac-GNR of N y = 11 carbon atoms wide (L y = 25.82Å). Except for some region around ε = ω/2, the shape of the LDoS is qualitatively the same as that of the non-interacting case (solid black), with a central peak related with the native edge states. This is expected as the e-ph interaction only becomes effective in the band crossing region at ω/2. Here, we observe a similar behavior as in the ε = 0 region, i.e. a depletion in the LDoS with a pronounced peak rising in its center.
To infer whether the peak at ω/2 survives far away from the border of the ribbon, we show in Fig. 3(b) the zero-phonon LDoS evaluated at the center of an infinite ac-GNR. The pronounced peaks at ε = 0 and ω/2 are no longer visible, and instead we can observe bandgap openings around these energies (blue arrows). For ω/2 this is the e-ph interaction induced bandgap, which was also predicted in vibrating CNTs. 22,23,41 As it is expressed in Eq. (8), the size of the vibration induced gap depends on the eigenmode we are looking at, which is observed in the LDoS maps of Figs. 3(c) and (d). For the chosen values ω = 5γ 0 and N y = 11, the band crossings occur for those eigenmodes fulfilling q < q max 5.29.
The minimum gap occurs for the mode with q closest to q max , which in this case is q = 5. Here, the effective interchain coupling between the phonon replicas [c.f. Eq. (A3)] is the smallest, and it increases for smaller q modes. The absence of peaks in the bulk LDoS makes us suspect that, as in the case of the native edge states at ε = 0, the peak at ω/2 is also related with states localized at the ribbon's border. Let us see, now, the behavior of the vibration induced peak as we move inside the ribbon. In Fig. 4(a) the zerophonon LDoS evaluated at half the phonon energy is shown as a function of the transversal layer number j [vertical lines of carbon atoms, see Fig. 1(a)]. The LDoS decays exponentially, such that for j ∼ 40 it becomes negligible. The way in which the peak decays is quite irregular, due to the superposition of the contributing tranversal modes (see inset plots in the figure).
To characterize the eigenenergy spectrum and the localized state wavefunctions we now consider a finite ac-GNR of N x = 200 carbon atoms long (L x = 211.58Å) and N y = 11 wide (L y = 25.82Å). In Fig. 4(b) we show the eigenenergy spectrum resolved in eigenmodes. To visualize such a spectrum as a perturbation of that in the non-interacting case, we weight each eigenstate over the n = 0 replica, i.e. the zero-phonon subspace. Given a k-eigenstate belonging to the eigenmode q, we can decompose it as the following superposition among the n replicas according to: and take its projection over the n = 0 subspace, i.e.
Although the number of replicas is infinite, the fact that we work in a perturbative regime allows us to truncate the full Fock space in a few replicas. In our case where ω = 5γ 0 and γ x γ 0 , the subspace associated with those replicas with n > 1 has negligible impact on the n = 0 replica. In Fig. 4(b), the red dots reveal how the valence band belonging to the n = 1 replica mixes with the conduction band of the n = 0 replica (black dots) for ε ∼ ω/2. Around this value, we can observe bandgap openings for those modes with q < q max , together with the presence of two degenerate midgap states per eigenmode (blue arrow).
In Fig. 4(c) we show one of the two midgap states for q = 4. As we move towards the center of the ribbon, the probability density oscillate according with the weights this particular qmode has on the sites composing the transversal layer. As it happens with the LDoS peak, the wavefunction also shows an exponential decay along the longitudinal direction, with a typical 'inverse gap' localization length. However, as the ribbon in this example has a finite length, the wavefunction has weight in the two borders. We can think of this state as bonding or antibonding combination of two states, |ψ L and |ψ R , localized at the left and the right border of the ribbon, respectively. For a ribbon with a large number of transversal layers (as it happens here), the overlap ψ L |ψ R can be neglected and one can consider |ψ L and |ψ R as linear combinations of the two degenerate midgap states.  We use a color scale to indicate the eigenstates weight with respect to the zero-phonon replica: Full weight (p q,k = 1) is in black, while zero weight (p q,k = 0) is in red. (c) Probability density |ψj| 2 (normalized to its maximum value and weighted over the n = 0 replica) as a function of the transversal layer j for one of the two eigenenergies at ε = ω/2 and eigenmode q = 4 in panel (b).

C. Topological origin of the edge states
As discussed before, when decreasing the phonon energy from ω = 6γ 0 on we will observe band crossings between different phonon replicas which, in turn, generate bandgap openings as new backscattering processes are introduced. This can be regarded as a band inversion process: In the crossing region, the valence band from the n = 1 replica happens to have less energy than that of the conduction band from the n = 0 replica. This band inversion is characteristic in topological phase transitions, together with the formation of localized midgap states. This strongly motivates the calculation of the Zak phase to infer about the topological nature of the vibration induced localized states. Although there is an infinite number of bands due to the structure of the Fock space, we can again truncate this by taking only those bands belonging to the n = 0 and n = 1 replicas. The corresponding topological invariant can be calculated as the integral of the Berry connection 60 with |u k,α the Bloch states belonging to the α-band and the integral taken over the first Brillouin zone. The bulk-boundary correspondence then allows us to characterize the existence of topological states with the Zak phase. Summing up Z α for all the bands with energy below a given gap yields the cumulative phase (modulo 2π), which indicates the existence (with cumulative phase π) or absence (zero cumulative phase) of topological midgap states. Although the Zak phase can be obtained analytically in the SSH model 51 and graphene nanoribbons, 59 we here proceed with a numerical calculation of the invariant. Computationally speaking, the Zak phase involves the calculation of wave-function amplitudes with some arbitrary gauge introduced by the diagonalization algorithm. In consequence, their possible outcomes, i.e. Z α = 0 or π, do not fully determine the band topology. However, one can infer its topology from the variation of Z α with respect to a reference case in which such a phase is known. In the previous sections, we concluded that for very high phonon energies, given by the condition ω > 6γ 0 , there are no band crossings, and in consequence the zero-phonon replica remains unperturbed. This high-frequency limit represents our reference case, where the Zak phase is well known for all band replicas.
We show in Table I the Zak phases corresponding to the four bands belonging to replicas n = 0 and n = 1 for the phonon energies ω = 5γ 0 and ω = 3γ 0 . As can be seen from the table, in this example the role of the interaction is to open a gap between bands 2 and 3, and adding a factor π to the band's cumulative phase. For ω = 5γ 0 this implies that modes with q = 1, . . . , 5 host vibration induced topological states between bands 2 and 3, while modes with q = 8, . . . , 11 host native topological states in the gap between bands 1 and 2 (related with the zero-phonon replica) and between bands 3 and 4 (related with the n = 1 phonon replica). For ω = 3γ 0 , the band crossing condition is fulfilled for q = 1, . . . , 9. For these modes, the Zak phase for the lowest energy band (valence band of the n = 0 replica) is equal to π, which means that native edge states are present. However, more interesting is the case of the q = 8, 9 modes, where the topological invariant for the second lowest energy band is equal to zero (cumulative Zak phase of π), which implies that these modes host both native and vibrational induced topological states.
In both cases we used γx = 0.1γ0. see next, this difference in the formation of the native and interaction-induced topological states brings with it important consequences when evaluating the robustness of such states against changes in the ribbon geometry and the introduction of several types of disorder.

D. Ribbon geometries and robustness against disorder
So far we have been discussing the effects of the e-ph interaction on a particular ribbon geometry where the eigenmodes remain decoupled even in the presence of the vibration. It therefore becomes natural to ask whether the topolog-ical states survive in other geometries and, in turn, if they are robust against different types of disorder which might couple these eigenmodes. In this section we provide an answer to these questions by analysing the LDoS along the ribbon edges for different geometries and by incorporating either vacancy or impurity disorder.
In Fig. 5 we show the LDoS evaluated around the edge region for different GNRs geometries. As in Fig. 3, we plot the quantity ln(1 + ρ) as to compensate the peaks height from the rest of the data. In these examples we considered randomlygenerated vacancy disorder of 0.5% over the complete sample. In all panels, we show the zero-phonon LDoS without disorder (solid black), a single disorder realization (shaded blue area), and an ensamble average over N = 200 realizations (red line). The insets in each panel show zoom regions around ε = 0 (bottom inset) and ε = ω/2 (top inset) to help visualization, while the schemes describe the type of ribbon geometry, characterized by longitudinal (cyan) and transversal (orange) borders. In all cases we used N x ∼ 1000 to ensure a ribbon length much larger than the localization lengths of all edge states.
Comparing the panels in Fig. 5 we first notice that, in absence of vacancy disorder, the LDoS peak in ε = 0 can be present or not depending on the ribbon geometry. This is illustrated by the black lines in the figure (see bottom inset in each panel), where panels (b) and (c) show no peak at this energy, while in panels (a) and (d) such a peak is certainly visible. The peak in ε = ω/2, on the other hand, is present in all panels. This important result allows for a clear distinction between the ε = 0 and the ε = ω/2 edge states. While the native states may appear or not depending on the particular ribbon geometry, the presence of topological states induced by the e-ph interaction is ensured by the crossing of conduction and valence bands belonging to different phonon replicas. For this reason we believe the native states are rather marginal: Although they admit a topological characterization through the Zak phase, their emergence is strictly determined by the D Ribbon geometries and robustness against disorder IV SUMMARY AND FINAL REMARKS.
When we include a 0.5% concentration of vacancies along the full sample, the main structure of the LDoS holds (blue lines), though it obviously displays a noisy pattern around the LDoS without disorder (black lines). Such a perturbative behavior can be tested by taking an ensamble average over several ribbon samples with the same vacancy concentration. The red lines show the LDoS averaged over N = 200 realizations, and superimpose the LDoS without disorder along almost the entire spectrum. However, looking closer at ε = 0 in Figs. 5(b) and (c) (see bottom insets) we notice that, in average, the zero energy peak returns when vacancy disorder is included. The fact that the average LDoS (red) in Fig. 5(b) shows a peak at ε = 0 while the trial LDoS for a single realization (blue) shows no peak, and that this peak is present in both cases in Fig. 5(c) indicates that this is an intermittent effect: In those cases where the LDoS shows no central peak without disorder, when including disorder this peak may appear for some vacancy configurations. In fact, these peaks have nothing to do with the native edge states appearing in Figs. 5(a) and (d), but these are related to localized states that surround the vacancies in the sample. 61 Thus depending on the presence (or absence) of vacancies nearby the region where the LDoS is being evaluated one can see (or not) a peak in the LDoS. In the considered examples of Fig. 5(b) and (c), the region where the LDoS was evaluated involves ∼100 carbon atoms, and since the vacancy concentration is 0.5%, one expects ∼0.5 vacancies in this region, so the chances of observing a peak in this region are one in two. Obviously, since the chances to have one or more localized states due to the presence of vacancies within the evaluation region grow with the vacancy concentration, we expect a simple relation between this quantity and the average height of the central peak. Having understood the role of disorder in the central peak of Figs. 5(b) and (c), we now observe that the shape and intensity of the central peak in Figs. 5(a) and (d) change little when including disorder, meaning that the native states, if present, are robust against moderate vacancy disorder. For the e-ph induced topological states we can arrive to the same conclusion, as the peaks in ω/2 are all the same, regardless of the vacancy disorder. We now investigate the role of impurity disorder on the ribbon's LDoS. This is modeled through a random variation of the on-site energies within the range [−W, W ]. This means that the pure electronic Hamiltonian in Eq. (2) is replaced bŷ with −W ≤ i ≤ W the random on-site energy. In Fig. 6 we show the LDoS around the e-ph band crossing point, centered at ε = ω/2, together with the LDoS around the Dirac point at ε = 0 (insets).
In these examples we evaluate the zero-phonon replica LDoS over 40 transversal layers of carbon atoms for both left and right borders, considered as mirror images each other. This was done for three ribbon geometries (see schemes in each panel), and we used W = 0.01γ 0 (red), 0.02γ 0 (blue), and 0.03γ 0 (green).
In these examples we calculated the average LDoS over 250 disorder realizations. Black lines exhibit the case without disorder (W = 0), while dotted orange lines (shaded area) illustrate the case of a single disorder realization for W = 0.03γ 0 . For the shown ribbon geometries, we can see that the e-ph induced LDoS peak at ω/2 (dotted orange) now splits out into several peaks lying within the bandgap region. To understand why this type of disorder produces such an effect, first notice that the LDoS peak at ω/2 without disorder (black lines) can be decomposed into several peaks, each one belonging to a topological edge state. The position of these peaks depends on the average of the on-site energies around the region where the probability density is finite. If we imagine the impurity disorder [first term in Eq. (12)] as a perturbation, and ψ k (r i ) represents the wavefunction amplitude of the (unperturbed) e-ph induced topological state k at position r i , then the energy ε k (and thus the peak position in the LDoS) will depend on the on-site energies as: which can be interpreted as the original position (i.e. without disorder), plus the k-state weighted average of the on-site energies. As for the considered ribbon sizes the topological edge states have finite weight over a small number of sites, the last term in the above equation may not vanish in general.
In fact, this quantity tends to increase with the degree of disorder W . This is reflected as a broadening of the averaged LDoS peaks (solid red, blue, and green) when W increases. Importantly, the area below the LDoS peak remains always constant, meaning that the number of topological states in the sample is independent of disorder. Of course, though W does not change the number of topological states, for larger W values these states may be located in energy regions outside the overall gap, thus difficulting a clear separation between localized and extended states. The same disorder-induced peak broadending can be observed in the central region around ε = 0 (see insets), but it is important to notice that, as in the vacancy disorder case, the native states can be present or not depending on the ribbon geometry. Additionally, due to the zoom factor, one can observe that the average peaks are not perfectly centered at ε = 0, but slightly shifted to the left. This is not related to disorder (the black line peaks are also shifted) but a second order effect in the coupling between the zero-and one-phonon replica bands.

IV. SUMMARY AND FINAL REMARKS.
In summary, we have shown that novel and robust states of topological origin form as a consequence of the electronphonon interaction in graphene nanoribbons. This study, based on a specific model for the electron-phonon interaction given by a stretching mode in graphene nanoribbons, serves as a proof of concept. The topological states form at the center of a bandgap (induced by the same interaction) located at half the phonon energy above the charge neutrality point. This was confirmed in several ribbon geometries and for vacancy and impurity disorder configurations. While both the native and the e-ph induced states were characterized through the Zak phase, and shown to be robust against disorder, the native states only appear in some specific ribbon geometries. Conversely, for a non-negligible e-ph interaction, the presence of the vibration induced topological states is guaranteed as long as the phonon energy does not exceed the typical band width, i.e. ω < 6γ 0 . Such a condition provides the required band inversion between the first two phonon replicas. Interestingly, this physics happens in our case even when the phonons do not break time-reversal symmetry, similar to a driven one-dimensional topological insulator. 54 In two dimensional systems, however, the e-ph induced bandgap finally closes for some particular mode k (see discussion on the N y → ∞ limit around Fig. 2) and one should break TRS to restore the gap (as it is required for light to induce Floquet topological states in the same material). In this sense, the recent observation of chiral phonons in two-dimensional materials 25 may open a promising way for studying electron-phonon induced topological phase transitions. the energy difference between the two replicas is ω. Let us suppose that ω is small enough such that the n = 0 conduction band and the n = 1 valence band cross. If k denotes the Bloch quasimomentum, there are two k-values (±k * ) where the band crossing occurs. If we focus on one of these points, say k * , we have two k-states at the same energy, each one belonging to one of the two replicas. This degeneration gets removed by the electron-phonon interaction, which in our case corresponds to the coupling between the replicas, and yields the bandgap opening. In this appendix we estimate the gap size produced by the e-ph interaction.
When using the eigenmode decomposition, the hopping term between sites depends on the mode q we are looking at, with a factor cos[qπ/(2N y + 1)], and q = 1, . . . , N y . We take the cosine argument as a continuous variable x, within the range 0 < x < π/2, and we simplify this analysis by truncating the full Fock space so we only keep the n = 0 and n = 1 replicas. Considering the bulk situation (i.e. an infite long dimer chain) we can use Bloch theorem and obtain the following Hamiltonian: where v (n) with a = 3a 0 /2 the unit cell length, and the bar standing for complex conjugation. Since we are only interested in the n = 0 conduction and n = 1 valence bands, we can reduce even more this Hamiltonian. To do so, we first diagonalize the 2 × 2 block matrices in the diagonal. As they commute each other, we can diagonalize the n = 0 block, and obtain the energies ε = ±|v (0) q |. Similarly, for the n = 1 block we obtain ε = ω ± |v (0) q |. The next step is to write the Bloch Hamiltonian in this new basis, such that it allows the proper truncatioñ where ∆ϕ = ϕ q . This Hamiltonian has the following eigenvalues (A4) With these expressions, we found the new eigenenergies for k close to the band crossing point. We can, in turn, particularize to the point in which this band crossing occurs and obtain an expression for the size of the gap ∆(x). Recalling that the band crossing takes place at ε = ω/2, we take |v (0) q | equal to this energy, and using the definitions for v (n) q given in Eq. (A2), we obtain ∆(x) = 2 9γ 2 x cos 2 x sin 2 ka 1 + 4 cos 2 x + 4 cos x cos ka . (A5) By finding the k-value for which |v (0) q | = ω/2 is satisfied, and defining the adimensional parameter η = ω/4γ 0 , we obtain that the size of the e-ph interaction induced bandgap is given by where η ± = η ± 1/2. The size of the gap will depend, therefore, on this parameter η and the particular eigenmode q (through the x-variable) we are looking at. The square root argument defines the band crossing regimes, as this quantity needs to be always positive. This yields the condition which is equivalent to Eq. (7) of the main text.