Exact diagonalization study of the tunable edge magnetism in graphene
Abstract
The tunable magnetism at graphene edges with lengths of up to 48 unit cells is analyzed by an exact diagonalization technique. For this we use a generalized interacting onedimensional model which can be tuned continuously from a limit describing graphene zigzag edge states with a ferromagnetic phase, to a limit equivalent to a Hubbard chain, which does not allow ferromagnetism. This analysis sheds light onto the question why the edge states have a ferromagnetic ground state, while a usual onedimensional metal does not. Essentially we find that there are two important features of edge states: (a) umklapp processes are completely forbidden for edge states; this allows a spinpolarized ground state. (b) the strong momentum dependence of the effective interaction vertex for edge states gives rise to a regime of partial spinpolarization and a second order phase transition between a standard paramagnetic Luttinger liquid and ferromagnetic Luttinger liquid.
pacs:
73.21.b,75.70.Rf,81.05.ue,73.22.PrI Introduction
Since it has first been isolated in the laboratory,(1) graphene, a twodimensional honeycomb lattice of carbon atoms,(2) attracts much attention. In fact, graphene has multiple amazing properties. To name only a few of them, it ranges among the mechanically strongest materials,(3) it shows a quantum Hall effect at room temperature,(4) and, due to its unusual Dirac band structure, it allows the study of relativistic quantum physics in a solid state environment.(5) Furthermore, its potential application as the basis of the next generation of electronic devices stimulated great efforts to gain experimental control as well as theoretical understanding of this astonishing material.
Usually a strong electron confinement increases the strength of electronelectron interactions by pushing the electrons close together. However, in spite of the extreme electron confinement to only one single layer of atoms, many experiments in graphene may be explained by assuming the electrons to be noninteracting. This is especially true for experiments probing the bulk properties of graphene, as the bulk density of states vanishes at the Fermi level, suppressing the manifestation of interaction effects. On the other hand, the properties of zigzag edges differ greatly from the bulk properties. Socalled edge states, i.e. onedimensional states with very small bandwidth, localized at these edges, give rise to a peak in the local density of states at the Fermi energy.(6) The enhanced density of states allows the electronelectron interaction to drive the zigzag edges to a ferromagnetic state with a magnetic moment localized at the edge.(7); (8); (9); (10); (11); (12) This phenomenon is known as edge magnetism.
At normal graphene edges the electronelectron interaction is so strong and the bandwidth of the edge states is so small that the spins of all electrons in the edge states are completely aligned. However, as has been proposed recently, graphene/graphane interfaces provide means to tune the bandwidth of the edge states to regimes in which the edge starts to depolarize and the edge magnetism is gradually suppressed until, for a critical edge state bandwidth, the magnetism disappears.(13)
In Ref. (13) it was argued that this interactioninduced magnetism can be understood on the basis of an effective model, describing the interacting onedimensional edge states only, while the bulk states are neglected. What at first glance appears to be a contradiction to the LiebMattis theorem,(14) stating that the ground state of interacting electrons in one dimension cannot be spin polarized, can be resolved by noting that the effective edge state model does not fulfill the prerequisites of the LiebMattis theorem.(13) The deeper reasons for the existence of a ferromagnetic ground state in a onedimensional interacting electron system, however, remained elusive. The present work is devoted to this issue.
In this paper we present a systematic exact diagonalization analysis of interacting edge states. Two striking features of edge states turn out to be most important for their magnetic properties: (a) the edge states exist only in a restricted part of the Brillouin zone and (b) the transverse edge state wave function has a strong characteristic momentum dependence. These features have consequences for the effective lowenergy theory, namely (a) no umklapp processes are allowed in the interaction Hamiltonian and (b) the interaction vertex acquires an unusually strong momentum dependence. In order to be able to study the consequences of these two features, we introduce a generalized model in which we add an artificial interaction term describing umklapp processes and allow the momentumdependence of the interaction vertex to be tuned from a momentumindependent vertex, as in usual metals, to the full momentumdependence, as it is found in edge states. Therefore, the generalized model can be tuned continuously from the limit in which it describes edge states to a limit which corresponds to usual onedimensional metals such as the Hubbard chain. We solve this generalized model for graphene zigzag edges of finite length unit cells (i.e. 12 nm) by exact diagonalization using the Lanczos method for the determination of the ground state of the effective model.(15); (16); (17); (18)
The paper is organized as follows. In Sec. II we review the direct model as it has been derived in Ref. (13) and introduce the more versatile generalized model with additional tunability. In Sec. III the exact diagonalization analysis of the generalized model is presented. Finally, the results are discussed in Sec. IV.
Ii Edge state models
In this section we introduce the models on which our analysis is based. The edge state model obtained from the direct projection of the honeycomb lattice Hubbard model to the Fock space spanned only by the edge states has been discussed in Ref. (13). This model will be called the direct model in the following. We identify two important features of the direct model: (a) the restriction of the Brillouin zone for the edge states and (b) the strong dependence of the transverse localization length on the momentum along the edge. After having analyzed the consequences of these features for the effective interaction vertex, we propose a generalized edge state model in which these features can be tuned. This allows us to investigate the impact of each of these edge state features on the magnetic properties. In particular, the generalized model can be tuned continuously from a Hubbard chain limit, i.e., a usual onedimensional metal without any ferromagnetic ground state, to the edge state limit with its ferromagnetic ground state.
ii.1 Direct derivation from the honeycomb model
We start from the simplest possible noninteracting tightbinding model of electrons in graphene
zigzag ribbons, taking into account only nearest neighbor hoppings of electrons , where runs over
nearest neighbor sites of a halfinfinite honeycomb lattice, is a collective site
index for the th unit cell and the sublattice (see Fig. 1), and
annihilates an electron at site with spin . Since we are exclusively
interested in the zero energy eigenstates of , the actual energy scale of
is unimportant so that we may drop it.
(1) 
where , is the momentum in direction (along the edge), and , with the number of unit cells in direction, i.e., along the edge. The dependent normalization constant can be interpreted as the weight of the edge state wave function right at the edge atoms where . It is easily seen that . As the edge state wave function is only nonzero on the sublattice we omit the sublattice index, setting it to .
The two most important features of the edge state wave function [Eq. (1)] are : (a) The edge state only exists for momenta . In the rest of the Brillouin zone the edge state wave function is not normalizable, as for these momenta. (b) In direction the edge state is sharply localized at the edge for , wereas for close to one of the Dirac points and , the wave function delocalizes into the bulk [see also Fig. 1(b)]. These two edge state properties are stable against adding more details, such as secondnearest neighbor hopping or various edge passivations, to the honeycomb Hamiltonian .(13); (20) The detailed analysis presented in this paper will clarify that the existence of edge magnetism and in particular its tunability are consequences of these two edge state properties.
The dependence of the localization length of has consequences for the edge states’ selfenergy as well as for their interaction vertex function . Neglecting the bulk state contributions,(13) the selfenergy correction due to a perturbation , which is invariant along the edge, is given by . Due to the delocalization of for near , edgelocalized perturbations lead to selfenergy corrections for which while . For sufficiently well behaved perturbations the selfenergy correction gives rise to a smooth edge state energy dispersion with a bandwidth . For a large class of these edgelocalized perturbations, the selfenergy correction approximately has the form
(2) 
with the dependent weight of the edge state wave function right at the edge. Eq. (2) expresses that an edge state which is more localized at the edge experiences a stronger selfenergy correction from an edgelocalized perturbation than an edge state which is delocalized into the bulk region. Examples of such perturbations are edge passivations, graphane termination,(20) or local interactions with the substrate. Note that the edge state bandwidth is experimentally tunable in various ways so that we consider as a free parameter. Therefore, the noninteracting part of the direct model (dm) edge state Hamiltonian is given by:
(3) 
where the sum is restricted such that only edge state operators with appear.
The effective interaction of the edge states, derived by projecting the Hubbard Hamiltonian on the twodimensional honeycomb lattice to the Fock space spanned by the edge states, reads(13)
(4) 
Together, we have the effective Hamiltonian . The interaction vertex is given by the overlap of the wave functions of all four fermions (with momenta ) participating in the interaction
(5) 
While the denominator, resulting from the geometric series over , turns out to lead only to unimportant quantitative corrections, the numerator of , which is the product of the wave function weights at the edge for each of the four fermion operators, leads to the momentumdependence of the interaction strength which is important for the stability of the weak edge magnetism. Essentially, the effective interaction becomes stronger the more localized the participating fermions are, i.e., the closer their momenta are to . If one or more of the momenta are close to the Dirac points , where the edge state wave functions delocalize into the bulk, the effective interaction is suppressed (see Fig. 1). Note that setting the denominator in Eq. (5) to unity corresponds to assuming that the Hubbard interaction is only present at the outermost line of carbon atoms right at the edge. Such an approximation has been used in Ref. (21). We find that this approximation is inessential for the edge magnetism, leading only to quantitative corrections.
An important consequence of the restriction of the summation in Eq. (4) is the absence of umklapp processes. As explained above, edge states only exist in one third of the Brillouin zone, i.e. for , so that no four fermion process with momentum exists. Indeed, within the restricted Brillouin zone, the process with the largest possible total momentum is , i.e. . Processes with larger total momentum leave the restricted Brillouin zone and are therefore suppressed, as they involve the overlap of edge states and bulk states, which is small. Also, most of the bulk states live in a different energy regime than the edge states.(13)
Thus, we have identified two properties of edge states which make them fundamentally different from usual onedimensional conductors:

