# Giant High-order Nonlinear and Nonreciprocal Electrical Transports Induced
by Valley Flipping in Bernal Bilayer Graphene

Yuelin Shao<EMAIL_ADDRESS>Beijing National Laboratory for Condensed
Matter Physics and Institute of Physics, Chinese Academy of Sciences, Beijing
100190, China  School of Physical Sciences, University of Chinese Academy of
Sciences, Beijing 100049, China  Xi Dai<EMAIL_ADDRESS>Department of Physics,
The Hongkong University of Science and Technology, Clear Water Bay, Kowloon
999077, Hong Kong, China

###### Abstract

We investigate the electrical transport properties of the mini-valley
polarized state proposed recently in slightly doped Bernal Bilayer Graphene
(BLG) in large electric displacement fields. By minimizing the Hartree-Fock
energy functional, we first confirm the appearance of mini-valley polarized
phase. At the low carrier doping regime, the 1-pocket state will be stabilized
where only one of the trigonal-wrapping-induced Fermi pockets near the atomic-
valley center is filled. Then we study the electrical transport of the
1-pocket state by solving the Boltzmann equation. We find that the valley
polarization could be easily flopped by an in-plane electrical field, which
will lead to hysteresis loop in the direct current (DC) $I-V$ curves. Such
irreversible current responses in the DC limit will directly induce strong
nonlinear and nonreciprocal alternating current (AC) responses, which has been
already observed in the recent experiments on BLG.

##  I Introduction

In the realm of two-dimensional materials, the valley degree of freedom is a
critical factor that influences their unique magnetic, transport, and optical
characteristics. This is in addition to the well-known degrees of freedom in
condensed matter, such as spin, orbital, and momentum. In certain systems,
like twisted bilayer graphene and transition metal dichalcogenides (TMDs), a
variety of valley polarization and hybridization states have been proposed[1,
2, 3, 4, 5, 6] and observed[7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]. These
states contribute to the emergence of exotic phenomena, including the quantum
anomalous Hall effect, correlated insulators, and orbital magnetism. Recently,
valley-polarized states have also been detected in non-twisted systems such as
Rhombohedral trilayer and Bernal bilayer graphene when a displacement field is
applied along the $z$-axis[18, 19, 20, 21, 22]. In this paper, we focus on the
Bernal bilayer graphene (BLG) whose atomic structure and Brillouin zone is
shown in Fig. 1(a). The displacement field first separates the two Dirac nodes
in energy and subsequently leads to the formation of a hybridization gap due
to inter-layer coupling[23, 24, 25], resulting in a “Mexican hat” type band
dispersion for each atomic valley $\pm K$. The trigonal warping effect in BLG
further creates three local maxima in the valence band (or minima in the
conduction band) near each atomic valley center[24], as indicated by
$Q_{1,2,3}$ in Figure 1(b). As a result, the highest valence band (HVB) of
gated BLG hosts six band maxima $\pm Q_{1,2,3}$ (the prefactor $\xi=\pm$ label
the atomic valley, the subscript $\alpha=1,2,3$ label the mini-valley around
the atomic valleys), which from now on are called all together as the valleys
thought out this paper for simplicity. When the filling factor of the system
is just slightly away from the charge neutrality point, quantum oscillation
experiments indicate an emergence of the half- and quarter-metal phases[19,
20, 21]. Starting from the quarter-metal phase where only one spin and atomic
valley is occupied, if the doping concentration is further reduced so that the
mini-valley $Q_{1,2,3}$ could also be viewed as a degree of freedom, a
polarized state to the mini-valley is further predicted[26]. Such kind of
valley-polarized states break the threefold rotation symmetry $C_{3}$ and has
been identified by the $C_{3}$-breaking pattern of angle-resolved second-
harmonic nonlinear transport measurements[27]. More recently, anomalous high-
order nonlinear and nonreciprocal transport signals have also been detected in
the slightly doped BLG[28], challenging traditional transport theories and
spurring the need for theoretical investigation from new vantage points.

In this paper, we first establish the presence of mini-valley-polarized phases
in gated BLG. This is deduced from the minimization of the system’s total
energy functional under Hartree-Fock approximation. The proximity of the three
valleys around the $K$ points makes the polarized states quite susceptible to
the transitions induced by current flow, potentially leading to significant
nonlinear and non-reciprocal transport behaviors. We then analyze the current
response of systems exhibiting ”Mexican hat” type band dispersion by employing
the Boltzmann equation, incorporating the effects of long-range Coulomb
interactions among carriers. Differentiating between relaxation processes
within the same mini-valley and across different valleys allowed us to derive
equations that effectively describe the dynamics within these systems. The
resulting phenomena from our analysis include: 1) The ease of flipping mini-
valley polarization via an electric current leading to hysteresis loops in the
direct current (DC) $I-V$ profile. 2) The presence of an irreversible current
response in the DC regime leads to significant nonlinear and nonreciprocal
alternating current (AC) responses at low frequencies, where the current
response remains adiabatic to the external drive. 3) The existence of a
critical external field intensity where the long-term behavior of the current
response becomes largely insensitive to the initial polarization direction,
leading to a near-symmetrical angular response in the system.

The rest of the paper will be organized as follows. In Section II, we first
confirm the existence of the mini-valley polarized ground states in the very
low hole doping regime by minimizing the Hartree Fock mean-field energy
functional. In Section III, we solve the Boltzmann equation for the mini-
valley polarized ground state with a DC electrical field, and a valley
flipping phenomenon is identified when the electrical field strength reaches
some critical value. By replacing the DC electrical field by a low frequency
sinusoidal AC field, the AC electrical transports are investigated in Section
IV. We find that the system enters into a steady valley oscillating state if
the AC electrical field amplitude is larger that the DC critical field for
valley flipping. At this time, the current responses in time domain show
discontinuous features which lead to very high order nonlinearity in frequency
domain responses. Besides, we also find that the valley oscillating state is
insensitive to the initial polarization direction when the electrical field
amplitude is large enough, and the current responses show a
$C_{3}$-persevering angular pattern although the initial mini-valley polarized
ground state breaks the threefold rotation symmetry. In Section V, symmetry
restrictions on the angular pattern of the current responses are discussed in
more details.

##  II The Valley-Polarized Ground States

Figure 1: (a) Real space structure and Brillouin zone of BLG. (b) Constant
energy surface of the HVB in the $K$ valley for inter-layer bias potential
$U=300\mathrm{meV}$. (c)(d) Detail view of the band structure along the tow
dotted red lines in (b). (e) DOS per spin and per atomic valley. The dashed
black lines in (b-e) mark the same energy level where the hole concentration
contributed by each pocket in (b) is about $0.1\times
10^{12}\mathrm{cm}^{-2}$.

The non-interacting band structures near the atomic valley center $\xi K$
($\xi=\pm$) could be described by the following $k\cdot p$ Hamiltonian[29, 30]
under atomic basis $\\{A_{1},B_{1},A_{2},B_{2}\\}$

$H_{\xi,\bm{k}}\approx\begin{bmatrix}-U/2&v_{0}\pi^{\dagger}&-v_{4}\pi^{\dagger}&v_{3}\pi\\\
v_{0}\pi&-U/2+\Delta^{\prime}&\gamma_{1}&-v_{4}\pi^{\dagger}\\\
-v_{4}\pi&\gamma_{1}&U/2+\Delta^{\prime}&v_{0}\pi^{\dagger}\\\
v_{3}\pi^{\dagger}&-v_{4}\pi&v_{0}\pi&U/2\end{bmatrix}$ (1)

where $\pi=\hbar(\xi k_{x}+ik_{y})$ and $U$ is the inter-layer potential arose
from gating. The velocities $v_{i}$ are defined by the nearest neighbor
hopping parameters $\gamma_{i}$ as $v_{i}=\sqrt{3}a\gamma_{i}/2\hbar$ and $a$
is the lattice constant (the hopping parameters in this paper is adopted from
Kuzmenko  _et al._ [31]). In this paper, we focus on the hole doping regime of
the gated BLG and the HVB dominates the low energy excitations. Taking
$U=300\mathrm{meV}$, the constant energy contour of the HVB near $K$ point is
plotted in Fig. 1(b) by the pseudo-color map. And the band structures near the
$-K$ valley are related with $K$ valley by time reversal symmetry. The detail
band structures along the two dotted red lines in Fig. 1(b) are shown in Fig.
1(c)(d). In Fig. 1(e), we also plot the density of states (DOS) per-valley and
per-spin and the Van Hove singularities are clearly shown. Due to the
trigonal-wrapping effect from the inter-layer coupling $\gamma_{3}$ between
$A_{1}$ and $B_{2}$ atoms, the non-interacting Fermi surfaces for very low
hole doping separate into three hole pockets as illustrated by the dashed
black lines in Fig. 1(b).

In the low doping regime, the quarter-metal phase of gated BLG is verified
both experimentally[19, 20] and theoretical[32, 26, 33]. We focus on the
density region where the quarter-metal state is stable and assume the
$(\uparrow,K)$ flavor is occupied, then the low energy physics could be
described by the HVB electron in $(\uparrow,K)$ flavor with creation operator
$c^{\dagger}_{\bm{k}}$. By projecting to the occupied flavor, the many-body
Hamiltonian could be approximately written as

