Interplay between the edge-state magnetism and long-range Coulomb interaction in zigzag graphene nanoribbons: quantum Monte Carlo study

# Interplay between the edge-state magnetism and long-range Coulomb interaction in zigzag graphene nanoribbons: quantum Monte Carlo study

## Abstract

We perform projective quantum Monte Carlo simulations of zigzag graphene nanoribbons within a realistic model with long-range Coulomb interactions. Increasing the relative strength of nonlocal interactions with respect to the on-site repulsion does not generate a phase transition but has a number of nontrivial effects. At the single-particle level we observe a marked enhancement of the Fermi velocity at the Dirac points. At the two-particle level, spin- and charge-density-wave fluctuations compete. As a consequence, the edge magnetic moment is reduced but the edge dispersion relation increases in the sense that the single-particle gap at momentum grows. We attribute this to nonlocal charge fluctuations which assist the spin fluctuations to generate the aforementioned gap. In contrast, the net result of the interaction-induced renormalization of different energy scales is a constant spin-wave velocity of the edge modes. However, since the particle-hole continuum is shifted to higher energies—due to the renormalization of the Fermi velocity—Landau damping is reduced. As a result, a roughly linear spin-wave-like mode at the edge spreads out through a larger part of the Brillouin zone.

## I Introduction

Quasi-one-dimensional graphene nanoribbons terminated by zigzag edges feature localized edge states at the Fermi energy. Nakada et al. (1996); Wakabayashi et al. (1999); Brey and Fertig (2006) Mean-field studies of the Hubbard model, with an on-site effective interaction only, Fujita et al. (1996); Fernández-Rossier (2008); Jung et al. (2009) as well as density-functional theory calculations, Son et al. (2006a); Pisani et al. (2007); Yazyev and Katsnelson (2008) predict spontaneously induced spin polarizations at the zigzag edges. The resulting finite dispersion of the low-energy electronic states is usually considered in the scanning tunneling spectroscopy as a signature of edge-state magnetism, Tao et al. (2011) despite merely a weak paramagnetic signal measured in the bulk graphene samples. Sepioni et al. (2010) The distinct low-energy electronic structure and the consequent transport properties of zigzag nanoribbons have been the subject of great interest. Wakabayashi et al. (2009); Yazyev (2010) Moreover, progress in atomic engineering makes it possible nowadays to build nanoribbons with atomically precise zigzag edges Ruffieux et al. (2016) and to verify an early promising proposal of their applicability in nanoelectronics. Son et al. (2006b); Rycerz et al. (2007)

The edge magnetism of zigzag graphene nanoribbons has also been extensively studied beyond the mean-field level using the density matrix renormalization group, Hikihara et al. (2003); Hagymási and Legeza (2016) exact diagonalization, Luitz et al. (2011) and quantum Monte Carlo (QMC) simulations. Feldner et al. (2011) A substantial ferromagnetic (FM) spin correlation length along the zigzag edge found in the latter study already for ribbons of moderate widths, i.e., with four zigzag chains, justifies the picture obtained within the mean-field type approximations. Feldner et al. (2011) In addition, an effective antiferromagnetic (AF) coupling between the opposite edges of the nanoribbon guarantees the spin-singlet nature of the ground state in accordance with Lieb’s theorem for bipartite lattices with a half-filled band. Lieb (1989) Thus, zigzag graphene nanoribbons constitute an interesting many-body system with a novel type of boundary-critical phenomena. Karimi and Affleck (2012)

However, in the semi-metallic state, as realized in graphene, the Coulomb interaction is only partially screened and remains long-range. Kotov et al. (2012) As shown within the leading-order perturbative renormalization group analysis González et al. (1994, 1999) and confirmed by more elaborated calculations, Barnes et al. (2014); Hofmann et al. (2014); Bauer et al. (2015); Sharma and Kopietz (2016) this results in a marginal Fermi liquid behavior and a logarithmic renormalization of the electron velocity at the Dirac points. The predicted renormalization of the Fermi velocity has been observed in suspended graphene using cyclotron resonance Elias et al. (2011) and has triggered intense theoretical studies of the effects driven by nonlocal Coulomb interactions in graphene and related materials. Araki and Semenoff (2012); Wu and Tremblay (2014); Smith and von Smekal (2014); Golor and Wessel (2015); Tang et al. (2015); Faye et al. (2015); Classen et al. (2016); Volpez et al. (2016); Katanin (2016); de la Peña et al. (2017); Tupitsyn and Prokof’ev (2017) A detailed analysis of this fundamental problem that includes the effect of both the short- and long-range parts of the Coulomb repulsion on the velocity renormalization can be found in Ref. Tang et al., 2017.

Given the long-range character of the Coulomb interaction in graphene, it stands to reason that there is a need to revisit the stability and dynamical signatures of spin-polarized edge states of zigzag graphene nanoribbons in the presence of nonlocal interactions. Only very recently, a step towards a realistic description of the nanoribbons beyond the mean-field treatment of Coulomb interactions Yamashiro et al. (2003); Wunsch et al. (2008); Jung (2011) was taken by Shi and Affleck. Shi and Affleck (2017) By neglecting the dynamics of low-lying bulk states, they derived an effective Hamiltonian with nonlocal interactions projected onto the Hilbert space of edge states. Within this framework, they were able to provide a sufficient condition under which the edge magnetism is stable against unscreened Coulomb interactions.

