Many-body effects, orbital mixing and cyclotron resonance in bilayer graphene

Graphene supports as charge carriers massless Dirac fermions which display unique and fascinating electronic properties. Recently, increasing attention is directed to bilayers NMMKF ; MF and fewlayers of graphene, where the physics and applications of graphene become far richer, with, e.g., a tunable band gap MF ; OBSHR in bilayer graphene.

In a magnetic field, graphene leads to a $``$relativistic” infinite tower of electron and hole Landau levels, with the lowest Landau level (LLL) consisting of a number of zero-energy levels. The presence of such zero-mode levels has a topological origin in the index of (the leading part of) the Dirac operators, or in the chiral anomaly NS . Monolayer graphene carries $2_{\rm valley}\times 2_{\rm spin}=4$ zero-mode levels. Bilayer graphene supports an octet of such levels, with an extra twofold degeneracy in Landau orbitals $n=0$ and $n=1$. In real samples, topological zero-energy levels evolve, by an interplay of general spin and valley splittings and Coulomb interactions, into a variety of pseudo-zero-mode (PZM) levels, or broken-symmetry quantum Hall states, as discussed theoretically BCNM ; KSpzm ; BCLM ; CLBM ; CLPBM ; NL ; MAF . Meanwhile it was noted that orbital degeneracy is also lifted by Coulomb interactions alone KS_Ls : Quantum fluctuations of the valence band split, like the Lamb shift in the hydrogen atom, the $n=0$ and $n=1$ levels appreciably. The orbital degeneracy and its lifting by the orbital Lamb shift are new features specific to the LLL in fewlayer graphene. Experimentally, partial or full resolution of the eightfold degeneracy of the LLL in bilayer graphene has been observed in transport FMY ; ZCZJ ; WAFM ; MEM ; VJB and capacitance MFW ; KFLH ; HLZW measurements, but the detailed structure of the LLL remains yet to be clarified.

Graphene supports cyclotron resonance (CR) in a variety of channels, both intraband and interband AbF . Interband resonance is specific to Dirac electrons, with a number of competing resonance channels simultaneously activated over a certain range of the total Landau-level filling factor $\nu$. CR provides a useful means to look into the many-body problem in graphene. In experiment CR has been explored for monolayer JHTWS ; DCNN ; HCJL and bilayer HJTS ; OFBB graphene. Many-body effects on CR IWF ; BM ; KS_CR ; RFG ; KS_CR_eh ; ST ; KS_CRnewGR ; SL were first detected rather indirectly by a comparison of some leading intraband and interband resonances JHTWS . Recently, Russell et al. RZTW have reported a direct signal of many-body effects on CR in high-mobility hBN-encapsulated monolayer graphene; see also ref. prepH . They observed significant variations of resonance energies over a finite range of filling factor $\nu$ under a fixed magnetic field $B$; the point lies in extracting genuine Coulombic effects as a function of $\nu$ while controlling the running of the velocity factor with $B$. The data show a profile nearly symmetric in $\nu$ and suggest a small band gap HuntYYY due to weak coupling to the hBN substrate.

No such $\nu$-dependent study of interband CR is yet available for bilayer graphene. The purpose of this paper is to examine, in anticipation of a future experiment, the detailed characteristics of the LLL in bilayer graphene, with full account taken of spin splitting, interlayer bias, weak electron-hole asymmetry and Coulomb interactions, and to show how they are detected by a $\nu$-dependent and fixed-$B$ survey of interband channels of CR. Particular attention is paid to electron-hole ($e$-$h$) conjugation and valley-interchange operations which relate the electron and hole bands within a valley or between the two valleys. We clarify how they govern the Landau-level and CR spectra.

In Sec. II we review some basic features of the effective theory of bilayer graphene in a magnetic field. In Sec. III we consider Coulombic corrections to Landau levels and CR. Such corrections inevitably contain ultraviolet divergences and, in Sec. IV, we carry out renormalization to extract observable many-body contributions. In Secs. V and VI, we examine how the PZM level spectra evolve with filling, note the general presence of orbital mixing and discuss how the associated many-body effects are revealed by interband CR. Section VII is devoted to a summary and discussion.

The electrons in bilayer graphene are described by four-component spinor fields on the four inequivalent sites $(A,B)$ and $(A^{\prime},B^{\prime})$ in the Bernal-stacked bottom and top layers of honeycomb lattices of carbon atoms. Their low-energy features are governed by the two Fermi points $K$ and $K^{\prime}$ in the Brillouin zone. The intralayer coupling $\gamma_{0}\equiv\gamma_{AB}\sim 3$  eV is related to the Fermi velocity $v=(\sqrt{3}/2)\,a_{\rm L}\gamma_{0}/\hbar\sim 10^{6}$ m/s ($a_{\rm L}=0.246$nm) in monolayer graphene. Interlayer hopping via the dimer coupling $\gamma_{1}\equiv\gamma_{A^{\prime}B}\sim 0.4$ eV makes the spectra quasi-parabolic MF in the low-energy branches $|\epsilon|<\gamma_{1}$.

The bilayer Hamiltonian with $\gamma_{0}$ and $\gamma_{1}$ leads to $e$-$h$ symmetric energy spectra. Infrared spectroscopy, however, detected weak asymmetry ZLBF ; LHJ_asym due to the $(A,B)$-sublattice energy difference $\Delta\sim 18$meV and the nonleading interlayer coupling $\gamma_{4}\equiv\gamma_{AA^{\prime}}\sim 0.04\gamma_{0}$.

