Excited spin density waves in zigzag graphene nanoribbons

In addition to the well-known anti-ferromagnetic and ferromagnetic edge states in zigzag graphene nanoribbons (GNR), we find that there also exist some excited spin density wave (ESDW) states, the energies of which are close to the anti-ferromagnetic state (ground state). We thus argue that these ESDW states may coexist in experiment. Our numerical results from the self-consistent mean-field method as well as the first-principles calculations indicate that the allowed ESDWs are commensurate; and their dispersion curves are linear in the long wave limit. The coupling of the two edge portions in an ESDW becomes very weak in a wide GNR, in which case these ESDW portions may have different phases. Finally our calculations in the non-equilibrium Green’s function theory combined with the Hubbard model show that these ESDW states can also exist in some localized (middle or terminal) region of GNR.


Introduction
Graphene, an important two-dimensional material, has attracted a lot of research interests in the past decade [1][2][3]. It has many novel properties such as Dirac-like electron behavior with the Klein tunneling effect [2,4]; the high opacity in the optical wavelength range [5]; the edge states on the zigzag graphene nanoribbons (GNR) [6]. For the GNRs, they have different energy gaps due to the widths and axis orientations and they may be metallic, semiconducting or insulated [7].
Among them, the magnetic properties of graphene are very interesting. It is well known that there exist spontaneous spin polarizations on the edges of zigzag GNR [8]. The two edge regions have either antiferromagnetic (AF) or ferromagnetic (FM) states: the later may be induced by the external magnetic field. The FM state often possesses a little higher energy than the AF state. Their energy differences decay as the increase of the ribbon width [8,9]. Based on the first-principles calculations or the tight-binding (TB) method with Hubbard model, many spintronic graphene structures are proposed. Some examples are the anti-dot lattices [10], the saw-tooth GNRs [11], the porous graphene sheets doped with nitrogen and boron atoms [12] and the dumb-bell-shaped GNR [13]. Their energy bands often exhibit some half-metallic property with potential applications in the spintronic devices. Besides the periodic systems, magnetic orders are also found in finite structures such as graphene nanoislands [14]. These ground states are in agreement with the Lieb's theorem [15].
Recently Magda et al have observed a semiconductor-to-metal transition in the GNR measurements as the increase of ribbon widths at the room temperature [16]. This results from the magnetic order change from the AF state to the FM state by the shift of the chemical potential.
However, most of these magnetic orders, existing on the zigzag edges, are either in the same direction (FM order) or the opposite direction (AF order) [17]. A very few works show more complicated magnetic states on one zigzag edge. Kim and Kim's calculations show that on the graphene-based spin valve devices, the zigzag GNR's magnetism may have opposite FM orders with a single domain wall (DM) in the middle region due to different magnetic fields in the leads [18]. With the first-principles calculations and non-equilibrium Green's function (NEGF) theory, the transmission currents under some bias voltages are obtained. Furthermore, the spin waves in zigzag GNR have been studied with the spin transverse dynamic susceptibility theory [19][20][21][22].
These studies are limited to simple and model systems. Recently Hagymasi et al use the density-matrix renormalization-group algorithm to study some excitations in zigzag GNR [23]. We also notice that some spin density wave (SDW) has been proposed in a zigzag GNR, which is induced by a periodically shear strain on the GNR [24].
In this study, based on the TB with Hubbard model and the first-principles calculations, we show that there exist some intrinsic spin variations in a uniform (unstrained) zigzag GNR. These states, corresponding to some excited spin-density-waves (ESDW), have very similar band gaps. We emphasize that this ESDW is not the common spin wave in traditional magnetic materials, which is described by the localized spin model, such as the Heisenberg model. They are the excitations of the AF/FM state, and are described by the Hubbard model with the kinetic energy and on-site Coulomb interaction. These states consist of various FM/AF orders or their combinations on two edges of a GNR supercell. In narrow GNRs, the magnetic orders on two edges are found to have strong coupling. Combined with the NEGF theory, we further find that in the local regions of a zigzag GNR there also exist similar SDW excitations. Our dispersion calculations show that the energies of these ESDW states may be arbitrarily close to the AF ground state, so they can partly contribute to the magnetizations of the insulated zigzag GNR in experimental measurements.

