# Inevitable First Order Phase Transitions in 3D Quantum Hall Systems

Kaiyuan Gu  Department of Physics, Princeton University, Princeton, New Jersey
08544, USA  Kai Torrens  Department of Physics, Harvard University, Cambridge,
Massachusetts 02138, USA  Biao Lian  Department of Physics, Princeton
University, Princeton, New Jersey 08544, USA

###### Abstract

Recent experiments suggest that low carrier density three-dimensional (3D)
metals ZrTe5 and HfTe5 exhibit the 3D quantum Hall (QH) effect with Hall
resistivity plateaus and a metal-insulator transition in strong magnetic
fields. The conventional 3D QH theory requires a fixed period charge density
wave (CDW), which is however not observed experimentally. We investigate
alternative non-CDW mechanisms by considering a 3D metal in strong magnetic
fields with electrons coupled to a boson (e.g., phonon) field. We show that
the model exhibits inevitable first order phase transitions at jumps of the
number of occupied Landau level bands, which do not involve CDW. These
transitions may drive the system into a phase separation state with
percolation transitions. We further show this can lead to Hall resistivity
quasi-plateaus similar to that observed experimentally, and can provide a
natural explanation for the metal-insulator transition.

The celebrated integer quantum Hall (QH) effect of two-dimensional (2D)
electron gas exhibits a quantized Hall conductance plateaus in units of
$e^{2}/h$ [1] in terms of electron charge $e$ and Planck constant
$h=2\pi\hbar$, with coefficient given by the topological Chern number which
equals the number of occupied Landau levels (LLs) [2, 3]. Its generalization
to the 3D QH effect in 3D metals, which is topologically equivalent to
infinite layers of 2D QH effects, is however extremely difficult due to the
gapless bands of 3D electrons in magnetic field. A gap opening is needed for
the 3D QH effect, which was first proposed to come from a fixed spatial period
charge density wave (CDW) along the magnetic field [4, 5, 6, 7]. Similarly,
the 3D quantum anomalous Hall effect which needs no magnetic field relies on
the underlying 3D periodic lattice [8, 9, 10, 11, 12].

Recently, quasi-quantized plateaus of Hall resistivity was observed in low
carrier density anisotropic semimetals ZrTe5 and HfTe5 [13, 14, 15, 16, 17,
18, 19, 20, 21], suggesting the potential realization of the 3D QH effect. A
magnetic field independent CDW induced 3D QH effect has been proposed as a
theoretical mechanism [15, 16, 22, 23], in which the longitudinal resistivity
approaches zero. However, the CDW mechanism is unsupported by other
experiments showing quasi-quantized Hall resistivity but finite longitudinal
resistivity [16, 17, 18], and transport and spectroscopy measurements found no
evidence of CDW [18, 19, 21]. Even if the conventional Fermi surface Peierls
instability CDW is present, its period is magnetic field dependent and cannot
explain the 3D QH effect.

In this letter, we investigate the instabilities other than CDW in such 3D
anisotropic metals in magnetic field. We show that when electrons interact
with some boson (e.g., phonon) field, without assuming CDW, a first order
phase transition (of uniform strain if the bosons are acoustic phonons) is
inevitable at the transition of number of occupied Landau level (LL) bands at
zero disorders. Increasing the disorders could weaken and eventually eliminate
the first order transitions. Such first order phase transitions may drive the
system into a phase separation state exhibiting percolation transitions, which
make the Hall resistivity deviate from the linear magnetic field dependence
and exhibit quasi Hall plateaus similar to the experimental observations in
ZrTe5 and HfTe5. Our percolation theory further provides a natural explanation
for the metal-insulator transition observed at large magnetic field [15],
which deserves future studies.

The model. Previous studies indicate that ZrTe5 and HfTe5 are highly
anisotropic low electron density metals containing both quadratic dispersion
electrons at $\Gamma$ point and massive Dirac electrons at M point [24, 16,
19, 25]. For universality, we construct our model as dimensionless via a
proper rescaling (see below Eq. 2), and assume a dimensionless magnetic field
$B$ is applied in the $z$ direction. This yields LLs in the $x$-$y$ plane, and
the $n$-th LL ($n\geq 0$, $n\in\mathbb{Z}$) has a 1D dispersion with respect
to the $z$ direction dimensionless momentum $k_{z}$ (Fig. 1(c)):

$\epsilon_{k_{z},n}(B)=\begin{cases}\frac{1}{2}k_{z}^{2}+B\Big{(}n+\frac{1}{2}\Big{)}\
,\quad\quad\ \text{(quadratic)}\\\
M\sqrt{M^{2}+2Bn+k_{z}^{2}}-M^{2},\text{(Dirac)}\end{cases}$ (1)

where $M$ is the dimensionless Dirac mass gap. For simplicity, we assume only
one type of electrons (quadratic or Dirac), and assume negligible spin Zeeman
splitting, so that spin only contributes a factor of $2$ of degeneracy to
electron numbers. Taking these simplifications or not does not qualitatively
alter our theory.

Focusing on the quasi-1D interacting physics, we assume the quasi-1D LL bands
in Eq. 1 have a dimensionless Hamiltonian per volume
$H=H_{\text{e}}+H_{\text{int}}+H_{\text{boson}}$, where

$\begin{split}&H_{\text{e}}=\frac{B}{2\pi
L_{z}}\sum_{n,k_{z}}\epsilon_{k_{z},n}(B)\ c_{n,k_{z}}^{\dagger}c_{n,k_{z}}\
,\\\ &H_{\text{int}}=-\frac{B}{2\pi L_{z}}\sum_{n}\int_{0}^{L_{z}}dz\
D\Delta(z)\ c_{n}^{\dagger}(z)c_{n}(z)\ ,\\\
&H_{\text{boson}}=\frac{1}{L_{z}}\int_{0}^{L_{z}}dz\
\bigg{[}\frac{1}{2}\Delta(z)^{2}+\frac{\beta}{4}\Delta(z)^{4}\bigg{]}\
,\end{split}$ (2)

where we have assumed a $z$-dependent mean field
$\Delta(z)=\langle\hat{\Delta}(\mathbf{r})\rangle$ of some dimensionless boson
field $\hat{\Delta}(\mathbf{r})$ coupled to the electrons. $c_{n,k_{z}}$ is
the annihilation operator of a representative electron orbital of the $n$-th
quasi-1D LL band (spin degenerate) satisfying
$\\{c_{m,k_{1}},c_{n,k_{2}}^{\dagger}\\}=\delta_{mn}\delta_{k_{1},k_{2}}$, and
$c_{n}(z)=\frac{1}{\sqrt{L_{z}}}\sum_{k_{z}}e^{ik_{z}z}c_{n,k_{z}}$, with
$L_{z}$ being the dimensionless $z$-direction system size. $H_{\text{e}}$ is
the electron kinetic energy, $H_{\text{int}}$ is the electron-boson
interaction with coupling constant $D>0$, and $H_{\text{boson}}$ is the boson
field energy, in which a quartic term with $\beta>0$ is needed for the total
Hamiltonian to be lower bounded.

The dimensionless model Eq. 2 is from its corresponding dimensionful model
rewritten in units of characteristic energy $E_{0}$ and lengths $a_{L}$,
$a_{z}$, defined by $E_{0}=\zeta
Ya_{L}^{2}a_{z}=\hbar^{2}/m_{L}a_{L}^{2}=\hbar^{2}/m_{z}a_{z}^{2}$ (see
supplementary material (SM) [26] Sec. I). Here $\zeta>0$ is a dimensionless
number which can be chosen freely, $Y$ is the coefficient of the kinetic
energy $\frac{Y}{2}\Delta^{\prime}(z^{\prime})^{2}$ of the dimensionful boson
mean field $\Delta^{\prime}(z^{\prime})$ (primed notations denote dimensionful
quantities), and $m_{z}$ and $m_{L}=\sqrt{m_{x}m_{y}}$ are the dimensionful
$z$-direction and in-plane geometric mean effective Newtonian masses of
electrons. For Dirac electrons with dimensionful Dirac mass $M^{\prime}$ and
velocities $v_{x,y,z}$, one has $m_{z}=M^{\prime}/v_{z}^{2}$,
$m_{L}=M^{\prime}/v_{x}v_{y}$, and the dimensionless Dirac mass
$M=\sqrt{M^{\prime}/E_{0}}$. The dimensionless quantities in Eq. 2 are then
given by $B=eB^{\prime}a_{L}^{2}/\hbar$, $D=\sqrt{2\zeta}D^{\prime}/E_{0}$,
$\beta=2\zeta\beta^{\prime}/Y$,
$\Delta(z)=\Delta^{\prime}(z^{\prime})/\sqrt{2\zeta}$, where $B^{\prime}$,
$D^{\prime}$ and $\beta^{\prime}$ are the dimensionful magnetic field,
coupling constant and quartic coefficient, respectively. The free parameter
$\zeta$ implies one of the dimensionless parameters in Eq. 2 is redundant, but
we keep $\zeta$ as a free parameter for later convenience.