The effective Hamiltonian with such intra- and inter-layer couplings is written as MF ; NCGP

	$\displaystyle H$	$\displaystyle=$	$\displaystyle\!\int\!d^{2}{\bf x}\,\Big{[}(\Psi^{K})^{{\dagger}}\,{\cal H}_{K}\Psi^{K}+(\Psi^{K^{\prime}})^{{\dagger}}\,{\cal H}_{K^{\prime}}\,\Psi^{K^{\prime}}\Big{]},$
	$\displaystyle{\cal H}_{K}$	$\displaystyle=$	$\displaystyle\!\!\left(\begin{array}[]{cccc}{1\over{2}}u&0&-v_{4}\,p^{{\dagger}}&v\,p^{{\dagger}}\\ 0&-{1\over{2}}u&v\,p&-v_{4}\,p\\ -v_{4}\,p&v\,p^{{\dagger}}&\Delta-{1\over{2}}u&\gamma_{1}\\ v\,p&-v_{4}\,p^{{\dagger}}&\gamma_{1}&\Delta+{1\over{2}}u\\ \end{array}\right),$		(5)

with $p=p_{x}+i\,p_{y}$, $p^{{\dagger}}=p_{x}-i\,p_{y}$ and $v_{4}\equiv(\gamma_{4}/\gamma_{0})\,v$. Here $\Psi^{K}=(\psi_{B^{\prime}},\psi_{A},\psi_{B},\psi_{A^{\prime}})^{\rm t}$ denotes the electron field in valley $K$, with $A$ and $B$ referring to the associated sublattices; $u$ stands for an interlayer bias, which opens a tunable valley gap OBSHR ; MF . We ignore the effect of trigonal warping $\propto\gamma_{3}\equiv\gamma_{AB^{\prime}}\sim 0.1$ eV which, in a strong magnetic field, causes only a negligibly small level shift. ${\cal H}_{K}$ is diagonal in the (suppressed) electron spin.

As a typical set of band parameters we adopt those by recent theoretical calculations JM ; fn ,

	$\displaystyle v$	$\displaystyle\approx$	$\displaystyle 0.845\times 10^{6}\,{\rm m/s},\ \gamma_{1}\approx 361\,{\rm meV},$
	$\displaystyle v_{4}/v$	$\displaystyle\equiv$	$\displaystyle\gamma_{4}/\gamma_{0}\approx 0.053,\ \triangle\approx 15\,{\rm meV},$		(6)

which are also used in an experimental analysis HLZW .

Let us note that ${\cal H}_{K}$ is unitarily equivalent to $-{\cal H}_{K}$ with the signs of $(u,\Delta,v_{4})$ reversed,

$$S^{{\dagger}}{\cal H}_{K}\,S=-{\cal H}_{K}|_{-u,-\Delta,-v_{4}},$$

(7)

with $S={\rm diag}(-1,1,-1,1)$; ${\cal O}|_{-u,-\Delta,-v_{4}}$ signifies setting $(u,\Delta,v_{4})\rightarrow(-u,-\Delta,-v_{4})$ in ${\cal O}$. The electron and hole bands are therefore intimately related within a valley.

The Hamiltonian ${\cal H}_{K^{\prime}}$ in another valley is given by ${\cal H}_{K}$ with $(v,v_{4},u)\rightarrow(-v,-v_{4},-u)$, and acts on a spinor of the form $\Psi^{K^{\prime}}=(\psi_{A},\psi_{B^{\prime}},\psi_{A^{\prime}},\psi_{B})^{\rm t}$. Actually, ${\cal H}_{K^{\prime}}$ is unitarily equivalent to ${\cal H}_{K}$ with the sign of $u$ reversed,

$${\cal H}_{K^{\prime}}={\cal U}^{{\dagger}}{\cal H}_{K}|_{-u}\,{\cal U},$$

(8)

with ${\cal U}={\rm diag}(1,1,-1,-1)$. In what follows we adopt ${\cal H}_{K}|_{-u}$ for ${\cal H}_{K^{\prime}}$ and simply pass to valley $K^{\prime}$ by reversing the sign of bias $u$ in valley-$K$ expressions. We also suppose, without loss of generality, that $u\geq 0$ in valley $K$; thus $u\leq 0$ refers to valley $K^{\prime}$.

The unitary equivalence

$${\cal H}_{K^{\prime}}\stackrel{{\scriptstyle\cal U}}{{\sim}}{\cal H}_{K}|_{-u}\stackrel{{\scriptstyle S}}{{\sim}}-{\cal H}_{K}|_{-\Delta,-v_{4}}$$

(9)

implies that the $e$-$h$ conjugation symmetry is kept exact only for $\Delta=v_{4}=0$ although, if $u\not=0$, it is apparently broken within each valley. This symmetry analysis also holds in the presence of Coulomb interactions.

Let us place bilayer graphene in a strong uniform magnetic field $B_{z}=B>0$ normal to the sample plane; we set, in ${\cal H}_{K}$, $p\rightarrow\Pi=p+eA$ with $A=A_{x}+iA_{y}=-By$, and rescale $\Pi=-(\sqrt{2}/\ell)\,Z^{{\dagger}}$ so that $[Z,Z^{{\dagger}}]=1$, where $\ell=1/\sqrt{eB}$ denotes the magnetic length.

