Chiral symmetry breaking and topological charge of zigzag graphene nanoribbons

Non-interacting zigzag graphene nanoribbons [1, 2] exhibit chiral symmetry and are topological insulators [3]. However, this symmetry breaks down in the presence of electron interactions, leading to uncertainty regarding the definition of a topological charge. Our paper aims to elucidate the formation of a topological charge in disorder-free interacting zigzag graphene nanoribbons. We will demonstrate that in the presence of such interactions, the relevant symmetry combines mirror and time-reversal operations. The breakdown of chiral symmetry leads us to utilize axial current instead of ordinary current for an accurate definition of a topological charge, mirroring the approach by Goldstone and Wilczek [4] in the context of solitons in polyacetylene [5]. This issue holds experimental significance due to the recent fabrication of atomically precise nanoribbons (Ruffieux et al [6], Kolmer et al [7], Houtsma et al [8]).

One way to define a topological charge is through the generation of mass (or gap), which breaks chiral symmetry [4]. This concept can be elucidated using a one-dimensional Dirac equation. In one spatial dimension (and generally for all odd spatial dimensions), the gapless (or massless in the particle physics context) Dirac fermions possess the property of chirality. The Lagrangian density of the massless Dirac fermion (see the appendix for details and conventions) is invariant under two kinds of global gauge transformations:

$$\begin{equation} \psi \to e^{i \Lambda} \psi \;\; \text{ordinary}, \quad \psi \to e^{i \gamma^5 \Lambda} \psi \;\;\text{axial or chiral}, \end{equation} \tag{ 1 }$$

where ψ is a Dirac spinor, and $\gamma^5 = \gamma^0 \gamma^1$. Chirality is defined to be the eigenvalue of γ⁵, which can take only ±1 values. According to classical Noether's theorem, these invariances imply the conservation of the corresponding currents (Fermi velocity may be multiplied to the spatial component):

$$\begin{equation} j^{\mu}_{\mathrm{ordinary}} = \bar{\psi} \gamma^\mu \psi, \quad j^{\mu 5}_{\mathrm{axial}} = \bar{\psi} \gamma^\mu \gamma^5 \psi. \end{equation} \tag{ 2 }$$

Owing to a special property of gamma matrices in two dimensions, the above two currents can be related by:

$$\begin{equation} j^{\mu 5} = - \epsilon^{\mu \nu} j_\nu, \end{equation} \tag{ 3 }$$

where ε^µν is the two-dimensional Levi–Civita symbol with $\epsilon^{01} = 1$. Expressing the Dirac spinor ψ in γ⁵-diagonal basis as $(\psi_R,\psi_L)^\textrm{t}$ ('t' is transpose), the explicit expressions of currents are given by (note that $j_1 = -j^1$, etc):

$$\begin{equation} j^0 = \psi^{\dagger}_R \psi_R + \psi^{\dagger}_L \psi_L = j^{1,5} , \quad j^1 = \psi^{\dagger}_R \psi_R - \psi^{\dagger}_L \psi_L = j^{0,5}. \end{equation} \tag{ 4 }$$

A mass term for Dirac fermion couples ψ_R and ψ_L . Two independent mass terms (in Hamiltonian) are allowed in the following way:

$$\begin{equation} m \left(\psi^{\dagger}_R \psi_L + \text{h.c.} \right) + i m ^{^{\prime}} \left( \psi^{\dagger}_R \psi_L - \text{h.c.} \right). \end{equation} \tag{ 5 }$$

The above mass terms break the chiral symmetry explicitly. From the viewpoint of a one-dimensional polyacetylene system, these masses correspond to the even and the odd components of dimerization. In the presence of the mass terms, the continuity equation for the axial current is modified as follows [4, 9]:

$$\begin{equation} \partial_\mu j^{\mu 5} = \frac{\partial j^1}{\partial t}+ \frac{\partial j^0}{\partial x} = - 2 \left(m- i m^{^{\prime}}\right) \psi_R^{\dagger} \psi_L - 2 \left(m + i m^{^{\prime}}\right) \psi_L^{\dagger} \psi_R \equiv \mathcal{M}. \end{equation} \tag{ 6 }$$

Equation (6) remains valid even if the mass parameters depend on coordinates. For slow adiabatic temporal changes, the time derivative of equation (6) can be neglected, and we have ($\rho = j^0$ denotes charge density).

$$\begin{equation} \frac{\partial\rho }{\partial x} = \mathcal{M}\left(x\right). \end{equation} \tag{ 7 }$$

An abrupt step-like change of the mass terms can give rise to a solitonic [5, 10] topological charge [4] accumulated at the soliton location, where the spectral weight of a soliton is equally divided between conduction and valence bands. More precisely, $\mathcal{M}(x)$ in equation (6) corresponds to the topological charge [11, 12] density after subtracting out the trivial vacuum contribution [4]. Then, the integration of equation (7) over some interval with endpoints (say, $[-\infty,\infty]$) gives the topological charge $Q_\textrm{top}$ associated with the solitonic configuration of mass parameters

$$\begin{equation} \rho\left(\infty\right)-\rho\left(-\infty\right) = \int_{-\infty}^{\infty} \textrm{d}x \; \mathcal{M}\left(x\right) = Q_\textrm{top}. \end{equation} \tag{ 8 }$$