$H=\sum_{\bm{k}}\varepsilon_{\bm{k}}^{0}c^{\dagger}_{\bm{k}}c_{\bm{k}}+\frac{1}{2\mathcal{S}}\sum_{\bm{k}\bm{k}^{\prime}\bm{q}}U_{\bm{k}\bm{k}^{\prime}}(\bm{q})c^{\dagger}_{\bm{k}}c^{\dagger}_{\bm{k}^{\prime}}c_{\bm{k}^{\prime}+\bm{q}}c_{\bm{k}-\bm{q}},$
(2)

where $\mathcal{S}$ is the 2D area of the system. The first term is the
kinetic part and the non-interacting dispersion $\varepsilon^{0}_{\bm{k}}$ is
shown in Fig. 1(b). The second term is the interacting part, where
$U_{\bm{k}\bm{k}^{\prime}}(\bm{q})\equiv
V(\bm{q})\langle\bm{k}|\bm{k}-\bm{q}\rangle\langle\bm{k}^{\prime}|\bm{k}^{\prime}+\bm{q}\rangle$
is the interaction matrix projected to HVB and $|\bm{k}\rangle$ is the
wavefunction of HVB electron in $(\uparrow,K_{+})$ flavor. To obtain the
previous form of interaction, only the electron-electron scattering with
momentum transfer much smaller than Brillouin zone, i.e. $q\ll G$, is
considered, which is a reasonable approximation for long range interaction. In
this paper, the interaction $V(\bm{q})$ is taken as the Coulomb one
$V_{col}(q)=2\pi/\epsilon q$ and the dielectric constant is set as
$\epsilon=5$.

Under Hartree-Fock approximation, the total energy per area is expressed as a
functional of the hole occupation function $n_{\bm{k}}\equiv 1-\langle
c^{\dagger}_{\bm{k}}c_{\bm{k}}\rangle$ as
$E_{tot}[n_{\bm{k}}]=E_{0}+E_{H}+E_{K}[n_{\bm{k}}]+E_{ex}[n_{\bm{k}}]$. The
vacuum energy of the fully filled HVB
$E_{0}=\mathcal{S}^{-1}\sum_{\bm{k}}\varepsilon^{0}_{\bm{k}}$ is a constant.
Besides, we assume the system is charge uniform and the Hartree energy $E_{H}$
is also a constant. Thus these two terms will be ignored in the following
text. The other two terms are kinetic energy $E_{K}$ and exchange energy
$E_{ex}$, which are separately written as

$\displaystyle
E_{K}[n_{\bm{k}}]=\frac{1}{\mathcal{S}}\sum_{\bm{k}}(-\varepsilon^{0}_{\bm{k}})n_{\bm{k}}$
(3a) $\displaystyle
E_{ex}[n_{\bm{k}}]=-\frac{1}{2\mathcal{S}^{2}}\sum_{\bm{k}\bm{k}^{\prime}}U_{\bm{k}\bm{k}^{\prime}}(\bm{k}-\bm{k}^{\prime})n_{\bm{k}}n_{\bm{k}^{\prime}}$
(3b)

For very low doping concentration $n_{h}$ where the three hole Fermi surfaces
as shown in Fig. 1(b) do not touch, we can divide the Brillouin zone into
three regions $\mathcal{Q}_{\alpha}$ ($\alpha=1,2,3$) as illustrated in Fig.
2(b-d) and treat the mini-valley as a new degree of freedom. Under this
assumption, we further approximate the full occupation function $n_{\bm{k}}$
as three separate fermi-Dirac distributions for different mini-valleys which
only depend on the corresponding valley occupation number $n_{\alpha}$. To be
specific, the occupation function is approximated by

$n_{\bm{k}}\approx
n^{0}_{\bm{k}}(n_{\alpha[\bm{k}]})=\Theta(\mu_{\alpha[\bm{k}]}+\varepsilon^{0}_{\bm{k}}),$
(4)

where $\alpha[\bm{k}]$ is the mini-valley index of the wavevector $\bm{k}$ and
$\Theta(\epsilon)$ is the Heaviside function (Fermi-Dirac distribution at zero
temperature). In Eq. (4), the parameter $\mu_{\alpha}$ is the valley-dependent
hole chemical potential determined from the hole concentration $n_{\alpha}$ by

$n_{\alpha}=\int_{\mathcal{Q}_{\alpha}}\frac{\textrm{d}^{2}k}{(2\pi)^{2}}\;n_{\bm{k}}=\int_{\mathcal{Q}_{\alpha}}\frac{\textrm{d}^{2}k}{(2\pi)^{2}}\;\Theta(\mu_{\alpha}+\varepsilon^{0}_{\bm{k}}).$
(5)

Then we can treat the three $n_{\alpha}$s as the only variational parameters
to greatly simplify the calculations. In Fig. 2(b-d), the approximated hole
occupation function gotten by Eq. (4) (shaded areas) are compared with those
from the full variational scheme (red lines), which reveal a significant
degree of agreement.

Under this approximation, the total energy is expressed as a function of
$n_{\alpha}$ as

$E_{tot}[n_{\bm{k}}]\approx E_{tot}(n_{1},n_{2},n_{3}).$ (6)

More details of the energy functional is given in Sec. I of the supplemental
materials (SM). By minimizing Eq. (6) with respect to $n_{\alpha}$ under
constrain $n_{1}+n_{2}+n_{3}=n_{h}$, we find that the only possible stable
sates (local minimums) are the $N$-pocket states ($N=1,2,3$), where there are
$N$ hole pockets equally occupied. And this is consistent with the previous
study by Dong  _et al._ [26]. In Fig. 2(a), the total energies per hole
$E_{tot}/n_{h}$ of the $N$-pocket states are plotted in the hole density range
where they are stable, and the ground state is the one with the lowest energy.
The typical ground state hole distributions of the $N$-pocket states are also
plotted in Fig. 2(b-d).

Figure 2: (a) Total energies per hole $E_{tot}/n_{h}$ of the $N$-pocket states
in the hole density range where they are stable. At fixed hole density, the
ground state is the one with the lowest energy. (b-d) Typical hole
distributions of the $N$-pocket ground states. The shaded areas are given by
the approximation Eq. (4) and the red lines are obtained from a full
variational scheme, which reveal a significant degree of agreement. The hole
density as taken as $n_{h}=0.06,0.16,0.27\times 10^{12}\mathrm{cm}^{-2}$
respectively.

##  III DC Transport: Valley Flipping and Current Hysteresis

Figure 3: (a) Starting from the $\mathcal{Q}_{1}$-polarized state $P_{0}$ at
zero electrical field strength, the steady hole distributions are smoothly
connected if the electrical fields vary in the blue region and its boundary
indicates the critical electrical field for valley flipping. (b)(c) The same
quantities as (a) except that the initial states at zero field strength are
polarized in the $\mathcal{Q}_{2}$ and $\mathcal{Q}_{3}$ valleys respectively.
(d-f) Steady hole distributions for points $P_{0-4}$ in (a). (g-i) Steady hole
distributions for points $P_{0-4}^{\prime}$ in (b). The blue and red color
boxed panels in (d-i) form a hysteresis loop when the electrical field varies
along the blue and red lines in (a) and (b). When generating these plots, the
parameters are taken as $n_{h}=0.06\times 10^{12}\;\mathrm{cm}^{-2}$,
$\tau_{i}=10\;\mathrm{ps}$ and $\tau_{x}/\tau_{i}=100$. The field unit is
defined as $E^{*}=\hbar\sqrt{3}Q/(2e\tau_{i})\approx
19.4\mathrm{mV}/\mathrm{\mu m}$. The direction of the electrical field
indicated by the blue and red lines is $\theta=7\pi/18$.

In this section, we investigate the DC electrical transport properties of the
valley polarized states. For simplicity, we focus on the low hole density
regime in this paper where only the $1$-pocket state is stable. The Boltzmann
equation for hole distribution function is

$\displaystyle\frac{\partial n_{\bm{k}}}{\partial t}=$
$\displaystyle-\frac{e}{\hbar}\bm{E}\cdot\nabla_{\bm{k}}n_{\bm{k}}$
$\displaystyle+\frac{2\pi
n_{imp}}{\hbar}\int\frac{\textrm{d}^{2}k^{\prime}}{(2\pi)^{2}}\;|T_{\bm{k}\bm{k}^{\prime}}|^{2}\delta(\tilde{\varepsilon}_{\bm{k}}-\tilde{\varepsilon}_{\bm{k}^{\prime}})(n_{\bm{k}^{\prime}}-n_{\bm{k}}).$
(7)