The eigenmodes of ${\cal H}_{K}$ are labeled by integers $N\equiv|n|=(0,1,2,\cdots)$ and plane waves with momentum $p_{x}$, and have the structure

$$\Psi^{n}=\Big{(}|N\!-\!2\rangle\,b^{\prime}_{n},|N\rangle\,c_{n},|N\!-\!1\rangle\,b_{n},|N\!-\!1\rangle\,c^{\prime}_{n}\Big{)}^{\rm t},$$

(10)

where only the orbital modes are shown using the standard harmonic-oscillator basis $\{|N\rangle\}$, with $|N\rangle=0$ for $N<0$. The coefficients ${\bf b}_{n}=(b^{\prime}_{n},c_{n},b_{n},c^{\prime}_{n})^{\rm t}$ for each $N$ are given by the normalized eigenvectors of the reduced Hamiltonian

$$\hat{\cal H}_{N}=\omega_{c}\left(\begin{array}[]{cccc}\mu&0&r\,\sqrt{N\!-\!1}&-\sqrt{N\!-\!1}\\ 0&-\mu&-\sqrt{N}&r\,\sqrt{N}\\ r\,\sqrt{N\!-\!1}&-\sqrt{N}&d-\mu&g\\ -\sqrt{N\!-\!1}&r\,\sqrt{N}&g&d+\mu\\ \end{array}\right),$$

(11)

where

$$\omega_{c}\equiv\sqrt{2}\,v/\ell\approx 36.3\times v[10^{6}{\rm m/s}]\,\sqrt{B[{\rm T}]}\ {\rm meV},$$

(12)

is the cyclotron energy for monolayer graphene; $g\equiv\gamma_{1}/\omega_{c}$, $\mu\equiv{1\over{2}}\,u/\omega_{c}$, $d\equiv\Delta/\omega_{c}$ and $r\equiv\gamma_{4}/\gamma_{0}$. Numerically, for the set of parameters in Eq. (6) at $B=20$ T,

$$\omega_{c}\approx 137\,{\rm meV},\ (g,d,r)\approx(2.63,0.109,0.053),$$

(13)

and $\mu\approx 0.073$ for $u=20$meV. In what follows, we employ this set (13) of parameters for numerical estimates at $B=20$ T.

The unitary equivalence in Eq. (9) reveals how the spectra $\epsilon_{n}$ and the associated eigenmodes ${\bf b}_{n}$ change via $e$-$h$ conjugation. Within each valley they read

	$\displaystyle\epsilon_{-n}=-\epsilon_{n}|_{-u,-\Delta,-v_{4}},$
	$\displaystyle(b^{\prime}_{-n},c_{-n},b_{-n},c^{\prime}_{-n})=(-b^{\prime}_{n},c_{n},-b_{n},c^{\prime}_{n})|_{-u,-\Delta,-v_{4}},\ \ \ \$		(14)

while, between the two valleys, they are related as

	$\displaystyle(\epsilon_{n};b^{\prime}_{n},c_{n},b_{n},c^{\prime}_{n})|^{K^{\prime}}=(\epsilon_{n};b^{\prime}_{n},c_{n},b_{n},c^{\prime}_{n})|^{K}_{-u}$
	$\displaystyle=(-\epsilon_{-n};b^{\prime}_{-n},c_{-n},-b_{-n},c^{\prime}_{-n})|^{K}_{-\Delta,-v_{4}}.$		(15)

Of our particular concern are the $n=0$ and $n=1$ modes. For $n=0$, $\hat{\cal H}_{N=0}$ has an obvious eigenvalue and eigenvector

$$\epsilon_{0}=-u/2=-\omega_{c}\mu,\ \ {\bf b}_{0}=(0,1,0,0)^{\rm t}.$$

(16)

For $N=1$, $\hat{\cal H}_{N}$ has three solutions $(\epsilon_{-1^{\prime}},\epsilon_{1},\epsilon_{1^{\prime}})$, with $\epsilon_{\pm 1^{\prime}}\sim\pm\gamma_{1}$ belonging to the higher branches.

The $n=1$ mode (with $\epsilon_{1}$) and the $n=0$ mode have zero energy for $(u,\Delta,r)=0$ and deviate from zero energy as $(u,\Delta,r)$ develop. These pseudo-zero modes (PZM) form the LLL in bilayer graphene. The eigenenergy and eigenvector of the $n=1$ mode are written as

	$\displaystyle\epsilon_{1}$	$\displaystyle=$	$\displaystyle\omega_{c}(-\mu+\kappa),$
	$\displaystyle{\bf b}_{1}$	$\displaystyle=$	$\displaystyle c_{1}\Big{(}0,1,-(g\kappa-r)/\tilde{g},(1+\kappa d-d^{2})/\tilde{g}\Big{)},\ \ \$
	$\displaystyle\tilde{g}$	$\displaystyle=$	$\displaystyle g+(d-\kappa)\,r.$		(17)

In what follows, we denote valley and $e$-$h$ symmetry breaking collectively as $X=(u,\Delta,v_{4})$ or $X=(\mu,d,r)$. As seen from the secular equation, the correction $\kappa\equiv\kappa(g;\mu,d,r)$ in $\epsilon_{1}$ is odd in $X$, and hence the normalization factor $c_{1}$ is even in $X$,

	$\displaystyle\kappa$	$\displaystyle\approx$	$\displaystyle{2\mu+d+2r\,g\over{g^{2}+1}}+O(X^{3}),$
	$\displaystyle c_{1}$	$\displaystyle\approx$	$\displaystyle 1/\sqrt{1+(1/g^{2})}+O(X^{2}).$		(18)