The boson field may originate from acoustic/optical phonons or other emergent
bosons in the system. In the case of acoustic phonons, the dimensionful boson
field
$\Delta^{\prime}(z^{\prime})=\langle\nabla^{\prime}\cdot\mathbf{u}\rangle$ is
the bulk strain, where $\mathbf{u}$ is the acoustic phonon field. Accordingly,
$Y$ and $D^{\prime}$ are the dimensionful bulk modulus and deformation
potential, respectively.

We also define the dimensionless 3D electron density
$n_{3D}=\frac{a_{z}a_{L}^{2}}{2}n_{3D}^{\prime}=\frac{B}{2\pi
L_{z}}\sum_{n,k_{z}}\langle c_{n,k_{z}}^{\dagger}c_{n,k_{z}}\rangle$,
dimensionless chemical potential $\mu=\mu^{\prime}/E_{0}$, and dimensionless
resistivity $\rho_{ij}=(e^{2}/ha_{z})\rho_{ij}^{\prime}$ ($i,j=x,y,z$), where
$n_{3D}^{\prime}$, $\mu^{\prime}$ and $\rho_{ij}^{\prime}$ are the
corresponding dimensionful quantities. Accordingly, the quasi-1D bands in Eq.
2 have an effective dimensionless 1D electron density
$n_{1D}=\frac{1}{L_{z}}\sum_{n,k_{z}}\langle
c_{n,k_{z}}^{\dagger}c_{n,k_{z}}\rangle=2\pi n_{3D}/B$.

Define $\Delta_{q}=\frac{1}{\sqrt{L_{z}}}\int dze^{-iqz}\Delta(z)$ as the
Fourier transform of the mean field. In the Peierls instability theory, if the
$n$-th quasi-1D LL band has Fermi momentum $k_{F,n}$, a nonzero CDW gap
$\Delta_{\pm 2k_{F,n}}$ will develop, resulting in an insulator. Since each
“state” of the quasi-1D bands is a LL with Chern number $1$, the total Chern
number of occupied states in the $x$-$y$ plane is
$C_{\text{tot}}=\sum_{n}2k_{F,n}L_{z}/2\pi=n_{1D}L_{z}=2\pi n_{3D}L_{z}/B$,
which leads to Hall and longitudinal resistivities

$\rho_{xy}=\frac{L_{z}}{C_{\text{tot}}}=\frac{B}{2\pi n_{3D}}\
,\quad\rho_{xx}=0\ .$ (3)

With $n_{3D}$ being a constant in the system, $\rho_{xy}\propto B$ gives no
plateaus. Thus, such Peierls CDWs, if existing, are irrelevant for
understanding the 3D Hall plateaus.

Figure 1: (a) $\frac{\partial\mathcal{F}}{\partial\Delta}$ of free energy
$\mathcal{F}$ at $B=1.9$ for $D=3$ and $\mu=1$, which is near the transition
between $N=1$ and $N=2$ phases. (b) $\mathcal{F}(\Delta)$ for different $B$,
where two local minima at $\Delta_{1}^{*}$ to $\Delta_{2}^{*}$ compete,
leading to a first order phase transition. (c) LL band dispersions of the two
ground states at equal free energy
$\mathcal{F}(\Delta_{1}^{*})=\mathcal{F}(\Delta_{2}^{*})$. The dashed line
denotes the Fermi level. (d) Ground state calculated at $B=2$ by keeping all
$\Delta_{q}$, which shows a major component at $q=0$ and small Peierls CDW
components at $q=\pm 2k_{F}$.

_First order phase transition._ As we will show below, the ground state of
model in Eq. 2 exhibits inevitable first order phase transitions when the
number of occupied quasi-1D LL bands jumps. We minimize the zero temperature
free energy $\mathcal{F}=\langle H\rangle-\mu n_{3D}$ at fixed chemical
potential $\mu$, allowing the electron density $n_{3D}$ to vary. The
minimization is done using the gradient descent numerical method with respect
to all $\Delta_{q}$ (see details in SM [26]). In the entire range of $B$ from
$0$ to the quantum limit, we find a predominantly large mean field component
$\Delta_{0}$ at $q=0$, as accompanied by small Peierls CDW components
$\Delta_{q}$ at $q=\pm 2k_{F,n}$, where $k_{F,n}$ is the Fermi momentum of the
$n$-th occupied quasi-1D band. An example is shown in Fig. 1(d). The small
Peierls CDWs agree with the Peierls theory prediction
$\Delta_{2k_{F}}\propto\exp(-2\pi^{2}k_{F}/B)$ (see SM [26]), which approaches
zero at large $k_{F}/B$.

Given the irrelevance of the small Peierls CDWs argued below Eq. 3, and the
lack of experimental evidence of CDWs [17, 18, 19, 21], hereafter we ignore
all the mean field components $\Delta_{q\neq 0}$, and only focus on the $q=0$
component $\Delta_{0}$. This amounts assuming a spatially uniform mean field
$\Delta(z)=\Delta_{0}/\sqrt{L_{z}}=\Delta$. By Eq. 2, this effectively shifts
the chemical potential of electrons from $\mu$ to $\mu+D\Delta$, and changes
the electron density $n_{3D}$ accordingly. In the $B\rightarrow 0$ limit,
where the electrons are 3D, for quadratic electron dispersion this yields a
free energy

$\lim_{B\rightarrow
0}\mathcal{F}(\Delta)=-\frac{2\sqrt{2}}{15\pi^{2}}(\mu+D\Delta)^{\frac{5}{2}}+\frac{1}{2}\Delta^{2}+\frac{\beta}{4}\Delta^{4}\
.$ (4)

Minimizing $\mathcal{F}$ would give a spontaneous nonzero $\Delta$.

At $B>0$, we define $N$ as the number of occupied (defined as having nonzero
electron density) quasi-1D LL bands, which decreases as $B$ increases. With
interaction $D>0$, assume the $B$ field is tuned such that the free energy
$\mathcal{F}$ is minimized at some $\Delta$ near
$\Delta_{N}=\frac{\epsilon_{N,0}(B)-\mu}{D}$. In this case, the $N$-th LL band
is occupied (unoccupied) when $\Delta>\Delta_{N}$ ($\Delta\leq\Delta_{N}$),
and has a Fermi momentum $k_{F,N}\propto\sqrt{\Delta-\Delta_{N}}$ if occupied.
In the vicinity of $\Delta_{N}$, this leads to a universal form of free energy
[26]:

$\mathcal{F}(\Delta)=\mathcal{F}_{N-1}(\Delta)-b(\Delta-\Delta_{N})^{\frac{3}{2}}\Theta(\Delta-\Delta_{N})\
,$ (5)

where $\mathcal{F}_{N-1}(\Delta)$ is a smooth function representing the free
energy of the first $N-1$ bands, $b=\sqrt{2}BD^{\frac{3}{2}}/3\pi^{2}$ for
quadratic dispersion, and $\Theta(x)$ is the Heaviside function which is $1$
for $x>0$ and $0$ for $x\leq 0$. This leads to
$\frac{\partial\mathcal{F}}{\partial\Delta}\simeq
c(\Delta-\Delta^{*}_{N-1})-\frac{3b}{2}\sqrt{\Delta-\Delta_{N}}\Theta(\Delta-\Delta_{N})$,
where $\Delta^{*}_{N-1}$ is the point $\mathcal{F}_{N-1}(\Delta)$ reaches its
minimum, and $c=\frac{\partial^{2}\mathcal{F}_{N-1}}{\partial\Delta^{2}}>0$ is
its second derivative. This implies
$\frac{\partial\mathcal{F}}{\partial\Delta}$ always increases (decreases) with
respect to $\Delta$ before (after) hitting $\Delta_{N}$, as shown in Fig.
1(a). Therefore, when $\Delta^{*}_{N-1}$ approaches $\Delta_{N}$ sufficiently
closely from below, two competing local minima $\Delta_{N-1}^{*}$ and
$\Delta_{N}^{*}$ with $\frac{\partial\mathcal{F}}{\partial\Delta}=0$
inevitably occur (Fig. 1(a),(b)), which have equal free energies
$\mathcal{F}(\Delta_{N-1}^{*})=\mathcal{F}(\Delta_{N}^{*})$ when the blue and
orange shaded areas in Fig. 1(a) are equal. Increasing $B$ thus leads to a
first order phase transition of $\Delta$ jumping from $\Delta_{N}^{*}$ to
$\Delta_{N-1}^{*}$, and the number of occupied bands jumping from $N$ to
$N-1$, as shown in Fig. 1(b)-(c).

