Effect of long-range interaction on graphene edge magnetism

TL;DR

This paper investigates how long-range Coulomb interactions affect edge magnetism in graphene, finding that ferromagnetism persists even with unscreened Coulomb forces, which is crucial for experimental realization.

Contribution

It introduces a model considering non-local Coulomb interactions on graphene edges and demonstrates that ferromagnetism remains stable under these conditions.

Findings

Ferromagnetism survives with unscreened Coulomb interactions

Long-range interactions do not destroy edge magnetism in graphene

The model includes electron, hole, and exciton excitations

Abstract

It has been proposed that interactions lead to ferromagnetism on a zigzag edge of a graphene sheet. While not yet directly studied experimentally, dramatically improving techniques for making and studying clean zigzag edges may soon make this possible. So far, most theoretical investigations of this claim have been based on mean field theories or more exact calculations using the Hubbard model. But long-range Coulomb interactions are unscreened in graphene so it is important to consider their effects. We study rather general non-local interactions, including of Coulomb $1/r$ form, using the technique of projection to a strongly interacting edge Hamiltonian, valid at first order in the interactions. The ground states as well as electron/hole and exciton excitations are studied in this model. Our results indicate that ferromagnetism survives with unscreened Coulomb interactions.

Equations

$$e_k^{†}=\frac{1}{\sqrt{2\pi }}\sum _{n\ge 0,m}{e}^{ik\frac{m}{2}}{g}_n\left(k\right)c_{m,n,A}^{†}\text{,}$$

$$g_n\left(k\right)\equiv \theta \left(k-\frac{2\pi }{3}\right)\theta \left(\frac{4\pi }{3}-k\right)\sqrt{1-4{\cos }^2\frac{k}{2}}{\left(-2\cos \frac{k}{2}\right)}^n$$

$$b_{k,k_y,s}^{†}$$

$$E_s\left(k,k_y\right)=st\sqrt{{\left(2\cos \frac{k}{2}\right)}^2+1+2\left(2\cos \frac{k}{2}\right)\cos k_y}$$

$$H_{int}=\frac{1}{2}\sum _{n,m}\sum _{{\delta }_m,{\delta }_n}{U}_{\left({\delta }_m,{\delta }_n\right)}\left(\sum _{\sigma =\pm }{c}_{m,n,\sigma }^{†}{c}_{m,n,\sigma }-1\right)\left(\sum _{{\sigma }^{\prime }=\pm }{c}_{m+{\delta }_m,n+{\delta }_n,{\sigma }^{\prime }}{c}_{m+{\delta }_m,n+{\delta }_n,{\sigma }^{\prime }}-1\right)\text{.}$$

$$U_{\left({\delta }_m,{\delta }_n\right)}=U_0\frac{d}{\sqrt{d^2+{\left|\vec{\delta }\right|}^2}}\text{,}$$

$$⟨b_{k,k_y,-,\sigma }^{†}{b}_{k^{\prime },k_y^{\prime },-,{\sigma }^{\prime }}⟩=⟨b_{k,k_y,+,\sigma }{b}_{k^{\prime },k_y^{\prime },+,{\sigma }^{\prime }}^{†}⟩={\delta }_{\sigma {\sigma }^{\prime }}\delta \left(k-k^{\prime }\right)\delta \left(k_y-k_y^{\prime }\right)\text{.}$$

$$H_{\text{int}}=\frac{1}{2}\sum _n\sum _{{\delta }_m,{\delta }_n}{U}_{\left({\delta }_m,{\delta }_n\right)}{\int }_{-\frac{2\pi }{3}}^{\frac{2\pi }{3}}\frac{dq}{2\pi }{e}^{iq\frac{{\delta }_m}{2}}{O}_{n+{\delta }_n}^{†}\left(q\right)O_n\left(q\right)\text{,}$$

$$O_n\left(q\right)\equiv \int dk{g}_n\left(k+q\right)g_n\left(k\right)\left[\sum _{\sigma =\pm }{e}_{k+q,\sigma }^{†}{e}_{k,\sigma }-\delta \left(q\right)\right]\text{.}$$

$$H=H_{\text{int}}+H_{\Delta }\text{, }{H}_{\Delta }=\Delta \sum _{\sigma =\pm }{\int }_{\frac{2\pi }{3}}^{\frac{4\pi }{3}}dk\left(2\cos k+1\right)e_{k,\sigma }^{†}{e}_{k,\sigma }\text{,}$$

$$H_{\text{int, Hubbard}}=\frac{1}{2}U\sum _{n=0}^{\infty }{\int }_{-\frac{2\pi }{3}}^{\frac{2\pi }{3}}\frac{dq}{2\pi }{O}_n^{†}\left(q\right)O_n\left(q\right)\text{;}$$

$$O_n\left(q\right)={\int }_0^{\pi }\frac{dK}{\pi }{O}_K\left(q\right)\sin K\left(n+1\right)\text{,}$$

$$H_{\text{int}}=\frac{1}{2}{\int }_{-\frac{2\pi }{3}}^{\frac{2\pi }{3}}\frac{dq}{2\pi }{\int }_0^{\pi }\frac{dK}{2\pi }\tilde{U}\left(k,q\right)O_K^{†}\left(q\right)O_K\left(q\right)\text{,}$$

$$\tilde{U}\left(K,q\right)\equiv \sum _{{\delta }_m,{\delta }_n}{U}_{\left({\delta }_m,{\delta }_n\right)}\cos \left(K{\delta }_n\right)\cos \frac{q{\delta }_m}{2}\text{.}$$

$$U_{\left(0,0\right)}\equiv U\text{, }{U}_{\left(\pm 2,0\right)}\equiv U_{2\parallel }\text{, }{U}_{\left(\pm 1,1\right)}=U_{\left(\pm 1,-1\right)}\equiv U_{2\angle }\text{,}$$

$$\tilde{U}\left(K,q\right)=U+2U_{2\parallel }\cos q+4U_{2\angle }\cos K\cos \frac{q}{2}\text{;}$$

$$\forall q\in \left[-\frac{2\pi }{3},\frac{2\pi }{3}\right]\text{, }U+2U_{2\parallel }\cos q\ge 4U_{2\angle }\cos \frac{q}{2}\text{.}$$

$${\int }_{\frac{2\pi }{3}}^{\frac{4\pi }{3}}dkf\left(k\right)e_{k,\uparrow }{e}_{k,\downarrow }^{†}|\text{FM}\uparrow ⟩\text{.}$$

$$E\left[f\right]=U{C}_0\left[f\right]+U_{2\parallel }{C}_{2\parallel }\left[f\right]+U_{2\angle }{C}_{2\angle }\left[f\right]\text{.}$$