Numerically, at $B=20$ T, $\kappa\approx 0.049+0.25\,\mu$ while $c_{1}\approx 0.934-0.007\,\mu-0.02\,\mu^{2}$ scarcely depends on $\mu$ [for small $\mu\sim O(0.1)$].

Interlayer bias $u\propto\mu$, if nonzero, breaks valley symmetry and shifts the PZM levels $n=\{0,1\}$ oppositely ($\propto\mp u/2$) in the two valleys. Accordingly, $(\epsilon_{0},\epsilon_{1})$ are nearly degenerate in each valley and are ordered so that $\epsilon_{0}^{K}\lesssim\epsilon_{1}^{K}<0<\epsilon_{1}^{K^{\prime}}\lesssim\epsilon_{0}^{K^{\prime}}$ for $\mu>0$.

For $N\geq 2$, $\hat{\cal H}_{N}$ has rank 4 and we denote the four branches of Landau-level spectra as $\epsilon_{-n^{\prime}}<\epsilon_{-n}<0<\epsilon_{n}<\epsilon_{n^{\prime}}$ (with $|\epsilon_{\pm n^{\prime}}|\gtrsim\gamma_{1}$) so that the index $(n,n^{\prime})$ reflects the sign of $\epsilon_{n}$. Let us denote $\epsilon_{n}$ in units of $\omega_{c}$,

$$\epsilon_{n}=\omega_{c}\,e_{n}\ {\rm and}\ \epsilon_{n^{\prime}}=\omega_{c}\,e_{n^{\prime}}.$$

(19)

For $X=(\mu,d,r)\rightarrow 0$, the (dimensionless) spectra $e_{n}\equiv e_{n}(g;\mu,d,r)$ of the lower branches take the form

$$e_{n}^{(0)}=s_{n}\sqrt{{2|n|(|n|-1)\over{g^{2}+a_{n}+\sqrt{D_{n}}}}}\sim O(1/g),$$

(20)

where $e_{n}^{(0)}\equiv e_{n}|_{X=0}$, $D_{n}=g^{4}+2a_{n}g^{2}+1$ and $a_{n}=2\,|n|-1$; $s_{n}\equiv{\rm sgn}[n]=\pm 1$ is the sign function. The full spectra, to first order in breaking $X$, are written as

$$e_{n}\approx e_{n}^{(0)}+{1\over{2}}(1-{g^{2}\over{\sqrt{D_{n}}}})\,d+{\mu+a_{n}\,r\,g\over{\sqrt{D_{n}}}}.$$

(21)

Actually, $e$-$h$ breaking $(d,r)$, listed in Eq. (13), has a sizable effect and shifts all the spectra, except for $\epsilon_{0}$, upwards appreciably [e.g., roughly by 10 % for $n=(2,3)$], as seen from Fig. 1, which, for later convenience, is presented in Sec. V. Unlike $e_{0}\sim e_{1}\sim O(\mu)$, $e_{n}$ depend on bias $u$ only weakly $\sim O(\mu/g^{2})$. The spectra $\epsilon_{n^{\prime}}$ of the higher branches $n^{\prime}=\pm 1^{\prime},\pm 2^{\prime},\dots$ show similar but partially different dependence on breaking $X$.

Each $|n|\geq 2$ is associated with a pair of electron and hole modes $\pm n$ (and $\pm n^{\prime}$). In contrast, the pseudo-zero modes $n=0$ and $n=1$ stand alone (per spin and valley) and are, in this sense, $e$-$h$ selfconjugate, with $\pm n\rightarrow n$ in Eqs. (14) and (15). Thus $\epsilon_{0}=-\epsilon_{0}|_{-X}=-\omega_{c}\mu$ and $\epsilon_{1}=-\epsilon_{1}|_{-X}$; that is, $\epsilon_{0}=\epsilon_{1}=0$ for $X\rightarrow 0$ and $\epsilon_{1}$ consists of odd powers of breaking $X$.

The Landau-level structure is made explicit by passing to the $|n,y_{0}\rangle$ basis (with $y_{0}\equiv\ell^{2}p_{x}$) via the expansion $\Psi^{a}_{\alpha}({\bf x})=\sum_{n,y_{0}}\langle{\bf x}|n,y_{0}\rangle\,\psi^{n;a}_{\alpha}(y_{0})$, where $n$ refers to the Landau level index, $\alpha\in(\uparrow,\downarrow)$ to the spin and $a\in(K,K^{\prime})$ to the valley. It is tacitly understood that the sum is taken over the higher branches $n^{\prime}=(\pm 1^{\prime},\pm 2^{\prime},\dots)$ as well. The one-body Hamiltonian is written as

$$H=\int dy_{0}\sum_{m,n}\sum_{a,\alpha}\psi_{\alpha}^{m;a{\dagger}}(y_{0})\,\epsilon_{n}^{a;\alpha}\,\psi_{\alpha}^{n;a}(y_{0}).$$

(22)