Such first order phase transitions are inevitable as long as the LL band
bottoms are quadratic at small $k_{z}$. To verify this, we calculate the mean
field $\Delta^{*}(\mu,B)$ minimizing the free energy $\mathcal{F}$ of model
Eq. 2 with appropriate parameters estimated for ZrTe5/HfTe5. Hereafter, we
choose the free parameter $\zeta$ to fix $a_{L}\simeq 25.7$nm equaling to the
magnetic length of $1$ Tesla field, such that the dimensionless $B$ is simply
the dimensionful magnetic field $B^{\prime}$ in units of Tesla. By assuming
the bulk strain from acoustic phonon field contributes a major part of the
boson field $\Delta$ [27, 28, 25, 15, 29], we estimate $D=3$ as a legistimate
coupling strength (see SM [26]), and set $\beta=0.2$ (which does not
sensitively affect the phase diagram). For Dirac band, $M=3$ is chosen based
on [22]. The experiments [15, 16, 17, 18] estimate a carrier density
$n_{3D}^{\prime}\sim 10^{17}$cm-3, which gives $n_{3D}\sim 0.3$. Fig. 2(a) and
(b) show $\Delta^{*}$ as a function of $B$ at fixed chemical potential $\mu=1$
for quadratic band and $\mu=1.7-0.4B$ for Dirac band (chosen to keep $n_{3D}$
around $0.3$), respectively, which shows sharp first order jumps at the
transitions of the number of occupied LL bands $N$ as expected. The jumps in
$\Delta^{*}$ is generically larger for smaller $N$ (i.e., larger $B$). Fig.
2(c) and (d) show the phase diagrams for quadratic and Dirac bands,
respectively, where $N$ labels the number of occupied quasi-1D LL bands.

Figure 2:  Fixing $D=3$, $\beta=0.2$, the ground state $\Delta^{*}$ with
respect to $B$ calculated for (a) quadratic band with $\mu=1$, and (b) Dirac
band with $M=3$ and $\mu=1.7-0.4B$. The phase diagram with the same parameters
for (c) quadratic band and (d) Dirac band, with phases labeled by number of
occupied LL bands $N$. The red dashed lines show the contour of $n_{3D}=0.3$.

Phase separation and percolation. In phase diagrams Fig. 2(c),(d) calculated
by fixing chemical potential $\mu$, the 3D electron density $n_{3D}$ jumps
across a first order phase boundary. In the physical material, $n_{3D}$ is
fixed, while $\mu$ is not. Therefore, the system should evolve with $B$ along
the constant $n_{3D}$ contour, as shown by the red dashed line in Fig. 2(c)
and (d), for which $n_{3D}=0.3$. In particular, when the contour lies on the
phase boundary between phases $N-1$ and $N$, the system will enter a phase
separation state of both phases with spatial fractions $p_{N-1}$ and
$p_{N}=1-p_{N-1}$, respectively, which is robust in 3D [30]. The fractions are
determined by

$n_{3D}^{(N-1)}p_{N-1}+n_{3D}^{(N)}p_{N}=n_{3D}\ ,$ (6)

where $n_{3D}^{(N-1)}$ and $n_{3D}^{(N)}$ are the densities of phases $N-1$
and $N$ adjacent to the phase boundary ($n_{3D}^{(N)}>n_{3D}^{(N-1)}$). Fig.
3(a) shows the fractions $p_{N}$ calculated along the red dashed $n_{3D}=0.3$
contour in Fig. 2(c).

In the presence of weak disorders and random fields (varying slowly
spatially), we expect the two phases $N-1$ and $N$ to form random shaped
domains in space, breaking translation symmetry in all directions. By the
percolation theory in 3D, the regions of phase $N$ can percolate to infinity
only if its spatial fraction $p_{N}>p_{c}$, and similarly for phase $N-1$,
where the threshold fraction $p_{c}\simeq 0.3$ [31, 32, 33]. The percolation
theory is known to be crucial for understanding the Hall plateaus of 2D QH
effect [34, 35, 36, 37, 38], where a partially filled LL will contribute a
Hall conductance $\frac{e^{2}}{h}$ if percolating. For the phase $N$ in 3D
here, the total number of occupied LLs in the $x$-$y$ plane is
$C_{\text{tot}}^{(N)}=2\pi n_{3D}^{(N)}L_{z}/B$. Therefore, phase $N$ has
$C_{\text{tot}}^{(N)}-C_{\text{tot}}^{(N-1)}$ more occupied LLs than phase
$N-1$. These LLs would contribute a (dimensionless) Hall conductivity
$\sigma_{xy}^{(1)}=f(p_{N})\frac{C_{\text{tot}}^{(N)}-C_{\text{tot}}^{(N-1)}}{L_{z}}$,
where $f(p_{N})>0$ only if $p_{N}>p_{c}$, namely, the regions of phase $N$
percolate, thus their chiral surface states can reach the boundaries of the
system, and $f(p_{N})\rightarrow 1$ when $p_{N}\rightarrow 1$. The rest
$C_{\text{tot}}^{(N-1)}$ occupied LLs possessed by both phases extend across
the entire system and contribute a Hall conductivity
$\sigma_{xy}^{(0)}=\frac{C_{\text{tot}}^{(N-1)}}{L_{z}}$. With Eq. 6, the
total Hall conductivity $\sigma_{xy}=\sigma_{xy}^{(0)}+\sigma_{xy}^{(1)}$ is
given by

$\sigma_{xy}=\frac{2\pi}{B}\left[n_{3D}+(n_{3D}^{(N)}-n_{3D}^{(N-1)})\big{(}f(p_{N})-p_{N}\big{)}\right]\
.$ (7)

Figure 3: Along the $n_{3D}=0.3$ contour in Fig. 2(c), (a) spatial fractions
$p_{N}$ of different phases $N$ up to $4$ (labels in the legend), and (b)
$\rho_{xy}$ (blue) and $\rho_{xx}$ (red) with respect to $B$. The (b) inset:
$\rho_{xx}$ (blue) and $\frac{Bd\rho_{xy}}{dB}$ (red) versus $1/B$.

Assuming Peierls CDW gaps are absent or smaller than the disorder strength,
all the phases with $N\geq 1$ are metallic. This is consistent with the lack
of experimental evidence of CDW [18, 19, 21], and finite longitudinal
resistivity $\rho_{xx}$ in experiments [16, 17, 18] (except for $\rho_{xx}$
approaching zero in [15]). The only exception is the $N=0$ phase, which has
$n_{3D}^{(0)}=0$ and is insulating. Without knowing the microscopic details,
we assume an approximate Drude model form (dimensionless) longitudinal
resistivity $\rho_{xx}$:

$\rho_{xx}=\frac{1}{4\pi
n_{3D}}\bigg{[}\Gamma_{0}(B)+\Gamma_{s}\big{(}p_{N}(1-p_{N})\big{)}\bigg{]}$
(8)

where $\Gamma_{0}(B)=\gamma_{0}(1+\chi_{0}B^{2})$ is the dimensionless
electron scattering rate typical for magneto-resistance of metals, and
$\Gamma_{s}(p_{N}(1-p_{N}))$ is the extra scattering rate due to numerous
domain walls when the system is in a phase separation state of phases $N$ and
$N-1$, which is peaked at $p_{N}=\frac{1}{2}$. Since the jump of $\Delta^{*}$
across the domain walls between the two phases increases with $B$ (Fig. 2(a)),
$\Gamma_{s}$ should increase with $B$, and we take an ansatz
$\Gamma_{s}(x)=\gamma_{1}Bx^{2}$. Here
$\Gamma_{\alpha}=\hbar\Gamma_{\alpha}^{\prime}/E_{0}$ in terms of dimensionful
scattering rates $\Gamma_{\alpha}^{\prime}$ ($\alpha=0,s$). The Hall
resistivity from Eqs. 7 and 8 is then
$\rho_{xy}=\frac{1+\sqrt{1-4\rho_{xx}^{2}\sigma_{xy}^{2}}}{2\sigma_{xy}}$.

To demonstrate the result, for Eq. 7 we take $f(p)=\frac{1}{2}[1+\tanh(6p-3)]$
which is approximately nonzero only when $p>p_{c}\simeq 0.3$, and we set
$\gamma_{0}=0.04$, $\chi_{0}=0.25$ and $\gamma_{1}=2$, which give $\rho_{xx}$
one order smaller than $\rho_{xy}$ as is the case in experiments [15, 16, 17,
18]. Fig. 3(b) shows $\rho_{xx}$ and $\rho_{xy}$ as a function of $B$ along
the constant $n_{3D}=0.3$ contour in Fig. 2(c). We find that $\rho_{xy}$
exhibits quasi-plateaus when the system is in a phase separation state with
one of the two phases not percolating. In contrast, when the system is
uniformly in one phase, $\rho_{xy}$ approaches the linear $B$ curve in Eq. 3.
The $\rho_{xy}$ plateaus are approximately equally spaced in $1/B$, as Fig.
3(b) inset shows. The $\rho_{xy}$ curve we obtained qualitatively resemble the
experimental observations especially of HfTe5 [17].

Lastly, the experiments revealed a metal-insulator transition at large
magnetic field [15, 16, 17, 18], which can also be captured by our theory.
Recall that the $N=0$ phase is insulating with $n_{3D}^{(0)}=0$. At
sufficiently large $B$ ($\gtrsim 3$ in Fig. 3), the system enters a phase
separation state with phases $N=1$ and $N=0$ coexisting, as shown in Fig. 2(c)
and (d). As $B$ increases, the electron density $n_{3D}^{(1)}$ of phase $N=1$
increases, thus its spatial fraction $p_{1}$ decreases (Fig. 3(a)). Assume
$p_{1}=p_{c}$ at magnetic field $B_{c}$. The $N=1$ phase is then percolating
(not percolating) when $B<B_{c}$ ($B>B_{c}$), thus the entire system is
metallic (insulating), leading to a metal-insulator transition at $B=B_{c}$.
Taking $p_{c}=0.3$, we find $B_{c}\simeq 8$ along the $n_{3D}=0.3$ contour in
Fig. 2(c), the order of which matches the experiments well
($B_{c}^{\prime}\simeq 6.7$T in [15]).