$$U_{\left(0,0\right)}\equiv U\text{, }{U}_{\left(0,\pm 2\right)}=U_{\left(\pm 3,1\right)}=U_{\left(\pm 3,-1\right)}\equiv U_5\text{,}$$

$$\tilde{U}\left(K,q\right)=U+2U_5\left(\cos 2K+2\cos K\cos \frac{3q}{2}\right)\text{,}$$

$$H_{\text{int}}$$

$$\Gamma \left(k,k^{\prime },q\right)$$

$${ϵ}_{p/h}\left(k\right)=\frac{1}{2}{\int }_{\frac{2\pi }{3}}^{\frac{4\pi }{3}}\frac{d{k}^{\prime }}{2\pi }\Gamma \left(k,k^{\prime },k^{\prime }-k\right)\pm \Delta \left(2\cos k+1\right)\text{.}$$

$${ϵ}_{p/h}\left(k\right)\approx \left(v\mp \sqrt{3}\Delta \right)\left(k-\frac{2\pi }{3}\right)\text{,}$$

$$v\equiv \frac{\sqrt{3}}{2}\sum _{{\delta }_m,{\delta }_n}{U}_{\left({\delta }_m,{\delta }_n\right)}{\int }_{\frac{2\pi }{3}}^{\frac{4\pi }{3}}\frac{d{k}^{\prime }}{2\pi }{\left(2\cos \frac{k^{\prime }}{2}\right)}^{|{\delta }_n|}\cos \left(k^{\prime }-\frac{2\pi }{3}\right)\frac{{\delta }_m}{2}\text{.}$$

$$v=\frac{\sqrt{3}}{6}U+\frac{3}{4\pi }{U}_{2\parallel }+\left(\frac{1}{\sqrt{3}}-\frac{3}{2\pi }\right)U_{2\angle }\text{.}$$

$${\int }_{\frac{2\pi }{3}}^{\frac{4\pi }{3}}\frac{d{k}^{\prime }}{2\pi }{\left(2\cos \frac{k^{\prime }}{2}\right)}^{|{\delta }_n|}\cos \left(k^{\prime }-\frac{2\pi }{3}\right)\frac{{\delta }_m}{2}$$

$$v$$

$${ϵ}_{p/h}\left(k\right)\approx \frac{2U_{0}d}{\pi }\left(k-\frac{2\pi }{3}\right)\ln \frac{\Lambda }{k-\frac{2\pi }{3}}\text{,}$$

Full text

I Introduction

Non-interacting graphene nanoribbons with zigzag edges are famous for hosting a nearly flat band of edge states.Fujita et al. (1996); Wakabayashi et al. (2010) In the presence of electron-electron interaction, the existence of edge magnetic orderYazyev (2010) has been predicted by a multitude of theoretical work using both analyticalFujita et al. (1996); Yamashiro et al. (2003); Jung et al. (2009); *PhysRevB.79.235433; Rhim and Moon (2009); Schmidt and Loss (2010); Jung (2011); Culchac et al. (2011); Schmidt (2012); Karimi and Affleck (2012) and numericalLuitz et al. (2011); Schmidt et al. (2013); Wunsch et al. (2008); Golor et al. (2013a, b, 2014); *PhysRevB.87.245431; Hagymási and Legeza (2016); Son et al. (2006); Pisani et al. (2007); Huang et al. (2013) techniques. The consensus emerging from these work is that edge states localized at the same edge are coupled ferromagnetically to form superspins, which then couple antiferromagnetically between edges. In addition to ground state properties, low-energy magnetic excitations in graphene nanoribbons have also attracted much theoretical attention.Wakabayashi et al. (1998); Yazyev and Katsnelson (2008); Rhim and Moon (2009); Culchac et al. (2011) A relatively large spin correlation length up to the order of micrometers has been found for a single zigzag edge; this is attributed to the large spin stiffness in this system, and boosts confidence in potential spintronics applications of graphene edge magnetism.Zhang et al. (2010) Although conclusive experimental evidence for edge magnetism is still lacking due to limited control over edge orientation, there has been significant progress in recent years towards the synthesis and characterization of zigzag edges.Magda et al. (2014); Ruffieux et al. (2015); Makarova et al. (2015)

A large number of theoretical studies on graphene edge magnetism represent the interaction by an on-site Hubbard term for simplicity. For the Hubbard model on a honeycomb lattice, arguments in support of edge magnetismKarimi and Affleck (2012) can be constructed based on Lieb’s theorem.Lieb (1989) The Coulomb interaction in pristine graphene on a non-metallic substrate is, nevertheless, poorly screened due to a vanishing density of states at the Dirac points.DiVincenzo and Mele (1984); Kotov et al. (2012) The influence of non-local components of the interaction has been investigated both in bulk grapheneWu and Tremblay (2014); Herbut (2006); Honerkamp (2008); Schüler et al. (2013) and in restricted geometries.Yamashiro et al. (2003); Jung (2011); Chacko et al. (2014); Zhu and Wang (2006); Wunsch et al. (2008) (By “non-local” we mean having a longer range than on-site.) However, many studies on graphene nanoribbons with non-local interactions have adopted a mean-field treatment, neglecting fluctuations whose role is especially important in low dimensions.Mermin and Wagner (1966) Exact diagonalization has been employed in other studies; despite the light it sheds on the nature of the ground states, correlations in manageably small systems are usually enhanced compared to the thermodynamic limit.

In the present work, we study the effect of long-range interactions on graphene edge ferromagnetism, in the limit of weak interactions but beyond the mean-field level. Focusing on a semi-infinite graphene sheet with a single zigzag edge, we find the effective Hamiltonian by projecting the interaction into the Hilbert space of edge states; we then propose a sufficient condition for the maximum spin ferromagnetic multiplet to be the half-filling ground states. Using exact diagonalization, we discuss the possible ground states for interactions in violation of this condition. The long-range Coulomb interaction is shown to satisfy the sufficient condition upon extrapolation to the limit of infinite long distance cutoff. We also examine the simplest low-energy excitations of the ferromagnetic ground states on a single edge. For short-range interactions, single-particle excitations and single-hole excitations have linear spectra $\propto v\delta k$ where $\left|\delta k\right|\ll 1$ is the distance from either Dirac point, with a slope $v$ controlled by the interaction strength. Spin-1 excitons have a small-momentum dispersion that is proportional to $vQ^{2}\ln Q$. For the long-range Coulomb interaction, $v\rightarrow\infty$, and the dispersion of single-particle or single-hole excitations near the Dirac points scales as $\delta k\ln\delta k$. Finally, for both short-range and Coulomb interactions, a sufficiently large particle-hole symmetry breaking term in the Hamiltonian can destabilize the ferromagnetic ground state.

II Model

We study a semi-infinite graphene sheet on the $xy$ plane, modeled by a honeycomb lattice which is terminated by an infinite zigzag edge (see Fig. 1). All carbon atoms reside in the half plane $y\geq 0$, and the outermost atoms on the zigzag edge (which belong to the $A$ hexagonal sublattice) lie on the $x$ axis. In units of the Bravais lattice constant $a=2.46$Å, it is convenient to represent the position of carbon atoms by $\vec{r}_{\left(m,n\right)}=\left(m/2\right)\hat{x}+\left(\sqrt{3}n/2\right)\hat{y}$ where $n\geq 0$. While $m$ is always an integer, note that $n$ is an integer only on the $A$ sublattice: for $A$ atoms $n$ and $m$ are both even or both odd, while for $B$ atoms $n+2/3$ and $m$ are both even or both odd.

The zero modes associated with the zigzag edge are given byFujita et al. (1996); Wakabayashi et al. (2010)

[TABLE]

where $k$ is the crystal momentum along the edge direction,

[TABLE]

describes the decay of the wave function into the bulk, and the $c$ operators obey the usual anticommutation relations $\left\{c_{m,n,A},c_{m^{\prime},n^{\prime},A}^{{\dagger}}\right\}=\left\{c_{m,n,B},c_{m^{\prime},n^{\prime},B}^{{\dagger}}\right\}=\delta_{mm^{\prime}}\delta_{nn^{\prime}}$, $\left\{c_{m,n,A},c_{m^{\prime},n^{\prime},B}^{{\dagger}}\right\}=0$. (We have temporarily suppressed the spin index.) These edge states exist only for $2\pi/3<k<4\pi/3$, i.e. in $1/3$ of the 1D Brillouin zone $0\leq k<2\pi$. The wave function is non-zero only on the $A$ sublattice, and is localized near the zigzag edge. The localization length $\xi_{k}=-\left[\ln\left|2\cos\left(k/2\right)\right|\right]^{-1}$ vanishes at $k=\pi$, and diverges near the Dirac points $k=2\pi/3$ and $k=4\pi/3$.

In addition to the edge states, we also have bulk states which are labeled by $k$, $k_{y}$ and $s$:

[TABLE]

Here the bulk dispersion relation is

[TABLE]

with nearest neighbor hopping strength $t$; $k_{y}$ is the crystal momentum perpendicular to the edge, $0\leq k_{y}\leq\pi$, and $s=\pm$ is a subband index. Near the Dirac points, where $\left(k,k_{y}\right)=\left(2\pi/3,\pi\right)+\left(\delta k,\delta k_{y}\right)$ or $\left(k,k_{y}\right)=\left(4\pi/3,0\right)+\left(\delta k,\delta k_{y}\right)$, $E_{s}\left(k,k_{y}\right)$ takes a Lorentz invariant form $E_{s}\left(k,k_{y}\right)=st\sqrt{\left(\delta k_{y}\right)^{2}+\left(3/4\right)\left(\delta k\right)^{2}}$. In this non-interacting model, at zero temperature and half-filling, the $s=-$ subband is completely filled and the $s=+$ subband is completely empty. While the edge states are half-filled, for the semi-infinite sheet we cannot ascertain which half is filled at this point, unless other ingredients–such as next-nearest-neighbor hopping, edge potential and electron-electron interaction–are present.

We now introduce a weak repulsive electron-electron interaction. The following extended Hubbard model manifestly respects $SU\left(2\right)$ spin symmetry, and also particle-hole symmetry at half-filling:

[TABLE]

$\left(\delta_{m},\delta_{n}\right)$ runs over all vectors $\vec{\delta}=\left(\delta_{m}/2\right)\hat{x}+\left(\sqrt{3}\delta_{n}/2\right)\hat{y}$ pointing from one lattice site to another; for instance, $U_{\left(0,0\right)}$ stands for the strength of the on-site Hubbard interaction, $U_{\left(0,2/3\right)}$ is the interaction between nearest neighbor sites (belonging to different sublattices) in the $y$ direction, $U_{\left(1,1/3\right)}$ is the interaction between nearest neighbor sites at $\pi/6$ angle with the $x$ direction, and $U_{\left(2,0\right)}$ is the interaction between next nearest neighbors (belonging to the same sublattice) in the $x$ direction. The sum over $n$ and $\delta_{n}$ is such that both $n\geq 0$ and $n+\delta_{n}\geq 0$. To lighten notations, we have suppressed the sublattice indices $A$ and $B$ in this expression, because they are uniquely determined by the position indices $\left(m,n\right)$ and $\left(m+\delta_{m},n+\delta_{n}\right)$.

In general $U_{\left(\delta_{m},\delta_{n}\right)}=U_{\left(-\delta_{m},-\delta_{n}\right)}$, but apart from this constraint $U$ can be an arbitrary function of $\delta_{m}$ and $\delta_{n}$. Nevertheless we further assume that $U$ obeys parity symmetry, $U_{\left(\delta_{m},\delta_{n}\right)}=U_{\left(-\delta_{m},\delta_{n}\right)}$. For the Hubbard model, $U_{\left(\delta_{m},\delta_{n}\right)}$ vanishes unless $\delta_{m}=\delta_{n}=0$. On the other hand, for the unscreened Coulomb interaction, $U_{\left(\delta_{m},\delta_{n}\right)}$ is inversely proportional to distance at large distances,Egger and Gogolin (1997); Zarea and Sandler (2007); *PhysRevLett.101.196804; Jung (2011)

[TABLE]

where $U_{0}$ is the on-site interaction, and the half-nearest-neighbor distance $d=1/\left(2\sqrt{3}\right)$ accounts for the finite spread of the carbon $\pi$ orbitals.

Assuming $U_{\left(\delta_{m},\delta_{n}\right)}\ll t$, we expect that the low-energy degrees of freedom are composed of the edge states $e_{k}$ with $2\pi/3<k<4\pi/3$, and the bulk states in the vicinity of the two Dirac points.Wunsch et al. (2008); Schmidt et al. (2013) As a first approximation at $O\left(U\right)$, we neglect the dynamics of the bulk states completely; they are assumed to be half-filled and not spin-polarized as in the non-interacting case.Schmidt and Loss (2010); Karimi and Affleck (2012) This approximation allows the projection of the interaction onto the Hilbert space of the edge states. More concretely, we invert Eqs. (1) and (3) to express the $c$ operators in terms of $e$ and $b$, then take the expectation values for pairs of $b$ operators using

