Symmetry characterization of the collective modes of the phase diagram of the $\nu=0$ quantum Hall state in graphene: Mean-field and spontaneously broken symmetries

TL;DR

This paper analyzes the phase diagram and collective modes of the $ u=0$ quantum Hall state in graphene, revealing symmetry-breaking mechanisms, collective mode structures, and the robustness of certain symmetries under various conditions.

Contribution

It provides a comprehensive classification of collective modes and demonstrates the persistence of $SO(5)$ symmetry in the mean-field phase diagram of graphene's quantum Hall states.

Findings

Phase transitions involve singlet orbital pseudospin modes independent of Coulomb strength.

Goldstone modes reveal valley-spin symmetry breaking mechanisms.

$SO(5)$ symmetry persists at phase boundaries under certain conditions.

Abstract

We devote this work to the study of the mean-field phase diagram of the $\nu=0$ quantum Hall state in bilayer graphene and the computation of the corresponding neutral collective modes, extending the results of recent works in the literature. Specifically, we provide a detailed classification of the complete orbital-valley-spin structure of the collective modes and show that phase transitions are characterized by singlet modes in orbital pseudospin, which are independent of the Coulomb strength and suffer strong many-body corrections from short-range interactions at low momentum. We describe the symmetry breaking mechanism for phase transitions in terms of the valley-spin structure of the Goldstone modes. For the remaining phase boundaries, we prove that the associated exact $SO(5)$ symmetry existing at zero Zeeman energy and interlayer voltage survives as a weaker mean-field symmetry of the Hartree-Fock equations. We extend the previous results for bilayer graphene to the monolayer scenario. Finally, we show that taking into account Landau level mixing through screening does not modify the physical picture explained above.

Equations

$$\hat{\psi }(\mathbf{x})$$

$${\hat{\psi }}_{\xi }(\mathbf{x})$$

$${\hat{H}}_0=\int {\mathrm{d}}^2\mathbf{x}{\hat{\psi }}^{†}(\mathbf{x})\left[H_B+{ϵ}_V{T}_{zz}-{ϵ}_Z{\sigma }_z\right]\hat{\psi }(\mathbf{x}),$$

$$H_B$$

$${\omega }_B$$

$$a_B=\frac{l_B}{\hslash }\frac{{\pi }_y+i{\pi }_x}{\sqrt{2}},l_B=\sqrt{\frac{\hslash }{e{B}_{\perp }}}=\frac{25.7}{\sqrt{B_{\perp }[\text{T}]}}\text{nm}.$$

$${\Psi }_{n,k}^0(\mathbf{x})$$

$${ϵ}_n$$

$$\Psi^{0}_{n,k}(\mathbf{x})=\frac{e^{iky}}{\sqrt{L_{y}}}\left[\begin{array}[]{c}0\\ \phi_{|n|}(x+kl^{2}_{B})\end{array}\right],~{}\epsilon_{n}=0$$

$$\hat{\psi }(\mathbf{x})=\text{sideset}{}^{\prime }{\sum }_{n=-\infty }^{\infty }\sum _{k,\alpha }{\Psi }_{n,k,\alpha }^0(\mathbf{x}){\hat{c}}_{n,k,\alpha }$$

$${\hat{H}}_C=\frac{1}{2}\int {\mathrm{d}}^2\mathbf{x}{\mathrm{d}}^2{\mathbf{x}}^{\prime }:[{\hat{\psi }}^{†}(\mathbf{x})\hat{\psi }(\mathbf{x})]V_0(\mathbf{x}-{\mathbf{x}}^{\prime })[{\hat{\psi }}^{†}({\mathbf{x}}^{\prime })\hat{\psi }({\mathbf{x}}^{\prime })]:$$

$${\hat{H}}_{sr}=\text{sideset}{}^{\prime }{\sum }_{i,j}\frac{1}{2}{g}_{ij}\int {\mathrm{d}}^2\mathbf{x}[{\hat{\psi }}^{†}(\mathbf{x})T_{ij}\hat{\psi }(\mathbf{x}){]}^2:$$

$$\frac{{ϵ}_Z}{\hslash {\omega }_B}=0.014\frac{B}{B_{\perp }}\ll 1$$

$$\frac{F_C}{\hslash {\omega }_B}=\frac{13.58}{\kappa \sqrt{B_{\perp }[\text{T}]}},F_C\equiv \frac{e_c^2}{\kappa l_B}$$

$$\frac{e_c^{2}d}{\kappa l_B^2\hslash {\omega }_B}\sim \frac{F_C}{\hslash {\omega }_B}\frac{d}{l_B}\sim \frac{0.1}{\kappa }\ll 1$$

$${\hat{H}}^{(0)}$$

$$\hat{\psi }(\mathbf{x})=\sum _{n=0,1}\sum _{k,\alpha }{\Psi }_{n,k,\alpha }^0(\mathbf{x}){\hat{c}}_{n,k,\alpha }$$

$${ϵ}_{n,\alpha }{\Psi }_{n,k,\alpha }(\mathbf{x})$$

$$|{\Psi }_0⟩=\prod _p{\hat{c}}_{0,p,a}^{†}{\hat{c}}_{1,p,a}^{†}{\hat{c}}_{0,p,b}^{†}{\hat{c}}_{1,p,b}^{†}|DS⟩,$$

$${ϵ}_{n,\alpha }{\chi }_{\alpha }$$

$$F_{00}=\sqrt{\frac{\pi }{2}}{F}_C,F_{01}=F_{10}=\frac{1}{2}{F}_{00},F_{11}=\frac{3}{4}{F}_{00}$$

$$F_0=\frac{3}{2}{F}_{00},F_1=\frac{5}{4}{F}_{00}$$

$${ϵ}_{n,(a,b)}=-\frac{F_n}{2}+{ϵ}_{(a,b)},{ϵ}_{n,(c,d)}=\frac{F_n}{2}+{ϵ}_{(c,d)}$$

$$E_{HF}$$

$$E(P)$$

$$V=({ϵ}_Z^{\ast },{ϵ}_V^{\ast })=(-2u_{\perp },u_z-u_{\perp }).$$

$${\chi }_a$$

$${\chi }_c$$

$${ϵ}_{a,b}$$

$$P$$

$${\chi }_a$$

$${\chi }_c$$

$${ϵ}_{a,b}$$

$$P$$

$${ϵ}_V-{ϵ}_Z=u_{\perp }+u_z$$

Full text

I Introduction

Since the discovery of the integer Klitzing et al. (1980); Tsui et al. (1982a) and fractional quantum Hall (QH) effects Tsui et al. (1982b); Stormer et al. (1983); Laughlin (1983), a large number of works have been devoted to the study of two-dimensional (2D) systems in the presence of strong perpendicular magnetic fields. Due to its chiral character, rich valley-spin structure and relatively large cyclotron frequency, graphene is a particulary interesting scenario to test QH features.

In this way, the study of QH states in graphene has been a hot topic of research in the last years McCann and Fal’ko (2006); Yang et al. (2006); Iyengar et al. (2007); Doretto and Smith (2007); Bychkov and Martinez (2008); Jung and MacDonald (2009); Abanin et al. (2009); Côté et al. (2010); Shizuya (2010); Gorbar et al. (2010, 2011, 2012); Kim et al. (2011); Shizuya (2011, 2012); Kharitonov (2012a, b, c, d); Lambert and Côté (2013); Tőke (2013); Wu et al. (2014); Knothe and Jolicoeur (2015); Murthy et al. (2016); Tikhonov et al. (2016); Knothe and Jolicoeur (2016); Jia et al. (2017). Among all the possible QH states, the $\nu=0$ QH state, corresponding to the charge neutrality point, has received special attention due to its intriguing strongly insulating behavior in both bilayer and monolayer graphene at zero in-plane magnetic field. In particular, the bilayer provides a more interesting scenario because of the richer structure of its zero-energy Landau level (LL) and the possibility of introducing an energy bias between the two valleys with the help of a perpendicular electric field. In fact, Ref. Kharitonov (2012b) presented a complete characterization of the mean-field phase diagram of the $\nu=0$ QH state in bilayer graphene, taking into account the most general spin symmetric interactions (including short-range valley/sublattice asymmetric interactions) and the introduction of a voltage between the two layers. The identified phases are ferromagnetic (F), canted anti-ferromagnetic (CAF), fully layer-polarized (FLP) and partially layer-polarized (PLP) Kharitonov (2012b). The resulting phase diagram is analog to that of monolayer graphene Kharitonov (2012a) due to the identical structure of the Hamiltonian governing short-range valley/sublattice asymmetric interactions.

Complementarily, the computation of the spectrum of the collective excitations in QH integer states has been also the subject of study of an important number of works Yang et al. (2006); Iyengar et al. (2007); Doretto and Smith (2007); Shizuya (2010); Tőke and Fal’ko (2011); Sári and Tőke (2013); Lambert and Côté (2013); Wu et al. (2014); Murthy et al. (2016). For instance, the spectrum of inter-LL collective excitations was obtained in Refs. Iyengar et al. (2007); Shizuya (2010); Sári and Tőke (2013). On the other hand, the spectrum of intra-LL excitations within the zero-energy Landau level was computed for the $\nu=0$ QH state Tőke and Fal’ko (2011) and the rest of integer QH states Lambert and Côté (2013) allowing only for long range and interlayer Coulomb interactions. In Ref. Wu et al. (2014), the dispersion relation of the modes for the Kekulé distortion (KD) and CAF phases of the $\nu=0$ QH state was computed for monolayer graphene using an effective low-energy model. Also within the $\nu=0$ QH state in the monolayer scenario, the bulk and edge collective modes of the CAF and F phases were obtained in Ref. Murthy et al. (2016).

We devote this work to the computation of the mean-field energies and the intra-LL collective modes of the whole phase diagram of the $\nu=0$ QH state of bilayer graphene, as presented in Ref. Kharitonov (2012b). Specifically, we compute the dispersion relation within the time-dependent Hartree-Fock approximation (TDHFA), considering the complete Hamiltonian of short-range interactions Kharitonov (2012b) and also including explicitly the interaction of the electrons with the filled Dirac sea Shizuya (2012), devoting special attention to the characterization of the rich orbital-valley-spin structure of the modes, induced by the short-range valley-asymmetric interactions, and to the study of phase transitions in terms of these modes. In this way, the work here presented provides a natural continuation to the recent literature on the field.

In particular, we show that only singlet modes in orbital pseudospin are strongly modified by many-body effects arising from short-range interactions while the remaining modes are almost unaffected by them due to the dominant character of long range Coulomb interactions. Indeed, as these orbital-singlet modes are independent of the Coulomb interaction strength at low momentum, they also represent the lowest-energy neutral excitations of the system, playing a crucial role in understanding the stability of the different phases. Regarding the valley-spin structure of the modes, we classify them according to the conserved valley and spin quantum numbers of the total Hamiltonian. We find that, while the FLP-PLP transition is governed by a spin-singlet mode, the F-CAF transition is governed by a valley-triplet one; this contrast arises due to the different nature of the spontaneous symmetry breaking mechanism of the CAF and PLP phases. Interestingly, we prove that the remaining phase boundaries, F-FLP and CAF-PLP, present a gapless mode arising from a mean-field symmetry inherited from the full exact $SO(5)$ symmetry existing at the same boundaries for zero Zeeman energy and zero interlayer voltage Wu et al. (2014). Moreover, we show that the CAF and PLP phases are able to present dynamical instabilities as a result of their spontaneously broken symmetries.

The strong analogy between the $\nu=0$ QH states in bilayer and monolayer graphene allows us to straightforwardly translate most of these results to the monolayer scenario, recovering essentially the same results of Ref. Murthy et al. (2016). We also study the effects of Landau level mixing by considering the screening effect of Coulomb interaction in the large-$N$ approximation Foster and Aleiner (2008); Basko and Aleiner (2008); Aleiner et al. (2007); Lemonik et al. (2010); Kharitonov (2012a), showing that the described orbital-valley-spin structure still holds, quantitatively (but not qualitatively) changing the dispersion relation of the modes. Finally, we relate the results presented in this work with experimental scenarios, including a computation of the mean-field transport gaps and a discussion on the collective modes detection.

We remark that the collective modes calculated in the present work are neutral and, therefore, topologically trivial Sondhi et al. (1993); Moon et al. (1995). The charged, topologically non-trivial excitations (skyrmions) have been studied in both monolayer Yang et al. (2006) and bilayer Abanin et al. (2009) graphene, and the effect of short-range interactions Wu et al. (2014) and screening Kharitonov (2012a) on the non-linear sigma-model stiffness coefficients has been addressed for monolayer graphene. Here we provide the corrections due to short range interactions to the stiffness coefficients in the bilayer case as well.

The article is arranged as follows: we first introduce the effective projected Hamiltonian considered in this work in Sec. II. We re-derive in Sec. III the phase diagram of Ref. Kharitonov (2012b) using a Hartree-Fock (HF) mean-field scheme and compute the corresponding mean-field energies and transport gaps. The dispersion relation of the different collective modes, computed within the TDHFA, is shown in Sec. IV. We translate the same calculations to monolayer graphene in Sec. V. Effects of LL mixing are studied in Sec. VI. A discussion on experimental features is presented in Sec. VII. Finally, the conclusions are drawn in Sec. VIII. Technical details about the diagonalization of the HF equations and the TDHFA are given in Appendices A-C.

II Effective Hamiltonian for the $\nu=0$ quantum Hall state in bilayer graphene

II.1 Low-energy Hamiltonian

In the first place, we briefly present the effective model used in this work for bilayer graphene, following Refs. McCann and Fal’ko (2006); Kharitonov (2012a, b, c); Sári and Tőke (2013), where the reader is referred for more details.

The effective dynamics at low energies can be described by a two-band model McCann and Fal’ko (2006), in which the field operator for electrons has eight components and reads:

[TABLE]

with $A,\tilde{B}$ the most far apart sublattices, $K,K^{\prime}$ the two valleys and $\xi$ the spin polarization. We note that the two sublattices are interchanged in the $K^{\prime}$ valley so, as usually done, we will refer to the corresponding subspace as $\bar{A}\bar{B}$ in order to avoid confusions. The eight components of the field operator then correspond to the total space $KK^{\prime}\otimes\bar{A}\bar{B}\otimes s$, $s$ being the spin space.

The corresponding effective Hamiltonian of the system is $\hat{H}=\hat{H}_{0}+\hat{H}_{C}+\hat{H}_{sr}$. After neglecting trigonal warping effects and other small corrections McCann and Fal’ko (2006); Kharitonov (2012c); Sári and Tőke (2013); Tőke (2013); Côté and Barrette (2013), the single-particle Hamiltonian $\hat{H}_{0}$ reads:

[TABLE]

where $T_{ij}=\tau_{i}^{KK^{\prime}}\otimes\tau_{j}^{\bar{A}\bar{B}}\otimes\hat{\mathbf{1}}^{s}$ and $i,j=0,x,y,z$ with $\tau_{i}$, $i=x,y,z$, the usual Pauli matrices and $\tau_{0}=\hat{\mathbf{1}}$ while $\sigma_{z}$ is the corresponding Pauli matrix in the spin space. In the following, the Pauli matrices in valley or sublattice space are denoted using the letter $\tau$ and the Pauli matrices in spin space are denoted using the letter $\sigma$.

The first term between square brackets in Eq. (13), $H_{B}$, is a $2\times 2$ matrix acting in the $\bar{A}\bar{B}$ subspace and corresponds to the kinetic energy:

[TABLE]

$m$ being the effective mass for which we take the experimental value $m=0.028m_{e}$, with $m_{e}$ the electron mass Mayorov et al. (2011), and $a_{B}$ the magnetic annihilation operator

[TABLE]

The operators $\pi_{i}=-i\hbar\partial_{i}+eA_{i}$, $i=x,y$, are the momentum components after the Peierls substitution Hofstadter (1976) and the magnitude $l_{B}$ is referred as the magnetic length.

The second term, $\epsilon_{V}T_{zz}$, arises from a voltage difference between the two layers, $\epsilon_{V}=Ea_{z}/2$, with $E$ the perpendicular electric field and $a_{z}\approx 0.35~{}\text{nm}$ the separation between the layers, while the third term takes into account the Zeeman effect, $\epsilon_{Z}=\mu_{B}B,B=\sqrt{B_{\parallel}^{2}+B_{\perp}^{2}}$. Here, we assume that the total magnetic field is not necessarily along the $z$ direction, i.e., it can present a parallel component to the graphene plane, $\mathbf{B}_{\parallel}$. The polarizations $\xi=\pm$ correspond to the spin components that are antiparallel (parallel) to the total magnetic field $\mathbf{B}$, respectively.

Taking into account that only the perpendicular magnetic field affects the orbital motion and using the Landau gauge, we can write the potential vector as $\mathbf{A}(\mathbf{x})=[0,B_{\perp}x,0]$. In this particular gauge, the eigenstates of $H_{B}$ are characterized by the following quantum numbers: the magnetic index $n$, which is an integer number and characterizes the energy of the corresponding Landau level, $\epsilon_{n}$; the momentum in the $y$-direction, $k$, and the polarization in the valley-spin ($KK^{\prime}\otimes s$) space, $\alpha$. Specifically, they are given by $\Psi^{0}_{n,k,\alpha}(\mathbf{x})=\Psi^{0}_{n,k}(\mathbf{x})\chi_{\alpha}$, where $\chi_{\alpha}$ is an arbitrary 4-component spinor in valley-spin space while the orbital wave function with components in the space $\bar{A}\bar{B}$ is

[TABLE]

for $|n|\neq 0,1$ and

[TABLE]

for the degenerate levels $|n|=0,1$ with zero energy (note that $n=\pm 1$ are indeed the same state). Hereafter, we refer to this manifold of states as the zero Landau level (ZLL). In the previous equations, $L_{y}$ is the length of the system in the $y$ direction and $\phi_{n}(x)$ is the usual harmonic oscillator wave function [see Eq. (115)]. We see that the kinetic energy is degenerate in $y$-momentum and valley-spin polarization; in particular, the degeneracy in $k$ for each magnetic level $n$ is $N_{B}=S/2\pi l^{2}_{B}$, $S$ being the total area of the system. Interestingly, the wave functions in the ZLL only have non-vanishing components in the subspace $KK^{\prime}\otimes\bar{B}\otimes s$.

Using the previous eigenfunctions, the field operator is decomposed as

[TABLE]

where ′ means that $n$ takes every integer value except $n=-1$.

With respect to the interacting part of the Hamiltonian, $\hat{H}_{C}$ represents the long range Coulomb interaction:

[TABLE]

Here, $:$ denotes normal ordering of the field operators and $V_{0}(\mathbf{x})=e^{2}_{c}/\kappa|\mathbf{x}|$ is the Coulomb potential, with $e^{2}_{c}\equiv e^{2}/4\pi\epsilon_{0}$ and $\kappa$ the dielectric constant of the environment. Finally, for the short-range interaction Hamiltonian, $\hat{H}_{sr}$, we consider the most general expression compatible with all the symmetries of the problem Kharitonov (2012a)

[TABLE]

The ′ in this sum denotes that we exclude the symmetric term $i=j=0$, already accounted by the long-range Coulomb interaction. These interactions are asymmetric in the valley and sublattice spaces. The origin of these short-range interactions are the Coulomb interaction between sublattice/valley spaces and the electron-phonon interactions, which we also treat as short-ranged Kharitonov (2012a). For shortness, we refer in the following to the long range Coulomb interactions as simply Coulomb interactions. Except for the kinetic energy term, this Hamiltonian is formally similar to that of monolayer graphene, see Sec. V.

The two-band model is expected to work quite well in a wide range of magnetic fields $1~{}\text{T}\lesssim B_{\perp}\lesssim 30~{}\text{T}$, specially for the ZLL Côté and Barrette (2013). For lower magnetic fields trigonal warping effects become important and for larger magnetic fields one has to use the complete four-band model, as the overlap of the wave function of the $n=1$ level with the ignored sublattices is not negligible. Other small corrections to the two-band Hamiltonian here considered break the degeneracy between the single-particle energies of the $n=0$ and $n=1$ levels Lambert and Côté (2013) but they are negligible compared to the analog Lamb shift, discussed in the next section.

II.2 Projection onto the zero Landau level

We now address the study of the $\nu=0$ QH state, which corresponds to a half-filling of the ZLL and complete filling of all LLs with $n\leq-2$. For that purpose, we make an estimation of the order of magnitude of the different terms in the Hamiltonian and we compare them with the typical energy difference between LLs, $\hbar\omega_{B}$. For the Zeeman term, we find

[TABLE]

for realistic values of the ratio $B/B_{\perp}$ while for the Coulomb interaction one has

[TABLE]

Usual values for the perpendicular magnetic field are $B_{\perp}\gtrsim 1~{}\text{T}$ and the highest available continuous magnetic field in the laboratory is $B_{\perp}\simeq 80~{}\text{T}$ Yang et al. (2016), which means that the strength of the Coulomb interaction verifies $F_{C}\gtrsim\hbar\omega_{B}$ when the environment is vacuum ($\kappa=1$). The interlayer voltage can in principle vary in a wide range Lambert and Côté (2013); Côté and Barrette (2013) but we keep its value sufficiently small here, $\epsilon_{V}\ll\hbar\omega_{B}$. Finally, a dimensional analysis of the short-range terms gives an estimation for the coupling constants $g_{ij}\sim e^{2}_{c}d/\kappa$ with $d\sim 0.1~{}\text{nm}$ the typical length scale of the lattice. Then, the energy scale associated to the short-range interactions is $\sim\frac{e^{2}_{c}d}{\kappa l^{2}_{B}}$, which is small compared to the LL separation as

[TABLE]

We note that the energy associated to the short-range interactions scales linearly with the magnetic field although for low values of the magnetic field the coupling constants themselves can be renormalized Kharitonov (2012a, b), see also Sec. VI.

As a result of the above analysis, the only term that is not small compared to the separation between LLs is that related with Coulomb interactions. We begin by treating Coulomb interactions as weak, $F_{C}\ll\hbar\omega_{B}$; for instance, by supposing a typical value $\kappa=5$ and a magnetic field $B\sim~{}20~{}T$, $F_{C}/\hbar\omega_{B}\sim 0.5$, which can be regarded as sufficiently small to treat it perturbatively. In Sec. VI, we address the usual situation $F_{C}\gtrsim\hbar\omega_{B}$ and explain how to deal with it.

Taking into account the previous considerations and in order to study the lowest energy excitations, we neglect LL mixing and restrict ourselves to the ZLL by projecting the full Hamiltonian into that subspace, as previously done in Refs. Tőke and Fal’ko (2011); Kharitonov (2012a, b); Lambert and Côté (2013). We remark [see Eq. (22) and ensuing discussion] that the states in the ZLL belong to the $KK^{\prime}\otimes\bar{B}\otimes s$ space, which means that they are localized, for each valley, in one sublattice or the other and correspondingly, in one layer or the other. Thus, within the ZLL, the sublattice degree of freedom becomes equivalent to the valley degree of freedom. The resulting effective Hamiltonian for the ZLL is:

[TABLE]

where $\tau_{i}\equiv\tau^{KK^{\prime}}_{i}$ and $g_{i}=g_{i0}+g_{iz},i=1,2,3$. As we restrict to the ZLL, the kinetic energy term is suppressed. We have neglected the symmetric short-range interaction, arising from the coupling $g_{0z}$, due to its smallness compared to the symmetric Coulomb interaction Kharitonov (2012a). Using symmetry considerations, it can be proven that $g_{x}=g_{y}\equiv g_{\perp}$ Aleiner et al. (2007); Kharitonov (2012a), so there are only two independent coupling constants $g_{\perp},g_{z}$. The potential $V_{DS}(\mathbf{x},\mathbf{x^{\prime}})$ represents the mean-field interaction of the ZLL with the (inert) Dirac sea compound by all the occupied states with $n\leq-2$ Shizuya (2012, 2013); Tőke (2013). As shown in Appendix A.2, this potential is diagonal within the ZLL; its explicit expression is given by Eq. (156). The projection onto the ZLL leads to a field operator of the form

[TABLE]

III Hartree-Fock equations and mean-field phase diagram

In order to obtain the mean-field phase diagram of the $\nu=0$ QH state at zero temperature, we use the Hartree-Fock (HF) approximation for the self-consistent single-particle wave functions, which we denote as $\Psi_{n,k,\alpha}$. The corresponding HF equations for the Hamiltonian (29) read (we refer the reader to Appendix A for all the technical details)

[TABLE]

where the indices $n,m=0,1$ label the magnetic levels, $k,p$ are the momenta in the $y$-direction and $\alpha,\beta$ represent the polarization in the valley-spin space. We have made explicit that we are dealing with an integer QH state, so every orbital $p$ is filled in the same way and then, the occupation number $\nu_{n,p,\alpha}$ of each state solely depends on $n$ and $\alpha$, $\nu_{n,p,\alpha}=\nu_{n,\alpha}$. In particular, for the $\nu=0$ QH state, only half of the ZLL is filled. Thus, for each value of the $y$-momentum $k$, only four states of the eightfold space, formed by the $0,1$ magnetic states and the valley-spin degrees of freedom, are occupied. As usual, the direct (Hartree) term for the Coulomb interaction is suppressed by the positive charge background.

An important result is that the orbital part of the self-consistent HF wave functions is equal to that of the non-interacting wave functions, given in Eq. (22); see Appendix A for the proof. Hence, the only remaining task is to specify the spinors $\chi_{\alpha}$. In order to minimize the dominant Coulomb interaction, the electrons occupy in the same way the valley-spin space for the two magnetic levels, i.e., $\nu_{0,a}=\nu_{0,b}=\nu_{1,a}=\nu_{1,b}=1$, with $\chi_{a,b}$ two orthogonal spinors Abanin et al. (2009); Kharitonov (2012b) so the mean-field ground state is of the form

[TABLE]

with $\ket{DS}$ the Dirac sea formed by all occupied Landau levels with $n\leq-2$.

The four remaining unoccupied states of the ZLL are characterized by the spinors $\chi_{c,d}$. In this way, the occupation number only depends on the valley-spin polarization $\alpha$, $\nu_{n,\alpha}=\nu_{\alpha}$, taking the values $\alpha=a,b,c,d$ that correspond to an orthonormal basis of the valley-spin space. These spinors are computed after projecting the HF equations into the orbital part of the wave functions, obtaining closed algebraic equations in valley-spin space:

[TABLE]

where $u_{i}=g_{i}/\pi l^{2}_{B}$ and $F_{n}=F_{n0}+F_{n1}$, with:

[TABLE]

so $F_{00}>F_{01}>F_{11}$ and then $F_{0}>F_{1}$ as

[TABLE]

The values of the factors $F_{nm}$ are obtained by inserting the Coulomb potential in Eq. (137).

The term $F_{n}/2$ in Eq. (33) is the analog of the Lamb shift Shizuya (2012, 2013), arising from the interaction of the ZLL with the Dirac sea after the proper regularization of the Hamiltonian (see Appendix A.2). The term $-F_{n}P$ arises from the exchange Coulomb interaction, where the matrix $P$ is the projector onto the subspace formed by $\chi_{a,b}$, $P=\chi_{a}\chi^{\dagger}_{a}+\chi_{b}\chi^{\dagger}_{b}$, and hence $P\chi_{a,b}=\chi_{a,b},~{}P\chi_{c,d}=0$. The remaining terms are those related with the short-range interactions and the valley-spin part of the single-particle Hamiltonian. We see that the sole dependence on the magnetic level is through the Coulomb and Lamb-shift terms, while the other contributions to the energy only depend on the spinor $\chi_{\alpha}$. Note that, although the interaction with the Dirac sea favors the filling of the magnetic level $n=1$, it is still more energetically favorable to occupy in the same way the levels $n=0,1$ due to the exchange interaction. Indeed, the mean-field state (32) is an exact eigenvalue of the effective Hamiltonian (29) when short-range interactions are neglected. As Coulomb interactions dominate over short-range interactions, this mean-field solution is expected to provide a very good approximation to the actual ground state Yang et al. (2006); Kharitonov (2012a), showing the robustness of the formalism here considered.

As a result of the previous discussion, the HF energies can be written as:

[TABLE]

with $\epsilon_{\alpha}$ depending only on the polarization $\alpha$. The energy of the total state (per wave vector state) is:

[TABLE]

The contribution from Coulomb interaction to the total energy turns out to be degenerate and does not depend on the specific form of the occupied spinors $\chi_{a,b}$. The actual ground state of the system is determined by comparing the energies corresponding to all possible solutions to the HF equations and selecting that with lowest energy $E(P)$, in the same fashion of Refs. Kharitonov (2012a, b). Hence, the corresponding mean-field phase diagram for the $\nu=0$ QH state is the same of those references and is represented in Fig. 1. The different possible phases are ferromagnetic (F), canted anti-ferromagnetic (CAF), fully layer-polarized (FLP) and partially layer-polarized (PLP). The expected phase for the $\nu=0$ QH state of bilayer graphene for $\epsilon_{V}=0$ and zero in-plane component of the magnetic field is the CAF phase Kharitonov (2012b, c); Maher et al. (2013), which implies that $u_{z}>-u_{\perp}>0$. The exact values of these short-range energies remain unknown but their order of magnitude is $u_{z},|u_{\perp}|\sim 0.1\hbar\omega_{B}$ Kharitonov (2012d); Young et al. (2014). In the following, we treat them as phenomenological inputs for the theory. All phase boundaries intersect at the critical point

[TABLE]

The complete phase diagram can be explored experimentally by manipulating the in-plane component of the magnetic field (which modifies the value of the total magnetic field $B$ and hence the value of $\epsilon_{Z}$) or the layer voltage (which modifies the value of $\epsilon_{V}$). This fact can be checked by looking at right Fig. 1, where we represent the phase diagram as a function of the Zeeman and voltage energies, $(\epsilon_{Z},\epsilon_{V})$, for fixed values of $u_{\perp},u_{z}$ such that $u_{z}>-u_{\perp}>0$.

The phase diagram here presented describes the $\nu=0$ QH state for $B_{\perp}\gtrsim 1$ T Kharitonov (2012b); the complementary phase diagram in the remaining limit of very low perpendicular magnetic fields $B_{\perp}\ll 1$ T is studied in Ref. Jia et al. (2017).

We now briefly describe all the phases and give the expressions for the filled and empty valley-spin spinors, the matrix $P$ and the corresponding valley-spin energies $\epsilon_{\alpha}$ for each phase.

III.1 Ferromagnetic phase

In the F phase, all the electrons have the spin aligned with the magnetic field. The complete set of solutions involves the following 4 spinors in valley-spin space, with eigenvalues:

[TABLE]

Here, $\ket{\pm s_{z}}$ denotes the state with spin polarization $\xi=\pm$, $I$ the $4\times 4$ identity matrix and $\ket{n_{z}}=\ket{K},\ket{-n_{z}}=\ket{K^{\prime}}$. The F phase is always a solution of the HF equations although it does not always correspond to the actual ground-state.

III.2 Full layer-polarized phase

The FLP phase is the equivalent of the F phase but in valley pseudospin, $KK^{\prime}$, which means that all the electrons are concentrated on one layer (due to the equivalency of valley-sublattice-layer in the ZLL in bilayer graphene). The corresponding spinors are now:

[TABLE]

In analogy to the F phase, it is always a solution to the HF equations but not necessarily the ground-state. By comparing the energies of the two phases, we obtain the boundary between the F and FLP phases:

[TABLE]

III.3 Canted anti-ferromagnetic phase

The previous phases would be the only possible phases if there were not short-range interactions, i.e., $u_{\perp}=u_{z}=0$. However, when taking into account these interactions, the system can exhibit canted anti-ferromagnetism or partially layer-polarization in order to minimize the interaction energy. In the CAF phase, we have that:

[TABLE]

where $\ket{s_{a,b}}$ are spin states with polarization given by the vectors $\mathbf{s}_{a,b}=[\pm\sin\theta_{S}\cos\phi_{S},\pm\sin\theta_{S}\sin\phi_{S},\cos\theta_{S}]$, with tilting angle $\cos\theta_{S}=-\epsilon_{Z}/2u_{\perp}$, and $\mathbf{s}_{\parallel}=[\cos\phi_{S},\sin\phi_{S},0]$. We choose the phases of these states in such a way that $\ket{s_{b}}=\sigma_{z}\ket{s_{a}}$ and $\ket{-s_{b}}=-\sigma_{z}\ket{-s_{a}}$ so the CAF phase continuously matches the F phase of Eq. (39) for $\theta_{S}=0$. As the azimuth $\phi_{S}$ is a free parameter, the CAF phase exhibits a $U(1)$ symmetry. When this solutions exists, it always has lower energy than the regular F phase, so the condition for the presence of the CAF phase is just

[TABLE]

with $\epsilon_{Zc}$ the critical Zeeman field.

III.4 Partially layer-polarized phase

In analogy with the relation between the FLP and F phases, the PLP is similar to the CAF phase but in valley space:

[TABLE]

with $\ket{\mathbf{n}}$ a state with valley-polarization given by $\mathbf{n}=[\sin\theta_{V}\cos\phi_{V},\sin\theta_{V}\sin\phi_{V},\cos\theta_{V}]$, $\cos\theta_{V}=\epsilon_{V}/(u_{z}-u_{\perp})$. The PLP phase also presents an $U(1)$ symmetry. In analogy with the F-CAF phase transition, whenever this solution exists, it has lower energy than the FLP phase so the existence condition for the PLP phase is

[TABLE]

with $\epsilon_{Vc}$ the critical voltage energy. On the other hand, the boundary between the PLP and the CAF phases is placed at

[TABLE]

III.5 Transport gap

An experimental magnitude of interest that can be obtained within the present mean-field computation is the transport gap $\Delta_{HF}$, defined as the energy difference between the lowest energy empty state and the last filled state. From the above results, and within the convention chosen here for the occupied and empty levels, it is immediate to show that in all phases

[TABLE]

in good agreement with Refs. Gorbar et al. (2011); Lambert and Côté (2013), where the only asymmetric interactions considered are interlayer Coulomb interactions.

IV Collective modes

IV.1 Preliminary remarks

We proceed to compute the neutral collective modes of the previous mean-field phases using the TDHFA. The general formalism of the TDHFA and its application to integer QH states is explained thoroughly in Appendix B while Appendix C provides all the technical details on the results presented in this section.

Due to the particular form of the $\nu=0$ QH state, in which the unoccupied levels have different valley-spin polarization with respect to the occupied levels, the excitations correspond to valley-spin waves. Within the TDHFA, the collective modes are obtained from the eigenvalues and eigenvectors of the matrix $\tilde{X}(\mathbf{k})$, given by Eq. (249). The wave vector $\mathbf{k}$ corresponds to the momentum of the so-called magnetoexciton Kallin and Halperin (1984), whose wave function is created by the action of the operator $\hat{M}^{\dagger}_{n\alpha n^{\prime}\alpha^{\prime}}(\mathbf{k})\equiv\hat{M}^{\dagger}_{n\lambda\xi n^{\prime}\lambda^{\prime}\xi^{\prime}}(\mathbf{k})$ on the mean-field ground state $\ket{\Psi_{0}}$, with

[TABLE]

where we have made explicit the dependence in valley $\lambda=K,K^{\prime}$ and spin $\xi=\pm$ indices of the total valley-spin polarization index $\alpha$. As well known from the general theory of integer QH effect, the collective modes are expressed in terms of linear combinations of magnetoexcitons. The resulting dispersion relation is isotropic and only depends on $k=|\mathbf{k}|$, as shown in Appendix C.

In order to classify these modes, we study their symmetries in the whole $8$-dimensional space $01\otimes KK^{\prime}\otimes s$, formed by the tensorial product of the magnetic levels of the ZLL and the valley-spin space. For that purpose, following the notation of Ref. Lambert and Côté (2013), we define the operators:

[TABLE]

which correspond to the components of the spin, valley and orbital pseudospin, respectively. In the above expression, $\mu_{i}$ are the corresponding Pauli matrices in the magnetic index, with $\mu_{z}=\pm 1$ for $n=1,0$.

We consider the behavior of the magnetoexcitons under transformations generated by the previous operators. For instance, for the spin operator, the commutator $[\hat{S}_{i},\hat{M}^{\dagger}_{n\lambda\xi n^{\prime}\lambda^{\prime}\xi^{\prime}}(\mathbf{k})]$ reads

[TABLE]

It is easy to prove that the matrices $G_{i}$ form a representation with spin $1/2\otimes 1/2=1\oplus 0$ of the Lie algebra of $SU(2)$. Thus, spin singlet and triplet magnetoexcitons can be constructed according to the value of the total spin $S$ and its $z$-component $S_{z}$, $\hat{M}^{\dagger}_{n\lambda n^{\prime}\lambda^{\prime},SS_{z}}(\mathbf{k})$. Specifically, the spin-triplet magnetoexcitons are given by