_Discussion_. We have shown that the 3D metal coupled with a boson field in
the disorderless limit exhibits inevitable first order phase transitions in
strong magnetic fields, in which the boson mean field $\Delta^{*}$ jumps when
the number of occupied LL bands $N$ jumps. Fig. 2(a),(b) shows that the jump
in $\Delta^{*}$ increases with $B$. Since the electrons feel a chemical
potential shift $D\Delta^{*}$, if the system has short-range correlated
disorders of energy scale $W$, we expect the first order phase transition to
disappear when the jump in $\Delta^{*}$ is smaller than $W/D$. So the first
order phase transition will only survive above a threshold magnetic field set
by $W$.

We then showed that phase separation and percolation at the first order phase
transitions give a possible non-CDW explanation of the Hall resistivity quasi-
plateaus observed in 3D QH experiments in ZrTe5/HfTe5 [15, 16, 17, 18].
Moreover, the observed large magnetic field metal-insulator transition may be
explained as the percolation transition of the phase separation of $N=0$ and
$N=1$ phases. An interesting future questions is to investigate the critical
exponents at such interacting percolation transitions. Lastly, if the boson
field $\Delta^{*}$ originates from bulk strain, our results suggest jumps of
the (dimensionful) bulk strain up to $10^{-4}$ in ZrTe5/HfTe5 at the first
order phase transitions [26], and it will be intriguing if this can be
measured.

###### Acknowledgements.

_Acknowledgments_. We thank Claudia Felser and Priscila Rosa for helpful
discussions. This work is supported by the National Science Foundation through
Princeton University’s Materials Research Science and Engineering Center
DMR-2011750, and the National Science Foundation under award DMR-2141966.
Additional support is provided by the Gordon and Betty Moore Foundation
through Grant GBMF8685 towards the Princeton theory program.

## References

  * Klitzing  _et al._ [1980] K. v. Klitzing, G. Dorda, and M. Pepper, New method for high-accuracy determination of the fine-structure constant based on quantized hall resistance, Phys. Rev. Lett. 45, 494 (1980).
  * Laughlin [1981] R. B. Laughlin, Quantized hall conductivity in two dimensions, Phys. Rev. B 23, 5632 (1981).
  * Thouless  _et al._ [1982] D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Quantized hall conductance in a two-dimensional periodic potential, Phys. Rev. Lett. 49, 405 (1982).
  * Halperin [1987] B. I. Halperin, Possible States for a Three-Dimensional Electron Gas in a Strong Magnetic Field, Japanese Journal of Applied Physics 26, 1913 (1987).
  * Montambaux and Kohmoto [1990] G. Montambaux and M. Kohmoto, Quantized hall effect in three dimensions, Phys. Rev. B 41, 11417 (1990).
  * Kohmoto  _et al._ [1992] M. Kohmoto, B. I. Halperin, and Y.-S. Wu, Diophantine equation for the three-dimensional quantum hall effect, Phys. Rev. B 45, 13488 (1992).
  * Koshino  _et al._ [2001] M. Koshino, H. Aoki, K. Kuroki, S. Kagoshima, and T. Osada, Hofstadter butterfly and integer quantum hall effect in three dimensions, Phys. Rev. Lett. 86, 1062 (2001).
  * Haldane [2004] F. D. M. Haldane, Berry curvature on the fermi surface: Anomalous hall effect as a topological fermi-liquid property, Phys. Rev. Lett. 93, 206602 (2004).
  * Burkov and Balents [2011] A. A. Burkov and L. Balents, Weyl semimetal in a topological insulator multilayer, Phys. Rev. Lett. 107, 127205 (2011).
  * Wang  _et al._ [2016] J. Wang, B. Lian, and S.-C. Zhang, Dynamical axion field in a magnetic topological insulator superlattice, Phys. Rev. B 93, 045115 (2016).
  * Lian and Zhang [2017] B. Lian and S.-C. Zhang, Weyl semimetal and topological phase transition in five dimensions, Phys. Rev. B 95, 235106 (2017).
  * Yin  _et al._ [2020] J.-X. Yin, W. Ma, T. A. Cochran, X. Xu, S. S. Zhang, H.-J. Tien, N. Shumiya, G. Cheng, K. Jiang, B. Lian, Z. Song, G. Chang, I. Belopolski, D. Multer, M. Litskevich, Z.-J. Cheng, X. P. Yang, B. Swidler, H. Zhou, H. Lin, T. Neupert, Z. Wang, N. Yao, T.-R. Chang, S. Jia, and M. Zahid Hasan, Quantum-limit chern topological magnetism in tbmn6sn6, Nature 583, 533–536 (2020).
  * Liu  _et al._ [2016] Y. Liu, X. Yuan, C. Zhang, Z. Jin, A. Narayan, C. Luo, Z. Chen, L. Yang, J. Zou, X. Wu, S. Sanvito, Z. Xia, L. Li, Z. Wang, and F. Xiu, Zeeman splitting and dynamical mass generation in Dirac semimetal ZrTe5, Nature Communications 7, 12516 (2016).
  * Wang  _et al._ [2018] H. Wang, H. Liu, Y. Li, Y. Liu, J. Wang, J. Liu, J.-Y. Dai, Y. Wang, L. Li, J. Yan, D. Mandrus, X. C. Xie, and J. Wang, Discovery of log-periodic oscillations in ultraquantum topological materials, Science Advances 4, 10.1126/sciadv.aau5096 (2018).
  * Tang  _et al._ [2019] F. Tang, Y. Ren, P. Wang, R. Zhong, J. Schneeloch, S. A. Yang, K. Yang, P. A. Lee, G. Gu, Z. Qiao, and L. Zhang, Three-dimensional quantum Hall effect and metal–insulator transition in ZrTe5, Nature 569, 537 (2019).
  * Wang  _et al._ [2020] P. Wang, Y. Ren, F. Tang, P. Wang, T. Hou, H. Zeng, L. Zhang, and Z. Qiao, Approaching three-dimensional quantum Hall effect in bulk $\mathrm{HfT}{\mathrm{e}}_{5}$, Phys. Rev. B 101, 161201 (2020).
  * Galeski  _et al._ [2020] S. Galeski, X. Zhao, R. Wawrzyńczak, T. Meng, T. Förster, P. M. Lozano, S. Honnali, N. Lamba, T. Ehmcke, A. Markou, Q. Li., G. Gu, W. Zhu, J. Wosnitza, C. Felser, G. F. Chen, and J. Gooth, Unconventional Hall response in the quantum limit of HfTe5, Nature Communications 11, 5926 (2020).
  * Galeski  _et al._ [2021] S. Galeski, T. Ehmcke, R. Wawrzyńczak, P. M. Lozano, K. Cho, A. Sharma, S. Das, F. Küster, P. Sessi, M. Brando, R. Küchler, A. Markou, M. König, P. Swekis, C. Felser, Y. Sassa, Q. Li, G. Gu, M. V. Zimmermann, O. Ivashko, D. I. Gorbunov, S. Zherlitsyn, T. Förster, S. S. P. Parkin, J. Wosnitza, T. Meng, and J. Gooth, Origin of the quasi-quantized Hall effect in ZrTe5, Nature Communications 12, 3197 (2021).
  * Tian  _et al._ [2021] Y. Tian, N. Ghassemi, and J. H. Ross, Gap-Opening Transition in Dirac Semimetal ZrTe5, Phys. Rev. Lett. 126, 236401 (2021).
  * Wu  _et al._ [2023] W. Wu, Z. Shi, Y. Du, Y. Wang, F. Qin, X. Meng, B. Liu, Y. Ma, Z. Yan, M. Ozerov, C. Zhang, H.-Z. Lu, J. Chu, and X. Yuan, Topological Lifshitz transition and one-dimensional Weyl mode in HfTe5, Nature Materials 22, 84 (2023).
  * Piva  _et al._ [2024] M. M. Piva, R. Wawrzyńczak, N. Kumar, L. O. Kutelak, G. A. Lombardi, R. D. dos Reis, C. Felser, and M. Nicklas, Importance of the semimetallic state for the quantum Hall effect in HfTe5, Phys. Rev. Mater. 8, L041202 (2024).
  * Qin  _et al._ [2020] F. Qin, S. Li, Z. Z. Du, C. M. Wang, W. Zhang, D. Yu, H.-Z. Lu, and X. C. Xie, Theory for the charge-density-wave mechanism of 3d quantum hall effect, Phys. Rev. Lett. 125, 206601 (2020).
  * Geng  _et al._ [2021] H. Geng, G. Y. Qi, L. Sheng, W. Chen, and D. Y. Xing, Theoretical study of the three-dimensional quantum hall effect in a periodic electron system, Phys. Rev. B 104, 205305 (2021).
  * Zhang  _et al._ [2017] Y. Zhang, C. Wang, L. Yu, G. Liu, A. Liang, J. Huang, S. Nie, X. Sun, Y. Zhang, B. Shen, J. Liu, H. Weng, L. Zhao, G. Chen, X. Jia, C. Hu, Y. Ding, W. Zhao, Q. Gao, C. Li, S. He, L. Zhao, F. Zhang, S. Zhang, F. Yang, Z. Wang, Q. Peng, X. Dai, Z. Fang, Z. Xu, C. Chen, and X. J. Zhou, Electronic evidence of temperature-induced lifshitz transition and topological nature in zrte5, Nature Communications 8, 10.1038/ncomms15512 (2017).
  * Wang [2021] C. Wang, Thermodynamically induced transport anomaly in dilute metals ${\mathrm{zrte}}_{5}$ and ${\mathrm{hfte}}_{5}$, Phys. Rev. Lett. 126, 126601 (2021).
  * [26] See Supplemental Material for details.
  * Fu  _et al._ [2020] B. Fu, H.-W. Wang, and S.-Q. Shen, Dirac polarons and resistivity anomaly in ${\mathrm{zrte}}_{5}$ and ${\mathrm{hfte}}_{5}$, Phys. Rev. Lett. 125, 256601 (2020).
  * Kamm  _et al._ [1985] G. N. Kamm, D. J. Gillespie, A. C. Ehrlich, T. J. Wieting, and F. Levy, Fermi surface, effective masses, and dingle temperatures of ${\mathrm{zrte}}_{5}$ as derived from the shubnikov–de haas effect, Phys. Rev. B 31, 7617 (1985).
  * Jain  _et al._ [2013] A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, and K. a. Persson, The Materials Project: A materials genome approach to accelerating materials innovation, APL Materials 1, 011002 (2013).
  * Imry and Ma [1975] Y. Imry and S.-k. Ma, Random-field instability of the ordered state of continuous symmetry, Phys. Rev. Lett. 35, 1399 (1975).
  * Isichenko [1992] M. B. Isichenko, Percolation, statistical topography, and transport in random media, Rev. Mod. Phys. 64, 961 (1992).
  * Chalker and Dohmen [1995] J. T. Chalker and A. Dohmen, Three-dimensional disordered conductors in a strong magnetic field: Surface states and quantum hall plateaus, Phys. Rev. Lett. 75, 4496 (1995).
  * Song  _et al._ [2021] Z.-D. Song, B. Lian, R. Queiroz, R. Ilan, B. A. Bernevig, and A. Stern, Delocalization transition of a disordered axion insulator, Phys. Rev. Lett. 127, 016602 (2021).
  * Chalker and Coddington [1988] J. T. Chalker and P. D. Coddington, Percolation, quantum tunnelling and the integer hall effect, Journal of Physics C: Solid State Physics 21, 2665 (1988).
  * Pruisken [1988] A. M. M. Pruisken, Universal singularities in the integral quantum hall effect, Phys. Rev. Lett. 61, 1297 (1988).
  * Huckestein and Kramer [1990] B. Huckestein and B. Kramer, One-parameter scaling in the lowest landau band: Precise determination of the critical behavior of the localization length, Phys. Rev. Lett. 64, 1437 (1990).
  * Huo and Bhatt [1992] Y. Huo and R. N. Bhatt, Current carrying states in the lowest landau level, Phys. Rev. Lett. 68, 1375 (1992).
  * Wang  _et al._ [2014] J. Wang, B. Lian, and S.-C. Zhang, Universal scaling of the quantum anomalous hall plateau transition, Phys. Rev. B 89, 085106 (2014).

