# Inhomogeneous magnetic catalysis on graphene’s honeycomb lattice

###### Abstract

We investigate the ordering instability of interacting (and for simplicity, spinless) fermions on graphene’s honeycomb lattice by numerically computing the Hartree self-consistent solution for the charge-density-wave order parameter in presence of both uniform and non-uniform magnetic fields. For a uniform field the overall behavior of the order parameter is found to be in accord with the continuum theory. In the inhomogeneous case, the spatial profile of the order parameter resembles qualitatively the form of the magnetic field itself, at least when the interaction is not overly strong. We find that right at the zero-field critical point of the infinite system the local order parameter scales as the square-root of the local strength of the magnetic field, apparently independently of the assumed field’s profile. The finite size effects on various parameters of interest, such as the critical interaction and the universal amplitude ratio of the interaction-induced gap to the Landau level energy at criticality are also addressed.

## I Introduction

Graphene has since its successful fabrication Geim () emerged as the prime electronic system of reduced dimensionality. Its structure can be described as two interpenetrating
triangular sublattices of carbon atoms, which together form a bipartite honeycomb lattice.
As a consequence of the lack of inversion symmetry around a site of honeycomb lattice, the valence
and the conduction bands touch each other at the six corners of the first Brillouin zone. At low
energies, one can linearize the electronic dispersion relation near those “Dirac” points. In
particular, at the filling one half when the conduction band is empty and the valence band is
filled, gapless quasiparticle excitations live in the vicinity of the Dirac points. In the continuum
limit, such excitations can then be described in terms of pseudo-relativistic massless Dirac fermions,
with the Fermi velocity playing the role of the velocity of light () castro ().

In its usual state graphene behaves like a semi-metal. A large overlap of the electronic wave
functions of the neighboring carbon atoms ( eV) protects such a phase against weak
electron-electron interactions. In the language of renormalization group, such a stability corresponds
to a large domain of attraction of the non-interacting Gaussian fixed point Igor (). On the other hand,
a strong enough interaction can bring on a Mott-Hubbard transition towards a gapped insulating phase
gonzales (); bitan (). For example, a sufficiently strong on-site Hubbard interaction or nearest-neighbor
Coulomb repulsion would turn the system into an insulator with the staggered pattern of either
average magnetization or density. For graphene Wehling (), eV, eV, whereas the critical values for insulation are Sorella (); Paiva () and Franz (). It appears that graphene lies safely on the semi-metallic side of possible
Mott transitions, but with the interactions which are nevertheless not too far from their critical values.
Currently, the interaction-to-bandwidth ratios that control the Mott transitions
in graphene are not easily tunable. However, subjecting the system to a finite magnetic flux quenches the kinetic energy and collapses the density
of states (DOS) onto a discrete set of Landau levels (LLs), and can this way “catalyze” the formation of ordered phases gusynin1 (). In presence of a magnetic field even an infinitesimal amount of the on-site Hubbard or the nearest-neighbor repulsion would turn the system at half filling into an ordered phase with either finite
Nel order or staggered density dima (); herbut4 (). Yet another and a qualitatively different insulator haldane () may result from a strong next-nearest-neighbor repulsion, which can induce a gapped phase with finite circulating currents between the sites on the same sublattice raghu (). Such a phase, in the spinless case, violates the time-reversal symmetry and represents an early example of a topological insulator. The same topological insulator may also be possible to catalyze by a fictitious magnetic field herbut3 () that would arise from specific deformations of the graphene sheet guinea (); levy (), for example.

The magnetic catalysis in the presence of uniform magnetic field is by now well understood. It has been proposed as a mechanism behind the formation of the Hall states in graphene at filling factors and herbut4 (); gusynin (), which become discernible at higher magnetic fields. Zhang () Sublinear scaling of the gap with the magnetic field, for example, strongly suggests that the electron-electron interactions are the cause of the gap at Zhang2 (); bitan2 (). In contrast, the behavior of interacting electrons in presence of an inhomogeneous magnetic field has not been studied much, although the issue of order parameter’s dependence on the local value of the magnetic field has been addressed analytically, for specific spatial profiles of the field Dunne (); Raya (). In this work we attempt to develop a more detailed understanding of the spatial variation and the field dependence of the order parameter when an inhomogeneous magnetic field penetrates through the system. The motivation for such a study comes in part from a closely related problem of interacting electrons in a pseudomagnetic field herbut3 (); guinea (), where field’s profile is typically non-uniform in space. On a methodological level, it seems also interesting to inquire how much of the catalysis mechanism remains in effect when the condition of uniformity of the magnetic field is relaxed.