[TABLE]

while the spin-singlet magnetoexciton reads

[TABLE]

Due to the formal analogy with the operators $\hat{L}_{i},\hat{O}_{i}$, we can also build similar singlet and triplet magnetoexcitons in valley and orbital pseudospin.

The importance of the previous considerations arises from the fact that the effective Hamiltonian (29) commutes with the operators $\hat{\mathbf{S}}^{2},\hat{\mathbf{L}}^{2},\hat{S}_{z},\hat{L}_{z}$, with

[TABLE]

so $S,L,S_{z},L_{z}$ are expected to be good quantum numbers that characterize the magnetoexcitons whenever $\ket{\Psi_{0}}$ is also an eigenstate of any of these operators. Thus, the collective modes can be labeled by a general index $\mu$, which represents a set of conserved quantum numbers in valley-spin space. Specifically, the F and FLP states have well-defined quantum numbers for $S,L,S_{z},L_{z}$ which, as discussed in the paragraph following Eq. (188) and in Appendix C, implies that they cannot exhibit dynamical instabilities (i.e., exploding modes with complex frequency and positive imaginary part), at least within the projected model considered in this work. On the other hand, the CAF and PLP states do not satisfy this property and, as shown in the next section, they are indeed able to display dynamical instabilities in some regions of the parameter space.

At $k=0$, apart from these valley-spin symmetries, the $z$ component of the orbital pseudospin, $O_{z}$, is also conserved by the Hamiltonian as it represents the angular momentum of the magnetoexciton Iyengar et al. (2007); Tőke and Fal’ko (2011); Sári and Tőke (2013). In fact, the total orbital pseudospin $O$ is another good quantum number in this limit. As a consequence, we can provide a complete classification of the modes at $k=0$ using their total orbital-valley-spin symmetry, denoting their frequencies as $\omega^{\mu}_{OO_{z}}$. As shown in Appendix C, the energy of the orbital-singlet modes $OO_{z}=00$ is completely independent of the Coulomb interaction strength. On the other hand, the orbital-triplet modes ($O=1$) depend on short-range interactions only through the valley-spin contributions of mean-field energies, but not through many-body corrections. Since Coulomb interactions are dominant, this structure implies that the triplet modes are shifted above by relatively large Coulomb energies and hence, the orbital-singlet modes are the lowest energy modes at $k=0$. In particular, the hierarchy $\omega^{\mu}_{00}<\omega^{\mu}_{11}=\omega^{\mu}_{1-1}<\omega^{\mu}_{10}$ is satisfied. It is worth noting that the triplet modes with $O_{z}=\pm 1$ are degenerate due to the analog Lamb shift.

For $k>0$, however, the orbital singlet and triplet magnetoexcitons are mixed and the previous classification is no longer valid. Then, in order to classify the orbital structure of the different magnetoexcitons for fixed $\mu$, we introduce the discrete index $N=0,1,2,3$ and denote the corresponding frequency as $\omega^{\mu}_{N}(k)$ so they respectively match the collective modes with orbital pseudospin $OO_{z}=00,11,1-1,10$ at $k=0$; hence, $\omega^{\mu}_{0}(k=0)=\omega^{\mu}_{00}$ and so on. The notation is chosen in such a way that, for fixed $\mu$, $\omega^{\mu}_{N}(k)<\omega^{\mu}_{N^{\prime}}(k)$ for $N<N^{\prime}$.

As the collective modes correspond to valley-spin waves, the behavior of the modes for low momentum is $\omega^{\mu}_{N}(k)\simeq\omega^{\mu}_{N}(0)+\rho^{\mu}_{N}k^{2}$, with $\rho^{\mu}_{N}$ the generalized valley-spin stiffness Kallin and Halperin (1984). However, this behavior is modified as $\omega^{\mu}_{0}(k)\simeq v_{G}k$ for the Goldstone modes of the phases with spontaneously broken symmetries.

At large momentum, $kl_{B}\gtrsim 1$, only Coulomb interactions are relevant since $C(\mathbf{k})$ decays as $\sim k^{-1}$ while $R(\mathbf{k})$ does it as $\sim e^{-\frac{(kl_{B})^{2}}{2}}$, where $C(\mathbf{k})$, $R(\mathbf{k})$ are the matrices that arise due to the many-body contributions from Coulomb and short-range interactions, respectively (check Appendix C for their precise definition). For $kl_{B}\gg 1$, the collective modes frequency approach asymptotically the mean-field particle-hole excitation energy $\hbar\omega\simeq\epsilon_{n,\alpha}-\epsilon_{n^{\prime},\alpha^{\prime}}$.

The above considerations imply that the modes $N=1,2,3$ are almost unaffected by short-range interactions since in both limits $kl_{B}\ll 1$, $kl_{B}\gg 1$ they solely depend on them through the valley-spin part of the mean-field contributions, which only provide a constant shift, and for $kl_{B}\sim 1$ Coulomb contributions are the dominant. Indeed, it is shown in Appendix C.2 that the stiffness coefficients of the $N=1$ modes depend very weakly on $u_{\perp},u_{z}$ while those of the $N=3$ modes are completely independent of them. Moreover, due to the analog Lamb shift, the $N=2$ modes can be computed explicitly as they correspond to a symmetric combination of the orbital $OO_{z}=11,1-1$ modes as given by Eqs. (264), (265), and it is seen that their dependence on $k$ only involves Coulomb interactions. As a result, we expect the dispersion relation of the $N=1,2,3$ modes to be quite independent of their valley-spin structure and then recover essentially the same results of Refs. Tőke and Fal’ko (2011); Lambert and Côté (2013), where only (intra- and inter-layer) Coulomb interactions are taken into account.

The situation is quite different for the $N=0$ modes, since at $k=0$ they correspond to the orbital-singlet modes which are independent of the Coulomb strength and the only modes that experience many-body corrections coming from short-range interactions. Hence, their behavior at low momentum greatly depends on the valley-spin structure of the modes. Besides, as they present a positive stiffness (see next section), $\omega^{\mu}_{0}(k)$ increases monotonically with $k$. Hence, the orbital-singlet modes at $k=0$ are a) the lowest energy excitations of the system and b) the most sensible modes to the valley-spin structure of interactions. We therefore conclude that the orbital-singlet modes are the most natural candidates to characterize the phase transitions. We note that the orbital-singlet modes are denoted as the “even” modes in Ref. Tőke and Fal’ko (2011).

IV.2 Collective modes: analysis of the results

In this section, we compute the dispersion relation for the collective modes of the different ground states obtained in Sec. III, see Eqs. (39),(40),(III.3),(44). For each phase, we discuss the symmetries of the excitations and of the ground state, give the explicit expression for the frequencies of the orbital-singlet modes $\omega^{\mu}_{00}$ in order to study the stability of the corresponding phase and the values of the associated stiffness, $\rho^{\mu}_{0}$, and plot the dispersion relation of the different collective modes, $\omega^{\mu}_{N}(k)$. We refer the reader to Appendix C.2 for the specific details on the calculation of the collective modes.

IV.2.1 Ferromagnetic phase

The ferromagnetic state has valley-pseudospin $L,L_{z}=0$ and spin $S,S_{z}=2N_{B}$, with $4N_{B}$ the total number of electrons in the ZLL. All magnetoexcitons carry spin $S=1,S_{z}=-1$ so we use the valley pseudospin in order to characterize the collective modes, $\omega^{\mu=LL_{z}}_{N}(k)$. The excitations of the ferromagnetic state consists of spin-flip excitations and full-flip excitations. The spin-flip excitations involve transitions between two different spin polarizations, keeping the same valley polarization, while in full-flip excitations both valley and spin indices are flipped at the same time. Specifically, the spin-flip excitations are characterized by the modes with $L=1,0$ and $L_{z}=0$. The energy of the orbital singlet associated to these modes, $\hbar\omega^{L0}_{00}$, is:

[TABLE]

The spin-singlet mode $\hbar\omega^{00}_{00}=2\epsilon_{Z}>0$ is the Larmor mode Murthy et al. (2016), while the spin-triplet mode characterizes the F-CAF phase transition, $\hbar\omega^{10}_{00}=2\epsilon_{Z}+4u_{\perp}=2(\epsilon_{Z}-\epsilon_{Zc})$; precisely, $\omega^{10}_{00}=0$ is the boundary [see Eq. (43)] between the two phases. The appearance of a such a gapless mode corresponds to the spontaneously broken $U(1)$ symmetry of the CAF phase. In the region where the CAF phase is present, the system is energetically unstable as $\omega^{10}_{00}<0$.

Regarding full-flip excitations, they are characterized by the triplet modes $L_{z}=\pm 1$. Indeed, their dispersion relation is the same, solely shifted by the layer voltage, $\omega^{1-1}_{N}(k)=\omega^{11}_{N}(k)+4\epsilon_{V}$, so we only need to compute $\omega^{11}_{N}(k)$ to characterize them. In particular, for the orbital-singlet mode,

[TABLE]

which is proportional to the difference between the mean-field energies of the FLP and F states. Hence, in analogy to the case of spin-flip excitations, $\omega^{11}_{00}=0$ is the boundary between the two phases.

The existence of such a gapless mode can be understood from the appearance of a mean-field symmetry right at the boundary between the FLP and the F phases. We define a mean-field symmetry as that which is not of the complete Hamiltonian but only of the mean-field HF solutions [see discussion after Eq. (203) for more details]. Specifically, when condition (41) is satisfied, the HF equations present a continuously degenerate mean-field ground state described by the parameters $t,\phi$ and characterized in valley-spin as

[TABLE]

with energy $E[P(t,\phi)]=E(P_{F})\cos^{2}t+E(P_{FLP})\sin^{2}t$, $P_{F,FLP}$ being the projectors of the F, FLP phases, respectively [see Eqs. (III),(39) and (40)], and the matrices $\Pi^{i}_{j}$ defined as $\Pi^{i}_{j}\equiv\frac{1}{2}\tau_{i}\otimes\sigma_{j}$. Thus, right at the phase boundary, the total mean-field energy of the state does not depend on the parameters $t,\phi$ and a continuous mean-field symmetry arises from this fact. In particular, the parameter $t$ describes the phase transition, with $t=0$ corresponding to the F state while the FLP state is obtained for $t=\frac{\pi}{2}$. The associated total state, denoted as $\ket{\Psi(t,\phi)}$, can be connected to the F state $\ket{F}$ by a continuous unitary transformation of the form:

[TABLE]

and, in terms of the projectors $P$, $P(t,\phi)=e^{G(t,\phi)}P_{F}e^{-G(t,\phi)}$, with $G(t,\phi)$ the matrix version of the operator $\hat{G}(t,\phi)$. Note that the azimuth $\phi$ simply arises from adding a trivial rotation in valley and/or spin space to $\hat{G}(t,0)$. The generator of Eq. (IV.2.1) can be further rewritten as

[TABLE]

The operators $\hat{\Pi}^{x}_{i},\hat{\Pi}^{y}_{i}$, along with $\hat{S}_{i},\hat{L}_{z}$, form a representation of $SO(5)$ that is an exact symmetry of the total Hamiltonian for $u_{\perp}+u_{z}=0$ and $\epsilon_{Z}=\epsilon_{V}=0$ Wu et al. (2014). Thus, the exotic $SO(5)$ symmetry existing at the F-FLP boundary for zero Zeeman and layer voltage terms survives as a weaker mean-field version at the same phase boundary in the realistic scenario $u_{z}+u_{\perp}\neq 0$. Remarkably, this gapless mode behaves as $\omega^{11}_{0}(k)\sim k^{2}$ for low momentum in contrast to the linear dependence found for the Goldstone modes of the phases with spontaneously broken symmetries.

The stiffness for the orbital-singlet modes is given, for dominant Coulomb interactions $F_{C}\gg u_{\perp},u_{z}$ [see Eq. (C.2.1) for the exact expression], by

[TABLE]

with $u^{\mu}$ a short-range coupling that depends on the valley-spin symmetry of the mode. Specifically, $u^{00}=-u_{z}-2u_{\perp}$, $u^{10}=-u_{z}+2u_{\perp}$ and $u^{11}=u^{1-1}=u_{z}$.

We plot in the upper left corner of Fig. 2 the dispersion relation of the full-flip (left panel) and spin-flip (right panel) excitations. In the case of spin-flip excitations, the $L=0$, $(L=1)$ branches are plotted in solid (dashed) line. As the corresponding mean-field energies are degenerate and the modes $N=1,2,3$ depend very weakly on the short-range interactions, the curves of both branches are very close to each other; indeed, this degeneracy becomes exact for the $N=2$ modes as explained in the previous section. The characteristic negative stiffness of the $N=1$ modes can also be observed in the plots [see Eq. (273)]. Only the frequencies $\omega^{00}_{0}(k)$, $\omega^{10}_{0}(k)$ are clearly distinguished due to the many-body corrections arising from short-range interactions at low momentum.

When comparing the plots of spin-flip and full-flip excitations, we see that the curves for the $N=1,2,3$ are extremely similar. Again, this can be explained by invoking the dominant character of Coulomb interactions and the weak dependence on short-range effects. As expected, the main differences are observed once more for the lowest energy modes $N=0$. Variations of the in-plane magnetic field or the interlayer voltage only affect through a linear energy offset in $\epsilon_{Z},\epsilon_{V}$ that depends on the valley-spin structure of each mode, as given by Eqs. (56), (57). Finally, we remark that variations of the different parameters only change quantitatively but not qualitatively the plots.

IV.2.2 Full layer-polarized phase

As this phase is the analog of the F phase in valley space, the results for the FLP phase are similar but interchanging $L,L_{z}$ and $S,S_{z}$. Specifically, the FLP state has spin $S,S_{z}=0$ and valley-pseudospin $L,L_{z}=2N_{B}$. All magnetoexcitons have valley pseudospin $L=1,L_{z}=-1$, so we use their spin to classify them, $\omega^{\mu=SS_{z}}_{N}(k)$. Since the interaction is spin-independent, the triplet modes present the same dispersion relation, only shifted by the Zeeman energy, $\omega^{1\pm 1}_{N}(k)=\omega^{10}_{N}(k)\mp 2\epsilon_{Z}$, hence we only need to compute $\omega^{11}_{N}(k)$ to characterize them.

The excitations of the FLP state correspond to valley-flip and full-flip excitations. Valley-flip excitations are characterized by the modes with $S=1,0$ and $S_{z}=0$. In particular, the energy of the orbital-spin-singlet mode is

[TABLE]

In analogy to the case of the F-CAF transition [see Eq. (56) and related discussion], $\omega^{00}_{00}=0$ is the boundary between the FLP and PLP phases; the gapless character of this mode arises from the broken $U(1)$ symmetry of the PLP phase.

Regarding the F-FLP boundary, we find that it is characterized by the mode

[TABLE]

which is just $-\hbar\omega^{S=1S_{z}=1}_{00}$ for the F state, see Eq. (57). Reasoning as before, the appearance of this gapless mode at the boundary between the two phases results from the residual mean-field symmetry shown in Eq. (58).

It is worth noting the subtle difference between the F-CAF and FLP-PLP transitions: while the former is governed by a valley triplet mode, the latter is governed by a spin singlet mode. This is due to the fact that, in the CAF phase, the occupied states $a,b$ correspond to two different spins $s_{a},s_{b}$ for the two valley polarizations $\pm n_{z}$, one being the other rotated $\pi$ in the $x-y$ spin plane. In contrast, in the PLP phase, all the occupied states present the same valley polarization, specified by the vector $\mathbf{n}$, independently of their spin. This relation can be observed in the $P$ matrices: the F-CAF transition is characterized by the appearance of an operator that goes as $\sim\sigma_{i}\tau_{z}$, which has valley pseudospin $L=1,L_{z}=0$, while the FLP-PLP transition only varies the orientation of the projector in valley space, which has zero spin $S,S_{z}=0$. On the other hand, the F-FLP transition is governed by triplet modes in both spin and valley spaces. This fact can be understood from the expression of the generator $\hat{G}(t,\phi)$ of the mean-field symmetry, Eq. (60), as it has $L=1,S=1$.

The stiffness coefficients of the orbital-singlet modes are also given by Eq. (61) but now $u^{00}=u_{z}+4u_{\perp}$ and $u^{11}=u_{z}$, noting that $\rho^{1-1}_{N}=\rho^{10}_{N}=\rho^{11}_{N}$.

In the upper right corner of Fig. 2, we represent the dispersion relation of the spin triplet (left panel) and singlet (right panel) excitations. The qualitative trends and properties of the curves are similar to those of the F state.

IV.2.3 Canted anti-ferromagnetic phase

The CAF state has only well defined value of the $z$ valley pseudospin, $L_{z}$, so we use this quantum number to characterize the magnetoexcitons. Due to this lack of symmetry, dynamical instabilities can appear as shown in Appendix C. The sector $L_{z}=0$ is characterized by two branches, labeled as $\omega_{N}^{\pm}(k)$, that correspond to generalizations of the spin-flip modes of the F phase; they match the $L=0,1$ modes for $\theta_{S}=0$, respectively. The energy for their orbital-singlet modes is:

[TABLE]

In analogy to the F state, the previous equation gives the Larmor mode $\hbar\omega^{+}_{00}=|4u_{\perp}|\cos\theta_{S}=2\epsilon_{Z}>0$ and the Goldstone mode $\hbar\omega^{-}_{00}=0$. The appearance of this Goldstone mode is expected because the CAF phase spontaneously breaks the $U(1)$ symmetry of the Hamiltonian (see discussion in Appendix B about Goldstone modes in the TDHFA for more details). The velocity of the Goldstone mode, given by $\omega^{-}_{0}(k)\simeq v_{G}k$ in the limit $k\rightarrow 0$, can be computed exactly and is given by Eq. (281). In particular, in the regime $F_{C}\gg u_{\perp},u_{z}$, it is well approximated by

[TABLE]

with $A=\epsilon_{Zc}\sin^{2}\theta_{S}$ and $\Delta=u_{z}-2u_{\perp}$. Near the phase transition, the velocity behaves as $v_{G}\sim\sqrt{A}\sim\sqrt{\epsilon_{Zc}-\epsilon_{Z}}$, in agreement with the critical behavior described in Ref. Murthy et al. (2016). Thus, the F-CAF transition is characterized in the CAF side by a vanishing group velocity instead of a closing energy gap. There are no complex-frequency modes in the $L_{z}=0$ sector.