Supplemental Material for “Inevitable First Order Phase Transitions in 3D
Quantum Hall Systems”

## I. Derivation from the Dimensionful Hamiltonian to the Dimensionless
Hamiltonian

Here we derive the dimensionless Hamiltonian in the main text Eqs. (1) and (2)
from the dimensionful Hamiltonian. Assume in the absence of magnetic field,
the metal has the following dimensionful 3D Hamiltonian
$H^{\prime}=H^{\prime}_{\text{e}}+H^{\prime}_{\text{int}}+H^{\prime}_{\text{Boson}}$
_per unit volume_ :

$\begin{split}&H^{\prime}_{\text{e}}=\frac{1}{L_{x}^{\prime}L_{y}L^{\prime}_{z}}\sum_{s=\uparrow,\downarrow}\int
d^{3}\mathbf{r}^{\prime}\
\psi_{s}^{\dagger}(\mathbf{r}^{\prime})\mathcal{H}^{\prime}(-i\nabla^{\prime})\psi_{s}(\mathbf{r}^{\prime})\
,\\\
&H^{\prime}_{\text{int}}=-\frac{1}{L_{x}^{\prime}L_{y}L^{\prime}_{z}}\sum_{s=\uparrow,\downarrow}\int
d^{3}\mathbf{r}^{\prime}\
D^{\prime}\hat{\Delta}^{\prime}(\mathbf{r}^{\prime})\psi_{s}^{\dagger}(\mathbf{r}^{\prime})\psi_{s}(\mathbf{r}^{\prime})\\\
&H^{\prime}_{\text{boson}}=\frac{1}{L_{x}^{\prime}L_{y}L^{\prime}_{z}}\int
d^{3}\mathbf{r}^{\prime}\
\bigg{\\{}\frac{Y}{2}\hat{\Delta}^{\prime}(\mathbf{r}^{\prime})^{2}+\frac{\beta^{\prime}}{4}\hat{\Delta}^{\prime}(\mathbf{r}^{\prime})^{4}+\frac{1}{2}\Big{[}\partial_{t}P(\hat{\Delta}^{\prime}(\mathbf{r}^{\prime}))\Big{]}^{2}\bigg{\\}}\end{split}$
(S1)

where all the primed notations are dimensionful physical quantities (meant to
be distinguished with the corresponding dimensionless quantities defined
later): $L_{\alpha}^{\prime}$, $k_{\alpha}^{\prime}$ ($\alpha=x,y,z$) are the
system size and momentum in the $\alpha$-th direction,
$\mathbf{r}^{\prime}=(x^{\prime},y^{\prime},z^{\prime})$ is the real space
position, and $s$ is the spin. $\psi_{s}(\mathbf{r}^{\prime})$ is the electron
field of spin $s$, while $\hat{\Delta}^{\prime}(\mathbf{r}^{\prime})$ is the
boson field coupled to it, with $D^{\prime}$, $Y$ and $\beta^{\prime}$ being
the corresponding dimensionful parameters.
$P(\hat{\Delta}^{\prime}(\mathbf{r}^{\prime}))$ is a certain function of the
boson field depending on the boson theory, and the term
$\frac{1}{2}\Big{[}\partial_{t}P(\hat{\Delta}^{\prime}(\mathbf{r}^{\prime}))\Big{]}^{2}$
is meant to represent the quantum dynamical term of the boson field
$\hat{\Delta}^{\prime}(\mathbf{r}^{\prime})$. The kinetic energy
$\epsilon^{\prime}(\mathbf{k})$ of electrons, namely, the eigenvalue of the
electron $\mathcal{H}^{\prime}(-i\nabla^{\prime})$ is assumed to take the
(spin-independent) form of either quadratic or a massive Dirac (conduction
band only) dispersion:

$\epsilon^{\prime}(\mathbf{k}^{\prime})=\begin{cases}\frac{\hbar^{2}k_{x}^{\prime
2}}{2m_{x}}+\frac{\hbar^{2}k_{y}^{\prime
2}}{2m_{y}}+\frac{\hbar^{2}k_{z}^{\prime 2}}{2m_{z}}\
,\qquad\qquad\quad\qquad\quad\quad\ \text{(quadratic)}\\\ \sqrt{M^{\prime
2}+\hbar^{2}(v_{x}^{2}k_{x}^{\prime 2}+v_{y}^{2}k_{y}^{\prime
2}+v_{z}^{2}k_{z}^{\prime 2})}-M^{\prime},\qquad\text{(Dirac)}\end{cases}$
(S2)

where $m_{\alpha}$ ($\alpha=x,y,z$) is the effective Newtonian mass in the
$\alpha$-th direction; $M^{\prime}$ is the Dirac mass (gap), and $v_{\alpha}$
($\alpha=x,y,z$) is the Dirac velocity in the $\alpha$-th direction. We have
set the band bottom to be at the zero energy.

We note that when the bosons are acoustic phonons, we have
$\hat{\Delta}^{\prime}(\mathbf{r}^{\prime})=\nabla^{\prime}\cdot\mathbf{u}(\mathbf{r}^{\prime})$
where $\mathbf{u}(\mathbf{r}^{\prime})$ is the acoustic phonon field. $Y$ and
$D^{\prime}$ are the bulk modulus and deformation potential of the material,
respectively.

When a (dimensionful) magnetic field $B^{\prime}$ is added in the
$z^{\prime}$-direction, Landau levels (LLs) are formed in the
$x^{\prime}$-$y^{\prime}$ plane, and $k_{z}^{\prime}\in
2\pi\mathbb{Z}/L_{z}^{\prime}$ in the $z$-direction remains a good quantum
number. For simplicity, we ignore the spin Zeeman splitting energy and assumes
the spin is degenerate, which does not affect the main message of our paper.
The spin degenerate single-electron kinetic energy of the $n$-th LL band is
then given by