A self-consistent profile of the local gap in the insulating phase is computed therefore numerically on a discrete lattice and at the level of Hartree approximation, with a special attention given to its spatial variation. We find that the system stills suffers a metal-insulator transition at weak nearest-neighbor interactions even in the presence of a inhomogeneous magnetic field. In the continuum, this phenomenon would be attributed to the delta-function density of states when the Fermi energy is at the Dirac point casher (). Interestingly, we find that the spatial profile of the interaction-induced gap (order parameter) in presence of a localized magnetic flux, although not matching exactly, still mimics closely the profile of the local strength of the magnetic field. Moreover, right at the zero-field metal-insulator quantum criticality, the local order parameter seems to vary very much like the square-root of the local magnetic field. This behavior is analogous to what we previously found analytically in a uniform magnetic field bitan2 (), and to what we also confirm here numerically (see below). Away from the critical point, at weak interactions the expectation value of the local order parameter reverts to a linear dependence of the local magnetic field, as one might expect.

For a uniform magnetic field, we in general find a very good agreement between the previous field theoretic results and our numerical calculations. We focus on the finite-ranged components of the Coulomb repulsion, and chose to keep only the simplest one, which acts between the nearest neighbors. The system at the filling one-half and in the magnetic field then develops a gap in the spectrum, even when the interaction is weak. Right at the metal-insulator quantum critical point the gap behaves as

(1) |

where is the LL energy. Here is the universal number, found to be in the largest
system considered here, within the the Hartree approximation. This is in satisfactory agreement with
the same quantity computed previously in the field-theoretic description bitan2 (), where we found it to
be , in the limit of infinite number of fermion components. These two procedures being equivalent, we indeed find that upon increasing system’s size the constant slowly approaches its value in the continuum.

The rest of the discussion is organized as follows. The system of free electrons on honeycomb lattice in presence of magnetic field is introduced in Sec. II. The Hartree mean-field theory of the electrons interacting via the
nearest-neighbor repulsion is outlined in Sec. III. In that section we also demonstrate the mechanism of magnetic catalysis in a finite size system. Sec. IV focuses on the scaling behavior of the gap with the strength of the
uniform magnetic field and the interaction. Sec. V is devoted to the spatial variation of the gap in presence of nonuniform magnetic field. In Sec. VI, we summarize the results and discuss some related issues.

## Ii Free fermions

Let us define the system of non-interacting electrons on honeycomb lattice in the presence of a uniform magnetic field. The tight-binding model with nearest-neighbor hopping is defined as

(2) |

where and are the usual fermionic annihilation and creation operators, respectively. Here we
omitted the spin degrees of freedom, for simplicity. denotes the sublattice generated by the linear
combination of basis vectors and , for example. The second sublattice is then at , with being either or , where is the lattice spacing. The magnetic
field may be introduced through the Peierls substitution , where is the usual flux-quantum, and counts the magnetic flux through each plaquette of the honeycomb lattice. In case of graphene, corresponds to a magnetic field T, with commonly assumed lattice constant,
. Such a high magnetic field corresponds to the magnetic length close to the lattice scale, . Therefore is equivalent to , where corresponds to T.
In Fig. [1] we have shown one way to introduce a uniform magnetic field on honeycomb lattice. By solving numerically the tight-binding model on a lattice with periodic boundary in x-direction and the field T, we clearly see the first few LLs as the well-separated energies where the
DOS is sharply peaked (black curve Fig. [2]). The energy spectrum is symmetric about
zero and the spacing among the LLs decreases with the LL index. The energy of the LLs varies as the square root of the
magnetic field, due to the relativistic nature of the quasi-particles (top curve Fig. [3]). We also found that the maximum energy of the free electron system is , in agreement with the previous results kohmoto (). It may be worth mentioning that in presence of the periodic boundary conditions the choice of gauge requires some care, and is chosen here so that only one out of the three bonds emanating from a site contributes to it. Such a choice is then equivalent to the Landau gauge in the continuum description.