The remaining collective modes correspond to the sectors with $L_{z}=\pm 1$, with frequencies $\omega^{1\pm 1}_{N}(k)$. In analogy to the F phase, they satisfy $\omega^{1-1}_{N}(k)=\omega^{11}_{N}(k)+4\epsilon_{V}$ so we just compute $\omega^{11}_{N}(k)$. We find that

[TABLE]

Interestingly, there are situations in which the quantity in the square root is negative and then a dynamical instability appears. It is straightforward to show that $\omega^{11}_{00}$ is purely real whenever

[TABLE]

is fulfilled. On the other hand, $\omega^{11}_{00}=0$ at the boundary between the PLP and CAF phases, i.e, whenever Eq. (46) is satisfied. In analogy to the F-FLP transition, this gapless mode can be understood in terms of a mean-field symmetry arising right at the phase boundary as the state $\ket{\Psi(t,\phi_{S},\phi_{V})}$, characterized by the projector

[TABLE]

with

[TABLE]

has an energy $E[P(t,\phi_{S},\phi_{V})]=E(P_{CAF})\cos^{2}t+E(P_{PLP})\sin^{2}t$, $P_{CAF,PLP}$ being the projector of a CAF, PLP state characterized by the angles $\theta_{S,V}$ and $\phi_{S,V}$. Therefore, as in the F-FLP case, the complete manifold of states is degenerate at the CAF-PLP boundary, giving rise to a continuous mean-field symmetry, with the parameter $t$ describing again the phase transition: the value $t=0$ yields the CAF phase while $t=\frac{\pi}{2}$ gives the corresponding PLP state at the other side of the phase transition. For simplicity, we do not give the spinors $\chi_{\alpha}(t,\phi_{S},\phi_{V})$ as their expression is quite cumbersome and unimportant for the present discussion.

A summary of the stability conditions for the CAF phase is given by

[TABLE]

Note that $\omega^{11}_{00}>0$ as long as we stay in the region where the ground state is the CAF phase.

With respect to the stiffness coefficients, we find in the regime $F_{C}\gg u_{\perp},u_{z}$ that [see Eq. (C.2.3) for the exact expression]

[TABLE]

with $u^{+}=-u_{z}-2u_{\perp}\cos^{2}\theta_{S}$ and $u^{11}=u_{z}+2u_{\perp}\sin^{2}\theta_{S}$. We note that the Goldstone mode has no stiffness coefficient due to its linear behavior and that $\rho^{1-1}_{N}=\rho^{11}_{N}$. Remarkably, the denominator of the above equation takes the value of the Larmor frequency for the $+$ mode and the resulting stiffness $\rho^{+}_{0}$ becomes relatively large due to the small value of the Zeeman energy.

In the lower left corner of Fig. 2 we plot the dispersion relation of the modes with $L_{z}=1$ (left panel) and $L_{z}=0$ (right panel). We see that a mode with zero energy that grows linearly with $k$ appears in the sector with $L_{z}=0$, corresponding to the Goldstone mode of the CAF state. Apart from this consideration, the qualitative form of the curves is similar to the previous cases. We note that, while a variation of the interlayer voltage only provides a trivial energy shift for the $L_{z}=\pm 1$ modes, the role of the in-plane magnetic field is much more critical as it determines the background mean-field state by tuning the canting angle $\theta_{S}$ through the Zeeman energy $\epsilon_{Z}$.

In left Fig. 3 we consider an experimentally unrealistic case with different values of the coupling constants, $-u_{\perp}>u_{z}>0$, in order to study the trends of the appearance of a complex-frequency mode in the CAF phase. The plot shows the real (solid curve) and the imaginary (dashed curve) parts of the unstable mode. The decrease (and eventual vanishing) of the imaginary part of the frequency is a consequence of the exponential decay of the many-body contributions from the short-range interactions, also giving rise to the term responsible for the appearance of the instability [see Appendix C for the details]. It is straightforward to show that the real part of the dynamical instability satisfies $\hbar\text{Re}[\omega^{11}_{0}(k)]=-2\epsilon_{V}$, as can be observed in the plot.

IV.2.4 Partially layer-polarized phase

In contrast with the CAF state, the PLP state has well defined quantum numbers in spin and valley-pseudospin. Indeed, the valley spin structure of the PLP phases is analog to that of the FLP phase but replacing $L_{z}$ by the component of the valley pseudospin along the direction of the vector $\mathbf{n}$, $L_{n}$. The excitations are also characterized by their spin, $\omega^{SS_{z}}_{N}(k)$, and the triplet modes satisfy $\omega^{1\pm 1}_{N}(k)=\omega^{10}_{N}(k)\mp 2\epsilon_{Z}$. However, the Hamiltonian does not commute with the operator $\hat{L}_{n}$ and because of this, dynamical instabilities can appear.

For the spin singlet, the orbital-singlet mode is gapless,

[TABLE]

as it is the Goldstone mode associated to the spontaneously broken $U(1)$ symmetry. Its velocity is computed in similar terms to that of the CAF phase, Eq. (65), but now $A=\epsilon_{Vc}\sin^{2}\theta_{V}$ and $\Delta=-u_{z}-4u_{\perp}$, presenting a similar critical behavior near the phase transition, $v_{G}\sim\sqrt{A}\sim\sqrt{\epsilon_{Vc}-\epsilon_{V}}$.

Considering the spin-triplet modes, the orbital-singlet frequency is

[TABLE]

It can be proven that $\omega^{11}_{00}$ is real if:

[TABLE]

A zero-frequency mode $\hbar\omega^{11}_{00}=0$ appears at the boundary between the PLP and CAF phases, arising from the mean-field symmetry already discussed for the CAF phase [see Eq. (IV.2.3) and related]. Hence, in analogy to the F-FLP case, the CAF-PLP transition is also governed by triplet modes. Putting all together as for the CAF state:

[TABLE]

In contrast to the CAF case, for the appearance of dynamical instabilities in the PLP state it is only required a sufficiently low value of the voltage; nevertheless, this instability occurs in a region of parameter space where the CAF phase is the ground state.

For the stiffness coefficients of the spin-triplet modes, Eq. (71) is still valid for sufficiently dominant Coulomb interactions, using the value $u^{11}=u_{z}-A$.

In the lower right corner of Fig. 2 we plot the dispersion relation of the spin triplet (left panel) and singlet (right panel) excitations. The trends are similar to the previous cases. In particular, in analogy with the CAF case, variations in the in-plane magnetic field only vary the trivial energy shift for the $S_{z}=\pm 1$ modes while the influence of the interlayer voltage is much more important since it controls the value of $\theta_{V}$.

In right Fig. 3 we study a case with sufficiently low value of the voltage so that the PLP state is dynamically unstable and hence, the triplet modes have complex frequency. The qualitative behavior resembles that of the unstable mode appearing for a CAF state, satisfying a similar relation for the real part, $\hbar\text{Re}[\omega^{11}_{0}(k)]=-2\epsilon_{Z}$.

Finally, in Fig. 4, we represent the velocity of the Goldstone modes for both the CAF and PLP phases as a function of the critical parameter $\delta_{c}$, defined as $\delta_{c}\equiv(\epsilon_{Zc,Vc}-\epsilon_{Z,V})/\epsilon_{Zc,Vc}$ for the CAF and PLP phases, respectively. For low in-plane magnetic field or interlayer voltage, the velocities saturate the limit $\theta_{S,V}\rightarrow\frac{\pi}{2}$. The agreement of the approximation (65) for dominant Coulomb interactions with the exact result of Eq. (281) is excellent (not shown). The corresponding plot for the energy gaps of the F and FLP phases is trivial as they are proportional to $-\delta_{c}$, see Eqs. (56), (62).

V Monolayer graphene

All of the previous work can be straightforwardly adapted to the $\nu=0$ QH state of monolayer graphene. We briefly revisit the content of Secs. II-IV, adjusting the results to the present case.

V.1 Effective Hamiltonian

We start by writing the effective Hamiltonian of monolayer graphene (see Ref. Kharitonov (2012a) for a more complete review). The corresponding field operator is given by:

[TABLE]

where now $A,B$ are the two sublattices of the same graphene layer. Once more, the two sublattices are interchanged in the $K^{\prime}$ valley so we denote the corresponding subspace as $\bar{A}\bar{B}$. The Hamiltonian is decomposed in the same fashion as in the bilayer case, $\hat{H}=\hat{H}_{0}+\hat{H}_{C}+\hat{H}_{sr}$, but the single-particle Hamiltonian $\hat{H}_{0}$ now reads

[TABLE]

(note that there is no layer voltage term). After an appropriated phase transformation, the matrix $H_{B}$ can be written as:

[TABLE]

with $v_{F}\simeq 10^{6}\text{m/s}$ the Fermi velocity. The corresponding eigenfunctions and eigenvalues are similar to those of the bilayer, $\Psi^{0}_{n,k,\alpha}(\mathbf{x})=\Psi^{0}_{n,k}(\mathbf{x})\chi_{\alpha}$, with

[TABLE]

for $n\neq 0$ and

[TABLE]

for the ZLL, that now corresponds to just the magnetic level $n=0$.

Due to the different scaling of the cyclotron frequency with the magnetic field, the Zeeman energy now satisfies

[TABLE]

while an order of magnitude analysis of the Coulomb interaction gives the dimensionless strength

[TABLE]

Note that it does not depend on the value of the magnetic field, in contrast to the case of bilayer graphene. The order of magnitude of short-range interactions is similar to the bilayer case and their typical energy also goes as $\sim F_{C}d/\hbar\omega_{B}l_{B}\ll 1$. Once more, for $\kappa=1$, Coulomb interactions are not weak. However, by taking $\kappa\sim 5$ we can get $F_{C}\sim 0.4\hbar\omega_{B}$, which can be regarded as a small value.

Under this assumption, we neglect LL mixing and project the Hamiltonian onto the ZLL. The states for the ZLL of monolayer graphene are also restricted to the $KK^{\prime}\otimes\bar{B}\otimes s$ subspace, which means that they are localized, for each valley, in one sublattice or the other and hence the sublattice degree of freedom becomes equivalent to the valley degree of freedom. Reasoning in the same fashion as in the bilayer case, the resulting effective Hamiltonian is formally equal to that of Eq. (29) but with $\epsilon_{V}=0$:

[TABLE]

The Dirac sea that creates the mean-field potential $V_{DS}(\mathbf{x},\mathbf{x^{\prime}})$ is composed now by all the occupied states with $n\leq-1$.

V.2 Hartree-Fock equations and mean-field phase diagram

The HF equations have the same form of Eq. (III). As the $\nu=0$ QH state corresponds to half-filling of the ZLL, the electrons occupy in the same way two orthogonal spinors $\chi_{a,b}$ in valley-spin space and leave empty the remaining orthogonal spinors $\chi_{c,d}$. Then, after projecting the HF equations into the orbital part of the wave functions, we get the algebraic equation

[TABLE]

where now $u_{i}=g_{i}/2\pi l^{2}_{B}$. The term $F_{00}/2$ arises from the interaction with the inert Dirac sea; however, in monolayer graphene, it is just a trivial energy shift as there is only one magnetic level in the ZLL.

The previous equation presents the same valley-spin structure of Eq. (33). As a consequence, the spinorial solutions $\chi_{\alpha}$ to the HF equations are identical to those of bilayer graphene and their mean-field energies are:

[TABLE]

Moreover, the total energy of the ground state (per wave vector state) is:

[TABLE]

with $E(P)$ given by Eq. (III). Therefore, we conclude that the mean-field phase diagram for the $\nu=0$ QH state in monolayer graphene is that of the bilayer for $\epsilon_{V}=0$ Kharitonov (2012a, b). In that case, the PLP phase of the bilayer changes to a fully interlayer coherent phase (ILC) since $\theta_{V}=\pi/2$ and the vector $\mathbf{n}$, pointing the polarization in the valley space, is fully contained in the $x-y$ plane, $\mathbf{n}=[\cos\phi_{V},\sin\phi_{V},0]$. For monolayer graphene, the equivalent of the ILC phase is the Kekulé distortion (KD) phase: since the valley degree of freedom is equivalent to the sublattice in the zero Landau level of monolayer graphene, this state corresponds to a coherent mixture of the two sublattices. For the same reason, the FLP phase is now a charge-density wave (CDW) phase. There are experimental evidences that the phase for the $\nu=0$ QH state of monolayer graphene for zero in-plane component of the magnetic field is also the CAF state Young et al. (2014), so $u_{z}>-u_{\perp}>0$ as in the bilayer case. However, as an analog of the layer voltage is lacking, in principle only the transition between CAF and F phases can be explored by changing the in-plane component of the magnetic field.

For the mean-field transport gap, we find that it satisfies a similar relation to Eq. (47),

[TABLE]

V.3 Collective modes

The TDHFA developed for computing the collective modes within the ZLL of the $\nu=0$ QH state of bilayer graphene can be straightforwardly translated to the monolayer. Indeed, the valley-spin structure of the excitations is the same as in the bilayer case, while the orbital structure is trivial as it corresponds to a one-dimensional subspace spanned by the magnetic level $n=0$. Hence, excitations are characterized just by their symmetry in valley-spin space, $\omega^{\mu}(k)$, with $\mu$ labeling the same set of conserved quantum numbers described in Sec. IV.2. Due to the simple structure in orbital space, the collective modes for monolayer graphene can be computed analytically and their explicit expression is given at the end of Appendix C. In particular, at $k=0$, $\omega^{\mu}(0)=\omega^{\mu}_{00}$ being $\omega^{\mu}_{00}$ the frequency of the corresponding orbital-singlet mode in the bilayer case. As a consequence, phase transitions are characterized in the same way.

We now plot the dispersion relation for every phase. The results are shown in Fig. 5. Note that for the study of the CDW and KD phases we have to use unrealistic values for the coupling constants. The qualitative trends of the several dispersion relations are similar to those of the $N=0$ orbital modes in bilayer graphene; in fact, at $k=0$, they have the same formal expression as explained above. When they exist, the dispersion relation of the dynamical instabilities of the CAF and KD phases is qualitatively similar to that of the bilayer case (not shown).

We note that the neutral excitations of the F and CDW phases were computed in Ref. Doretto and Smith (2007) using a bosonization approach and a different type of short-range interactions, nevertheless finding similar analytical expressions for the frequency of the modes. In addition, Ref. Murthy et al. (2016) obtained the dispersion relation for the F and CAF phases using the same valley-asymmetric interactions here considered but replacing the long-range Coulomb interaction by an effective short-range one in order to simulate the valley-spin stiffness of the waves.

VI Effects of Landau-level mixing

In this section, we consider the usual case where $F_{C}\gtrsim\hbar\omega_{B}$ and LL mixing cannot be neglected. One possible way to take it into account is through a static screening of the Coulomb interaction in the large-$N$ approximation Foster and Aleiner (2008); Basko and Aleiner (2008); Aleiner et al. (2007); Lemonik et al. (2010); Kharitonov (2012a). Other approaches consider dynamical screening within the same approximation Gorbar et al. (2010, 2011, 2012) or allow for LL mixing in the TDHFA formalism Tőke (2013). As we are mainly interested in the low-energy modes describing the phase transitions, static screening is expected to provide a good approximation in this limit; on the other hand, screening by LL mixing is not expected to describe correctly the dispersion relation near $k=0$ Tőke (2013), which is precisely the most interesting region for our purposes.

In the large-$N$ approximation, the effective Coulomb interaction potential is obtained using a RPA-type screening

[TABLE]

where $V_{0}(\mathbf{k})$ is the Fourier transform of the bare Coulomb interaction and $\Pi^{0}(\mathbf{k},\omega)$ is the non-interacting polarization for the $\nu=0$ QH state,

[TABLE]

with $n_{k},n_{l}$ taking values for all integers except $-1$, as in Eq. (23). In the same fashion, $D^{0}_{nk\alpha_{k}n_{l}\alpha_{l}}(\omega)$ is the non-interacting two-particle propagator [see Eq. (225) and related discussion for more details], with the non-interacting values for the energies and occupation numbers. Here, $A^{(2)}_{nn^{\prime}}(\mathbf{k})$ is the bilayer graphene magnetic form factor:

[TABLE]

$A_{nn^{\prime}}(\mathbf{k})$ being the usual magnetic form factor, given by Eq. (A.1). The polarization $\Pi^{0}(\mathbf{k},\omega)$ is the Fourier transform of the non-interacting density-density correlation function and is obtained following an analog calculation to that described in Appendix B.2 but using the bare vertex and allowing for LL mixing.

After neglecting small corrections due to the Zeeman effect and the layer voltage, the free static polarization reads:

[TABLE]

with $N=4$ the number of the valley-spin components of the field, which is expected to be a sufficiently large value to provide a good approximation Basko and Aleiner (2008); Lemonik et al. (2010); Kharitonov (2012a). The function $f(x)$ is dimensionless as $|A^{(2)}_{n_{k}n_{l}}(\mathbf{k})|^{2}$ only depends on the momentum through the quantity $kl_{B}$. The effective screened Coulomb interaction can be then rewritten as:

[TABLE]

where we have made explicit the rotational invariance of all the magnitudes. For $F_{C}\ll\hbar\omega_{B}$, the screened interaction reduces to the bare Coulomb interaction used in the previous calculations.

The dimensionless function $f(x)$ is plotted in Fig. 6, in perfect agreement with the results of Refs. Gorbar et al. (2010, 2012). For small $x$, it behaves as $f(x)\simeq\alpha x^{2}$,

[TABLE]

and then $\bar{V}(k)\simeq V_{0}(k)$ for low momentum. On the other hand, $f(x)\sim 1$ for $x\gtrsim 1$ so for relevant momentum $kl_{B}\sim 1$ the screened potential satisfies:

[TABLE]

The dimensionless strength of the screened Coulomb interaction is then of order $N^{-1}=0.25\ll 1$ and thus we can now safely neglect LL mixing and restrict once more ourselves to the ZLL.

The other effect that results from allowing LL mixing is the renormalization of the coupling constants Kharitonov (2012a, b), that amounts to replace the bare coupling constants by their renormalized values, $\bar{g}_{\perp},\bar{g}_{z}$. However, this process is only strong for low values of the magnetic field Kharitonov (2012b) so we do not study its effect in detail here, treating again the renormalized coupling constants as inputs for the computations.

Taking into account the previous observations, we consider the following effective Hamiltonian for the ZLL

[TABLE]

instead of that of Eq. (29). We see that the result of considering LL mixing is reduced to replacing the Coulomb interaction and the coupling constants by their effective values $\bar{V}$ and $\bar{g}_{i}$. Note that, since the excitations within the ZLL change valley or spin, the direct contribution of the Coulomb interaction vanishes and then there is no double-counting in this approximation (see Appendix B for more details about the diagrammatic formalism of the TDHFA).

