Zero-modes and global antiferromagnetism in strained graphene

A novel magnetic ground state is reported for the Hubbard Hamiltonian in strained graphene. When the chemical potential lies close to the Dirac point, the ground state exhibits locally both the N\'{e}el and ferromagnetic orders, even for weak Hubbard interaction. Whereas the N\'{e}el order parameter remains of the same sign in the entire system, the magnetization at the boundary takes the opposite sign from the bulk. The total magnetization this way vanishes, and the magnetic ground state is globally only an antiferromagnet. This peculiar ordering stems from the nature of the strain-induced single particle zero-energy states, which have support on one sublattice of the honeycomb lattice in the bulk, and on the other sublattice near the boundary of a finite system. We support our claim with the self-consistent numerical calculation of the order parameters, as well as by the Monte Carlo simulations of the Hubbard model in both uniformly and non-uniformly strained honeycomb lattice. The present result is contrasted with the magnetic ground state of the same Hubbard model in the presence of a true magnetic field (and for vanishing Zeeman coupling), which is exclusively N\'{e}el ordered, with zero local magnetization everywhere in the system.


I. INTRODUCTION
Despite the manifest stability of the Dirac fermions in graphene against the effects of Coulomb interaction, 1 the nature of the possible broken symmetry phases at strong coupling continues to be an issue of fundamental importance [2][3][4][5][6][7] . The minimal onsite Hubbard interaction, for example, when sufficiently strong, is believed to produce the antiferromagnetic Néel ground state 2,4,[8][9][10][11] . The universality class of the semimetal-to-Néel insulator quantum phase transition can be captured by the effective Gross-Neveu-Yukawa field theory 12 , and studied systematically near the upper critical (spatial) dimension of three. Recent quantum Monte Carlo simulations of the Hubbard model in a half-filled honeycomb lattice suggest a direct transition from the Dirac semimetal to the Néel state, with the critical exponents in reasonable agreement with the field-theoretic predictions. 9,10 The correlated phases of graphene, Néel state included, unfortunately may be laying at too strong a coupling to be realized in graphene in its pristine state, even when placed in vacuum. 13 The deformability of the graphene membrane may facilitate a different path towards the realization of some of the symmetry-broken phases at weaker couplings. 14 Arguments along this line and in favor of the topological quantum anomalous (spin) Hall insulator at * Current affiliation weak finite-range repulsion [14][15][16] , or of unconventional superconductors 17 at weak attraction, have recently been put forward. The physical reason behind this electromechanical phenomenon is the generic appearance of the single-particle zero-energy states in graphene under strain. Upon bulging a graphene flake, the quasirelativistic low-energy electronic degrees of freedom couple to a finite, time-reversal-symmetric, magnetic-like field. 18 Such an, and not necessarily uniform, axial magnetic field, similarly to the true magnetic field, falls under the jurisdiction of index theorems, 19 and as such brings a finite number of states close to the Dirac point. This creates an ideal situation for the fermions to form various particle-hole or particle-particle condensates at weak interactions. The single-particle states at zero energy, responsible for the weak-coupling instabilities, are in the axial case, however, special: 14 normalizability forces them to live exclusively on one of the two sublattices of the honeycomb lattice. The remaining states in the zero-energy subspace, which would be discarded from the spectrum as non-normalizable in an infinite system, in a finite system are found near the boundary, and to be living on the opposite sublattice.
Due to this inextricable correlation between the realspace and sublattice degrees of freedom within the zeroenergy subspace the ground state of the Hubbard model in strained graphene is non-trivial. While it seems natural to expect that a short-range interaction such as Hubbard's would lead to a spin-polarized "Hund" ferromagnetic state in the presence of a flat band in strained graphene, we show here, first via a self-consistent numerical calculations, that the ground state of the Hubbard model in this system is more interesting: while it is locally displaying the expected ferromagnetic ordering, the sign of the magnetization varies in precisely such a way so that the total space-integrated magnetization in fact vanishes. Since the magnetization is tied to the zero-modes, however, its support, as well as its sign, also switches between the two sublattices when traversing the system from its bulk to the boundary, so that the Néel order parameter is also finite, and of the same sign everywhere. The appearance of this unusual magnetic ordering is also supported by a quantum Monte Carlo calculation on the Hubbard model on a half-filled strained honeycomb lattice. We name this unconventional magnetic ordering the global (edge-compensated) antiferromagnet. Recent experimental progress in realizing the axial magnetic field in real and artificial graphene, 20-22 offers hope for the detection of this unusual ground state.
The paper is organized as follows. In the next section we propose a specific modulation of the nearest-neighbor hopping amplitude that introduces a finite axial magnetic field in a graphene flake. We also discuss and compare the zero-modes in the presence of the strain-induced axial, and true magnetic fields. In Sec. III we introduce the on-site Hubbard interaction, and discuss possible magnetic ground states in strained graphene. In Sec. IV we present the self-consistent Hatree solution for the magnetic ground state in the presence of (roughly) uniform and well localized axial magnetic fields. We contrast our results with the magnetic ground state in graphene in a true (time-reversal symmetry breaking) magnetic field in Sec. V. Quantum Monte Carlo simulation of the Hubbard model in strained graphene is presented in Sec. VI. We summarize our findings in Sec. VII.