Next, we consider still non-interacting electrons on the honeycomb lattice, but now subject to an inhomogeneous magnetic field. Numerically diagonalizing the free electron Hamiltonian, the DOS is found to be a smooth
function of energy, and peaked only at the zero energy (red curve in Fig. [2]). There we considered a lattice, with a open boundary and total flux . We assumed the field to be uniform in the x-direction, and bell-shaped in y-direction, with the maximum at the center. The number of near zero energy states is proportional to the total flux of the magnetic field enclosed by the system.casher () The maximum energy in the free electron spectrum is found to be .

## Iii Interactions and magnetic catalysis

Next, we turn on the short-range electron-electron interaction. The Hamiltonian in the presence of only the nearest-neighbor Coulomb repulsion is given by

(3) |

where stands for the summation over the nearest-neighbor sites, is the total number of electrons and is the chemical potential. After the usual Hartree decomposition the effective single-particle Hamiltonian for interacting electrons becomes

(4) | |||||

where counts the self-consistent site-dependent average electron density on sub-lattice B(A). Let us measure these relative to the uniform density at half-filling by defining

(5) |

The positive quantities determine the local charge-density-wave order parameter (OP). Both and will be functions of position with the constraint that the system is precisely at half filling,

(6) |

We also choose the value of .

We have computed the (Hartree) self-consistent solutions for the OPs, for different values of the flux and for a variety of interaction strengths (), at . Consider a lattice with a periodic boundary in the direction and let us conveniently define the local OP as,

(7) |

where is either one of the two nearest-neighbors to the site , on the same row in x-direction. this way measures the order parameter
in a unit cell. On the other hand, the OP will be averaged over the points connected by the symmetry,
when we considered a quasi-circular system with an open boundary. In the presence of a uniform magnetic field,
is found to be uniform in the bulk of the system. However, becomes position dependent and
proportional to the local field when the system is subject to a inhomogeneous field.

Before we proceed, it is worth pausing to establish a practical definition of the “critical interaction” associated with the metal-insulator transition in the finite size system like ours. For a sufficiently strong nearest-neighbor interaction, fermions reside only on one sublattice, with the other one completely empty. The system is then deep in the insulating phase. As the interaction is weakened, the size of the order parameter decreases and we numerically find the system to go through a well-defined transition into the semi-metallic phase, where the
order parameters are zero in the entire system. Right above that particular interaction which we call critical, there is a finite, but slightly inhomogeneous staggered density everywhere in the system. We will designate that value as the critical interaction corresponding to the metal-insulator transition. The described scenario is quite generic and occurs both in the absence and presence of magnetic fields, which also may be either uniform or nonuniform. The observed non-analytic behavior in a finite system, however, is clearly an artifact of the Hartree approximation, i. e. of the requirement of the self-consistency of the solution.

We computed the variation of the critical interaction defined this way
with the lattice size and the geometry, for a particular
magnetic flux () through each plaquette of the lattice. As may be seen from the Fig. [4], for a small system size the critical interaction is large even in the presence
of a magnetic field, and also depends on the geometry of the lattice, as well as on the choice of the gauge. Upon
increasing the size of the system, the value of the critical interaction decreases and appears to approach zero
in the thermodynamic limit, in agreement with the results obtained in the continuum theory. A typical distribution of
the OP in a lattice with periodic boundary is shown in inset of Fig. [4]. From that one can
conclude that upon decreasing , the size of the order parameter decreases in the entire system, both
in the bulk and at the edge. However a finite gap in the spectrum exists even at a rather weak interaction, (bottom curve). Hence, in the presence of a uniform magnetic field, we expect that a large system would find itself in a gapped insulating phase even at an infinitesimal interaction, as found in the continuum theory.

On the other hand, when the system is exposed to a localized flux of magnetic field
OP develops a local expectation value (see Fig. [5]). The local OP is found to be proportional to the local magnetic field. As one enters the regime of weaker interaction OP decreases both in the bulk and the edge of the system. Yet we managed to observe finite expectation value of the OP in the entire system even at the smallest interaction considered here, . Therefore, the system can also find itself in a ordered phase at weak interaction when an inhomogeneous flux of the magnetic field pierces through it. This phenomenon can be attributed to the finite density of state at zero energy, where the chemical potential lies at the filling one-half.

Besides the condensation in the bulk of the system, the OP acquires spikes near
the edges of the system, in presence of both uniform (Inset of Fig. [4]) and non-uniform (Fig. [5]) field. The spikes in the OP near the edge of the system arise from the finite size effect. Such edge effects die out as one increases
the system size. Moreover, with the increasing magnetic field, those spikes also dissolve and give rise to a uniform condensation throughout the system, at sufficiently large magnetic field. These effects on the OP are demonstrated in Fig. [6]. We exhibited the finite-size effects in presence of a uniform flux only, but the result is qualitatively the same in presence of a localized flux as well.

## Iv Scaling in uniform magnetic field

We now investigate the dependence of the gap on the magnetic field and the interaction .
First we take the field to be uniform, and still consider the lattice with the periodic boundary in the direction.
The functional dependence of , defined in Eq. 7, on the magnetic fields and interactions is shown in Fig. [7]. Here, is computed on a lattice and we considered OP only in the bulk of the system. With the parametrization as in Sec. II, the lowest value of the magnetic field ( T), considered here, is about four times larger than the current highest constant laboratory magnetic field. However, upon using a larger system one can get down to a more realistic strength of the field.

For sufficiently small interactions, order parameter ( or ) varies almost linearly with (bottom
curve in Fig. [7]). As one increases the strength of interaction, there is a crossover to a sub-linear dependence
of the mass () on bitan2 (). In particular, right at , we find the best overall fit of the mass to the magnetic field. We therefore designate that interaction to be the critical interaction at . When we computed the ratio of the first LL energy to the interaction induced gap at , it came out to be a universal number , independent of the magnetic field (inset Fig. [3]). The value of the number is in satisfactory agreement with the same quantity previously calculated in the continuum description and in the large-N limit bitan2 (). There we obtained , with the difference between the two values that can be attributed to the finite size (Inset Fig. [3]).

We also found a similar dependence of the mass on interactions and magnetic fields, computed on a quasi-circular lattice with open boundary, preserving the symmetry of a hexagon. A spatial variation of the interaction induced gap in the presence of a uniform magnetic field is shown in Fig. [8] (black curves). In that case, we found the
ratio of the first LL energy to the interaction induced gap to be , in a system of lattice points. Such a particular choice of lattice turns out to be useful when one imposes rotationally symmetric inhomogeneous magnetic field. By considering a graphene sheet with open boundary, we computed the OP in two different gauges, equivalent to
and and it turned out to be gauge independent, as expected.

## V Interacting fermions in inhomogeneous field

Next, we consider spinless interacting fermions on a honeycomb lattice subject to a inhomogeneous magnetic field in more detail. It was previously shown that, for a specific realization of the inhomogeneous magnetic field and in the limit of a large magnetic flux, the order parameter in the insulating phase computed within the zero-energy manifold matches exactly the local profile of the magnetic field Dunne (). Here we determine the order parameter self-consistently (at ) and on the honeycomb lattice, and include the contributions from all the states into account. We will consider two specific configurations of spatially modulated magnetic field. Localized field in one direction,
in our case, but extended in the orthogonal direction, and rotationally symmetric localized field with the maximum strength at the center. We imposed the field of type on a lattice with periodic boundary in direction. On
the other hand, a quasi-circular lattice with open boundary, preserving the symmetry of a hexagon, is exposed to a localized field of type . As mentioned previously, even in presence of a non-uniform field, there is a large (and in the continuum limit, infinite) DOS at zero energy. Therefore, right at the filling one half, one expects that even a weak interaction can place the system into a insulating phase.