The first term on the right side of the equation describes the driving effect
by a DC electrical field and the second term is the impurity scattering, where
$n_{imp}$ is the impurity concentration, $T_{\bm{k}\bm{k}^{\prime}}$ is the
impurity scattering matrix and
$\tilde{\varepsilon}_{\bm{k}}\equiv\partial_{n_{\bm{k}}}E_{tot}[n_{\bm{k}}]$
is the quasi-particle energy of the holes. Considering Eq. (3)(4), the quai-
particle energy in the $\mathcal{Q}_{\alpha}$ valley could be approximated as
$\tilde{\varepsilon}_{\bm{k}}\approx-\varepsilon^{0}_{\bm{k}}+\Sigma_{\bm{k}}$
where
$\Sigma_{\bm{k}}\approx-{\mathcal{S}}^{-1}\sum_{\bm{k}^{\prime}}U_{\bm{k}\bm{k}^{\prime}}(\bm{k}-\bm{k}^{\prime})n^{0}_{\bm{k}^{\prime}}$
is the self-energy of holes due to interaction. Near the mini-valley center,
the self-energy could be viewed as constants
$\Sigma_{\bm{k}}\approx\Sigma_{\alpha[\bm{k}]}$ for different valleys, where
$\Sigma_{\alpha}$ is obtained by
$\Sigma_{\alpha}\approx\partial_{n_{\alpha}}E_{ex}(n_{1},n_{2},n_{3})$. As for
the scattering matrix, we also assume that $|T_{\bm{k}\bm{k}^{\prime}}|$ to be
constants $T_{i}$ and $T_{x}$ for the intra- and inter-valley scattering
process, which satisfies $T_{i}\gg T_{x}$. For charge impurity with long range
impurity potential in real space, the scattering matrix drops quickly with
respect to the transferred momentum $|\bm{k}-\bm{k}^{\prime}|$, and $T_{i}\gg
T_{x}$ wil be a good approximation. Under the above assumption, the intra-
valley relaxation process is considered to be much faster than the inter-
valley process. Therefore, we can approximately separate the fast and slow
processes by ignoring the inter-valley relaxation process when dealing with
the intra-valley dynamics. Assume the hole concentration in
$\mathcal{Q}_{\alpha}$ valley is $n_{\alpha}$ at time $t$, by adopting the
relaxation time approximation, the Boltzmann equation in
$\mathcal{Q}_{\alpha}$ valley becomes

$\frac{\partial n_{\bm{k}}}{\partial
t}\approx-\frac{e}{\hbar}\bm{E}\cdot\nabla_{\bm{k}}n_{\bm{k}}-\frac{n_{\bm{k}}-n_{\bm{k}}^{0}(n_{\alpha})}{\tau_{i}},$
(8)

where $n^{0}_{\bm{k}}(n_{\alpha})$ is the Fermi-Dirac distributions given by
Eq. (4), $\tau_{i}$ is the intra-valley relaxation time[34] defined as
$\tau_{i}^{-1}\equiv 2\pi n_{imp}|T_{i}|^{2}g/\hbar$ and $g$ is the DOS at the
hole Fermi energy. For 2D system, the DOS per mini-valley is nearly a constant
and will be approximated by its value at the band edge which is about
$g\approx 0.078\;\mathrm{nm}^{-2}\mathrm{eV}^{-1}$ as shown in Fig. 1(e).
Using the Chamber’s formula[35], Eq. (8) could be analytically solved as

$n_{\bm{k}}(n_{1},n_{2},n_{3})=\frac{1}{\tau_{i}}\int_{-\infty}^{0}\textrm{d}t^{\prime}\textrm{e}^{t^{\prime}/\tau_{i}}n^{0}_{\bm{k}+e\bm{E}t^{\prime}/\hbar}(n_{\alpha[\bm{k}+e\bm{E}t^{\prime}/\hbar]})$
(9)

To further include the inter-valley process, we integrate Eq. (7) over
$\mathcal{Q}_{\alpha}$ and the Boltzmann equation could be formally written as

$\frac{\partial n_{\alpha}}{\partial t}=-\left(\frac{\partial
n_{\alpha}}{\partial t}\right)_{d}+\left(\frac{\partial n_{\alpha}}{\partial
t}\right)_{s}.$ (10)

On the right side of Eq. (10), the first term is the field driven term which
is written as

$\displaystyle\left(\frac{\partial n_{\alpha}}{\partial t}\right)_{d}=$
$\displaystyle\int_{\mathcal{Q}_{\alpha}}\frac{\textrm{d}^{2}k}{(2\pi)^{2}}\;\frac{e}{\hbar}\bm{E}\cdot\nabla_{\bm{k}}n_{\bm{k}}(n_{1},n_{2},n_{3})$
$\displaystyle=$
$\displaystyle\int_{\mathcal{Q}_{\alpha}}\frac{\textrm{d}^{2}k}{(2\pi)^{2}}\;\nabla_{\bm{k}}\cdot[\dot{\bm{k}}n_{\bm{k}}(n_{1},n_{2},n_{3})].$
(11)

Notice that $\dot{\bm{k}}=e\bm{E}/\hbar$ is the velocity of in $k$-space, we
can define the current flows in $k$-space as
$\bm{\mathcal{J}}_{\bm{k}}\equiv\dot{\bm{k}}n_{\bm{k}}$. Using Gauss’ law, Eq.
(11) could be further written as

$\left(\frac{\partial n_{\alpha}}{\partial
t}\right)_{d}=\frac{1}{(2\pi)^{2}}\oint_{\partial\mathcal{Q}_{\alpha}}\bm{\mathcal{J}}_{\bm{k}}\cdot(\textrm{d}\bm{k}\times\hat{z}).$
(12)

In other words, the field driving term could be interpreted as changing rate
of the particle numbers in area $\mathcal{Q}_{\alpha}$ caused by the $k$-space
current flow $\bm{\mathcal{J}}_{\bm{k}}$. The second term on the right side of
Eq. (10) is the changing rate in area $\mathcal{Q}_{\alpha}$ due to impurity
scattering. According to the derivation in Sec. II of the SM, this term could
be approximated by

$\displaystyle\left(\frac{\partial n_{\alpha}}{\partial t}\right)_{s}\approx$
$\displaystyle\frac{2\pi
n_{imp}|T_{x}|^{2}}{\hbar}\sum_{\beta\neq\alpha}\left[\Theta(\Sigma_{\beta}+\mu_{\beta}-\Sigma_{\alpha}-\mu_{\alpha})\right.$
$\displaystyle\left.\times
g^{2}\min(\mu_{\beta},\Sigma_{\beta}+\mu_{\beta}-\Sigma_{\alpha}-\mu_{\alpha})-\alpha\leftrightarrow\beta\right],$
(13)

As an example, we will use Fig. 3(g) to explain the physical meaning of Eq.
(13). The inter-valley scattering from the $\mathcal{Q}_{1}$ valley to
$\mathcal{Q}_{2}$ valley can only occur when
$\mu_{1}+\Sigma_{1}>\mu_{2}+\Sigma_{2}$ is satisfied. For a scattering channel
with energy $\varepsilon$, the energy should satisfy
$\Sigma_{1}<\varepsilon<\Sigma_{1}+\mu_{1}$ to make sure there are occupied
initial states in $\mathcal{Q}_{1}$ valley. Besides, the energy should also
satisfy $\varepsilon>\Sigma_{2}+\mu_{2}$ to make sure there are empty final
states in $\mathcal{Q}_{2}$ valley. Thus the energy of the scattering channel
should satisfy
$\max(\Sigma_{1},\Sigma_{2}+\mu_{2})<\varepsilon<\Sigma_{1}+\mu_{1}$ and the
total number of scattering channels are given by
$g^{2}\min(\mu_{1},\Sigma_{1}+\mu_{1}-\Sigma_{2}-\mu_{2})$. It’s convenient to
define the inter-valley relaxation time as $\tau_{x}^{-1}\equiv 2\pi
n_{imp}|T_{x}|^{2}g/\hbar$. Then the ratio between these two relaxation times
is $\tau_{x}/\tau_{i}=|T_{i}|^{2}/|T_{x}|^{2}\gg 1$, which means the intra-
valley relaxation process is much faster than the inter-valley one. By
separating the fast relaxation process in the same mini-valley and the slow
process across different valleys, the original Boltzmann equation Eq. (7) is
decoupled into the two equations Eq. (8)(10)describing the intra and inter
mini-valley dynamics respectively and will be the starting point of the
further study. From the above two equations, the steady state solution can be
obtained under particular electric field, where the two terms on the right
side of Eq. (10) cancel each other.

The 1-pocket state is threefold degenerate which could be polarized to
$\mathcal{Q}_{\alpha=1,2,3}$ valleys with equal possibility. As a consequence,
the steady solution of the Boltzmann equation ($\partial_{t}n^{S}_{\bm{k}}=0$)
is not only a function of the static electrical field, but also depends on the
initial valley polarization. For example, starting from the
$\mathcal{Q}_{1}$-polarized state as illustrated by Fig. 3(e) ($P_{0}$ point
in Fig. 3(a)), the final steady sates are evolved smoothly if the electrical
fields only vary in the blue region in Fig. 3(a) (the radial coordinate
represents electrical field strength and the angular coordinate represents the
direction). When the electrical field is applied along the
$\overrightarrow{P_{0}P_{1}P_{2}}$ direction in Fig. 3(a), the holes will
accelerate as $\hbar\dot{\bm{k}}=e\bm{E}$. Such movements will be balanced by
the intra-valley relaxation process, which results into a shifted steady hole
distributions as shown in Fig. 3(d) just like what happens in normal metals.
When the electrical field reverses to the $\overrightarrow{P_{0}P_{3}}$
direction before reaching the boundary of the blue region in Fig. 3(a), holes
will be driven into the $\mathcal{Q}_{2}$ valley which is balanced by the
inter-valley relaxation process and the typical steady hole distributions are
shown in Fig. 3(f). With the increase of the electrical field strength, the
steady hole concentration in $\mathcal{Q}_{2}$ increases and there must be a
critical point ($P_{4}$ in Fig. 3(a)) where the inter-valley relaxation can
not balance the field driving effect and a valley flipping occurs. After that,
the steady state will become the $\mathcal{Q}_{2}$-polarized state illustrated
by Fig. 3(h) ($P_{0}^{\prime}$ point in Fig. 3(b)). The cases for electrical
fields varying along the red line in Fig. 3(b) are similar and the
corresponding steady hole distributions for points $P_{0-4}^{\prime}$ are
shown in Fig. 3(g-i). Due to the valley flipping phenomenon, the steady hole
distribution will form a hysteresis loop when electrical field strength takes
values along the blue and red lines in Fig. 3(a)(b). Between red and blue
boxed panels in Fig. 3(d-i), the blue and red arrows indicate smooth
reversible transitions while the black dashed arrow lines means irreversible.