The determination of this topological charge relies on the topological aspects of the mass parameters, defined as $\mathcal{M}(x) = \frac{\partial \theta}{\partial x}$, where the order parameter $\theta(x)$ represents a kink:

$$\begin{equation} \theta\left(x\right) = \begin{cases} -1/2 & \text{for }x<0\\ +1/2 & \text{for }x>0. \end{cases} \end{equation} \tag{ 9 }$$

In polyacetylene $\theta(x)$ represents dimerization order [5]. However, this topological charge is not necessarily quantized in general circumstances [4]. Then, by equation (8), the induced fermion number can be identified with the topological charge $Q_\textrm{top}$ of mass parameters, which is in turn determined by their boundary behavior [10].

Our goal in this paper is to try to understand the formation of a topological charge in disorder-free interacting zigzag graphene nanoribbons in view of the above statement. First, let us give a brief introduction to non-interacting graphene nanoribbons. The band structure of a non-interacting zigzag graphene nanoribbon is shown in figure 1. Zigzag graphene nanoribbons [2, 13–15] consist of A and B sublattices, as seen in the top figure of figure 2. The A/B sublattice degree of freedom is referred to as the chirality. The chiral transformation [16, 17] on electron operators is implemented by

$$\begin{equation} \widehat{\mathsf{a}}_i\rightarrow \widehat{\mathsf{a}}_i, \quad \widehat{\mathsf{b}}_i\rightarrow -\widehat{\mathsf{b}}_i, \end{equation} \tag{ 10 }$$

where $\widehat{\mathsf{a}}_i$ and $\widehat{\mathsf{b}}_i$ are the electron destruction operators of A and B sublattice sites which are similar to the Dirac fermion $\psi_R, \psi_L$ under the action of γ⁵ (see the appendix). Under the transformation equation (10), the nearest-neighbor hopping term between A–B sublattices of the tight-binding Hamiltonian of graphene nanoribbon without on-site interaction changes sign, and this is called the chiral symmetry. This chiral symmetry implies a particle-hole symmetry within the energy spectrum and ensures the existence of zero-energy soliton [18, 19] edge states. Thus, non-interacting graphene zigzag nanoribbons are deemed a symmetry-protected topological insulator [20, 21].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Now, the inclusion of on-site repulsion interactions among electrons explicitly breaks chiral symmetry, as the density operator does not alter its sign under the transformation described by equation (10). This interaction introduces an energy gap, Δ. Despite this energy gap, the persistence of edge states is evident, as depicted in figure 2. Furthermore, interacting zigzag graphene nanoribbons manifest as short-range entangled symmetry-protected topological insulators [21]. (Here, we exclude consideration of disordered zigzag ribbons that exhibit long-range entangled [22, 23] topological order [21].) Consequently, it is natural to inquire about the topological origin of these edge states in the presence of on-site interactions. (In this paper, we do not consider armchair ribbons since the physical origin of a gap may differ.)

We aim to address the aforementioned questions by investigating the interacting nanoribbons within the framework of the Hartree–Fock (HF) approach and employing an anomalous continuity equation similar to equation (6). The self-consistent HF fields break chiral symmetry. However, the self-consistent fields of the opposite zigzag edges exhibit antisymmetry under mirror symmetry operation. The relevant symmetry combines mirror and time-reversal operations, enabling the formation of a topological charge induced by the mass on ferromagnetic zigzag edges in undoped ribbons. This mechanism is analogous to the equation in equation (8) and is of a non-perturbative nature. This topological charge persists even in the low-doped region where the ground state displays edge spin density waves. It appears as a nearly zero energy edge mode. Our work demonstrates that interacting graphene nanoribbons retain their status as symmetry-protected topological insulators. We believe that the conventional classification scheme of topological insulators [24, 25] does not directly apply to zigzag graphene nanoribbons because interacting zigzag ribbons display gapful edge excitations, whereas the conventional classification applies to gapless excitations. Moreover, the quasi-one-dimensionality of zigzag ribbons is important as it spatially separates the opposite zigzag edges of different chiralities.

Let us consider the narrowest armchair nanoribbon with zigzag edges and the narrowest zigzag nanoribbon, as illustrated in the top figure of figure 2. The Hamiltonian of the tight-binding Hubbard model of these nanoribbons [26–29] consists of the hopping term with hopping amplitude t and the Hubbard term for on-site Coulomb repulsion U

$$\begin{equation} \widehat{\mathsf{H}} = \widehat{\mathsf{H}}_{t} + \widehat{\mathsf{H}}_{U}. \end{equation} \tag{ 11 }$$

The second quantized hopping Hamiltonian is (h.c. is the Hermitian conjugate, and $\sigma = \uparrow,\downarrow$ is spin)

$$\begin{equation} \widehat{\mathsf{H}}_{t} = - t \sum_{< \mathbf{R}, \mathbf{R}^{^{\prime}} > \sigma} \left( \widehat{\mathsf{c}}^{\dagger}_{\mathbf{R} \sigma } \widehat{\mathsf{c}}_{\mathbf{R}^{^{\prime}} \sigma} + \mathrm{h.c.} \right), \end{equation} \tag{ 12 }$$