This kind of effective renormalized Hamiltonian has already been used for studying the $\nu=0$ state or fractional QH effect in the ZLL Kharitonov (2012a, d); Abanin et al. (2013). As the valley-spin structure of the Hamiltonian is preserved, the solutions to the corresponding HF equations are the same as before and thus the phase diagram is that of Fig. 1 but replacing $u_{i}$ by $\bar{u}_{i}$. The same is true for the HF energies, obtained from the screened values $\bar{F}_{nm}$ [see Eq. (164)], finding numerically that they satisfy the same relation $\bar{F}_{00}>\bar{F}_{10}>\bar{F}_{11}$ and then $\bar{F}_{0}>\bar{F}_{1}$, so the equation for the transport gap is given by the adapted version of Eq. (47).

Regarding the collective modes, using the same reasoning, we find that all the orbital-valley-spin classification of Sec. IV still holds, including the orbital pseudospin structure at $k=0$. Hence, only the quantitative values of the dispersion relation change but not the qualitative features and the resulting physical conclusions. In particular, the structure of the Goldstone modes and of the dynamical instabilities is preserved.

We now plot the numerical results for the computation of the dispersion relation using the screened Coulomb potential of Eq. (107) for a Coulomb interaction strength $F_{C}=4\hbar\omega_{B}$ and keeping the same values for the rest of parameters (including $\bar{u}_{i}=u_{i}$) in order to observe the effects of screening. In Fig. 7, we plot the dispersion relation of the modes for every phase and in Fig. 8 we represent the real and imaginary parts of the frequency of a dynamically unstable mode. We observe that the qualitative structure of the curves is similar to the unscreened case.

All the calculations developed in this section can be straightforwardly adapted to monolayer graphene, obtaining similar conclusions.

VII Experimental remarks

We make here some experimental remarks about the magnitudes computed in this work. First, we analyze the behavior of the transport gap, a quantity that can be measured, for instance, through local compressibility measurements Martin et al. (2010), electronic transport measurements Zhao et al. (2010) or using the bias as a spectroscopic tool Velasco Jr et al. (2012); Velasco et al. (2014); Shi et al. (2016). As revealed by Eq. (47), the transport gap has a contribution from Coulomb interactions, independent of the phase, and another contribution which actually depends on the valley-spin structure of the state. Using Eqs. (107), (137), it is straightforward to show that the screened value of the Coulomb contribution satisfies $\bar{F}_{1}=\hbar\omega_{B}g\left(\frac{F_{C}}{\hbar\omega_{B}}\right)$, with $g(x)$ a dimensionless function that asymptotically behaves as [check Eq. (166) and related discussion for more details]

[TABLE]

Hence, for low $F_{C}/\hbar\omega_{B}$ we recover the result for bare Coulomb interactions $\bar{F}_{1}\simeq F_{1}\propto\sqrt{B_{\perp}}$ while for strongly screened Coulomb interactions, due to the mild behavior of the logarithm, we are close to a linear behavior $\bar{F}_{1}\sim B_{\perp}$, already predicted in Ref. Gorbar et al. (2012). On the other hand, the short-range energies also scale linearly with the perpendicular magnetic field, $u_{\perp},u_{z}\sim B_{\perp}$, at least for moderate fields $B_{\perp}\gtrsim 2~{}\text{T}$ Kharitonov (2012b). Hence, the gap is expected to grow linearly with $B_{\perp}$, as confirmed by experiments Martin et al. (2010); Velasco Jr et al. (2012); only for very high fields and large values of the dielectric constant $\kappa$ the unscreened regime $\Delta_{HF}\sim\sqrt{B_{\perp}}$ can be reached (note that even in this regime the Coulomb contribution is still dominant over those of short-range interactions). In fact, this behavior has been reported for transport gaps of other integer and fractional QH states Zhao et al. (2010); Shi et al. (2016).

Interestingly, more information about the short-range energies can be obtained by comparing the transport gap of two different phases. Suppose that the transport gap for a F state with $\epsilon_{V}=0$ and a relatively high Zeeman energy $\epsilon_{Z,F}$ is measured and the process is repeated, keeping constant $B_{\perp}$, for a FLP state with voltage energy $\epsilon_{V,FLP}$ and negligible Zeeman term $\epsilon_{Z}\simeq 0$. Then, as the Coulomb contribution is expected to be the same, by subtracting both gaps one finds

[TABLE]

For sufficiently high in-plane magnetic field for the F phase, $\epsilon_{Z,F}$ becomes independent of $B_{\perp}$ and hence one can obtain the value of $u_{z}+u_{\perp}$ from a linear fit of $\Delta_{HF,\text{F}}-\Delta_{HF,\text{FLP}}$ versus $B_{\perp}$. We note that, in order to determine the transport gap, this magnetoexciton gap has to be compared with the skyrmionic gap. According to Ref. Tőke and Fal’ko (2011), skyrmions are expected to be energetically favorable only for very large magnetic fields $B_{\perp}\gtrsim 30~{}\text{T}$.

Regarding the detection of the collective modes, it is well known that their frequencies are the poles of the response functions to external fields $\hat{H}_{\rm{ext}}(\mathbf{x},t)$ (see Appendix B for a detailed explanation). We now present a study of the required valley-spin structure of the operator $\hat{H}_{\rm{ext}}$ in order to detect the various modes. For instance, a magnetic field in the $i$ direction couples to the spin density, proportional to the operator $\hat{S}_{i}$ while a layer voltage difference couples to the charge-density difference between the layers, proportional to the $\hat{L}_{z}$ operator. A recent promising work has shown that the layer polarization can be measured in detail in capacitive measurements, arising as a powerful tool to characterize quantum Hall states in bilayer graphene samples Hunt et al. (2016). Unfortunately, the operators $\hat{L}_{x,y}$, corresponding to the effective ladder operators between the layers in the ZLL, are not so easily translated to real physical observables. Compound operators of the form $\hat{S}_{i}\hat{L}_{z}$ can be interpreted as a spin-density difference between the layers although it is not clear how they could be measured. Note that in the projected model considered in this work, the modes do not couple to density perturbations (proportional to a scalar in valley-spin space) since the excitations are valley-spin waves that change necessarily one of these two quantum numbers Murthy et al. (2016).

In this way, for the F phase, the modes with $L=0$ couple to just the $\hat{S}_{-}$ operator, with $\hat{S}_{\pm}=\hat{S}_{x}\pm i\hat{S}_{y}$, while the $L=1,L_{z}=0,\pm 1$ modes couple to the compound operator $\hat{S}_{-}\hat{L}_{z,\pm}$, respectively. Thus, in principle, only the $L=0$ modes can be detected using their coupling to a magnetic field in the $x,y$ direction. For the FLP phase, the discussion is exactly the same but changing the role of spin and layer operators, which makes much harder their detection as they all couple to the $\hat{L}_{-}$ operator.

On the other hand, for the phases with spontaneously broken symmetries, the situation is much easier due to the lower number of conserved quantities. For instance, for the CAF phase, in the sector $L_{z}=0$, the modes $+$ couple to the $\hat{S}_{\pm}$ operators while the modes $-$, including the Goldstone mode, couple to $\hat{S}_{z}$: as pointed out in Ref. Murthy et al. (2016), the coupling of the CAF phase with the $\hat{S}_{z}$ operator could be used in principle to unambiguously distinguish this phase from the others. The modes with $L=1,L_{z}=\pm 1$ couple to the $\hat{S}_{z}\hat{L}_{\pm}$ operators. For the PLP phase, the modes with $S=0$ couple to just the $\hat{L}_{z}$ operator while the modes with $S=1,S_{z}=0,\pm 1$ couple to the operators $\hat{L}_{z}\hat{S}_{z,\pm}$, respectively. Then, in a similar way to the CAF phase, the coupling of the modes to a layer voltage provides a unique feature of the PLP phase that can be used to clearly characterize it.

An alternative way to study the collective-mode spectra through measurements of thermal transport has been proposed in a very recent work Pientka et al. (2017).

Regarding the detection of the dynamical instabilities, from Eqs. (IV.2.3), (IV.2.4) it is seen that the conditions for the appearance of unstable modes for both CAF and PLP phases are

[TABLE]

Both conditions of instability are depicted in Fig. 9. Interestingly, the saturation of the two previous inequalities correspond to the boundary between the PLP and CAF phases for $\epsilon_{V}=0$ ($\epsilon_{Z}=0$). The first condition is unlikely to be achieved for the expected values of the coupling constants. The second condition is at least compatible with those values, but corresponds to a region where the PLP is not the actual ground state, as can be seen in the right plot. Naturally, the system is not dynamically unstable when starting from an equilibrium state. However, these potentially unstable dynamics could be explored in out of equilibrium situations.

For monolayer graphene, it is well known that the transport gap is governed by skyrmions Yang et al. (2006), scaling as $\sim\sqrt{B_{\perp}}$ for both bare and screened interactions Kharitonov (2012a), in agreement with the experimental results Jiang et al. (2007). The same scaling occurs for the mean-field transport gap since the screened Coulomb contribution also behaves in monolayer graphene as $\bar{F}_{00}=\hbar\omega_{B}h\left(\frac{F_{C}}{\hbar\omega_{B}}\right)$, with $h$ some dimensionless function. However, as now $F_{C}/\hbar\omega_{B}$ is independent of $B_{\perp}$, the gap always scales as $\sim\hbar\omega_{B}\sim\sqrt{B_{\perp}}$ for any value of $B_{\perp}$. With respect to the detection of the collective modes, as the valley-spin structure is formally equal to the bilayer scenario, the above discussion still holds although, unfortunately, there is not a direct equivalent of the layer voltage term. Alternative candidates for the job are distortions in the lattice Abanin et al. (2007); Fuchs and Lederer (2007) but they seem to be much less controllable from an experimental point of view.

VIII Conclusions

In this work, we have studied the $\nu=0$ QH state within a mean-field Hartree-Fock approach projected onto the zero-energy Landau level. We have reproduced the mean-field phase diagram of the $\nu=0$ QH state of Ref. Kharitonov (2012a) by solving the self-consistent Hartree-Fock equations and computed the different mean-field energies and transport gaps.

After that, using the time-dependent Hartree-Fock approximation, we have computed the corresponding collective modes, that in the limit of zero momentum can be obtained analytically. In particular, we have provided a complete classification of the rich orbital-valley-spin structure of the modes and their symmetries. As a result, we have identified the singlet modes in orbital pseudospin at zero momentum as those characterizing the phase transitions since they are the lowest energy modes and the only ones coupled to the many-body contributions of short-range interactions at low momentum.

With respect to the valley-spin structure of the modes, at the boundary between the ferromagnetic and the fully layer-polarized phases, there is a gapless mode, resulting from a mean-field symmetry that can be regarded as a weak extension of the complete $SO(5)$ symmetry described in Ref. Wu et al. (2014); the same phenomenon appears at the boundary between the canted anti-ferromagnetic and the partially layer-polarized phases. On the other hand, Goldstone modes appear at the boundaries between the ferromagnetic and canted anti-ferromagnetic phases and the fully and partially layer-polarized phases due to the spontaneous symmetry breakings. It is remarkable that, while the former Goldstone mode corresponds to a triplet mode in valley pseudospin, the latter is a spin-singlet mode. This contrast arises due to the different nature of both symmetry breaking mechanisms. A complementary study of the phase transitions could be provided in terms of the edge modes, extending the work presented in Refs. Murthy et al. (2016); Tikhonov et al. (2016) for the ferromagnetic-canted anti-ferromagnetic transition in monolayer graphene.

Another interesting feature of the phases with spontaneously broken symmetries is that they are able to present dynamical instabilities within this projected model. The study of such unstable dynamics in out of equilibrium scenarios is of particular interest as both the in-plane magnetic field and the interlayer voltage are easily manipulable in the laboratory. For instance, one can try to explore the occurrence of dynamical instabilities for a state initially prepared in the phases with spontaneously broken symmetries by applying a sudden quench in the external fields de Nova et al. .

The performed calculations for bilayer graphene can be straightforwardly adapted to monolayer graphene, due to the formal analogy of the Hamiltonian. We have found that most of the conclusions of the previous paragraphs still hold in the monolayer scenario. We have also analyzed the effects of LL mixing and we have accounted them by screening the long-range Coulomb interaction and by renormalizing the coupling constants of the short-range interactions. The resulting effective Hamiltonian is formally similar to that previously considered and hence, the phase diagram and the collective modes present the same structure discussed above.

Finally, we have made some experimental remarks about the transport gap and the detection of the instabilities. We have found that, within the current detection schemes, one can only aim at characterizing a restricted number of modes. Specifically, among all the modes describing the phase transitions, only those related to the F-CAF and FLP-PLP transitions couple to physical observables in the side of the spontaneously broken symmetry phase. Hence, the detection of the collective modes motivates the study and design of new possible experimental methods able to provide effective observables that couple to valley-pseudospin waves and measure its dispersion relation. In this sense, the recent work of Ref. Hunt et al. (2016) provides a promising line of research as it shows measurements of the layer polarization of the states.

Acknowledgements.

We would like to specially thank E. Demler for his collaboration at the beginning the project. We also thank E. Shimsoni, P. Tikhonov, L. Brey, F. Guinea, F. Sols and T. Stauber for useful discussions. This work has been supported by MINECO (Spain) through grants FIS2010-21372 and FIS2013-41716-P, Comunidad de Madrid through grant MICROSERES-CM (S2009/TIC-1476) and also partially by the Israel Science Foundation.

Appendix A Solutions of the Hartree-Fock equations

We compute in this Appendix the solutions to the Hartree-Fock equations for the $\nu=0$ QH state. We start by reviewing the case of the 2D electron gas (2DEG), as many of the developed techniques are also used for the actual computation for graphene, presented later.

A.1 Basic results of the 2DEG

As a first step, we define the magnetic form factors

[TABLE]

where the functions $\phi_{n}(x)$ are the usual harmonic oscillator wave functions

[TABLE]

In order to compute the explicit expression of the magnetic form factors, we rewrite Eq. (114) as:

[TABLE]

and we recall that $P_{x}$, the momentum operator in the $x$-direction, is the generator of translations in the $x$ direction. Rewriting the previous expression in terms of the usual harmonic oscillator destruction operator $a$ in the $x$-direction yields

[TABLE]

Then, after some straightforward manipulations, we obtain

[TABLE]

with $L^{m}_{n}(x)$ a generalized Laguerre polynomial. From the above relations it follows that $A_{nn^{\prime}}\left(\mathbf{k}\right)=A^{*}_{n^{\prime}n}\left(-\mathbf{k}\right)$. It is worth studying the dependence on $\mathbf{k}$ of the magnetic form factors. If we switch to the following polar coordinates, $(k_{x},k_{y})=k(\sin\varphi_{\mathbf{k}},\cos\varphi_{\mathbf{k}})$, we find that

[TABLE]

As in practice we restrict to the ZLL, we give the explicit expression of the magnetic form factors involved in the calculations:

[TABLE]

Interestingly, in the context of the Wigner function, the magnetic form factors are the Moyal functions of the harmonic oscillator Schleich (2001).

For simplicity, we first review the usual case of integer QH states in the 2DEG, where the field operator only has two components corresponding to the spin polarizations $\xi=\pm$. The extension to the graphene scenario is discussed in the next subsection. The non-interacting eigenfunctions of the 2DEG are:

[TABLE]

where $\phi_{n,k}(\mathbf{x})$ are the orbital wave functions of the magnetic levels and $\chi_{\xi}$ is the spin wave-function with polarization $\xi$. The non-interacting eigenvalues of these wave functions are

[TABLE]

with $\omega_{B}$ and $\epsilon_{Z}$ the corresponding cyclotron and Zeeman energy.

We now suppose that the electrons interact through an arbitrary scalar potential $V$ (as it is, for instance, the Coulomb interaction $V_{0}$). Its spatial matrix elements are given in terms of the magnetic form factors:

[TABLE]

The HF equations associated to the previous interacting potential are (see Appendix B.1.1 for a derivation of the HF equations)

[TABLE]

with $\nu_{m,\xi^{\prime}}$ the occupation number of each LL and the components of the vector $\mathbf{\pi}=(\pi_{x},\pi_{y})$ given after Eq. (17). We analyze now the spatial structure of the previous equation. For that purpose, we define two associated mean-field potentials, denoted as $V^{H,F}$, that take into account the Hartree (direct) and Fock (exchange) contributions from the potential $V$, respectively. First, we address the non-self-consistent (NSC) problem, where the mean-field potentials are created by the bare (non-interacting) wave functions $\phi^{0}_{n,k,\xi}$. In that case, the above HF equations are:

[TABLE]

where $V_{m}^{H,F}$ are the contributions from each magnetic level to the mean-field Hartree and Fock potentials:

[TABLE]

The function $K_{n}(\mathbf{x},\mathbf{x^{\prime}})\equiv K_{nn}(\mathbf{x},\mathbf{x^{\prime}})$ is the spatial part of the Green’s function, with

[TABLE]

We can compute explicitly $K_{nn^{\prime}}(\mathbf{x},\mathbf{x^{\prime}})$ by transforming the discrete sum over $p$ into an integral and by making the transformation $x_{p}=X+pl^{2}_{B}$, with the center of mass and relative coordinates defined as $\mathbf{R}\equiv(\mathbf{x}+\mathbf{x^{\prime}})/2$ and $\Delta\mathbf{x}\equiv\mathbf{x}-\mathbf{x^{\prime}}$. The result is

[TABLE]

In particular,

[TABLE]

with $\Delta r=|\Delta\mathbf{x}|$. Since $K_{nn^{\prime}}(\mathbf{x},\mathbf{x})=\frac{\delta_{nn^{\prime}}}{2\pi l^{2}_{B}}$, $(2\pi l^{2}_{B})^{-1}$ being the homogeneous density corresponding to a completely filled magnetic level, the Hartree contribution is uniform:

[TABLE]

$V(0)$ being the Fourier transform of the potential evaluated at $\mathbf{k}=0$. In the usual case of the Coulomb interaction, the Hartree potential is canceled by the positive charge background. For a short-range interaction, proportional to $\delta(\mathbf{x}-\mathbf{x^{\prime}})$, the spatial part of the Hartree potential is equal to the Fock potential.

With respect to the Fock potential, for general purposes, we consider the following matrix elements

[TABLE]