II. AXIAL MAGNETIC FIELD AND ZERO-MODES
Before delving into the effects of the electron-electron interactions, let us set the stage by reviewing the physics of zero-modes in the continuum and on the lattice, and by introducing a specific realization of the axial magnetic fields on a finite honeycomb lattice. As well known, the low energy degrees of freedom in graphene may be collected into an 8-component Dirac spinor Ψ , and σ =↑, ↓ are the two projections of electron spin along the z-axis. 4 u σ is the fermionic annihilation operator on A-sublattice, generated by the linear combination of basis vectors a 1 = ( √ 3, −1)a, a 2 = (0, 1)a, where a is the lattice spacing, and v † σ is the fermionic creation operator on B-sublattice, where B = A + c, with c = (1/ √ 3, 1)(a/2). Two inequivalent Dirac points, or valleys, may be chosen as at ± K, . When | q| ≪ | K|, the non-interacting low-energy Hamiltonian with only the nearest-neighbor hopping takes the relativistically invariant form H D = σ 0 ⊗ iγ 0 γ jqj , with the first matrix acting on spin, and the second on sublattice and valley indices.
The coupling of the time-reversal-symmetric axial magnetic field b( x) to the low-energy Dirac fermions then reads as 14,23-25 The axial vector potential here is a i ( x) = ǫ ij ∂ j χ( x), and therefore b( x) = ∂ 2 χ( x). The mutually anticommuting 4-dimensional γ-matrices may be represented as (σ 0 , σ) are the usual two-dimensional Pauli matrices. In this representation the time-reversal symmetry operator is (anti-linear) I K = σ 2 ⊗ iγ 1 γ 5 K, where K is the complex conjugation. 4 A random distribution of the axial gauge field in graphene, for example, results from the presence of ripples, with the net axial flux as zero. If graphene is deliberately buckled, on the other hand, a finite total flux of the axial field may be introduced, which by index theorem would bring a finite number of states at zero energy. These special zero-modes can be written as The number of zero energy states labeled by index n equals the total axial flux through the system. The matrix γ 0 in the exponent changes sign between two sublattices, whereas the function χ( x) is a monotonic function at a large distance | x| from the location of the axial flux. The normalizable zero-energy states therefore must reside only on one of the two sublattices, which we will call the sublattice A. On the finite honeycomb lattice, these are the bound states in the interior of the system, where the flux is located. The remaining non-normalizable zero energy states, with the support on the sublattice B, on the other hand, in the continuum increase exponentially towards the infinity. On a finite lattice, however, this is tantamount to their localization near the boundary of the system. Eq. (1) suggests an introduction of an axial magnetic field on a lattice via the following modification of the nearest-neighbor hopping amplitude (t) where i ∈ A, j ∈ B. 16,26 Hereafter we set t = 1 for convenience. Let us define a quantity, say R, which counts the minimal number of bonds required to reach a particular site in the system from the central hexagon (centered at (j x , j y ) = (0, 0)) in Fig. 1. For all the six sites on the central hexagon of the system R = 0, for example. We then assign the parameter χ(R), such that for all the sites with same R, χ(R) is same, as shown in Fig. 1, and the hopping between nearest-neighbor sites are modified according to Eq. (3). Otherwise, χ(R) increases monotonically from the center towards the boundary of the system. As a result, a finite axial magnetic field is introduced in the system. For example, if χ(R) ∼ R 2 the system experiences a roughly uniform axial field, whereas a bell-shaped localized axial flux around the center of the system can be realized by setting χ(R) ∼ log R. This configuration of χ(R) is slightly different than in the previous work 16 , with the advantage that the increase of the band width in the presence of axial fields can be somewhat better controlled here. In a finite strained honeycomb lattice there are of course no states at precise zero energy. Nevertheless, irrespective of the spatial profile of the axial magnetic field, there are always a finite number of states in the close vicinity of the zero energy, 16 which bear the signature of the axial magnetic field. Lets consider two such states | ± δ , with the energies ±δ, where |δ| ∼ 0. The symmetric |S and the anti-symmetric |AS combinations of these two states live on A and B sub-lattice, respectively. In the presence of axial magnetic fields, |S lives inside the bulk and corresponds to the normalizable zero energy states in Eq. (2). The |AS lives near the boundary of a finite system, and is to be identified as the nonnormalizable state in the continuum description. Hence, the bulk and the boundary of a finite graphene system become inequivalent in their sublattice structure in the presence of axial magnetic fields.
It is worth contrasting the structure of the zero energy subspaces in the presence of true vs. axial magnetic field. Coupling of Dirac fermions to the true magnetic field involves the same Eq. (1), except for the replacement of the matrix iγ 3 γ 5 by the 4-dimensional unit matrix, and for taking the matrix M (χ) to be exp Likewise, the zero energy states in the true magnetic field still have the form of Ψ 0,n in Eq. (2), but with the matrix γ 0 in the exponent replaced by iγ 3 γ 5 . As a consequence, the normalizable zero-modes now reside on both sublattices. In a finite graphene system in a true magnetic field, the near-zero-energy states, residing in the bulk, similarly populate equally both sublattices. The non-normalizable zero-modes in the continuum, located near the boundary in a finite system, are also shared equally between the two sublattices. Hence, in stark contrast to the axial field, in the presence of true magnetic field, the interior and the boundary of a finite honeycomb lattice have the same unresolved sublattice structure.

