Ground-state atlas of a three-dimensional semimetal in the quantum limit

TL;DR

This paper provides a comprehensive theoretical phase diagram of a three-dimensional semimetal in the quantum limit, revealing various ground states influenced by electron interactions and screening effects, with implications for experimental observations.

Contribution

It introduces an unbiased theoretical framework combining parquet RG and mean-field methods to map ground states of 3D semimetals under high magnetic fields, considering both repulsive and attractive interactions.

Findings

Ground state varies from excitonic insulator to spin density wave depending on screening.

Topological excitonic insulator supports massless Dirac fermions on the surface.

Surface SdH oscillations exhibit a $ oot H_{\perp}$ dependence.

Abstract

An interplay between electron correlation and reduced dimensionality due to the Landau quantization gives rise to exotic electronic phases in three-dimensional semimetals under high magnetic field. Using an unbiased theoretical method, we clarify for the first time comprehensive ground-state phase diagrams of a three-dimensional semimetal with a pair of electron and hole pockets in the quantum limit. For the electron interaction, we consider either screened Coulomb repulsive interaction or an attractive electron-electron interaction mediated by a screened electron-phonon coupling, where a screening length is generally given by a dimensionless constant times magnetic length $l$. By solving the parquet RG equation numerically and employing a mean-field argument, we construct comprehensive ground-state phase diagrams of the semimetal in the quantum limit for these two cases, as a function of the Fermi wave length and the screening length (both normalized by $l$). In the repulsive interaction case, the ground state is either excitonic insulator (EI) in strong screening regime or Ising-type spin density wave in weak screening regime. In the attractive interaction case, the ground state is either EI that breaks the translational symmetries (strong screening regime), topological EI, charge Wigner crystal (intermediate screening regime), plain charge density wave or possible non-Fermi liquid (weak screening regime). We show that the topological EI supports a single copy of massless Dirac fermion at its side surface, and thereby exhibit a $\sqrt{H_{\perp}}$-type surface Shubnikov-de Haas (SdH) oscillation in in-plane surface transports as a function of a canted magnetic field $H_{\perp}$. Armed with these theoretical knowledge, we discuss implications of recent transport experiments on graphite under the high field.

Equations

$$V(q_{\perp },p_z\simeq 0)$$

$$V(q_{\perp },p_z\simeq 2k_F)$$

$$\displaystyle V_{\rm eff}(q_{\perp},p_{z})=-\frac{\rho_{0}}{Mc^{2}}\bigg{(}\frac{4\pi Ze^{2}l^{2}}{\varepsilon}\frac{1}{(q^{2}_{\perp}+p^{2}_{z})l^{2}+\frac{1}{A}e^{-\frac{1}{2}q^{2}_{\perp}l^{2}}}\bigg{)}^{2}$$

$${\hat{H}}_T=\int d^3\mathbf{r}{\hat{h}}_0(\mathbf{r})+\frac{1}{2}\int d^3\mathbf{r}{d}^3{\mathbf{r}}^{\prime }\hat{\rho }(\mathbf{r})\hat{\rho }({\mathbf{r}}^{\prime })V(\mathbf{r}-{\mathbf{r}}^{\prime })$$

$${\hat{h}}_0(\mathbf{r})\equiv \sum _{\tau =+(\uparrow ),-(\downarrow )}$$

$$\displaystyle\hskip 2.84544pt\bigg{[}\Psi^{\dagger}_{e,\tau}({\bm{r}})\Big{\{}-E_{g}+\frac{(-\hbar^{2}\nabla^{2}_{z}+{\bm{\pi}}^{2})}{2m_{e}}-H_{\rm z}\tau\Big{\}}\Psi_{e,\tau}({\bm{r}})$$

$$\displaystyle\hskip 8.5359pt+\Psi^{\dagger}_{h,\tau}({\bm{r}})\Big{\{}E_{g}-\frac{(-\hbar^{2}\nabla^{2}_{z}+{\bm{\pi}}^{2})}{2m_{h}}-H_{\rm z}\tau\Big{\}}\Psi_{h,\tau}({\bm{r}})\bigg{]}$$