where $\widehat{\mathsf{c}}_{\mathbf{R},\sigma}$ is the electron destruction operator at site R with spin σ and $\lt \mathbf{R}, \mathbf{R}^{^{\prime}} \gt$ denotes a nearest-neighbor pair. One of the pair belongs to A-sublattice, while the other is in B-sublattice. The Hubbard term is given by

$$\begin{equation} \widehat{\mathsf{H}}_{U} = U \sum_{\mathbf{R} } \widehat{\mathsf{n}}_{\mathbf{R} \uparrow} \widehat{\mathsf{n}}_{\mathbf{R} \downarrow}, \end{equation} \tag{ 13 }$$

where $\widehat{\mathsf{n}}_{\mathbf{R} \sigma} = \widehat{\mathsf{c}}^{\dagger}_{\mathbf{R} \sigma} \widehat{\mathsf{c}}_{\mathbf{R} \sigma}$ is the electron number operator at site R with spin σ. As discussed in section 1, $\widehat{\mathsf{H}}_t$ posssses the chiral symmetry, while the $\widehat{\mathsf{H}}_{U}$ explicitly breaks it.

Applying HF mean-field theory to equation (13) at half-filling, a quadratic mean-field Hamiltonian $\widehat{\mathsf{H}}_{\text{HF}}$ can be obtained:

$$\begin{equation} \widehat{\mathsf{H}}_{\text{HF}} = \widehat{\mathsf{H}}_{t} + \widehat{\mathsf{H}}_{U,\textrm{HF}}, \quad \widehat{\mathsf{H}}_{U,\textrm{HF}} = \sum_{\mathbf{R} \sigma} \, \Delta_{\mathbf{R} \sigma} \widehat{\mathsf{n}}_{\mathbf{R} \sigma}, \end{equation} \tag{ 14 }$$

where the gap parameter $\Delta_{\mathbf{R} \sigma}$, which corresponds to the mass parameter of equation (6), is defined by

$$\begin{equation} \Delta_{\mathbf{R} \sigma} = U \left( \langle \widehat{\mathsf{n}}_{\mathbf{R} \bar{\sigma}} \rangle -\frac{1}{2} \right),\quad \bar{\sigma} = \text{the opposite spin of } \; \sigma. \end{equation} \tag{ 15 }$$

A constant term has been added to $\widehat{\mathsf{H}}_{\text{HF}}$ to make the Fermi energy to be zero. The site gap $\Delta_{\mathbf{R} \overline{\sigma}}$ associated with spin $\overline{\sigma}$ can be understood as the staggered potential energy for an electron with the opposite spin σ. This potential is staggered since it changes sign between A- and B-sites, as depicted in figure 4. It is important to note that the mean-field HF Hamiltonian can be split into the spin-up and the spin-down components [30]:

$$\begin{equation} \widehat{\mathsf{H}}_{\text{HF}} = \widehat{\mathsf{H}}_{{\text{HF}} \uparrow} + \widehat{\mathsf{H}}_{{\text{HF}} \downarrow}. \end{equation} \tag{ 16 }$$

The band structure without the Hubbard interaction $\widehat{\mathsf{H}}_U$ taken into account was shown in figure 1. The nearly zero-energy edge states exist close to the Brillouin zone boundary, and they exhibit non-magnetic behavior (the net spin polarization of a zigzag edge is zero). The existence of these zero energy chiral edge states crucially rests on the chiral symmetry. However, a self-consistent solution of the HF mean-field Hamiltonian reveals that for ribbon widths shorter than 100 Å, the ground state is antiferromagnetic in the undoped case. Local density approximation calculations [27, 31, 32] also yield similar results. In this state, site spins of sublattice A (B) are spin-up (-down), and the state is doubly degenerate with an energy gap. The chiral edge states do exist, and these chiral edge states are now oppositely polarized (edge antiferromagnetism) [26]. Each zigzag edge displays ferromagnetism but the opposite zigzag edges with different chirality are antiferromagnetically coupled (see the top figure of figure 2). The band structures are displayed in the second figure in figure 2. The HF results qualitatively agree with those of density matrix renormalization methods [29], both in doped and undoped zigzag ribbons. In the low-doped regime, it is found that edge spin density waves develop in contrast to the ferromagnetic edges observed in the undoped case. Thus, doping is a singular effect. A staggered potential may be applied to graphene nanoribbons [33]. In this case, mirror and time reversal symmetries are broken, as is their combined operation. The HF method shows that the charge of each zigzag edge may assume non-integer values instead of integer values.

Ribbons exhibit a symmetry based on the combined operation of parity transformation (mirror reflection) and time reversal transformations, as described by the relationship:

$$\begin{equation} \langle \widehat{\mathsf{n}}_{\mathbf{R} \sigma} \rangle = \langle \widehat{\mathsf{n}}_{\overline{\mathbf{R}} \bar{\sigma}} \rangle. \end{equation} \tag{ 17 }$$