[TABLE]

After some algebra, we find

[TABLE]

where the sum over $\left(\delta_{m},\delta_{n}\right)$ is now limited to vectors on one of the sublattices; recalling that edge states only exist on the $A$ sublattice, $\delta_{m}$ and $\delta_{n}$ are now both even or both odd. Again $n\geq 0$ and $n+\delta_{n}\geq 0$. $O_{n}\left(q\right)$ is bilinear in $e$,

[TABLE]

$q$ measures the momentum difference between two edge states, so the operator $O_{n}\left(q\right)$ is nontrivial only when $\left|q\right|<2\pi/3$. Note that $O_{n}\left(q\right)$ annihilates all members of the fully polarized ferromagnetic multiplet at half-filling for any $n$ and $q$, which means the ferromagnetic multiplet states are always eigenstates of $H_{\text{int}}$ with zero energy.

Due to the constraint on the $\left(\delta_{m},\delta_{n}\right)$ summation, many terms in the interaction (most notably the nearest-neighbor interaction) do not enter the projected effective Hamiltonian in the edge state subspace, Eq. (8). Although the authors of Ref. Yamashiro et al., 2003 predict a charge-polarized ground state when the nearest-neighbor interaction prevails over the on-site interaction, our picture is consistent with their weak interaction limit, where the charge-polarized state always has a higher energy and the nearest-neighbor interaction is unimportant.

Just as Eq. (5), Eq. (8) manifestly respects $SU\left(2\right)$ symmetry and particle-hole symmetry at half-filling. In particular, the particle-hole transformation $c_{m,n,\sigma}\rightarrow c_{m,n,\sigma}^{{\dagger}}$ corresponds to $e_{k,\sigma}\rightarrow e_{2\pi-k,\sigma}^{{\dagger}}$ and $O_{n}\left(q\right)\rightarrow-O_{n}\left(q\right)$ in the edge state subspace. (The form $e_{k,\sigma}\rightarrow e_{k,\sigma}^{{\dagger}}$ previously suggested in the Hubbard modelKarimi and Affleck (2012) is the combination of a particle-hole transformation and a parity transformation.) The particle-hole symmetry is broken by either a weak next-nearest neighbor hopping $\left|t_{2}\right|\ll t$ in the bulk, or a weak potential localized at the edge $\left|V_{e}\right|\ll t$; the latter can arise, for example, at a graphene-graphane interface.Schmidt and Loss (2010) When $\Delta=t_{2}-V_{e}\neq 0$, a dispersion develops for the edge states:

[TABLE]

assuming the Fermi energy is fixed at the new Dirac point $\epsilon_{F}=3t_{2}$.Schmidt and Loss (2010); Karimi and Affleck (2012)

In the remainer of this paper we analyze the edge state Hamiltonian given by Eq. (10) at half-filling.

III Ground state

We first study the ground state of the particle-hole symmetric Hamiltonian Eq. (8), keeping $\Delta=0$.

For the projected Hubbard model, it has been proven in Ref. Karimi and Affleck, 2012 that the fully polarized ferromagnetic multiplet states with maximum total spin are the unique ground states. In the Hubbard case, Eq. (8) becomes

[TABLE]

It is obvious that $H_{\text{int, Hubbard}}$ is positive semi-definite. Since the ferromagnetic multiplet states are always zero energy eigenstates, they must belong to the ground state manifold of $H_{\text{int, Hubbard}}$. Furthermore, it is also possible to show that they are the only states annihilated by $O_{n}\left(q\right)$ for any $n$ and $q$, and therefore the unique ground states of $H_{\text{int, Hubbard}}$.Karimi and Affleck (2012) We emphasize again that the proof rests on the positive semi-definiteness of the Hamiltonian.

Let us explore the extent to which the proof outlined above can be generalized in our extended Hubbard model. In analogy to a semi-infinite tight-binding chain, through the following transformation

[TABLE]

the generic interaction Hamiltonian Eq. (8) can be formally diagonalized:

[TABLE]

where

[TABLE]

The spectrum of $\tilde{U}\left(K,q\right)$ does not give the spectrum of the interacting problem because $O_{K}\left(q\right)$ does not obey simple commutation relations. Nevertheless, if $\tilde{U}\left(K,q\right)$ is positive semi-definite for $0\leq K\leq\pi$ and $-2\pi/3\leq q\leq 2\pi/3$, we can borrow the arguments from the case of the Hubbard model, and show that the ferromagnetic multiplet states are the unique ground states of Eq. (8) at half filling. [That a state is annihilated by all $O_{n}\left(q\right)$ is equivalent to it being annihilated by all $O_{K}\left(q\right)$.] The positive semi-definiteness of $\tilde{U}\left(K,q\right)$ is thus a sufficient condition for ferromagnetic ground states.

As a simple example, we consider the model with only on-site and next-nearest-neighbor interactions:

[TABLE]

and $U_{\left(\delta_{m},\delta_{n}\right)}=0$ for other $\left(\delta_{m},\delta_{n}\right)$. (The nearest-neighbor interactions drop out, as remarked in Section II.) In the next-nearest-neighbor interaction we have introduced an anisotropy between the direction parallel to the edge ($U_{2\parallel}$) and the directions at an angle of $\pi/3$ with the edge ($U_{2\angle}$). While such anisotropy is not necessarily realistic, we shall see that $U_{2\parallel}$ and $U_{2\angle}$ have very different effects on edge magnetism.

For this model,

[TABLE]

as $\cos q/2>0$, the minimum of $\tilde{U}$ with respect to $K$ is obtained at $K=\pi$. The positive semi-definiteness condition of $\tilde{U}\left(K,q\right)$ is therefore equivalent to

[TABLE]

This is a sufficient condition for the ground states to be ferromagnetic in the model specified by Eq. (15). It requires that neither $U_{2\parallel}$ nor $U_{2\angle}$ should be greater than $U$. In particular, Eq. (17) becomes $U_{2\angle}\leq U/4$ when $U_{2\parallel}=0$, and $U_{2\parallel}\leq U$ when $U_{2\angle}=0$; in the isotropic case $U_{2\parallel}=U_{2\angle}\equiv U_{2}$,  Eq. (17) is reduced to $U_{2}\leq U/3$.