$$\displaystyle(\pi_{x},\pi_{y})=\left\{\begin{array}[]{c}(-i\hbar\partial_{x},-i\hbar\partial_{y}+eHx)\\ (-i\hbar\partial_{x}-eHy,-i\hbar\partial_{y})\\ \end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{c}E^{e}_{n,\tau}(k_{z})=-E_{g}+\frac{\hbar^{2}p^{2}_{z}}{2m_{e}}-H_{z}\tau+\big{(}n+\frac{1}{2}\big{)}\hbar\omega_{e}\\ E^{h}_{n,\tau}(k_{z})=E_{g}-\frac{\hbar^{2}p^{2}_{z}}{2m_{h}}-H_{z}\tau-\big{(}n+\frac{1}{2}\big{)}\hbar\omega_{h}\\ \end{array}\right.$$

$$\int d^3\mathbf{r}{\hat{h}}_0(\mathbf{r})-\mu \hat{N}=2\pi l\sum _{\sigma =\pm }{\int }_{-\Lambda }^{+\Lambda }\frac{dp}{2\pi }\sigma v_{F}p$$

$$\displaystyle\hskip 5.69046pt\int d(lQ)\Big{\{}e^{\dagger}_{\sigma}(Q,p)e_{\sigma}(Q,p)-h^{\dagger}_{\sigma}(Q,p)h_{\sigma}(Q,p)\Big{\}}$$

$${\Psi }_{e,\uparrow }(\mathbf{r})$$

$$+{\int }_{-\Lambda }^{+\Lambda }\frac{dp}{2\pi }{e}^{i(-k_F+p)z}\int d(l{q}_x){\varphi }_{q_x}(x,y)e_{-}(q_x,p),$$

$${\Psi }_{h,\downarrow }(\mathbf{r})$$

$$+{\int }_{-\Lambda }^{+\Lambda }\frac{dp}{2\pi }{e}^{i(-k_F+p)z}\int d(l{k}_y){\psi }_{k_y}(x,y)h_{-}(k_y,p),$$

$${\psi }_{k_y}(x,y)$$

$${\varphi }_{q_x}(x,y)$$

$${\int }_{-\infty }^{\infty }\frac{d(l{k}_y)}{\sqrt{2\pi }}{e}^{i{q}_x{k}_y{l}^2}{\psi }_{k_y}(x,y)={\varphi }_{q_x}(x,y),$$

$${\int }_{-\infty }^{\infty }\frac{d(l{q}_x)}{\sqrt{2\pi }}{e}^{-i{q}_x{k}_y{l}^2}{\varphi }_{q_x}(x,y)={\psi }_{k_y}(x,y).$$

$$\frac{1}{2}\int d^3\mathbf{r}{d}^3{\mathbf{r}}^{\prime }\hat{\rho }(\mathbf{r})\hat{\rho }({\mathbf{r}}^{\prime })V(\mathbf{r}-{\mathbf{r}}^{\prime })\equiv H_1+H_2,$$

$$H_1\equiv \frac{1}{2}\sum _{\mu ,\nu =e_{+},e_{-},h_{+},h_{-}}\int \frac{dp}{2\pi }\int d(l{Q}_1)\int d(l{Q}_1^{\prime })\int d(l{Q}_2^{\prime })\int d(l{Q}_2)$$

$${\Gamma }_{\mu \nu }(Q_1,Q_1^{\prime },Q_2^{\prime },Q_2;I_0)\int \frac{d{p}_2}{2\pi }{a}_{\mu }^{†}(Q_1^{\prime },p_2+p)a_{\mu }(Q_2^{\prime },p_2)\int \frac{d{p}_1}{2\pi }{a}_{\nu }^{†}(Q_1,p_1-p)a_{\nu }(Q_2,p_1),$$

$$H_2\equiv \sum _{\mu ,\nu =e,h}\int \frac{dp}{2\pi }\int d(l{Q}_1)\int d(l{Q}_1^{\prime })\int d(l{Q}_2^{\prime })\int d(l{Q}_2)$$

$${\Phi }_{\mu \nu }^{+-}(Q_1,Q_1^{\prime },Q_2^{\prime },Q_2;I_{2k_F})\int \frac{d{p}_2}{2\pi }{a}_{{\mu }_{+}}^{†}(Q_1^{\prime },p_2+p)a_{{\mu }_{-}}(Q_2^{\prime },p_2)\int \frac{d{p}_1}{2\pi }{a}_{{\nu }_{-}}^{†}(Q_1,p_1-p)a_{{\nu }_{+}}(Q_2,p_1),$$

$$\displaystyle\left\{\begin{array}[]{c}a_{e_{+}}(Q,p)\equiv e_{+}(k_{y},p),\\ a_{e_{-}}(Q,p)\equiv e_{-}(q_{x},p),\\ a_{h_{+}}(Q,p)\equiv h_{+}(q_{x},p),\\ a_{h_{-}}(Q,p)\equiv h_{-}(k_{y},p).\\ \end{array}\right.$$

$$I_0(q_x,k_y)$$

$$I_{2k_F}(q_x,k_y)$$

$${\Gamma }_{e_{+}{e}_{+}}(k_1,k_1^{\prime },k_2^{\prime },k_2;I_0)=\cdots ={\Gamma }_{h_{-}{e}_{+}}(k_1,k_1^{\prime },k_2^{\prime },k_2;I_0)$$

$$=\delta (k_1+k_1^{\prime }-k_2-k_2^{\prime })\int d{q}_x$$

$$I_0(q_x,-k_1+k_2)e^{-i\frac{1}{2}{q}_x(k_1+k_2-k_1^{\prime }-k_2^{\prime })l^2},$$

$${\Gamma }_{e_{-}{e}_{-}}(q_1,q_1^{\prime },q_2^{\prime },q_2;I_0)=\cdots ={\Gamma }_{h_{+}{e}_{-}}(q_1,q_1^{\prime },q_2^{\prime },q_2;I_0)$$

$$=\delta (q_1+q_1^{\prime }-q_2-q_2^{\prime })\int d{k}_y$$

$$I_0(q_1-q_2,k_y)e^{i\frac{1}{2}{k}_y(q_1+q_2-q_1^{\prime }-q_2^{\prime })l^2},$$

$${\Gamma }_{e_{+}{e}_{-}}(q_1,k_1,k_2,q_2;I_0)=\cdots ={\Gamma }_{h_{-}{e}_{-}}(q_1,k_1,k_2,q_2;I_0)$$

$$=e^{i{\mathbf{k}}_1\wedge {\mathbf{k}}_2{l}^2}{I}_0(q_1-q_2,-k_1+k_2),$$

$${\Gamma }_{e_{-}{e}_{+}}(k_1,q_1,q_2,k_2;I_0)=\cdots ={\Gamma }_{e_{-}{h}_{-}}(k_1,q_1,q_2,k_2;I_0)$$

$$=e^{i{\mathbf{k}}_1\wedge {\mathbf{k}}_2{l}^2}{I}_0(q_1-q_2,-k_1+k_2),$$

$$\displaystyle{\bm{k}}_{1}\wedge{\bm{k}}_{2}\equiv\left(\begin{array}[]{c}k_{1}\\ q_{1}\\ \end{array}\right)\wedge\left(\begin{array}[]{c}k_{2}\\ q_{2}\\ \end{array}\right)=k_{1}q_{2}-q_{1}k_{2}.$$

$${\Phi }_{ee}^{+-}(q_1,k_1,q_2,k_2;I_{2k_F})={\Phi }_{hh}^{+-}(k_1,q_1,k_2,q_2;I_{2k_F})$$

$$=l^2{e}^{i{\mathbf{k}}_1\wedge {\mathbf{k}}_2{l}^2}\int \frac{d{q}_{x}d{k}_y}{2\pi }{I}_{2k_F}(q_x,k_y)e^{i(k_1-k_2)q_x{l}^2-i(q_1-q_2)k_y{l}^2},$$

$${\Phi }_{eh}^{+-}(k_2,k_1,q_1,q_2;I_{2k_F})=l^2{e}^{i(k_1{q}_1+k_2{q}_2)l^2}\int \frac{d{q}_{x}d{k}_y}{2\pi }$$

$$I_{2k_F}(q_x,k_y)e^{i{q}_x{k}_y{l}^2+i{q}_x(k_1-k_2)l^2+i{k}_y(q_1-q_2)l^2},$$

$${\Phi }_{he}^{+-}(q_1,q_2,k_2,k_1;I_{2k_F})=l^2{e}^{-i(k_1{q}_1+k_2{q}_2)l^2}\int \frac{d{q}_{x}d{k}_y}{2\pi }$$

$$I_{2k_F}(q_x,k_y)e^{-i{q}_x{k}_y{l}^2+i{q}_x(k_1-k_2)l^2+i{k}_y(q_1-q_2)l^2}.$$

$$\int d(lQ)d(l{Q}^{\prime }){\Gamma }_{\mu \nu }(Q,Q_2,Q_1,Q^{\prime };f){\Gamma }_{\nu \lambda }(Q_2^{\prime },Q^{\prime },Q,Q_1^{\prime };g)$$

Full text

††preprint: APS/123-QED

I introduction

One of the fundamental challenge in condensed matter physics is a realization of three-dimensional unconventional electronic phases in the quantum limit halperin87 . Recent experimental discoveries of Dirac, Weyl and nodal Dirac semimetal materials weyl-dirac18 lead to growing research interests on novel quantum transports and quantum phase transitions under high magnetic field in these new compounds tang19 ; zhang19 ; fujioka19 ; zhang17 ; ramshaw18 as well as celebrated semimetal compounds such as bismuth zhu12 ; zhu17a and graphite fauque13 ; zhu19 ; arnold17 ; leboeuf17 . In fact, these semimetallic compounds under the high field often exhibit low-temperature metal-insulator transitions within wide ranges of the field tang19 ; fujioka19 ; yaguchi09 ; fauque13 ; akiba15 ; arnold17 ; leboeuf17 ; zhu17b ; taen18 ; zhu19 .

The quantizing effect of the strong magnetic field confines electrons into cyclotron motions in the Landau levels, while the electron’s kinetic energy along the field direction remains unaffected. This leaves the system with pristine one-dimensional momentum-energy dispersions along the field, making the system extremely sensitive to various instabilities gruner00 . Previous theories proposed a number of spontaneous symmetry broken (SSB) phases as well as non-Fermi liquid phaseyakovenko93 . The SSB phases proposed include charge-density wave halperin87 ; macdonald87 ; celli65 ; lee69 ; fukuyama78 ; yoshioka81 ; tesanovic87 ; zhang17 , three-dimensional quantum Hall halperin87 ; balents96 ; bernevig07 ; zhang17 , charge-Wigner crystal kleppmann75 ; fukuyama78 ; biagini01 ; tsai02a ; tsai02b ; alicea09 , spin-density wavetakahashi94 ; takada98 ; yaguchi09 ; pan18 , excitonic insulator fenton68 ; jerome67 ; pan18 ; song17 ; pan18 , valley-density wave, and three-dimensional topological excitonic insulator pan18 . Recent theoretical efforts on Dirac and Weyl semimetal models can be found in Ref. song17, ; trescher17, . In spite of these efforts during last decades, indentities of the low-temperature insulating phases in the experiments are still veiled in mystery due to a lack of comprehensive microscopic theory based on an unbiased theoretical method.

An electronic state of the prototypical semimetal materials can be captured by a pair of electron and hole band. Under the magnetic field ($\parallel$ $z$), the electron/hole’s motions in the $xy$ plane are confined into clockwise/anticlockwise cyclotron orbits around the field respectively. The counter-propagating cyclotron motion in the $xy$ plane inspires an ‘electron-hole’ analogy of the two-dimensional quantum spin Hall physics qi11 ; hasan10 . The one-dimensional dispersions along $z$ of the electron and hole bands go across the Fermi level at several Fermi points in the Brillouin zone. The experimental Hall conductivity measurements conclude that the relevant semimetal material within the relevant field regime yaguchi09 ; fauque13 ; akiba15 ; arnold17 ; leboeuf17 ; zhu17b ; taen18 ; zhu19 are in the charge neutrality region, where electron and hole densities compensate with each other completely uji98 ; kopelevich09 ; kumar10 ; akiba15 ; pan18 ; zhu19 . Thereby, to uncover the identities of the low-temperature insulator phases in the experiments yaguchi09 ; fauque13 ; akiba15 ; arnold17 ; leboeuf17 ; zhu17b ; taen18 ; zhu19 ; liang19 , it is vital to understand a ground-state phase diagram of a microscopic Hamiltonian for the semimetal material in the quantum limit at their charge neutrality point.

In this paper, using an unbiased theoretical method, we clarify the comprehensive ground-state phase diagrams of a prototype model for a three-dimensional semimetal under the magnetic field $H$ ($\parallel\!\ z$). The semimetal model has a pair of electron pocket with $\uparrow$ spin and hole pocket with $\downarrow$ spin. We study two limiting cases at the charge neutrality point, (i) the model with screened Coulomb interaction and (ii) the model with an effective attractive interaction mediated by the screened electron-phonon interaction. We rederive parquet renormalization group (RG) equations, that was originally given by Brazovskii  abrikosov70 ; brazovskii72 ; yakovenko93 , solve numerically the parquet RG equations, and complete the ground-state electronic phase diagrams for these two cases. The ground-state phase diagrams thus obtained exhibit rich phase diagram structures as a function of three important length scales in the system, ($\alpha$) magnetic length $l\equiv\sqrt{\hbar/eH}$ , ($\beta$) Fermi wave length (along the field) $2\pi/(2k_{F})$, and ($\gamma$) screening length of the long-range Coulomb interaction, $\sqrt{A}l$ [see Eq. (1) for the definition of a dimensionless parameter $A$, see Eq. (57) for its RPA evaluation]. In the repulsive interaction case, we found that the ground state is either Ising-type spin density wave phase or excitonic insulator (EI) phase (Fig. 1). In the attractive interaction case, we found that the ground state is either charge Wigner crystal phase, charge density wave phase, possible non-Fermi liquid yakovenko93 phase, EI phase that breaks the translational symmetry within the $xy$ plane or three-dimensional topological EI phase with a topological surface massless Dirac state pan18 (Fig. 2). We show in the paper that a canted magnetic field $H_{\perp}$ splits the surface Dirac state into surface Landau levels (sLL). The result suggests that when in-plane transport in the topological EI phase is dominated by surface transport, the in-plane resistivity must show a $\sqrt{H_{\perp}}$-type surface Shubnikov-de Haas (SdH) oscillation under the canted magnetic field.

I.1 highlight of the paper

In the next section, we introduce a model Hamiltonian for the semimetal material in the quantum limit, where a pair of the electron pocket with $\uparrow$ spin and hole pocket with $\downarrow$ spin go across the Fermi level at $p_{z}=\pm k_{F}$ at the charge neutrality point [$p_{z}$ is momentum along the field direction]. By using the RPA approximation, we discuss in Sec. III how the long-range Coulomb interaction $V({\bm{r}})\equiv e^{2}/{\varepsilon r}$ is screened by a density fluctuation at $p_{z}=0$ and a fast mode of the density fluctuation at $p_{z}=2k_{F}$. The screened Coulomb interactions thus obtained take the following forms in the momentum space tsvelik18 ,

[TABLE]

where $q_{\perp}$ is momentum within the $xy$ plane. Here two dimensionless parameters $1/A$ and $1/A^{\prime}$ are nothing but bare polarization functions associated with the density fluctuation at $p_{z}=0$ and the fast mode of the density fluctuation at $p_{z}=2k_{F}$ respectively. Eq. (1) especially indicates that the screening length of the long-range Coulomb interaction is given by the magnetic length times the dimensionless parameter $\sqrt{A}$; larger/smaller $A$ stands for the weak/strong screening respectively.

The polarization function associated with the slow density fluctuation at $p_{z}=2k_{F}$ as well as the polarization function for an excitonic fluctuation have logarithmic singularities abrikosov70 ; brazovskii72 . These singularities indicate several competing instabilities in the semimetal model at lower temperature. To identify the most dominant instability in the model within a controlled theoretical framework, we rederive in Sec. IV the one-loop parquet renormalization group (RG) equations, where the slow $2k_{F}$ density fluctuation as well as other fluctuations with the logarithmic singularity are taken into account on the same footing brazovskii72 ; zheleznyak97 . Using the screened Coulomb interaction [Eqs. (1,2)] as an ‘initial’ interaction form for the RG equations, we solve the parquet RG equations numerically and identify the most relevant fluctuations (instabilities) at lower temperature for different values of $A$, $A^{\prime}$ and $2k_{F}l$ [an overall factor of $V(q_{\perp},p_{z})$, $(4\pi e^{2}l^{2})/\varepsilon$, can be absorbed into a RG scale change at the one-loop RG equation; it does not change the phase diagram]. By combining mean-field arguments with the numerical RG solutions, we construct in Sec. V a comprehensive ground-state phase diagram in the presence of the repulsive Coulomb interaction (Fig. 1).

As the complimentary aspect for the semimetal in the quantum limit, we also study in Sec. VI an effect of electron-phonon interaction in the semimetal model under high magnetic field. Thereby, we employ a correspondence between an electron-phonon coupled system and a system with an electron-electron interaction, and adopt an effective attractive electron-electron interaction fetter03 ; mahan00 . The effective interaction is mediated by the screened Coulomb interaction between electron and (acoustic) phonon, and thereby it takes the following form in the momentum space,

[TABLE]

Here $Z$ and $c$ are an electron valence of positively charged nucleus ion and a sound velocity of the acoustic phonon respectively, $\rho_{0}$ and $M$ are the density of the charged nucleus ions and a mass of the ion respectively. To clarify possible instabilities in the semimetal model in the presence of the electron-phonon interaction, we use this effective attractive interaction [Eq. (3)] as an initial interaction form for the parquet RG equations, and solve the RG equations numerically for different values of $A$, $A^{\prime}$ and $2k_{F}l$. The numerical solutions in combination with mean-field arguments gives out a comprehensive ground-state phase diagram in the presence of the effective attractive interaction (Fig. 2).

The two ground-state phase diagrams thus obtained accommodate a rich variety of electronic phases as a function of $2k_{F}l$ and $A$. In the repulsive Coulomb interaction case (Sec. V), the ground state (GS) for strong screening regime [$A\leq 3$ for $2k_{F}l\simeq 1$] is an EI phase with a spatially even-parity excitonic pairing between electron and hole at the same momentum $p_{z}$. Since the pairing is between the electron with $\uparrow$ spin and the hole with $\downarrow$ spin and it is between the electron and hole at the same spatial location within the $xy$ plane, the excitonic pairing field leads to a long-ranged ferro-type order of a $XY$ component of the spin-1 moment. For the weak screening regime [$A\geq 3$ for $2k_{F}l\simeq 1$], the GS is a plain superposition of a density wave (DW) of the electron band with $\uparrow$ spin and a DW of the hole band with $\downarrow$ spin, that have the $2\pi/(2k_{F})$ spatial pitch along the field. Due to the Coulomb interaction between the two DWs, a relative phase between the two DWs is locked to $\pi$. Such superposition leads to an Ising-type spin density wave without any charge density modulation; the spatial pitch of the Ising-type antiferromagnetic order is $2\pi/(2k_{F})$.

In the attractive interaction case (Sec. VI), the GS for a strong screening regime [$A\leq 0.3$ for $2k_{F}l\simeq 1$] is either one of two distinct EI phases or the charge Wigner crystal. In one of the EI phases, the excitonic pairing is the spatially even-parity excitonic pairing between electron and hole at the same momentum $p_{z}$, but it is between the electron and hole at the different spatial location within the $xy$ plane [the field $\parallel\!\ z$]. As a result, the EI phase forms a two-dimensional texture of the $XY$ component of the spin-1 moment, breaking the translational symmetries within the $xy$ plane. The charge Wigner crystal phase breaks both the translational symmetries along the field and within the $xy$ plane by a three-dimensional texture of the charge density. The GS for the weak screening regime [$A\geq 0.3$ for $2k_{F}l\simeq 1$] is either a plain charge density wave phase with the $2\pi/(2k_{F})$ spatial pitch along the field [$2k_{F}l<1$ for $A\simeq 1$] or a possible non-Fermi liquid phase [$2k_{F}l>1$ for $A\simeq 1$].

The other EI phase found in the strong screening regime is a three-dimensional topological band insulator pan18 ; the EI phase supports a single copy of massless surface Dirac fermion state at its side surface [side surface is parallel to the field; $zx$ and $yz$ planes with the field along $z$]. The EI phase in the bulk is characterized by a spatially odd-parity pairing between electron and hole at the same momentum $p_{z}$ and at the same spatial location within the $xy$ plane. As a result, the EI phase does not break any translational symmetries. Besides, it has no local $XY$ component of the spin-1 moment, since the odd parity leads to a cancellation between $p_{z}$ and $-p_{z}$. Meanwhile, the odd-parity excitonic pairing in the bulk reconstructs a surface chiral Fermi arc state of the electron band with $\uparrow$ spin and that of the hole band with $\downarrow$ spin into the massless surface Dirac state with a helical spin texture. According to the so-called ‘periodic table’ of non-interacting topological insulator and topological crystalline insulator schnyder08 ; kitaev09 ; shiozaki14 ; kruthoff17 , the EI phase can be classified as topological ‘magnetic crystalline’ insulator shiozaki14 ; sato19 , where the massless nature of the surface Dirac fermion is protected by a magnetic point group symmetry $C_{2,\perp}T$ [$C_{2,\perp}$ denotes a $\pi$ rotation that changes $z$ to $-z$, and $T$ is the time reversal]. To give a physical characterization to the topological EI phase, we show in Sec. VII that a canted magnetic field $H_{\perp}$ splits the massless surface Dirac state into surface Landau levels (sLL), whose energy spacing is proportional to $\sqrt{H_{\perp}\Delta_{0}}$ [$\Delta_{0}$ is a strength of the excitonic pairing]. The result suggests that the longitudinal electric surface transport in the $xy$ direction can show a $\sqrt{H_{\perp}}$-type surface Shubnikov-de Haas (SdH) oscillation under the canted magnetic field. Based on these finding, we give a brief summary and discussion on the semimetal experiments in Sec. VIII.

II model Hamiltonian

An interacting electron model with a pair of electron pocket and hole pocket under high magnetic field $H$ ($\parallel z$) is considered;

[TABLE]

with ${\bm{\pi}}\equiv(\pi_{x},\pi_{y})$. $-E_{g}$ and $E_{g}$ are charge state energies of an electron-type and hole-type bands at the $\Gamma$ point. $m_{e}$ and $m_{h}$ are effective masses of the electron and hole band respectively, $m_{e}>0$ and $m_{h}>0$. $\tau=\pm$ refers to the spin $1/2$ degree of freedom ($+\equiv\!\ \uparrow$, $-\equiv\!\ \downarrow$). $H_{z}$ denotes the Zeeman field and we assume that the $g$ factor is isotropic in spin and same for electron and hole band.

An electron density $\hat{\rho}({\bm{r}})$ is a sum of the density of the electron band and that of the hole band, $\rho({\bm{r}})\equiv\sum_{a=e,h}\sum_{\tau}\Psi^{\dagger}_{a,\tau}({\bm{r}})\Psi_{a,\tau}({\bm{r}})$. In this paper, we consider as the electron correlation $V({\bm{r}})$ either repulsive Coulomb interaction (Secs. III,V) or an effective attractive electron-electron interaction mediated by the screened electron-phonon interaction (Sec. VI).

Due to the Landau quantization, the kinetic energy within a plane perpendicular to the field is quenched, where the electron band and hole band form a sequence of the Landau levels respectively;

[TABLE]

with the cyclotron frequency $\hbar\omega_{e,h}\equiv eH/m_{e,h}$. we consider the charge neutrality region in the quantum limit, where only the lowest Landau levels ($n=0$) with $\uparrow$-spin electron and $\downarrow$-spin hole bands go across the Fermi level at the same Fermi points ($p_{z}=\pm k_{F}$) , while all the others Landau levels leave the Fermi level (Fig. 3). For simplicity, we assume that the electron mass and hole mass are same, $m_{e}=m_{h}\equiv m$, $\omega_{e}=\omega_{h}\equiv\omega$.

We linearize the kinetic energy along the field direction around the two Fermi points. This leads to the following low-energy Hamiltonian for $\hat{h}_{0}$,

[TABLE]

where the magnetic length $l\equiv\sqrt{\hbar/eH}$ and $v_{F}=\frac{\hbar^{2}k_{F}}{m}$. $\sigma=\pm$ distinguishes two Fermi points, $k_{z}=\pm k_{F}$. $Q$ in Eq. (12) denotes momentum within the $xy$ plane. For the ‘right-mover’ fermion with positive velocity along the field, $e_{+}(Q,p)$ and $h_{-}(Q,p)$, we use a Landau gauge with an eigenstate localized along $x$-direction ($x$-gauge); $Q$ is momentum along $y$-direction, $Q=k_{y}$. For the ‘left-mover’ fermion with the negative velocity, $e_{-}(Q,p)$ and $h_{+}(Q,p)$, we use the Landau gauge with an eigenstate localized along $y$ ($y$-gauge); $Q$ is the momentum along $x$, $Q=q_{x}$ yakovenko93 . To be more specific, electron-band and hole-band creation operators in Eq. (5) are expanded as,

[TABLE]

with the eigenstates in the LLL,

[TABLE]

For later convenience, note that the $x$-gauge eigenstates and $y$-gauge eigenstates are transformed to each other by a Fourier transformation;

[TABLE]

By substituting Eqs. (13,14) into the electron-electron interaction in Eq. (4), we obtain,

[TABLE]

[TABLE]

with a notation of

[TABLE]

$H_{1}$ is a sum of all the interactions that carry the zero momentum along the field (Fig. 4(a)), while $H_{2}$ is a sum of all the interactions that carry the $2k_{F}$ momentum along the field (Fig. 4(c)). The respective interaction potentials are given by functionals of following two ‘bare’ functions of the in-plane momentum $(q_{x},k_{y})$,

[TABLE]

with $V(q_{x},k_{y},p_{z})\equiv\frac{4\pi e^{2}}{\varepsilon(q^{2}_{x}+k^{2}_{y}+p^{2}_{z})}$. Specifically, $\Gamma_{\mu\nu}$ for $(\mu,\nu)=(e_{+},e_{+})$, $(h_{-},h_{-})$, $(e_{+},h_{-})$, $(h_{-},e_{+})$, $\Gamma_{\mu\nu}$ for $(\mu,\nu)=(e_{-},e_{-})$, $(h_{+},h_{+})$, $(e_{-},h_{+})$, $(h_{+},e_{-})$, $\Gamma_{\mu\nu}$ for $(\mu,\nu)=(e_{+},e_{-})$, $(h_{-},h_{+})$, $(e_{+},h_{+})$, $(h_{-},e_{-})$, and $\Gamma_{\mu\nu}$ for $(\mu,\nu)=(e_{-},e_{+})$, $(h_{+},h_{-})$, $(h_{+},e_{+})$, $(e_{-},h_{-})$ are given by the same functionals of $I_{0}(q_{x},k_{y})$ respectively;

[TABLE]

with

[TABLE]

$\Phi^{+-}_{\mu\nu}$ for $(\mu,\nu)=(e,e),(h,h)$, $\Phi^{+-}_{\mu\nu}$ for $(\mu,\nu)=(e,h)$, and $\Phi^{+-}_{\mu\nu}$ for $(\mu,\nu)=(h,e)$ are given by the following functionals of $I_{2k_{F}}(q_{x},k_{y})$,

[TABLE]

For later convenience, note that these functionals can be regarded homomorphic mappings of the functions. Namely, a product between two functionals of functions $f$ and $g$ is a functional of $fg$,

[TABLE]

for any $\mu,\nu,\lambda=e_{+},e_{-},h_{+},h_{-}$. Here a summation over $\nu$ is not taken in their left-hand sides. Similarly,

[TABLE]

for any $\mu,\nu,\lambda=e,h$. These homomorphic natures of the functionals are useful in the next section.

In the following, we consider Eqs. (12,19,20,21) as the prototype model Hamiltonian for a semimetal in the quantum limit at the charge neutrality point.

III screened Coulomb interaction

The interaction potentials in $H_{1}$ are screened by low-energy density fluctuations at $p_{z}=0$, while the interaction potentials in $H_{2}$ are screened by the $2k_{F}$ density fluctuations. The respective screened interaction comprises of a sum of the bare interaction part and an effective interaction mediated by the density fluctuations. Using the random phase approximation (Fig. 4(b,d)) with a help of the homomorphic nature of the interaction potentials, Eqs. (41,42), we can show that the screened forms for the interaction potentials $\Gamma_{\mu\nu}$ and $\Phi^{+-}_{\mu\nu}$ take exactly the same forms as their respectively bare forms in Eqs. (20,21), except that their arguments, $I_{0}(q_{x},k_{y})$ and $I_{2k_{F}}(q_{x},k_{y})$, are replaced by their screened counterparts, $\overline{I}_{0}(q_{x},k_{y})$ and $\overline{I}_{2k_{F}}(q_{x},k_{y})$, respectively,

[TABLE]

Here $\Pi_{0,\lambda}(\omega)$ denotes a bare polarization function for the $p_{z}=0$ density fluctuation in the right/left-mover electron band ($\lambda=e_{+/-}$) or left/right-mover hole band ($\lambda=h_{+/-}$). It is given by

[TABLE]

with the Heaviside step function $\theta(x)$. From Eq. (12), $\varepsilon_{\lambda,p}-\varepsilon_{F}\equiv v_{F}p$ for $\lambda=e_{+},h_{-}$ and $\varepsilon_{\lambda,p}-\varepsilon_{F}\equiv-v_{F}p$ for $\lambda=e_{-},h_{+}$. Noting that $p$ is much smaller than $k_{F}$, one obtains the bare polarization function as

[TABLE]

for $\lambda=e_{+},e_{-},h_{+},h_{-}$. $\Pi^{-+}_{0,\lambda}(\omega)$ is the bare polarization function for the $2k_{F}$ density fluctuation within the electron band ($\lambda=e$) or hole band ($\lambda=h$),

[TABLE]

From Eq. (12), $\varepsilon_{\lambda_{\pm},p}-\varepsilon_{F}=\pm v_{F}p$ for $\lambda=e$ and $\varepsilon_{\lambda_{\pm},p}-\varepsilon_{F}=\mp v_{F}p$ for $\lambda=h$.

The polarization function for $2k_{F}$ density fluctuation can be divided into low-energy (slow) $2k_{F}$ density fluctuation part and high-energy (fast) $2k_{F}$ density fluctuation part,

[TABLE]

Here $\Lambda_{\rm UV}$ is a momentum cut off associated with the Brillouin zone boundary, while $\Lambda$ separates the slow mode from the fast mode, $\Lambda<\Lambda_{\rm UV}$. The slow $2k_{F}$ density fluctuation part ($|p_{1}|<\Lambda$) leads to the logarithmic singularity at $p=0$ zheleznyak97 ,

[TABLE]

for any $\lambda=e,h$. The singularity is a seed of the $2k_{F}$ DW instability in each pocket gruner00 ; abrikosov70 ; brazovskii72 ; yakovenko93 . More generally, ‘bubble’ Feynman diagrams composed of two single-particle Green functions with opposite sign of the Fermi velocities have the same logarithmic singularity both in particle-hole and in particle-particle channels abrikosov70 ; brazovskii72 . These logarithmic singularities suggest that the ground state of the electron-hole model at the charge neutrality point has several competing instabilities at lower temperature. To clarify the most dominant instability precisely, we thus take into account the slow $2k_{F}$ density fluctuation in the framework of parquet RG equations and include it as well as the other low-energy fluctuations with the logarithmic singularity on the equal footing (see the next section). Therefore, to avoid the double counting of the slow $2k_{F}$ density fluctuation part, we include in Eq. (44) only the high-energy (fast) $2k_{F}$ density fluctuation part first. This determines a form of the screened interaction. The screened interaction thus obtained is then used for an ‘initial’ interaction potential for the RG equations (Fig. 5). Finally, the low-energy (slow) $2k_{F}$ density fluctuation as well as other dominant low-energy fluctuations shall be included sequentially in the framework of the RG procedure (Fig. 6). For example, in Fig. 6, the last two Feynman diagrams in the right-hand side of the first and second lines as well as the first two Feynman diagrams in the right-hand side of the fourth line represents the inclusions of the slow $2k_{F}$ density fluctuations. The fast $2k_{F}$ density fluctuation part, Eq. (50), takes a constant finite value at $p=0$,

[TABLE]

for $\lambda=e,h$.

To summarize, we will use the following form of the screened interaction potentials as the initial interaction forms for the later parquet RG studies;

[TABLE]

with $B\equiv 4k^{2}_{F}l^{2}$ and

[TABLE]

For simplicity, we take $\ln(\Lambda_{\rm UV}/\Lambda)=2$ henceforth and identify $A^{\prime}$ with $A$. For reminder, the functionals $\Gamma_{\mu\nu}$ ($\mu,\nu=e_{+},e_{-},h_{+},h_{-}$) and $\Phi^{+-}_{\mu\nu}$ ($\mu,\nu=e,h$) in Eqs. (53,54) are defined in Eqs. (29,30,31,32,38, 39,40) respectively.

IV parquet renormalization group equation

The polarization function for the slow $2k_{F}$ density fluctuation has the logarithmic singularity [Eq. (51)]. More generally, all the ‘bubble’ diagrams composed of the two Green functions with opposite sign of the Fermi velocities have the same logarithmic singularity in both particle-hole and particle-particle channels. The presence of the logarithmic singularities in several distinct channels means competing ground-state instabilities in the semimetal model. To reveal the ground-state phase diagram of the model precisely, we thus include all the relevant fluctuations with the logarithmic singularity on the equal footing.

To this end, we derive in this section the parquet renormalization group (RG) equations abrikosov70 ; brazovskii72 ; yakovenko93 , where consecutive integration of the higher-energy fermionic degree of freedom renormalizes the interaction potentials among the lower-energy fermions. The renormalization gives rise to either enhancement, suppression or convergence of the interaction potentials. By identifying the most divergent potentials among the others, we shall tell the dominant ground-state instability in the model.

The one-loop parquet RG equations can be derived by a standard momentum shell renormalization. Thereby, we begin with a partition function of the interacting fermion model;

[TABLE]

[TABLE]

Here $Q=k_{y}$ or $q_{x}$ is rescaled by the magnetic length $l$; $Q_{\rm new}\equiv Q_{\rm old}l$. Besides, we used the following notations,

[TABLE]

For the repulsive Coulomb interaction case, the interaction potentials in $S_{1}$ are given by either some of $\overline{H}_{1}$ or $\overline{H}_{2}$ or their combination from Eqs. (53,54,55,56) (see Fig. 5). Namely, $W_{b}({\bm{k}})$ and $W_{d}({\bm{k}})$ are given by a sum of Eq. (31) and Eq. (38) with $I_{0}$ and $I_{2k_{F}}$ replaced by $\overline{I}_{0}$ and $\overline{I}_{2k_{F}}$. $W_{e}({\bm{k}})$ and $W_{g}({\bm{k}})$ are given by Eq. (31) and by Eq. (40) respectively with $\overline{I}_{0}$ and $\overline{I}_{2k_{F}}$. To be more specific, we consider the following set of the screened interaction as the initial interaction forms of the RG equations,

[TABLE]

with ${\bm{k}}\equiv(k_{y},q_{x})$ and Eqs. (55,56). $``\cdots"$ in Eq. (60) denotes those interaction parts that are not renormalized by others and do not renormalize others at the level of the one-loop RG equations. Such interaction parts are irrelevant in the framework of the one-loop RG analyses; we thus omit them henceforth.

Following the standard momentum shell renormalization process (Appendix B), we first decompose the fermionic field into fast mode ($\Lambda-d\Lambda<|p|<\Lambda$) and slow mode ($|p|<\Lambda-d\Lambda$) in the momentum along the field. The integration of the fast mode in the partition function leads to the renormalizations of the interaction potentials among the slow modes. This gives out a set of coupled RG equations for the interaction potentials,

[TABLE]

with ${\bm{k}}\equiv(k,q)$, ${\bm{k}}^{\prime}\equiv(k^{\prime},q^{\prime})$, and $d{\bm{k}}^{\prime}\equiv dk^{\prime}dq^{\prime}$. $d\xi$ denotes a RG scale,

[TABLE]

In order to solve the coupled RG equations numerically, we put them in the dual-space representation by the Fourier transform of $W_{\mu}({\bm{k}})$,

[TABLE]

for $\mu=b,d,e,g$ with

[TABLE]

${\bm{r}}\equiv(r_{x},r_{y})$ and ${\bm{k}}\equiv(k,q)$. It turns out that parquet RG equations for $F_{\mu}({\bm{r}})$ ($\mu=b,d,e$) and $\tilde{F}_{g}({\bm{r}})$ as well as their initial forms respect the following O(2) symmetry and real-valued-ness;

[TABLE]

with $r\equiv|{\bm{r}}|$ for arbitrary $\theta\in(0,2\pi]$. Using this symmetry, we can finally reach O(2)-symmetric parquet RG equations for $\Gamma_{\mu}(r)$ ($\mu=b,d,e$) and $\Gamma_{g}(r)$ as follows brazovskii72 ,

[TABLE]

with

[TABLE]

and Bessel function $J_{2m}(r)$ (integer $m$). $\overline{\Gamma}_{g}(r)$ is a Hankel transform of $\Gamma_{g}(r)$;

[TABLE]

The initial forms of $\Gamma_{\mu}(r)$ ($\mu=b,d,e,g$) and $\overline{\Gamma}_{g}(r)$ are obtained from Eqs. (61,62,63) as follows,

[TABLE]

with

[TABLE]

Note that the overall factor of $\overline{I}_{0}(r)$ and $\overline{I}_{2k_{F}}(r)$, $4\pi e^{2}l^{2}/\varepsilon$, can be absorbed into a redefinition of the RG scale change $\xi$; it does not alter the ground-state phase diagram. Only two dimensionless parameters, $A$ and $B\equiv 4k^{2}_{F}l^{2}$ in Eqs. (88,89), play vital role in a determination of the ground-state phase diagram.

We solved Eqs. (77,78,79,80) numerically, with $\Gamma_{\mu}(r)$ ($\mu=b,d,e,g$) and $\overline{\Gamma}_{g}(r)$ at the initial RG scale $(\xi=0)$ being given by Eqs. (84,85,86, 87,88,89). By doing so, we numerically observed that in the two-dimensional $A$-$B$ space, either a set of $\Gamma_{b}(r,\xi)=\Gamma_{d}(r,\xi)$ and $\Gamma_{g}(r,\xi)$ or a set of $\Gamma_{e}(r,\xi)$ and $\overline{\Gamma}_{g}(r,\xi)$ show divergences at certain values of $r$ and $\xi$ (see the next section). The divergence indicates a certain type of pairing instabilities in the ground state. To identify the favored pairings and natures of resulting symmetry-broken phases, we rewrite the interaction potentials in Eq. (60) in the same basis of the Landau gauge. For example, we put $e_{-}$ and $h_{+}$ as well as $e_{+}$ and $h_{-}$ in the basis of the $x$-gauge eigenstates by using Eqs. (17,18). This leads to

[TABLE]

where $k_{1}$, $\overline{k}_{1}$, $k_{2}$ and $\overline{k}_{2}$ are the momentum along the $y$-direction. Note that due to the translational symmetry along the $y$ direction in the $x$-gauge, all the interaction potentials preserve a center of mass in the momentum. From Eqs. (70,72,73), one can readily see that $\Phi_{b}(k,k^{\prime})$, $\Phi_{e}(k,k^{\prime})$ and $\Phi_{g}(k,k^{\prime})$ in Eq. (90) are given by $\Gamma_{b}(r)=\Gamma_{d}(r)$, $\Gamma_{e}(r)$, $\Gamma_{g}(r)$ and $\overline{\Gamma}_{g}(r)$ as follows,

[TABLE]

One could also rewrite $e_{+}$ and $h_{-}$ as well as $e_{-}$ and $h_{+}$ in the basis of the $y$-gauge eigenstates. Of course, this leads to the same conclusions as we will reach in the $x$-gauge eigenstates (see the following two sections).

V Ground-state phase diagram in the presence of repulsive Coulomb interaction

The parquet RG equations have a dual structure; $\Gamma_{b}(r)=\Gamma_{d}(r)$ and $\Gamma_{g}(r)$ couple with each other exactly in the same way as $\Gamma_{e}(r)$ and $\overline{\Gamma}_{g}(r)$ do, and $\Gamma_{g}(r)$ and $\overline{\Gamma}_{g}(r)$ are Fourier transforms of the other [Eqs. (82,83)]. In the case of the repulsive interaction, this dual structure in the RG equations leads to a ground-state competition between the excitonic insulator phase brazovskii72 and Ising-type spin density wave phase (Fig. 1). The numerical solution of the RG equations shows that in the two-dimensional $A$-$B$ parameter space, either a set of $\Gamma_{e}(r)$ and $\overline{\Gamma}_{g}(r)$ or a set of $\Gamma_{b}(r)$ and $\Gamma_{g}(r)$ diverge at a certain critical RG scale, $\xi=\xi_{c}$.

V.1 strong screening region

To understand the phase diagram qualitatively, let us keep only those terms in the parquet RG equations that couple the functions locally in the radial coordinate $r$; brazovskii72 ; yakovenko93 ; tsai02a ; alicea09

[TABLE]

and

[TABLE]

When the RG scale is near (but below) the critical RG scale, $\xi\lesssim\xi_{c}$, the local terms become leading order than those terms neglected, and the approximation and solutions below are justified. Without the constraint between $\Gamma_{g}(r)$ and $\overline{\Gamma}_{g}(r)$ (Eqs. (82),(83)), the approximate RG equations can be solved and the solutions are determined only by the initial forms of the interaction potentials,

[TABLE]

and

[TABLE]

With Eqs. (84,85,86,87,88,89) at the initial RG scale ($\xi=0$), $\Gamma_{e}(r,0)+\overline{\Gamma}_{g}(r,0)$ takes the largest positive value at $r=0$ among the other three at any $r$,

[TABLE]

Thus, the approximate solution dictates that positive $\Gamma_{e}(r,\xi)$ and positive $\overline{\Gamma}_{g}(r,\xi)$ diverge at $r=0$ simultaneously on the renormalization as brazovskii72 ,

[TABLE]

Fig. 7 demonstrates how the four interaction potentials, $\Gamma_{b}(r,\xi)=\Gamma_{d}(r,\xi)$, $\Gamma_{e}(r,\xi)$, $\Gamma_{g}(r,\xi)$ and $\overline{\Gamma}_{g}(r,\xi)$, change their forms under the parquet RG equations, Eqs.(77,78,79,80), in the strong screening regime.

When $\Gamma_{e}(r=0)$ and $\overline{\Gamma}_{g}(r=0)$ dominate over the others, an excitonic pairing is formed between electron and hole bands at the same Fermi point and at the same spatial coordinate within the $xy$ plane,

[TABLE]

Namely, the asymptotic forms of $\Gamma_{e}(r,\xi_{c})$ and $\overline{\Gamma}_{g}(r,\xi_{c})$ make the following scattering channels in Eq. (90) to be dominant among the others,

[TABLE]

for arbitrary $\overline{k}_{2}-\overline{k}_{1}$ (Eq. (106)) and arbitrary $k_{2}-\overline{k}_{2}$ (Eq. (107)) respectively. These scatterings favor electron-hole pairings at the same Fermi points and at the same two-dimensional space coordinates within the $xy$ plane;

[TABLE]

Note that the relative U(1) phase between the excitonic pairing field at the right Fermi point and that at the left Fermi point is locked to be zero by the positively large $\Phi_{g}(k_{1}-\overline{k}_{2}=0,k_{2}-\overline{k}_{2})$.

The excitonic pairing between the electron band with $\uparrow$ spin and the hole band with $\downarrow$ spin results in a ferro-type order of an $XY$ component of the spin-1 moment;

[TABLE]

The ferro-type order breaks the U(1) spin rotational symmetry around the magnetic field. However, detailed microscopic magnetism of the excitonic phase depends on atomic orbitals (localized Wannier orbitals) that form the electron band and the hole band.

V.2 weak screening region

When the screening length is longer than the magnetic length [$A\geq 3$ for $2k_{F}l\simeq 1$], the numerical solution shows that $\Gamma_{b}(r)$ and $\Gamma_{g}(r)$ diverges at $r=0$ as;

[TABLE]

The divergence identifies the relevant scattering channels in Eq. (90) as,

[TABLE]

for any $\overline{k}_{2}-\overline{k}_{1}$, and

[TABLE]

for any $k_{1}-\overline{k}_{2}$. These scattering channels cause an instability to a charge density wave of the electron band and that of the hole band,

[TABLE]

Both density waves share the same spatial pitch ($2\pi/2k_{F}$) along the field direction. The relative U(1) phase between the electron-band density wave and hole-hand density wave is locked to be $\pi$ by the positively large $\Phi_{g}(k_{1}-\overline{k}_{2},k_{2}-\overline{k}_{2}=0)$. Due to the $\pi$ phase shift, the ground state in the weak screening region is accompanied by Ising-type spin density wave that preserves the U(1) spin rotational symmetry around the magnetic field (Fig. 8).

VI Ground-state phase diagram in the presence of effective attractive interaction

In the previous section, we have studied how the repulsive Coulomb interaction leads to the low-temperature instability in the semimetal under high magnetic field. As the complimentary aspect, we consider in this section an effect of another relevant many-body interaction; electron-(acoustic) phonon interaction. We employ an argument based on an equivalence between an electron-phonon coupled system and a system with an electron-electron interaction, and adopt the following effective attractive electron-electron interaction;

[TABLE]

Here $U_{\rm eff}({\bm{q}})$ is the Fourier transform of a screened Coulomb potential between electron and (longitudinal acoustic) phonon. $Z$ and $c$ are an electron valence of positively charged nucleus ion and a sound velocity of the acoustic phonon respectively, $\rho_{0}$ and $M$ is the density of the charged nuclei, and a mass of the charged nucleus ion. Within the random phase approximation, (the square of) the screening length ‘$A$’ in $U_{\rm eff}({\bm{q}})$ was calculated in the previous section [Eq. (57,58)]. Here we consider the case with $A=A^{\prime}$ for simplicity.

Using Eqs. (111,112,113,114) as the effective electron-electron interaction, we study low-temperature instabilities in semimetal under high magnetic field in the presence of the electron-phonon coupling. To this end, we solve numerically the same parquet RG equations as in the previous section, while we use the following set as the interaction forms at the initial RG scale ($\xi=0$);

[TABLE]

with

[TABLE]

Again, the overall factor of $\overline{I}^{\prime}_{0}(r)$ and $\overline{I}^{\prime}_{2k_{F}}(r)$ does not play any role in a determination of the phase diagram within the one-loop RG analyses. Only the two dimensionless parameters $A$ and $B$ play the crucial role.

VI.1 intermediate screening region

Fig. 2 is a phase diagram obtained by the numerical solutions. In an intermediate screening region ($A\simeq 10^{-1}$), the ground state shows an instability toward a charge Wigner crystal phase, where $\Gamma_{b}(r)$ and $\Gamma_{g}(r)$ diverges at nonzero $r$ ($r=r_{c}\neq 0$) at a certain critical RG scale ($\xi=\xi_{c}$) as,

[TABLE]

Substituting this into Eqs. (91,93), one can see that dominant scattering channels in Eq. (90) take the following asymptotic forms,

[TABLE]

for any $\theta\in[0,\pi)$, and for any $\overline{k}_{2}-\overline{k}_{1}$ (Eq. (122)) and any $k_{1}-\overline{k}_{2}$ (Eq. (123)), respectively. The scattering channels induce $2k_{F}$ density-wave pairings within the electron band and hole hand. The induced density-wave pairings generally connect different two-dimensional coordinates within the $xy$ plane,

[TABLE]

Due to the coordinate-dependent ($k$-dependent) cosine functions in Eqs. (122,123), the action $S_{1}(\theta)$ is fully minimized by the pairing fields that have coordinate-dependent phases,

[TABLE]

Such pairings lead to the density waves in the electron and hole bands, that break the translational symmetries within the $xy$ plane,

[TABLE]

The density wave of the electron band with $\uparrow$ spin and the density wave of the hole band with the $\downarrow$ spin have the same phase; the superpose of these two is nothing but the charge density wave without any spin texture. The spatial pitches within the $xy$ plane and along the field direction is $2\pi l/r_{c}$ and $2\pi/(2k_{F})$ respectively.

The ‘propagation’ direction of the density wave within the $xy$ plane is specified by $\theta$, that can take any value in $[0,\pi)$ according to Eqs. (122,123,124). The ground state is generally a superposition of the density waves with different propagation directions within the $xy$ plane. One of the most plausible superposition is a symmetric superposition,

[TABLE]

with ${\bm{r}}_{\perp}=(x,y)$, ${\bm{n}}_{1}=(1,0)$, ${\bm{n}}_{2}=(-\frac{1}{2},\frac{\sqrt{3}}{2})$, and ${\bm{n}}_{3}=(-\frac{1}{2},-\frac{\sqrt{3}}{2})$ [or its O(2) rotation within the $xy$ plane]. This leads to a triangle lattice of the charge density within the $xy$ plane (charge Wigner crystal; Fig. 9(a)).

VI.2 strong screening region

In a strong screening region ($A<10^{-2}$), the ground state exhibits an instability to an excitonic phase, where $\Gamma_{e}(r)$ and $\overline{\Gamma}_{g}(r)$ diverge at nonzero $r$ ($r=r_{c}\neq 0$) at the critical RG scale ($\xi=\xi_{c}$) as;

[TABLE]

The divergence gives rise to the following forms of the dominant scattering channels in Eq. (90),

[TABLE]

for any $\theta$, and for any $\overline{k}_{2}-\overline{k}_{1}$ (Eq. (128)) and any $k_{2}-\overline{k}_{2}$ (Eq. (129)) respectively. These scattering channels mediate the excitonic pairings between different spatial coordinate within the $xy$ plane;

[TABLE]

Namely, the action with Eqs. (128,129) is minimized by the excitonic pairing within the same Fermi points but between different spatial coordinates within the $xy$ plane. The pairing fields thus determined have the coordinate-dependent phase factors,

[TABLE]

Such excitonic pairings leads to a density wave of the $XY$ component of the spin-1 moment, that breaks the translational symmetry within the $xy$ plane;

[TABLE]

The propagation direction of the $XY$ spin density wave is characterized by the arbitrary phase $\theta$; the ground state takes a form of the superposition of the density waves over different propagation directions within the $xy$ plane. From an analogy of the charge Wigner crystal phase, one of the possible spatial structures of the $XY$ spin moment is the symmetric superposition,

[TABLE]

with ${\bm{n}}_{1}=(1,0)$, ${\bm{n}}_{2}=(-\frac{1}{2},\frac{\sqrt{3}}{2})$, and ${\bm{n}}_{3}=(-\frac{1}{2},-\frac{\sqrt{3}}{2})$. This results in a two-dimensional vortex lattice, where vortices of the $XY$ spin with $\pm 1$ charges form a two-dimensional honeycomb structure (Fig. 9(b,c)).

In actual semimetal compounds, the emergent two-dimensional lattice structure of the $XY$ spin moment as well as the charge density wave must be extremely sensitive to actual crystal symmetry of underlying lattice structure in each compound. Especially, $k$-$p$ expansion around a high symmetric momentum line (parallel to the magnetic field $\parallel z$) often gives rise to an anisotropy in the effective mass or effective velocity within a plane perpendicular to the field ($xy$ plane). The anisotropy reduces the in-plane O(2) symmetry in the model dictated by Eqs. (4,5,8) down to a discrete rotational symmetry around the field. For example, in the case of the graphite, the relevant electron and hole pockets around the zone boundary lines ($H$-$K$-$H$ and $H^{\prime}$-$K^{\prime}$-$H^{\prime}$) respect a Z3 discrete rotational symmetry, reflecting the graphite crystal structure. Speaking symmetry, the triangle lattice structure of the charge density (Fig. 9(a)) as well as the two-dimensional vortex lattice structure of the $XY$ spin moment (Fig. 9c) is compatible with this Z3 discrete rotational symmetry.

VI.3 weak screening region

In the weak screening region ($A\geq 1$), the phase diagram is covered by either charge density wave (smaller $k_{F}l$ region) or possible non-Fermi liquid (larger $k_{F}l$ region). In the charge density wave phase, $\Gamma_{b}(r)$ and $\Gamma_{g}(r)$ diverges at $r=0$ at a certain critical RG scale as,

[TABLE]

Equivalently, the effective pontentials in Eq. (90) will be dominated by the following scattering channels,

[TABLE]

for any $\overline{k}_{2}-\overline{k}_{1}$ and any $k_{1}-\overline{k}_{2}$ respectively. As in the previous section, the scatterings give rise to the $2k_{F}$ density wave of the electron band and that of the hole band. The relative U(1) phase between the two density waves is locked to be zero by the negatively large $\Phi_{g}(k_{1}-\overline{k}_{2},k_{2}-\overline{k}_{2}=0)$,

[TABLE]

The resulting ground state has a simple charge density modulation along the field direction (without any spin texture), whose spatial pitch is $2\pi/(2k_{F})$.

When the spatial pitch of the charge density modulation becomes shorter than the magnetic length ($1/(2k_{F}l)\lesssim 1$), the density wave undergoes a phase transition, and the ground state becomes a critical phase. In the critical phase, $\Gamma_{g}(r)$ as well as $\overline{\Gamma}_{g}(r)$ get renormalized to the zero at any $r$. Since $\Gamma_{g}(r)\equiv\overline{\Gamma}_{g}(r)\equiv 0$, the coupled parquet RG equations are decoupled into two RG equations,

[TABLE]

for $\mu=e,b$. Being attractive, both $\Gamma_{b}(r)$ and $\Gamma_{e}(r)$ converge to universal functions of $r$. The universal functions are solutions of the decoupled RG equation at larger RG scale, $\xi\gg\xi_{1}$ where $\xi_{1}$ is a certain short-range cutoff of the RG scale. The functions have a ‘self-similar’ structure (Fig. 10) yakovenko93 ,

[TABLE]

Yakovenko previously discovered this critical phase in single band model under the magnetic field and characterized this critical phase as marginal Fermi liquid phase, where the renormalization factor of the electron Green function vanishes in the large $\xi$ limit yakovenko93 .

VI.4 topological excitonic insulator

The numerical solutions also found a three-dimensional topological excitonic insulator phase between the charge Wigner crystal phase and the excitonic insulator phase with the $XY$-spin vortex lattice. Thereby, $\Gamma_{e}(r)$ and $\overline{\Gamma}_{g}(r)$ diverge at $r=0$ as

[TABLE]

The divergence chooses Eqs. (106,107) as the dominant scattering channels in Eq. (90), while the sign of $\Phi_{g}(k_{1}-\overline{k}_{2}=0,k_{2}-\overline{k}_{2})$ is negative. Such scattering channels lead to a formation of a ‘(spatially) odd-parity’ excitonic pairing that connects the same spatial coordinate within the $xy$ plane;

[TABLE]

Due to the opposite sign between the two pairings at the right and left Fermi points, the $XY$ components that come from these two Fermi points cancel each other. The phase has no local $XY$ component of the spin-1 moment.

As shown by the author previously, the excitonic insulator phase can be regarded as a topological band insulator that has a single copy of $(2+1)$D massless Dirac surface fermion at its side surface [side surface is along the field direction; $zx$ plane or $yz$ plane] pan18 . The emergence of the surface state results from the odd-parity excitonic pairing in the bulk and is a direct consequence of a ${\bm{Z}_{2}}$ topological integer defined in a bulk mean-field electronic Hamiltonian.

To explain this, note first that the bulk mean-field Hamiltonian takes a form of a sum of ‘one-dimensional’ Hamiltonian, as the excitonic pairing connects the same two-dimensional spatial coordinate within the $xy$ plane, $H_{\rm mf}\equiv\int dk_{y}H_{\rm 1D}(k_{y})$ with;

[TABLE]

[TABLE]

and $\mu_{0}\equiv E_{g}+H_{z}-\frac{1}{2}\hbar\omega$. Here we went back to Eq. (11) and wrote down explicitly the whole $k_{z}$-dependence of the kinetic energy along the field. Besides, we put a confining potential $V_{c}(x)$ in the $x$ gauge (Landau gauge) with $x=k_{y}l^{2}$. $V_{c}(x)$ is zero in the bulk region ($|x|<L/2$) and it becomes positively large in the vacuum region ($|x|>L/2$). $\Delta(k_{z})$ is the excitonic pairing potential. From Eq. (108) and Eq. (139), the potential is an odd function of the momentum along the field,

[TABLE]

When the U(1) phase in Eq. (139) is spatially uniform, one can absorb the phase into a relative gauge between the electron and hole bands. For a fixed $\theta$, one defines a winding number for the one-dimensional mean-field Hamiltonian heeger88 ; wen89 ,

[TABLE]

with

[TABLE]

Note that in the bulk region ($V_{c}(k_{y}l^{2})=0$), $N_{3}(k_{z},k_{y})$ is negative for $|k_{z}|<k_{F}$ and positive otherwise. Thus, the winding number must be an odd integer ($\pm 1$), since $N_{1}(k_{z},k_{y})$ is an odd function in $k_{z}$. Meanwhile, in the vacuum region ($V_{c}(k_{y}l^{2})=+\infty$), the confining potential becomes positively large, so that the winding number is always zero; the electron band and the hole band are ‘re-inverted’ and $N_{3}$ is always positive for all $k_{z}$.

The odd-even difference in the winding number in the one-dimensional mean-field Hamiltonians causes the emergence of the surface state in a boundary region between bulk and vacuum. Namely, by regarding $k_{y}l^{2}$ as a ‘parameter’ of one-dimensional electronic system, one can expect that the one-dimensional topological band insulator ‘phase’ with the odd integer winding ($|k_{y}l^{2}|<L/2$) and one-dimensional trivial band insulator ‘phase’ with the zero winding ($|k_{y}l^{2}|>L/2$) must be intervened by a one-dimensional topological ‘critical point’, that comes at the boundaries ($|k_{y}l^{2}|\simeq L/2$). In fact, since the excitonic pairing is spatially odd, the critical point is generally described by the (1+1)D massless Dirac fermion with a linear dispersion along the momentum $k_{z}$ at $k_{z}=0$. Besides, $M(k_{z}=0,k_{y})$ changes its sign at $|k_{y}l^{2}|\simeq L/2$. Thus, the mean-field Hamiltonian forms the $(2+1)$D massless Dirac Hamiltonian in the $k_{z}$-$k_{y}$ plane around $k_{z}=0$ and $|k_{y}l^{2}|\simeq L/2$;

[TABLE]

Note that the massless surface Dirac fermion has helical velocities in any directions within the side surface ($yz$ plane in the $x$-gauge). It has a helical velocity not only along the field direction ($\parallel z$) but also along the perpendicular direction ($\parallel y$). From these observations, the excitonic insulator phase with the odd-parity excitonic pairing can be regarded as a three-dimensional topological band insulator in the quantum limit.

In the next section, we will describe the effect of the tilted magnetic field on the topological surface state on the side surface.

VII effect of titled magnetic field on topological surface states

When the magnetic field is tilted from the $z$-axis to $Z$-axis with $Z\equiv-x\sin\theta+z\cos\theta$ ($0\geq\theta<\pi$), the excitonic pairing in the bulk remains intact; the three-dimensional semimetal model has a spatially isotropic effective mass [Eqs. (4,5)]. Meanwhile, the $(2+1)$D massless surface Dirac fermion on the side surface ($yz$ plane) forms Landau levels due to a finite out-of-surface component of the magnetic field. Equivalently, we can consider the same situation by tilting the ‘side’ surface from the $yz$ plane to the $yZ$ plane, and keep the field along the $z$ axis (Fig. 11).

Specifically, we add in Eqs. (5,8) a confining potential $V_{c}(X)$ that depends only on $X\equiv x\cos\theta+z\sin\theta$. For simplicity, we take the system is translationally symmetric along the $y$-direction, so that $-i\hbar\partial_{y}$ in Eqs. (5,8) is replaced by $\hbar k_{y}$. This gives out

[TABLE]

Here we assume that the odd-parity excitonic pairing is linear in $k_{z}$, $\Delta(k_{z})\equiv\Delta_{0}k_{z}$. The confining potential $V_{c}(X)$ takes a constant value, $V_{c}(X)=-E_{g}-H_{z}$, for those $X$ in the bulk region. $V_{c}(X)$ becomes increasingly large for those $X$ in the vacuum region. In the following, we obtain the eigenstates and eigenvalues of this mean-field Hamiltonian, that are localized at the boundary region. To this end, we Taylor-expand $V_{c}(X)$ around the boundary and keep only up to the linear term in the spatial coordinate,

[TABLE]

with positive $V_{1}$. We define the boundary, $L$, such that for $\theta=0$, $M(k_{y},k_{z}=0)$ changes its sign at $X=x=-k_{y}l^{2}-\frac{V_{1}}{m\omega^{2}}=\frac{L}{2}$. This definition of $L$ gives out,

[TABLE]

The Taylor expansion will be a priori justified, provided that the potential varies in space much slower than the magnetic length; $l\partial_{X}V_{c}\ll\hbar\omega$ (see below).

Under a proper basis change of the 2 by 2 Pauli matrices, $\hat{h}_{\rm 1D}(k_{y})$ thus given can be expressed in terms of raising and lower operators,

[TABLE]

with $\beta\equiv\sqrt{2\Delta_{0}V_{1}\sin\theta}$. The raising and lower operators, $[a,a^{\dagger}]=[b,b^{\dagger}]=1$, are defined in the following way,

[TABLE]

with $l_{\perp}\equiv\sqrt{\frac{\Delta_{0}}{V_{1}\sin\theta}}$, $\tilde{x}\equiv x-x_{0}$, $\tilde{z}\equiv z-z_{0}$ and

[TABLE]

The raising (lowering) operators, $a^{\dagger}$ ($a$) and $b^{\dagger}$ ($b$), have ladders of number states, $|n\rangle_{a}$, $|n\rangle_{b}$,

[TABLE]

These number states are functions only of $z$ and $x$. They are localized around $z=z_{0}$ and $x=x_{0}$ with localization length $l_{\perp}$ and $l$ respectively. Especially, $|0\rangle_{a}$ is given by the Airy function.

$\hat{h}_{\rm 1D}(k_{y})$ thus given has following set of eigenstates and eigenvalues;

[TABLE]

with $n\geq 1$ and $E_{n}\equiv\sqrt{2\Delta_{0}\sin\theta H|n|}$. The $k_{y}$ dependence is encoded into $x_{0}$ and $z_{0}$ in the number states. Irrespective of $k_{y}$, the eigenstates are localized around $X=\frac{L}{2}$ along the $X$-direction;

[TABLE]

Here $\lambda$ is a characteristic length scale with which the confining potential varies in space around the boundary, $\lambda V_{1}\equiv\hbar\omega$. Provided that $\lambda\gg l$, the eigenstates with different $k_{y}$ are all localized at $X=L/2$. The localized feature of the eigenstates a priori justifies the Taylor expansion of $V_{c}(X)$ around $X=L/2$ in $\hat{h}_{\rm 1D}(k_{y})$.

The eigenstates with different $k_{y}$ are energetically degenerate in each sLL and they are localized at different locations along the $Z$-axis,

[TABLE]

Accordingly, the degeneracy at each surface Landau level is proportional to an area of the side surface and the out-of-surface component of the magnetic field,

[TABLE]

In conclusion, the $(2+1)$D massless surface Dirac state in the topological excitonic insulator under the tilted magnetic field forms a sequence of the surface Landau levels, $E_{n}={\rm sgn}(n)\sqrt{2\Delta_{0}H_{\perp}|n|}$ ($n=\cdots,-1,0,1,\cdots$). Each surface Landau level has an degeneracy of $L_{y}L_{Z}\sin\theta/(2\pi l^{2})$, where $H_{\perp}$ is the out-of-surface component of the magnetic field, $H_{\perp}\equiv H\sin\theta$.

VIII summary and discussion on experiment

In this paper, we clarify comprehensive ground-state phase diagrams of a three-dimensional semimetal model in the quantum limit. The semimetal model has a pair of electron and hole pocket. We study two limiting cases at the charge neutrality point, (i) the model with screened Coulomb interaction and (ii) the model with an effective attractive interaction mediated by the screened electron-phonon interaction. The results show rich phase diagram structures as a function of the Fermi wavelength and the screening length (normalized by the magnetic length). In the repulsive interaction case, we found that an Ising-type spin density wave phase / excitonic insulator phase with ferro-type order of $XY$ spin moment is stabilized in the weak / strong screening regime respectively. In the attractive interaction case, we found that the plain charge density phase or a possible non-Fermi liquid phase is stabilized for weak screening regime, while from the intermediate to strong screening regimes, the ground state is dominated by charge Wigner crystal phase, topological excitonic insulator phase, and excitonic insulator with a two-dimensional vortex lattice of the $XY$ component of the spin-$1$ moment.

The topological excitonic insulator phase in the attractive interaction case is an three-dimensional interaction-driven topological band insulator in the quantum limit. Thereby, the odd-parity excitonic pairing in the bulk gives rise to a single copy of the $(2+1)$D massless surface Dirac fermion state at those surfaces parallel to the magnetic field. We show that when an in-plane transport is dominated by the surface transport through the $(2+1)$D massless Dirac state, the in-plane resistivity must show a $\sqrt{H_{\perp}}$-type surface SdH oscillation under canted magnetic field $H_{\perp}$.

Recently, a comprehensive resistivity measurement in graphite under high magnetic field has been carried out up to 90 T zhu19 . The graphite under the high magnetic field (the field $\perp$ the graphene plane) exhibits consecutive metal-insulator transitions as well as insulator-metal re-entrant transition at low temperature in an electric resistivity along the out-of-plane (field) direction yaguchi09 ; yaguchi98a ; fauque13 ; akiba15 ; arnold17 ; zhu19 . Experimentally, there exists two insulating phases, one insulating phase in the range of 30 T $<\!\ H\!\ <$ 53 T and the other in the range of 53 T $<\!\ H\!\ <$ 75 T. The recent experiment shows that the resistivity within the in-plane direction shows unusual ‘metallic’ behaviour in the second ‘insulating’ phase at 53 T $<\!\ H\!\ <$ 75 T yaguchi09 ; fauque13 ; akiba15 ; zhu19 .

The semimetal model studied in this paper can be applicable to the latter field regime (53 T $<\!\ H\!\ <$ 75 T), where an electron pocket locates around the $K$ point and a hole pocket locates around the $H$ point in the first Brillouin zone of the graphite under $H$. Previously, the authors argued that the low-$T$ insulator phase in the range of 53 T $<\!\ H\!\ <$ 75 T could be the topological EI phase, where the unusual metallic in-plane resistivity was attributed to the surface transport through the $(2+1)$D massless Dirac states pan18 .

A relevant electronic energy band in graphite under the high magnetic field has a band width of $40$ meV, the lattice constant along the out-of-plane direction $c_{0}$ is 6.7 $\AA$, and the relative permittivity $\varepsilon$ in graphite is from 9 to 16. We assume that $k_{F}=\pi/(6c_{0})$ for $H=64$ T. For $H=64$ T, this gives out $\hbar v_{F}=-\partial t\cos(kc_{0})/\partial k|_{k=k_{F}}=tc_{0}/2=10$ meV $\times 6.7\AA$ with $t=20$ meV, $\log_{10}B=2\log_{10}(2k_{F}l)\simeq 1.4$. From Eq. (57), $1/A=(4e^{2})/(\hbar\pi v_{F}\varepsilon)\simeq 22$ and $\log_{10}A\simeq-1.34$. The sound velocity in graphite along the $c_{0}$ axis is around 500 m/s. Carbon is 12 amu heavy ($M=12$ amu), and the density of carbon atom in graphite is $\rho_{0}=0.12...\AA^{-3}$. From low-carrier density feature in graphite in the zero field (at most $n=10^{18}$ cm*-3*), we take $Z=10^{-5}$. For $H=64$ T with $1/A=22$, this set of parameters give a ratio between the overall factor of the effective attractive interaction mediated by the screened electron-phonon interaction and that of the screened repulsive Coulomb interaction. It turns out that the effective attractive interaction is much smaller than the screened repulsive Coulomb interaction,

[TABLE]

The small value of the effective attractive interaction is mainly because of tiny electron valence of charged nucleus ion (carbon atom), $Z$. The tiny $Z$ even overcomes very large screening length, e.g. $l_{\rm scr}\equiv\sqrt{A}l=6.8\AA$ for $H=64$ T. Thereby, if we simply add these two interactions with the opposite signs at the initial RG scale, the repulsive interaction clearly dominates over the attractive interaction. This would be the case even if we used 100 times larger value of $Z$ than the value given above. From this observation, we consider in the following the case with only the repulsive Coulomb interaction.

In Fig. 1, the parameter point with $(\log_{10}A,\log_{10}B)\simeq(-1.34,1.4)$ corresponds to the EI phase with the broken U(1) spin rotational symmetry around the field direction. Please note that the excitonic pairing in the EI phase in Fig. 1 has the spatially even parity [see Eqs. (104,105)], and thereby it is non-topological EI instead of the topological EI. In fact, the non-topological EI phase seems to be consistent with the second ‘insulating’ phase in a recent graphite experiment. The recent in-plane resistivity experiment under the canted magnetic field does not show any SdH oscillation as a function of the canted component of the magnetic field zhu19 , unlike the expectation from the surface transport in the topological EI phase.

Depending on other factors, the excitonic pairing in the non-topological EI phase in the repulsive interaction case could be formed between electron band and hole band at different spatial coordinates within the $xy$ plane, as was the case in Sec. VIB and Fig. 9(c). For example, a ratio between $A$ and $A^{\prime}$ may not be $1$. $A$ itself could be smaller by several factors than the value given above, due to additional screening from the higher LLs and from the other valley in graphite.

The excitonic pairings between different spatial coordinates within the $xy$ plane could induce coherent carrier transports within the plane. Since the excitonic pairing is between $\uparrow$ spin electron-type band and $\downarrow$ spin hole-type band, the transport must be free from pinning effect due to charged impurities shayegan97 . Thereby, we can expect that such EI phase with broken translational symmetries within the $xy$ plane may give a simple theory explanation for the in-plane metallic bulk-transport behaviour in the second ‘insulating’ phase of 53 T $<H<$ 75 T in the graphite experiment yaguchi09 ; fauque13 ; akiba15 ; zhu19 . In fact, the recent transport experiment up to 90 T shows that the in-plane resistivity in the second ‘insulating’ phase is nearly constant in the field zhu19 . This observation is consistent with the two-dimensional $XY$-spin vortex lattice shown in Fig. 9(c) whose lattice constant is proportional to the magnetic length $l$. Since the lattice constant within the $xy$ plane is scaled by $l$, an Aharonov-Bohm (AB) flux that penetrates through a unit cell of the two-dimensional spin vortex lattice is independent of the field. This results in an absence of any SdH-like oscillation in the in-plane transport inside the second ‘insulating’ phase. Nonetheless, for further understandings of the unusual transport in graphite as well as the re-entrant insulator-metal transition, we need further theoretical studies and relevant results will be discussed elsewhere.

ACKNOWLEDGMENTS

RS acknowledge helpful and enlighting discussions with Zengwei Zhu, Miguel A. Cazalilla, Tomi Ohtsuki, Alexei Tsvelik, and Masatoshi Sato. The work was supported by NBRP of China (Grant No. 2014CB920901, Grant No. 2015CB921104 and Grant No. 2017A040215).

Appendix A RPA screening

The interaction potentials that carry the zero momentum along the field, i.e. $\Gamma_{\mu\nu}$ in Eq. (20), are screened by low-energy density fluctuations within each branch (‘right-mover’ or ‘left-mover’ branch) of the electron-type band or hole-type band. The screened interaction comprises of a sum of the bare interaction part and an effective interaction mediated by the density fluctuations fetter03 . According to the linear response theory, the effective interaction part is given by a retarded correlation functions between the density fluctuations;

[TABLE]

with $\mu,\nu,\lambda,\psi=e_{+},e_{-},h_{+},h_{-}$ and Eq. (26). The first term in the right-hand side is the bare interaction part and second term is the effective interaction part. The retarded correlation function $D^{R}_{\lambda\psi}(-p^{\prime},p,\omega)$ is obtained from a time-ordered correlation function in the static limit, $D^{R}_{\lambda\psi}(-p^{\prime},p,\omega=0)=D^{T}_{\lambda\psi}(-p^{\prime},p,\omega=0)$ with;

[TABLE]

Here $|\Psi_{0}\rangle$ is a many-body ground-state wavefunction and (real-)time dependence of the operator is in the Heisenberg picture. $\delta\hat{T}_{\mu}(Q_{1},Q_{2},q)$ is the density fluctuation operator within every branch $\mu=e_{+},e_{-},h_{+},h_{-}$,

[TABLE]

According to the Feynman-Dyson perturbation theory fetter03 ; mahan00 , the time-ordered correlation function is given by a proper part of the polarization function. The RPA approximates the proper part by its lowest order in the electron correlation. This gives out

[TABLE]

where the summation over $\rho=e_{+},e_{-},h_{+},h_{-}$ is omitted in the right hand side. In the static limit, the bare polarization function $\Pi_{0,\lambda}(\omega)$ for $\lambda=e_{+},e_{-},h_{+},h_{-}$ is given by Eqs. (45,46). In terms of the homomorphic nature of the interaction potential functional, Eq. (41), Eq. (179) with the RPA correlation function Eq. (180) reduces to

[TABLE]

where $\overline{I}_{0}(q_{x},k_{y})$ is given by Eq. (43).

The interaction potentials that carry $2k_{F}$ momentum along the field, $\Phi^{+-}_{\mu\nu}$ in Eq. (21), are also screened by $2k_{F}$ density fluctuations. As above, the screened interaction is characterized by the retarded density correlation function between the $2k_{F}$ density fluctuation operators;

[TABLE]

In the static limit ($\omega=0$), the retarded correlation function $D^{R,-+}_{\lambda\psi}(-p^{\prime},p,\omega)$ is identical to the corresponding time-ordered correlation function;

[TABLE]

$\delta\hat{S}^{\pm}_{\mu}$ is the $\pm 2k_{F}$ density fluctuation operator within electron pocket ($\mu=e$) or hole pocket ($\mu=h$),

[TABLE]

Within the RPA, the time-ordered correlation function is given by a bare polarization function that carries $2k_{F}$ momentum;

[TABLE]

where the polarization function at $p_{z}=2k_{F}$, $\Pi^{-+}_{0,\lambda}(\omega=0)$, is given by Eq. (47). In terms of the homomorphic relation, Eq. (42), Eq. (182) with Eq. (183) reduces to

[TABLE]

where $\overline{I}_{2k_{F}}(q_{x},k_{y})$ is given by Eq. (44).

Appendix B derivation of parquet RG equation

A derivation of the one-loop parquet RG equation can be implemented by a standard momentum shell renormalization. Thereby, we begin with a partition function of the interacting fermion model, Eqs. (59,60), and decompose the fermionic field into fast mode ($e_{\pm,>}$, $h_{\pm,>}$) and slow mode ($e_{\pm,<}$, $h_{\pm,<}$) in the momentum space

[TABLE]

with $\Lambda^{\prime}\equiv\Lambda e^{-\ln b}$. The integration of the fast mode in the partition function leads to a renormalization of the effective action for the slow mode,

[TABLE]

where

[TABLE]

and

[TABLE]

$S_{1,<}$ is the interaction part that is comprised only of the slow modes. $S_{1,>}$ is the other part of the interaction term that contains the fast modes. $\langle S_{1,>}\rangle_{0,>}$ in Eq. (191) renormalizes the Fermi velocity of the electron and hole pocket. Due to a particle-hole symmetry that exchanges the electron ahd hole bands ($m_{e}=m_{h}$), the renormalization of the Fermi velocity of the electron band and that of the hole band are identical to each other at the charge neutrality point. At the one-loop level of the renormalization group (RG), the Fermi velocity renormalization can be always absorbed into a scale change of the RG (see Eq. (202)). Thereby, we do not keep track of the Fermi velocity renormalization from $\langle S_{1,>}\rangle_{0,>}$ in Eq. (191).

$\langle S^{2}_{1,>}\rangle_{0,>}-\langle S_{1,>}\rangle^{2}_{0,>}$ gives rise to a renormalization of the interaction potentials. To calculate the renormalization, we have only to consider the following part of $S_{1,>}$,

[TABLE]

(the others do not contribute to the renormalization of the interaction potentials at the one-loop level RG).

$\langle S^{2}_{1,>}\rangle_{0,>}-\langle S_{1,>}\rangle^{2}_{0,>}$ in Eq. (191) gives out the following one-loop renormalization to the interaction potentials in Eqs. (60);

[TABLE]

with $d\Lambda\equiv\Lambda\ln b$ ${\bm{k}}\equiv(k,q)$, ${\bm{k}}^{\prime}\equiv(k^{\prime},q^{\prime})$ and $d{\bm{k}}^{\prime}\equiv dk^{\prime}dq^{\prime}$. After the integration of the fast modes, we scale the momentum along the field ($p$), single-particle frequency ($\omega$) and the field operators ($e_{\sigma}$ and $h_{\sigma}$) as

[TABLE]

This (tree-level) scale change keeps $S_{0,<}$ as well as $S_{1,<}$ to be invariant, while putting $\Lambda^{\prime}$ in $S_{0,<}$ and $S_{1,<}$ back to $\Lambda$. Accordingly, Eqs. (193,194,195,196) lead to the following one-loop renormalization group equations for the interaction potentials,

[TABLE]

with

[TABLE]

Note that the above one-loop RG equations as well as the initial forms of the interaction potentials, Eqs. (61,62,63), respect the following symmetries,

[TABLE]

Using these symmetries, the RG equations can be also written in Eqs. (64,65,66,67).

Consider the Fourier transform of $W_{\mu}({\bm{k}})$,

[TABLE]

for $\mu=b,d,e,g$ with

[TABLE]

and ${\bm{r}}\equiv(r_{x},r_{y})$, ${\bm{k}}\equiv(k,q)$. In terms of this dual representation, Eqs. (64,65,66,67) reduce to

[TABLE]

From Eqs. (61,62,63), the initial function forms for $F_{\mu}({\bm{r}})$ ($\mu=b,d,e$) and $\tilde{F}_{g}({\bm{r}})$ are as follows,

[TABLE]

with ${\bm{k}}\equiv(k,q)$. These initial forms as well as the RG equations in the dual space respect the following symmetries,

[TABLE]

Accordingly, Eqs. (206,207,208) can be rewritten into more symmetric forms,

[TABLE]

The RG equations thus obtained as well as the initial forms have the following O(2) symmetry;

[TABLE]

with $r\equiv|{\bm{r}}|$ for $\mu=b,d,e$ and arbitrary $\theta\in(0,2\pi]$. Utilizing this symmetry, we can reduce Eqs. (214,215,216) into the RG equations for $\Gamma_{\mu}(r)$ ($\mu=b,d,e$) and $\Gamma_{g}(r)$, Eqs. (77,78,79).