The charge density $\rho_{-{\bf p}}=\int d^{2}{\bf x}\,e^{i{\bf p\cdot x}}\,\rho$ with $\rho=(\Psi^{K})^{{\dagger}}\Psi^{K}+(\Psi^{K^{\prime}})^{{\dagger}}\Psi^{K^{\prime}}$ is thereby written as KS_CR

	$\displaystyle\rho_{-{\bf p}}$	$\displaystyle=$	$\displaystyle\gamma_{\bf p}\sum_{m,n=-\infty}^{\infty}\sum_{a,\alpha}g^{mn;a}_{\bf p}\,R^{mn;aa}_{\alpha\alpha;\bf-p},$
	$\displaystyle R^{mn;ab}_{\alpha\beta;{\bf-p}}$	$\displaystyle\equiv$	$\displaystyle\int dy_{0}\,{\psi^{m;a}_{\alpha}}^{{\dagger}}(y_{0})\,e^{i{\bf p\cdot r}}\,\psi^{n;b}_{\beta}(y_{0}),$		(23)

where $\gamma_{\bf p}=e^{-\ell^{2}{\bf p}^{2}/4}$; ${\bf r}=(i\ell^{2}\partial/\partial y_{0},y_{0})$ stands for the center coordinate with uncertainty $[r_{x},r_{y}]=i\ell^{2}$. The charge operators $R^{mn;aa}_{\alpha\alpha;\bf p}$ obey the $W_{\infty}$ algebra GMP .

The coefficient matrix $g^{mn;a}_{\bf p}\equiv g^{mn}_{\bf p}|^{a}$ in valley $a\in(K,K^{\prime})$ is constructed from the eigenvectors ${\bf b}_{n}|^{a}$,

	$\displaystyle g^{mn}_{\bf p}$	$\displaystyle=$	$\displaystyle c_{m}\,c_{n}\,f_{\bf p}^{|m|,|n|}+b^{\prime}_{m}\,b^{\prime}_{n}\,f_{\bf p}^{|m|-2,|n|-2}$		(24)
			$\displaystyle+(b_{m}\,b_{n}+c^{\prime}_{m}\,c^{\prime}_{n})\,f_{\bf p}^{|m|-1,|n|-1},$

where

$$f^{mn}_{\bf p}=\sqrt{{n!\over{m!}}}\,\Big{(}{-\ell p^{{\dagger}}\over{\sqrt{2}}}\Big{)}^{m-n}\,L^{(m-n)}_{n}\Big{(}{1\over{2}}\ell^{2}{\bf p}^{2}\Big{)}$$

(25)

for $m\geq n\geq 0$, and $f^{nm}_{\bf p}=(f^{mn}_{\bf-p})^{{\dagger}}$; $f^{mn}_{\bf p}=0$ for $m<0$ or $n<0$; $p=p_{x}\!+i\,p_{y}$. In view of Eqs. (14) and  (15), $g^{mn}_{\bf p}$ have the following property under $e$-$h$ conjugation,

	$\displaystyle g^{-m,-n;a}_{\bf p}=g^{m,n;a}_{\bf p}|_{-\mu,-d,-r},$
	$\displaystyle g^{mn}_{\bf p}|^{K^{\prime}}=g^{m,n}_{\bf p}|^{K}_{-\mu}=g^{-m,-n}_{\bf p}|^{K}_{-d,-r}.$		(26)

For $n,m\in(0,1)$, $g^{m,n;a}_{\bf p}=g^{m,n;a}_{\bf p}|_{-X}$ are even functions of breaking $X$. They actually take simple form

	$\displaystyle g^{00}_{\bf p}$	$\displaystyle=$	$\displaystyle 1,\ \ g^{11}_{\bf p}=1-c_{1}^{2}\,\textstyle{1\over{2}}\ell^{2}{\bf p}^{2},$
	$\displaystyle g^{01}_{\bf p}$	$\displaystyle=$	$\displaystyle c_{1}\ell\,p/\sqrt{2},\ \ g^{10}_{\bf p}=-c_{1}\ell\,p^{{\dagger}}/\sqrt{2},$		(27)

in each valley, with $c_{1}|^{K^{\prime}}=c_{1}|^{K}_{-\mu}$.

The Coulomb interaction $V={1\over{2}}\sum_{\bf p}v_{\bf p}\,{:\!\rho_{\bf-p}\,\rho_{\bf p}\!:}$ is written as

$$V={1\over{2}}\sum_{\bf p}v_{\bf p}\,\gamma_{\bf p}^{2}\,g^{jk;a}_{\bf p}\,g^{mn;b}_{\bf-p}:\!R^{jk;aa}_{\alpha\alpha;{\bf-p}}\,R^{mn;bb}_{\beta\beta;{\bf p}}\!:,$$

(28)

with the potential $v_{\bf p}=2\pi\alpha_{e}/(\epsilon_{\rm b}|{\bf p}|)$, $\alpha_{e}\equiv e^{2}/(4\pi\epsilon_{0})$ and the substrate dielectric constant $\epsilon_{\rm b}$; $\sum_{\bf p}\equiv\int d^{2}{\bf p}/(2\pi)^{2}$; $:\ :$ denotes normal ordering. For simplicity, we ignore a small effect of interlayer separation. Here and from now on, we suppress summations over levels $n$, valleys $a$ and spins $\alpha$, with the convention that the sum is taken over repeated indices. The one-body Hamiltonian $H$ is thereby written as

$$H=\epsilon_{n}^{a;\alpha}\,R^{nn;aa}_{\alpha\alpha;{\bf p=0}}-\mu_{\rm Z}\,(T_{3})_{\alpha\beta}R^{nn;aa}_{\alpha\beta;{\bf p=0}}.$$

(29)