III. MAGNETIC GROUND STATES WITH HUBBARD INTERACTION
Let us focus next on the effect of electron-electron interaction in strained graphene, with the chemical potential (µ) tuned to the charge neutrality point. Here we consider only the onsite Hubbard interaction (U ) among the fermions. The standard Hubbard Hamiltonian is with the interaction Hamiltonian given by and where n i,σ is the fermionic number operator with spin projection σ =↑, ↓ at site i. The total number of electrons in the system is N . The charge-neutrality in the system is maintained through the constraint µ = 0. The usual Hartree decomposition of the onsite interaction leads to an effective single particle Hamiltonian One can rewrite n x,σ as where σ = (+, −) respectively represents (↑, ↓) projections of electron's spin along z-axis. δ x,σ (r) corresponds to the site-dependent local deviation of the electronic density from the uniform background, and is to be determined self-consistently in a finite honeycomb lattice in presence of the axial magnetic field representing strain. We here always take the system to be quasi-circular, and r = 1, 2, · · · , n is a discrete variable, representing the n th ring of around center of the system. In the presence of ordering (δ = 0) the overall charge neutrality of the system is achieved through the constraint r σ=± in addition to µ = 0. Notice that the Fock decoupling of the Hubbard term would yield expectation values of the components of the local magnetization that are orthogonal to the chosen (z) direction, which we neglect. Depending on the relative sings of the site parameters δs, one can realize two different magnetic ground states: (i) If the ground state configuration is such that all δ A/B,↑/↓ > 0 in the above equation, then this would correspond to an antiferromagnetic phase, whereas (ii) a ground state with δ A,↑/↓ > 0, but δ B,↑/↓ < 0 would be identified as a ferromagnetic state. Before we proceed with the numerical simulation of the Hubbard model, it is worth pausing to compare these two magnetic phases in strained graphene. Magnetization on the two triangular sublattices points in the opposite directions when the ground state is antiferromagnetic, whereas in a ferromagnet the two sublattice magnetizations are aligned. Recall that all the zero-modes in strained graphene are localized on one sublattice in the bulk, and on the other sublattice near the boundary of a finite system. Therefore, both the Néel and the ferromagnetic orders in strained graphene give rise to a finite local magnetization everywhere in the system. However, globally the two states may be distinguished.
Due to the spatial separation of the zero-modes that are localized on different sublattices, the magnetization in the state with an equal sign of the Néel order parameter in the entire system would need to change sign as one approaches the boundary from the bulk, so that the total magnetization could in fact vanish. In the truly ferromagnetic state, on the other hand, the sign of the magnetization in the bulk and the boundary would be the same, so that the system would develop an overall finite magnetization. We will show shortly through a detailed numerical calculation that by taking into account the entire zero subspace in a finite system, the magnetic ground state in strained graphene is uniquely determined to be of the first variety, i. e, a global antiferromagnet with zero total magnetization. Nevertheless, the continuum picture already provides a valuable insight into the nature of the competition between the Néel and the ferromagnetic states. The order parameters of these two states read as N = Ψ † ( σ ⊗ γ 0 )Ψ , and F = Ψ † ( σ ⊗ I 4 )Ψ , respectively. Both states split the zero-energy subspace in strained graphene, and open a gap at the Dirac points. However, the matrix appearing in the Néel order parameter anticommutes with the Dirac Hamiltonian, and as such it corresponds to a chiral-symmetry-breaking massterm for the Dirac quasi-particles. In the ferromagnetic order parameter, in contrast, the corresponding matrix commutes with the Dirac Hamiltonian. As a result, besides the splitting of the zero-energy subspace common to both orders, the Néel order at the mean-field level pushes down in energy all the filled states below the chemical potential. In contrast, the ferromagnetic order parameter splits all the energy levels equally, half up and half down in energy, and therefore only lowers the energy of the half-filled zero-energy subspace. By spontaneously developing the Néel order the system can thus more efficiently minimize the ground state energy.