It is natural to wonder whether the fully polarized ferromagnetic multiplet remains the ground states of Eq. (15) at half filling when the sufficient condition Eq. (17) is violated. To answer the question we perform exact diagonalization on Eq. (15). Assuming a system size of $L$ unit cells along the edge, the number of different edge state momenta allowed is approximately$\ N=L/3$. It is convenient to take advantage of the good quantum numbers of the Hamiltonian, namely the $z$ component of the total spin $S_{z}$ and also the total momentum $Q$ along the edge direction.Wunsch et al. (2008) We measure $Q$ relative to the fully polarized state $\left|\text{FM}\uparrow\right\rangle$ where every edge state is singly occupied by a spin-up electron; for this state $S_{z}=N/2$ and $Q=0$.

In Fig. 2 we plot the ferromagnetic phase boundary for Eq. (15) on the $U_{2\parallel}$-$U_{2\angle}$ plane, obtained from exact diagonalization. For comparison we also show the region where the sufficient condition Eq. (17) is satisfied. In most of the parameter space, we find that the ground states at half filling are uniquely given by the $\left(N+1\right)$-fold degenerate ferromagnetic multiplet with $S_{z}=-N/2$, $-N/2+1$, …, $N/2$, and $Q=0$. In particular, the ground states are always ferromagnetic in the isotropic case $U_{2\parallel}=U_{2\angle}$. However, in the region above the phase boundary where $U_{2\angle}$ is relatively large compared to both $U$ and $U_{2\parallel}$, the ground states are not part of the ferromagnetic multiplet, but rather form a negative-energy manifold with a lower degeneracy and a lower total spin. For fixed $U$ and $U_{2\parallel}$, the degeneracy is reduced as $U_{2\angle}$ gradually increases, and eventually for sufficiently large $U_{2\angle}$ the ground state becomes a non-degenerate singlet state in the $S_{z}=0$ sector.

In Fig. 3, choosing a fixed $U_{2\parallel}/U$, we plot $E_{\text{GS}}\left(S_{z}\right)$ (the ground state energy in the sector labeled by $S_{z}$) as a function of $\left|S_{z}\right|$ for different $U_{2\angle}/U$ outside of the ferromagnetic regime. We observe that $E_{\text{GS}}\left(S_{z}\right)$ is a monotonically increasing function of $\left|S_{z}\right|$ in general, and becomes a strictly increasing function of $\left|S_{z}\right|$ if U_{2\angle}\is sufficiently large. This property of $E_{\text{GS}}\left(S_{z}\right)$ allows us to determine the ferromagnetic phase boundary in Fig. 2 by calculating $E_{\text{GS}}\left(S_{z}=N/2-1\right)$, which for a given $N$ is considerably less numerically intensive than $E_{\text{GS}}\left(S_{z}=0\right)$. Reasonably accurate estimates of the phase boundary can then be made through a variational calculation. We can characterize an arbitrary $Q=0$ state in the $S_{z}=N/2-1$ sector by

[TABLE]

The ferromagnetic state in this sector corresponds to $f\left(k\right)=1$, i.e. an equal-weighted superposition of all states where every edge state is singly occupied. The energy expectation value as a functional of $f$ is a linear combination of $U$, $U_{2\parallel}$ and $U_{2\angle}$:

[TABLE]

If $E\left[f\right]<0$, the ground states cannot be the ferromagnetic multiplet whose energy is always zero. For $f\left(k\right)\propto\sin k$, $C_{0}=0.100$, $C_{2\parallel}=0.0964$ and $C_{2\angle}=-0.0730$; for $f\left(k\right)\propto k-\pi$, $C_{0}=0.0946$, $C_{2\parallel}=0.0887$ and $C_{2\angle}=-0.0687$. For these two trial wave functions, the trajectories above which $E\left[f\right]<0$ are plotted in Fig. 2; both trajectories are very close to the ferromagnetic phase boundary obtained from exact diagonalization.

It should also be cautioned that anisotropy is not necessary to stabilize non-ferromagnetic ground states. For instance, we can also study an isotropic interaction consisting of an on-site term and six fifth-nearest-neighbor terms (or equivalently, next-nearest-neighbor terms on the same sublattice):

[TABLE]

and $U_{\left(\delta_{m},\delta_{n}\right)}=0$ for other $\left(\delta_{m},\delta_{n}\right)$. For this model

[TABLE]

so our sufficient condition for ferromagnetism becomes $U_{5}\leq U/3$. In a system with $N=720$, exact diagonalization shows that a non-ferromagnetic ground state appears when $U_{5}>80.48U$, i.e. when the non-local $U_{5}$ term is far stronger than the on-site interaction.

Our exact diagonalization results for both models indicate that while ferromagnetism is favored by the on-site interaction, it may be destabilized by sufficiently strong non-local interactions. This is in agreement with the findings of Ref. Schüler et al., 2013 that the effective on-site part of the interaction in bulk graphene is reduced by a weighted average of non-local interactions.

We now investigate whether the unscreened Coulomb interaction, Eq. (6), satisfies the sufficient condition for ferromagnetism. To this end, we introduce a long-distance cutoff $R$, and minimize $\tilde{U}\left(K,q\right)$ for the interaction that is given by Eq. (6) for $\left|\vec{\delta}\right|\leq R$ but vanishes for $\left|\vec{\delta}\right|>R$. In Fig. 4 we show $\tilde{U}_{\min}$, the minimum of $\tilde{U}\left(K,q\right)$ for $0\leq K\leq\pi$ and $-2\pi/3\leq q\leq 2\pi/3$, as a function of $R$ for $R\leq 500$. While $\tilde{U}_{\min}$ oscillates wildly, its lower envelope is an increasing function of $R$, and $\tilde{U}_{\min}$ does not go below $0.2U_{0}$ for $50\leq R\leq 500$. This strongly implies that $\tilde{U}_{\min}$ remains positive as $R\rightarrow\infty$, and provides evidence that the ferromagnetic multiplet states are the unique ground states for the unscreened Coulomb interaction.

A remark is in order about the short-distance cutoff $d=1/\left(2\sqrt{3}\right)$ in Eq. (6). If $d$ is treated as a tunable parameter of our model, then the observation that $\tilde{U}_{\min}\left(R\rightarrow\infty\right)>0$ is only valid when $d\lesssim 1$. If $d$ is close to $1$, $\tilde{U}_{\min}$ oscillates around zero even for $R$ up to $500$. Nevertheless, as shown in the next-nearest-neighbor model and the fifth-nearest-neighbor model, violation of the sufficient condition for ferromagnetism $\tilde{U}_{\min}\geq 0$ is not an indication of ground states being non-ferromagnetic. Indeed, we have verified in the $S_{z}=N/2-1$ sector that the ground states remain ferromagnetic for $R$ up to $20$ and $d$ up to $10$.

IV Low-energy excitations

In this section we discuss the low-energy single-particle, single-hole and particle-hole excitations of the ferromagnetic ground state, and also the effect of the particle-hole symmetry breaking term $\Delta$.