Hubbard model and band structure calculation
In the unrestricted Hartree-Fock model with the mean-field approximation [25][26][27], the Hamiltonian is written as The first two terms are the hopping energy of the TB model with two spin components. 〈i,j〉 means the hopping terms exist only between the nearest neighboring sites and t=−2.7 eV. The last term is also called the Hubbard term, resulting from the Coulomb repulsion of two on-site electrons with different spins. Here we choose U=2.0 eV for this value which can best reproduce the first-principles calculation results for zigzag GNR [26,27]. á ñ  n i and á ñ  n i are the average electron densities at site i with different spins. In a periodic system the band structure is obtained from the eigenvalue calculation [28], where N 0 is the number of atoms in a unit cell, k is the Bloch wavevector, n=−1, 0, 1, corresponds to three neighboring unit cells (in 1D case). d i j , =1 (i=j) or 0 (i≠j), comes from the orthogonal relation of two bases functions in the TB model. Equation (2) may be rewritten in a matrix form where  ( ) H k is from the combination of the Hamiltonians (in equation (1)) in the neighboring unit cells with the Bloch factor ⋅ e . k R i n When considering the Hubbard term, we may double the dimension of the Hamiltonian and the eigen-basis with two spin components. Instead of doubling the dimension, here we use two sets of eigenvalue equations corresponding to two spin cases. From equation (1) the diagonal elements of the each Hamiltonian consists of the electron densities with opposite spin, which makes the two equations coupled to each other. This is a selfconsistent (SC) problem and has to be solved with the numerical iterations. In equation (1), á ñ   ( ) n i is the average electron density on each site, which comes from all eigenstates in the occupied orbitals and averaged in the first Brillouin zone, as shown below [25] 'occ' means the summation is utilized only for the occupied orbitals (for zero temperature); p R is the volume of the first Brillouin zone (1D case).   ( ) is the eigenvector of ith site for the pth eigenvalue with the Bloch wavevector K. In our work all the bands (with spins) are half-filled expect the FM state.

NEGF theory
Besides the periodic GNR, we also investigate the non-periodic GNR with the NEGF theory [29]. The Hubbard term exists both in the device and the leads, where I is the unitary matrix with the dimension of the device; S S S = + r L r R r is the retarded self-energy, which is given by the coupling Hamiltonians and the surface Green's function a g r (see equation (7), a = L R , , which means the left or right lead.) a g r is obtained from the recursive Green's function method [30]. The imaginary part of ) f E , is the Fermi function with the chemical potential m. Similar to the band calculations before, we use two sets of matrix equations from equation (6). They are coupled to each other by the electron density   ( ) n i with opposite spins. This is because that in equations (6) and (7)   ( ) n i appears in the Hamiltonians and self-energies, and they also result from   ( ) G r, and the integral (equation (9)). So this is also a self-consistent problem and has to be solved iteratively. If we want to calculate the semi-infinite system, we may set the self-energy of one lead (such as right lead) to be zero, which means S S = .
r L r Then the calculated   ( ) G r, and   ( ) n i correspond to a semi-infinite system. The integral in equation (9) can be effective evaluated by the residue theorem [31,32]. The details can be seen in appendix A and [33,34].

Results and discussions
3.1. Spin polarizations of the FM-type and AF-type ESDW In this part we use the band structure calculation to study the ESDW states in the GNR supercells. The TB model with the Hubbard term and the density-functional-based-tight-binding (DFTB) model [35] are employed in the calculations. We also calculate the Hubbard energies and total energies of these ESDW for their dispersion relation. Some hybrid ESDW states in long GNR supercells are also investigated.
In the paper of Kim and Kim [17], a GNR system with two FM orders is proposed, which is induced by two opposite magnetic fields in the leads. These two FM regions are flipped to each other through a single domain wall in the middle.
Here in our SC calculations we further find that there exists such 'FM flip state' in a periodic GNR system. To study the complicated states in the periodic system, a supercell of GNR is used, as shown in figure 1. With the proper initial condition of á ñ   ( ) n i in the SC calculation as described in section 2, we find there exists some stable state with the flipped FM orders in the periodic GNR regions. Figure 2 for this state in a GNR supercell (M=12, N=8). In this state, the GNR edge on the left (or right) side has the same FM order, and this order flips into the opposite direction on the other side. From the 3D plot we see that the FM order changes gradually like a sinusoidal function. In the middle part the magnetic order is very small, which can be regarded as a domain wall. So this state is a type of ESDW state. With the similar approach, we also find another SDW excitation in the same GNR supercell as shown in figure 2(b). This state has the AF order: opposite spin polarizations exist on the two edges of GNR, and these AF orders gradually change and flip in another period as an ESDW. Now we use the first-principles calculations with the DFTB model to confirm these novel spin polarization states. This model is an approximation of density-functional-theory (DFT) method derived from the secondorder expansion of Kohn-Sham energy around the reference charge density [35][36][37]. In this model, the KS matrix elements of the DFT Hamiltonian and the overlap integrals of non-orthogonal orbitals are numerically computed and stored in Slater-Koster [38] tables as function of the inter-element distance. A free software package: DFTB+ [35][36][37] developed by the Th. Frauenheim's group is used in our paper. A supercell is used in the band calculations and hydrogen atoms are added on the edge carbon atoms. The vacuum layers of 12 nm in the ribbon plane and 14 nm in the normal direction is applied to avoid the interaction between GNR in different unit cells. In the SC iteration we choose the Broyden mixing scheme and the energy tolerance is set to 5 × 10 −6 eV. The collinear spin-polarization algorithm with the PBE exchange-correlation functional is chosen in our calculations. In the band calculation 100 k-points are used in the first Brillouin zone and the temperature is set as 300 K.
To obtain these ESDW states, we initially set some spin configuration which is close to the final spin polarization. Then we do the SC calculation by this DFTB code. It is noted that in common cases DFT calculates for the ground state. However DFT can also deal with some metastable states of GNR such as the FM state [9] and the domain wall configuration [18], which has the energies slightly higher than the ground state. We will investigate these metastable states and their energies in the following sections. Figure 3(a) shows the spin polarizations for the FM-type ESDW state in a GNR supercell (M=12; N=10) from the DFTB calculation. To obtain this state, a proper initial spin configuration is set in the input file of DFTB+package. We see that the DFTB+method can generate a similar state, compared to that from the TB +Hubbard method (figure 2(a)). The only difference is that there is some phase shift for the domain wall. But since the system is periodic, it is also reasonable. We also obtain the AF-type ESDW state by the DFTB+software, as shown in figure 3(b), which is similar to the state obtained from our TB+Hubbard model ( figure 2(b)).

Band structures of ESDW states
We also calculate the band structures of these AF/FM-type ESDW states from the DFTB model and our TB+Hubbard method, as shown in figures 4(a) and (b). In these bands two spin components are degenerate. We see that the bands from the DFTB calculation are similar to those from the TB+Hubbard calculation. From the insets of figures 4(a) and (b), we see the bands also have the similar properties by these two methods: above the Fermi level the first two bands of the FM-type ESDW are separated; and these bands of the AF-type ESDW are connected at the Brillouin zone boundary (k * R = π). These properties also remain the same in figure 4(c). We will interpret the physical origin in detail later. There still some difference of shapes in these bands from DFTB and TB+Hubbard methods. For example, the two DFTB bands above the Fermi level have a larger separation than the TB bands. This comes from the different methods of the two models.
In figures 4(a) and (b), both AF-type and FM-type ESDW states have a band gap near the Fermi level. This gap comes from the Coulomb repulsion on the zigzag edges as in the pure AF state. We here explain this with  more details. In a simple zigzag GNR there exist a very flat band near the Brillouin boundary even without the Hubbard term [7]. The electrons in the flat band have nearly zero group velocity and their kinetic energy is almost quenched. Thus the Coulomb interaction plays a major role to form some many-body states, such as FM or AF state [39]. Furthermore, these flat bands results in a very large density of states (DOS). We also check that for the zigzag edge atoms in the GNR. Their LDOS near the Fermi level also approach to infinity. So an infinitesimal small on-site repulsion (with the opposite spin on one edge atom) will lead a spin separation either in the vertical direction (which forms the AF state) or in the horizontal direction periodically (which forms the ESDW states).
We further investigate the bands for two ESDW states and the pure AF state in a GNR supercell (figure 4(c)). The band of the AF state comes from the simple Brillouin zone folding of a single GNR system [9]. For the FMand AF-type ESDW states, one or two gaps are opened at the boundary or center of the Brillouin zone. Similar to the charge density waves, the two degenerate states in the Fermi surface (here are the Fermi points in 1D case) are related by a nesting vector. They have interaction with each other to lead to a gap which lowers the total energy of the system. The nesting vector is formed by the new periodic potential field, which results from the spin condensation in these ESDW states [40]. We find that all these ESDW state are commensurate with the GNR lattice.
In figure 4(c) we notice that besides some gaps at k * R=0 or k * R =π, for the AF-type ESDW the gap between the first and second bands (above the Fermi level) at k * R=π disappears. And the gap between the third and fourth bands at k * R=π also disappears. This seems not be consistent with the previous statement about the gap opening mechanism due to the spin condensation. However, considering the spatial distribution of spin polarizations (figure 2), we find positive and negative spin polarizations on the two edge atoms. The Hubbard terms on the edges have different signs for the AF-type and FM-type ESDWs. These diagonal terms provide additional potentials to affect the gaps (see figure 4(d)). Appendix B gives the detailed demonstration.
We also calculate some other GNR supercell systems. The results are shown in table 1. We see that the gap of the FM-type ESDW state is usually slightly smaller than that of the AF-type ESDW state. In a wider GNR, the gap difference between two ESDW states is smaller than that in a narrow GNR.