The case of the Fock potential in Eq. (A.1) is obtained for $m=m^{\prime}$. To compute the previous integral, we switch to the center of mass and relative coordinates defined before Eq. (129) and integrate along the center of mass coordinates $X,Y$, obtaining:

[TABLE]

Now, we assume the typical situation where the potential $V(\mathbf{x})$ is rotationally invariant (as the Coulomb potential, for example). The polar dependence described in Eq. (119) gives the result

[TABLE]

with $F_{nmm^{\prime}}$ the on-shell value of the previous integral. Specifically, for the Fock potential, $m=m^{\prime}$, we get

[TABLE]

Note that $F_{mn}=F_{nm}$. Amazingly, the Hartree and Fock potentials of Eq. (A.1) are diagonal in the magnetic base of the eigenfunctions of Eq. (121), which implies that the self-consistent orbital wave functions are indeed the same as the non-interacting ones. Moreover, from Eq. (125), it is easy to check that this result also holds for the spin part. Thus, the complete self-consistent wave functions are equal to the non-interacting ones.

It is useful to reproduce the previous results in Fourier space by writing the matrix elements of the Hartree and Fock potentials in terms of the matrix elements of the potential $V$ in Eq. (123). For instance, for the Hartree potential we trivially find:

[TABLE]

On the other hand, for the generalized Fock potential of Eq. (132), we obtain

[TABLE]

Again, if we consider that the potential is rotationally invariant, its Fourier transform is also rotationally invariant and hence we recover the result of Eq. (133). In particular, for $m=m^{\prime}$ we find

[TABLE]

Equations (134) and (137) are related to each other through the identity:

[TABLE]

which can be proven by inserting the definition of the magnetic form factors, Eq. (114).

A.2 Diagonalization of the Hartree-Fock equations in graphene

After the previous training, we obtain the HF solutions for the $\nu=0$ QH in bilayer graphene (the monolayer case is just a trivial extension of this calculation). We start by computing the eigenfunctions of the single-particle Hamiltonian (13). Its matrix elements, $\braket{nk\alpha}{H_{0}}{n^{\prime}k^{\prime}\alpha^{\prime}}$, with $\braket{\mathbf{x}}{nk\alpha}=\Psi^{0}_{n,k,\alpha}(\mathbf{x})$ the wave functions of Eqs. (20)-(22), are given by:

[TABLE]

with $(\tau_{i})_{\alpha\alpha^{\prime}}=\chi^{\dagger}_{\alpha}\tau_{i}\chi_{\alpha^{\prime}}$. We see that the previous Hamiltonian is diagonal within the ZLL while it mixes the LLs $\pm n$ due to the layer voltage. The reason is that, outside the ZLL, the wave functions are not localized on one specific layer and thus they are able to experiment the effect of the voltage. However, as long as $\epsilon_{V}\ll\hbar\omega_{B}$, it is a good approximation to consider that $H_{0}$ is diagonal in every LL so their corresponding eigenfunctions are still given by $\Psi^{0}_{n,p,\alpha}(\mathbf{x})$. The associated eigenvalues $H_{0}\ket{nk\alpha}=\hbar\omega^{0}_{n,\alpha}\ket{nk\alpha}$ are then

[TABLE]

where $\lambda,\xi$ label the valley and spin polarizations, $\lambda=\pm 1$ corresponding to $K,K^{\prime}$ valleys. The non-interacting spinors are given by all possible orthogonal combinations of valley-spin polarizations, $\chi_{\alpha}^{0}=\ket{\pm n_{z}}\otimes\ket{\pm s_{z}}$; see Eq. (39) for the notation of the wave-functions in valley-spin space. For the $\nu=0$ QH state, all the states with $n\leq-2$ are filled and empty for $n\geq 2$ while the ZLL is half-filled. As the non-interacting energy only depends on the polarization in valley-spin space, the two magnetic levels are filled in the same way in valley-spin space. Thus, following the notation of the main text, we label the occupied non-interacting spinors as $\chi^{0}_{a,b}$ and the empty ones as $\chi^{0}_{c,d}$. In particular, for $\epsilon_{V}<\epsilon_{Z}$, the ZLL is filled in the same F configuration of Eq. (39) and for $\epsilon_{V}>\epsilon_{Z}$ it is filled in the FLP configuration of Eq. (40). Hence, in the non-interacting problem, these are the only two possible phases for the $\nu=0$ QH state.

We now address the effect of the interactions, given by Eqs. (24), (25). The resulting HF equations for the $\nu=0$ QH state are:

[TABLE]

where the Hartree term of the Coulomb potential is canceled by the positive charge background. As in the 2DEG case, we first consider the NSC problem in order to understand the structure of the equations.

A.2.1 Non-self consistent problem

The spatial part of the non-interacting Green’s function is given by a $2\times 2$ matrix in the $\bar{A}\bar{B}$ subspace

[TABLE]

Operating in the same fashion as in the 2DEG case, we obtain:

[TABLE]

In particular, for $\mathbf{x}=\mathbf{x^{\prime}}$, $K^{(2)}_{n}(\mathbf{x},\mathbf{x^{\prime}})=(4\pi l^{2}_{B})^{-1}\text{diag}[1,1]$ for $|n|\neq 0,1$ and $K^{(2)}_{n}(\mathbf{x},\mathbf{x^{\prime}})=(2\pi l^{2}_{B})^{-1}P_{\bar{B}}$ for $n=0,1$, $P_{\bar{B}}=\text{diag}[0,1]$ being the projector onto the subspace $\bar{B}$.

The matrix $K^{(2)}_{n}(\mathbf{x},\mathbf{x^{\prime}})$ is key to understand the orbital structure of the HF equations. With that finality, we write the NSC version of the HF equations (A.2) as:

[TABLE]

and separate the different contributions to the non-self-consistent mean-field Hamiltonian $H_{NSCHF}$:

[TABLE]

where $H_{0}$ is just the non-interacting Hamiltonian. The term $H_{NSCDS}$ represents the NSC mean-field potential created by the Dirac sea, formed by all the non-interacting states with $n\leq-2$,

[TABLE]

Since all valley-spin polarizations are occupied in the Dirac sea, the Hartree term of the short-range interactions vanishes and the corresponding Fock term is just an scalar. The last term of Eq. (150) corresponds to the NSC mean-field potential created by the filled states of the ZLL:

[TABLE]

with $P^{0}=\chi^{0}_{a}\chi^{0\dagger}_{a}+\chi^{0}_{b}\chi^{0\dagger}_{b}$ the analog of the matrix $P$ [defined after Eq. (33)] for the non-interacting spinors, $P_{\bar{A}}=\text{diag}[1,0]$ the projector onto the subspace $\bar{A}$ and $v_{i}=(g_{ix}+g_{iy})/\pi l^{2}_{B}$. The matrix elements of the Fock contribution of the Coulomb potential are

[TABLE]

with $C_{nm},D_{nm}$ some coefficients that can be expressed in terms of the values of the generalized matrix elements of Eq. (133); their particular expression is not interesting for the current purposes. As we see, the Fock term of the Coulomb potential mixes $n$ with $-n$ outside the ZLL and it is diagonal within the ZLL, in the same fashion of the non-interacting Hamiltonian. It is easy to show that this behavior also persists in the matrix elements of the short-range terms; therefore, we conclude that the NSC Hamiltonian $H_{NSCHF}$ only couples the LLs $\pm n$ and leaves invariant the ZLL. Moreover, even when the mean-field potentials are self-consistently renormalized, this structure is preserved. Hence, the orbital part of the self-consistent wave functions $\Psi_{n,k,\alpha}$ is constructed from the non-interacting orbital wave functions $\Psi^{0}_{\pm n,k}$. In the following, we neglect the coupling between $\pm n$ LLs as we are assuming that interactions are weak enough to neglect LL mixing. A nice discussion about the coupling of the $\pm n$ LLs for strong Coulomb interactions is presented in Ref. Tőke (2013).

A.2.2 Self-consistent problem

We now switch to the actual self-consistent problem. We rewrite Eq. (A.2) as previously done for the NSC problem

[TABLE]

and separate the different contributions along the same lines of Eq. (150)

[TABLE]

As we are neglecting LL mixing, we can regard the Dirac sea as frozen and hence the self-consistent HF mean-field potential created by the Dirac sea, $H_{DS}$, is the same as the NSC one of Eq. (A.2.1), $H_{DS}=H_{NSCDS}$. With respect to the self-consistent mean-field potential created by the filled states of the ZLL, Eq. (152), we only have to replace the NSC projector $P^{0}$ by the self-consistent one, $P$, as the filling of the ZLL is of the same form [see discussion after Eq. (III)]. Hence, we focus on the ZLL in order to determine the self-consistent valley-spin wave functions of the ZLL.

For that purpose, since $H_{HF}$ leaves invariant the ZLL, we only need to take the components of $H_{HF}$ in the sublattice $\bar{B}$, obtaining the projected HF Hamiltonian $H^{(0)}_{HF}$ for the ZLL, $H^{(0)}_{HF}=P_{\bar{B}}H_{HF}P_{\bar{B}}=H^{(0)}_{0}+H^{(0)}_{DS}+H^{(0)}_{ZLL}$. $H^{(0)}_{0}$ corresponds to just the layer-voltage and Zeeman terms, $H^{(0)}_{0}=-\epsilon_{V}\tau_{z}-\epsilon_{Z}\sigma_{z}$, as the kinetic energy is zero. The interaction of the ZLL with the inert Dirac sea can be taken into account by a scalar non-local potential $V_{DS}(\mathbf{x},\mathbf{x^{\prime}})$,

[TABLE]

while the self-consistent interaction inside the ZLL gives rise to

[TABLE]

Note that, after projecting, the interaction with the Dirac sea corresponds to a one-body operator.

Then, the resulting HF equation for the ZLL

[TABLE]

is completely equivalent to the HF equation (III) derived from the effective Hamiltonian (29) and gives

[TABLE]

where we have reabsorbed in the Hamiltonian the infinite energy shift provided by the sum of the term $\bar{u}\delta(\mathbf{x}-\mathbf{x^{\prime}})$ in the expression of $V_{DS}(\mathbf{x},\mathbf{x^{\prime}})$. Since the mean-field HF potentials are invariant under unitary transformations that leave invariant the subspace formed by the filled states [see Eq. (172) and related discussion], the self-consistent problem is still diagonal in the orbital part as the two magnetic levels $n=0,1$ are filled exactly in the same way. Thus, the self-consistent wave functions have the form $\Psi_{n,k,\alpha}(\mathbf{x})=\Psi^{0}_{n,k}(\mathbf{x})\chi_{\alpha}$.

In order to arrive at an equation for the spinor $\chi_{\alpha}$, we multiply by $\Psi^{0,\dagger}_{n,k}(\mathbf{x})$ and integrate in Eq. (A.2.2), obtaining:

[TABLE]

with $F_{nm}$ the eigenvalues of the Fock potential associated to the Coulomb interaction, Eq. (134). The contribution from the mean-field interaction with the Dirac sea, which is the origin of the analog of the Lamb shift, corresponds to the first term at the r.h.s. of the above equation. Following the regularization procedure of Ref. Shizuya (2012), we rearrange the associated series as:

[TABLE]

The completeness relation

[TABLE]

can be proven from the definition of the magnetic form factors Eq. (114) and the completeness relation of the harmonic oscillator wave functions. Then, using Eqs. (137), (162), we find that

[TABLE]

which represents a constant (infinite) energy shift that can be absorbed in the Hamiltonian, as that arising from the short-range interactions [see discussion below Eq. (156)]. After this regularization, and by defining $F_{n}\equiv F_{n0}+F_{n1}$ as in the main text, we finally get Eq. (33).

The formalism developed in this section is also valid for screened Coulomb interactions if we just replace $V_{0}(k)$ by $\bar{V}(k)$ of Eq. (107). The expression for the screened Fock energies $\bar{F}_{nm}$ then takes the form:

[TABLE]

with $g_{nm}(x)$ the dimensionless functions obtained from the integral

[TABLE]

In particular, the dimensionless function $g(x)$ associated to the screened Coulomb contribution to the transport gap, $\bar{F}_{1}=\bar{F}_{10}+\bar{F}_{11}$ [see Eq. (111)], reads:

[TABLE]

For small $x$, one recovers the unscreened result for $g(x)$ while for large $x$, $g(x)$ diverges logarithmically since $f(z)\simeq\alpha z^{2}$ for small $z$.

Appendix B Time-dependent Hartree-Fock approximation

We review in this Appendix the basic theory of the TDHFA. First, we consider a variational approach for the wave function Thouless (1961) and then we connect the results with a diagrammatic calculation of the correlation functions Kallin and Halperin (1984); Wang et al. (2002), which gives identical results for the computation of the collective modes. Alternative approaches can also be seen in Refs. Giuliani and Vignale (2005); Negele and Orland (2008).

B.1 Variational formalism

B.1.1 Stationary situation

We start by reviewing the stationary HF equations in the usual case of a system of $N$ fermions governed by a second-quantization Hamiltonian of the form Fetter and Walecka (2003):

[TABLE]

where the eigenfunctions of the single-particle Hamiltonian $H_{sp}$ are known. The indices $l,k,j,m$ label the states of an orthonormal basis of the single-particle Hilbert space. The stationary HF equations are obtained by looking for a Slater determinant that minimizes the expectation value of the Hamiltonian

[TABLE]

$\lambda=1\ldots N$ being $N$ orthogonal states. Now, we consider small particle-hole perturbations around the single-particle states in the so-called Thouless parametrization Thouless (1961):

[TABLE]

where the sum in $\Lambda$ runs over all the unoccupied states. Along this section, we will use lower Greek indices $\lambda,\sigma$ for labeling occupied levels and upper Greek indices $\Lambda,\Sigma$ for empty levels, while Latin indices $l,k,j,m$ will label arbitrary (filled or empty) states.

The condition for $\ket{\Psi_{0}}$ to be an extreme of the energy, $\bra{\Psi_{0}}\hat{H}\ket{\Psi_{0}}$, reads:

[TABLE]

which implies that $H^{HF}_{\Lambda\lambda}=0$, with

[TABLE]

Thus, if $\ket{\Psi_{0}}$ is an extreme, the mean-field Hartree-Fock Hamiltonian $\hat{H}^{HF}$ created by $\ket{\Psi_{0}}$ cannot connect occupied and empty states. Moreover, as $\hat{H}^{HF}$ is invariant under unitary transformations that leave invariant the subspace spanned by the occupied and empty states, we can choose an orthonormal basis of the Hilbert space that diagonalizes $\hat{H}^{HF}$ so

[TABLE]

which corresponds to the usual HF equations for the self-consistent wave functions written in terms of matrix elements. The Hartree-Fock Hamiltonian is then

[TABLE]

The total energy of the state $\ket{\Psi_{0}}$ is:

[TABLE]

In order to check if $\ket{\Psi_{0}}$ is a true minimum of the expectation value of $\hat{H}$, we consider fluctuations around $\ket{\Psi_{0}}$ [as given by Eq. (B.1.1)] and keeping only terms up to second order in the coefficients $w_{\Lambda\lambda}$, we get the quadratic form

[TABLE]

with $\mathbf{W}$ being a column vector containing the coefficients $w_{\Lambda\lambda},w^{*}_{\Lambda\lambda}$,

[TABLE]

The elements of the matrix $X$ are

[TABLE]

Then, if the matrix $X$ is positive definite, $\ket{\Psi_{0}}$ corresponds to a true local minimum of the expectation value of the Hamiltonian and hence it is an energetically stable solution. We will refer to the elements $X_{\lambda\Lambda,\sigma\Sigma}$, $X_{\Lambda\lambda,\Sigma\sigma}$ as normal and to the elements $X_{\lambda\Lambda,\Sigma\sigma}$, $X_{\Lambda\lambda,\sigma\Sigma}$ as anomalous since they are not given by a matrix elements between two particle-hole excitations but rather between the ground state and an excited state with two pairs of particle-holes.

For notational convenience, we rewrite the matrix $X$ in terms of arbitrary states $kl,jm$ as:

[TABLE]

where $\nu_{k},\nu_{l}$ represent the number occupation of the state $k,l$. We keep in mind that the only valid matrix elements of $X$ in Eq. (178) are those with the pair index $kl$ corresponding to one level filled and one level empty, so $\nu_{k}-\nu_{l}=\pm 1$ if $k$ is filled (empty) and $l$ is empty (filled); this also applies for the pair $jm$. In the same fashion, we rewrite the components of the vector $\mathbf{W}$ as $W_{kl}=w_{lk}$ for $\nu_{k}-\nu_{l}=1$ and $W_{kl}=w^{*}_{kl}$ for $\nu_{k}-\nu_{l}=-1$.

Using this notation, we can rewrite the quadratic form of Eq. (B.1.1) in a more compact way:

[TABLE]

where ′ denotes now that we only sum over the proper values of the pair indices $kl,jm$.

B.1.2 Time-dependent situation

The Schrödinger equation for the time evolution of the wave function can also be derived from a variational principle. In particular, a solution of the Schrödinger equation, $\ket{\Psi(t)}$, must be an extreme of the functional

[TABLE]

The equations for the TDHFA arise when imposing a time-dependent solution of the same form of Eq. (B.1.1),

[TABLE]

In that case, we find

[TABLE]

where we have only kept terms up to second order in $w_{\Lambda\lambda}$. The resulting equation of motion for $w_{\Lambda\lambda}$ is:

[TABLE]

Combining this equation of motion with that of $f(t)$, it is straightforward to show that:

[TABLE]

Thus, $f^{*}\frac{df}{dt}$ is quadratic in the coefficients $w_{\Lambda\lambda}$ and since we must consistently keep only the lowest order terms for the equation of motion for $w_{\Lambda\lambda}$, Eq. (183) reduces to the linear $f$-independent equation

[TABLE]

If we now perform a linear expansion in modes for $w_{\Lambda\lambda}$,

[TABLE]

$\gamma_{n}$ being the amplitude of each mode, we finally arrive at the equations for the TDHFA

[TABLE]

In matrix form, the previous equation simply reads as

[TABLE]

with $N,A$ denoting the normal, anomalous sectors of the matrix $X$. As $N$ is an Hermitian matrix, whenever the anomalous elements are non-zero, Eq. (188) is a non-Hermitian eigenvalue equation and thus, it can present complex eigenvalues. We also remark that the $u$ and $v$ components are mixed only for non-vanishing anomalous elements. It is worth noting the strong formal analogy of the above equations with the bosonic Bogoliubov-de Gennes equations Pitaevskii and Stringari (2003); Pethick and Smith (2008).

For convenience, we rewrite the above equations using the compact notation developed in Eqs. (178),(179). In this way, Eq. (185) can be put as