It is simplest to consider the excitations from the maximum $S_{z}$ state $\left|\text{FM}\uparrow\right\rangle$. We can rewrite the projected Hamiltonian of Eq. (10) in a form which explicitly annihilates $\left|\text{FM}\uparrow\right\rangle$:

[TABLE]

where the domains of integration are such that all edge states have momenta between $2\pi/3$ and $4\pi/3$, the interaction kernel is

[TABLE]

and the energy to create one single spin-down electron or one single spin-up hole is

[TABLE]

As noted in Refs. Luitz et al., 2011; Karimi and Affleck, 2012, the interaction $\Gamma\left(k,k^{\prime},q\right)$ is strongly momentum-dependent. For both Hubbard and Coulomb interactions, $\Gamma\left(k,k^{\prime},q\right)$ is positive so that spin-down electrons attract spin-down holes, which favors the formation of bound states between the two. The third term in Eq. (22) generally gives rise to interaction between edge states with the same spin orientation, although for the Hubbard model it vanishes due to an additional symmetry of the kernel, $\Gamma\left(k,k^{\prime},q\right)=\Gamma\left(k,k^{\prime},k^{\prime}-k-q\right)$.

IV.1 Single particle and single hole excitations

We first examine the eigenstates deviating slightly from half-filling, namely the single particle excitations and single hole excitations. They are represented by $e_{k,\downarrow}^{{\dagger}}\left|\text{FM}\uparrow\right\rangle$ [of energy $\epsilon_{p}\left(k\right)$] and $e_{k,\uparrow}\left|\text{FM}\uparrow\right\rangle$ [of energy $\epsilon_{h}\left(k\right)$] respectively. Using the definitions Eqs. (24) and (23) and the fact that $\delta_{n}+\delta_{m}$ is even, it is easy to show that $\epsilon_{p/h}\left(k\right)=\epsilon_{p/h}\left(2\pi-k\right)$, so we may focus on $2\pi/3\leq k\leq\pi$.

Near the Dirac point $0<k-2\pi/3\ll 1$, we can expand Eq. (24) to obtain

[TABLE]

where the velocity $v$ depends only on the interactions:

[TABLE]

Since the $k^{\prime}$ integral is finite, $v$ is finite for any short-range interaction. Eq. (25) shows that, as in the projected Hubbard model,Karimi and Affleck (2012) the single-particle and single-hole excitations are generally gapless at the Dirac points for a single zigzag edge.

For the next-nearest-neighbor model Eq. (15), $v$ is always positive:

[TABLE]

Nevertheless, $v$ may become negative for certain strongly non-local interactions. An example is the term with $\delta_{m}=4$ and $\delta_{n}=0$, which gives a coefficient of $-3/\left(16\pi\right)$. The ferromagnetic ground state will be unstable against the creation of electrons or holes near the Dirac points in the case of $v<0$, or more generally $v<\sqrt{3}\left|\Delta\right|$ when the particle-hole symmetry breaking term $\Delta$ is nonzero.

The case of unscreened Coulomb interaction Eq. (6) is especially interesting. In this case the low-energy behavior of $\epsilon_{p/h}\left(k\right)$ is controlled by the long range part of $U_{\left(\delta_{m},\delta_{n}\right)}$. When $\left|\delta_{n}\right|\gg 1$ or $\left|\delta_{m}\right|\gg 1$, the $k^{\prime}$ integral is dominated by $k^{\prime}$ near the Dirac points, and we findPereira et al. (2006)

[TABLE]

Approximating the sum over $\delta_{n}$ and $\delta_{m}$ by integrals over $x=\delta_{m}/2$ and $y=\sqrt{3}\delta_{n}/2$, and discarding the subleading contribution from the oscillating term, we see that $v\propto\ln R$ where $R$ is the long-distance cutoff:

[TABLE]

As $R\rightarrow\infty$, the only other large distance scale in the problem is given by the inverse distance to the Dirac points, which should therefore replace $R$ as the distance cutoff. In other words, for the unscreened Coulomb interaction, $\epsilon_{p/h}$ has the following behavior for $0<k-2\pi/3\ll 1$:

[TABLE]

where $\Lambda\ll 1$ is a momentum cutoff. This behavior is not affected by the particle-hole symmetry breaking term $\Delta$, which merely shifts $\Lambda$.

In Fig. 5 we plot $\epsilon_{p/h}\left(k\right)/\left(k-2\pi/3\right)$ versus $\ln\left(k-2\pi/3\right)$ at $0<k-2\pi/3\ll 1$ for the Coulomb interaction with various $R$, and show how the logarithmic divergence in Eq. (30) is cut off at low energies by $R$. We also plot the velocity $v$ given by Eq. (26) as a function of $\ln R$ in Fig. 6. These results suggest that the Coulomb interaction produces a divergent “Fermi velocity” for edge modes near the Dirac points, a behavior reminiscent of the marginal Fermi liquid in bulk graphene with Coulomb interaction.González et al. (1999)

It is also useful to consider $k=\pi$, since this is where $\epsilon_{p/h}\left(k\right)$ obtains its maximum for the Hubbard interaction and the Coulomb interaction, in the absence of particle-hole symmetry breaking. At $k=\pi$ Eq. (24) is again greatly simplified:

[TABLE]

where the sum is over even $\delta_{m}$ with $\delta_{m}\geq 4$.

$\epsilon_{p/h}\left(\pi\right)$ is also finite for any short-range interaction. Interestingly, $\epsilon_{p/h}\left(\pi\right)$ depends on $U_{\left(\delta_{m},\delta_{n}\right)}$ only if $\delta_{n}=0$: it is, for instance, independent of $U_{2\angle}$ in the next-nearest-neighbor model Eq. (15). For the Coulomb interaction Eq. (6), the $\delta_{m}$ sum turns out to be convergent, and we find

[TABLE]

where $U_{0}$ is the on-site interaction strength. Therefore, when $\Delta>\Delta_{c}\left(R\rightarrow\infty\right)=0.189U_{0}$ for the unscreened Coulomb interaction, the maximum $S_{z}$ state $\left|\text{FM}\uparrow\right\rangle$ becomes unstable towards the creation of a spin-down electron at $k=\pi$, e.g. by absorption from the bulk. [For Hubbard interaction with strength $U$, the condition is $\Delta>\Delta_{c}\left(R=0\right)=\left(\sqrt{3}/\left(2\pi\right)-1/6\right)U\approx 0.109U$.]Karimi and Affleck (2012) Similarly, when $\Delta<-\Delta_{c}\left(R\rightarrow\infty\right)$ there is an instability towards the creation of a spin-up hole at $k=\pi$.