This relation is illustrated in figure 3. Note that $\overline{\mathbf{R}}$ is the mirror symmetric site of R and $\bar{\sigma}$ is the opposite spin of σ. Here, the relevant mirror axis of the zigzag nanoribbon is the horizontal x-axis and that of the armchair nanoribbon the vertical y-axis through the center of the ribbon. This symmetry aligns with the existence of a spin-up edge state on one zigzag edge and the corresponding spin-down edge state on the opposite zigzag edge.

Zoom In Zoom Out Reset image size

Then the definition of the gap parameter equation (15) immediately implies that

$$\begin{equation} \Delta_{\mathbf{R} \sigma} = U \left( \langle \widehat{\mathsf{n}}_{\mathbf{R} \bar{\sigma}} \rangle-\frac{1}{2}\right) = U \left( \langle \widehat{\mathsf{n}}_{\overline{\mathbf{R}} \sigma} \rangle-\frac{1}{2} \right) = \Delta_{\overline{\mathbf{R}} \bar{\sigma}}. \end{equation} \tag{ 18 }$$

Note that equation (18) is valid even for the low-doped case as long as equation (17) is valid. It is important to note that this mirror symmetry is accompanied by spin reversal and that the symmetry in the occupation numbers involves a spin reversal but not in the site gap values. At half-filling (undoped case), the condition $\langle \widehat{\mathsf{n}}_{\mathbf{R} \sigma} \rangle + \langle \widehat{\mathsf{n}}_{\mathbf{R} \bar{\sigma}} \rangle = 1$ being combined with equation (17) yields

$$\begin{equation} \Delta_{\mathbf{R} \sigma} = - \Delta_{\overline{\mathbf{R}} \sigma}, \end{equation} \tag{ 19 }$$

which is demonstrated in figure 4.

Zoom In Zoom Out Reset image size

Let us consider zigzag ribbons in the low doping limit. It turns out that doping has a singular effect [29]: the grounds state displays an edge spin density wave, as shown in figure 5. This is in contrast to the uniform edge ferromagnetism of the half-filled case. Note that in the low-doped region, the electron spins are collinear to a very good approximation, as confirmed by the density matrix renormalization group approach in the matrix product representation. This implies that our Hamiltonian at half-filling may still be applicable [29].

Zoom In Zoom Out Reset image size

To understand the role of the doping we write the filling condition as

$$\begin{equation} \langle \widehat{\mathsf{n}}_{\mathbf{R} \sigma} \rangle + \langle \widehat{\mathsf{n}}_{\mathbf{R} \bar{\sigma}} \rangle = 1+\delta, \end{equation} \tag{ 20 }$$

where δ is the doping level. Then we have

$$\begin{equation} \Delta_{\mathbf{R} \sigma} = U \left( \langle \widehat{\mathsf{n}}_{\mathbf{R} \bar{\sigma}} \rangle-\frac{1}{2} \right) = U \left(1+\delta - \langle \widehat{\mathsf{n}}_{\mathbf{R} \sigma} \rangle-\frac{1}{2} \right) = - \Delta_{\mathbf{R} \bar{\sigma}}+ U \delta = - \Delta_{\overline{\mathbf{R}} \sigma}+U \delta. \end{equation} \tag{ 21 }$$

Equation (21) shows that the antisymmetry of the gap parameter $\Delta_{\mathbf{R} \sigma}$ observed at half-filling in figure 4 is disrupted when moving away from the half-filling state relative to the mirror axis.

The total current operator of the tight-binding Hamiltonian equation (11) is given by

$$\begin{equation} \mathbf{J}_{\text{tot}} = i t \sum_{ \lt \mathbf{R}, \mathbf{R}^{^{\prime}}>, \sigma} \left(\mathbf{R}-\mathbf{R}^{^{\prime}}\right) \, \left(\widehat{\mathsf{c}}^{\dagger}_{\mathbf{R}, \sigma} \widehat{\mathsf{c}}_{\mathbf{R}^{^{\prime}}, \sigma} - \text{h.c.} \right), \end{equation} \tag{ 22 }$$

which is defined such that the ordinary continuity equation is satisfied. Remember, in one-dimensional Dirac fermions, the axial current can be defined, where the roles of charge density ρ and spatial charge current density j interchange between ordinary and axial current (as shown in equation (2)). Additionally, the Frölich equation governing the (unpinned) sliding charge density wave [34] is expressed as:

$$\begin{equation} \nabla \rho + \frac{1}{v_F^2} \frac{\partial j}{\partial t} = \kappa E. \end{equation} \tag{ 23 }$$

Here, κ represents compressibility, v_F signifies Fermi velocity, and E denotes an electric field. In one-dimensional systems, equation (23) essentially embodies the well-known axial anomaly [9]. These insights lead us to consider the time derivative of the spatial current, as outlined in equation (22). This approach is motivated by the pronounced manifestation of chiral (axial) symmetry breaking effects, notably evident in the continuity equation for the axial current (equation (6)).

The time derivative can be found from the commutator with $\widehat{\mathsf{H}}_{\text{HF}}$ equation (14):