The aim of this paper is to validate those conclusions by performing lattice simulations of the full Hamiltonian with the unscreened tail of the electron-electron repulsion. To this end, we carry out systematic studies of zigzag graphene nanoribbons with up to lattice sites using the unbiased QMC technique capable of dealing with long-range density-density interactions. Hohenadler et al. (2014) Our results lend further support for the robustness of the edge magnetism. We also find that nonlocal interactions play an important role in the renormalization of the dispersion of edge states. Since the screening length in graphene can be efficiently engineered by substrate modification, Hwang et al. (2012) deposition of a zigzag nanoribbon on substrates with a different dielectric constant should allow one to tune accordingly the low-energy electronic states.

The paper is organized as follows: In Sec. II, we define the model and introduce the projective QMC method. In Sec. III, we focus on static properties of the system and address the behavior of equal-time spin and charge correlation functions upon increasing nonlocal Coulomb interactions. In Sec. IV, we analyze the corresponding evolution of single- and two-particle excitation spectra. Section V provides a summary and the conclusions.

## Ii Model and the QMC method

We consider the following Hamiltonian:

 ^H=−t∑⟨ij⟩,σ^c†iσ^cjσ+14∑ijVij(^ni−1)(^nj−1), (1)

where creates an electron with spin at lattice site and is the local particle number operator. The first term in Eq. (1) describes the kinetic energy of the electrons with the hopping amplitude between nearest-neighbor sites on a finite-size half-filled honeycomb lattice with open boundary conditions at the zigzag edges. The width of the graphene nanoribbon is defined by the number of zigzag chains with 2 carbon atoms each and its length is measured in units of the lattice vector where is the carbon-carbon bond length. The second term is a long-range density-density interaction with matrix elements,

 Vij={2U,if |i−j|=0αUa|i−j|,if |i−j|>0. (2)

Here, is the minimal distance between two lattice sites and while determines the relative strength of nonlocal interactions with respect to the on-site repulsion . In particular, for , one recovers the Hubbard interaction:

 ^HU=U2∑i(^ni−1)2. (3)

To address the ground-state properties of the Hamiltonian (1), we use a projective (zero temperature) QMC algorithm based on the imaginary-time evolution of a trial wave function to the ground state ,Assaad and Evertz (2008)

 ⟨Ψ0|^O|Ψ0⟩⟨Ψ0|Ψ0⟩=limΘ→∞⟨ΨT|e−Θ^H/2^Oe−Θ^H/2|ΨT⟩⟨ΨT|e−Θ^H|ΨT⟩. (4)

It relies on a Trotter-Suzuki decomposition of the imaginary-time propagation,

 e−Θ^H=Lτ∏i=1e−Δτ^Hte−Δτ^HV+O(Δτ2), (5)

into steps of size . Numerical results presented in Secs. III and IV were obtained with the projection parameter , chosen sufficiently large to guarantee the convergence to the ground state , and with an imaginary time discretization of .

The simulations at finite were performed by means of a recently implemented generalization of the QMC method which allows one to handle long-range density-density interactions. Hohenadler et al. (2014) The approach we use, was introduced and used in the context of hybrid QMC methods. Brower et al. (); Ulybyshev et al. (2013) It is based on the observation that the quartic nonlocal interaction term can be recast into a quadratic one by coupling the local electron density operator to an auxiliary continuous scalar field, provided the interaction matrix is positive definite. This renders it possible to integrate out fermionic degrees of freedom at the expense of a determinant yielding an action solely dependent on the scalar field. The integration over the field configurations is carried out stochastically using the Monte Carlo method with a sequential Metropolis updating scheme. Finally, a standard discrete Hubbard-Stratonovich transformation was employed in the limit of a purely local Hubbard interaction. Assaad and Evertz (2008) In order to ensure the explicit conservation of the SU(2) spin-rotation symmetry at half-filling, we have opted for an Ising spin auxiliary field coupling to the local charge operator at every space-time point,

 e−ΔτU(^ni−1)2=12∑siτ=±1eiλsiτ(^ni−1), (6)

with coupling constant given by .

## Iii Static spin and charge correlations

Nonlocal Coulomb interactions can make charge fluctuations as large as the spin fluctuations. Zhang and Callaway (1989); Ohta et al. (1994); Aichhorn et al. (2004); Davoudi and Tremblay (2007); van Loon et al. (2014); Huang et al. (2014); Terletska et al. (2017) This contrasts with systems with the repulsive on-site interaction only where the spin fluctuations dominate. In the extreme case, increasing the range of Coulomb interactions can result in a phase transition from a uniform state favoring singly occupied sites and hence, local moment formation, to the charge-ordered state with alternating doubly occupied and empty sites.

On the one hand, starting from the 1D and strong-coupling limits, the critical nearest-neighbor Coulomb repulsion driving the spin-density wave to a charge-density-wave phase transition along the 1D edge is larger than in the bulk graphene with the coordination number 3. Thus, this oversimplified view implies robustness of the edge magnetism. On the other hand, Coulomb interaction couples edge and bulk states and the bulk corrections to the edge physics become more important with increasing ribbon width Koop and Schmidt (2015) This coupling makes the issue of the stability of magnetic edge correlations less obvious and might result in a complex many-body state with intertwined spin and charge fluctuations both contributing to the spectral properties of the nanoribbons.

Given this background, we perform QMC studies with the purpose of elucidating to what extent enhanced charge fluctuations affect the edge magnetism upon increasing the strength of nonlocal interactions controlled by . The QMC simulations have been done at a fixed value of the on-site repulsion , well below the transition from the semimetal to the bulk AF insulator at in the honeycomb Hubbard model, Sorella et al. (2012); Assaad and Herbut (2013) and up to where charge and spin correlations are nearly degenerate in the classical limit. Notably, as estimated from previous QMC simulations, already increases substantially the critical interaction for the aforementioned magnetic transition signifying the importance of nonlocal interactions. Hohenadler et al. (2014)

First, we study properties of the ground state of the system. To this end, we calculate the static real-space spin-spin correlation function,

 S(r)=43⟨^Sr⋅^S0⟩, (7)

and charge-charge correlation function,

 N(r)=⟨^nr^n0⟩−⟨^nr⟩⟨^n0⟩. (8)

The evolution of the real-space correlation pattern which stems from the QMC simulations of our largest nanoribbon with width and length , measured relative to the central site of the edge (indicated by the arrow), is illustrated in Fig. 1. For a purely local interaction () one finds a clear evidence of quasi-long range FM spin correlations along the edges, see Fig. 1(a). The spin-spin correlations are antialigned to each other on the opposite ribbon edges and quickly decay into the inner sites. The corresponding charge-charge correlations shown in Fig. 1(b) are slightly suppressed around the origin with respect to the noninteracting system. Such behavior is expected for a fluid of charge carriers with contact interactions where charged particles avoid each other at short distances. In the case of a long-range interaction with , numerical data shown in Fig. 1(c,d) indicate the competition between short-range charge correlations, driven by the nonlocal interactions, and AF spin correlations promoted by the on-site repulsion. Consequently, the spin correlations fall off into the bulk sites more quickly as compared to the Hubbard model, see Fig. 1(c). Nevertheless, the tendency towards the extended spin polarization along the edges remains dominant over the sublattice charge fluctuations observed in , see Fig. 1(d).

As summarized in Table 1 for ribbons of width and 8, the enhancement of charge fluctuations upon increasing is also signaled by a systematic growth of the average double occupancy in the system,

 d=1N∑i⟨^ni↑^ni↓⟩. (9)

The increase of can be equally well understood as resulting from the tendency of nonlocal Coulomb interactions to renormalize downward the effective on-site interaction. Schüler et al. (2013) Note that enhanced error bars at the smallest nonzero follow from particularly low acceptance rates within the adopted sequential Metropolis updating scheme of the time slices. However, the efficiency of the implemented sampling of the scalar field improves at larger which is reflected in better data quality.

In order to address the edge magnetism in the presence of long-range Coulomb interactions more quantitatively, we plot in Fig. 2 the spin-spin correlation function as a function of the distance along the edge. Following up Ref. Feldner et al., 2011, we consider a function,

 Sfit(r)∝r−ηe−r/ξ+(L−r)−ηe−(L−r)/ξ, (10)

combining both exponential and power-law behaviors. An important conclusion of the QMC studies in Ref. Feldner et al., 2011 reached on systems up to was that the correlation length increases with the width of ribbons such that already the data for the ribbons can be equally well fitted assuming an infinite spin correlation length. By fitting to our QMC data obtained on shorter systems with , we find exceedingly large magnetic correlation length even for the narrowest ribbon with . However, as summarized in Table 2, the actual change of exponent with increasing strength of the nonlocal part of Coulomb interactions is the width-dependent.

In the case of the narrow ribbons, finite enhances the interedge AF coupling and brings the system closer to the regime of even-leg spin ladders as signals by a larger with respect to exponent and thus a more rapid exponential-like decay of spin correlations. White et al. (1994); Dagotto and Rice (1996) The tendency to form rung singlets becomes less important for larger interedge distances. Already for the ribbons, the dominant effect is the reduction of short-range spin correlations by competing charge-density fluctuations while leaving the long-wavelength physics nearly unaltered. This is seen as an effective reduction of the exponent with , cf. Table 2. Same behavior is also resolved for the ribbons.

The special character of the ribbon with the relatively strongest reduction of spin correlations upon increasing is further confirmed by comparing in Fig. 3(a) the evolution of the magnetic moment at the largest distance along the edge,

 mL/2=√43⟨^SL/2⋅^S0⟩, (11)

with that obtained for wider and 8 ribbons. In contrast, the reduction of the local moment at the edge,

 m0=43⟨^S20⟩=1−2⟨^n0↑^n0↓⟩, (12)

is very much the same for all the ribbon geometries studied in this work, see Fig. 3(b).

Equation (12) implies that a decrease in the local moment is equivalent to the enhancement of double occupancy and reflects growing charge fluctuations in the system. Thus, we study next the evolution of spatial dependence of double occupancy upon increasing . To this end, we plot in Fig. 4 double occupancy for the widest ribbon at the edge as well as at its nearest-neighbor (NN) and next-nearest-neighbor (NNN) sites in the unit cell. For comparison, we also show double occupancy in the bulk. On the one hand, due to the tendency to local moment formation, the edge clearly stands out from the rest through the strongest reduction of double occupancy, see Fig. 4(a). On the other hand, double occupancy of both the NN and NNN sites is essentially indistinguishable from the bulk and becomes at our largest only slightly reduced compared to its noninteracting value 0.25, see Fig. 4(b).

## Iv Dynamical spectral functions

Given the tendency of nonlocal Coulomb interactions to renormalize downwards the effective on-site interaction, Schüler et al. (2013) the reduction of spin correlations along the edge by nonlocal interactions revealed in Sec. III can be mimicked within a simplified model with a purely local interaction. However, an important question is how accurately the actual physical situation in the system can be accounted for by such an effective model?

To address this issue let us now discuss the single- and two-particle excitation spectra of the system. Using the translational invariance of the system along the edge direction, it is convenient to introduce a momentum-space representation of the fermion operator,

 ^cℓqσ=1√LL∑j=1eiq|a1|j^ci(j,ℓ)σ, (13)

where the lattice site is given by the unit cell index and the intracell index . The knowledge of the momentum-resolved Green’s function,

 Gℓ(q,τ)=12∑σ⟨^cℓqσ(τ)^c†ℓqσ(0)⟩, (14)

allows one to obtain the single-particle spectral function on the real frequency axis by the inversion of,

 Gℓ(q,τ)=1π∫dωe−τωAℓ(q,ω), (15)

by means of the stochastic analytic continuation method. Beach (2004) Similarly, imaginary-time-displaced spin-spin and charge-charge correlation functions allow one to extract the dynamical spin and charge structure factors by analytic continuation of,

 43⟨^Sℓ(q,τ)⋅^Sℓ(−q,0)⟩ =1π∫dωe−τωSℓ(q,ω), (16) ⟨^Nℓ(q,τ)⋅^Nℓ(−q,0)⟩ =1π∫dωe−τωCℓ(q,ω), (17)

where,

 ^Sℓ(q) =1√LL∑j=1eiq|a1|j^Si(j,ℓ), (18) ^Nℓ(q) Missing or unrecognized delimiter for \left (19)

and is the average filling level. Below we discuss dynamical spectral functions along the ribbon edge and within the center (bulk) obtained for our largest ribbon of width .

### iv.1 Single-particle spectral function

The essential feature of the tight-binding () band structure of a semi-infinite graphene sheet with a zigzag edge is a pair of flat zero-energy bands in the range between the two Dirac points. Electron interactions modify the single-particle spectrum and lead to distinct spectral features related to the onset of edge magnetism. We illustrate this effect in Fig. 5. It compares a momentum resolved single-particle spectral function along the zigzag edge in the limit of local () Hubbard interaction [Fig. 5(a)] with that obtained with nonlocal interactions [Fig. 5(b)]. Note that the considered value is chosen below the critical value for the Mott transition on the honeycomb lattice Sorella et al. (2012); Assaad and Herbut (2013) Hence, on our limited-size lattice, the gap around and wavevectors is a finite size effect. Moreover, the observed in Ref. Hohenadler et al., 2014 increase of the critical value upon going from a Hubbard to a long-range interaction, ensures the gapless nature of edge modes such that the system remains a semimetal.

As apparent, in the presence of the Hubbard interaction [Fig. 5(a)], the initially flat edge-localized states acquire a finite dispersion reaching its maximum at the momentum. The position of the single-particle peak at the wavevector is proportional to both the strength of the Hubbard interaction and the magnetic moment at the edge sites, Feldner et al. (2011) . As shown in Fig. 5(b), switching on nonlocal interactions renormalizes further the dispersion of edge states and gives rise to the enhancement of its bandwidth. While it is suggestive of the many-body renormalization of the electron velocity at low energies found in bulk graphene, Elias et al. (2011) we note that much larger system sizes than those studied in this work are required to unambiguously resolve the logarithmically divergent Fermi velocity for edge modes near the Dirac points. Shi and Affleck (2017)

In Fig. 6 we display the corresponding single-particle spectrum in the bulk . A complete suppression of the low-energy spectral weight for wavevectors can be traced back to the localized nature of the edge states which quickly decay into the inner sites. Furthermore, upon going from a Hubbard [Fig. 6(a)] to a long-range interaction [Fig. 6(b)], one observes a substantial widening of the single-particle spectrum.

To get more insight into this effect and, in particular, to ascertain whether the observed renormalization of the single-particle spectral function at the ribbon edge and in the bulk can be traced to the same origin, we focus on the momentum. This wavevector is of major relevance to magnetic properties of the system: it corresponds to the most localized states along the edge and thus—according to the Stoner criterion—plays a dominant role in the spontaneous appearance of FM spin correlations along the edges.

A detailed evolution of and upon increasing for a ribbon of width is summarized in Fig. 7. We also plot the spectral data from the QMC simulations on narrower and ribbons. First of all, the observed transfer of spectral weight towards higher frequencies turns out to be robust against changes in the nanoribbon width excluding the possibility that it simply stems from a finite size effect enhanced by long-range Coulomb interactions. Furthermore, the striking difference between both panels is a steady shift in the peak position in which persists even at our largest . Thus, one may conclude that the shift originates from the long-range tail of Coulomb interactions; the resultant renormalization of the electron velocity at the Dirac points produces a rigid shift of the electronic structure in the bulk towards higher energies, see Fig. 6.

In contrast, the corresponding transfer of spectral weight in the single-particle spectrum at the edge seems to saturate with increasing . This indicates that it stems from short-range interactions. Such an interpretation is further supported by the results of Ref. Shi and Affleck, 2017 reporting a significant change in the position of the low-energy peak in the presence of Coulomb interactions with a short-range real-space cutoff. However, recalling that in the Hubbard model, , the simultaneous shift of towards higher frequencies and the reduction of the magnetic moment cannot be consistently reproduced in the effective model with a renormalized on-site interaction thus pointing out the intricate nature of the dynamics of edge states.

### iv.2 Dynamical charge structure factor

In order to resolve this puzzle, we turn now our attention to the charge dynamics. From the equal-time data in Fig 1(d), revealing enhanced short-range charge correlations at finite , one might expect to resolve the concomitant redistribution of charge spectral weight in the dynamical charge structure factor along the zigzag edge . The comparison of obtained in the Hubbard model [Fig. 8(a)] with that calculated with long-range interactions [Fig. 8(b)] does indeed reveal the emergence of broad incoherent low-energy excitations. The observed softening near and wavevectors stems from scattering between the Dirac points and reflects growing charge fluctuations in the system. Therefore, it is natural to ascribe the shift of in Fig. 7(a) to short-range charge-density-wave fluctuations. In this scenario, the quasiparticle scattering off charge fluctuations adds up to scattering off spin fluctuations and effectively reinforces the edge magnetism induced single-particle gap at the wavevector. Finally, we also note that the quickly vanishing low-energy charge weight with increasing momentum in Fig. 8(b) is consistent with the expected increase of the velocity of long-wavelength charge excitations in the presence of nonlocal interactions.

Given the semimetallic nature of the bulk, it is expected to retrace soft features seen in in the charge excitation spectrum in the bulk. However, the corresponding dynamical charge structure factor in Fig. 9 appears to be gapfull. The apparent contradiction is resolved by comparing the single-particle spectra along the edge (Fig. 5) and in the bulk (Fig. 6). The presence of a weakly dispersive quasiparticle band along the edge directly affects the low-energy scattering between the Dirac points. The resulting large phase space for particle-hole excitations gives rise to an extended particle-hole continuum in the dynamical charge structure factor along the edge , see Fig. 8(b). In contrast, a pointlike Fermi surface with vanishing density of states in the bulk strongly restricts decay channels for low-energy quasiparticles due to phase space limitations. Consequently, the intensity of the particle-hole continuum quickly decreases in the low-energy limit making it very difficult to resolve the corresponding charge excitations in , see Fig. 9(b).

### iv.3 Dynamical spin structure factor

Finally, let us address the effect of nonlocal Coulomb interactions on the spin dynamics. Figure 10(a) reports the dynamical spin structure factor along the edge in the Hubbard model. As apparent, a prominent signature of quasi-long range FM spin correlations is the roughly linear low-energy spin-wave-like mode which merges into a particle-hole continuum at higher frequencies. Considering the FM nature of spin correlations along the edge, it might—at first glance—be surprising to see that the spin-wave dispersion is not quadratic in the long-wavelength limit. However, the linear dispersion is a consequence of AF spin correlations between the opposite edges. Wakabayashi et al. (1998); Culchac et al. (2011) This behavior has a close analogy with spin excitations in -type Heisenberg antiferromagnets. They are composed of one-dimensional chains with a FM superexchange embedded in a higher-dimensional system with AF couplings along the other directions. In such a system, the AF interactions dominate the actual form of spin-wave excitations: one finds linear Goldstone modes irrespective of the cubic direction while the FM interactions modify merely the spin-wave velocity Raczkowski and Oleś (2002)

A comparison with Fig. 10(b) reveals that the magnetic spectrum is directly affected by nonlocal interactions with spin-wave-like excitations continuing to disperse up to larger wavevectors. We attribute this to the transfer of spectral weight in the single-particle spectrum at the edge in Fig. 7(a) which reduces Landau damping at low frequencies. In contrast, despite the interaction-induced renormalization of the electron velocity at the Dirac points, which should in turn increase the spin-wave velocity, the spin spectral function remains essentially identical for small momenta , see the inset in Fig. 10(b). This can be understood as resulting from the simultaneous reduction of the effective on-site repulsion by nonlocal interactions (or equivalently the reduction of the magnetic moment at the edge in Fig. 3(a) by enhanced charge fluctuations) such that both renormalization effects effectively cancel each other out.

Figure 11 displays the dynamical spin structure factor along the bulk sites . In contrast to the edge spin spectrum featuring a low-energy spin-wave-like mode, the bulk spin spectral function is broad and corresponds to the particle-hole continuum. In addition to the low-energy spectral weight in the long-wave limit, the spin excitations seem to also soften on approaching and wavevectors. These two momenta correspond to wavevectors connecting two Dirac cones in the single-particle spectral function , see Fig. 6. The softening at finite momenta bears a striking similarity to the continuum of excitations with soft features around and wavevectors resolved in the dynamical charge structure factor along the zigzag edge , see Fig. 8(b). The similarity between the low-energy spin and charge excitations at finite momenta confirms that they stems from scattering between the Dirac points.

## V Summary

By performing projective QMC simulations of a realistic model with long-range Coulomb interactions, we have revisited the stability of spin-polarized edge states of graphene nanoribbons with up to lattice sites. The tendency towards the extended FM spin polarization along the edges promoted by the on-site interaction remains dominant over the competing short-range charge correlations favored by long-range Coulomb interactions. Our results agree with the findings of Shi and Affleck using an effective edge model, Shi and Affleck (2017) providing further evidence in support of edge magnetism for moderate values of nonlocal interactions.

While the reduction of spin correlations along the edge can be understood as resulting from the tendency of nonlocal Coulomb interactions to renormalize downwards the effective on-site interaction, the evolution of the single-particle spectrum cannot be accounted for within a simplified model with a purely local Hubbard interaction. Indeed, on the basis of our QMC simulations we have found that spin and charge fluctuations are intrinsically mixed and both are of vital importance to the actual form of the dispersion relation of edge states. In particular, increasing the relative strength of nonlocal interactions with respect to the on-site repulsion tunes the single-particle gap at the momentum. Thus, given that the screening length in graphene can be efficiently engineered by substrate modification, Hwang et al. (2012) deposition of a zigzag nanoribbon on substrates with a different dielectric constant should allow one to modify accordingly the low-energy electronic states.

Our final remark concerns doped zigzag nanoribbons. It has been proposed within the Hartee-Fock approximation of the Hubbard model that the extra charge will accommodate at the edges forming topological defects, i.e., domain walls, across which there is a change in the phase of the edge spins. López-Sancho and Brey (2017) Given the tendency for charge topological solitons in doped nanoribbons, the inclusion of the long-range part of the Coulomb interaction shall drive the extra charge apart, possibly affecting the charge concentration in the domain walls and thus the distance between solitons. A more quantitative discussion requires a separate analysis which constitutes an interesting subject for future.

###### Acknowledgements.
This work was supported by the German Research Foundation (DFG) through SFB 1170 ToCoTronics. The authors gratefully acknowledge the computing time granted by the John von Neumann Institute for Computing (NIC) and provided on the supercomputer JURECA Jülich Supercomputing Centre (2016) at Jülich Supercomputing Centre (JSC).

### References

1. K. Nakada, M. Fujita, G. Dresselhaus,  and M. S. Dresselhaus, Edge state in graphene ribbons: Nanometer size effect and edge shape dependence, Phys. Rev. B 54, 17954 (1996).
2. K. Wakabayashi, M. Fujita, H. Ajiki,  and M. Sigrist, Electronic and magnetic properties of nanographite ribbons, Phys. Rev. B 59, 8271 (1999).
3. L. Brey and H. A. Fertig, Electronic states of graphene nanoribbons studied with the Dirac equation, Phys. Rev. B 73, 235411 (2006).
4. M. Fujita, K. Wakabayashi, K. Nakada,  and K. Kusakabe, Peculiar Localized State at Zigzag Graphite Edge, J. Phys. Soc. Jpn. 65, 1920 (1996).
5. J. Fernández-Rossier, Prediction of hidden multiferroic order in graphene zigzag ribbons, Phys. Rev. B 77, 075430 (2008).
6. J. Jung, T. Pereg-Barnea,  and A. H. MacDonald, Theory of Interedge Superexchange in Zigzag Edge Magnetism, Phys. Rev. Lett. 102, 227205 (2009).
7. Y.-W. Son, M. L. Cohen,  and S. G. Louie, Energy Gaps in Graphene Nanoribbons, Phys. Rev. Lett. 97, 216803 (2006a).
8. L. Pisani, J. A. Chan, B. Montanari,  and N. M. Harrison, Electronic structure and magnetic properties of graphitic ribbons, Phys. Rev. B 75, 064418 (2007).
9. O. V. Yazyev and M. I. Katsnelson, Magnetic Correlations at Graphene Edges: Basis for Novel Spintronics Devices, Phys. Rev. Lett. 100, 047209 (2008).
10. C. Tao, L. Jiao, O. V. Yazyev, Y.-C. Chen, J. Feng, X. Zhang, R. B. Capaz, J. M. Tour, A. Zettl, S. G. Louie, H. Dai,  and M. F. Crommie, Spatially resolving edge states of chiral graphene nanoribbons, Nat. Phys. 7, 1745 (2011).
11. M. Sepioni, R. R. Nair, S. Rablen, J. Narayanan, F. Tuna, R. Winpenny, A. K. Geim,  and I. V. Grigorieva, Limits on Intrinsic Magnetism in Graphene, Phys. Rev. Lett. 105, 207205 (2010).
12. K. Wakabayashi, Y. Takane, M. Yamamoto,  and M. Sigrist, Electronic transport properties of graphene nanoribbons, New J. Phys. 11, 095016 (2009).
13. O. V. Yazyev, Emergence of magnetism in graphene materials and nanostructures, Rep. Prog. Phys. 73, 056501 (2010).
14. P. Ruffieux, S. Wang, B. Yang, C. Sánchez-Sánchez, J. Liu, T. Dienel, L. Talirz, P. Shinde, C. A. Pignedoli, D. Passerone, T. Dumslaff, X. Feng, K. Müllen,  and R. Fasel, On-surface synthesis of graphene nanoribbons with zigzag edge topology, Nature (London) 531, 489 (2016).
15. Y.-W. Son, M. L. Cohen,  and S. G. Louie, Half-metallic graphene nanoribbons, Nature (London) 444, 347 (2006b).
16. A. Rycerz, J. Tworzydlo,  and C. W. J. Beenakker, Valley filter and valley valve in graphene, Nat. Phys. 3, 172 (2007).
17. T. Hikihara, X. Hu, H.-H. Lin,  and C.-Y. Mou, Ground-state properties of nanographite systems with zigzag edges, Phys. Rev. B 68, 035432 (2003).
18. I. Hagymási and O. Legeza, Entanglement, excitations, and correlation effects in narrow zigzag graphene nanoribbons, Phys. Rev. B 94, 165147 (2016).
19. D. J. Luitz, F. F. Assaad,  and M. J. Schmidt, Exact diagonalization study of the tunable edge magnetism in graphene, Phys. Rev. B 83, 195432 (2011).
20. H. Feldner, Z. Y. Meng, T. C. Lang, F. F. Assaad, S. Wessel,  and A. Honecker, Dynamical Signatures of Edge-State Magnetism on Graphene Nanoribbons, Phys. Rev. Lett. 106, 226401 (2011).
21. E. H. Lieb, Two Theorems on the Hubbard Model, Phys. Rev. Lett. 62, 1201 (1989).
22. H. Karimi and I. Affleck, Towards a rigorous proof of magnetism on the edges of graphene nanoribbons, Phys. Rev. B 86, 115446 (2012).
23. V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea,  and A. H. Castro Neto, Electron-Electron Interactions in Graphene: Current Status and Perspectives, Rev. Mod. Phys. 84, 1067 (2012).
24. J. González, F. Guinea,  and M. Vozmediano, Non-Fermi liquid behavior of electrons in the half-filled honeycomb lattice (A renormalization group approach), Nucl. Phys. B 424, 595 (1994).
25. J. González, F. Guinea,  and M. A. H. Vozmediano, Marginal-Fermi-liquid behavior from two-dimensional Coulomb interaction, Phys. Rev. B 59, R2474 (1999).
26. E. Barnes, E. H. Hwang, R. E. Throckmorton,  and S. Das Sarma, Effective field theory, three-loop perturbative expansion, and their experimental implications in graphene many-body effects, Phys. Rev. B 89, 235431 (2014).
27. J. Hofmann, E. Barnes,  and S. Das Sarma, Why Does Graphene Behave as a Weakly Interacting System? Phys. Rev. Lett. 113, 105502 (2014).
28. C. Bauer, A. Rückriegel, A. Sharma,  and P. Kopietz, Nonperturbative renormalization group calculation of quasiparticle velocity and dielectric function of graphene, Phys. Rev. B 92, 121409 (2015).
29. A. Sharma and P. Kopietz, Multilogarithmic velocity renormalization in graphene, Phys. Rev. B 93, 235425 (2016).
30. D. C. Elias, R. V. Gorbachev, A. S. Mayorov, S. V. Morozov, A. A. Zhukov, P. Blake, L. A. Ponomarenko, I. V. Grigorieva, K. S. Novoselov, F. Guinea,  and A. K. Geim, Dirac cones reshaped by interaction effects in suspended graphene, Nat. Phys. 7, 701 (2011).
31. Y. Araki and G. W. Semenoff, Spin versus charge-density-wave order in graphenelike systems, Phys. Rev. B 86, 121402 (2012).
32. W. Wu and A.-M. S. Tremblay, Phase diagram and Fermi liquid properties of the extended Hubbard model on the honeycomb lattice, Phys. Rev. B 89, 205128 (2014).
33. D. Smith and L. von Smekal, Monte Carlo simulation of the tight-binding model of graphene with partially screened Coulomb interactions, Phys. Rev. B 89, 195429 (2014).
34. M. Golor and S. Wessel, Nonlocal density interactions in auxiliary-field quantum Monte Carlo simulations: Application to the square lattice bilayer and honeycomb lattice, Phys. Rev. B 92, 195154 (2015).
35. H.-K. Tang, E. Laksono, J. N. B. Rodrigues, P. Sengupta, F. F. Assaad,  and S. Adam, Interaction-Driven Metal-Insulator Transition in Strained Graphene, Phys. Rev. Lett. 115, 186602 (2015).
36. J. P. L. Faye, P. Sahebsara,  and D. Sénéchal, Chiral triplet superconductivity on the graphene lattice, Phys. Rev. B 92, 085121 (2015).
37. L. Classen, I. F. Herbut, L. Janssen,  and M. M. Scherer, Competition of density waves and quantum multicritical behavior in Dirac materials from functional renormalization, Phys. Rev. B 93, 125119 (2016).
38. Y. Volpez, D. D. Scherer,  and M. M. Scherer, Electronic instabilities of the extended Hubbard model on the honeycomb lattice from functional renormalization, Phys. Rev. B 94, 165107 (2016).
39. A. Katanin, Effect of vertex corrections on the possibility of chiral symmetry breaking induced by long-range Coulomb repulsion in graphene, Phys. Rev. B 93, 035132 (2016).
40. D. S. de la Peña, J. Lichtenstein,  and C. Honerkamp, Competing electronic instabilities of extended Hubbard models on the honeycomb lattice: A functional renormalization group calculation with high-wave-vector resolution, Phys. Rev. B 95, 085143 (2017).
41. I. S. Tupitsyn and N. V. Prokof’ev, Stability of Dirac Liquids with Strong Coulomb Interaction, Phys. Rev. Lett. 118, 026403 (2017).
42. H.-K. Tang, J. N. Leaw, J. N. B. Rodrigues, I. F. Herbut, P. Sengupta, F. F. Assaad,  and S. Adam, The role of electron-electron interactions in two dimensional Dirac fermions, Preprint  (2017).
43. A. Yamashiro, Y. Shimoi, K. Harigaya,  and K. Wakabayashi, Spin- and charge-polarized states in nanographene ribbons with zigzag edges, Phys. Rev. B 68, 193410 (2003).
44. B. Wunsch, T. Stauber, F. Sols,  and F. Guinea, Interactions and Magnetism in Graphene Boundary States, Phys. Rev. Lett. 101, 036803 (2008).
45. J. Jung, Nonlocal exchange effects in zigzag-edge magnetism of neutral graphene nanoribbons, Phys. Rev. B 83, 165415 (2011).
46. Z. Shi and I. Affleck, Effect of long-range interaction on graphene edge magnetism, Phys. Rev. B 95, 195420 (2017).
47. M. Hohenadler, F. Parisen Toldin, I. F. Herbut,  and F. F. Assaad, Phase diagram of the Kane-Mele-Coulomb model, Phys. Rev. B 90, 085146 (2014).
48. C. Hwang, D. A. Siegel, S.-K. Mo, W. Regan, A. Ismach, Y. Zhang, A. Zettl,  and A. Lanzara, Fermi velocity engineering in graphene by substrate modification, Sci. Rep. 2, 590 (2012).
49. F. Assaad and H. Evertz, in Computational Many-Particle Physics, Lecture Notes in Physics, Vol. 739, edited by H. Fehske, R. Schneider,  and A. Weiße (Springer, Berlin Heidelberg, 2008) pp. 277–356.
50. R. Brower, C. Rebbi,  and D. Schaich, Hybrid Monte Carlo simulation on the graphene hexagonal lattice, PoS(Lattice 2011)056 (arXiv:1204.5424) .
51. M. V. Ulybyshev, P. V. Buividovich, M. I. Katsnelson,  and M. I. Polikarpov, Monte Carlo Study of the Semimetal-Insulator Phase Transition in Monolayer Graphene with a Realistic Interelectron Interaction Potential, Phys. Rev. Lett. 111, 056801 (2013).
52. Y. Zhang and J. Callaway, Extended Hubbard model in two dimensions, Phys. Rev. B 39, 9397 (1989).
53. Y. Ohta, K. Tsutsui, W. Koshibae,  and S. Maekawa, Exact-diagonalization study of the Hubbard model with nearest-neighbor repulsion, Phys. Rev. B 50, 13594 (1994).
54. M. Aichhorn, H. G. Evertz, W. von der Linden,  and M. Potthoff, Charge ordering in extended Hubbard models: Variational cluster approach, Phys. Rev. B 70, 235107 (2004).
55. B. Davoudi and A.-M. S. Tremblay, Non-perturbative treatment of charge and spin fluctuations in the two-dimensional extended Hubbard model: Extended two-particle self-consistent approach, Phys. Rev. B 76, 085115 (2007).
56. E. G. C. P. van Loon, A. I. Lichtenstein, M. I. Katsnelson, O. Parcollet,  and H. Hafermann, Beyond extended dynamical mean-field theory: Dual boson approach to the two-dimensional extended Hubbard model, Phys. Rev. B 90, 235135 (2014).
57. L. Huang, T. Ayral, S. Biermann,  and P. Werner, Extended dynamical mean-field study of the Hubbard model with long-range interactions, Phys. Rev. B 90, 195114 (2014).
58. H. Terletska, T. Chen,  and E. Gull, Charge ordering and correlation effects in the extended Hubbard model, Phys. Rev. B 95, 115149 (2017).
59. C. Koop and M. J. Schmidt, Effective spin theories for edge magnetism in graphene zigzag ribbons, Phys. Rev. B 92, 125416 (2015).
60. S. Sorella, Y. Otsuka,  and S. Yunoki, Absence of a Spin Liquid Phase in the Hubbard Model on the Honeycomb Lattice, Sci. Rep. 2, 992 (2012).
61. F. F. Assaad and I. F. Herbut, Pinning the Order: The Nature of Quantum Criticality in the Hubbard Model on Honeycomb Lattice, Phys. Rev. X 3, 031010 (2013).
62. M. Schüler, M. Rösner, T. O. Wehling, A. I. Lichtenstein,  and M. I. Katsnelson, Optimal Hubbard Models for Materials with Nonlocal Coulomb Interactions: Graphene, Silicene, and Benzene, Phys. Rev. Lett. 111, 036601 (2013).
63. S. R. White, R. M. Noack,  and D. J. Scalapino, Resonating Valence Bond Theory of Coupled Heisenberg Chains, Phys. Rev. Lett. 73, 886 (1994).
64. E. Dagotto and T. M. Rice, Surprises on the Way from One- to Two-Dimensional Quantum Magnets: The Ladder Materials, Science 271, 618 (1996).
65. K. S. D. Beach, Identifying the maximum entropy method as a special limit of stochastic analytic continuation, arXiv e-prints  (2004), cond-mat/0403055 .
66. A. Schneider and H. Wierstorf, Gnuplot-Colorbrewer: Colorbrewer Color Schemes For Gnuplot, Zenodo  (2014), 10.5281/zenodo.10282.
67. K. Wakabayashi, M. Sigrist,  and M. Fujita, Spin Wave Mode of Edge-Localized Magnetic States in Nanographite Zigzag Ribbons, J. Phys. Soc. Jpn. 67, 2089 (1998).
68. F. J. Culchac, A. Latgé,  and A. T. Costa, Spin waves in zigzag graphene nanoribbons and the stability of edge ferromagnetism, New J. Phys. 13, 033028 (2011).
69. M. Raczkowski and A. M. Oleś, Spin stiffness and quantum fluctuations in C-type and A-type antiferromagnets, Phys. Rev. B 66, 094431 (2002).
70. M. P. López-Sancho and L. Brey, Charged Topological Solitons in Zigzag Graphene Nanoribbons, arXiv e-prints  (2017), 1705.07753 .
71. Jülich Supercomputing Centre, JURECA: General-purpose supercomputer at Jülich Supercomputing Centre, Journal of large-scale research facilities 2 (2016), 10.17815/jlsrf-2-121.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters