Chiral symmetry and its manifestation in optical responses in graphene: interaction and multilayers

As to the electron–electron interaction, we take the two-body interaction V_ijn_in_j between sites 〈ij〉, which can include long-range interactions. Note that

$$\begin{eqnarray} &&n_{i } n_{j } = c_{i } ^\dagger c_{i } c_{j } ^\dagger c_{j } = c_{i } ^\dagger c_{j } ^\dagger c_{j }c_{i } = -1+ n_{i } +n_{j } + c_{i }c_{j } c_{j } ^\dagger c_{i } ^\dagger . \end{eqnarray} \tag{ 21 }$$

We can express the interaction in an electron–hole symmetric form, at half-filling, as

$$\begin{eqnarray} &&\fl {\cal H}_{\rm int} = \sum_{ i j \rangle} V_{ij} \left( n_i- \frac {1}{2} \right) \left( n_j- \frac {1}{2} \right) = \frac {1}{2} \sum_{\langle ij \rangle} V_{ij} \bigg[ c_i ^\dagger c_j ^\dagger c_j c_i + (c\rightleftarrows c ^\dagger ) \bigg]+{\rm const}. \end{eqnarray} \tag{ 22 }$$

We assume that the electron–electron interaction is sufficiently weak as compared with the Landau gap between the n = 0 Landau level and those for n = ± 1. When the filling factor of the n = 0 LL is ν < 1, the ground state is then given by the configurations of the n = 0 LL states, while the filled states below them can be regarded as a 'Dirac sea'. Since the one-particle spectrum in graphene is given by the -zero modes with chiral pairs as discussed above, the kinetic energy from the -zero modes is negligible, so that the ground state is given perturbatively (i.e. for the interaction energy that is smaller than the Landau gap) as

$$\begin{eqnarray} &&|\Psi \rangle = \sum_{i_1,\ldots, i_M\subset\{1,\ldots, M\}} C_{i_1,\ldots, i_M} d_{i_1} ^\dagger \cdots d_{i_M} ^\dagger | D_< \rangle ,\quad | D_< \rangle = \prod_{j=1}^{M'} {d_<}_{j} ^\dagger | 0 \rangle , \end{eqnarray} \tag{ 23 }$$

where |D_<〉 is the filled Dirac sea (d_i|D_<〉 = 0 for all i's), {i₁,...,i_M} a subset of the -zero modes and M is the number of particles. The free-fermion Hamiltonian ${\cal H}_0$ is written as

$$\begin{eqnarray} &&{\cal H}_0 = c ^\dagger \psi_T \,{\rm diag}\, ({\cal E},- \Lambda ,\Lambda )\psi_T ^\dagger c = {d}^\dagger {\cal E} d -{d_<}^\dagger \Lambda d_< +{d_>}^\dagger \Lambda d_>, \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} &&\left[ \begin{array}{@{}c@{}} {d }\\ {d_<}\\ {d_>}\end{array}\right] =\psi _T ^\dagger c = \left[ \begin{array}{@{}c@{}} {\psi ^\dagger c }\\ {\varphi ^\dagger c }\\ {\varphi^\dagger \Gamma c }\end{array}\right], \quad c = \psi _T\left[ \begin{array}{@{}c@{}} {d }\\ {d_<}\\ {d_>}\end{array}\right] =(\psi ,\varphi,\Gamma \varphi ) \left[ \begin{array}{@{}c@{}} {d }\\ {d_<}\\ {d_>}\end{array}\right] , \end{eqnarray} \tag{ 25 }$$

where d^† = (d^†₁,...,d_M), d_<^† = (d_<1^†,...,d_<M'^†) and d_>^† = (d_>1^†,...,d_>M'^†) are, respectively, fermion creation operators for the -zero mode, the negative-energy states and positive ones with {d_ℓ,d^†_ℓ'} = δ_ℓℓ', {d_<ℓ,d^†_<ℓ'} = δ_ℓℓ', {d_>ℓ,d^†_>ℓ'} = δ_ℓℓ', etc.

Now let us take the chiral basis for the -zero modes to describe the interaction between the particles. We assume that the perturbative ground state within the configuration of the n = 0 Landau level is determined by the projected interaction $\widetilde {\cal H}_{\mathrm { int}}$. The projected Hamiltonian is defined as

$$\begin{eqnarray} &&\fl \widetilde{\cal H}_{\rm int} = \frac {1}{2} \sum_{ij}V_{ij} [ \tilde {c}_{i } ^\dagger \tilde {c}_{j } ^\dagger \tilde {c}_{j } \tilde {c}_{i } + (c\rightleftarrows c ^\dagger )] = \frac {1}{2} \sum_{ij}V_{ij} (d ^\dagger \psi^\dagger )_i (d ^\dagger \psi^\dagger )_j (\psi d)_j (\psi d)_i + {\rm C.c.}, \end{eqnarray} \tag{ 26 }$$

where C.c. is a charge conjugation defined below, and $\tilde {c}_{i}$ is a projected fermion defined as

$$\begin{eqnarray} &&\tilde{c} \equiv \Psi _T\left[ \begin{array}{@{}c@{}} {d }\\ {0}\\ {0}\end{array}\right] =(\psi ,\varphi,\Gamma \varphi ) \left[ \begin{array}{@{}c@{}} {d }\\ {0}\\ {0}\end{array}\right]=\psi d =\psi \psi ^\dagger c = P c, \end{eqnarray} \tag{ 27 }$$

with P = ψψ^† = P^† being the projection to the -zero modes. These $\widetilde {c}_i$'s are not canonical fermion operators, since

$$\begin{eqnarray} &&\{{\tilde{c}}_i ,\tilde{c}_j ^\dagger \} = \{ (P )_{ii'}c_{i'} , c_{j'} ^\dagger (P^\dagger )_{j'j} \} =(P P ^\dagger )_{ij} =(P^2 )_{ij}=(P )_{ij}. \end{eqnarray} \tag{ 28 }$$

Note that the projected Hamiltonian is clearly positive semi-definite.

Let us take a chiral basis to define canonical fermions d_i for the -zero modes, defined as

$$\begin{eqnarray} &&{d} = \left[ \begin{array}{@{}c@{}} {{d}_+ }\\ {{d}_- }\end{array}\right]= \psi_\Gamma ^\dagger c =\left[ \begin{array}{@{}cc@{}} {\psi_\bullet ^\dagger } & {O }\\ {O } & {\psi_\circ ^\dagger } \end{array}\right]c = \left[ \begin{array}{@{}c@{}} {\psi_ \bullet ^\dagger c_\bullet }\\ {\psi_ \circ ^\dagger c_\circ }\end{array}\right], \end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray} &&{d}_+ = \left[ \begin{array}{@{}c@{}} {d_{1+}}\\ {\vdots}\\ {d_{M_++}}\end{array}\right], \quad {d}_-= \left[ \begin{array}{@{}c@{}} {d_{1-}}\\ {\vdots}\\ {d_{M_--}}\end{array}\right], \end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray} &&\tilde{c} = \left[ \begin{array}{@{\,}c@{}} { \widetilde {c} _\bullet }\\ { \widetilde {c} _\circ }\end{array}\right] = \left[ \begin{array}{@{}cc@{}} {\psi_\bullet } & {0 }\\ {0 } & {\psi_\circ }\end{array}\right] d = \left[ \begin{array}{@{}c@{}} {\psi_\bullet {d}_+}\\ {\psi_\circ {d}_-} \end{array}\right]. \end{eqnarray} \tag{ 31 }$$

We can define a chiral transformation ${\cal U}_\theta$ as

$$\begin{eqnarray} &&{\cal O} \mapsto {\cal U}_\theta ^\dagger {\cal O}{\cal U}_\theta, \quad {\cal U}_\theta = \mathrm{e} ^{\mathrm{i}\theta {\cal G}}, \end{eqnarray} \tag{ 32 }$$

where the generator of the transformation, the chirality ${\cal G}$, is given as

$$\begin{eqnarray} &&{\cal G} = \tilde{c}^\dagger \Gamma \tilde{c} = \sum_{i\in \bullet }\tilde {c} _i ^\dagger \tilde {c} _i - \sum_{i\in \circ }\tilde {c} _i ^\dagger \tilde {c} _i = {d}^\dagger \psi_\Gamma ^\dagger \Gamma \psi_\Gamma {d} = {d}^\dagger \Gamma _0 {d} = d_ + ^\dagger d_+-d_- ^\dagger d_-. \end{eqnarray} \tag{ 33 }$$

The chiral transformation operates as

$$\begin{eqnarray} &&d \mapsto {d}_\theta = {\cal U}_\theta ^\dagger \, {d} \, {\cal U}_\theta = {\rm e}^{{\rm i} \theta {\cal G}} \, {d} \, {\rm e}^{-{\rm i} \theta {\cal G}} = {\rm e}^{-{\rm i}\theta \Gamma _0} {d} =\left[ \begin{array}{@{}c@{}} {{\rm e}^{-{\rm i}\theta }{d}_+ }\\ {{\rm e}^{{\rm i}\theta}{d}_- }\end{array}\right] , \end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray} \tilde {c}_i \mapsto \left\{ \begin{array}{@{}cc@{}} {{\rm e}^{-{\rm i}\theta }\tilde {c}_i,} & {i\in \bullet, }\\ {{\rm e}^{+{\rm i}\theta }\tilde {c}_i, } & {i\in \circ. }\end{array}\right. \end{eqnarray} \tag{ 35 }$$

Then it is clear that the projected interaction is invariant under the chiral transformation as

$$\begin{eqnarray} &&\widetilde{\cal H}_{\rm int} \mapsto {\rm e}^{-{\rm i}\theta {\cal G}} \widetilde{\cal H}_{\rm int} \,{\rm e}^{+{\rm i}\theta {\cal G}} = \widetilde{\cal H}_{\rm int}. \end{eqnarray} \tag{ 36 }$$

For an infinitesimal transformation this implies a conservation law,

$$\begin{eqnarray} &&[{\cal G}, \widetilde{\cal H}_{\rm int}] = 0. \end{eqnarray} \tag{ 37 }$$

Let us next define a charge conjugation which is anti-unitary as

$$\begin{eqnarray} &&O\mapsto {\cal A}_C ^\dagger O {\cal A}_C, \end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray} &&{\cal A}_C = K \prod_{\ell=1}^M (d_{\ell} +d_{\ell} ^\dagger ) \prod_{\ell=1}^{M'} (d_{<\ell} +{d}_{<\ell} ^\dagger ) (d_{>\ell} +{d}_{>\ell} ^\dagger ), \end{eqnarray} \tag{ 39 }$$

where K is complex-conjugation with K² = 1. Its operation is, for example,

$$\begin{eqnarray} &&{\cal A}_C ^\dagger \, d_\ell \, {\cal A}_C =(-)^{N-1} d_\ell ^\dagger,\quad {\cal A}_C ^\dagger \, d_\ell ^\dagger \, {\cal A}_C =(-)^{N-1} d_\ell. \end{eqnarray} \tag{ 40 }$$

Then the invariance of the interaction Hamiltonian under the charge conjugation,

$$\begin{eqnarray} &&{\cal A}_C ^\dagger \widetilde{\cal H}_{\rm int} {\cal A}_C = \widetilde{\cal H}_{\rm int}, \end{eqnarray} \tag{ 41 }$$

trivially follows.

In this section, we restrict ourselves to considering a repulsive interaction that only acts between the • and ○ sites. The projected interaction in this case is written as

$$\begin{eqnarray} &&\fl \widetilde{\cal H}_{\rm int}= \frac {1}{2} \sum_{ i_\bullet j_\circ } V_{ i_\bullet j_\circ } ( \tilde {c}_{i \bullet} ^\dagger \tilde {c}_{j \circ} ^\dagger \tilde {c}_{j \circ} \tilde {c}_{i \bullet} + {\rm C.c.}) \end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray} &&\fl\hphantom{\widetilde{\cal H}_{\rm int}} = \frac {1}{2} \sum_{ i_\bullet j_\circ} V_{i_\bullet j_\circ } \bigg[ ( {d}_+ ^\dagger \psi_\bullet ^\dagger )_i ( {d}_- ^\dagger \psi_\circ ^\dagger )_j ( \psi_\circ {d}_- )_j ( \psi_\bullet {d}_+ )_i + ( \psi_\bullet {d}_+ )_i ( \psi_\circ {d}_- )_j ( {d}_- ^\dagger \psi_\circ ^\dagger )_j ( {d}_+ ^\dagger \psi_\bullet ^\dagger )_i \bigg], \end{eqnarray} \tag{ 43 }$$

where V_i•j○ > 0. Here we allow any (i.e. including long-ranged) repulsive interactions between the • and ○ sites.

Now we can readily see that the two states

$$\begin{eqnarray} &&|G _\pm \rangle = d_{1\pm} ^\dagger \cdots d_{M_\pm \pm} ^\dagger | D_< \rangle , \end{eqnarray} \tag{ 44 }$$

which have maximum and minimum chiralities, respectively, are special, since each of them has vanishing interaction energy with $\widetilde {\cal H}_{\mathrm { int}} |G_\pm \rangle = 0$. Namely, |G₊〉 has all the d_ℓ+ states filled with (d^†₊ψ^†_•)_i|G₊〉 = 0. Also, none of the d_ℓ− states are occupied, so that we have (ψ_○d₋)_j|G₊〉 = 0. Since the interaction Hamiltonian $\widetilde {\cal H}_{\mathrm { int}}$ is semi-positive definite, these two states are thus the ground states of $\widetilde {\cal H}_{\mathrm { int}}$ with a degeneracy of two.

As to the chiral transformation, we have

$$\begin{eqnarray} &&|G_\pm \rangle \mapsto {\cal U}_\theta ^\dagger |G_\pm \rangle = \mathrm{e}^{\mp \mathrm{i} M_\pm\theta} | G_\pm \rangle , \end{eqnarray} \tag{ 45 }$$

where ${\cal U}_\theta$ is the chiral operation defined above. This property follows, since |G_±〉 are eigenstates of the chirality ${\cal G}$ as

$$\begin{eqnarray} &&{\cal G} | G_\pm \rangle = \pm M_\pm | G_\pm \rangle. \end{eqnarray} \tag{ 46 }$$

Thus we have shown that the ground states are the chiral condensates, where the chirality is macroscopic (M_± ∼ N).

As to the charge conjugation, it operates as

$$\begin{eqnarray} &&| G_\pm \rangle \mapsto d_{1\mp} ^\dagger \cdots d_{M_\mp \mp} ^\dagger | D_> \rangle, \end{eqnarray} \tag{ 47 }$$

where |D_>〉 is the positive Dirac sea up to the sign.

The Hall conductivity of the doubly degenerate chiral condensates can be calculated with the Niu–Thouless–Wu formula as

$$\begin{eqnarray} &&\sigma _{xy} = \frac {e^2}{h} \frac {1}{N_{\rm D}} C, \quad C = \frac {1}{2\pi \mathrm{i}} \int {\rm Tr}_2\, F,\quad F = \mathrm{d}A + A^2,\quad A = \Psi ^\dagger \,\mathrm{d} \Psi, \end{eqnarray} \tag{ 48 }$$

where Ψ = (|G₊〉,|G₋〉) is the chiral-condensed ground states with the degeneracy of N_D = 2.

The doublet, being degenerate, can be mixed, but diagonalization of the chirality ${\cal G}$ acts to fix the gauge. Here we assume a finite energy gap above the ground state multiplet. Then the Chern number of the doublet becomes well defined, and is given by the sum as

$$\begin{eqnarray} &&C = C_++C_-, \end{eqnarray} \tag{ 49 }$$

where C_± is the Chern number for |G_±〉. Further, we can decompose the condensate as

$$\begin{eqnarray} &&C_\pm = C_{\psi_\pm }+C_{{\rm D}_<}, \end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray} &&C_{\psi_\pm } = \frac {1}{2\pi\mathrm{ i} } \int {\rm Tr}_{M_\pm}\, \mathrm{d} \psi_\pm ^\dagger \,\mathrm{d} \psi_\pm, \quad C_{{\rm D}_<}=\frac {1}{2\pi \mathrm{i} } \int \langle \mathrm{d} D_< | \mathrm{d} D_< \rangle , \end{eqnarray} \tag{ 51 }$$

where C_D< is the Chern number of the filled Dirac sea. Since the chiral operator Γ is unitary, the Chern number of the positive Dirac sea is the same as that of the negative Dirac sea as

$$\begin{eqnarray} &&C_{{\rm D}_>}=\frac {1}{2\pi \mathrm{i} } \int \langle\mathrm{ d} D_> |\mathrm{ d} D_> \rangle = C_{{\rm D}_<}\equiv C_{{\rm D}}. \end{eqnarray} \tag{ 52 }$$

The charge conjugation is anti-unitary, which implies that

$$\begin{eqnarray} &&C_{{\rm D}_>} + C_{\psi_-} = -(C_{{\rm D}_<} + C_{\psi_+} ), \end{eqnarray} \tag{ 53 }$$

and we end up with the total Chern number of the doubly degenerate chiral condensate as

$$\begin{eqnarray} &&C = \big(C_{{\rm D}_<} + C_{\psi_-} \big)+ \big(C_{{\rm D}_<} + C_{\psi_+} \big)= C_{\psi_-} + C_{\psi_+} + 2C_{\rm D} =0. \end{eqnarray} \tag{ 54 }$$

Namely, the exact chiral condensate for the half-filling (i.e. ν = 0) constructed from the chiral doublet for the projected interaction is a Hall insulator with the topological degeneracy of 2. Physical behaviors such as the boundary effects of the edge states and correlation functions are discussed in a recent paper [16].

The discussion here is for spinless fermions. An extension to chiral condensates for spinful cases can be important for the graphene if the state is spin unpolarized [13]. Also, optical responses of the chiral condensates are among the interesting future problems. These will be published elsewhere based on the present formulation.

Although the low-energy effective Hamiltonian has been discussed previously [17–19], the chiral symmetric model has several advantages. One is that the discussion of low-energy physics is much clearer since the derivation of the low-energy Hamiltonian is substantially simplified, especially as the behaviors of the wavefunctions are very clearly described. Also, we can treat the tight-binding Hamiltonian with a weak magnetic field directly since the zero energy is special for the chiral symmetric model and the physics near the zero energy, which is our main interest, is accessible by the sparse matrix technique as described below. This is essential to describe the physics of multilayers such as trigonal warping which is only relevant for a very weak magnetic field. The tight-binding model with such a very weak magnetic field is only accessible by using the chiral symmetry and the sparse matrix technique.

Let us start with the Hamiltonian of bilayer graphene in a magnetic field as a simple generalization of the monolayer model [6]. With j = (j₁,j₂) denoting two-dimensional coordinates, e_1,2 the unit translations and φ the total flux per hexagon in units of the flux quantum φ₀ = h/e (see figure 1), the Hamiltonian is written as

$$\begin{eqnarray} &&\fl {\cal H}^{\rm bilayer} = H_d + H_u + H_{ud}, \end{eqnarray} \tag{ 55 }$$

$$\begin{eqnarray} &&\fl {\cal H}_d = \gamma _0\sum_j\bigg[ d_\bullet ^\dagger ( j ) d_\circ (j )+ {\rm e}^{{\rm i}2\pi \varphi j_1 } d_\bullet ^\dagger (j ) d_\circ (j -e _2) + {\rm e}^{-{\rm i}2\pi \varphi (j_1 + \frac {3}{6} )} d_\bullet ^\dagger (j +e _1) d_\circ (j ) + {\rm h.c.} \bigg], \end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray} &&\fl {\cal H}_u = {\cal H}_d (d\to u,j_1\to j_1+2/6), \end{eqnarray} \tag{ 57 }$$

$$\begin{eqnarray} &&\fl {\cal H}_{ud}= \sum_j\bigg\{ \gamma _1 u_\bullet ^\dagger (j ) d_ \circ (j ) + \gamma _3 \bigg[ {\rm e}^{{\rm i}2\pi \varphi (j_1 +\frac {1}{6} )} d_\bullet ^\dagger (j) + d_\bullet ^\dagger (j+ e_1 ) \nonumber \\ &&+ {\rm e}^{-{\rm i}2\pi \varphi (j_1 +\frac {4}{6} )}d_\bullet ^\dagger (j+ e_1 -e_2 ) \bigg] u_ \circ (j-e_2 ) + {\rm h.c.} \bigg\}, \end{eqnarray} \tag{ 58 }$$

where we have used the Peierls phases to describe the effects of the magnetic field. Then the lattice periodicity is enlarged by the magnetic field. One needs a very large system for the systems with a weak magnetic field. Here 'u(d)' denote up (down) layers, γ₀ is the nearest-neighbor hopping within each layer and γ₁ is the inter-layer hopping perpendicular to the layer. The hopping γ₃, the second largest inter-layer hopping, is along an oblique direction, and causes trigonal deformation of the energy dispersion (trigonal warping) [9] in zero magnetic field. In a finite magnetic field, the hopping acquires Peierls phases appearing in the Hamiltonian above, again due to the oblique directions of the hopping.

Zoom In Zoom Out Reset image size

As to the trilayer, we have

$$\begin{equation} {\cal H}^{\rm trilayer} = {\cal H}_d + {\cal H}_u +{\cal H}_t + {\cal H}_{ud} + {\cal H}_{ut}^{ABA, ABC}, \end{equation} \tag{ 59 }$$

$$\begin{equation} {\cal H}_t = {\cal H}_d(d\to t,j_1\to j_1+4/6), \end{equation} \tag{ 60 }$$

$$\begin{equation} H_{ut}^{ABC} ={\cal H}_{ud}(d\to u, u\to t, j_1\to j_1+2/6), \end{equation} \tag{ 61 }$$

$$\begin{equation} H_{ut}^{ABA}: H_{ud} (d\to t, u \to u). \end{equation} \tag{ 62 }$$

Here we have added 't' (top layer) on top of the d,u (down, up) layers, and ABA (ABC) stands for the trilayer graphene with the ABA (ABC) stacking.

Now let us Fourier-transform d_α,u_α,t_α, first along the e₂ direction as

$$\begin{eqnarray} &&d_\alpha(j ) = \int _{-\pi }^{\pi} \frac {\mathrm{d}k_2}{2\pi} {\rm e}^{{\rm i}k_2 j_2} d_\alpha(j_1,k_2),\;{\rm etc}, \end{eqnarray} \tag{ 63 }$$

where the momentum is defined as

$$\begin{eqnarray*} &&{k} = (k_x,k_y)=k_1 \tilde {e}_1+ k_2 \tilde {e}_2,\quad \tilde {e} _i\cdot {e}_j= \delta_{ij} . \end{eqnarray*}$$

Then we have

$$\begin{eqnarray} &&\fl {\cal H}^{b} = \int_{-\pi}^\pi \frac {\mathrm{d}k_2}{2\pi} {\cal H}^{1{\rm D},b}(k_2), \quad {\cal H}^{t} = \int_{-\pi}^\pi \frac {\mathrm{d}k_2}{2\pi} {\cal H}^{1{\rm D},t}(k_2), \end{eqnarray} \tag{ 64 }$$

$$\begin{eqnarray} &&\fl {\cal H} ^{1{\rm D},b}(k_2) = {\cal H}_d^{1{\rm D}}(k_2) + {\cal H}_u^{1{\rm D}}(k_2) + {\cal H}_{du}^{1{\rm D}}(k_2), \end{eqnarray} \tag{ 65 }$$

$$\begin{eqnarray} &&\fl {\cal H} ^{1{\rm D},t}(k_2) = {\cal H}_d^{1{\rm D}}(k_2) + {\cal H}_u^{1{\rm D}}(k_2) + {\cal H}_t^{1{\rm D}}(k_2) + {\cal H}_{du}^{1{\rm D}}(k_2) + {\cal H}_{ut}^{1{\rm D},ABC(ABA)}(k_2), \end{eqnarray} \tag{ 66 }$$

$$\begin{eqnarray} &&\fl {\cal H}_d^{1{\rm D}}(k_2) = \gamma _0\sum_{j_1}\bigg[\! \big( 1\! + {\rm e}^{{\rm i}(2\pi \varphi j_1 -k_2)} \big) d_\bullet ^\dagger (j_1,k_2 ) d_\circ (j_1,k_2) + {\rm e}^{-{\rm i} 2\pi\varphi(j_1+ \frac {3}{6} )} d_\bullet ^\dagger (j_1\! +\!1,k_2) d_\circ (j_1 ,k_2) + {\rm h.c.}\! \bigg],\nonumber\\ \end{eqnarray} \tag{ 67 }$$

$$\begin{eqnarray} &&\fl {\cal H}_u^{1{\rm D}}(k_2) = {\cal H}_d^{1{\rm D}}(k_2) : d\to u, j_1\to j_1+2/6, \end{eqnarray} \tag{ 68 }$$

$$\begin{eqnarray} &&\fl {\cal H}_t^{1{\rm D}}(k_2) = {\cal H}_d^{1{\rm D}}(k_2) : d\to t, j_1\to j_1+4/6, \end{eqnarray} \tag{ 69 }$$

$$\begin{eqnarray} &&\fl {\cal H}_{du}^{\rm 1D}(k_2) = \sum_{j_1}\bigg\{ \gamma _1 u _ \bullet ^\dagger (j_1,k_2 ) d_ \circ (j_1,k_2 ) + \gamma _3 \bigg[ {\rm e}^{{\rm i}(-k_2+2\pi \varphi (j_1 +\frac {1}{6} ))} d_\bullet ^\dagger (j_1,k_2) \nonumber \\ &&+ (\mathrm{ e}^{-k_2} + {\rm e}^{-{\rm i} 2\pi \varphi (j_1+ \frac {4}{6}) } ) d_\bullet ^\dagger (j_1+1,k_2) \bigg] u_ \circ (j_1,k_2) + {\rm h.c.} \bigg\}, \end{eqnarray} \tag{ 70 }$$

$$\begin{eqnarray} &&\fl {\cal H}_{ut}^{1{\rm D},ABC}(k_2) ={\cal H}_{du}^{1{\rm D}}(k_2) : d\to u,u\to t,j_1\to j_1+2/6, \end{eqnarray} \tag{ 71 }$$

$$\begin{eqnarray} &&\fl {\cal H}_{ut}^{1{\rm D},ABA}(k_2) ={\cal H}_{du}^{1{\rm D}}(k_2) : d\to t,u\to u. \end{eqnarray} \tag{ 72 }$$

In figure 2, we have shown the spectrum of bilayer graphene on a cylinder. Here the translation symmetry along e₂ makes the wave number k₂ a good quantum number, so that the dispersion against k₂ is depicted. Due to the chiral symmetry, states near the zero energy, which are our focus, are obtained by diagonalizing the sparse matrices D^†D or DD^†, which are semi-positive definite (see equation (4)) only near the minimum energies. Numerically, it is substantially easy compared with diagonalizing H directly.

Zoom In Zoom Out Reset image size

The 'bulk-edge correspondence', which dictates that bulk topological features such as the quantized Hall conductivity have a one-to-one correspondence to edge properties, allows us to obtain the Hall conductivity of the system by counting the number of edge modes. The numerical results with edges for a weak magnetic field are given in figure 2. Counting the number of branches, the Hall conductances in units of e²/h are given in the figure.

For zero magnetic field, we can Fourier transform u_α,d_α,t_α along the e₁ as

$$\begin{eqnarray} &&u_\alpha(k_1 ,k_2 ) = \int _{-\pi }^{\pi} \frac {\mathrm{d}k_1}{2\pi} {\rm e}^{{\rm i}k_1 j_1} u_\alpha, d_\alpha, t_\alpha (j_1,k_2),\quad {\rm etc}. \end{eqnarray} \tag{ 73 }$$

The fully Fourier-transformed bilayer Hamiltonian is

$$\begin{eqnarray} &&\fl {\cal H}^{\rm bilayer} = \int \frac {\mathrm{d}^2k}{(2\pi)^2} {\cal H}({k} ), \end{eqnarray} \tag{ 74 }$$

$$\begin{eqnarray} &&\fl {\cal H}({k} ) = \gamma _0 ( 1 + {\rm e}^{-{\rm i}k_2} + {\rm e}^{-{\rm i} k_1} ) (d_\bullet ^\dagger d_\circ +u_\bullet ^\dagger u_\circ ) +\gamma _1 u _ \bullet ^\dagger d_ \circ + \gamma _3 ( {\rm e}^{-{\rm i} k_1} + {\rm e}^{-{\rm i}k_2} + {\rm e}^{-{\rm i}(k_1+k_2)} ) d_\bullet ^\dagger u_ \circ + {\rm h.c.} \end{eqnarray} \tag{ 75 }$$

$$\begin{eqnarray} &&\fl \hphantom{{\cal H}({k} )} = ( u_\bullet ^\dagger , d_\bullet ^\dagger , u_\circ ^\dagger , d_\circ ^\dagger ) H^L_{b} \left[ \begin{array}{@{}c@{}} {u_\bullet }\\ {d_\bullet }\\ {u_\circ }\\ {d_\circ}\end{array}\right], \end{eqnarray} \tag{ 76 }$$

$$\begin{eqnarray}\fl H_b^L = \left[ \begin{array}{@{}cc@{}} {O} & {D }\\ {D ^\dagger } & {O}\end{array}\right],\quad D = \left[ \begin{array}{@{}cc@{}} {\Delta } & {\gamma _1}\\ {\gamma _3'} & {\Delta }\end{array}\right], \end{eqnarray} \tag{ 77 }$$

where

$$\begin{eqnarray} &&\Delta(k) = \gamma _0(1+{\rm e}^{-{\rm i} k_1}+{\rm e}^{-{\rm i} k_2}), \end{eqnarray} \tag{ 78 }$$

$$\begin{eqnarray} &&\gamma _3' = \gamma _3({\rm e}^{-{\rm i} k_1}+{\rm e}^{-{\rm i} k_2}+{\rm e}^{-{\rm i} k_1}\,{\rm e}^{-{\rm i} k_2}). \end{eqnarray} \tag{ 79 }$$

For the trilayer, we have

$$\begin{eqnarray} &&{\cal H}^{\rm trilayer} = \int \frac {\mathrm{d}^2k}{(2\pi)^2} {\cal H}({k} ), \end{eqnarray} \tag{ 80 }$$

$$\begin{eqnarray} &&{\cal H}({k} ) = \gamma _0 ( 1 + {\rm e}^{-{\rm i}k_2} + {\rm e}^{-{\rm i} k_1} ) (d_\bullet ^\dagger d_\circ +u_\bullet ^\dagger u_\circ +t_\bullet ^\dagger t_\circ ) + {\rm h.c.} \end{eqnarray} \tag{ 81 }$$

$$\begin{eqnarray} &&({\rm for} \ ABC)\quad +\gamma _1 (t _ \bullet ^\dagger u_ \circ +u _ \bullet ^\dagger d_ \circ ) + \gamma _3' (u_\bullet ^\dagger t_ \circ +d_\bullet ^\dagger u_ \circ ) + {\rm h.c.} \end{eqnarray} \tag{ 82 }$$

$$\begin{eqnarray} &&({\rm for} \ ABA)\quad +\gamma _1 (u _ \bullet ^\dagger t_ \circ +u _ \bullet ^\dagger d_ \circ ) + \gamma _3' ( t_\bullet ^\dagger u_ \circ + d_\bullet ^\dagger u_ \circ ) + {\rm h.c.} \end{eqnarray} \tag{ 83 }$$

$$\begin{eqnarray} &&= ( t_\bullet ^\dagger , u_\bullet ^\dagger , d_\bullet ^\dagger , t_\circ ^\dagger , u_\circ ^\dagger , d_\circ ^\dagger ) H^L_{b} \left[ \begin{array}{@{}c@{}} {t_\bullet }\\ {u_\bullet }\\ {d_\bullet }\\ {t_\circ }\\ {u_\circ }\\ {d_\circ}\end{array}\right], \end{eqnarray} \tag{ 84 }$$

$$\begin{eqnarray} H_t^L = \left[ \begin{array}{@{}cc@{}} {O} & {D^{ABC,ABA} }\\ {{D^{ABC,ABA}} ^\dagger } & {O}\end{array}\right], \end{eqnarray} \tag{ 85 }$$

$$\begin{eqnarray} D^{ABC} = \left[ \begin{array}{@{}ccc@{}} {\Delta} & {\gamma _1} & { 0}\\ {\gamma' _3} & {\Delta} & {\gamma _1}\\ {0}&{\gamma' _3}&{\Delta}\end{array}\right],\quad D^{ABA} = \left[ \begin{array}{@{}ccc@{}} {\Delta} & {\gamma' _3} & { 0}\\ {\gamma _1} & {\Delta} & {\gamma _1}\\ {0} & {\gamma' _3} & {\Delta}\end{array}\right]. \end{eqnarray} \tag{ 86 }$$

Further, we can extend the above argument to consider the stacked system with a general number (p) of ABC layers. In the present representation, the Hamiltonian simplifies to

$$\begin{eqnarray} H_p=\left[ \begin{array}{@{}cc@{}} {O} & {D_p}\\ {D_p ^\dagger } & {O}\end{array}\right],\quad D_p = \Delta {\bf 1} + \gamma_1 J_p+ \gamma _3' \skew3\tilde {J}_p , \end{eqnarray} \tag{ 87 }$$

where 1 is a p-dimensional unit matrix and

$$\begin{eqnarray*} J_p= \left [\begin{array}{@{}ccc@{}} 0 & 1 & 0\\ 0 & 0 & \ddots \\ \vdots & \ddots & \ddots \end{array} \right]. \end{eqnarray*}$$

4.2.1. Low-energy Hamiltonians around K and K'

The Dirac points (K and K' points) for the monolayer graphene are specified by e^−ik0₁ = ω and e^−ik0₂ = ω², ω³ = 1, ω ≠ 1, at which Δ vanishes along with γ₃' → 0. Around the Dirac points we can expand them as

$$\begin{eqnarray} &&\Delta = - \mathrm{i}\gamma _0 \omega \xi,\quad \gamma _3' = \mathrm{i} \gamma _3 \omega ^2 \xi^*, \quad \xi= \delta k_1 + \omega \delta k_2, \end{eqnarray}\noindent \tag{ 88 }$$

where δk_i = k_i − k⁰_i. As to multilayer graphene, gapless momenta are determined by $\det \,D$'s. For the bilayer, it is given as

$$\begin{eqnarray} &&\det D^{AB}= \Delta ^2-\gamma _1 \gamma _3'= - \gamma_0^2 \omega ^2\left( \xi^2 - \frac {\mathrm{i} \gamma _1 \gamma _3}{\gamma _0^2} \xi^*\right). \end{eqnarray} \tag{ 89 }$$

This vanishes at four points, ξ = 0,ξ₀,ωξ₀,ω²ξ₀, where ξ₀ = (γ₁γ₃/γ²₀)e^iπ/6. Namely, each Dirac point proliferates into four.

For the trilayer, it is

$$\begin{eqnarray} &&\det D^{ABC} =\det D^{ABA} =\Delta ^3-2 \Delta \gamma _1 \gamma _3'. \end{eqnarray}\noindent \tag{ 90 }$$

This gives Dirac cones at $\xi = \sqrt {2}\xi _0, \sqrt {2}\xi _0\omega , \sqrt {2}\xi _0\omega ^2$ and double zero at ξ = 0. One may discuss the case of p-layers as well.

The Dirac points in the chiral-symmetric system at E = 0 are specified by $\det D=0$, which determines the gap closing momenta. Since $\det D D ^\dagger =\det D ^\dagger D=|\det D|^2$ is a product of the energies of all branches and the high-energy ones are considered as constants near the gap closing momentum, the low-energy dispersion near E = 0 is given by $\epsilon (k) \propto \pm |\det D|$. In figure 3, we have shown $|\det D|$ and Arg $\det D$ for a finite value of γ₃ in the bilayer system. From this plot, we can see that the chiralities of the four Dirac points are +1,−1,−1,−1 or −1,+1,+1,+1. This characterization of the Dirac cones by the phase of $\det D$ is only possibly using the special form of the chiral symmetric Hamiltonian.

Zoom In Zoom Out Reset image size

Now let us discuss the low-energy Hamiltonian in the absence of the trigonal warping (γ₃ = 0). In the ABC-stacked p-layered graphene (including bilayer), D_p is triangular

$$\begin{eqnarray} D_p (k)= \gamma _1 \left[ \begin{array}{@{}cccc@{}} z(k) & 1 & 0 & \\ 0 & z(k) & 1 & \ddots \\ & \ddots & \ddots & \ddots \end{array} \right], \quad z(k) = \Delta(k) /\gamma _1, \end{eqnarray} \tag{ 91 }$$

$$\begin{eqnarray} &&\det D_p (k)= \Delta ^p, \end{eqnarray} \tag{ 92 }$$

where taking a continuous limit around K and K' points implies that we are assuming |z| ≪ 1. Precisely at the K and K' points, we have D^†_pD_p = γ²₁diag (0,1,1,...), i.e. there exists only one low-energy mode for D^†_pD_p while the others have γ²₁. Then the chiral symmetry (i.e. H_p is composed of D_p as the off-diagonal block) implies that H_p has two low-energy modes, ±(k), with $\epsilon (k) = |\det D_p(k)|/\gamma _1^{p-1} = \gamma _1|z(k)|^p$. Here we can take the chiral basis introduced in previous sections to expand the low-energy doublet as

$$\begin{eqnarray} \psi=(\psi_+,\psi_- )= \left[ \begin{array}{@{}cc@{}} {\psi_\bullet } & {0}\\ {0} & {\psi_\circ }\end{array}\right],\quad D_p D_p ^\dagger \psi_\bullet =\epsilon^2 \psi_\bullet , \quad D_p ^\dagger D_p \psi_\circ =\epsilon^2 \psi_\circ . \end{eqnarray} \tag{ 93 }$$

ψ_•,○ are asymptotically (|z| ≪ 1) normalized and given as

$$\begin{eqnarray} &&\psi_\bullet = \left[ \begin{array}{@{}c@{}} {(-z^*)^{p-1} }\\ {\vdots}\\ {-z^*}\\ {1}\end{array}\right],\quad \psi_\circ = \left[ \begin{array}{@{}c@{}} {1}\\ {-z}\\ {\vdots}\\ {(-z)^{p-1} }\end{array}\right] \end{eqnarray} \tag{ 94 }$$

with $\psi _\bullet ^\dagger \psi _\bullet =\psi _\circ ^\dagger \psi _\circ =(1-|z|^p)(1-|z|^2) ^{-1} =1+{\cal O}(|z|^2)$. They are consistent with equation (93) up to the errors arising from the normalization of ψ^†_•ψ_• and ψ^†_○ψ_○. Then projecting out the high-energy sectors, we have an effective Hamiltonian for the low-energy doublet formed by the chiral basis ψ (including the monolayer case) as

$$\begin{eqnarray} H_p^{\rm eff}=\psi ^\dagger H_p \psi = -\gamma _1\left[ \begin{array}{@{}cc@{}} {0} & {(-z)^p}\\ {(-z^*)^p} & {0}\end{array}\right] = (-\gamma _1)^{-(p-1)} \left[ \begin{array}{@{}cc@{}} {0} & {\Delta ^p}\\ {{\Delta^*} ^p} & {0}\end{array}\right]. \end{eqnarray} \tag{ 95 }$$

It is given in [19] with a unitary transformation. Here the derivation is substantially simplified with the projection operator exploiting the chiral symmetry.

As to the ABA trilayer, (D^ABA)^†D^ABA has two low-energy modes around the K and K' points, which implies that the low-energy modes of H^ABA₃ are spanned by a four-dimensional orthonormalized chiral basis as

$$\begin{eqnarray} \fl \psi=(\psi_+,\psi_- )= \left[ \begin{array}{@{}cc@{}} \psi_\bullet & 0 \\ 0 & \psi_\circ \end{array} \right] = \left[ \begin{array}{@{}cccc@{}} {\psi_\bullet^{(1)} } & {\psi_\bullet^{(2)} } & {0} & {0}\\ {0} & {0} & {\psi_\circ^ {(1)} } & {\psi_\circ ^{(2)} } \end{array} \right], \end{eqnarray} \tag{ 96 }$$

$$\begin{eqnarray} &&\fl \psi ^\dagger \psi= {\rm diag} ( 1,1,1,1)+{\cal O}(|z|), \end{eqnarray} \tag{ 97 }$$

$$\begin{eqnarray} &&\fl D_p D_p ^\dagger \psi_\bullet^{(i)} =\epsilon_ i^2 \psi_\bullet^{(i)} , \quad D_p ^\dagger D_p \psi_\circ^{(i)} =\epsilon_i^2 \psi_\circ ^{(i)} ,\quad \epsilon_1 = |\Delta |,\quad \epsilon_2 = \gamma _1 ^{-1} |\Delta |^2/2, \end{eqnarray} \tag{ 98 }$$

where we have, up to the leading order,

$$\begin{eqnarray} \psi_\bullet^{(1)} &=& \frac {1}{\sqrt{2}} \left[ \begin{array}{@{}c@{}} {-1}\\ {0}\\ {1}\end{array}\right], \quad \psi_\bullet^{(2)} = \frac {1}{2} \left[ \begin{array}{@{}c@{}} {-z^*}\\ {2}\\ {-z^*}\end{array}\right], \nonumber \\\\ \psi_\circ ^{(1)} &=& \frac {1}{\sqrt{2}} \left[ \begin{array}{@{}c@{}} {-1}\\ {0}\\ {1}\end{array}\right], \quad \psi_\circ ^{(2)} = \frac {1}{\sqrt{2}} \left[ \begin{array}{@{}c@{}} {1}\\ {-z}\\ {1}\end{array}\right].\nonumber \end{eqnarray} \tag{ 99 }$$

Now we have

$$\begin{eqnarray*} \psi ^\dagger H_3^{ABA} \psi = \left[ \begin{array}{@{}cc@{}} {O} & {\psi_\bullet ^\dagger D \psi_ \circ }\\ {\psi_\circ ^\dagger D ^\dagger \psi_\bullet } & {O}\end{array}\right] \end{eqnarray*}$$

and

$$\begin{eqnarray*} \psi_\bullet ^\dagger D \psi_ \circ =\gamma _1\left[ \begin{array}{@{}cc@{}} {z} & {0}\\ {0} & {-z^2/\sqrt{2}}\end{array}\right], \end{eqnarray*}$$

so that we end up with a simple decomposition,

$$\begin{eqnarray} H_3^{ABA:{\rm eff}}=\gamma _1\left[ \begin{array}{@{}cc@{}} {0} & {z ^\dagger }\\ {z} & {0}\end{array}\right]\oplus \frac {\gamma _1}{\sqrt{2}} \left[ \begin{array}{@{}cc@{}} {0} & {(-z ^\dagger )^2}\\ {-z^2} & {0}\end{array}\right] = H_1 \oplus H_2/\sqrt{2}. \end{eqnarray} \tag{ 100 }$$

4.2.2. Landau levels

Around the Dirac cones, we can expand the effective Hamiltonian, as in the monolayer case, as

$$\begin{eqnarray} &&\fl H_1\to H_1^{\rm eff} = ( \sigma \cdot X) \delta k_x + (\sigma \cdot Y )\delta k_y, \quad X = \left[ \begin{array}{@{}c@{}} { {\rm Re}\,\Delta_x }\\ {-{\rm Im}\,\Delta_x }\\ {0}\end{array}\right],\quad Y = \left[ \begin{array}{@{}c@{}} {{\rm Re}\,\Delta _y}\\ {-{\rm Im}\,\Delta _y}\\ {0}\end{array}\right], \end{eqnarray} \tag{ 101 }$$

where $\Delta _\alpha =\frac {\partial \Delta }{\partial k_\alpha }\big |_{k=K,K'}$. Then we have

$$\begin{eqnarray} &&\fl (H_1^{\rm eff})^2 =(\hbar c \overline{\delta k})^2,\quad \overline{\delta k}^2 \equiv [\delta k_x,\delta k_y] \Xi \left[ \begin{array}{@{}c@{}} {\delta k_x}\\ {\delta k_y}\end{array}\right],\quad \Xi= \frac {1}{|X\times Y|} \left[ \begin{array}{@{}cc@{}} {X\cdot X} & {Y\cdot X}\\ {X\cdot Y} & {Y\cdot Y}\end{array}\right], \end{eqnarray} \tag{ 102 }$$

where c² = |X × Y |/ℏ² is an effective 'light velocity' and $\overline {\delta k}$ is the averaged momentum near the Dirac cones with $\det \Xi =(|X|^2|Y|^2-|X\cdot Y|^2)/|X\times Y|^2=1$. Note that

$$\begin{eqnarray} &&X\times Y = \left[ \begin{array}{@{}c@{}} {0}\\ {0}\\ {{\rm Im}\, \Delta _x \Delta _y^*}\end{array}\right] =\left[ \begin{array}{@{}c@{}} {0}\\ {0}\\ {\chi \hbar ^2 c^2}\end{array}\right], \end{eqnarray} \tag{ 103 }$$

where χ = sgn Im Δ_xΔ*_y is the chirality of the Dirac cones.

In a magnetic field, the Hamiltonian H^C₁ in the continuum limit is obtained by replacing ℏδk_α → π_α = p_α − eA_α = −iℏ∂_α − eA_α where [π_x,π_y] = iℏeB = i(ℏ/ℓ_B)² with $\ell _B=\sqrt {\hbar /eB}$ (eB > 0). As to the ABC stacked p-layered graphene (including monolayer and bilayer), the effective continuum Hamiltonian is obtained as the extension of McCann–Fal'ko [17] as

$$\begin{eqnarray} H_p^{\rm eff} \to H_p^{\mathrm{C}} \equiv (- \gamma _1)^{-(p-1)} \left[ \begin{array}{@{}cc@{}} {0} & {(\Pi ^\dagger)^p }\\ {\Pi^p} & {0}\end{array}\right], \end{eqnarray} \tag{ 104 }$$

$$\begin{eqnarray} &&\Delta \to \Pi^\dagger \equiv \hbar ^{-1} (\Delta _x \pi_x+\Delta _y \pi_y). \end{eqnarray} \tag{ 105 }$$

This is another derivation of the pth layer effective Hamiltonian discussed by Koshino–McCann [19] for the chiral symmetric case. Here the commutation relation of Π is given as

$$\begin{eqnarray} &&[\Pi,\Pi ^\dagger ] =2\hbar ^{-2} {\rm Im\,} \Delta _x \Delta _y^* = 2(c\hbar/\ell_B)^2 \chi. \end{eqnarray} \tag{ 106 }$$

Then one can define a canonical boson operator a ([a,a^†] = 1) as

$$\begin{eqnarray} \Pi /(\sqrt{2}c\hbar/\ell_B) =\left\{ \begin{array}{@{}cc@{}} {a} & {(\chi=+1)},\\ {a ^\dagger} & {(\chi=-1)}\end{array}\right. . \end{eqnarray} \tag{ 107 }$$

For the positive chirality χ = + 1, one has

$$\begin{eqnarray} \fl (H_p^{\mathrm{C}} )^2= N_p^2 \left[ \begin{array}{@{}cc@{}} {(a ^\dagger )^p a ^p} & {0}\\ {0} & {a ^p(a ^\dagger )^p }\end{array}\right] = N_p ^2 \left[ \begin{array}{@{}cc@{}} {\prod_{i=1}^p(n-i+1)} & {0}\\ {0} & {\prod_{i=1}^p(n+i)}\end{array}\right], \end{eqnarray} \tag{ 108 }$$

where $N_p=\gamma _1(\gamma _1 ^{-1} c \sqrt {2e\hbar B})^{p}$, n = a^†a, and

$$\begin{eqnarray} &&(a ^\dagger )^2 a ^2 =a ^\dagger n a =(n-1)n,\ (a ^\dagger )^3 a ^3 =(n-2)(n-1) n,\ \ldots, \end{eqnarray} \tag{ 109 }$$

$$\begin{eqnarray} &&a ^2 (a ^\dagger )^2 = (n+2)(n+1), \ldots. \end{eqnarray} \tag{ 110 }$$

Then we have a series of non-equally spaced Landau levels

$$\begin{equation} E_n = \pm N_p\sqrt{n(n-1)(n-2)\cdots(n-p+1)}. \end{equation} \tag{ 111 }$$

The χ = −1 case follows trivially.

In this section, let us focus on the Lifshitz transition caused by the trigonal warping term γ₃, by directly calculating the energy levels on a honeycomb lattice without using the low-energy effective model [17, 18]. We stress here that our numerical calculations are performed for the tight-binding model on a honeycomb lattice with the Peierls phases for the magnetic field which require diagonalization of large matrices for systems with a weak magnetic field. This becomes only possible by the chiral symmetry using the sparse matrix technique. It is enough to diagonalize the sparse matrix D^†D (see equation (4)) near the zero energies. The results are compared with the discussion of the Lifshitz transition using the effective low-energy model in [17, 18]. However, the crossover phenomenon from the lattice model to the low-energy effective Hamiltonian caused by γ₁ is only possible by the present direct lattice calculation (see figure 4). The Landau levels of the bilayer graphene as a function of γ₁ are shown in figure 4. One can observe the crossover, that is, the γ₁ produces bonding–antibonding splitting, so that the low-energy sector is projected out from the high-energy one.

Zoom In Zoom Out Reset image size

Figure 5, on the other hand, displays the Landau levels as a function of the γ₃. One can see that the Landau levels at γ₃ = 0 are adiabatically connected to two classes of Landau levels at γ₃ ≠ 0, which correspond to those of two kinds of four Dirac cones with different Fermi velocities. To observe the phenomenon at the moderate value of the flux, we take an artificially large value for γ₁/γ₀ = 0.5 in the left panel of figure 5. In the right panel, the same plot for a much smaller, realistic γ₁ and the magnetic field is shown. Numerical calculations for such a weak magnetic field become possible only by imposing the chiral symmetry. Then the Lifshitz transition with realistic parameters is clearly demonstrated using the tight-binding model.

Zoom In Zoom Out Reset image size

We start with the optical conductivities for the bilayer system with γ₀,γ₁ terms (equation (104) with p = 2). In figures 6(a) and (b), we show optical longitudinal and Hall conductivities (σ_xx(_F,ω),σ_xy(_F,ω)) plotted against the Fermi energy _F and the frequency ω. Energy and frequency are normalized with the cyclotron frequency ω_c = 2v²eB/γ₁ for bilayer parabolic bands. If we label LLs with the Landau index n and an electron/hole band index s = ± , the selection rule for bilayer graphene is described with intra-band transitions (n,s) → (n + 1,s) and inter-band transitions (n,−s) → (n + 1,s). For both σ_xx(_F,ω),σ_xy(_F,ω), intra-band transitions (n,s) → (n + 1,s) occur around

$$\begin{eqnarray*} &&\omega_{\rm intra} \sim \omega_{\rm c}=2 v^2 e B/\gamma_1, \end{eqnarray*}$$

with the Fermi velocity of the Dirac cone v (=c in section 4), since LLs are almost equally spaced, unlike in the monolayer case where LL energy $\propto \sqrt {n}$ is not equally separated. The inter-band transitions across the band-touching point occur around

$$\begin{eqnarray*} &&\hbar\omega_{\rm inter} \simeq 2 |\epsilon_{\rm F}|, \end{eqnarray*}$$

for large enough n. When both transitions (n,−) → (n + 1,+) and (n + 1,−) → (n,+) are allowed, these contributions for optical conductivities add up for σ_xx(_F,ω) while they cancel out for σ_xy(_F,ω) due to the sign of current matrix elements, so that we have inter-band transitions wide in _F for σ_xx(_F,ω) while they are seen in σ_xy(_F,ω) only for a region of Fermi energy where one of the transitions is allowed.

Zoom In Zoom Out Reset image size

For ABA stacked trilayer, the effective Hamiltonian around K₊/K₋ points is given by a 6 × 6 matrix (the dimension being 2 sublattices × 3 layers) as [25–28]

$$\begin{equation} H_{ABA}= \left( \begin{array}{@{}cccccc@{}} 0 & v \pi^\dagger & 0 & v_3 \pi & \gamma_2/2 & 0\\ v \pi & \Delta' & \gamma_1 & 0 & 0 & \gamma_5/2 \\ 0 & \gamma_1 & \Delta' & v \pi^\dagger & 0 & \gamma_1\\ v_3 \pi^\dagger & 0 & v \pi & 0 & v_3 \pi^\dagger & 0 \\ \gamma_2/2 & 0 & 0 & v_3 \pi & 0 & v \pi^\dagger \\ 0 & \gamma_5/2 & \gamma_1 & 0 & v \pi & \Delta' \end{array} \right), \end{equation} \tag{ 113 }$$

with π = π_x + iπ_y, a velocity v₃ associated with γ₃ as v₃ = vγ₃/γ₀, Δ' the on-site energy difference between the atoms with and without vertical bond γ₁, and γ₂ (γ₅) the next-nearest interlayer hoppings between A1 and A3 (B1 and B3).

First, we only retain γ₀,γ₁ terms in equation (113), in which case the Hamiltonian is chiral-symmetric. In figures 7(a) and (b), the optical longitudinal and Hall conductivities σ_xx(_F,ω),σ_xy(_F,ω) are plotted against the Fermi energy _F and frequency ω. The magnetic field is chosen so that an interlayer hopping energy is γ₁/ℏω_c = 5, with the monolayer cyclotron energy $\omega _{\mathrm { c}}= \sqrt {2} v/\ell =v \sqrt {2eB/\hbar }$. We can see a monolayer contribution (Dirac cyclotron frequency =ω_c) and a bilayer contribution with a cyclotron energy, $\hbar \omega _{\rm {bilayer}}/ \sqrt {2},$ both of which show intra-band and inter-band transitions. For moderate magnetic fields B ∼ 1 T, the cyclotron frequency for the monolayer is much larger than that for the bilayer, ω_c ≫ ω_bilayer. In this case, σ_xx is an even function with respect to _F = 0, while σ_xy is odd due to the electron–hole symmetry, a consequence of the chiral symmetry. The jump of σ_xy at _F = 0 is related to the chiral zero modes.

Zoom In Zoom Out Reset image size

Now, if we include all the hopping terms in equation (113), the chiral symmetry is broken, and we no longer have the chiral protection for zero modes (three zero-energy LLs). We have then a massive Dirac band plus gapped bilayer bands with energy shifts. As to hopping parameters, we adopt the values for graphite, γ₀ = 3.2 eV, γ₁ = 0.39 eV, γ₃ = 0.32 eV, γ₂ = −0.020 eV, γ₅ = 0.038 eV, Δ' = 0.050 eV [9, 29].

Figures 8(a) and (b) depict the result for the optical longitudinal σ_xx(_F,ω) and Hall conductivity σ_xy(_F,ω) plotted against the Fermi energy _F and frequency ω for ABA-stacked trilayer graphene, including all the hoppings. We see contributions from monolayer-like Dirac LLs (labeled with M) and those from bilayer LLs (B), both of which exhibit intra-band and inter-band transitions. Since the Dirac cone is massive due to the chiral symmetry breaking and the zero-energy LL for the monolayer band (M0) is situated at the bottom of the conduction band for K₊ valley and the top of the valence band for K₋, M0 → M(1,+) resonance occurs at an energy lower than M(1,−) → M0 for K₊, and vice versa for K₋ (figure 8(c)). A cancellation of resonances in σ_xy, due to the opposite signs in current matrices, occurs between M(1,−) → M0 for K₊ and M0 → M(1,+) for K₋ for a region of Fermi energy between M0(K₊) and M0(K₋), while this is not the case with σ_xx. For bilayer contributions, satellites appear as $B(n,\pm ) \leftrightarrow B(n+1+3m,\pm )$ and $B(n,\pm ) \leftrightarrow B(n+2+3m,\pm )$, since the trigonal warping term mixes n LLs and n + 3 LLs [30, 31]. The resonance frequency for intra-band transition within the conduction band is larger than that within the valence band, which is a consequence of the electron–hole asymmetry in the bilayer bands, triggered by a breaking of the chiral symmetry. A deviation in the cyclotron mass for electron and hole bands prevents complete cancellation between B(n,−) → B(n + 1,+) and B(n + 1,−) → B(n,+) transitions, which results in small inter-band transitions in a wide region of Fermi energy.

Zoom In Zoom Out Reset image size

The low-energy effective Hamiltonian of ABC-stacked trilayer graphene in a 2 × 2 matrix equation (104) is a cubic form in the momentum, if we neglect hopping terms but γ₀,γ₁. For this effective Hamiltonian, LL energy, equation (111) shows a magnetic-field dependence $\propto B^{\frac {3}{2}}$, resulting in a smaller LL spacing compared to the single-layer LL $\propto B^{\frac {1}{2}}$ and bilayer LL ∝B for weak magnetic fields.

Now we turn to a result for the optical conductivities in ABC trilayer for an interlayer hopping energy γ₁/ℏω_c = 5, with monolayer cyclotron frequency $\omega _{\mathrm { c}}= v \sqrt {2eB/\hbar }$. Figures 9(a) and (b) show the optical longitudinal and Hall conductivity σ_xy(_F,ω) plotted against the Fermi energy _F and frequency ω for ABC-stacked trilayer graphene in strong magnetic fields. In contrast to ABA stacking (figure 7), we only see a single sequence of intra-band and inter-band transitions with a much smaller cyclotron energy than ℏω_c due to the dependence on magnetic fields $\propto B^{\frac 3 2}$ due to the LL energy structure (equation (111) with p = 3). The inter-band transition occurs at ℏω ∼ 2_F for the same reason as in the bilayer, while the intra-band transition energy grows with increasing LL index n as

$$\begin{eqnarray*} &&\omega \sim (n+1)^{3/2}-n^{3/2} \sim n^{\frac 1 2}, \end{eqnarray*}$$

which explains why the intra-band resonance energy increases with n in figure 9(a), while it decreases with n in the case of monolayer graphene [22]. We can then predict in general that, for ABC p-layered graphene, the intra-band transition occurs at

$$\begin{eqnarray*} &&\omega\sim n^{p/2-1}, \end{eqnarray*}$$

and the inter-band transition at ω ∼ 2_F. Thus, we end up with intra-band transition energies that exhibit different behaviors with the Landau index n as

$$\begin{eqnarray*} &&\propto n^{-1/2} \;{\rm (monolayer)},\quad {\rm const. \; (bilayer)},\quad \propto n^{1/2} \;(ABC\;{\rm trilayer}), \end{eqnarray*}$$

while the inter-band transition energies are qualitatively the same with ∼2_F.

Zoom In Zoom Out Reset image size

Finally, let us mention the effects of electron–electron interactions on the optical responses. It is argued that the many-body gap is open for bilayer and trilayer graphene [32, 33]. This many-body gap causes a splitting of n = 0 LL and the optical transitions associated with n = 0 LL would also be split as we have seen in ABA trilayer with all hoppings. Another interesting problem is the effects of the inter-valley plasmon and possible CDW [34], which can cause a similar splitting of optical transitions or low-energy modifications of optical spectra.