[TABLE]

with the elements of the matrix $\tilde{T}$ given by $\tilde{T}_{kl,jm}=(\nu_{k}-\nu_{l})\delta_{kj}\delta_{lm}$. The expansion in modes of Eq. (186) reads:

[TABLE]

and the equations of the TDHFA, Eq. (B.1.2), are simply

[TABLE]

Exploiting the analogies with the bosonic BdG equations, we find some interesting properties of the TDHFA equations. For instance, if $\mathbf{Z}_{n}$ is a mode with frequency $\omega_{n}$, then the conjugate $\bar{\mathbf{Z}}_{n}$ is as well a mode with eigenvalue $-\omega^{*}_{n}$.

[TABLE]

The TDHFA equations also have an associated Klein-Gordon type scalar product given by:

[TABLE]

In particular, if $\mathbf{Z}_{n},\mathbf{Z}_{m}$ are two eigenmodes of Eq. (190), then

[TABLE]

from which follows that two modes with different eigenvalues are orthogonal according to this scalar product. Another important property is that the previous scalar product is not positive definite so the norm of a given solution $\mathbf{Z}_{n}$, defined as $(\mathbf{Z}_{n}|\mathbf{Z}_{n})$, can be positive, negative or zero. In particular, the norm of the conjugate $\bar{\mathbf{Z}}_{n}$ has the opposite sign to that of $\mathbf{Z}_{n}$. It is immediately deduced from Eq. (199) that any mode with complex frequency has zero norm. Also, as $\tilde{T}\tilde{X}^{\dagger}\tilde{T}=\tilde{X}$, the complex-frequency modes appear in pairs $(\omega_{n},\omega_{n}^{*})$.

Interestingly, we can further relate the energetic stability of the HF solution $\ket{\Psi_{0}}$ (discussed at the end of Sec. B.1.1) with the frequencies of the collective modes computed in the TDHFA through

[TABLE]

If the state $\ket{\Psi_{0}}$ truly minimizes the expectation value of the Hamiltonian, the matrix $X$ is positive definite and then there are no dynamical instabilities since that would imply the presence of a mode with zero norm, $(\mathbf{Z}_{n}|\mathbf{Z}_{n})=0$. Therefore, the presence of dynamical instabilities is only possible if the state is energetically unstable.

Also, if $\ket{\Psi_{0}}$ is energetically stable, the modes with positive (negative) frequency have positive (negative) normalization. Thus, the presence of a mode with positive (negative) normalization and negative (positive) frequency reveals that the system is energetically unstable. Hence, in practice, we only need to compute the modes with positive norm since the rest of the modes are obtained through the conjugation $\mathbf{Z}_{n}\rightarrow\bar{\mathbf{Z}}_{n}$. Indeed, they are the two sides of the same coin as they both give rise to the same wave function.

More interesting properties appear in the presence of a continuous symmetry in the Hamiltonian. A continuous symmetry arises when the Hamiltonian is invariant under a continuous transformation of the form $\hat{U}^{\dagger}(\varphi)\hat{H}\hat{U}(\varphi)=\hat{H}$, with $\hat{U}(\varphi)=e^{-i\varphi\hat{L}}$, $\hat{L}$ being a Hermitian operator of the form

[TABLE]

As well known, the presence of such a symmetry gives rise to the appearance of a gapless Goldstone mode when the symmetry is spontaneously broken by the ground state. This result is recovered in the TDHFA by taking into account that the ground state is a solution of the HF equations. In that case, if we rewrite Eq. (170) as

[TABLE]

and expand the exponential for small $\varphi$, we find that a mode with zero frequency arises, $X\mathbf{Z}_{G}=0$, with

[TABLE]

$L_{\Lambda\lambda}$ being the matrix elements of $\hat{L}$ that connect filled and empty states. This gapless mode disappears whenever all the $L_{\Lambda\lambda}$ are zero; in that case, the unitary transformation given by $\hat{U}(\varphi)$ leaves invariant the state $\ket{\Psi_{0}}$ so the symmetry is not spontaneously broken.

The above argument also predicts the existence of such a gapless mode whenever the system presents a continuous mean-field symmetry, by which we denote a continuous transformation that is not a exact symmetry of the Hamiltonian, $\hat{U}^{\dagger}(\varphi)\hat{H}\hat{U}(\varphi)\neq\hat{H}$, but it satisfies Eq. (202), which means that $\hat{U}(\varphi)\ket{\Psi_{0}}$ is also a mean-field ground state with the same energy of $\ket{\Psi_{0}}$.

Using the developed formalism, it is straightforward to compute the response function to a small perturbation introduced by an external field $H_{\rm{ext}}(t)$ within the TDHFA. Adding $H_{\rm{ext}}$ to the equilibrium Hamiltonian $\hat{H}$ in the functional of Eq. (180), using the same ansatz of Eq. (181) and retaining only the lowest order terms in the amplitude of the external field gives an inhomogeneous version of Eq. (189)

[TABLE]

Thus, to first order, the system only responds to the particle-hole matrix elements of the external field, perturbing only the expectation value of particle-hole operators. Then, if we consider the change to lowest order in the expectation value of a particular observable $A$ due to the coupling with an external field $H_{\rm{ext}}=B$, $\Delta A(t)=\sideset{}{{}^{\prime}}{\sum}_{kl}A_{kl}W_{kl}(t)$, we find by taking the Fourier transform in the above expression that

[TABLE]

with $\mathbf{A},\mathbf{B}$ vectors containing the particle-hole matrix elements of $A,B$, $\mathbf{A}_{kl}=A^{*}_{kl}$ and the same for $\mathbf{B}$. The function $\chi(\omega)$ is called the mother of all response functions Giuliani and Vignale (2005). Its poles are given by the condition

[TABLE]

which after a trivial transformation gives the TDHFA eigenvalue equation $\det~{}[\tilde{X}-\hbar\omega]=0$, completely equivalent to Eq. (196). Hence, the poles of the different response functions are indeed the frequencies of the collective modes, as well known from the general theory of the response function Fetter and Walecka (2003).

B.1.3 TDHFA for the $\nu=0$ quantum Hall state in graphene

We now apply the previous formalism to the study of the $\nu=0$ QH state in graphene, although many of the results presented proceed from the general theory of integer QH states. In that case, we take the Hamiltonian of Eq. (167) as the effective Hamiltonian of Eq. (29). The states of the Hilbert space are labeled by the dummy index $k$, which represents all the quantum numbers $(n_{k},p_{k},\alpha_{k})$ at the same time, with $n_{k}$ the magnetic number, $p_{k}$ the $y$-momentum and $\alpha_{k}$ the polarization in valley-spin space. In particular, in our ZLL projected Hamitonian, $n_{k}=0,1$ and $\alpha_{k}=a,b,c,d$. The single-particle Hamiltonian $H_{sp}$ corresponds to the layer voltage, the Zeeman term and the mean-field interaction with the Dirac sea while the total effective potential corresponds to the Coulomb and short-range interactions. Since we are restricting to the ZLL, where the orbital part of the wave functions is proportional to the magnetic wave function $\phi_{n}$, the matrix elements of the total effective interaction potential [$V_{lk,jm}$ in Eq. (167)] can be straightforwardly calculated with the help of Eq. (123) as

[TABLE]

The corresponding HF equations (172) read in terms of the self-consistent wave functions as Eq. (III). As discussed in the main text and in Appendix A.1, the orbital part of the self-consistent wave functions is equal to that of the non-interacting wave functions and the ground state $\ket{\Psi_{0}}$ is given by Eq. (32).

With respect to the TDHFA equations (196), since the occupation number does not depend on the $y$-momentum for integer QH states, we perform the following transformation in momentum space for an arbitrary vector $\mathbf{Z}$:

[TABLE]

so we get rid of the momentum coordinates and obtain a discrete matrix equation only in the magnetic and valley-spin indices:

[TABLE]

$\mathbf{Z}(\mathbf{k})$ being a vector with components $Z_{kl}(\mathbf{k})$ and $\tilde{T}$ the correspondent version of the same matrix in only magnetic and valley-spin indices. The elements of the matrix $X(\mathbf{k})$, $X_{kl,jm}(\mathbf{k})\equiv X_{n_{k}\alpha_{k}n_{l}\alpha_{l},n_{j}\alpha_{j}n_{m}\alpha_{m}}(\mathbf{k})$, are computed from

[TABLE]

which gives:

[TABLE]

The energies $U^{RPA},U^{LAD}$ result from direct (exchange) interactions; see also discussion after Eq. (230). The matrix $X(\mathbf{k})$ satisfies the properties:

[TABLE]

Following Ref. Wang et al. (2002), we will refer to the matrix $X(\mathbf{k})$ as the dispersion matrix.

As a result, the transformation given by Eq. (210) diagonalizes the TDHFA equations so the wave vector $\mathbf{k}$ becomes a good quantum number that represents the momentum of the magnetoexciton wave function, Eq. (48). By defining $\hat{M}^{\dagger}_{kl}\equiv\hat{M}^{\dagger}_{n_{k}\alpha_{k}n_{l}\alpha_{l}}(\mathbf{k})$, the matrix elements of $X(\mathbf{k})$ can be easily expressed in terms of magnetoexciton matrix elements as

[TABLE]

with $X_{\Lambda\lambda,\sigma\Sigma}(\mathbf{k})$, $X_{\lambda\Lambda,\Sigma\sigma}(\mathbf{k})$ the anomalous matrix elements.

The wave functions associated to the collective modes computed from Eq. (209), $\ket{\Psi(\mathbf{k},t)}=\hat{M}^{\dagger}(\mathbf{k},t)\ket{\Psi_{0}}$, are characterized by the particle-hole operators $\hat{M}^{\dagger}(\mathbf{k},t)$, which are given in terms of linear combinations of magnetoexcitons

[TABLE]

as can be seen from Eqs. (186), (B.1.3).

If all the $v_{\Lambda\lambda}$ coefficients are zero, the excitation is a pure magnetoexciton with momentum $\mathbf{k}$ while for non-vanishing $v_{\Lambda\lambda}$ the excitations are combinations of magnetoexcitons with momentum $\pm\mathbf{k}$. This can be understood from the fact that, as the $v$ components arise from the anomalous matrix elements (connecting the ground state with a state with two magnetoexcitons created), in order to conserve total momentum the two magnetoexcitons must sum total momentum zero, see Eq. (B.1.3).

With respect to the response function, after inverting relation (210), one gets from Eq. (205) that

[TABLE]

with the elements $\mathbf{A}_{kl}(\mathbf{k}),\mathbf{B}_{kl}(\mathbf{k},\omega)$ defined as in Eq. (B.1.3). The mother of all response functions is now diagonalized in momentum and reads $\chi(\mathbf{k},\omega)=\left[\hbar\omega\tilde{T}-X(\mathbf{k})\right]^{-1}$. Of special interest is the case where the operators $A,B$ are of the form:

[TABLE]

with $\theta_{A}$ some operator in valley-spin space, as they characterize the response of the system to the introduction of a perturbation in the charge, spin or interlayer density. After some manipulations, it is seen that

[TABLE]

B.2 Diagrammatic formalism

Following Refs. Kallin and Halperin (1984); Wang et al. (2002), we rederive the results of the previous section using a diagrammatic expansion. For definiteness, we restrict from the beginning to the particular case of the $\nu=0$ QH state (although the extrapolation of the results for general integer QH states is immediate). The first step is to compute the self-consistent Green’s function, whose equation is given by the diagrammatic representation of the first line of Fig. 10. As well known Fetter and Walecka (2003), this equation leads to the HF equations for the self-consistent wave functions, Eq. (III), in terms of which the Green’s functions is written:

[TABLE]

with $\eta=0^{+}$ and $\hbar\omega_{k}=\epsilon_{k}$ the HF energy. We follow here the notation introduced in Sec. B.1.3 so the dummy index $k$ labels all the quantum numbers $(n_{k},p_{k},\alpha_{k})$.

Since the corresponding formalism is closely related to that of the response function Fetter and Walecka (2003), the collective-mode frequencies can also obtained from the poles of the set of correlation functions:

[TABLE]

where $x=(\mathbf{x},t)$, $T$ denotes time-ordering, all the expectation values are evaluated in the ground state of the system in the Heisenberg picture and $\Delta\hat{O}=\hat{O}-\braket{\hat{O}}$ denotes the fluctuations of an operator around its mean value.

The Fourier transform of the correlation function gives:

[TABLE]

after taking into account that the correlation function is invariant under time translations (due to the fact that the Hamiltonian is time independent) and also under spatial translations (this last property is shown explicitly later).

We now proceed to compute the previous correlation functions within the diagrammatic version of the TDHFA, in which the equation for $\Pi_{A}(\mathbf{k},\omega)$ reads

[TABLE]

where we have introduced the dressed vertex function $\Lambda^{A}_{kl}(\mathbf{k^{\prime}},\omega)$. In the above equation, we have made use of the identity $\sum_{\mathbf{k^{\prime}}}e^{i\mathbf{k^{\prime}}\mathbf{x^{\prime}}}=S\delta(\mathbf{x^{\prime}})$. We remark that the index $k$ in Eq. (222) labels the quantum numbers $(n_{k},p_{k},\alpha_{k})$ and not momentum. The expression for the non-interacting vertex $\Lambda^{A(0)}_{kl}(\mathbf{k},\omega)$ is:

[TABLE]

with $\theta_{\alpha_{l}\alpha_{k}}$ defined in Eq. (B.1.3).

In the TDHFA, the Dyson’s equation for the dressed vertex function is given by:

[TABLE]

with $V_{lk,jm}$ the matrix elements of the total interaction potential, see Eq. (207).

The diagrams for the total bare interaction, the correlation function and the vertex equation in the TDHFA are shown in Fig. 10. Hence, within the TDHFA, the dressed vertex function is the bare vertex (first term) plus the series corresponding to bubble (second term) and ladder (third term) diagrams. The RPA approximation would correspond to just keep only the bubble diagrams. By defining the two-particle propagator

[TABLE]

we simplify Eqs. (222), (224) as

[TABLE]

We note that the two-particle propagator $D_{kl}(\omega)$ does not depend on the momenta $p_{k},p_{l}$ and it is only non-zero whenever the pair index $kl$ represents one filled level and one empty.

Further simplifications can be made by summing over all the momenta

[TABLE]

The function $L^{A}_{kl}$ only depends in $k,l$ through the magnetic level and the valley-spin polarization so, thanks to these manipulations, we get rid of the momentum indices and simplify the vertex equation. For instance, for the bare vertex, $\Lambda^{A(0)}_{kl}(\mathbf{k^{\prime}},\omega)$, one finds

[TABLE]

while for the ladder and bubble diagrams, after using the expression for the matrix elements given in Eq. (207),

[TABLE]

with $W_{lk,jm}(\mathbf{k})$ given in Eq. (211). From the above expression, it is easy to understand the notation of the quantities $U^{RPA},U^{LAD}$: they take into account the energy contribution from bubble and ladder diagrams, respectively Wang et al. (2002). Since we have performed the summation in all the momentum dummy variables, from now on we use the index $k$ to label just the quantum numbers $(n_{k},\alpha_{k})$.

After inserting the expression for $\Lambda^{A}_{kl}(\mathbf{k},\omega)$ of Eq. (226) in Eq. (228) and using Eqs. (229)-(230), one realizes that $L^{A}_{kl}(\mathbf{k},\mathbf{k^{\prime}},\omega)$ is diagonal in $(\mathbf{k},\mathbf{k^{\prime}})$. Hence, after defining a new dressed vertex function:

[TABLE]

the original integral equation (224) for the dressed vertex is transformed into a discrete matrix equation:

[TABLE]

where we remark that the labels $l,k,j,m$ now only represent the pair index corresponding to the magnetic level $n$ and the valley-spin polarization $\alpha$. The correlation function then reads:

[TABLE]

The $\delta_{\mathbf{k},\mathbf{k^{\prime}}}$ factor appearing in Eq. (231) ensures the translational invariance of the correlation function previously assumed. In order to solve the matrix equation (B.2) we define:

[TABLE]

and then

[TABLE]

with $\mathbf{a}_{kl}$ given by Eq. (B.1.3). According to Eqs. (225), (234) and (235), the only valid matrix elements of $X_{kl,jm}(\mathbf{k},\omega)$ are those where the index pair $kl$ correspond to one level occupied and the other level empty (the same goes for $jm$). Therefore, $X_{kl,jm}(\mathbf{k},\omega=0)=X_{kl,jm}(\mathbf{k})$ is exactly the same of Eq. (211).

Indeed, after rewriting Eq. (233) as

[TABLE]

and following the discussion leading to Eq. (206), we find that the poles of the correlation function also give the collective-mode energies, as they are given by the condition $\det~{}[\tilde{X}(\mathbf{k})-\hbar\omega]=0$, which is the same eigenvalue problem of Eq. (209). In particular, note that Eqs. (B.1.3), (237) are equivalent as $X^{-1}(\mathbf{k},\omega)=-\chi(\mathbf{k},\omega)$, revealing the close link between the correlation and the response functions mentioned at the beginning of the section.

Although the diagrammatic formalism may be more robust from a fundamental point of view, the form of the TDHFA equations obtained in Sec. B.1.3 provides a simpler and clearer physical picture of the TDHFA since it allows us to understand the collective modes in terms of wave functions.

Appendix C Computation of the dispersion relation in the TDHFA

We devote this section to the computation of the solutions to the TDHFA equations for the $\nu=0$ QH state within the projected model considered in this work.

C.1 Analytical results for the dispersion matrix

We start by computing analytically the elements of the dispersion matrix, $X_{kl,jm}(\mathbf{k})$, given by Eq. (211). In particular, we focus on the non-trivial many-body contribution arising from interactions, $W_{lk,jm}(\mathbf{k})$. In order to give its explicit expression, we separate the Coulomb from the short-ranged terms

[TABLE]

The Coulomb term, $W^{C}_{lk,jm}(\mathbf{k})$, has no direct (RPA) contribution as it vanishes when considering the proper elements of the dispersion matrix because in our model the empty levels have different valley-spin polarization from that of the filled ones. Thus, we only have to consider the exchange (ladder) contribution, which yields:

[TABLE]

These energies are given in terms of modified Bessel functions, as seen by computing the element $U^{C}_{00,00}(\mathbf{k})$. After changing to polar coordinates, taking into account that only the real part of the complex exponential survives, performing the integral over the radial coordinate and using the integral representation of the modified Bessel functions,

[TABLE]

one finds that:

[TABLE]