IV. SELF-CONSISTENT CALCULATION
We now present the results of the self-consistent calculation of the magnetic order parameters with onsite Hubbard interaction (U ) in strained graphene. We numerically search for the self-consistent solution of δs with two different initial ansatz for the magnetic ordering (i) when δ A/B,↑/↓ > 0, which can be idenified as antiferromagnet ordering, and (ii) δ A,↑/↓ > 0, δ B,↑/↓ < 0, which corresponds to a ferromagnetic order. Here, all the δs are kept as a function of position, and we always maintain the overall charge-neutrality of the system. In the Hartree self-consistent calculation, electronic spin is treated as Ising variable, σ → σ 3 , and the effect of fluctuation is neglected for a moment. Later we treat the fermionic spin as O(3) vectors in a separate quantum Monte Carlo simulation of the Hubbard model in strained graphene, which explicitly takes into account the effect of the fluctuations. We here search for the self-consistent solution of the magnetic order specifically for weak Hubbard interactions, U ≪ U c , where U c ≈ 3.8t is the zero axial field critical strength of the onsite interaction for antiferromagnetic ordering 10 . Due to the presence of a finite density of state near the zero energy in strained graphene, ordering takes place even for onsite interaction as weak as U/t = 0.1, irrespective of the spatial profile of the axial magnetic field, and a spectral gap opens up at the Dirac points. The resulting magnetic ground states is insensitive to the initial ansatz (i) or (ii), and thus our self-consistent analysis can be considered as variational calculation. To further explore the nature of such magnetic ground state, we define two local order parameters as which correspond to local Néel and ferromagnetic order parameters, respectively. δs are as defined in Eq. (7). We here obtain the numerical results in a quasi-circular honeycomb lattice of 600 sites or r max = 10. In such a system the self-consistent solutions for all the δs are essentially without any finite size effects, for all values of the sub-critical Hubbard interactions, down to U = 0.1, and for both the uniform and localized axial fields (see the captions of Fig. 2 and Fig. 3 for details of these parameters). Self-consistent solutions of N (r) and M (r) in the presence of roughly uniform axial magnetic fields and a wide range of sub-critical on-site interaction are shown in Fig. 2 (upper panel) and Fig. 2 (lower panel), respectively. From the spatial variation of these two order parameters we see that the local antiferromagnetic order parameter N (r) is off the same sign in the bulk as well as in the boundary of the system, whereas M (r) near the boundary is of the opposite sign from the bulk. The total magnetization for all chosen values of the interaction and uniform axial fields is zero within the numerical accuracy ∼ 10 −12 . Therefore, the magnetic ground state is indeed an antiferromagnet, as one would anticipate from the continuum description of this problem. The Néel order in strained graphene is different from the conventional one on the honeycomb lattice 10 as it also carries a finite local magnetization everywhere in the system. We dub this unconventional magnetic ground state global (edge-compensated) antiferromagnet.
We also obtained the self-consistent solution for the magnetic ground state when the graphene flake is subject to non-uniform axial fields, localized in the vicinity of the center of the system. The spatial variation of the local order parameters N (r) and M (r), for a wide range of sub-critical interactions, and total flux of localized axial magnetic field is shown in Fig. 3 (upper panel), and Fig. 3 (lower panel), respectively, obtained in a quasi-circular system of 600 sites (r max = 10). From Figs. 2 and 3, it is clear that the nature of the magnetic ground state is insensitive to the exact profile of the axial field, and it is always the global antiferromagnet. However, the location where the local magnetization changes its sign depends on the profile of the axial field. In the presence of uniform axial field, magnetization flips its sign only very close to the boundary of the system, whereas in the presence of localized fields the same change occurs in the middle of the system. The difference in the position of the domain wall between the two signs of the magnetization is tied to the spatial distribution of the zero energy states in the system. In the uniform axial field the zero-modes on A-sublattice are distributed over a large portion of the bulk, whereas they are highly localized in the deep interior when the axial field is peaked near the center of the system. Therefore, the local magnetization suffers a change in sign roughly where the zero energy states on A-sublattice loose their support. The magnitude of both of the order parameters N (r) and M (r) increases with the strength of the interaction and of the axial field, as shown in Figs. 2 and 3.
The spatial distribution of the zero-modes also dictates spatial variation of the Néel order parameter. Note that in the presence of uniform axial fields, the Néel order parameter N (r) is more or less uniform in the bulk (r ≤ 6) of the system (see Fig. 2, the upper panel). The abrupt increment in N (r) near the boundary of the system arises from the existence of zero-energy states on the B-sublattice in that region. When the axial field assumes a spatially localized profile, around the center of the system, the Néel order parameter predominantly develops in the region where the axial flux penetrates the graphene flake; see Fig. 3 (the upper panel). The existence of the zero-energy states on the sites of the Bsublattice near the boundary leads to the spikes in N (r) even when the axial field is inhomogeneous. Therefore, the spatial modulation of the Néel order follows closely the profile of the axial magnetic field, resembling in this regard the spatial variation of the quantum anomalous Hall insulator in strained graphene with next-nearestneighbor interaction.