Figure 4: (a)(b)(c)(d) Hysteresis curves for $p^{S}$, $j_{N,\parallel}$ and
$j_{N,\perp}$ and $j_{H}$ as functions of the electric field strength $E$
along the blue and red lines in Fig. 3(a)(b). The critical field strengths for
the valley flipping are about $E_{c}^{1}/E^{*}\approx 0.1693$ and
$E_{c}^{2}/E^{*}\approx 0.1691$. The red and blue solid lines are reversible
while the black dashed arrow lines are one-way through. The direction of the
electrical field and the different current components are also indicated in
the insets respectively.

To see the hysteretic effect more clearly, we define the polarization strength
of the steady hole concentrations between $\mathcal{Q}_{1}$ and
$\mathcal{Q}_{2}$ valleys as $p^{S}\equiv(n^{S}_{2}-n^{S}_{1})/n_{h}$. Then
the hysteresis curve of $p^{S}$ along the blue and red lines in Fig. 3(a)(b)
is plotted in Fig. 4(a). The points $P_{0-4}$ and $P^{\prime}_{1-4}$ are also
marked in Fig. 4(a) accordingly. The hysteresis of the distribution functions
and hole concentrations will directly lead to the current hysteresis. With the
steady hole distributions $n_{\bm{k}}^{S}$ solved from the equations Eq.
(8)(10), the steady electrical current density is calculated as

$\displaystyle\bm{j}(n^{S}_{\bm{k}})=$
$\displaystyle\frac{e}{\hbar}\int\frac{\textrm{d}^{2}k}{(2\pi)^{2}}\;(\nabla_{\bm{k}}\tilde{\varepsilon}_{\bm{k}}+e\bm{E}\times\hat{z}\Omega_{\bm{k}})n^{S}_{\bm{k}}$
$\displaystyle\approx$
$\displaystyle\frac{e}{\hbar}\int\frac{\textrm{d}^{2}k}{(2\pi)^{2}}\;(-\nabla_{\bm{k}}{\varepsilon}^{0}_{\bm{k}}+e\bm{E}\times\hat{z}\Omega_{\bm{k}})n^{S}_{\bm{k}}.$
(14)

In the above expression, the first term is the normal current $\bm{j}_{N}$ and
the second term is the anomalous Hall current $\bm{j}_{H}$ due to the non-zero
Berry curvature $\Omega_{\bm{k}}$ of HVB[36, 37]. Decomposed to the
longitudinal direction of the electrical field $\hat{r}$ and transverse
direction $\hat{\theta}\equiv\hat{z}\times\hat{r}$, the normal current has
both longitudinal and transverse components
$\bm{j}_{N}=j_{N,\parallel}\hat{r}+j_{N,\perp}\hat{\theta}$, while the Hall
current is always along the transverse direction
$\bm{j}_{H}=j_{H}\hat{\theta}$. In Fig. 4(b-d), the hysteresis curves for the
different current components $j_{N,\parallel}$, $j_{N,\perp}$ and $j_{H}$ are
also plotted.

In Fig. 3 and 4, the hole density is taken as $n_{h}=0.06\times
10^{12}\;\mathrm{cm}^{-2}$ and the relaxation times are taken as
$\tau_{i}=10\;\mathrm{ps}$ and $\tau_{x}/\tau_{i}=100$. The electrical field
unit is defined as $E^{*}=\hbar\sqrt{3}Q/(2e\tau_{i})\approx
19.4\mathrm{mV}/\mathrm{\mu m}$, where $Q=|\bm{Q}_{\alpha}-\bm{K}_{+}|$. The
current unit is defined by
$j^{*}=e^{2}\tau_{i}n_{h}E^{*}/\sqrt{m_{r}m_{\varphi}}\approx
3.33\;\mathrm{nA}/\mathrm{nm}$ where $m_{r}$ and $m_{\varphi}$ are effective
masses at the mini-valley center in radial and angular direction respectively.

##  IV AC Transport: Valley Oscillation And High-Order Nonlinearity

In this section, we replace the DC field discussed in last section by a
sinusoidal AC field
$\tilde{\bm{E}}(t)=\tilde{E}(\hat{x}\cos\theta+\hat{y}\sin\theta)\sin(2\pi
t/T)$, which is more close to the realistic experimental setups. We assume
that the period $T$ of the external field is much longer than the inter-valley
relaxation time $\tau_{x}$, which is always true for the experimental
measurements below MHz. Then an adiabatic approximation could be safely
implemented.

Figure 5: (a) Current responses in the time domain. The inset shows the
direction of the AC electrical field and different current components. The
longitudinal and transverse normal currents $j_{N,\parallel}$, $j_{N,\perp}$
(left axis) are encoded by blue and red colors respectively, while the Hall
current $j_{H}$ (right axis) is represented by gray color. The field amplitude
is taken as $\tilde{E}/E^{*}=0.2$, and the steady state is an oscillating
state between $\mathcal{Q}_{1}$ and $\mathcal{Q}_{2}$ valleys. (b) The
corresponding current responses in the frequency domain. The results for even
($n=2m$) and odd ($n=2m+1$) frequencies are represented by dots and circles
respectively. Both the horizontal and vertical axes are drawn on logarithmic
scales.

From the hysteresis curve in Fig. 4(a), it’s easy to see that the final state
is an oscillating state between $\mathcal{Q}_{1}$ and $\mathcal{Q}_{2}$
valleys if the AC electrical field amplitudes satisfies
$\tilde{E}>E_{c}^{1,2}$. Taking $\tilde{E}/E^{*}=0.2$ and $\theta=7\pi/18$,
the current responses during one period are directly got by tracing the
hysteresis curves in Fig. 4(b-d) and plotted in Fig. 5(a). Due to the valley
flipping phenomenon, the current responses in the time domain show
discontinuous features when the electrical field passes the critical values.
By defining the current response at the $n$-th harmonic frequency as

$j^{(n)}\equiv\frac{1}{T}\int_{0}^{T}\textrm{d}t\;j(t)\textrm{e}^{in\omega
t},$ (15)

the current responses in the frequency domain are also calculated and plotted
in Fig. 5(b). As a result of the discontinuities in the time domain, the
current responses in the frequency domain decay in a power-law manner, i.e.
$\log j^{(n)}\sim-a\log n+b$ and $j^{(n)}\sim n^{-a}$. This power-law decay
leads to a very high-order nonlinearity, which is beyond the perturbation
description. For the perfect step-like discontinuities shown in Fig. 5(a), the
power exponent is fitted as $a=1$. While in realistic experiments, a larger
exponent $a>1$ is expected, since the discontinuities are always smoothed out
by inhomogeneity of the sample.

Figure 6: The line plots represent the critical AC electrical field amplitudes
$\tilde{E}_{c}^{\alpha}$ for reaching a valley-oscillating state, where the
superscripts $\alpha=1,2,3$ indicate the initial valley polarization
$\mathcal{Q}_{\alpha}$. The filled regions represent the valley-oscillating
states which are independent of the initial polarization. For example, the
purple colored region labeled by
$\mathcal{Q}_{1}\leftrightarrow\mathcal{Q}_{2}$ means that the final state is
an oscillating state between $\mathcal{Q}_{1}$ and $\mathcal{Q}_{2}$ valleys
when the electrical field amplitude is taken in this region. The maximum value
of all the critical field amplitude is
$\tilde{E}_{M}\equiv\max(\tilde{E}^{\alpha}_{c})\approx 0.3177E^{*}$ and the
minimum value is $\tilde{E}_{m}\equiv\min(\tilde{E}^{\alpha}_{c})\approx
0.1589E^{*}$

If Fig. 6, we also plot the critical amplitudes $\tilde{E}_{c}^{\alpha}$ for
reaching the valley-oscillating states (the superscripts $\alpha=1,2,3$
indicate the initial valley polarization $\mathcal{Q}_{\alpha}$). An
interestingly condition is when the the electrical amplitude is larger than
all of the three critical amplitudes $\tilde{E}>\tilde{E}_{c}^{\alpha=1,2,3}$.
If this is satisfied, the final valley-oscillating state will be independent
of the initial polarization. Using the oscillating state between
$\mathcal{Q}_{1}$ and $\mathcal{Q}_{2}$ for example (the purple region labeled
by $\mathcal{Q}_{1}\leftrightarrow\mathcal{Q}_{2}$ in Fig. 6), it’s not
surprising to reach such an oscillating state if the initial state is
polarized in $\mathcal{Q}_{1}$ or $\mathcal{Q}_{2}$ valleys as indicated by
the hysteresis loop Fig. 3(d-i) and Fig. 4(a). If the initial state is
polarized in the $\mathcal{Q}_{3}$ valley, the system will be first flipped
into the $\mathcal{Q}_{1}$-polarized state if the electrical field is applied
along the line
$\overrightarrow{P_{0}^{\prime\prime}P_{1}^{\prime\prime}P_{2}}$ in Fig. 3(c).
After that, the system also enters into a steady oscillating state between the
$\mathcal{Q}_{1}$ and $\mathcal{Q}_{2}$ valleys.