$$\begin{equation} \frac{\partial \mathbf{J}_{\text{tot}}}{\partial t} = \frac{i}{\hbar}\left[\widehat{\mathsf{H}}_{U,\text{HF}},\mathbf{J}_{\text{tot}}\right] +\frac{i}{\hbar}\left[\widehat{\mathsf{H}}_{t},\mathbf{J}_{\text{tot}}\right]. \end{equation} \tag{ 24 }$$

This expression comprises two contributions: one from the HF decomposed Hubbard interaction and the other from the hopping term. Taking the HF ground state expectation value of equation (24) and considering the temporal adiabatic limit, $(-i \hbar) \frac{\partial \langle \mathbf{J}_{\text{tot}} \rangle}{\partial t} \approx 0$ we find :

$$\begin{equation} \left(-i \hbar\right) \frac{\partial \langle \mathbf{J}_{\text{tot}} \rangle } {\partial t} = \Big\langle \left[\widehat{\mathsf{H}}_{U,\text{HF}},\mathbf{J}_{\text{tot}}\right] \Big\rangle + \Big\langle \left[\widehat{\mathsf{H}}_{t},\mathbf{J}_{\text{tot}}\right] \Big\rangle\approx 0. \end{equation} \tag{ 25 }$$

Then, we explore the implications of equation (25) in light of the (integrated) anomalous continuity equation (equation (6)) and the accumulation of the topological charge (equation (8)).

The computation of the commutators are straightforward. For each spin projection, we find

$$\begin{equation} \left[\widehat{\mathsf{H}}_{U,\text{HF}},\mathbf{J}_{\text{tot},\sigma}\right] = \left(+i\right) t \sum_{< \mathbf{R}, \mathbf{R}^{^{\prime}} > } \left(\mathbf{R}-\mathbf{R}^{^{\prime}}\right) \left( \Delta_{\mathbf{R} \sigma}-\Delta_{\mathbf{R}^{^{\prime}} \sigma} \right) \left(\widehat{\mathsf{c}}^{\dagger}_{\mathbf{R}, \sigma} \widehat{\mathsf{c}}_{\mathbf{R}^{^{\prime}} \sigma}+ \text{h.c.} \right), \end{equation} \tag{ 26 }$$

and

$$\begin{equation} \left[\widehat{\mathsf{H}}_{t},\mathbf{J}_{\text{tot},\sigma}\right] = \left(-i\right) t^2\sum_{ \lt \mathbf{R}_1,\mathbf{R}>} \sum_{ \lt \mathbf{R}, \mathbf{R}^{^{\prime}}>} \left(\mathbf{R}-\mathbf{R}^{^{\prime}}\right) \left( \widehat{\mathsf{c}}^{\dagger}_{\mathbf{R}_1 \sigma} \widehat{\mathsf{c}}_{\mathbf{R}^{^{\prime}} \sigma}+\text{h.c.} \right). \end{equation} \tag{ 27 }$$

It's important to clarify that in equation (27), R₁ and $\mathbf{R}^{^{\prime}}$ are not a nearest neighbor pair, see figure 6. To elaborate further, if we consider a fixed $\mathbf{R}^{^{\prime}}$, then R must be a nearest neighbor in the other sublattice, thereby determining the bond vector $\mathbf{R} - \mathbf{R}^{^{\prime}}$. Subsequently, R₁ should be a nearest neighbor of R, placing it in the same sublattice as $\mathbf{R}^{^{\prime}}$. Now, considering the proximity of $\mathbf{R}^{^{\prime}}$ to R, R₁ can be $\mathbf{R}^{^{\prime}}$, resulting in $\widehat{\mathsf{c}}^{\dagger}\textrm{ger}_{\mathbf{R}^{^{\prime}}\sigma} \widehat{\mathsf{c}}_{\mathbf{R}^{^{\prime}} \sigma} + \text{h.c.}$, effectively representing two times the number operator at $\mathbf{R}^{^{\prime}}$. However, it is crucial to note that R₁ might not necessarily be the same as $\mathbf{R}^{^{\prime}}$; in such cases, R₁ becomes a next-nearest neighbor site from $\mathbf{R}^{^{\prime}}$ while still being in the same sublattice. Thus, equation (27) can be expressed as

$$\begin{equation} \left[\widehat{\mathsf{H}}_{t},\mathbf{J}_{\text{tot},\sigma}\right] = \left(-i\right) 2 t^2 \sum_{\mathbf{R}^{^{\prime}}} \, \sum_{ \lt \mathbf{R}, \mathbf{R}^{^{\prime}}>} \left(\mathbf{R}-\mathbf{R}^{^{\prime}}\right) \widehat{\mathsf{n}}_{\mathbf{R}^{^{\prime}} \sigma} + \text{next-nearest bond terms}. \end{equation} \tag{ 28 }$$

Zoom In Zoom Out Reset image size

Some examples of the next-nearest bond terms can be visualized using figure 7, such as $\langle\widehat{\mathsf{c}}^{\dagger}_{l+1,3,\uparrow} \widehat{\mathsf{c}}_{l,1,\uparrow}\rangle$. We note that in the framework of the anomalous continuity equation treated here, no bonds longer than the next-nearest ones appear.

Zoom In Zoom Out Reset image size