Energies of ESDW states
In the TB calculations we find that the involution of the Hubbard term would also result in the change of the kinetic energy. Here we emphasize the energy related to the Hubbard term, called the Hubbard energy ( ) E . Hubbard It is from the following Hubbard Hamiltonian (equation (10)), and often calculated by the mean-field approximation (equation (11) The calculated E Hubbard is then averaged in the first Brillouin zone with the Brillouin sampling number N k =400. The large sampling is necessary for a converged result in the SC calculation. Here the energy is scaled by the number of carbon atoms in the GNR supercell. So it reflects the intensive spin excitation energy and it is not proportional to the supercell size. Figures 5(a) and (b) show the Hubbard energies and total energies (sum of the eigen-energies of the occupied states) of the AF-type ESDW states in different lengths of supercells (with the width parameter M=12). The horizontal axis is the wavevector q (q corresponds to the supercell period: = p q , R 2 where = R N a 3 , R is the supercell length and it is also the wavelength of ESDW). They are the dispersion curves of these ESDW modes. The curve is quite linear near the Brillouin zone center (with very large supercell periods) and then it deviates at large wavevectors. This curve also agrees well with the work in literatures [18,20]. We find that in some short supercells (for example, M=12, N=7 or less), the AF-type or FM-type ESDW states does not exist in our SC calculations. In other words, they are unstable modes. In the paper of Culchac et al they noticed that some highenergy spin excitation modes in GNR are strongly damped and impossible to observe [19]. Our calculations give similar results. Here we notice that the ESDW states on the two edges are not the standard harmonic (sinusoidal) modes. We will see this more apparently in very long or wide supercells later (see figures 6 and 7). For these itinerant electron systems from the Hubbard model, Tian has proved that the collective spin excitations are gapless from the AF ground state, based on the Lieb's theorem [41]. It is consistent with our linear dispersion ESDW curve. This means that there exist the very long-wavelength ESDW state with an energy arbitrarily close to that of the AF ground state. So we believe these ESDW can coexist with the AF state in GNR at any temperature. Figure 5(c) shows the maximum magnetization and Hubbard energies of AF-type ESDW states in different widths (M) of supercells (with the same length N=12). The magnetization of the AF states is also plotted for comparison. For the AF states the magnetization increases monotonously [9]. For the ESDW states the magnetization increases with small widths (M <= 14) and then oscillates with large widths (M>14). We check the edge spin distributions and find that when M is larger than 14, the eigenmode is not sinusoidal, which consists of high harmonic modes as those in figure 6. Thus their magnetization behavior differs from the sinusoidal AF-type ESDW modes with monotonous increasing magnetization on edges. The magnetization of these hybrid ESDW modes are often much higher than those of the sinusoidal AF-type ESDW modes with a smaller width. Figure 5(d) shows an energy evolution from the AF state (ground state) to the AF-type ESDW state. As we stated before, either the AF/FM states or the ESDW states have the local minimum of energy in the spin configuration space. Here we give a calculation for this in a GNR supercell system. The detailed mathematical formulas are in the appendix C. We change the eigenvectors from the AF state to the AF-type ESDW state, with all the eigenvectors transformed linearly from one state to another, except one component left to adjust for maintaining the normalization of these eigenvectors. The parameter i x reflects this process (the points i x =0 and 50 correspond to the two metastable states). It is easy to see that the two local minimums correspond to the AF and AF-type ESDW states. In others words, we give a numerical verification for the ESDW to be a metastable state in the spin configuration space. Table 1 shows the band gaps, Hubbard energies and total energies in different ESDW states. We see that the pure AF state has a lower energy, compared to the pure FM state. This is aggregable with others' work [8,16]. In table 1 we see that the AF-type ESDW often has a lower energy than the FM-type ESDW. The gaps in these ESDW states are much smaller than that of the pure AF state (about 0.36 eV). When the GNR supercell becomes wider (M changes from 12 to 16), the energy difference between the two ESDW states becomes smaller (with the same supercell length: N=8). We also observe the energy gaps of the two ESDW states tend to the same value in a very wide supercell. This means in a wider GNR system the SDW excitations on the two edges have less interaction with each other. In the limit of a very wide GNR, they behave as two independent modes.

Hybrid ESDW in long GNR supercells
With these self-consistent calculations, we also discover a lot of hybrid ESDW states in the GNR supercells. For example, beginning with different initial spin configurations, we obtain the hybrid ESDW state 1 (figure 6(a)) and state 2( figure 6(b)). We see that in the hybrid state 1, the left part of upper edge has an oscillation of a small period (N/2); and the right part of the upper edge has an oscillation of a large period (N). The oscillation in the lower edge is similar but the order is reversed. In the hybrid state 2, the upper edge has the large-period oscillation and the lower edge has the small-period oscillation. Figures 6(c) and (d) show the magnetization on the two edges for the above hybrid ESDW states (upper parts). Then we employ a discrete Fourier transformation (DFT) for a harmonic mode analysis, as shown in the same figures (lower parts). We see that in the hybrid ESDW state 1 there are about 4 modes which contribute to this spin order. In the hybrid ESDW state 2, the harmonic mode 1 (with the period of N) and mode 3 (with the period of N/3) are two main contributions to the lower-edge spin order; and the harmonic mode 2 (with the period of N/2) is mainly attributed to the upper-edge spin order. Figure 6(e) shows the third ESDW state: the upper edge has the pure FM order and the lower edge has a spin oscillation with a period of the half GNR length. From this state we can conclude that the two spin orders are almost independent on the upper and lower edges.

Out-of-phase ESDW in wide GNR supercells
For the wide GNR systems, we find that the ESDW portions on the two zigzag edges may have different phases, as shown in figure 7. We obtained these ESDW states by setting the proper initial spin configurations in the selfconsistent calculations. Figure 7(a) shows the TB result for the AF-type ESDW state with the same phases on two edges in a wider GNR (M=32, N=10). The DFTB calculation also gives similar in-phase results. Figure 7(b) shows the TB result for the ESDW state with different phases on two edges. In calculations we find for some wide GNR, such as this case (M=32, N=10), the mean-field Hubbard model does not give a stable iteration result for the out-of-phase or in-phase ESDW. We choose the Hubbard model without the mean-field approximation (equation (10)) to reach a convergent SC result. Figure 7(c) shows the similar DFTB result for the ESDW state with different phases on two edges. However, both by the TB and the DFTB models we find that in narrow GNR systems, even if the initial spin configuration is set to be out of phase, the final self-consistent results always give an in-phase ESDW portions on two edges return. This means in the narrow GNR, the two edge portions of ESDW have strong coupling which coordinates their phases.
The calculated Hubbard energies from TB method or total energies from the DFTB method are listed in table 2. The data show that the energies for the out-of-phase ESDW state and the in-phase ESDW state are almost the same (degenerate). From these results we conclude that in a wide GNR the ESDW portions on two edges have a much weaker coupling with each other, so they may have different phases.

Localized ESDW in non-periodic systems
Besides the periodic GNR systems, we also use the NEGF theory coupled with the TB+Hubbard model to calculate some non-periodic GNR systems. We find that these interesting spin polarization states also may exist in some local regions of a uniform GNR. The system consists of a device and two leads, as shown in figure 1. The calculation process is given in section 2.
Firstly a benchmark NEGF calculation is done for a uniform GNR. We calculate the transmission spectrum and band structure for a FM state in such GNR system ( figure 8(a), in the transmission the device size is: M=12, N=8 and the two leads sizes are M=12, N=2). The transmission spectrum agrees well with the FM band structure (also see [16]). In some energy regions the transmissions are different for two spins in a FM state. The integer transmission plateaus indicate this system is a uniform periodic system.
We then calculate the Green's function and the DOS of the device by using equations (6)- (8) and do the energy integral with equation (9) to obtain the electron density. The contour integral method is used in the integration (see appendix A). After obtaining the electron density, we put them into the Hubbard term and calculate the Green's function again for the iterations until a SC solution is reached. In the beginning of the SC calculation, we set some initial electron densities with spins in the Hubbard term. Different initial spin configurations, for example, the FM-type or AF-type ESDW in the middle of the device are set in the beginning of the SC calculation. To calculate the self-energy, electron densities with spin configurations (such as AF order and FM order) are set in the Hubbard terms of the lead Hamiltonians similarly.
Two SC solutions are shown in the following figure. Figure 8(b) is the FM-type SDW excited in the middle of the device, while the leads have the opposite FM orders. Figure 8(c) is the AF-type SDW excited in the middle of the device while the two leads have the opposite AF orders. We also observe that the transmission spectra of these two states are almost the same and without any spin separation.
We also find that with the initial sinusoidal spin configuration, the FM-domain-wall state is reached after the SC calculation (left panel of figure 8(d)). This is similar to the results in the literatures [18,42]. Furthermore, we also obtain some hybrid ESDW state within this SC NEGF method, as shown in the right panel of figure 8(d).
To evaluate the total energy of these local ESDW states in the device, we use the following formula   energy gap of GNR (AF state has a much larger gap). This will influence the absorption spectrum, as well as the electronic behaviors when it is used in the graphene transistor devices. However, since the energy of ESDW is higher than AF state (ground state), the occupation probability of ESDW is usually smaller. Finally we show the excited SDW states in terminal of semi-infinite GNR systems. The calculation process is very similar with that of the lead-device-lead system above, except one lead is not coupled to the device.
Our results are shown in the following. Figures 9(a) and (b) are two states with a half period of ESDW (AF or FM order) in the terminal region, and the magnetic orders in the inner region tends to the same (AF or FM order). The polarization near the terminal side is observed to be larger. These states can be explained as the surface ESDW excitations in semi-infinite ribbons from the original AF or FM order.
Similarly, with some other initial spin configurations, we find out another two hybrid ESDW surface states, as shown in figures 9(c) and (d). These two states begin with the FM (or AF) order near the terminal side and then transforms into the AF (or FM) order in the inner region.

Conclusions
We have systematically investigated the excited SDW modes in the zigzag GNR. Many interesting spin polarization states have been found, such as the AF-type and FM-type ESDW states, and some hybrid ESDW states. The first-principles calculation has been employed to confirm these states. We have calculated the Hubbard energy and the band structures of these states. All these novel states can be viewed as some spin excitations in the configuration space of magnetic orders.
The periods of these ESDW states are commensurate. In the Brillouin zone boundary or center the bands open some gaps due to the interaction of two degenerate states related by nesting vectors. Two simple models have been used to explain the disappeared gaps at the Brillouin boundary for the AF-type and FM-type ESDWs. The dispersion curve of the AF-type ESDW states is linear in the long wavelength limit. In a wide GNR the excited spin polarizations on two edges tend to be weakly coupled. The energy difference between the AF-type and FM-type ESDW states becomes small as the width of GNR is increased. We also have utilized the Lagrange multiplier method and numerical calculations to analyze the metastable states for these ESDWs. The edge  portions of ESDW in a wide GNR may possess different phases, while in a narrow GNR only the same phase exists. In very wide or long GNR supercells the ESDW modes are not sinusoidal and they have some high-order harmonic components. The amplitude of these ESDW states ranges from about 0.2 m B to 0.3 m .

B
Finally we have utilized the NEGF theory with the Hubbard model to calculate the SDW excitations in nonperiodic zigzag GNRs. Some interesting magnetic states such as the middle-part and terminal-part SDW excitation states have been investigated. The energy and transmission gap of these localized ESDW states have also been calculated by this method.
We believe these ESDW states are intrinsic spin excitations in a uniform zigzag GNR without any lattice deformation. They can coexist with the AF ground state in GNR at any temperature. These ESDW may be excited by some circular polarized light with angular momentum and detected by some technique such as neutron scattering. The coexistence of these states may influence the absorption spectrum, as well as the electronic behaviors of GNR, which will be studied in the future. Appendix B. Simple chain models for the band gap openning mechanism in ESDW systems To demonstrate the band gap behaviors at the Brillouin zone boundary in figure 4(c), we propose some atomchain models with one orbital for each atom. Firstly we adopt a double-chain model with 4 atoms in a unit cell (see the first two insets of figure 4(d)). According to the statements in the context, the Hubbard terms from the ESDW correspond to the positive or negative diagonal contributions in the Hamiltonian. So here we add two perturbation potentials (d 1 and d 3 on atom 1 and 3) on the two edges to simulate the different spin polarization in the AF/FM ESDW states. With the TB model, the eigenvalue calculation for the band structure leads to the following secular equation  1 3 the four eigen-energies are easily obtained as , , 2 For the FM-type ESDW, we set d d d = = , 1 3 the four eigen-energies are , , .So we see the gap for the AFtype ESDW is very small: are applied to simulate the Hubbard potentials on one edge. The calculated result (the right panel of figure 4(d)) clearly shows that due to the positive and negative potentials on the two segments of one edge, the third and fourth bands meet with each other at k * R=π without any gap.