In Fig. 7 we plot $\epsilon_{p}\left(k\right)$ versus $k$ for $2\pi/3\leq k\leq 4\pi/3$ for the Coulomb interaction with different values of long-distance cutoff $R$, both when $\Delta=0$ and when $\Delta=\Delta_{c}\left(R\right)$ so that $\epsilon_{p}\left(\pi\right)$ vanishes. Notice that for the Coulomb interaction $\epsilon_{p}\left(k\right)>0$ for $0<k-2\pi/3\ll 1$ even when $\Delta=\Delta_{c}\left(R\right)$; that is, as we increase $\left|\Delta\right|$, single particle or single hole creation energy becomes negative at $k=\pi$ sooner than it does near the Dirac points.

IV.2 1-particle-1-hole sector

We turn to the half-filled sector with $N-1$ spin-up electrons and $1$ spin-down electron, so that $S_{z}=N/2-1$. This sector hosts $1$ spin-down electron and $1$ spin-up hole relative to the $\left|\text{FM}\uparrow\right\rangle$ state, and accommodates the excitations that would be seen as magnons in an effective spin model.

Let the total momentum relative to $\left|\text{FM}\uparrow\right\rangle$ be $Q$, and without loss of generality we assume $Q\geq 0$. Denoting an eigenstate by

[TABLE]

we obtain the following Schroedinger’s equation:

[TABLE]

where $E$ is the energy eigenvalue. The ferromagnetic state in the 1-particle-1-hole sector, $f\left(k;Q=0\right)=1$, is obviously a zero-energy solution.

It is possible for $f\left(k,Q\right)$ to have a $\delta$-function peak at $k=k_{0}$. In this case the solution to Eq. (34) is part of the 1-particle-1-hole continuum, and has an energy $E=\epsilon_{h}\left(k_{0}\right)+\epsilon_{p}\left(k_{0}+Q\right)$. Another possibility is having $E<\epsilon_{h}\left(k\right)+\epsilon_{p}\left(k+Q\right)$ for any $k$, in which case $f\left(k;Q\right)$ does not have any $\delta$-function peaks, and the solution is a particle-hole bound state, or an exciton. Since it reduces $S_{z}$ by $1$, it can also be viewed as a magnon in an effective spin model.

For any short range interaction, we can show that the exciton energy has the following $Q\ll 1$ behavior:

[TABLE]

where $v$ is the velocity Eq. (26) that also appears in the single particle dispersion, and $\Lambda^{\prime}\ll 1$ is again a momentum cutoff. The inverse exciton mass, or the spin stiffness of the ferromagnetic zigzag edge, is therefore logarithmically divergent. The derivation of Eq. (35) is sketched in Appendix A, where we see the divergence arises due to the linear behavior of $\epsilon_{p/h}\left(k\right)$ near the Dirac points. This divergence is possibly related to the large spin stiffness found by Refs. Yazyev and Katsnelson, 2008; Rhim and Moon, 2009 for $U$ comparable to $t$.

Although similar exciton dispersions have been previously reported in carbon nanotubes,Jiang et al. (2007); *PhysRevLett.106.136805; *PhysRevLett.109.187403 in contrast to Eq. (35) they originate from the long-range nature of the Coulomb interaction. In fact, since in the Coulomb interaction with a long-distance cutoff $R$ we have $v\propto\ln R$, we expect that Eq. (35) is modified to $E\left(Q\right)\propto Q^{2}\ln^{2}Q$ for $R\rightarrow\infty$; that is, the spin stiffness is even more divergent than a logarithm for the unscreened Coulomb interaction. Fig. 8 shows $E\left(Q\right)/Q^{2}$ plotted against $\ln Q$ at $0<Q\ll 1$ for some values of $R$ and $\Delta=0$, where $E\left(Q\right)$ is found by solving Eq. (34) numerically via Chebyshev series expansion.Abramowitz and Stegun (2012)

It is also helpful to examine the effect of $\Delta$ on the exciton dispersion, taking as an example the Coulomb interaction with a long-distance cutoff $R$. As depicted in Fig. 9, when $\left|\Delta\right|=\Delta_{c}\left(R\right)$ so that $\epsilon_{p}\left(\pi\right)=0$, the exciton dispersion $E\left(Q\right)$ calculated numerically also approximately vanishes at $Q=\pm\pi/3$, and the exciton wave function strongly favors the state with a spin-down electron at $\pi$ and a spin-up hole at either Dirac point. For $\left|\Delta\right|>\Delta_{c}\left(R\right)$, in parallel with the Hubbard case,Karimi and Affleck (2012) $E\left(\pm\pi/3\right)$ becomes negative which indicates that the ground state at half-filling is no longer maximally spin polarized; instead, the edge states near $\pi$ become more likely to be doubly occupied and the edge states near the Dirac points become more likely to be unoccupied.

Finally, we mention that in the 2-particle-2-hole sector, the excitons in the 1-particle-1-hole sector can form an additional bound state below the exciton continuum. Nevertheless, for both the Hubbard and the Coulomb interactions with $\left|\Delta\right|<\Delta_{c}\left(R\right)$, we find numerically that the bottom of the 2-particle-2-hole bound state dispersion remains positive; we thus conjecture that the ferromagnetic ground state is stable for $\left|\Delta\right|$ up to $\Delta_{c}\left(R\right)$. We also mention that the bound state picture provides an intuitive explanation for the non-ferromagnetic regime in Fig. 3: for $M<N/2$, we can usually form an $M$-particle-$M$-hole bound state with a non-negative binding energy, i.e. with an energy lower than or equal to the sum of energies of an $\left(M-1\right)$-particle-$\left(M-1\right)$-hole bound state and a 1-particle-1-hole bound state. Therefore, if the 1-particle-1-hole ground state has a negative energy as happens for sufficiently large $U_{2\angle}$, then as $M$ increases and $\left|S_{z}\right|$ decreases, the ground state energy in the $S_{z}$ sector either stays the same or decreases.

V Discussion and conclusions