It is important to emphasize that these results describe two-dimensional electron motion, covering the charge movement both along and perpendicular to the ribbon. In the following section we apply the above results to the zigzag nanoribbons depicted in the top figure of figure 2.

In consideration of the anomalous continuity equation of the axial current (equation (6)), we must compute the expectation value of the gap term (equation (26)) and the terms originating from the hopping Hamiltonian (equation (28)) concerning the HF ground state. Recalling that the HF Hamiltonian is split into the spin-up and the spin-down components (see equation (16)) the expectation values can be computed for each spin component separately.

Now it is crucial to note that there is a vertical mirror axis through the center of the zigzag ribbon (see figure 4). Then the expectation values of the bond operators of equations (26) and (28) should respect this mirror symmetry. We hasten to comment that this mirror symmetry operation does not involve spin reversal. However, the mirror symmetry employed in the discussion of the gap parameter in section 3 does involve spin reversal. Therefore, these two mirror symmetry operations are distinct from each other.

This vertical mirror symmetry implies that the horizontal x-component of the expectation values of the operators of equations (26) and (28) vanish identically owing to the oddness of the horizontal component of the bond vector $\mathbf{R}-\mathbf{R}^{^{\prime}}$ of equations (26) and (28) under the mirror reflection, so that the topological charges can be accumulated only along the vertical direction: the upper and the lower chiral zigzag edges. Below $l,m$ denote the unit cell and the position within the unit cell, respectively (see figure 7).

For simplicity, we will evaluate the axial current for the smallest-width zigzag nanoribbons. Subsequently, we will validate this result for larger widths. This model of zigzag nanoribbons will clearly reveal the main physics. Firstly, regarding the expectation value of the hopping Hamiltonian in equation (28), it comprises two distinct components. Initially, we focus on the first term, the contribution stemming from the density operator terms. We observe that the contributions from m = 2 and m = 3 sites (refer to figure 7) vanish. This occurs due to the vertical components of the nearest bond vector emanating from the site, summing up to zero. Thus, only the upper and lower edge sites contribute to the summation for the vertical y-component in the first term of equation (28). This yields:

$$\begin{equation} \left[\widehat{\mathsf{H}}_{t},\mathbf{J}_{\text{tot},\sigma}\right] = \left(-2 i \right) \Bigg\langle t^2 \sum_{\mathbf{R}^{^{\prime}}} \, \sum_{ \lt \mathbf{R}, \mathbf{R}^{^{\prime}}>} \left(\mathbf{R}-\mathbf{R}^{^{\prime}}\right)_y \widehat{\mathsf{n}}_{\mathbf{R}^{^{\prime}} \sigma} \Bigg\rangle = 2 i t^2 a \sum_l \left( \langle \widehat{\mathsf{n}}_{l,1, \sigma} \rangle - \langle \widehat{\mathsf{n}}_{l,4, \sigma} \rangle \right). \end{equation} \tag{ 29 }$$

For edge states that are well-localized on the zigzag edges this result is exactly the structure anticipated from equation (8), namely only the boundary contributions remain. Now, we estimate the second term of equation (28), consisting of next-nearest bond contributions. It may be ignored. This is because the next-nearest bond terms are small when the distance between these bonds, $|\mathbf{R}-\mathbf{R}^{^{\prime}}|$, exceeds the correlation length ξ. Note that the one-particle correlation function, expressed as:

$$\begin{equation} C_{\sigma}\left(\mathbf{R},\mathbf{R}^{^{\prime}}\right) = \langle c^+_{\mathbf{R} \sigma} c_{\mathbf{R}^{^{\prime}} \sigma} \rangle \end{equation} \tag{ 30 }$$

remains significant only within the correlation length. As U increases, Δ also rises, leading to a decrease in the correlation length $\xi \sim \hbar v_F/\Delta$ (where v_F represents the velocity of Dirac electrons).

The expectation value of the vertical component of the gap parameter contribution equation (26) is found to be

$$\begin{align} \langle [ \widehat{\mathsf{H}}_{U,\text{HF}}, \mathbf{J}_{\textrm{tot},\sigma} ] \rangle \big \vert_y = & -i t a \, \sum_l \, \Big[ \frac{1}{2} \left(\Delta_{l,2,\sigma}-\Delta_{l,1,\sigma} \right) \langle \widehat{\mathsf{c}}^{\dagger}_{l,2,\sigma} \widehat{\mathsf{c}}_{l,1,\sigma} +\text{h.c.} \rangle \nonumber \\ &+ \left(\Delta_{l,3,\sigma}-\Delta_{l,2,\sigma} \right) \langle \widehat{\mathsf{c}}^{\dagger}_{l,3,\sigma} \widehat{\mathsf{c}}_{l,2,\sigma} +\text{h.c.} \rangle +\frac{1}{2} \left(\Delta_{l,4,\sigma}-\Delta_{l,3,\sigma} \right) \langle \widehat{\mathsf{c}}^{\dagger}_{l,4,\sigma} \widehat{\mathsf{c}}_{l,3,\sigma} +\text{h.c.} \rangle \nonumber \\ &+\frac{1}{2}\left(\Delta_{l+1,2,\sigma}-\Delta_{l,1,\sigma}\right) \langle \widehat{\mathsf{c}}^{\dagger}_{l+1,2,\sigma} \widehat{\mathsf{c}}_{l,1,\sigma} +\text{h.c.} \rangle \nonumber \\ &+ \frac{1}{2} \left(\Delta_{l,4,\sigma}-\Delta_{l+1,3,\sigma}\right) \langle \widehat{\mathsf{c}}^{\dagger}_{l,4,\sigma} \widehat{\mathsf{c}}_{l+1,3,\sigma} +\text{h.c.} \rangle \Big]. \end{align} \tag{ 31 }$$