Here the Zeeman term $\mu_{\rm Z}\equiv g^{*}\mu_{\rm B}B$ is introduced via the spin matrix $T_{3}=\sigma_{3}/2$.

In this section we study Coulombic contributions to the Landau-level and associated CR spectra. In graphene, unlike conventional quantum Hall systems, the electrons and holes are always subject to quantum fluctuations of the infinitely-deep filled valence band (or the Dirac sea), which, being strong, lead to ultraviolet divergences. One first has to handle them properly.

The Coulomb direct interaction leads to a divergent self-energy $\propto v_{\bf p\rightarrow 0}$, which, as usual, is removed when a neutralizing background is taken into account. The exchange interaction, on the other hand, gives rise to corrections to level spectra $\epsilon_{n}^{a;\alpha}$ of the form

$$\Delta\epsilon_{n}^{a;\alpha}=-\sum_{\bf p}v_{\bf p}\gamma_{\bf p}^{2}\,\sum_{m}\nu_{m}^{a;\alpha}\,|g^{nm;a}_{\bf p}|^{2}.$$

(30)

Here $0\leq\nu_{n}^{a;\alpha}\leq 1$ stands for the filling fraction of the $(n,a,\alpha)$ level. The exchange interaction preserves the spin and valley $(\alpha,a)$. Accordingly, from now on, we suppress them and mainly refer to valley $K$.

The self-energies $\Delta\epsilon_{n}$ involve a sum over infinitely many filled levels $m$ in the valence band. Their structure is clarified if one notes the completeness relation

$$\sum_{k=-\infty}^{\infty}|g^{nk}_{\bf p}|^{2}=1/\gamma_{\bf p}^{2}=e^{{1\over{2}}\ell^{2}{\bf p}^{2}}.$$

(31)

Actually, this follows from the fact KS_LWGS that $G^{nk}_{\bf p}\equiv\gamma_{\bf p}\,g^{nk}_{\bf p}$ are ($W_{\infty}$) unitary matrices that obey the composition law $G_{\bf p}\,G_{\bf q}=e^{i{1\over{2}}\ell^{2}{\bf p}\times{\bf q}}\,G_{\bf p+q}$ with ${\bf p}\times{\bf q}=p_{x}q_{y}-p_{y}q_{x}$. The half-infinite sum in $\Delta\epsilon_{n}$ is then cast in the form

	$\displaystyle\gamma_{\bf p}^{2}\!\sum_{k\leq-2}|g^{nk}_{\bf p}|^{2}\!$	$\displaystyle=$	$\displaystyle\!{\textstyle{1\over{2}}}-{\textstyle{1\over{2}}}\,F_{n}(z)-{\textstyle{1\over{2}}}\gamma_{\bf p}^{2}\,(|g^{n0}_{\bf p}|^{2}+|g^{n1}_{\bf p}|^{2}),$
	$\displaystyle F_{n}(z)$	$\displaystyle\equiv$	$\displaystyle\gamma_{\bf p}^{2}\sum_{k=2}^{\infty}\{|g^{nk}_{\bf p}|^{2}-|g^{n,-k}_{\bf p}|^{2}\},$		(32)

where $z={1\over{2}}\ell^{2}{\bf p}^{2}$. Here again it is tacitly understood that the sum is taken over $k^{\prime}\leq-1$ or $k^{\prime}\geq 1$ as well.

The self-energies $\Delta\epsilon_{n}$ are thereby rewritten as

$$\Delta\epsilon_{n}=\sum_{\bf p}v_{\bf p}\Big{[}\!-{\textstyle{1\over{2}}}+{\textstyle{1\over{2}}}\,F_{n}(z)-\sum_{k}\nu[k]\,\gamma_{\bf p}^{2}|g^{nk}_{\bf p}|^{2}\Big{]}.$$

(33)

Here the last term with the $``$electron-hole” filling factor,

	$\displaystyle\nu[k]$	$\displaystyle=$	$\displaystyle\nu_{k}\,\theta_{(k\geq 2)}-(1-\nu_{k})\theta_{(k\leq-2)}$		(34)
			$\displaystyle+(\nu_{1}-{\textstyle{1\over{2}}})\,\delta_{k,1}+(\nu_{0}-{\textstyle{1\over{2}}})\,\delta_{k,0},$

where $\theta_{(k\geq 2)}=1$ for $k\geq 2$ and $\theta_{(k\geq 2)}=0$ otherwise, etc., refers to a finite number of filled electron or hole levels around the PZM sector $n=\{0,1\}$.

Now $F_{n}(z)\equiv F_{n}(z;g,\mu,d,r)$ involve a sum over infinitely many levels. $F_{n}(z)\propto 1/\sqrt{z}$ for $z\rightarrow\infty$ and yield ultraviolet divergences upon integration over ${\bf p}$ with $v_{\bf p}$. In view of Eq. (26), $F_{\pm n}(z)$ are related in each valley or between the valleys,

	$\displaystyle F_{-n}(z)$	$\displaystyle=$	$\displaystyle-F_{n}(z)|_{-\mu,-d,-r},$
	$\displaystyle F_{n}(z)|^{K^{\prime}}$	$\displaystyle=$	$\displaystyle F_{n}(z)|^{K}_{-\mu}=-F_{-n}(z)|^{K}_{-d,-r}.$		(35)

For $n=(0,1)$, $F_{n}(z)=-F_{n}(z)|_{-X}$; $F_{0}(z)$ and $F_{1}(z)$ are odd in breaking $X=(u,d,r)$. The PZM sector thus has no divergence for $X=0$, $(F_{0}(z),F_{1}(z))|_{X\rightarrow 0}\rightarrow 0$.

