Competing charge-density wave, magnetic, and topological ground states at and near Dirac points in graphene in axial magnetic fields

In the presence of axial magnetic fields that can be realized in deliberately buckled monolayer graphene, quasirelativistic Dirac fermions may find themselves in a variety of broken-symmetry phases even for weak repulsive interactions. Through a detailed Hartree self-consistent numerical calculation in finite strained graphene with cylindrical and open boundaries, we establish the possibility of realizing a charge-density wave order for the spinless fermions in the presence of weak nearest-neighbor repulsion. Such an instability gives rise to a staggered pattern of average fermionic density between bulk and boundary of the system as well as among two sublattices of honeycomb lattice, due to the spatial separation of the zero-energy states localized on opposite sublattices. Although with fermions spin restored, an unconventional magnetic order driven by the onsite repulsion possibly leads to the dominant instability at the Dirac point, the proposed charge-density wave order can nevertheless be realized at finite doping, which is always accompanied by a finite ferromagnetic moment. Additionally, the charge-density wave phase supports a quantized charge or spin Hall conductivity when its formation away from the Dirac point is further preceded by the appearance of topological anomalous or spin Hall insulator, respectively. The topological orders in strained graphene can be supported by weak second-neighbor repulsion, for example. Therefore, depending on the relative strength of various short-range components of the Coulomb interaction, several broken-symmetry phases can be realized within the zero-energy manifold in strained graphene.

Figure 1

Proposed modulation of the NN hopping amplitude in honeycomb lattice with cylindrical boundary, yielding a finite axial magnetic field. Red (black) corresponds to stronger (weaker) bond. Along the blue bonds, hopping is unchanged ($t=1$). Red/black bonds experience a gradient along the axis of the cylinder (not shown explicitly). $A$ and $B$ are the two sublattices of honeycomb lattice, and $j$ represents the row index.

Figure 2

Plots of LDOS in noninteracting system on sites of $A$ (red) and $B$ (black) sublattices, at the bottom (top, left), near the center (top, right), at the center (bottom, left), and top (bottom, right) row of the cylinder, for $q=0.02$. $E$ is measured in the units of $t(=1)$. We average the density of states over an interval $\Delta E=0.1$, yielding the broadening in LDOS.

Figure 3

Plots of LDOS on $A$ (red), and $B$ (black) sublattices for $V_1=0.2$ and $q=0.02$. Choices of the location of sites on two sublattices in different plots are same as in Fig. 2. Here as well $E$ is measured in units of $t$, and $\Delta E=0.1$.

Figure 4

Spatial variation of the local CDW order ${\delta }_R$, defined in Eq. (10), for $q=0.02$, and $V_1=0.2$ (blue), $0.4$ (red), $0.6$ (black). Inset: plot of the same quantity, but for $V_1=1.0$ (brown), $1.25$ (red), $1.5$ (black), $2.0$ (blue), $2.5$ (green). Here, $j$ is the row index of the cylinder shown in Fig. 1.

Figure 5

Left: spatial variation of ${\delta }_R\left(r\right)$ in a quasicircular honeycomb lattice with open boundary for $q=0.0125$. Strength of the NN interaction reads as $V_1=0.2$, $0.3$, $0.4$, $0.5$, $0.6$, $0.76$, $0.9$, from bottom to top. Right: scaling of the average CDW order in the bulk of the system, $\langle {\delta }_R\rangle ={\sum }_{r=1}^6{\delta }_R\left(r\right)/6$, with the axial magnetic field ($b$). The strength of the NN interaction reads as $V_1=0.76$, $0.6$, $0.5$, $0.4$, $0.3$, $0.2$ from top to bottom. The bottom (top) solid line shows almost linear ($\sqrt{b}$) dependence of the gap with $b$ for $V_1=0.2[0.76={\left(V_1\right)}_c]$. Here, $b$ is measured in units of $b_0\sim {10}^4$ T, an axial field associated with the lattice spacing $c\sim 3\text{Å}$.

Figure 6

Left: finite-size scaling of the local CDW order parameter near the center of the cylinder (circles), at the top edge (triangles), for $V_1=0.2$ (green), $0.4$ (black), and $q=0.02$. Here, $j_{\max }$ corresponds to the maximal number of row in cylindrical graphene system. Right: finite-size scaling of the same quantity in honeycomb lattice with open boundary, at the center (red) and middle (black) of the system for $V_1=0.3$ (circles), $0.4$ (triangles), and $q=0.015$.

Figure 7

Top: splitting of the zero-energy states in strained graphene (for spinless fermions), when (a) $V_1\gg V_2$, (b) $V_1\sim V_2$, and (c) $V_2\gg V_1$. States localized on the $A$ ($B$) sublattice are shown in blue (red), and $\pm$ correspond to their valley index $(\pm \vec{Q})$. Bottom: a representative mean-field phase diagram of the $V_1-V_2$ model in strained graphene, when the chemical potential is pinned at the charge-neutrality point.

Figure 8

Interaction-driven splittings of the topological zero-energy manifold in strained graphene. Here, blue and red correspond to the zero-energy states, localized on $A$ and $B$ sublattices, respectively, and $\pm$ to the ones localized near the Dirac valleys at $\pm \vec{Q}$. Projections of electron's spin along the $z$ axis are represented by $\uparrow ,\downarrow$. (a) Leading instability at the Dirac point driven by global antiferromagnetic order. (b) Subleading splittings driven by the (i) QSHI, (ii) QAHI, and (iii) CDW, when the chemical potential $\mu \to {\mu }_1$. Splittings shown in panel (c) are followed by the ones shown in panel (b), and all the degeneracies of the zero-energy subband are removed as $\mu \to {\mu }_2$.