V. HUBBARD MODEL IN TRUE MAGNETIC FIELD
The conical dispersion of the Dirac quasi-particles quenches into a set of well separated Landau levels, when the graphene flake is subject to a uniform true magnetic field. Although the discrete quantization of the spectrum smears out if the magnetic field is spatially modulated, a finite number of states at zero energy states always persists, irrespective of the profile of the magnetic fields. 19 Finite density of states near the Dirac points, here as well, catalyzes the effect of electron-electron interactions. To compare the present unconventional magnetic ground state of the Hubbard model in strained graphene with the one in the presence of a true magnetic field, we perform the same numerical self-consistent analysis in a finite honeycomb lattice, placed in uniform and nonuniform true magnetic fields. The orbital effect of the true magnetic field is included by incorporating the Peierls phase into the nearest-neighbor hopping amplitudes. 27 We here omit the single-particle Zeeman coupling of electron's spin with the magnetic field. 28 In true magnetic field the zero energy states are found on both sublattices in the bulk, as well as near the boundary of the system. In this situation, only the antiferromagnetic ansatz (i) leads to a finite gap at the Dirac point.
From the self consistent solution of the magnetic ground state, obtained in a finite honeycomb system of 384 sites (r max = 8), we compute the two local order parameters N (r) and M (r), defined in Eq. (9). The results are presented in Fig. 4, when the field is uniform (left) and nonuniform (once again a bell-shaped one, peaked around the center of the system) (right). The local Néel order parameter N (r) again follows the profile of the true magnetic field, resembling in this regard the spatial distribution of the charge-density-wave order, obtained previously for spinless fermions in honeycomb lattice with the nearest-neighbor interaction. 27 Due to the existence of a finite density of states at the Dirac point, the antiferromagnetic ordering can be found for the interaction as weak as U = 0.4. 29 On the other hand, M (r) is zero everywhere in the system. The dramatic difference in the magnetic ground states between the axial and the true magnetic fields arises entirely as a consequence of the distinct structure of the zero-modes. The antiferromagnetic ground state we find in the presence of true magnetic field and at weak interaction (U ≪ U c ) is the exact replica of the one in graphene at strong interaction (U ≥ 3.8t) and in zero field.