In the presence of a non-uniform field, but with a finite total magnetic flux, the system develops a gap in the spectrum even at sub-critical interactions () (Fig. [5]). The spatially variation of the order parameter in presence
of a inhomogeneous magnetic field of type is depicted in Fig. [9]. We considered the OP
only far from the edges of the system, and normalized the OP as well the magnetic fields with respect to their
maximum values at the center of the system. The order parameter appears to follow the spatial profile of the magnetic field, and to depend on its local strength for . Near the zero-field criticality the profile of the OP follows the magnetic field’s more closely, whereas at large interactions the effect of inhomogeneous magnetic field becomes irrelevant, leading to uniform condensation. In Fig. [8] (red curves), we exhibited the spatial distribution of the OP in presence of a non-uniform magnetic field, applied on a lattice that preserves the symmetry of a hexagon. Our computation yields an interaction induced OP as a function of space qualitatively similar to the assumed profile of the magnetic field itself.

Let us now turn to functional dependence of the OP on the magnetic field and interaction, when the former is
space dependent. At sufficiently weak interactions, the size of the gap at different region of the bulk of the system varies almost precisely linearly with the local strength of magnetic field. The linear dependence of the local OP at with the local magnetic field is shown in Fig. [10]. As the interaction is increased there is a
crossover to a sub-linear dependence of the mass on the local magnetic field. Situation is quite similar to the one
in presence of a uniform field. Right at the zero-field criticality , the local OPs in the entire bulk of the system varies as , independent of the position (blue dots in Fig. [3]). This suggests that the OP in the insulating phase may be a universal function of the local magnetic field, independent of its spatial distribution.

## Vi Summary and Discussion

In the present work we systematically studied magnetic catalysis for the spinless interacting electrons on a honeycomb lattice of finite extension, both for uniform and spatially modulated magnetic fields. In presence of the magnetic field, either uniform or non-uniform, the semimetal-insulator transition takes place at weak interaction in a large system. We here considered only the nearest-neighbor component () of the Coulomb interaction, and omitted its long-ranged tail gonzales (); semenoff () for simplicity. We computed the self-consistent Hartree solution of the interaction-induced gap while keeping the system at the filling one half and presented a scaling behavior of the interaction-induced order parameter (or a gap ) with the magnetic field and interaction, at . At weak interaction we observed a linear variation of the interaction induced local OP with the local magnetic field. With increase in the strength of the interaction we find a crossover from linear to a sub-linear dependence of the mass on the local magnetic field. A perfect dependence of the OP emerges when the system is tuned to be precisely at the zero-field criticality, which we identified to be at . This is close to the value found analytically Franz ().

In our analysis we have considered only the nearest-neighbor hopping amplitude , while neglecting the next-nearest-neighbor hopping , which in graphene, for example, is finite but rather small. The main effect of a finite is the violation of the perfect particle-hole symmetry of the free electron spectrum. On the basis of continuum theory we expect, however, that the inclusion of would not change our results in a significant way, once the chemical potential is adjusted so that the central (formerly zero-energy) LL is half filled. A more detailed analysis is left for future study.

If we were to restore the spin of electrons, we would need to include a finite on-site Hubbard interaction as well. In absence of a magnetic field, anti-ferromagnetic (AF) ground state, is energetically favored for a large on-site Hubbard interaction when the chemical potential is at the Dirac point bitan (). The presence of magnetic field stabilizes such ground state even at an infinitesimal on-site interaction herbut4 (); Meng (). We therefore expect that the system would develop a local expectation value of the Nel order parameter when in a inhomogeneous magnetic field, if . If , on the other hand, the system would decrease the energy more by forming a charge density wave order of the type we considered here. At zero magnetic field the two quantum phase transitions belong to distinct Gross-Neveu universality classes vladimir (); subir (). A scaling behavior in a similar model has also been studied recently in presence of both real and pseudo magnetic fields.bitanown ()

## Vii Acknowledgement

This work is supported by the NSERC of Canada. B. Roy thanks Peter Smith and Kamran Kaveh for useful discussions, and to Vladimir Juričić, Kelly Cheng and Payam Mousavi for their help with the manuscript.

## References