The critical amplitudes in Fig. 6 have both maximum value
$\tilde{E}_{M}=\max(\tilde{E}_{c}^{\alpha})\approx 0.3177E^{*}$ and minimum
value $\tilde{E}_{m}=\min(\tilde{E}_{c}^{\alpha})\approx 0.1589E^{*}$ and the
angular patterns of the current responses will behave differently for
$\tilde{E}<\tilde{E}_{m}$ and $\tilde{E}>\tilde{E}_{M}$. By assuming a
$\mathcal{Q}_{1}$-polarized initial state, the second harmonic currents
$j^{(2)}$ of the final periodic steady state are calculated and their norms
$|j^{(2)}|$ are plotted in Fig. 7 as functions of $\theta$ in polar
coordinates. Different rows mean different AC electrical field amplitudes as
indicated in the left axis while different columns represent different current
components as indicated in the top axis. These plots are colored red and blue
according to the sign of their real part $\mathrm{Re}(j^{(2)})$. All the plots
in Fig. 7 are antisymmetric about the inversion operation, this arises from
the additional tempo-spatial symmetry of the sinusoidal AC electrical field
$\tilde{\bm{E}}(t;\theta+\pi)=\tilde{\bm{E}}(t+T/2;\theta)$ as discussed in
Sec. IV of the SM. In the first row, the field amplitude is taken as
$\tilde{E}=0.05E^{*}<\tilde{E}_{m}$ and no valley flipping occurs. Then the
current responses show $C_{3}$-breaking patterns as has already been observed
by Lin  _et al._ [27]. In the third row, a much larger electrical field
amplitudes $\tilde{E}=0.35E^{*}>\tilde{E}_{M}$ is taken and the broken $C_{3}$
symmetry recovers. This is not surprising since the final states are valley-
oscillating states in any electrical field direction as indicated by Fig. 6
and the system will lose its memory about the initial polarization. As
discussed in Sec. III of the SM, the critical field strength for valley
flipping decreases with hole concentration $n_{h}$, suggesting the
$C_{3}$-persevering patterns are more likely to appear near the phase boundary
of the 1-pocket state, which was also observed by Lin  _et al._ [27] but might
be wrongly interpreted as an unpolarized state in their paper. At the
intermediate field amplitudes, for example $\tilde{E}=0.2E^{*}$ in the second
row, the angular patterns are similar to the first row in the directions with
no valley flipping and similar to the third row in the directions where valley
oscillation occurs.

Figure 7: Current responses at the second harmonic frequency $j^{(2)}$ as
functions of electrical field direction $\theta$ in polar coordinates.
According to the sign of $\mathrm{Re}(j^{(2)})$, the plots are colored red and
blue respectively. The initial valley polarization is also marked by the
$\mathcal{Q}_{1}$ labels. Different rows mean different AC electrical field
amplitudes while different columns represent different current components.

##  V Discussion

When the spin and atomic valley degeneracy is lifted in the quarter-metal
phase, the inversion and time reversal symmetries are both broken, which
ensures the existence of the anomalous Hall currents as illustrated by the
gray line in Fig. 5(a). At this time, the approximated point group is $D_{3d}$
which contains the threefold rotations $C_{3}$ and mirror symmetries
$M_{\theta}$ with respect to the mirror lines $\theta=0,\;\pm\pi/3$. If the
system is further condensed to the 1-pocket state, for example the
$\mathcal{Q}_{1}$-polarized state, the rotation symmetries are broken and only
one mirror symmetry $M_{\pi/3}$ survives. As a consequence of the low
symmetry, all the three electrical current components $j_{N,\parallel}$,
$j_{N,\perp}$ and $j_{H}$ usually exist at the same time and they have both
even and odd frequency responses in the frequency domain as shown in Fig.
5(b). To distinguish the two kind of currents in the transverse direction, the
surviving mirror symmetry could be utilized. As discussed in Sec. IV of the
SM, the mirror symmetry of the $\mathcal{Q}_{1}$-polarized state requires the
angular patterns of the current responses to satisfy

$\displaystyle j_{N,\parallel}(2\pi/3-\theta)=j_{N,\parallel}(\theta)$ (16a)
$\displaystyle j_{N,\perp}(2\pi/3-\theta)=-j_{N,\perp}(\theta)$ (16b)
$\displaystyle j_{H}(2\pi/3-\theta)=j_{H}(\theta).$ (16c)

And the odd and even parities of the normal current $j_{N,\perp}$ and Hall
current $j_{H}$ under mirror operation are clearly shown in Fig. 7 by the last
two columns.

Another concerning problem is spatial domains where the holes populate
different mini-valleys, which make the system nonuniform. However, as shown in
our paper, the valley polarizations in different domains could be aligned by a
large enough electrical field or equivalently current flow. And the discussion
of the uniform case in our paper is sufficient to describe realistic
experiments.

###### Acknowledgements.

We thank Prof. Guangyu Zhang for their inspiring experimental observations and
helpful discussions. This work was fully supported by a fellowship award from
the Research Grants Council of the Hong Kong Special Administrative Region,
China (Project No. HKUST SRFS2324-6S01)