Due to the restricted Brillouin zone, umklapp processes are forbidden.

The transverse localization of the edge state wave function gives rise to a tunable band width . Furthermore, the interaction vertex becomes weaker if the momenta of the participating fermions approach a Dirac point, i.e. .
We will show that these properties are the basis for the magnetism at graphene edges.
ii.2 Generalized model
We now introduce a generalized model in which the different aspects of the effective electronelectron interaction, found in the previous subsection, may be tuned independently. For this, we map the edge state operators which correspond to right (left) moving modes for (), to fermionic operators in which specifies the direction of motion and , i.e. and . The direction of motion corresponds to when used in formulas. Note that the zero point of has been shifted so that corresponds to for right and left movers, respectively (see Fig. 2).
For the noninteracting part of the generalized edge state model we assume a linear spectrum with slope
(6) 
This linearization of the selfenergy [Eq. (2)] only leads to inessential quantitative corrections (see also Appendix B).
The partitioning into left and rightmovers always gives rise to four terms in the interaction part of the Hamiltonian involving different combinations of left and rightmoving modes. Conventionally, these terms are called processes (see Fig. 2 and Ref. (22)). correspond to forward scattering, involving processes that scatter only between modes with the same direction of motion, refers to backward scattering, and are the umklapp terms which are forbidden in edge states. Note that, unlike in usual ology,(22) we may not assume that the coupling constants for the individual processes are constant. The momentumdependence of the must be taken into account.
The two forward scattering processes and may be merged together into one
Hamiltonian
(7) 
where enforces the normal order(22) of the operator . The primed sum is restricted such that for all momentum arguments in the electron operators. In order to be able to change the amplitude of the momentumdependence of the interaction vertex [see Eq. (5)], we introduce the factors
(8) 
from which we build the interaction vertex for the generalized model. The factor quantifies the momentumdependence. For the interaction is momentumindependent. This limit corresponds to usual onedimensional Hubbard chains. For the interaction goes to zero if at least one of the fermions is close to the upper band edge (). This corresponds to the direct model, where the trigonometric term under the square root in has been replaced by a linear approximation. The differences between the generalized model in the edge state limit and the direct model only lead to quantitative renormalizations of the critical point, as shown in Appendix B. The essential property of the interaction vertex is that it approaches zero if one of the fermion momenta gets close to the Dirac points. This feature is present in the direct and in the generalized model with .
The form of the backscattering Hamiltonian is similar to . However the scattering takes place between left and rightmovers
(9) 
We have introduced the additional parameter which allows us to tune the overall strength of the processes relative to the processes. corresponds to the physical backscattering strength which is required by SU(2) invariance.(22) Nevertheless, we will investigate the consequences of a suppression of backscattering since this will be important for a bosonization analysis of the generalized model which will be presented in an upcoming paper.(24)
As already pointed out, an important feature of edge states is the absence of umklapp processes in the effective electronelectron interaction. However, in order to be able to compare the edge state model to a Hubbard chain, we add an artificial umklapp process with relative strength to the Hamiltonian of the generalized model
(10) 
Varying between 1 (Hubbard chain limit) and 0 (edge state limit) allows us to investigate the consequences of the presence of umklapp processes for onedimensional ferromagnetism.
Alltogether, the four parameters , , , and define the phase space of the generalized model
(11) 
The following limits of this model may be identified:

Edge state limit: the generalized model with the parameters , , and , is a good approximation of the direct model.

Hubbard chain limit: for , , and , the generalized model essentially describes a onedimensional Hubbard chain. The only difference is the assumption of a linearized singleparticle spectrum instead of the dispersion.
Note that it is important to work in the kspace formulation because it is difficult to control the umklapp scattering or the momentum dependence of the interaction vertex in a real space formulation. One reason for this is that an interaction vertex with a nontrivial dependence does not transform to a real space interaction of the form but to a complicated nonlocal interaction. This also hampers the application of DMRG methods to this problem.
Iii Exact diagonalization
The ground state of the generalized model is calculated for finite sized zigzag edges up to by
the Lanczos exact diagonalization method
The Hamiltonian [Eq. (11)] conserves the numbers of upspin and downspin electrons, so that is block diagonal in the subspaces, which we define by the total spinpolarization in direction
(12) 
The total number of electrons is kept constant. This corresponds to halffilling. Note, however, that the filling is physically relevant only if umklapp scattering is present (i.e. ). For the edge states in which we are finally interested, umklapp scattering is forbidden so that the filling is irrelevant as it only leads to quantitative renormalizations of the interaction strength and the Fermi velocity. In the following, we determine the ground state of in each subspace separately.
Note that by the definition of the subspaces we have chosen a spin quantization axis. The Hamiltonian , however, is SU(2) symmetric if the backscattering is at its physical strength . Furthermore, since we are dealing with finite systems, there will be no spontaneous rotational symmetry breaking. Thus, edge magnetism will become manifest in a fold ground state degeneracy, corresponding to a high spin () state. For instance, if in a system with electrons the lowest energy states in the subspaces are the degenerate ground states, then the spins of two electrons point into the same direction, building an super spin. Because the SU(2) symmetry of the individual electron spins is not broken, also this composite super spin has full rotational symmetry. The quantum numbers of the degenerate subspaces then correspond to the magnetization of this composite spin system. Note that the spinorbit interaction lowers the symmetry of the super spin, as it breaks the SU(2) invariance of the individual electron spins which form the super spin.
For practical reasons, we extract the total spin quantum number of the ground state from its ground state degeneracy , which is obtained from the subspace ground state energies. We have checked that this is equivalent to calculating the total spin of the ground state directly.
iii.1 Hubbard chain vs. edge states
First we study the crossover from a usual Hubbard chain to interacting edge states. As explained above, the generalized model can be tuned continuously between these two limiting cases by means of the parameters and . We postpone the analysis of backscattering to the next subsection and set here.
It is most instructive to begin with the Hubbard chain limit of the generalized model, which is
characterized by the full umklapp process strength and a suppressed
momentum dependence of the interaction . With this parameter set the direct model
resembles a onedimensional metal with a linear single particle dispersion instead of a cosdispersion.
Next, the generalized model is tuned away from the Hubbard chain limit by suppressing the umklapp scattering . Suppressed umklapp scattering is one of the properties of edge states which makes them fundamentally different from usual onedimensional metals. In the inset of Fig. 3, the lowest eigen energies of the and the (full spinpolarization) subspaces are shown as is reduced from 1 to 0 in steps of 0.2. For any there is a nonzero critical value for below which the lowest energy states of these two subspaces and also for all in between (not shown in the inset of Fig. 3) are equal. This corresponds to a high spin state of size in which the spins of all electrons point into the same direction. The critical point at which the transition between and takes place depends on the umklapp scattering strength
(13) 
For we find for the exponent . Obviously, the absence of umklapp scattering allows a high spin ground state. However, for and , the system instantly jumps from zero polarization to the maximal possible polarization at the critical point . This is a first order phase transition. For the case of completely suppressed umklapp scattering this is shown in Fig. 4 (a), where the lowest eigen energies of all subspaces are plotted as a function of : the subspace contains the nondegenerate ground state until at the critical point the lowest energy eigenstates of all subspaces form the degenerate ground state; no intermediate regime of exists in which there is only a degeneracy between some of the subspaces.
The reason for the instant jump in the total spin is as follows: once the Stoner criterion is met, the interaction energy gain associated with developing a certain spin polarization is larger than the corresponding kinetic energy penalty . Unlike in two or three dimensions, however, for onedimensional systems with momentumindependent interactions, has no minimum, i.e. , for all . Thus, the system instantly ’flows’ to the highest possible polarization , once the Stoner criterion is met. This is a rather common feature of onedimensional systems with a constant interaction vertex (such as the Hubbard interaction) and can easily be observed in a variational calculation of the ground state properties (see Appendix A).
The momentumdependence of the interaction vertex () reduces the interaction energy
gain as the spinpolarization becomes larger. This is because for larger , the Fermi level of the
spinup rightmovers is shifted to higher momenta where the interaction is suppressed by the factors [see Eq. (7)]. Similarly, for the spinup leftmovers, the Fermi level is then
shifted to smaller momenta, where the suppress the interaction.
If the total spin is plotted as a function of , decreases from to 0 in
steps. These steps correspond to the positions , where the degree of the
ground state degeneracy changes, indicated by arrows in Fig. 4. For , there is one step where the spinpolarization
jumps from to , while for the maximal , there are
steps at each of which the spinpolarization is decreased by .
Note that the nonzero found in the exact diagonalization reveals a weakness of the fermionic meanfield theory in which this minimum momentum dependence, above which a WEM regime appears, is zero (see Appendix A). A nonzero means that the small momentumdependencies which always follow from a dependence of the Bloch wave functions in usual onedimensional conductors on the momentum are not necessarily sufficient to stabilize the weak edge magnetism; the momentumdependence of the interaction vertex must be sufficiently strong for this.
In the limit , which cannot be accessed within exact diagonalization, of course, becomes a smooth function of . We approximate this smooth function by a power law
(14) 
Fermionic meanfield theory (see Appendix A) predicts . Because the exact diagonalization study is limited to small systems , it is difficult to obtain a decent estimate of the exact exponent within this work. Fitting the edge state limit of the generalized model to the center of the plateaus of the results of the exact diagonalization, we obtain , dependent on how many plateaus are included in the fit. The critical point for this fit is obtained by extrapolating the rightmost step from the data sets to .
Interestingly, not only affects the order of the transition but also the critical . This is also not correctly predicted by the meanfield approach (see Appendix A), which, independently of , finds to be the critical point. For small and the exact diagonalization gives
(15) 
For the maximal , the position of the leftmost step can be calculated by exact diagonalization for very large systems.(18) We performed calculations for system sizes up to in order to extrapolate this step position. Within the limits of the accuracy of this extrapolation, the critical point between the SEM and the WEM regime coincides with the meanfield prediction (see Appendix A). This extrapolation to the thermodynamic limit, in combination with the extrapolation of the critical point between the WEM and the LL regime, is a strong evidence for the existence of the WEM phase for in the thermodynamic limit.
For completeness we note that our exact diagonalization analysis shows that a SEM phase also exists in the general model with umklapp scattering if . However even for there is no weak edge magnetism phase between the Luttinger liquid and the saturated edge magnetism as long as .
iii.2 The relevance of backscattering
The backscattering Hamiltonian is important for the SU(2) invariance of the Hamiltonian. It is easily seen that only for the SU(2) symmetry is preserved. At real graphene edges, of course, the backscattering cannot be tuned experimentally. Nevertheless it is interesting to study the consequences of a suppression of since in a bosonization treatment of the generalized model translates to a sineGordon term which is difficult to analyze. Therefore, some insight into the relevance of is helpful from a theoretical point of view. Following the philosophy of the previous subsection, we restrict the discussion to the spin polarization properties of the ground state. The analysis of more complicated observables such as spinspin correlation functions is beyond the scope of this work and will be discussed in another paper.
Fig. 7 compares the lowest eigen energies of the subspaces from calculations with and without backscattering. The most striking feature of the suppression of backscattering is the lifting of the ground state degeneracy in the SEM regime. This effect is easily understood by noting that the very reason for the ground state degeneracy in the case was the SU(2) symmetry, which, however, is broken for . Interestingly, the lifting of the degeneracy is such that the lowest energy states of the subspaces with highest , in the ground space for , form the ground state for . This means that suppressing backscattering introduces an Ising anisotropy along the spin quantization axis chosen in the definition of the model.
Apart from this degeneracy lifting, the evolution of the ground state properties with is very similar for calculations with and without backscattering. The positions of the highest spinpolarization steps are practically unchanged. Only at the steps close to the phase transition between the LL and the WEM regime, a deviation of the results from the results can be observed. A handwaving explanation of this behavior can be given in terms of the bosonization analysis of the WEM regime in Ref. (13). As soon as the Fermi levels for the upspin electrons and the downspin electrons are split, the backscattering process for electrons right at the Fermi surface is forbidden because it is not momentumconserving. Thus, in order to conserve momentum, the electrons are forced to scatter to higher energies if there is a nonzero spinpolarization. This mechanism suppresses the backscattering. In the bosonization language the backscattering Hamiltonian acquires a spatially oscillating phase which makes the corresponding operator irrelevant in the renormalization group. Thus in the WEM regime, not too close to the critical point, is suppressed and does not give an important contribution.
Close to the critical point, however, Fig. 7 indicates that becomes more important. This observation is consistent with the qualitative bosonization argument: Close to the critical point the phase oscillations in the bosonic backscattering Hamiltonian get slower until they completely disappear at the critical point.
Iv Discussion
On the basis of a generalized class of effective models for onedimensional interacting electrons we have studied the magnetic properties of a graphene zigzag edge. Using exact diagonalization we confirmed the existence of three phases within these models, namely the saturated edge magnetism phase which is present at normal graphene edges, the Luttinger liquid phase which appears for edge states with strongly enhanced bandwidth, and an intermediate regime of weak edge magnetism. The latter phase is a realization of a ferromagnetic Luttinger liquid, a onedimensional itinerant ferromagnet. We presented evidence that the transition between the Luttinger liquid and the weak edge magnetism phase becomes a second order quantum phase transition in the limit of long edges.
Beyond the identification of the magnetic properties of edge states, we examined the question why electrons in onedimensional edge states have such a rich phase diagram with two types of ferromagnetic ground states, while usual onedimensional electrons do not show any ferromagnetism. In view of the LiebMattis theorem,(14) which actually forbids a spinpolarized ground state for interacting electrons in one dimension, this question becomes even more pressing.
A closer inspection of the edge state model, which was derived directly from the graphene crystal structure,(13) revealed two unusual features of edge states which cannot be found in other onedimensional electronic systems, such as quantum wires, for instance. These are (a) the total absence of umklapp processes in the effective electronelectron interaction, independently of the filling factor, and (b) a strong momentum dependence of the effective interaction vertex. Each of these features precludes the applicability of the LiebMattis theorem. The momentumdependence gives rise to a complicated nonlocal interaction, which cannot be written as , with an electron density operator, as it is required for the LiebMattis theorem. And even in the limit of momentumindependent interactions ( in the general model) the suppressed umklapp scattering makes the reformulation as a densitydensity interaction in real space impossible. In order to further track down the particular consequences of these special features for the magnetic properties of graphene edges, we replenished the direct model with a tunable umklapp scattering term and replaced the interaction vertex by a generalized vertex function in which the momentumdependence can be switched on and off.
The study of this generalized model, which can be tuned continuously between its edge state limit and a regime in which it describes normal onedimensional metals, revealed the significance of the two edge state features: the absence of umklapp processes is responsible for the existence of a spinpolarized ground state, and the strong momentum dependence of the interaction vertex stabilizes a regime of weak edge magnetism and gives rise to a second order phase transition between the paramagnetic Luttinger liquid and the ferromagnetic Luttinger liquid.
It is interesting to note that the stabilization of the weak edge magnetism phase seems to be very robust against changes in the details of the interaction vertex function. Apparently it is only important that the vertex is suppressed as one of the four momenta of the participating fermions gets close to one of the Dirac points. The exact functional form of this suppression, however, seems to be irrelevant, since the qualitative behavior of the spinpolarization did not depend on whether we used the interaction vertex of the direct model or the interaction vertex of the generalized model with maximal momentum dependence. These two vertex functions have in common that they vanish as one of the momenta approaches a Dirac point. However, their functional forms are very different.
Finally we note that onedimensional itinerant magnetism has also been studied in Hubbard chains with an additional second neighbor hopping,(30); (31) showing that it is indeed possible to define onedimensional models which, at first sight, seem to comply with the prerequisites of the LiebMattis theorem, but nevertheless have a high spin ground state. The physical picture behind the model discussed in Refs. (30); (31), however, is much different from the present work. Interestingly, the sign of the hopping amplitude to the nearest neighbor must be different from the sign of the nextnearest neighbor hopping for the system to have a ferromagnetic ground state.
Acknowledgements.
D.J.L. and F.F.A. acknowledge financial support from the DFG for grant AS120/43. M.J.S. acknowledges financial support from the Swiss NSF and from the NCCR QSIT.Appendix A Variational analysis of the generalized model
We calculate the magnetic ground state properties of the generalized model within a fermionic meanfield approximation. It is assumed that only the averages are nonzero, so that the umklapp Hamiltonian and the backscattering Hamiltonian drop out of the meanfield treatment. The resulting noninteracting Hamiltonian is diagonal in the momentum , in the direction of motion and in the spin projection, so that the meanfield theory is equivalent to a variational ansatz based on the trial wave function
(16) 
with an asymmetric occupation of spinup and spindown states. The variational parameter is related to the spindependent Fermi levels by
(17) 
and to the spinpolarization , used in Sect. III, by . For finite size systems, as discussed in the main part of this paper, the Fermi level cannot be varied continuously so that also is a discrete variable in this case. However, within meanfield theory it is easy to perform the calculations in the thermodynamic limit, so that we will consider to be a continuous variable and interpret it as the magnetization order parameter.
The variational energy is easily calculated from the Hamiltonian in Eq. (11)
(18) 
For , the minimum of is at , while for the meanfield ground state has a finite magnetization
(19) 
Note that by definition the magnetization cannot become larger than 1. From Eq. (19) it becomes obvious that a nonzero momentumdependence is required to stabilize the regime of weak edge magnetism. For , the magnetization would jump from 0 to 1 at the critical point .
The existence of the weak edge magnetism can be traced back to the term in Eq. (18) which is generated by the momentumdependence of the interaction vertex. In dimensions higher than one, such terms emerge also from momentumindependent interactions or directly from the kinetic energy, so that at least on the meanfield level is required for the stabilization of weak ferromagnetism only in one dimension.
Appendix B Exact diagonalization of the direct model
The direct model Hamiltonian defined by Eqs. (3  4) and the edge state limit of the generalized model Hamiltonian [Eq. (11) with , and ] are not exactly equal, as the general model linearizes the single particle dispersion and replaces the factors by the approximation . Nevertheless, the most important properties of graphene edge states, i.e. the momentumdependence of the interaction vertex and the absence of umklapp scattering, are properly described by both, the direct model and the general model in the edge state limit.
In this appendix, we check that the direct model has qualitatively the same magnetic properties as the general model in the edge state limit. In Fig. 8, we present the spinpolarization as a function of , obtained from the exact diagonalization of the direct model. The bandwidth parameter of the direct model corresponds to the Fermi velocity of the general model.
Clearly, for larger , which corresponds to the parameter in the generalized model, we obtain a Luttinger liquid phase with a ground state of total spin . An intermediate regime with weak edge magnetism exists, where the total spin of the ground state is not maximal. As in the exact diagonalization analysis of the general model in the main text, only some of the lowest eigen energies in different subspaces are degenerate and form the ground state. For small , the saturated edge magnetism phase is reached and the spin of the ground state is maximal, i.e., the lowest eigen energies in all subspaces are degenerate.
Note that the energy of the fully spin polarized eigenstate of the direct model Hamiltonian has a finite slope (see Fig. 8). This is because the direct model lacks a symmetry of of the generalized model leading to (cf. Fig. 4). As only the degeneracy of the lowest eigen energies are important, but not their absolute values, this difference does not have any physical consequences.
Footnotes
 It is only important to assume that this energy scale is large enough so that the restriction to the zero energy sector of is well justified. The validity of this assumption has been checked in Ref. (13).
 The equal strength of the and processes is due to the symmetry property , which is a consequence of the origin of in the 2D honeycomb Hubbard model (see Ref. (13)).
 An edge with unit cells in length corresponds to only kspace points in the reduced Brillouin zone in which the edge states are defined.
 We have checked that the linearization of the single particle spectrum does not change the results qualitatively.
 The point corresponds to a pathologic potential in Ref. (14), as it can be reached by . At this point, all sectors with total spin are degenerate.
 Note that the increase in the interaction vertex for the spindown electrons with lowered Fermi level is overcompensated by the suppression due to the higher Fermi level of the spinup electrons, so that in total the interaction is reduced as grows.
 This decrease of in is due to the contribution of the left and rightmoving branch to the spinpolarization: each branch contributes one spin flip.