VI. MONTE CARLO CALCULATION
Axial magnetic fields, resulting from the modification of the hopping matrix elements between the nearest neighbors according to Eq. (3), do not break the particlehole symmetry. Thereby, auxiliary field quantum Monte Carlo simulations do not suffer from the infamous minus sign problem and accurate simulations on large lattices can be carried out. Here, we have used the projective zero temperature approach based on the equation, in which the ground state is filtered out of a single Slater determinant by propagating along the imaginary time axis. It is beyond the scope of this paper to go into the details of the implementation and the readers are referred to Ref. 30 for an overview of the algorithm. Let us however comment on some aspects of our implementation. The axial field does not break SU (2)-spin rotation symmetry and we have found it important to impose this symmetry by opting for a discrete Hubbard-Stratonovitch transformation coupling to the charge: with cos(α) = e −∆τ U/2 . We have furthermore used a symmetric Trotter decomposition with ∆τ t = 0.1, and the trial wave function corresponds to the ground state of the non-interacting Hamiltonian. For this choice of the trial wave function projection parameters in the range Θt = 40 − 60 suffice to guarantee convergence to the ground state within the quoted accuracy.
Since the Monte Carlo simulations do not break SU (2) spin symmetry we have to rely on spin-spin correlations to detect the global edge-compensated antiferromagnetic state. In Fig. 5 (upper panel) we consider 600-site flakes at U/t = 2. This choice of U/t places us well below U c ≈ 3.8t at which the transition to the antiferromagnetic Mott insulating state occurs in the absence of axial field. 10 The reference site i i i R = (−1, −1/ √ 3)/2 is chosen to belong to the sub-lattice which hosts the (normalizable) zero-energy modes. In the very close vicinity of i i i R antiferromagnetic correlations are apparent and they rapidly give way to dominant ferromagnetic correlations. For the uniform, χ(R) = qR 2 , axial field the extent of the dominant ferromagnetic correlations is considerably larger than for the localized one, χ(R) = q log(R). Finally strong antiferromagnetic correlations are present at the edge of the flake. Hence the overall features present in the quantum Monte Carlo calculations support the mean-field picture of the global edge-compensated antiferromagnetic spin structure.
A more precise measure for the global edgecompensated antiferromagnetic state can be obtained by considering In the absence of axial field (q = 0) only short range antiferromagnetic correlations are present and S(R) quickly decays to zero. At finite values of the axial field S(R) images the global edge-compensated antiferromagnetic spin structure. For the localized field the zero energy modes are localized around the center and S(R) quickly approaches a plateau value before being compensated by the edge antiferromagnetic correlations. In contrast, for the uniform field S(R) builds up as a function of distance before being again compensated by the edge magnetism. In Fig. 6 we consider larger values of U/t = 4. In the absence of the axial field, this choice of the Hubbard interaction places us in the antiferromagnetic Mott insulting phase. As a consequence, and in comparison to the U/t = 2 case, the integrated spin-spin correlations of Eq. (12) show small fluctuations up to large distances. It is interesting to note that even starting from the Mott insulating state, the axial field leads to the same reorganization of the spin-spin correlations as observed at weak couplings. One will nevertheless observe substantial antiferromagnetic oscillations superimposed on the edge-compensated antiferromagnetic spin structure.