## References

  * Liu  _et al._ [2019] J. Liu, Z. Ma, J. Gao, and X. Dai, Quantum Valley Hall Effect, Orbital Magnetism, and Anomalous Hall Effect in Twisted Multilayer Graphene Systems, Physical Review X 9, 031021 (2019).
  * Kang and Vafek [2019] J. Kang and O. Vafek, Strong Coupling Phases of Partially Filled Twisted Bilayer Graphene Narrow Bands, Physical Review Letters 122, 246401 (2019).
  * Zhang  _et al._ [2020] Y. Zhang, K. Jiang, Z. Wang, and F. Zhang, Correlated insulating phases of twisted bilayer graphene at commensurate filling fractions: A Hartree-Fock study, Physical Review B 102, 035136 (2020).
  * Xie and MacDonald [2020] M. Xie and A. H. MacDonald, Nature of the Correlated Insulator States in Twisted Bilayer Graphene, Physical Review Letters 124, 097601 (2020).
  * Liu and Dai [2021] J. Liu and X. Dai, Theories for the correlated insulating states and quantum anomalous Hall effect phenomena in twisted bilayer graphene, Physical Review B 103, 035427 (2021).
  * Devakul  _et al._ [2021] T. Devakul, V. Crépel, Y. Zhang, and L. Fu, Magic in twisted transition metal dichalcogenide bilayers, Nature Communications 12, 6730 (2021).
  * Cao  _et al._ [2018] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature 556, 80 (2018).
  * Sharpe  _et al._ [2019] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene, Science 365, 605 (2019).
  * Jiang  _et al._ [2019] Y. Jiang, X. Lai, K. Watanabe, T. Taniguchi, K. Haule, J. Mao, and E. Y. Andrei, Charge order and broken rotational symmetry in magic-angle twisted bilayer graphene, Nature 573, 91 (2019).
  * Choi  _et al._ [2019] Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora, R. Polski, Y. Zhang, H. Ren, J. Alicea, G. Refael, F. von Oppen, K. Watanabe, T. Taniguchi, and S. Nadj-Perge, Electronic correlations in twisted bilayer graphene near the magic angle, Nature Physics 15, 1174 (2019).
  * Cao  _et al._ [2020] Y. Cao, D. Rodan-Legrain, O. Rubies-Bigorda, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Tunable correlated states and spin-polarized phases in twisted bilayer–bilayer graphene, Nature 583, 215 (2020).
  * Serlin  _et al._ [2020] M. Serlin, C. L. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. F. Young, Intrinsic quantized anomalous Hall effect in a moiré heterostructure, Science 367, 900 (2020).
  * Chen  _et al._ [2020] G. Chen, A. L. Sharpe, E. J. Fox, Y.-H. Zhang, S. Wang, L. Jiang, B. Lyu, H. Li, K. Watanabe, T. Taniguchi, Z. Shi, T. Senthil, D. Goldhaber-Gordon, Y. Zhang, and F. Wang, Tunable correlated Chern insulator and ferromagnetism in a moiré superlattice, Nature 579, 56 (2020).
  * Polshyn  _et al._ [2020] H. Polshyn, J. Zhu, M. A. Kumar, Y. Zhang, F. Yang, C. L. Tschirhart, M. Serlin, K. Watanabe, T. Taniguchi, A. H. MacDonald, and A. F. Young, Electrical switching of magnetic order in an orbital Chern insulator, Nature 588, 66 (2020).
  * Tschirhart  _et al._ [2021] C. L. Tschirhart, M. Serlin, H. Polshyn, A. Shragai, Z. Xia, J. Zhu, Y. Zhang, K. Watanabe, T. Taniguchi, M. E. Huber, and A. F. Young, Imaging orbital ferromagnetism in a moiré Chern insulator, Science 372, 1323 (2021).
  * Choi  _et al._ [2021] Y. Choi, H. Kim, Y. Peng, A. Thomson, C. Lewandowski, R. Polski, Y. Zhang, H. S. Arora, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, Correlation-driven topological phases in magic-angle twisted bilayer graphene, Nature 589, 536 (2021).
  * Pierce  _et al._ [2021] A. T. Pierce, Y. Xie, J. M. Park, E. Khalaf, S. H. Lee, Y. Cao, D. E. Parker, P. R. Forrester, S. Chen, K. Watanabe, T. Taniguchi, A. Vishwanath, P. Jarillo-Herrero, and A. Yacoby, Unconventional sequence of correlated Chern insulators in magic-angle twisted bilayer graphene, Nature Physics 17, 1210 (2021).
  * Zhou  _et al._ [2021] H. Zhou, T. Xie, A. Ghazaryan, T. Holder, J. R. Ehrets, E. M. Spanton, T. Taniguchi, K. Watanabe, E. Berg, M. Serbyn, and A. F. Young, Half- and quarter-metals in rhombohedral trilayer graphene, Nature 598, 429 (2021).
  * de la Barrera  _et al._ [2022] S. C. de la Barrera, S. Aronson, Z. Zheng, K. Watanabe, T. Taniguchi, Q. Ma, P. Jarillo-Herrero, and R. Ashoori, Cascade of isospin phase transitions in Bernal-stacked bilayer graphene at zero magnetic field, Nature Physics 18, 771 (2022).
  * Seiler  _et al._ [2022] A. M. Seiler, F. R. Geisenhof, F. Winterer, K. Watanabe, T. Taniguchi, T. Xu, F. Zhang, and R. T. Weitz, Quantum cascade of correlated phases in trigonally warped bilayer graphene, Nature 608, 298 (2022).
  * Zhou  _et al._ [2022] H. Zhou, L. Holleis, Y. Saito, L. Cohen, W. Huynh, C. L. Patterson, F. Yang, T. Taniguchi, K. Watanabe, and A. F. Young, Isospin magnetism and spin-polarized superconductivity in Bernal bilayer graphene, Science 375, 774 (2022).
  * Han  _et al._ [2023] X. Han, Q. Liu, R. Niu, Z. Qu, Z. Wang, Z. Li, C. Han, K. Watanabe, T. Taniguchi, Z. Gan, and J. Lu, Spontaneous isospin polarization and quantum Hall ferromagnetism in a rhombohedral trilayer graphene superlattice, Chinese Physics B 32, 117201 (2023).
  * McCann [2006] E. McCann, Asymmetry gap in the electronic band structure of bilayer graphene, Physical Review B 74, 161403 (2006).
  * McCann and Fal’ko [2006] E. McCann and V. I. Fal’ko, Landau-Level Degeneracy and Quantum Hall Effect in a Graphite Bilayer, Physical Review Letters 96, 086805 (2006).
  * Oostinga  _et al._ [2008] J. B. Oostinga, H. B. Heersche, X. Liu, A. F. Morpurgo, and L. M. K. Vandersypen, Gate-induced insulating state in bilayer graphene devices, Nature Materials 7, 151 (2008).
  * Dong  _et al._ [2023] Z. Dong, M. Davydova, O. Ogunnaike, and L. Levitov, Isospin- and momentum-polarized orders in bilayer graphene, Physical Review B 107, 075108 (2023).
  * Lin  _et al._ [2023] J.-X. Lin, Y. Wang, N. J. Zhang, K. Watanabe, T. Taniguchi, L. Fu, and J. I. A. Li, Spontaneous momentum polarization and diodicity in Bernal bilayer graphene (2023), arxiv:2302.04261 [cond-mat] .
  * Zan  _et al._ [2024] X. Zan, W. Yang, L. Du, and G. Zhang, Private Communication (2024).
  * McCann and Koshino [2013] E. McCann and M. Koshino, The electronic properties of bilayer graphene, Reports on Progress in Physics 76, 056503 (2013).
  * Jung and MacDonald [2014] J. Jung and A. H. MacDonald, Accurate tight-binding models for the $\pi$ bands of bilayer graphene, Physical Review B 89, 035405 (2014).
  * Kuzmenko  _et al._ [2009] A. B. Kuzmenko, I. Crassee, D. van der Marel, P. Blake, and K. S. Novoselov, Determination of the gate-tunable band gap and tight-binding parameters in bilayer graphene using infrared spectroscopy, Physical Review B 80, 165406 (2009).
  * Xie and Das Sarma [2023] M. Xie and S. Das Sarma, Flavor symmetry breaking in spin-orbit coupled bilayer graphene, Physical Review B 107, L201119 (2023).
  * Jang  _et al._ [2023] Y. Jang, Y. Park, J. Jung, and H. Min, Chirality and correlations in the spontaneous spin-valley polarization of rhombohedral multilayer graphene, Physical Review B 108, L041101 (2023).
  * Mahan [2000] G. D. Mahan,  _Many-Particle Physics_ (Springer Science & Business Media, 2000).
  * Budd [1962] H. Budd, Chamber’s Solution of the Boltzmann Equation, Physical Review 127, 4 (1962).
  * Sundaram and Niu [1999] G. Sundaram and Q. Niu, Wave-packet dynamics in slowly perturbed crystals: Gradient corrections and Berry-phase effects, Physical Review B 59, 14915 (1999).
  * Xiao  _et al._ [2010] D. Xiao, M.-C. Chang, and Q. Niu, Berry phase effects on electronic properties, Reviews of Modern Physics 82, 1959 (2010).

Supplemental Materials: Giant High-order Nonlinear and Nonreciprocal
Electrical Transports Induced by Valley Flipping in Bernal Bilayer Graphene

##  I The Exchange Energy Functional

To minimize $E_{tot}(n_{1},n_{2},n_{3})$ in the main text, the exchange energy
functional is further approximated by its series expansions with respect to
$n_{\alpha=1,2,3}$. For convenience, we define the intra-valley exchange
energy $E_{ex}^{i}$ and inter-valley exchange $E_{ex}^{x}$ as

$\displaystyle E_{ex}^{i}=$
$\displaystyle-\frac{1}{2\mathcal{S}^{2}}\sum_{\alpha}\sum_{\bm{k}\bm{k}^{\prime}\in\mathcal{Q}_{\alpha}}U_{\bm{k}\bm{k}^{\prime}}(\bm{k}-\bm{k}^{\prime})n_{\bm{k}}n_{\bm{k}^{\prime}}$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{\alpha}\int\frac{\textrm{d}^{2}\delta
k\textrm{d}^{2}\delta
k^{\prime}}{(2\pi)^{4}}\;|\langle\bm{Q}_{\alpha}+\delta\bm{k}|\bm{Q}_{\alpha}+\delta\bm{k}^{\prime}\rangle|^{2}V_{col}(\delta\bm{k}-\delta\bm{k}^{\prime})n^{0}_{\alpha,\bm{Q}_{\alpha}+\delta\bm{k}}n^{0}_{\alpha,\bm{Q}_{\alpha}+\delta\bm{k}^{\prime}}$
(S1) $\displaystyle E_{ex}^{x}=$
$\displaystyle-\frac{1}{2\mathcal{S}^{2}}\sum_{\alpha\neq\beta}\sum_{\bm{k}\in\mathcal{Q}_{\alpha},\bm{k}^{\prime}\in\mathcal{Q}_{\beta}}U_{\bm{k}\bm{k}^{\prime}}(\bm{k}-\bm{k}^{\prime})n_{\bm{k}}n_{\bm{k}^{\prime}}$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{\alpha\neq\beta}\int\frac{\textrm{d}^{2}\delta
k\textrm{d}^{2}\delta
k^{\prime}}{(2\pi)^{4}}\;|\langle\bm{Q}_{\alpha}+\delta\bm{k}|\bm{Q}_{\beta}+\delta\bm{k}^{\prime}\rangle|^{2}V_{col}(\delta\bm{k}-\delta\bm{k}^{\prime}+\bm{Q}_{\alpha}-\bm{Q}_{\beta})n^{0}_{\alpha,\bm{Q}_{\alpha}+\delta\bm{k}}n^{0}_{\beta,\bm{Q}_{\beta}+\delta\bm{k}^{\prime}}$
(S2)

where $\delta k\ll Q_{\alpha}$ is measured from the mini-valley center
$\bm{Q}_{\alpha}$. Due to the $1/q$ divergence of the Coulomb interaction
$V_{col}(q)$ in Eq. (S1), the lowest order term when expanding with respect to
$n_{\alpha}$ is $n_{\alpha}^{3/2}$, which is similar to the conventional 2D
electron gas. Then up to $n_{h}^{3}$, the intra-valley exchange energy is
expanded as

$E_{ex}^{i}\approx-\sum_{\alpha}(A_{1}n_{\alpha}^{3/2}+A_{2}n_{\alpha}^{5/2})+\mathcal{O}(n_{h}^{7/2})$
(S3)

As for the intra-valley term Eq. (S2), since $\delta k\ll Q_{\alpha}$ is
satisfied in the low doping regime, and
$\delta\bm{k}-\delta\bm{k}^{\prime}+\bm{Q}_{\alpha}-\bm{Q}_{\beta}$ can never
be zero, there will be no singularities in the integration and the expansion
of $E_{ex}^{x}$ with respect to $n_{\alpha}$ is the conventional Taylor
expansion. Then up to $n_{h}^{3}$, the inter-valley exchange energy should
have the expansion form

$E_{ex}^{x}\approx-\sum_{\alpha\neq\beta}(B_{1}n_{\alpha}n_{\beta}+B_{2}n_{\alpha}^{2}n_{\beta})+\mathcal{O}(n_{h}^{4})$
(S4)

By fitting the numerical results, the expansion coefficients are gotten and
listed in Tab. 1

Table 1: Expansion coefficients of the exchange energy functional fitted from numerical results. $A_{1}$ | $A_{2}$ | $B_{1}$ | $B_{2}$
---|---|---|---
$0.429\;\mathrm{eVnm}$ | $-16.0\;\mathrm{eVnm^{3}}$ | $0.366\;\mathrm{eVnm^{2}}$ | $41.1\;\mathrm{eVnm^{4}}$