Let us eliminate from $\Delta\epsilon_{n}$ a constant, $\sum_{\bf p}v_{\bf p}\,[-{1\over{2}}+{1\over{2}}F_{0}(z)|_{\mu=0}]$ common to all levels, which is safely done by adjusting zero of energy. Via such regularization, self-energies $\Delta\epsilon_{n}$ are cast in the form

$$\Delta\epsilon_{n}\stackrel{{\scriptstyle\rm reg}}{{=}}{\textstyle{1\over{2}}}\sum_{\bf p}v_{\bf p}\,{\cal F}_{n}(z)-\sum_{k}\nu[k]\sum_{\bf p}v_{\bf p}\,\gamma_{\bf p}^{2}|g^{nk}_{\bf p}|^{2}.$$

(36)

Note that

$${\cal F}_{n}(z)\equiv F_{n}(z)-F_{0}(z)|_{\mu=0}.$$

(37)

inherit the conjugation property of $F_{n}(z)$ in Eq. (35).

A key property of the electron-hole filling factor $\nu[k]$ defined in Eq. (34) is that it is odd under $e$-$h$ conjugation: $\nu[k]$ changes sign upon interchanging the electron and hole levels $(k\rightarrow-k)$ by replacing, in $\nu[-k]$, $\nu_{-k}\rightarrow 1-\nu_{k}$. Noting this feature and rewriting Eq. (36) then allows us to relate, as done earlier KS_CRnewGR for monolayer graphene, $\Delta\epsilon_{\pm n}$ in a valley or between the valleys. In terms of the full spectra to $O(V)$,

$$\hat{\epsilon}_{n}=\epsilon_{n}+\Delta\epsilon_{n},$$

(38)

the result is

$$\hat{\epsilon}_{n}^{K}|^{N_{\rm f}=m}=-\hat{\epsilon}_{-n}^{K}|^{N_{\rm f}=2-m}_{-\mu,-d,-r}=-\hat{\epsilon}_{-n}^{K^{\prime}}|^{N_{\rm f}=2-m}_{-d,-r}.$$

(39)

Here $N_{\rm f}$ specifies filling of the relevant valley; note that $\Delta\epsilon_{n}$ depend on filling. We assign $N_{\rm f}=0$ to the empty PZM sector (of a given spin and valley) with levels below it ($n\leq-2$) all filled, or with the $``$uppermost filled level” $n_{\rm f}=-2$; accordingly, $N_{\rm f}=(-1,0;3,4)$, e.g., refer to $n_{\rm f}=(-3,-2;2,3)$. $N_{\rm f}=1$ refers to the PZM sector with one filled level and $N_{\rm f}=2$ to the filled sector. For definiteness, in what follows, we focus on cases of integer filling, in which a distinct band gap is present, and take $N_{\rm f}$ in Eq. (39) to be integers. Let us now imagine, e.g., valley $K$ with $N_{\rm f}=m$. Interchanging electrons and holes (according to $\nu_{-k}\rightarrow 1-\nu_{k}$) yields a configuration with $N_{\rm f}=2-m$, and Eq. (39) implies that the spectra of such $e$-$h$ conjugate configurations are intimately related. A typical $e$-$h$ conjugate pair are the empty PZM sector ($N_{\rm f}=0$) and the filled one ($N_{\rm f}=2$).

Some care is needed in handling the $N_{\rm f}=1$ configuration, which is not uniquely fixed, e.g., $N_{\rm f}=1$ with $n_{\rm f}=0$ or $n_{\rm f}=1$ or with any linear combination of $\{0,1\}$. Note also that $e$-$h$ conjugation reverses the level layout of the PZM sector, e.g., $(0,1)\rightarrow(1,0)$. As a result, if, e.g., $\hat{\epsilon}_{0}<\hat{\epsilon}_{1}$ in a valley, Eq. (39) reads

$$(\hat{\epsilon}_{0},\hat{\epsilon}_{1})|^{N_{\rm f}=1,n_{\rm f}=0}=(-\hat{\epsilon}_{0},-\hat{\epsilon}_{1})|^{N_{\rm f}=1,n_{\rm f}=1}_{-X},$$

(40)

which means that the filled $n=0$ level turns into the filled $n=1$ level in the conjugate configuration. In Sec. V, we study the PZM sector over a continuous range $0\leq N_{\rm f}\leq 2$ and encounter an interesting case of mixed $\{0,1\}$ levels at $N_{\rm f}=1$.