$\epsilon_{k_{z}^{\prime},n}^{\prime}(B^{\prime})=\begin{cases}\frac{\hbar^{2}k_{z}^{\prime
2}}{2m_{z}}+\frac{\hbar eB^{\prime}}{m_{L}}\Big{(}n+\frac{1}{2}\Big{)}\
,\qquad\qquad\qquad\qquad\ \text{(quadratic)}\\\ \sqrt{M^{\prime 2}+2\hbar
v_{x}v_{y}eB^{\prime}n+\hbar^{2}v_{z}^{2}k_{z}^{\prime 2}}-M^{\prime}\
,\qquad\text{(Dirac)}\end{cases}$ (S3)

where $m_{L}=\sqrt{m_{x}m_{y}}$ is the geometric mean Newtonian mass in the
$x^{\prime}$-$y^{\prime}$ plane, and $n\geq 0$. As we did in the main text, we
first impose translational symmetry in the $x^{\prime}$ and $y^{\prime}$
directions, and assume the boson field develops a time-independent mean field
$\langle\hat{\Delta}^{\prime}(\mathbf{r}^{\prime})\rangle=\Delta^{\prime}(z^{\prime})$.
We define $c^{\dagger}_{n,k_{z}^{\prime}}$ and $c_{n,k_{z}^{\prime}}$ as the
creation and annihilation operators of a representative electron eigenstate
orbital of the $n$-th LL with $z$-direction momentum $k_{z}^{\prime}$ and a
particular spin (e.g., the eigenstate orbital at the origin of the
$x^{\prime}$-$y^{\prime}$ plane with spin up), which satisfies
$\\{c_{n,k_{z}^{\prime}},c^{\dagger}_{n,k_{z}^{\prime\prime}}\\}=\delta_{n,n^{\prime}}\delta_{k_{z}^{\prime},k_{z}^{\prime\prime}}$.
Such an orbital has an area $2\pi l_{B}^{\prime 2}=2\pi\hbar/eB^{\prime}$ in
the $x^{\prime}$-$y^{\prime}$ plane, where
$l_{B}^{\prime}=\sqrt{\hbar/eB^{\prime}}$ is the magnetic length. Accordingly,
for each given $n$ and $k_{z}$, the electron states have a LL degeneracy
$N_{L}=2L_{x}^{\prime}L_{y}^{\prime}/2\pi l_{B}^{\prime
2}=2L_{x}^{\prime}L_{y}^{\prime}eB^{\prime}/2\pi\hbar$ in the
$x^{\prime}$-$y^{\prime}$ plane, where the factor $2$ comes from spin
degeneracy. The translational symmetry in the $x^{\prime}$ and $y^{\prime}$
directions requires each (spin degenerate) LL to be uniformly
occupied/unoccupied in the $x^{\prime}$-$y^{\prime}$ plane, in which the
number of electrons is effectively equal to
$N_{L}c^{\dagger}_{n,k_{z}^{\prime}}c_{n,k_{z}^{\prime}}$. Thus, we can
effectively rewrite the dimensionful Hamiltonian per unit volume in the mean
field approximation as

$\begin{split}&H^{\prime}_{\text{e}}=\frac{N_{L}}{L_{x}^{\prime}L_{y}L^{\prime}_{z}}\sum_{n,k^{\prime}_{z}}\epsilon_{k_{z}^{\prime},n}^{\prime}(B^{\prime})\
c_{n,k^{\prime}_{z}}^{\dagger}c_{n,k^{\prime}_{z}}=\frac{2eB^{\prime}}{2\pi\hbar
L^{\prime}_{z}}\sum_{n,k^{\prime}_{z}}\epsilon_{k_{z}^{\prime},n}^{\prime}(B^{\prime})\
c_{n,k^{\prime}_{z}}^{\dagger}c_{n,k^{\prime}_{z}}\ ,\\\
&H^{\prime}_{\text{int}}=-\frac{2eB^{\prime}}{2\pi\hbar
L^{\prime}_{z}}\sum_{n}\int_{0}^{L_{z}^{\prime}}dz^{\prime}\
D^{\prime}\Delta^{\prime}(z^{\prime})\
c_{n}^{\dagger}(z^{\prime})c_{n}(z^{\prime})\ ,\\\
&H^{\prime}_{\text{boson}}=\frac{1}{L^{\prime}_{z}}\int_{0}^{L_{z}^{\prime}}dz^{\prime}\
\bigg{[}\frac{Y}{2}\Delta^{\prime}(z^{\prime})^{2}+\frac{\beta^{\prime}}{4}\Delta^{\prime}(z^{\prime})^{4}\bigg{]}\
,\end{split}$ (S4)

where we have defined the 1D Fourier transformed electron operator
$c_{n}^{\dagger}(z^{\prime})=\frac{1}{\sqrt{L_{z}^{\prime}}}\sum_{k_{z}^{\prime}}e^{ik_{z}^{\prime}z^{\prime}}c_{n,k^{\prime}_{z}}^{\dagger}$
for the $n$-th LL.

We now redefine the Hamiltonian model into a dimensionless form by rescaling.
We define a characteristic energy $E_{0}$ and characteristic lengths $a_{L}$
in the $x$-$y$ plane and $a_{z}$ in the $z$ direction by requiring:

$E_{0}=\zeta
Ya_{L}^{2}a_{z}=\frac{\hbar^{2}}{m_{L}a_{L}^{2}}=\frac{\hbar^{2}}{m_{z}a_{z}^{2}}\
,$ (S5)

where $\zeta>0$ is a dimensionless positive number one is free to choose
(which will be chosen properly later). For quadratic electron band,
$m_{L}=\sqrt{m_{x}m_{y}}$ and $m_{z}$ are simply the effective Newtonian
masses. For Dirac electron band, the Newtonian masses are defined as
$m_{z}=M^{\prime}/v_{z}^{2}$ and $m_{L}=M^{\prime}/v_{x}v_{y}$, since in the
small momentum limit the dispersion in Eq. S2 takes the quadratic form
$\epsilon^{\prime}(\mathbf{k}^{\prime})\simeq\frac{\hbar^{2}}{2M^{\prime}}(v_{x}^{2}k_{x}^{\prime
2}+v_{y}^{2}k_{y}^{\prime 2}+v_{z}^{2}k_{z}^{\prime 2})$. Solving Eq. S5
yields the following expressions:

$a_{L}=\left(\frac{\hbar^{2}m_{z}^{1/2}}{\zeta Ym_{L}^{3/2}}\right)^{1/5}\
,\qquad a_{z}=\left(\frac{\hbar^{2}m_{L}}{\zeta Ym_{z}^{2}}\right)^{1/5}\
,\qquad
E_{0}=\left(\frac{\hbar^{6}\zeta^{2}Y^{2}}{m_{L}^{2}m_{z}}\right)^{1/5}\ .$
(S6)

We then define the following dimensionless quantities (which do not have prime
in notation):

$\begin{split}&H=\frac{H^{\prime}}{2\zeta Y}\ ,\quad
L_{z}=\frac{L_{z}^{\prime}}{a_{z}}\ ,\quad z=\frac{z^{\prime}}{a_{z}}\ ,\quad
k_{z}=k_{z}^{\prime}a_{z}\ ,\quad B=\frac{a_{L}^{2}}{l_{B}^{\prime
2}}=\frac{ea_{L}^{2}}{\hbar}B^{\prime}\
,\quad\epsilon_{k_{z},n}(B)=\frac{\epsilon_{k_{z}^{\prime},n}^{\prime}(B^{\prime})}{E_{0}}\
,\ M=\sqrt{\frac{M^{\prime}}{E_{0}}}\ ,\\\
&c^{\dagger}_{n,k_{z}}=c^{\dagger}_{n,k_{z}^{\prime}}\ ,\quad
c_{n}^{\dagger}(z)=\frac{1}{\sqrt{L_{z}}}\sum_{k_{z}}e^{ik_{z}z}c_{n,k_{z}}^{\dagger}=\sqrt{a_{z}}c_{n}^{\dagger}(z^{\prime})\
,\quad\Delta(z)=\frac{\Delta^{\prime}(z^{\prime})}{\sqrt{2\zeta}}\
,\quad\beta=\frac{2\zeta\beta^{\prime}}{Y}\ ,\quad
D=\frac{\sqrt{2\zeta}D^{\prime}}{E_{0}}\ .\end{split}$ (S7)

In addition, assume $\mu^{\prime}$ and $n_{3D}^{\prime}$ are the dimensionful
physical chemical potential and 3D electron density, respectively. We define
the dimensionless chemical potential $\mu$ and dimensionless electron density
$n_{3D}$ as

$\mu=\frac{\mu^{\prime}}{E_{0}}\ ,\qquad
n_{3D}=\frac{n_{3D}^{\prime}a_{L}^{2}a_{z}}{2}\ ,$ (S8)

After such a rescaling, the dimensionless electron dispersions becomes