##  II The Dynamics Equation of Valley Concentrations

When integrate over the $\mathcal{Q}_{\alpha}$ valley, the Boltzmann equation
becomes

$\displaystyle\frac{\partial n_{\alpha}}{\partial t}=$
$\displaystyle-\int_{\mathcal{Q}_{\alpha}}\frac{\textrm{d}^{2}k}{(2\pi)^{2}}\;\frac{e}{\hbar}\bm{E}\cdot\nabla_{\bm{k}}n_{\bm{k}}+\frac{2\pi
n_{imp}}{\hbar}\int_{\mathcal{Q}_{\alpha}}\frac{\textrm{d}^{2}k}{(2\pi)^{2}}\int\frac{\textrm{d}^{2}k^{\prime}}{(2\pi)^{2}}\;|T_{\bm{k}\bm{k}^{\prime}}|^{2}\delta(\tilde{\varepsilon}_{\bm{k}}-\tilde{\varepsilon}_{\bm{k}^{\prime}})(n_{\bm{k}^{\prime}}-n_{\bm{k}})$
$\displaystyle=$
$\displaystyle-\int_{\mathcal{Q}_{\alpha}}\frac{\textrm{d}^{2}k}{(2\pi)^{2}}\;\frac{e}{\hbar}\bm{E}\cdot\nabla_{\bm{k}}n_{\bm{k}}+\frac{2\pi
n_{imp}|T_{i}|^{2}}{\hbar}\int_{\mathcal{Q}_{\alpha}}\frac{\textrm{d}^{2}k}{(2\pi)^{2}}\int_{\mathcal{Q}_{\alpha}}\frac{\textrm{d}^{2}k^{\prime}}{(2\pi)^{2}}\;\delta(\tilde{\varepsilon}_{\bm{k}}-\tilde{\varepsilon}_{\bm{k}^{\prime}})(n_{\bm{k}^{\prime}}-n_{\bm{k}})$
$\displaystyle+\frac{2\pi
n_{imp}|T_{x}|^{2}}{\hbar}\int_{\mathcal{Q}_{\alpha}}\frac{\textrm{d}^{2}k}{(2\pi)^{2}}\int_{\bar{\mathcal{Q}}_{\alpha}}\frac{\textrm{d}^{2}k^{\prime}}{(2\pi)^{2}}\;\delta(\tilde{\varepsilon}_{\bm{k}}-\tilde{\varepsilon}_{\bm{k}^{\prime}})(n_{\bm{k}^{\prime}}-n_{\bm{k}})$
(S5)

where $\bar{\mathcal{Q}}_{\alpha}$ means regions exclude
$\mathcal{Q}_{\alpha}$. The three terms describe the field driving, intra- and
inter-valley scattering processes respectively. The intra-valley term change
sign when switching the indices $\bm{k}$ and $\bm{k}^{\prime}$ which means the
integral is zero. This reflects the fact that the intra-valley scattering
process will not effect the hole concentrations in each valley. To further
simplify the inter-valley term, we approximate the distribution $n_{\bm{k}}$
in the integrand by $n^{0}_{\bm{k}}$ and the last term becomes

$\displaystyle\frac{2\pi
n_{imp}|T_{x}|^{2}}{\hbar}\int_{\mathcal{Q}_{\alpha}}\frac{\textrm{d}^{2}k}{(2\pi)^{2}}\int_{\bar{\mathcal{Q}}_{\alpha}}\frac{\textrm{d}^{2}k^{\prime}}{(2\pi)^{2}}\;\delta(\tilde{\varepsilon}_{\bm{k}}-\tilde{\varepsilon}_{\bm{k}^{\prime}})(n^{0}_{\bm{k}^{\prime}}-n^{0}_{\bm{k}})$
$\displaystyle\approx$ $\displaystyle\frac{2\pi
n_{imp}|T_{x}|^{2}g^{2}}{\hbar}\sum_{\beta\neq\alpha}\int\textrm{d}\varepsilon\int\textrm{d}\varepsilon^{\prime}\;\delta(\varepsilon+\Sigma_{\alpha}-\varepsilon^{\prime}-\Sigma_{\beta})[\Theta(\mu_{\beta}-\varepsilon^{\prime})-\Theta(\mu_{\alpha}-\varepsilon)]$
$\displaystyle=$ $\displaystyle\frac{2\pi
n_{imp}|T_{x}|^{2}g^{2}}{\hbar}\sum_{\beta\neq\alpha}\left[\Theta(\Sigma_{\beta}+\mu_{\beta}-\Sigma_{\alpha}-\mu_{\alpha})\min(\mu_{\beta},\Sigma_{\beta}+\mu_{\beta}-\Sigma_{\alpha}-\mu_{\alpha})-\alpha\leftrightarrow\beta\right].$
(S6)

##  III Dependence of the Critical Field Strength on Hole Concentration and
Relaxation Times

In this part, we investigate the dependence of the critical field strength for
valley flipping on hole concentration $n_{h}$ and relaxation times
$\tau_{i,x}$. Without losing generality, we will assume that the electrical
field is applied in the direction $\pi/3<\theta<\pi$ and the initial state is
$\mathcal{Q}_{1}$-polarized. Other cases could by easily got by symmetry
operations.

we first consider the zero-doping limit $n_{h}\to 0$. At this time, the
dynamics equation of valley concentrations could be much simplified and
analytically solved. In the zero doping limit, an electrical field in the
direction $\pi/3<\theta<\pi$ can only drive holes from the $\mathcal{Q}_{1}$
to the $\mathcal{Q}_{2}$ valley and we could define the driving rate as

$\displaystyle\eta\equiv$ $\displaystyle\lim_{n_{h}\to
0}-\frac{1}{n_{1}}\int_{\mathcal{Q}_{2}}\frac{\textrm{d}^{2}k}{(2\pi)^{2}}\;\frac{e}{\hbar}\bm{E}\cdot\nabla_{\bm{k}}n_{\bm{k}}$
$\displaystyle=$ $\displaystyle\lim_{n_{h}\to
0}-\frac{1}{n_{1}}\int_{\mathcal{Q}_{2}}\frac{\textrm{d}^{2}k}{(2\pi)^{2}}\;\int_{-\infty}^{0}\textrm{d}t^{\prime}\;\frac{\textrm{e}^{t^{\prime}/\tau_{i}}}{\tau_{i}}\frac{e}{\hbar}\bm{E}\cdot\nabla_{\bm{k}}n^{0}_{1,\bm{k}+e\bm{E}t^{\prime}/\hbar}(n_{1})$
$\displaystyle=$
$\displaystyle-\int_{-\infty}^{0}\textrm{d}t^{\prime}\;\frac{\textrm{e}^{t^{\prime}/\tau_{i}}}{\tau_{i}}\partial_{t^{\prime}}\int_{\mathcal{Q}_{2}}{\textrm{d}^{2}k}\;\delta(\bm{k}+e\bm{E}t^{\prime}/\hbar-\bm{Q}_{1})$
$\displaystyle=$
$\displaystyle\int_{-\infty}^{0}\textrm{d}t^{\prime}\;{\textrm{e}^{t^{\prime}/\tau_{i}}}\int_{0}^{\infty}{\textrm{d}k_{y}}\;\delta(k_{y}+eEt^{\prime}\sin\theta/\hbar+\sqrt{3}Q/2)\int_{-\infty}^{k_{y}/\sqrt{3}}{\textrm{d}k_{x}}\;\delta(k_{x}+eEt^{\prime}\cos\theta/\hbar+Q/2)$
$\displaystyle=$
$\displaystyle\int_{-\infty}^{0}\textrm{d}t^{\prime}\;{\textrm{e}^{t^{\prime}/\tau_{i}}}\int_{0}^{\infty}{\textrm{d}k_{y}}\;\delta(k_{y}+eEt^{\prime}\sin\theta/\hbar+\sqrt{3}Q/2)\Theta(k_{y}/\sqrt{3}+eEt^{\prime}\cos\theta/\hbar+Q/2)$
$\displaystyle=$
$\displaystyle\int_{-\infty}^{0}\textrm{d}t^{\prime}\;{\textrm{e}^{t^{\prime}/\tau_{i}}}\int_{\max(-\sqrt{3}eEt^{\prime}\cos\theta/\hbar-\sqrt{3}Q/2,0)}^{\infty}{\textrm{d}k_{y}}\;\delta(k_{y}+eEt^{\prime}\sin\theta/\hbar+\sqrt{3}Q/2)$
$\displaystyle=$
$\displaystyle\int_{-\infty}^{0}\textrm{d}t^{\prime}\;{\textrm{e}^{t^{\prime}/\tau_{i}}}\Theta[-\max(eEt^{\prime}\sin\theta/\hbar+\sqrt{3}Q/2,eEt^{\prime}(\sin\theta-\sqrt{3}\cos\theta)/\hbar)]$
$\displaystyle=$ $\displaystyle\int_{-\infty}^{-\sqrt{3}\hbar
Q/(2eE\sin\theta)}\textrm{d}t^{\prime}\;\textrm{e}^{t^{\prime}/\tau_{i}}$
$\displaystyle=$ $\displaystyle\frac{\textrm{e}^{-\sqrt{3}\hbar
Q/(2eE\tau_{i}\sin\theta)}}{\tau_{i}}$ (S7)

In the zero-doping limit, according to Eq. (S3) and (S4), we also have