Let us next study CR, namely, optical interlevel transitions at zero momentum transfer, with the selection rule AbF $\Delta|n|=\pm 1$ for graphene, i.e., (i) intraband channels $\{n\leftarrow n-1,-(n-1)\leftarrow-n\}$ and (ii) interband channels $T_{n}\equiv\{n\leftarrow-(n-1),\ n-1\leftarrow-n\}$ for $n=1,2,\cdots$; $T_{1}=\{1\!\leftrightarrow\!0\}$, $T_{2}=\{2\leftarrow 1,\ 1\leftarrow-2\}$, $T_{3}=\{3\leftarrow-2,\ 2\leftarrow-3\}$, etc. Interband CR is specific to Dirac electrons. Consider now CR from level $j$ to level $n$ for each (valley, spin)=$(a,\alpha)$ channel and denote the associated excitation energy as $\epsilon_{\rm exc}^{n\leftarrow j}=\epsilon_{n}-\epsilon_{j}+\Delta\epsilon^{n,j}$; CR preserves the valley and spin so that $(a,\alpha)$ will be suppressed from now on. In the mean-field treatment KS_CR ; KH ; MOG , the corrections

	$\displaystyle\Delta\epsilon^{n,j}$	$\displaystyle=$	$\displaystyle\Delta\epsilon_{n}-\Delta\epsilon_{j}-(\nu_{j}-\nu_{n})\,{\cal A}_{n,j},$
	$\displaystyle{\cal A}_{n,j}$	$\displaystyle=$	$\displaystyle\sum_{\bf p}v_{\bf p}\gamma_{\bf p}^{2}\,g^{nn}_{\bf-p}\,g^{jj}_{\bf p},$		(41)

consist of the self-energy difference and the Coulombic attraction ${\cal A}_{n,j}$ between the excited electron and created hole. The full CR spectra thus differ from the dressed level gaps by attraction energies ${\cal A}_{n,j}$,

$$\epsilon_{\rm exc}^{n\leftarrow j}=\hat{\epsilon}_{n}-\hat{\epsilon}_{j}-(\nu_{j}-\nu_{n})\,{\cal A}_{n,j}.$$

(42)

It will be clear now that, via $e$-$h$ conjugation, $\epsilon_{\rm exc}^{n\leftarrow-j}$ turns into $\epsilon_{\rm exc}^{j\leftarrow-n}$. In particular, the interband CR channels $T_{n}\equiv\{n\leftarrow-(n-1),\ n-1\leftarrow-n\}$ are intimately related,

	$\displaystyle\epsilon_{\rm exc}^{n\leftarrow-(n-1)}|^{K;N_{\rm f}=m}\!$	$\displaystyle=$	$\displaystyle\epsilon_{\rm exc}^{n-1\leftarrow-n}|^{K;N_{\rm f}=2-m}_{-\mu,-d,-r},$		(43)
		$\displaystyle=$	$\displaystyle\epsilon_{\rm exc}^{n-1\leftarrow-n}|^{K^{\prime};N_{\rm f}=2-m}_{-d,-r};$

we suppose $N_{\rm f}\not=1$ here. Such an $``$$e$-$h$ conjugate” character of $T_{n}$ was noticed earlier KS_CRnewGR for monolayer graphene, for which $e$-$h$ conjugation is an exact symmetry. Here for bilayer graphene $e$-$h$ conjugation, though not a symmetry, is an operation that tells us how the CR spectra deviate from the symmetric ones by breaking $(\mu,d,r)$. Experimentally the conjugate channels of a given set $T_{n}$ are observable as competing signals over a finite range of filling factor $\nu$ and are indistinguishable unless polarized light is used.

The self-energies $\Delta\epsilon_{n}$ are afflicted with ultraviolet divergences. In this section we consider how to extract physically observable information out of them. To this end, one first defines renormalized parameters by setting $v=Z_{v}v^{\rm ren}=v^{\rm ren}+\delta v$, $\gamma_{1}=\gamma_{1}^{\rm ren}+\delta\gamma_{1}$, $\Delta=\Delta^{\rm ren}+\delta\Delta$, etc., with counterterms $\delta v=(Z_{v}-1)\,v^{\rm ren}\sim O(\alpha_{e})$, $\delta\gamma_{1}$, etc. The one-body Hamiltonian $H=H^{\rm ren}+\delta_{\rm ct}H$ is then divided into the portion $H^{\rm ren}$ involving only the renormalized quantities and the counterterms $\delta_{\rm ct}H$. One starts with $H^{\rm ren}$, calculates the Coulombic corrections and encounters divergences. If one could remove them by adjusting $\delta_{\rm ct}H$, renormalization is done properly.

Fortunately, the renormalization procedure to $O(V)\sim O(\alpha_{e})$ for bilayer graphene in a magnetic field $B$ has been formulated earlier KS_CR_eh . The key point is that, since $B$ simply acts as a long-wavelength cutoff $\sim\ell$, the short-distance structure of the present theory is known from its free-space ($B=0$) version: (i) The divergence associated with velocity $v$ is the same as that for monolayer graphene,

$$\delta v=-(\alpha_{e}/8\epsilon_{b})\log(\ell\Lambda)^{2},$$

(44)

with a momentum cutoff $\Lambda$. (ii) $v_{4}$ and $u$ remain finite,

$$\delta v_{4}=\delta u=0.$$

(45)

We thus take $r=v_{4}/v^{\rm ren}$ and $\mu=(u/2)/\omega_{c}^{\rm ren}$ finite; $\omega_{c}^{\rm ren}\equiv\sqrt{2}v^{\rm ren}/\ell$. (iii) $\gamma_{1}$ and $\Delta$, though mixed under renormalization, are also governed by $\delta v$,

	$\displaystyle\delta\gamma_{1}$	$\displaystyle=$	$\displaystyle(\gamma_{1}^{\rm ren}+r\,\Delta^{\rm ren})\,h(r)\,\delta v/v^{\rm ren},$
	$\displaystyle\delta\Delta$	$\displaystyle=$	$\displaystyle 2\,(\Delta^{\rm ren}+r\,\gamma_{1}^{\rm ren})\,h(r)\,\delta v/v^{\rm ren},$		(46)

where $h(r)=1/(1-r^{2})$.