$\epsilon_{k_{z},n}(B)=\begin{cases}\frac{1}{2}k_{z}^{2}+B\Big{(}n+\frac{1}{2}\Big{)}\
,\quad\quad\qquad\qquad\text{(quadratic)}\\\
M\sqrt{M^{2}+2Bn+k_{z}^{2}}-M^{2},\qquad\ \text{(Dirac)}\end{cases}$ (S9)

where for the Dirac band the dimensionless mass is $M=M^{\prime}a_{z}/\hbar
v_{z}=\sqrt{M^{\prime}/E_{0}}$ as define in Eq. S7. The dimensionless
Hamiltonian per unit volume $H=H_{\text{e}}+H_{\text{int}}+H_{\text{boson}}$
after the rescaling takes the form

$\begin{split}&H_{\text{e}}=\frac{B}{2\pi
L_{z}}\sum_{n,k_{z}}\epsilon_{k_{z},n}(B)\ c_{n,k_{z}}^{\dagger}c_{n,k_{z}}\
,\\\ &H_{\text{int}}=-\frac{B}{2\pi L_{z}}\sum_{n}\int_{0}^{L_{z}}dz\
D\Delta(z)\ c_{n}^{\dagger}(z)c_{n}(z)\ ,\\\
&H_{\text{boson}}=\frac{1}{L_{z}}\int_{0}^{L_{z}}dz\
\bigg{[}\frac{1}{2}\Delta(z)^{2}+\frac{\beta}{4}\Delta(z)^{4}\bigg{]}\
,\end{split}$ (S10)

as given in the main text. The dimensionless 3D electron density in Eq. S8
satisfies

$n_{3D}=\frac{B}{2\pi L_{z}}\sum_{n,k_{z}}\ \langle
c_{n,k_{z}}^{\dagger}c_{n,k_{z}}\rangle\ .$ (S11)

In particular, we note that the form of the dimensionless model in Eq. S10 is
independent of the choice of the dimensionless number $\zeta$ in Eq. S5. This
implies that the three dimensionless parameters $B$, $D$ and $\beta$ are not
independent: we can set one of the parameters to $1$ by choosing a number
proper $\zeta$. However, we keep all the three parameters for convenience.

We now take the estimated physical parameters for ZrTe5 and HfTe5, and
estimate the dimensionless parameters here, which we use for numerical
calculations. For this purpose, we assume the boson field
$\Delta^{\prime}=\langle\nabla^{\prime}\cdot\mathbf{u}(\mathbf{r}^{\prime})\rangle$
is from acoustic phonon field $\mathbf{u}(\mathbf{r}^{\prime})$, and thus
$D^{\prime}$ is the deformation potential, and $Y$ is the bulk modulus. In
principle, the other boson fields such as optical phonons can also contribute,
which we assume are of similar order of magnitudes. From the literature [27,
28, 25, 15, 22, 29], we adopt the following order estimations of the physical
parameters:

$D^{\prime}\sim 10\text{eV}\ ,\quad Y\sim 5\times 10^{9}\text{Pa}\simeq
3.1\times 10^{28}\text{eV}\cdot\text{m}^{-3}\ ,\quad m_{L}\simeq 0.2m_{e}\
,\quad m_{z}\simeq 2m_{e}\ ,\quad M^{\prime}\simeq 5\text{meV}\ ,$ (S12)

where $m_{e}\simeq 0.91\times 10^{30}$kg is the bare electron mass, which
yields $\hbar^{2}/m_{e}\simeq 1.22\times 10^{-38}$J$\cdot$m${}^{2}=7.63\times
10^{-20}$eV$\cdot$m2. Allowing $\zeta$ to be a free parameter, this yields

$a_{L}\simeq 0.52\zeta^{-1/5}\text{nm}\ .$ (S13)

For convenience, we choose $\zeta$ to set $a_{L}$ to be equal to the magnetic
length $\sqrt{\hbar\times 1\text{T}/e}=25.7$nm for a 1 Tesla magnetic field,
which gives:

$a_{L}\simeq 0.52\zeta^{-1/5}\text{nm}\simeq 25.7\text{nm}\
,\qquad\rightarrow\qquad\zeta\simeq 3.4\times 10^{-9}\ .$ (S14)

In this way, the dimensionless magnetic field $B$ is simply the dimensionful
magnetic field $B^{\prime}$ measured in Tesla, which we consider as a
convenient choice. Under this choice, we derive the characteristic scales in
Eq. S6 and dimensionless parameters in Eq. S7 to be

$a_{L}\simeq 25.7\text{nm}\ ,\quad a_{z}\simeq 8.1\text{nm}\ ,\quad
E_{0}\simeq 0.5\text{meV}\ ,\quad D=\frac{\sqrt{2\xi}D^{\prime}}{E_{0}}\sim 2\
,\quad M=\sqrt{\frac{M^{\prime}}{E_{0}}}\sim 3\ .$ (S15)

In the experiment, the 3D electron density is around $n_{3D}^{\prime}\sim
10^{17}$cm${}^{-3}=10^{4}$nm-3, which gives the dimensionless density

$n_{3D}=\frac{n_{3D}^{\prime}a_{L}^{2}a_{z}}{2}\sim 0.3\ .$ (S16)

Considering that the coupling to the other boson fields may further enhance
the effective coupling constant $D$, we choose

$D=3\ ,\qquad n_{3D}=0.3$ (S17)

for our calculations. The parameter $\beta$ of the quartic boson field term is
not reported in literature, and is solely for the purpose to lower bound the
Hamiltonian. For small $\beta$, our calculations are not sensitive to $\beta$.
The value of $\beta=0.2$ chosen in our calculations belong to this range.

An additional note: from the above derivations, $E_{0}\simeq 0.5$meV gives the
magnitude of the dimensionful LL spacings at physical magnetic field
$B^{\prime}\simeq 1$T. At the same time, the spin Zeeman splitting is on the
order of $\mu_{B}B^{\prime}\simeq 0.1$meV for $B^{\prime}\simeq 1$T. This
justifies our simplifying approximation that the spin Zeeman splitting is
ignorable compared to the orbital LL spacing.

## II. Numerical Calculation of ground states with respect to all $\Delta_{q}$

In this section, we show the numerical calculation of the ground states of the
system by optimizing the free energy $\mathcal{F}(\Delta_{q})=\langle
H\rangle-\mu n_{3D}$ at fixed $\mu$ with respect to all $\Delta_{q}$, where
$\Delta_{q}$ are the Fourier components of
$\Delta(z)=\frac{1}{\sqrt{L_{z}}}\sum_{q}\Delta_{q}e^{iqz}$. We rewrite the
dimensionless Hamiltonian (Eq. S10) in momentum space:

$\begin{split}&H_{\text{e}}=\frac{B}{2\pi
L_{z}}\sum_{n,k_{z}}\epsilon_{k_{z},n}(B)\ c_{n,k_{z}}^{\dagger}c_{n,k_{z}}\
,\\\ &H_{\text{int}}=-\frac{B}{2\pi
L_{z}}\sum_{n,k_{z},q}\frac{D\Delta_{q}}{\sqrt{L_{z}}}c_{n,k_{z}+q}^{\dagger}c_{n,k_{z}}\\\
&H_{\text{boson}}=\frac{1}{L_{z}}\bigg{[}\frac{1}{2}\sum_{q}\Delta_{q}\Delta_{-q}+\frac{\beta}{4L_{z}}\sum_{q_{1},q_{2},q_{3}}\Delta_{q_{1}}\Delta_{q_{2}}\Delta_{q_{3}}\Delta_{-(q_{1}+q_{2}+q_{3})}\bigg{]}\end{split}$
(S18)

We allow all $\Delta_{q}$ to be complex in our calculations, which are
required to satisfy $\Delta_{q}=\overline{\Delta}_{-q}$, where $\overline{c}$
stands for the complex conjugate of $c$.

Denote $\Psi_{n}=(\\{c_{n,k_{z}}\\})^{T}$ as the base of $n^{th}$ Landau band,
then the Hamiltonian subtracting $\mu\hat{n}_{3D}$ ($\hat{n}_{3D}$ is the 3D
electron density operator) can be written as

$H-\mu\hat{n}_{3D}=\frac{B}{2\pi
L_{z}}\sum_{n}\Psi_{n}^{\dagger}\big{(}\mathcal{H}_{e}+\sum_{q}\mathcal{H}_{q}-\mu\big{)}\Psi_{n}+\sum_{q}\Delta_{q}\Omega_{q}$
(S19)

where $\mathcal{H}_{e}$ is a diagonal matrix for electron energy term
$\varepsilon_{0}$, and $\mathcal{H}_{q}$ accounts for the electron-boson mean
field coupling, including the uniform term at $q=0$, and

$\Omega_{q}=\frac{1}{2L_{z}}\Delta_{-q}+\frac{\beta}{4L_{z}^{2}}\sum_{q_{1}+q_{2}+q_{3}=-q}\Delta_{q_{1}}\Delta_{q_{2}}\Delta_{q_{3}}\
.$ (S20)