with $F_{00}$ given by Eq. (34).

The other elements $U^{C}_{n_{j}n_{k},n_{l}n_{m}}(\mathbf{k})$ can also be computed analytically. Instead of using the complicated expression found in the tables, we simply take into account that $A_{n_{l}n_{m}}(\mathbf{q})$ is a polynomial in $q_{x},q_{y}$ multiplied by $e^{-\frac{(ql_{B})^{2}}{4}}$. Then,

[TABLE]

with $P_{n_{j}n_{k},n_{l}n_{m}}(q_{x},q_{y})$ some polynomial. For $U^{C}_{00,00}(\mathbf{k})$, $P_{00,00}(q_{x},q_{y})=1$. Therefore, from the usual properties of Fourier transforms, one has that

[TABLE]

Since the Bessel functions are solutions of a second order differential equation, their higher order derivatives can always be put in terms of themselves and their first derivatives. For $I_{0}(x)$, $\frac{dI_{0}}{dx}=I_{1}(x)$, which means that the elements $U^{C}_{n_{j}n_{k},n_{l}n_{m}}(\mathbf{k})$ are expressed through combinations of $I_{0},I_{1}$ and polynomials in $k_{x},k_{y}$. For instance, from Eq. (120), we obtain

[TABLE]

and so on. For large $x$, a saddle-point approximation of the integral in Eq. (240) gives $I_{0}(x)\sim\frac{e^{x}}{\sqrt{x}}$. Then, for large wave vector $k$, $kl_{B}\gg 1$, $U^{C}_{00,00}(\mathbf{k})\sim(kl_{B})^{-1}$. As the other elements are obtained from derivatives of $U^{C}_{00,00}(\mathbf{k})$, the total matrix $U^{C}$ decays at least as $(kl_{B})^{-1}$ for large $k$.

We now turn our attention to the short-range interactions. Using the property (138), one finds that the orbital part of RPA and ladder contributions is the same and then:

[TABLE]

Thus, all the valley-spin structure of the short-range interactions is captured by the effective couplings $u^{\alpha_{k}\alpha_{l},\alpha_{j}\alpha_{m}}$, which present the following properties:

[TABLE]

In order to simplify the notation, we define a matrix $C(\mathbf{k})$ containing the exchange Coulomb energies with the following index ordering in orbital space:

[TABLE]

We also define an analog matrix $R(\mathbf{k})$ for the short-range interactions

[TABLE]

Then, after taking into account that the only valid matrix elements of $X_{kl,jm}(\mathbf{k})$ are those with the pair indices $kl,jm$ corresponding to one level filled and one level empty and that the occupation number for the $\nu=0$ QH state only depends on the valley-spin polarization, it is seen that the matrix $\tilde{X}(\mathbf{k})$ of the TDHFA equations (209) is a $32\times 32$ matrix of the form:

[TABLE]

where the $4\times 4$ matrices $Y^{\alpha_{k}\alpha_{l},\alpha_{j}\alpha_{m}}(\mathbf{k})$ are the building blocks of $\tilde{X}(\mathbf{k})$ and represent the different valley-spin sectors. As discussed after Eq. (200), computing just the modes with positive norm is sufficient for characterizing the collective modes.

After some straightforward manipulations, the matrices $Y^{\alpha_{k}\alpha_{l},\alpha_{j}\alpha_{m}}(\mathbf{k})$ can be written as:

[TABLE]

We remind that the occupied states ($\nu_{\alpha}=1$) correspond to polarizations $\alpha=a,b$ while the empty ones ($\nu_{\alpha}=0$) correspond to $\alpha=c,d$ and that $I$ is the $4\times 4$ identity matrix. Between the square brackets of Eq. (250), $\Delta^{\alpha_{k}\alpha_{l}}$ represents the valley-spin part of the mean-field energy gap [see Eq. (36)] while the matrix $F(\mathbf{k})$ contains all the terms involving Coulomb interactions, both mean-field and many-body contributions. At the same time, they are all multiplied by a diagonal tensor with respect to the valley-spin indices, so the different valley-spin sectors of $\tilde{X}$ are only connected through terms arising from short-range interactions, proportional to $u^{\alpha_{k}\alpha_{l},\alpha_{j}\alpha_{m}}$. This fact simplifies notably the calculations since many of the effective couplings $u^{\alpha_{k}\alpha_{l},\alpha_{j}\alpha_{m}}$ are related to each other through Eq. (246) and also some of them vanish due to symmetry considerations (see next section) and consequently we do not need to take into account the full $32\times 32$ problem. In particular, using Eq. (246), it is shown that the only independent and non-vanishing element in the anomalous sector [the off-diagonal boxes in Eq. (249)] is $u^{ac,db}$.

Some analytical properties of the dispersion relation can be obtained from the above results. First, as explained in the main text, at $\mathbf{k}=0$ the total orbital pseudospin $O$ and its $z$-component $O_{z}$ are good quantum numbers. We can check this fact by making a unitary transformation in the matrix $\tilde{X}$ that switches to the singlet and triplet magnetoexciton base, $\hat{M}^{\dagger}_{OO_{z},\alpha\alpha^{\prime}}(\mathbf{k})$,

[TABLE]

This transformation yields a simplified form at $\mathbf{k}=0$ for the matrices $Y^{\alpha_{k}\alpha_{l},\alpha_{j}\alpha_{m}}$,

[TABLE]

where the order of the magnetic indices for rows and columns corresponds to $OO_{z}=00,11,1-1,10$. All matrices are now diagonal in this new basis. In particular, the energy of the singlet $OO_{z}=00$ mode is independent of the Coulomb strength and is the only one that presents a non-trivial structure in valley-spin due to many-body contributions from short-range interactions, represented by the matrix $R$. This contrasts to the case of orbital-triplet modes, where the many-body corrections due to short-range interactions vanish and the collective-mode frequencies are immediately obtained from the diagonal of $\Delta^{\alpha_{k}\alpha_{l}}I+F$.

For $\mathbf{k}\neq 0$, we show that the resulting dispersion relation is isotropic. The key point is the polar structure of the magnetic form factors, shown in Eq. (119). As $(k_{x},k_{y})=k(\sin\varphi_{\mathbf{k}},\cos\varphi_{\mathbf{k}})$,

[TABLE]

with $R_{n_{k}n_{l},n_{j}n_{m}}(k)$ a function that only depends on $k=|\mathbf{k}|$ [not to be confused with the subindex $k$ that labels the quantum numbers $(n_{k},\alpha_{k})$]. For the Coulomb interaction, we switch to the following polar coordinates in the integral of Eq. (C.1), $\mathbf{q}=q(\sin\varphi_{\mathbf{q}},\cos\varphi_{\mathbf{q}})$. This gives

[TABLE]

with $w_{n_{j}n_{k},n_{l}n_{m}}(q)$ some function that only depends on $q$. The polar integral gives

[TABLE]

This implies that the matrix $C(\mathbf{k})$ has the same dependence on the polar angle of Eq. (253). On the other hand, for the diagonal elements $n_{k}n_{l}=n_{j}n_{m}$ and the factor $e^{i(n_{l}-n_{k})\varphi_{\mathbf{k}}}e^{i(n_{j}-n_{m})\varphi_{\mathbf{k}}}$ becomes the unity. Hence, $\tilde{X}_{kl,jm}(\mathbf{k})$ is of the form

[TABLE]

with $X_{kl,jm}(k)$ depending only on $k=|\mathbf{k}|$. Then, after making an appropriated phase transformation in the magnetic indices we get rid of the dependence on the polar angle of the momentum and obtain a matrix $\tilde{X}_{kl,jm}(k)$ which explicitly depends solely on $k=|\mathbf{k}|$. Therefore, the resulting dispersion relation is isotropic, $\omega(\mathbf{k})=\omega(k)$.

We remark that the previous proof only relies on the fact that the Coulomb interaction is rotationally invariant and not on its particular form. Thus, this result applies for any rotationally invariant interaction potential; in particular, it holds for the screened interaction considered in Sec. VI.

For completeness, we give the expression of the rotationally invariant matrices $F(k),R(k)$ that arise from the matrix $\tilde{X}(k)$ in Eq. (256). For Coulomb interactions,

[TABLE]

and $F(k)$ obtained through Eq. (250), while for short-range interactions

[TABLE]

Note that the matrix $R(k)$ is a matrix of rank $1$ as it is proportional to the projector on the vector $a(k)$. Interestingly, from the above expressions it is immediate to check that the vector $[0,1,1,0]^{T}$ is an eigenvector of the matrices $F(k),R(k)$; in particular, it is orthogonal to $a(k)$. Hence, the following linear combination of magnetoexcitons with $O_{z}=\pm 1$ and $\nu_{\alpha_{k}}-\nu_{\alpha_{l}}=1$,

[TABLE]

is always an eigenmode of the TDHFA equations, which describes the $N=2$ orbital modes, with frequency

[TABLE]

Note that the only dependence on the short-range interactions of this mode is contained in the valley-spin contribution of the single-particle gap, $\Delta^{\alpha_{k}\alpha_{l}}$, so the dependence in $k$ solely involves Coulomb interactions.

Finally, we note that considering screened Coulomb interactions only amounts to replace the matrix $F(k)$ by its screened version, $\bar{F}(k)$.

C.2 Computation of the collective modes in bilayer graphene

We now discuss the details of the computation of the collective modes for the different states of Eqs. (39),(40),(III.3), (44). For all the phases, we give the independent non-vanishing coefficients $u^{\alpha_{k}\alpha_{l},\alpha_{j}\alpha_{m}}$ that characterize the valley-spin structure, discuss the symmetries of the modes and the associated eigenvalue problem and compute the stiffness coefficients for all the modes.

C.2.1 Ferromagnetic phase

[TABLE]

The problem is diagonalized in valley-spin space by considering magnetoexcitons with well defined valley pseudospin numbers, $\hat{M}^{\dagger}_{nn^{\prime},LL_{z}}(\mathbf{k})$:

[TABLE]

Since the ground state has well-defined quantum number for $S,S_{z},L,L_{z}$ and all magnetoexcitons have $S=1,S_{z}=-1$, the anomalous matrix element vanish due to spin conservation and hence dynamical instabilities cannot appear in the ferromagnetic state. As a consequence, for a fixed value of $LL_{z}$, the frequency and the orbital structure of the modes (corresponding to the $N=0,1,2,3$ orbital modes) are computed from the $4\times 4$ eigenvalue equation

[TABLE]

where we have made use of the isotropic expressions (256)-(263). All the matrices $Y^{\mu=LL_{z}}(k)$ present the same structure:

[TABLE]

where the valley-spin gaps $\Delta^{\mu}$ and effective coupling energies $u^{\mu}$ are

[TABLE]

The previous relations imply that the dispersion relation for the triplet modes with $L_{z}=\pm 1$ is the same, only shifted by the layer voltage, $\omega^{1-1}(k)=\omega^{11}(k)+4\epsilon_{V}$, and that $Z_{1-1}(k)=Z_{11}(k)$. The wave function of each mode is created by the operator

[TABLE]

Using standard techniques it is straightforward to compute from Eq. (269) the stiffness coefficients of the different modes:

[TABLE]

The expression for the stiffness of the $N=2$ modes can also be obtained directly from Eq. (265). Note that the stiffness coefficients of the modes $N=2,3$ are independent of short-range interactions and hence of the valley-spin symmetry of the mode; this result holds for all the phases. Remarkably, $\rho^{\mu}_{1}<0$ for every mode as the leading contribution goes as

[TABLE]

C.2.2 Full layer-polarized phase

[TABLE]

As this phase is the analog of the F phase in valley space, the results for the FLP phase are formally analog to those of the F phase and they are obtained by replacing $LL_{z}$ by $SS_{z}$ in Eqs. (267)-(C.2.1). As mentioned in the main text, all triplet modes present the same dispersion relation, only shifted by the Zeeman energy, $\omega^{1\pm 1}(k)=\omega^{10}(k)\mp 2\epsilon_{Z}$ and $Z_{11}(k)=Z_{10}(k)=Z_{1-1}(k)$.

C.2.3 Canted anti-ferromagnetic phase

[TABLE]

Since the CAF state has only well defined value of the $z$ valley pseudospin, $L_{z}=0$, the anomalous element $u^{ac,db}$ becomes non-zero, so the collective modes are combinations of magnetoexcitons with $\pm\mathbf{k}$, see Eq. (215) and related discussion. For the same reason, complex-frequency modes can appear.

The eigenvalue problem is split in sectors with fixed $L_{z}$. First, we consider the subspace with $L_{z}=0$, where we define:

[TABLE]

Note that, due to the different spin orientations of $s_{a},s_{b}$, these modes have no well-defined value of $L$; only in the limit $\theta_{S}=0$ they match the proper $L=0,1$ modes. In this new basis, the eigenvalue problem is diagonalized in valley-spin space and reduced to two independent $8\times 8$ blocks:

[TABLE]

where the reader is advised to pay attention to the interchange of magnetic indices in the $v$ components of the vector $Z$. The matrices $Y^{\pm}(k)$ are given by Eq. (269) with the values $\Delta^{\pm}=\Delta^{ac}$ and $u^{\pm}=u^{ac,ac}\pm u^{ac,bd}$ and the matrix $A(k)$ characterizes the anomalous matrix elements of $\tilde{X}(k)$, $A(k)\equiv u^{ac,db}R(k)$.

The velocity of the Goldstone mode is computed, after some tedious but straightforward algebra [note that the eigenvalue problem of Eq. (277) at $k=0$ gives a non-diagonal Jordan canonical form for the orbital-singlet modes rather than the usual diagonal matrix], as

[TABLE]

The wave function of the modes is created by the operator:

[TABLE]

For the components with $L_{z}=\pm 1$, the eigenvalue problem reads:

[TABLE]

where the expressions for $Y^{1\pm 1}(k)$ are the same as in the F state, and the operator associated to the wave function of the collective modes is:

[TABLE]

which we see that mixes magnetoexcitons with $\pm\mathbf{k},~{}L_{z}=\pm 1$ in order to conserve both total momentum and total value of $L_{z}$ in the anomalous matrix element. As in the F phase, $\omega^{1-1}_{N}(k)=\omega^{11}_{N}(k)+4\epsilon_{V}$ and $Z_{1-1}(k)=Z_{11}(k)$.

The expression of the stiffness coefficients for the orbital-singlet modes (when $\mu$ does not correspond to the Goldstone one) is:

[TABLE]

For the $N=1$ modes, the expression varies if there is a Goldstone mode:

[TABLE]

We note that the leading contribution in the Coulomb interaction strength for $\rho^{\mu}_{1}$ still satisfies Eq. (273).

C.2.4 Partially layer-polarized phase

[TABLE]

The valley-spin structure of the modes is similar to that of the FLP phase and the problem is diagonalized by considering modes with well defined spin number $S,S_{z}$. The main difference is that, as the Hamiltonian does not commute with the operator $\hat{L}_{n}$, the anomalous matrix element is non-zero.

For the modes with $S_{z}=0$, we find that the collective modes are given by:

[TABLE]

The corresponding wave functions are those of Eq. (C.2.3) but changing $\pm$ by $SS_{z}=00,10$, respectively. The exact expression for the velocity of the Goldstone mode, $\omega^{00}_{0}(k)\simeq v_{G}k$, is also given by Eq. (281). For the modes with $S_{z}=\pm 1$, the results are formally analog to those of Eqs. (C.2.3), (286), but with $SS_{z}$ playing the role of $LL_{z}$. The values of $\Delta^{SS_{z}},~{}u^{SS_{z}}$ characterizing the matrices $Y^{SS_{z}}(k)$ are computed in the same way as in the FLP case, also finding that $\omega^{1\pm 1}(k)=\omega^{10}(k)\mp 2\epsilon_{Z}$ and $Z_{11}(k)=Z_{10}(k)=Z_{1-1}(k)$.

The expression of the stiffness coefficient for the orbital-singlet mode $\rho^{10}_{0}$ is given by Eq. (C.2.3) and those of $\rho^{10}_{1},\rho^{00}_{1}$ by the first and second line of Eq. (288), respectively.

C.3 Computation of the collective modes in monolayer graphene

The TDHFA previously developed for computing the collective modes of the $\nu=0$ QH state of bilayer graphene is also valid for the monolayer. In particular, the dispersion relation is also obtained from the eigenvalues of the matrix $\tilde{X}(\mathbf{k})$ of Eq. (249) but now there is only one possible value for the magnetic index, $n=0$. Thus, as the structure in the valley-spin space is the same as for the bilayer problem, the building blocks of the dispersion matrix $\tilde{X}$, $Y^{\alpha_{k}\alpha_{l},\alpha_{j}\alpha_{m}}(\mathbf{k})$, are now scalars instead of $4\times 4$ matrices:

[TABLE]

where we have made explicit the isotropy of the dispersion relation. All the coefficients $\Delta^{\alpha_{k}\alpha_{l}}$ and $u^{\alpha_{k}\alpha_{l},\alpha_{j}\alpha_{m}}$ have the same value as in the bilayer case but with $\epsilon_{V}=0$; the same holds for the magnitudes $\Delta^{\mu}$ and $u^{\mu}$. Hence, the dispersion relation is computed in the same way as in the previous section but replacing all the matrices $Y$ by the corresponding scalars, notably simplifying the calculations.

We now give the explicit analytical results for the dispersion relation of the collective modes $\omega^{\mu}(k)$ and the associated stiffness $\rho^{\mu}$ for every phase, characterized by the same valley-spin symmetries $\mu$ as in the bilayer case. Note that the components $u_{nn^{\prime}},v_{nn^{\prime}}$ of the magnetoexciton wave functions are reduced now to scalars as the orbital structure is trivial for monolayer graphene.

C.3.1 Ferromagnetic and charge-density wave phases

[TABLE]

As $\epsilon_{V}=0$, for the F phase $\omega^{1-1}(k)=\omega^{11}(k)$. Note the strong similarity in the formula for the stiffness with that of bilayer graphene for dominant Coulomb interactions, Eq. (61).

C.3.2 Canted anti-ferromagnetic and Kekulé distortion phases

[TABLE]

Similarly, $\omega^{1-1}(k)=\omega^{11}(k)$ in the CAF phase.

The velocity of the Goldstone modes is simply:

[TABLE]

with $\Delta^{G}=\Delta^{-}$ for the CAF phase and $\Delta^{G}=\Delta^{00}$ for the KD phase, while the stiffness coefficients are:

[TABLE]

Both results can be compared with those for the bilayer for dominant Coulomb interactions, Eqs. (65), (71).