- (1) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, A. A. Firsov, Science 306, 666 (2004).
- (2) For a review, see A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
- (3) I. F. Herbut, Physics 2, 57 (2009).
- (4) J. Gonzales, F. Guinea, M. A. H. Vozmediano, Nucl. Phys. B424, 595 (1994); Phys. Rev. B 59, 2474 (1999).
- (5) I. F. Herbut, Phys. Rev. Lett. 97, 146401 (2006); I. F. Herbut, V. Juričić, B. Roy, Phys. Rev. B 79, 085116 (2009), and references therein.
- (6) T. O. Wehling, E. Sasioglu, C. Friedrich, A. I. Lichtenstein, M. I. Katsnelson, S. Blügel, arXiv:1101.4007.
- (7) S. Sorella and E. Tosatti, Europhys. Lett. 19, 699 (1992).
- (8) T. Paiva, R.T. Scallettar, W. Zheng, R. R. P. Singh, J. Oitmaa, Phys. Rev. B 72, 085123 (2005).
- (9) C. Weeks and M. Franz, Phys. Rev. B 81, 085105 (2010).
- (10) V. P. Gusynin, V. A. Miransky, and I. A. Shovkovy, Phys. Rev. Lett. 73, 3499 (1994); Phys. Rev. D 52, 4718 (1995).
- (11) D. V. Khveshchenko, Phys. Rev. Lett. 87, 246802 (2001); H. Leal and D. V. Khveshchenko, Nucl. Phys. B B687, 323 (2004).
- (12) I. F. Herbut, Phys. Rev. B 75, 165411 (2007); Phys. Rev. B 76, 085432 (2007).
- (13) F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).
- (14) S. Raghu, Xiao-Liang Qi, C. Honerkamp, S.-C. Zhang, Phys. Rev. Lett. 100, 156401 (2008).
- (15) I. F. Herbut, Phys. Rev. B 78, 205433 (2008).
- (16) F. Guinea, M. I. Katsnelson, A. K. Geim, Nat. Phys. 6, 30 (2010).
- (17) N. Levy, S. A. Burke, K. L. Meaker, M. Panlasigui, A. Zettl, F. Guinea, A. H. Castro Neto, M. F. Crommie, Science 329, 544 (2010).
- (18) V. P. Gusynin, V. A. Miransky, S. G. Sharapov, and I. A. Shovkovy, Phys. Rev. B 74, 195429 (2006); E. V. Gorbar, V. P. Gusynin, V. A. Miransky, I. A. Shovkovy, Phys. Rev. B 78 085437 (2008).
- (19) Y. Zhang, Z. Jiang, J. P. Small, M. S. Purewal, Y. W. Tan, M. Fazlollahi, J. D. Chudow, J. A. Jaszczak, H. L. Stormer and P. Kim, Phys. Rev. Lett. 96, 136806 (2006).
- (20) Z. Jiang, Y. Zhang, H. L. Stormer, and P. Kim, Phys. Rev. Lett. 99, 106802 (2007).
- (21) I. F. Herbut, B. Roy, Phys. Rev. B 77, 245438 (2008).
- (22) G. Dunne and T. Hall, Phys. Rev. D 53, 2220 (1996).
- (23) A. Raya and E. Reyes, Phys. Rev. D 82, 016004 (2010).
- (24) Y. Aharonov and A. Casher, Phys. Rev. A 19, 2461 (1979).
- (25) Y. Hasegawa, M. Kohmoto, Phys. Rev. B 74, 155415 (2006); D. R. Hofstadter, Phys. Rev. B. 14, 2239 (1976).
- (26) V. Juričić, I. F. Herbut, G. W. Semenoff, Phys. Rev. B, 80, 081405 (R) (2009); V. Juričić, O. Vafek, I. F. Herbut, Phys. Rev. B 82, 235402 (2010), and references therein.
- (27) See however, Z. Y. Meng, T. C. Lang, S. Wessel, F. F. Assaad and A. Muramatsu, Nature 464, 847 (2010) about the possible appearance of a spin-liquid phase at an intermediate strength of Hubbard .
- (28) I. F. Herbut, V. Juričić, O. Vafek, Phys. Rev. B 80, 075432 (2009).
- (29) S. Sachdev, arXiv:1012.0299.
- (30) B. Roy, arxiv:1012.2109.