In our numerical calculation, We sample $k_{z}$ with the standard spacing
$2\pi/L_{z}$. For each configuration $\Delta_{q}$, we diagonalize
$\mathcal{H}=\mathcal{H}_{e}+\sum_{q}\mathcal{H}_{q}$, and calculate the total
energy by summing all the occupied electron states with eigenvalues smaller
than $\mu$ for each Landau band $n$, after which we get the free energy
$\mathcal{F}(\Delta_{q})$. Then, we iterate using the gradient descent method
to minimize the the free energy for each iteration: in the $j$-th iteration
step, we update the boson mean field by gradient descent method,
$\Delta_{q}^{(j)}=\Delta_{q}^{(j-1)}-\gamma\frac{\delta\mathcal{F}(\Delta_{q}^{(j-1)})}{\delta\overline{\Delta}_{q}}$
with a suitable length of step $\gamma$, which would converge to the ground
state or some local minimum. By choosing sufficiently many initial
configurations for iteration, we select the lowest final state as the ground
state.

For quadratic Landau bands with $L_{z}=100$, $\mu=1$, $B=2$, $D=3$ and
$\beta=0.2$, the result is shown in the main text Fig. 1(d), where both
$\Delta_{0}$ and the Peierls CDW gap $\Delta_{2k_{F}}$ are found to be
nonzero, but $|\Delta_{0}|$ is much larger than $|\Delta_{2k_{F,n}}|$. For
different $B$, the amplitudes of the two peaks are shown in Fig. S1. We note
that the CDW gaps $\Delta_{2k_{F,n}}$ are generically much smaller than
$\Delta_{0}$ except when $k_{F,n}$ is very small, namely, close to the first
order phase transitions. However, when $k_{F,n}$ is very small,
$\Delta_{2k_{F,n}}$ has almost zero momentum $q=2k_{F,n}\approx 0$, which is
not quite distinguishable from $\Delta_{0}$, so we consider them playing the
same role as giving an almost uniform mean field $\Delta(z)=\Delta$. This
justifies our approximation of keeping only $\Delta_{0}$ in the main text.

Figure S1: Comparison of energy shift by $\Delta_{0}$ and the CDW gap for
different $B$.

We can also further verify the above results and the smallness of the Peierls
CDW gaps by analytical calculations in the quantum limit. Assume there is only
one LL band occupied, and only $\Delta_{0}$ and $\Delta_{2k_{F}}$ are non-
zero. We also neglect the $\beta$ quartic term of the boson energy. This goes
back to the conventional Peierls theory. Then, for a fixed $\Delta_{0}$, the
free energy as a function of $\Delta_{2k_{F}}$ is

$\mathcal{F}(\Delta_{2k_{F}})\approx-\frac{B}{2\pi}\int_{-k_{F}}^{k_{F}}\frac{dk}{2\pi}\sqrt{\Big{(}\varepsilon_{0}(k)-\mu-\frac{D\Delta_{0}}{\sqrt{L_{z}}}\Big{)}^{2}+\bigg{(}\frac{D\Delta_{2k_{F}}}{\sqrt{L_{z}}}\bigg{)}^{2}}+\frac{1}{L_{z}}|\Delta_{2k_{F}}|^{2}\
.$ (S21)

Minimizing the free energy by
$\partial\mathcal{F}(\Delta_{2k_{F}})/\partial\Delta_{2k_{F}}=0$ gives

$\Delta_{2k_{F}}\approx\frac{k_{F}^{2}}{2\sinh\big{(}\frac{2\pi^{2}k_{F}}{B}\big{)}}$
(S22)

which decreases exponentially when $B$ is small,
$\Delta_{2k_{F}}\propto\exp(-2\pi^{2}k_{F}/B)\rightarrow 0$ for sufficiently
large $k_{F}$. When $B$ is large, this CDW gap can be eventually comparable to
$\Delta_{0}$.

## III. Free energy with respect to uniform $\Delta$

Here we derive the free energy $\mathcal{F}(\Delta)$ of the dimensionless
Hamiltonian (at fixed chemical potential $\mu$) by assuming uniform boson mean
field $\Delta(z)=\Delta$, based on which we calculate the phase diagrams in
the main text by minimizing $\mathcal{F}(\Delta)$.

For constant $\Delta(z)=\Delta=\Delta_{0}/\sqrt{L_{z}}$, the dimensionless
free energy of the system is

$\mathcal{F}(\Delta,B)=\frac{B}{2\pi
L_{z}}\sum_{n\in\text{occ},k_{z}}\big{[}\epsilon_{k_{z},n}(B)-D\Delta-\mu\big{]}+\frac{1}{2}\Delta^{2}+\frac{\beta}{4}\Delta^{4}$
(S23)

where the summation in $n$ and $k_{z}$ is within the range
$\epsilon_{k_{z},n}(B)\leq D\Delta+\mu$ (i.e., occupied bands). Define
$k_{F,n}\geq 0$ as the Fermi momentum of the $n$-th LL band satisfying
$\epsilon_{k_{F,n},n}(B)=D\Delta+\mu$. When $L_{z}\rightarrow\infty$, we have

$\mathcal{F}(\Delta,B)=\frac{B}{2\pi}\sum_{n\in\text{occ}}\int_{-k_{F,n}}^{k_{F,n}}\frac{dk_{z}}{2\pi}\
\big{[}\epsilon_{k_{z},n}(B)-D\Delta-\mu\big{]}+\frac{1}{2}\Delta^{2}+\frac{\beta}{4}\Delta^{4}\
,$ (S24)

where $n\in\text{occ}$ runs over all the occupied bands (with
$\epsilon_{0,n}(B)\leq D\Delta+\mu$).

(i) For the quadratic band where
$\epsilon_{k_{z},n}(B)=\frac{1}{2}k_{z}^{2}+B(n+\frac{1}{2})$, by substituting
into variable $E=\frac{1}{2}k_{z}^{2}$ and integration by part,

$\begin{split}\mathcal{F}(\Delta,B)&=-\frac{B}{2\pi}\sum_{n\in\text{occ}}\frac{2}{2\pi}\int_{0}^{\mu+D\Delta-B(n+1/2)}\sqrt{2E}dE+\frac{1}{2}\Delta^{2}+\frac{\beta}{4}\Delta^{4}\\\
&=-\frac{\sqrt{2}}{3\pi^{2}}B\sum_{n\in\text{occ}}\bigg{[}\mu+D\Delta-B\big{(}n+\frac{1}{2}\big{)}\bigg{]}^{3/2}+\frac{1}{2}\Delta^{2}+\frac{\beta}{4}\Delta^{4}\
,\end{split}$ (S25)

and the corresponding density $n_{3D}$ is

$n_{3D}(\Delta,B)=-\frac{\partial\mathcal{F}(\Delta,B)}{\partial\mu}=\frac{\sqrt{2}}{2\pi^{2}}B\sum_{n\in\text{occ}}\bigg{[}\mu+D\Delta-B\big{(}n+\frac{1}{2}\big{)}\bigg{]}^{1/2}\
.$ (S26)

When $B\rightarrow 0$,

$\lim_{B\rightarrow
0}\mathcal{F}(\Delta,B)=-\frac{2\sqrt{2}}{15\pi^{2}}(\mu+D\Delta)^{5/2}+\frac{1}{2}\Delta^{2}+\frac{\beta}{4}\Delta^{4}\
.$ (S27)

and

$\lim_{B\rightarrow
0}n_{3D}(\Delta,B)=\frac{\sqrt{2}}{3\pi^{2}}(\mu+D\Delta)^{3/2}$ (S28)

(ii) For the Dirac band where
$\epsilon_{k_{z},n}(B)=M\sqrt{M^{2}+2Bn+k_{z}^{2}}-M^{2}$, similarly we have

$\begin{split}\mathcal{F}(\Delta,B)&=-\frac{B}{2\pi}\sum_{n\in\text{occ}}\frac{2}{2\pi}\int_{M(\sqrt{M^{2}+2Bn}-M)-D\Delta}^{\mu}\sqrt{\big{(}\frac{E}{M}+M\big{)}-M^{2}-2Bn}\
dE+\frac{1}{2}\Delta^{2}+\frac{\beta}{4}\Delta^{4}\\\
&=-\frac{B}{2\pi^{2}}\sum_{n\in\text{occ}}\int_{\sqrt{M^{2}+2Bn}}^{\frac{\mu+D\Delta}{M}+M}\sqrt{E^{\prime
2}-M^{2}-2Bn}\ dE^{\prime}+\frac{1}{2}\Delta^{2}+\frac{\beta}{4}\Delta^{4}\\\
&=-\frac{B}{4\pi^{2}}\sum_{n\in\text{occ}}(M^{2}+2Bn)\bigg{[}x\sqrt{x^{2}-1}-\ln(x+\sqrt{x^{2}-1})\bigg{]}+\frac{1}{2}\Delta^{2}+\frac{\beta}{4}\Delta^{4}\
,\\\ \end{split}$ (S29)

where

$x=\frac{\frac{\mu+D\Delta}{M}+M}{\sqrt{M^{2}+2Bn}}\ ,$ (S30)

and the 3D density is

$n_{3D}(\Delta,B)=\frac{B}{2\pi^{2}}\sum_{n\in\text{occ}}\sqrt{M^{2}+2Bn}\sqrt{x^{2}-1}\
.$ (S31)