The site numbers of the unit cell used in this equation are defined in figure 7. For a ribbon with a large number of unit-cells, the expectation values of the nearest-neighbor bond operators will be identical between unit cells. Next let us focus m = 3 site. In the presence of translational symmetry the expectation values of bond operator with m = 4 site in the lth and the l − 1 cells are identical. Then, the sum of the nearest-neighbor bond expectation values connected to m = 3 site in the lth unit cell takes the following form:

$$\begin{equation} \text{sum} = \Delta_{l,3,\sigma} \left( \langle \widehat{\mathsf{c}}^{\dagger}_{l,3,\sigma} \widehat{\mathsf{c}}_{l,2,\sigma} +\text{h.c.} \rangle - 2 \times \frac{1}{2} \langle \widehat{\mathsf{c}}^{\dagger}_{l,4,\sigma} \widehat{\mathsf{c}}_{l,3,\sigma} +\text{h.c.} \rangle \right). \end{equation} \tag{ 32 }$$

Our numerical evaluation shows that, for $U\gt\gt 2t$, the correlation function values satisfy:

$$\begin{equation} \langle \widehat{\mathsf{c}}^{\dagger}_{l,1,\sigma} \widehat{\mathsf{c}}_{l,2,\sigma} \rangle \approx \langle \widehat{\mathsf{c}}^{\dagger}_{l,2,\sigma} \widehat{\mathsf{c}}_{l,3,\sigma} \rangle \approx \langle \widehat{\mathsf{c}}^{\dagger}_{l,3,\sigma} \widehat{\mathsf{c}}_{l,4,\sigma} \rangle . \end{equation} \tag{ 33 }$$

We thus find that the difference between the two expectation values of bonds in equation (32) is small. This means that the contribution of the gap parameter $\Delta_{l,3,\sigma}$ will be negligible. The same argument can be applied to the m = 2 site, too.

Therefore, the result equation (31) simplifies to

$$\begin{equation} \langle \left[ \widehat{\mathsf{H}}_{U,\text{HF}}, \mathbf{J}_{\textrm{tot},\sigma} \right] \rangle \big \vert_y \approx - i t a \mathcal{B} \left( \sum_{l} \left( \Delta_{l,4,\sigma}-\Delta_{l,1,\sigma} \right)\right), \end{equation} \tag{ 34 }$$

where $\mathcal{B}$ denotes the expectation value of the nearest-neighbor bond operator (as given in equation (33)) plus its hermitian conjugate. The simplified result equation (34) tells us that the contributions from the upper and the lower zigzag edges are dominant. Now employing the temporal adiabatic condition equation (25),

$$\begin{equation} \Big\langle \left[\widehat{\mathsf{H}}_{t},\mathbf{J}_{\text{tot},\sigma}\right] \Big\rangle \big \vert_y = - \Big\langle \left[ \widehat{\mathsf{H}}_{U,\text{HF}}, \mathbf{J}_{\textrm{tot},\sigma} \right] \Big\rangle \Big \vert_y \end{equation} \tag{ 35 }$$

we arrive at

$$\begin{equation} \sum_l \left( \langle \widehat{\mathsf{n}}_{l,1, \sigma} \rangle - \langle \widehat{\mathsf{n}}_{l,4, \sigma} \rangle \right) \approx \frac{\mathcal{B}}{2t} \sum_{l} \left( \Delta_{l,4,\sigma}-\Delta_{l,1,\sigma} \right). \end{equation} \tag{ 36 }$$

We have numerically tested this result for various widths, as depicted in figure 8. As expected, the agreement is excellent for $U\gtrsim 3t$. When the terms dependent on-site gap values from the next-nearest bonds are added to the right side of equation (36), there is an almost exact agreement, even for small values of $U/t$. The above result explicitly demonstrates that the edge occupations of the zigzag nanoribbons are determined by the values of the gap parameters in the vicinity of the edge. This allows us to recognize the edge occupations of the zigzag nanoribbons as a topological charge linked to the solitonic behavior of the gap parameters. This identification parallels equation (8), which associates the induced fermion numbers with the topological charge of a soliton. When a topological charge is well-defined, i.e. when equation (36) is satisfied, it implies that the corresponding edge state is highly localized at the boundary sites.

Zoom In Zoom Out Reset image size

So far we have investigated the topological charges of zigzag ribbons. It turns out that an interacting rectangular armchair ribbon with two zigzag edges and two armchair edges may also support edge states on its two zigzag edges. As the ribbon width increases more topological edge states are formed in the gap, as shown in figure 9. When a staggered potential is present [33], mirror and time reversal symmetries are broken, as is their combined operation. The HF method shows that the charge of each zigzag edge may assume non-integer values, i.e. they are not quantized.