References
 K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 (2004)
 A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (Jan 2009)
 C. Lee, X. Weil, J. W. Kysar, and J. Hone, Science 321, 385 (2008)
 Y. Zhang, Y.W. Tan, H. L. Stormer, and P. Kim, Nature 438, 201 (Nov 2005), ISSN 00280836, http://dx.doi.org/10.1038/nature04235
 K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature 438, 197 (Nov 2005), ISSN 00280836
 M. Fujita, K. Wakabayashi, K. Nakada, and K. Kusakabe, Journal of the Physical Society of Japan 65, 1920 (1996), http://jpsj.ipap.jp/link?JPSJ/65/1920/
 Y.W. Son, M. L. Cohen, and S. G. Louie, Phys. Rev. Lett. 97, 216803 (Nov 2006)
 Y.W. Son, M. L. Cohen, and S. G. Louie, Nature 444, 347 (2007)
 J. Jung and A. H. MacDonald, Phys. Rev. B 79, 235433 (Jun 2009)
 H. Feldner, Z. Y. Meng, A. Honecker, D. Cabra, S. Wessel, and F. F. Assaad, Phys. Rev. B 81, 115416 (Mar 2010)
 H. Feldner, Z. Y. Meng, T. C. Lang, F. F. Assaad, S. Wessel, and A. Honecker, “Dynamical signatures of edgestate magnetism on graphene nanoribbons,” (2011), arXiv:1101.1882
 T. Hikihara, X. Hu, H.H. Lin, and C.Y. Mou, Phys. Rev. B 68, 035432 (Jul 2003)
 M. J. Schmidt and D. Loss, Phys. Rev. B 82, 085422 (Aug 2010)
 E. Lieb and D. Mattis, Phys. Rev. 125, 164 (Jan 1962)
 C. Lanczos, J. Res. Nat. Bur. Stand. 45, 255 (1950)
 E. Dagotto, Rev. Mod. Phys. 66, 763 (Jul 1994)
 D. Sénéchal, “An introduction to quantum cluster methods,” (2008), arXiv:0806.2690v2
 At filling, and , we could access edges of length with a dimension of the corresponding Hilbert space of after exploiting all symmetries. For nearly maximal and filling, we were able to calculate the groundstate energy of an edge with length with a Hilbert space dimension of only . In this case we were limited only by the 64 bit limit of the internal storage of the basis vectors and longer edges are in principle accessible by a modification of our code.
 It is only important to assume that this energy scale is large enough so that the restriction to the zero energy sector of is well justified. The validity of this assumption has been checked in Ref. \rev@citealpnumtem_schmidt_loss_2010.
 M. J. Schmidt and D. Loss, Phys. Rev. B 81, 165439 (Apr 2010)
 M. Hohenadler, T. C. Lang, and F. F. Assaad(2011), arXiv:1011.5063
 T. Giamarchi, Quantum Physics in One Dimension (Oxford Univ. Press, 2003)
 The equal strength of the and processes is due to the symmetry property , which is a consequence of the origin of in the 2D honeycomb Hubbard model (see Ref. \rev@citealpnumtem_schmidt_loss_2010).
 M. J. Schmidt(2011)
 An edge with unit cells in length corresponds to only kspace points in the reduced Brillouin zone in which the edge states are defined.
 We have checked that the linearization of the single particle spectrum does not change the results qualitatively.
 The point corresponds to a pathologic potential in Ref. \rev@citealpnumlieb_mattis, as it can be reached by . At this point, all sectors with total spin are degenerate.
 Note that the increase in the interaction vertex for the spindown electrons with lowered Fermi level is overcompensated by the suppression due to the higher Fermi level of the spinup electrons, so that in total the interaction is reduced as grows.
 This decrease of in is due to the contribution of the left and rightmoving branch to the spinpolarization: each branch contributes one spin flip.
 S. Daul, Eur. Phys. J. B 14, 649 (2000)
 S. Daul and R. M. Noack, Phys. Rev. B 58, 2635 (1998)