In our effective model Eq. (8) at $O\left(U\right)$, we have ignored the dynamics of low-lying bulk degrees of freedom near the Dirac points, so an obvious issue is whether this approximation is justified. For the on-site Hubbard interaction, the answer is partly given by Refs. Karimi and Affleck, 2012; Koop and Schmidt, 2015, where effective Hamiltonians are found to $O\left(U^{2}/t\right)$ by integrating out the bulk states and neglecting retardation. While Ref. Karimi and Affleck, 2012 finds that the $O\left(U^{2}/t\right)$ correction to the Hamiltonian has a $q^{2}\ln q$ behavior for small momentum transfer $q$, such behavior does not necessarily hint at a breakdown of the perturbation theory, as logarithms also appear at $O\left(U\right)$, e.g. in the exciton dispersion Eq. (35). Ref. Koop and Schmidt, 2015 further shows that, as far as the effective spin model is concerned, the interaction strengths are only weakly modified by the bulk states even for $U$ comparable to $t$. In other words, there is no evidence that the perturbation theory in $U/t$ is divergent. However, while a weak Hubbard interaction is known to be irrelevant in the bulk, a weak Coulomb interaction is marginally irrelevant and may lead to further logarithmic corrections.González et al. (1999); Kotov et al. (2012) It therefore remains an open question whether integrating out the bulk states at $O\left(U^{2}/t\right)$ qualitatively changes the physics of the $O\left(U\right)$ edge model for the unscreened Coulomb interaction.

Another problem that we have not discussed so far is the inter-edge coupling in realistic graphene nanoribbons. We now consider a ribbon of large but finite width $W\gg 1$ with two zigzag edges, whose overall ground state is antiferromagnetic. The inter-edge coupling originates in part from the direct interaction between opposite edges, which is significant even at the first order in interaction if it is long-ranged [$O\left(U_{0}/W\right)$ in the Coulomb case]. Inter-edge coupling is also mediated by bulk states, which is second order in interaction and is $O\left(U^{2}/\left(tW^{2}\right)\right)$ in the Hubbard case.Karimi and Affleck (2012) Yet another source is the hopping amplitude between edge states of opposite edges, which exists even in the absence of interactions and leads to an energy gap exponentially small in $W$. For wide ribbons $W\gg 1$, it is well known that the edge states are no longer strictly localized near one edge when their momenta are within $O\left(1/W\right)$ of the Dirac points. The hopping amplitude at momentum $k$ thus grows rapidly as $k$ approaches the Dirac points, eventually reaching $O\left(t/W\right)$.Brey and Fertig (2006); Jung et al. (2009) Under our assumption $U\ll t$, this is actually a much larger energy scale than that of the direct inter-edge Coulomb interaction or that of the bulk-mediated inter-edge interaction. Thus it is not justified to ignore the inter-edge hopping amplitude near the Dirac points in the effective model for a nanoribbon. In fact, at the mean-field level, it is exactly the part of the Brillouin zone near the Dirac points that contributes the most to the inter-edge superexchange interaction,Jung et al. (2009); Jung (2011) and the spin wave dispersion becomes linear for small momenta once the inter-edge coupling is taken into account.Culchac et al. (2011) Although an effective edge model incorporating the inter-edge hoppingSchmidt et al. (2013) is often much less analytically accessible beyond the mean-field level, we hope further insight on the effect of Coulomb interaction in finite width nanoribbons can be gleaned from exact diagonalization.

In conclusion, we have investigated the effects of long-range interactions on the zigzag edge states of a semi-infinite graphene sheet. By projecting the interaction onto the edge state subspace, we obtain an effective model for which the states in the maximally polarized ferromagnetic multiplet are zero energy eigenstates. A sufficient condition is found for the ferromagnetic multiplet to be the ground states, and we present evidence that the unscreened Coulomb interaction satisfies this condition, which implies that its ground states are ferromagnetic. In cases where the sufficient condition is not met, exact diagonalization results indicate that the ground state can be non-ferromagnetic, provided that certain non-local components of the interaction are sufficiently strong. We also discuss the single-particle excitations, single-hole excitations and spin-1 excitons of the maximum $S_{z}$ ground state. For short range interactions the single-particle and single-hole excitations have linear dispersions near the Dirac points, as described in Eq. (25). The slope $v$ also governs the exciton energy at small momenta, Eq. (35), which shows a $vQ^{2}\ln Q$ behavior. For the unscreened Coulomb interaction $v$ becomes logarithmically divergent as a function of the long-distance cutoff, corresponding to a $\delta k\ln\delta k$ behavior where $\delta k\ll 1$ is the distance from either of the Dirac points. The edge states acquire a dispersion due to a particle-hole symmetry breaking perturbation $\Delta$; the ferromagnetic ground state can be destroyed if $\left|\Delta\right|$ is large enough.

Acknowledgements.

This work was supported in part by NSERC of Canada, Discovery Grant 04033-2016 and the Canadian Institute for Advanced Research. ZS would like to acknowledge helpful discussions with Emilian Nica.

Appendix A Exciton dispersion at small momenta for a single zigzag edge with

short-range interactions

For simplicity, we illustrate the derivation of Eq. (35) with the Hubbard interaction $U_{\left(0,0\right)}=U$. Generalization to non-local interactions is tedious but straightforward; it is briefly discussed at the end of this appendix.

Expanding the denominator of the kernel $\Gamma$ in Eq. (34), we can isolate the $k$ dependence of $f\left(k;Q\right)$:

[TABLE]

where $\Phi$’s are independent of $k$, and are defined as

[TABLE]

Inserting Eq. (36) into Eq. (37), we obtain an infinite number of linear equations satisfied by $\Phi$:

[TABLE]

For $Q\ll 1$, the integrand on the right-hand side can be expanded to $O\left(E\right)$ and $O\left(Q^{2}\right)$.

[TABLE]

In the second and the third lines we have approximated $\Phi_{l^{\prime}}\left(Q\right)\approx\Phi_{l^{\prime}}\left(0\right)$, assuming that $\Phi_{l}\left(Q\right)$ is well-behaved at $Q=0$ and any difference is $O\left(Q\right)$. In the third line we have retained the most singular contribution at $O\left(Q^{2}\right)$, which are from the vicinity of the Dirac points (hence the factor of $2$), as the remaining terms contain no infrared divergence.

Using Eq. (37) and recalling that the $Q=0$ solution is $f\left(k;0\right)=1$, we have

[TABLE]

and

[TABLE]

therefore

[TABLE]

Now, we multiply the entire expression by $\left[1-4\cos^{2}\left(k/2\right)\right]\left[4\cos^{2}\left(k/2\right)\right]^{l}U$, then sum over $l$ and integrate over $k$. The left hand side then cancels the first term on the right hand side, and using $v=U/\left(2\sqrt{3}\right)$, we are left with Eq. (35).

In the presence of non-local interactions, one needs to assign three more indices to $\Phi$, namely $\delta_{m}$, $\delta_{n}$ and $\alpha=1$, $2$ [corresponding to the two terms in the third line of Eq. (23)]. All three indices should be summed over in Eq. (38), and subsequently in Eqs. (40), (41) and (42).