VII. SUMMARY AND DISCUSSION
To summarize, we proposed a specific modulation of the nearest-neighbor hopping amplitudes in honeycomb lattice that captures the coupling of the low-energy Dirac quasi-particles to the (time-reversal-symmetric) axial magnetic fields. Due to the presence of the axial magnetic field a finite number of states appears at (near) zero energy, which in turn enhances the effect of electron-electron interaction. Various orderings can take place this way in strained graphene even at weak interactions. 14-17 .
In this paper we considered only the onsite Hubbard interaction between the fermions, and studied the nature of the magnetic ground state in strained graphene. Due to the special structure of the zero-energy states, which are supported by one sublattice in the bulk of the system, and by the other one near the boundary, the magnetic ground state in strained graphene lacks any analog in pristine graphene, or in graphene in true magnetic fields. Through the numerical self-consistent Hartree calcula-tion, and a separate quantum Monte Carlo simulation, we established that the magnetic ground state gives rise to both antiferromagnetic and ferromagnetic orders, locally, everywhere in the system. Although the antiferromagnetic order parameter is of the same sign in the entire system, the magnetization changes its sign near the boundary, so that the total magnetization is actually zero. Such ordering takes place even for weak Hubbard interactions. We named the magnetic ground state in strained graphene edge-compensated antiferromagnet.
In contrast, the ground state of the Hubbard model on honeycomb lattice subject to true magnetic fields (and with the Zeeman coupling ignored 28 ) is the conventional Néel antiferromagnet. Through the self-consistent calculation we show that such Néel ordering takes place again for weak interactions, and the order parameter closely resembles the profile of the true magnetic fields. The magnetization is in this case, however, identically zero everywhere in the system. the grant number AS120/9-1 and AS120/10-1 (Forschergruppe FOR 1807). IFH has been supported by the NSERC of Canada. We thank the Jülich Supercomputing Centre and the Leibniz-Rechenzentrum in Münich for generous allocation of CPU time.