$\displaystyle\lim_{n_{h}\to 0}\Sigma_{\alpha}=\lim_{n_{h}\to
0}\partial_{n_{\alpha}}E_{ex}=-3A_{1}n_{\alpha}^{1/2}/2+\mathcal{O}(n_{h}^{3/2})$
(S8a) $\displaystyle\lim_{n_{h}\to
0}\mu_{\alpha}=n_{\alpha}/g+\mathcal{O}(n_{h}^{2})$ (S8b)

Thus the inter-valley scattering term becomes

$\displaystyle\lim_{n_{h}\to
0}\frac{g}{\tau_{x}}\sum_{\beta\neq\alpha}[\Theta(\tilde{\mu}_{\beta}-\tilde{\mu}_{\alpha})\min(\mu_{\beta},\tilde{\mu}_{\beta}-\tilde{\mu}_{\alpha})-\alpha\leftrightarrow\beta]$
$\displaystyle=$ $\displaystyle\lim_{n_{h}\to
0}\frac{g}{\tau_{i}}\sum_{\beta\neq\alpha}[\Theta(\sqrt{n_{\alpha}}-\sqrt{n_{\beta}})\min(n_{\beta}/g,3A_{1}(\sqrt{n_{\alpha}}-\sqrt{n_{\beta}})/2)-\alpha\leftrightarrow\beta]$
$\displaystyle=$
$\displaystyle\frac{1}{\tau_{i}}\sum_{\beta\neq\alpha}\left[\Theta(n_{\alpha}-n_{\beta})n_{\beta}-\alpha\leftrightarrow\beta\right]$
(S9)

In summary, the dynamics equation for $n_{\alpha}$ in the zero-doping limit is
finally simplified to

$\displaystyle\dot{n}_{1}=-\eta
n_{1}+\frac{1}{\tau_{x}}[\Theta(n_{1}-n_{2})n_{2}-1\leftrightarrow 2]$
$\displaystyle\dot{n}_{2}=\eta
n_{1}+\frac{1}{\tau_{x}}[\Theta(n_{2}-n_{1})n_{1}-1\leftrightarrow 2]$

Define the valley polarization strength as $p\equiv(n_{2}-n_{1})/n_{h}$, then
the dynamics equation of $p$ is

$\displaystyle\dot{p}=\eta(1-p)+\tau_{x}^{-1}\left[\Theta(p)(1-p)-\Theta(-p)(1+p)\right]=(\eta+\tau_{x}^{-1}\mathrm{sgn}\;p)-(\eta+\tau_{x}^{-1})p$
(S10)

The steady polarization is gotten by $\dot{p}^{S}=0$, which gives

$\displaystyle p^{S}=\frac{\eta-\tau_{x}^{-1}}{\eta+\tau_{x}^{-1}},\;p^{S}\leq
0\;$ (S11a) $\displaystyle p^{S}=1,\;p^{S}>0.$ (S11b)

The solution is discontinuous at $p^{S}=0$, which gives the critical
electrical field for valley flipping as,

$\frac{E_{c}}{E^{*}}=\frac{1}{\sin\theta\ln(\tau_{x}/\tau_{i})}$ (S12)

where $E^{*}=\sqrt{3}\hbar Q/(2e\tau_{i})$ is the electrical field unit.

At $\theta=\pi/2$, the critical field strength in the zero doping limit Eq.
(S12) is plotted as a function of $\tau_{x}/\tau_{i}$ in Fig. S1(a) by the
blue line. Besides, numerically results by solving the Boltzmann equation at
finite doping concentrations are also plotted for comparisons. By taking the
relaxation times as $\tau_{x}/\tau_{i}=100$ (the gray dashed in Fig. S1(a)),
the critical field strength is also plotted as a function of hole doping
strength $n_{h}$ in Fig. S1(b).

Figure S1: (a) Critical field for valley flipping as a function of
$\tau_{x}/\tau_{i}$ at fixed electrical field direction $\theta=\pi/2$ and the
initial state is $\mathcal{Q}_{1}$-polarized. The blue line for the zero
doping limit is directly got by Eq. (S12) and the orange and red lines at
finite doping concentrations are numerically solved from the Boltzmann
equation. (b) The critical field strength as a function of hole concentrations
at $\tau_{x}/\tau_{i}=100$ (the dashed gray line in (a)).

##  IV Symmetry Restrictions on the Angular Pattern of Current Responses

Assume the system is $\mathcal{Q}_{1}$-polarized and the surviving mirror
symmetry is $M=M_{\pi/3}$ with respect to the mirror line $\theta=\pi/3$.
After the mirror operation, the electrical field and steady hole distribution
transform to $M\bm{E}$ and $n^{S}_{M\bm{k}}$ respectively, and the current
becomes

$\displaystyle\bm{j}^{\prime}=$
$\displaystyle\frac{e}{\hbar}\int\frac{\textrm{d}^{2}k}{(2\pi)^{2}}\;[-\nabla_{\bm{k}}{\varepsilon}^{0}_{\bm{k}}+e(M\bm{E})\times\hat{z}\Omega_{\bm{k}}]n^{S}_{M\bm{k}}$
$\displaystyle=$
$\displaystyle\frac{e}{\hbar}\int\frac{\textrm{d}^{2}k}{(2\pi)^{2}}\;M[-\nabla_{\bm{k}}{\varepsilon}^{0}_{\bm{k}}+e(M\bm{E})\times\hat{z}\Omega_{\bm{k}}]n^{S}_{\bm{k}}$
$\displaystyle=$ $\displaystyle M\bm{j}_{N}-\bm{j}_{H}.$ (S13)

We can see that the normal current $\bm{j}^{\prime}_{N}=M\bm{j}_{N}$
transforms along with the mirror operation while the Hall current
$\bm{j}^{\prime}_{H}=-\bm{j}_{H}$ only changes sign. This reflects the fact
that the Hall current is a pseudo-vector. Besides, after the mirror operation,
the unit vectors in longitudinal and transverse directions also change to
$\hat{r}^{\prime}=M\hat{r}$ and
$\hat{\theta}^{\prime}=\hat{z}\times(M\hat{r})=-M(\hat{z}\times\hat{r})=-M\hat{\theta}$.
Thus the different current components transform as
$j^{\prime}_{N,\parallel}\equiv\hat{r}^{\prime}\cdot\bm{j}^{\prime}_{N}=j_{N,\parallel}$,
$j^{\prime}_{N,\perp}\equiv\hat{\theta}^{\prime}\cdot\bm{j}^{\prime}_{N}=-j_{N,\perp}$
and $j^{\prime}_{H}\equiv\hat{\theta}^{\prime}\cdot\bm{j}^{\prime}_{H}=j_{H}$
respectively. In other words, the electrical current as a function of $\theta$
should satisfy

$\displaystyle j_{N,\parallel}(2\pi/3-\theta)=j_{N,\parallel}(\theta)$ (S14a)
$\displaystyle j_{N,\perp}(2\pi/3-\theta)=-j_{N,\perp}(\theta)$ (S14b)
$\displaystyle j_{H}(2\pi/3-\theta)=j_{H}(\theta).$ (S14c)

In the AC case, since the applied AC electrical field satisfies
$\tilde{\bm{E}}(t;\theta+\pi)=\tilde{\bm{E}}(t+T/2;\theta)$, there is an
additional tempo-spatial symmetry $\bm{j}(t;\theta+\pi)=\bm{j}(t+T/2;\theta)$
which reads

$\displaystyle\bm{j}^{(n)}(\theta+\pi)=$
$\displaystyle\frac{1}{T}\int_{0}^{T}\textrm{d}t\;\bm{j}(t;\theta+\pi)\textrm{e}^{in\omega
t}$ $\displaystyle=$
$\displaystyle\frac{1}{T}\int_{0}^{T}\textrm{d}t\;\bm{j}(t+T/2;\theta)\textrm{e}^{in\omega
t}$ $\displaystyle=$
$\displaystyle\frac{1}{T}\int_{-T/2}^{T/2}\textrm{d}t\;\bm{j}(t;\theta)\textrm{e}^{in\omega
t}(-1)^{n}$ $\displaystyle=$ $\displaystyle(-1)^{n}\bm{j}^{(n)}(\theta)$ (S15)

in the frequency domain. When the electrical field reverses in the second half
period ($t\to t+T/2$, $\theta\to\theta-\pi$), the unit vectors in the
longitudinal and transverse direction also reverse, i.e.,
$\hat{r}\to\hat{r}^{\prime\prime}=-\hat{r}$,
$\hat{\theta}\to\hat{\theta}^{\prime\prime}=-\hat{\theta}$. So the different
current components should satisfy

$\displaystyle
j^{(n)}_{N,\parallel}(\theta+\pi)=\hat{r}^{\prime\prime}\cdot\bm{j}^{(n)}(\theta+\pi)=(-1)^{n+1}j^{(n)}_{N,\parallel}(\theta)$
(S16a) $\displaystyle
j^{(n)}_{N,\perp}(\theta+\pi)=\hat{\theta}^{\prime\prime}\cdot\bm{j}^{(n)}(\theta+\pi)=(-1)^{n+1}j^{(n)}_{N,\perp}(\theta)$
(S16b) $\displaystyle
j^{(n)}_{H}(\theta+\pi)=\hat{\theta}^{\prime\prime}\cdot\bm{j}^{(n)}(\theta+\pi)=(-1)^{n+1}j^{(n)}_{H}(\theta)$
(S16c)