Zoom In Zoom Out Reset image size

We investigate whether ribbons in the low doping limit, characterized by the edge spin density ground state [29], possess topological edge states. Our numerical studies reveal that within the low-doped region, the ribbon can host nearly zero-energy gap states with a topological charge, as illustrated in figure 10. However, it is important to note that in the undoped case, the edge states have energies around $E\approx \pm\Delta/2$, as depicted in figure 2.

Zoom In Zoom Out Reset image size

Our findings suggest that chiral symmetry breaks due to on-site repulsion among electrons. However, an abrupt spatial change in the gap at the zigzag edges enables the emergence of topological charges with energies around $\pm\Delta/2$, leading to the formation of a symmetry-protected topological insulator. This symmetry is dependent on combined operations involving parity (mirror reflection) and time reversal transformations. The importance of this combined symmetry may be further illustrated by considering a staggered potential [33]. In the presence of such a potential, both mirror and time reversal symmetries are broken, as is their combined operation. The HF method shows that the charge of each zigzag edge may assume non-quantized charges as opposed to integer charges.

Are topological edge states robust in the presence of disorder? The introduction of disorder induces a topological phase transition from a symmetry-protected insulator to a topologically ordered insulator (see [23, 30]). Disorder in zigzag ribbons is a singular perturbation that gives rise to instantons. It also leads to the emergence of new edge states exhibiting e/2 fractional charges, which are protected by topological order, despite the presence of broken symmetries. But these new edge states are mixed chiral edge states, i.e. they are bonding and antibonding states of the chiral edge states. So, in this sense, the topological edge states are robust in the presence of disorder.

The boundary states (edge states) of our system are gapful. Therefore, it is not clear to what extent the classification scheme of references [24, 25] can be applied to our system. Moreover, the quasi-one-dimensionality of zigzag ribbons is important as it spatially separates the opposite zigzag edges of different chiralities. Disorder couples these edges and produce $e/2$ fractional charges. We believe that zigzag graphene nanoribbons are outside the conventional classification scheme. It would be interesting to expand the classification scheme to zigzag ribbons.

In low-doped zigzag ribbons, the ground state forms an edge spin density wave rather than a ferromagnetic edge. This state supports nearly zero-energy zigzag edge states, distinct from the $\pm\Delta/2$ energies, exhibiting a gap-induced topological charge.

We hope that our work will stimulate the study of atomically precise nanoribbons recently fabricated [6, 7] from the perspective of symmetry-protected topological insulators. It would be intriguing to observe topological charges of nearly zero energy states in low-doped ribbons that exhibit the edge spin-density wave ground state. Scanning tunneling microscopy [35] could be employed for this purpose.

All data that support the findings of this study are included within the article (and any supplementary files).

For a precise description let us introduce Dirac Hamiltonian. The the classical Lagrangian density of the gapless Dirac fermion in the presence U(1) gauge potential A_µ is given by (Einstein summation convention assumed)

$$\begin{equation} \mathcal{L} = \bar{\psi} i \gamma^\mu \left(\partial_\mu + i e A_\mu \right) \psi, \quad \bar{\psi} = \psi^{\dagger} \gamma^0, \end{equation} \tag{ 37 }$$

where the spacetime indices $\mu,\nu$ take the value $0,1$ for time t and space x, respectively. The gamma matrices are defined by

$$\begin{align} &\left\{\gamma^\mu,\gamma^\nu \right\} = 2 g^{\mu \nu} = 2 \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}, \quad \mu = 0,1 \end{align} \tag{ 38 }$$

$$\begin{align} \gamma^0 & = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix}, \quad \gamma^1 = \begin{pmatrix} 0 & i \\ i & 0 \end{pmatrix} = i \sigma_x, \quad \gamma^5 = \gamma^0 \gamma^1 = \begin{pmatrix} 1 & 0 \\ 0& -1 \end{pmatrix}. \end{align} \tag{ 39 }$$

The representation of gamma matrices where γ⁵ is diagonal is called the chiral representation, so equation (38) is a chiral representation.

Writing ($R,L$ denotes the handedness )

$$\begin{equation} \psi = \begin{pmatrix} \psi_R \\ \psi_L \end{pmatrix}, \end{equation} \tag{ 40 }$$

the Lagrange equations of motion in the absence of the gauge field yield

$$\begin{equation} \left(\partial_0 + \partial_1 \right) \psi_R = 0, \quad \left(\partial_0 - \partial_1 \right) \psi_L = 0, \end{equation} \tag{ 41 }$$

which implies that the right handed fermion is right moving, and vice versa.

The eigenvalue of γ⁵ is defined to be the chirality, and the following properties of γ⁵ play an important role:

$$\begin{equation} \left\{\gamma^5,\gamma^\mu \right\} = 0, \quad \gamma^\mu \gamma^5 = -\epsilon^{\mu \nu} \gamma_\nu, \quad \epsilon^{01} = +1 \end{equation} \tag{ 42 }$$

where ε^µν is a totally antisymmetric tensor and $\gamma_0 = \gamma^0, \gamma_1 = -\gamma^1$.